Dispersion relation for hadronic light-by-light scattering: theoretical foundations

In this paper we make a further step towards a dispersive description of the hadronic light-by-light (HLbL) tensor, which should ultimately lead to a data-driven evaluation of its contribution to $(g-2)_\mu$. We first provide a Lorentz decomposition of the HLbL tensor performed according to the general recipe by Bardeen, Tung, and Tarrach, generalizing and extending our previous approach, which was constructed in terms of a basis of helicity amplitudes. Such a tensor decomposition has several advantages: the role of gauge invariance and crossing symmetry becomes fully transparent; the scalar coefficient functions are free of kinematic singularities and zeros, and thus fulfill a Mandelstam double-dispersive representation; and the explicit relation for the HLbL contribution to $(g-2)_\mu$ in terms of the coefficient functions simplifies substantially. We demonstrate explicitly that the dispersive approach defines both the pion-pole and the pion-loop contribution unambiguously and in a model-independent way. The pion loop, dispersively defined as pion-box topology, is proven to coincide exactly with the one-loop scalar QED amplitude, multiplied by the appropriate pion vector form factors.


Introduction
The anomalous magnetic moment of the muon, a µ = (g − 2) µ /2, is one of the very rare quantities in particle physics where a significant discrepancy between its experimental determination and the Standard-Model evaluation still persists. A robust, reliable estimate of the theoretical uncertainties is thus mandatory, especially in view of the projected accuracy of the forthcoming (g − 2) µ experiments at FNAL [1] and J-PARC [2], to settle whether this discrepancy can be attributed to New Physics.
The uncertainty on the theory side is dominated by hadronic contributions, see e.g. [3][4][5], the largest error arising from hadronic vacuum polarization (HVP). Since HVP is straightforwardly related to the total e + e − hadronic cross section via a dispersion relation, the evaluation of this contribution is expected to become more accurate within the next few years [6] thanks to upcoming improvements in the experimental input, 1 although reducing further the present sub-percent accuracy will be challenging. This implies that hadronic light-by-light (HLbL) scattering is likely to soon dominate the theory uncertainty in a µ . 2 Using previous approaches [11][12][13][14][15][16][17][18][19][20][21][22], systematic errors are difficult, if not impossible, to quantify. A novel strategy is required to go beyond the state of the art, to avoid or at least reduce model dependence to a minimum, make a solid estimate of theoretical errors, and possibly reduce them. Lattice QCD is a natural candidate to achieve this goal, but it is not clear yet if or when this method will become competitive [23][24][25]. Alternatively, as we recently argued [26], it is possible to establish a rigorous framework based on dispersion theory that directly links a HLbL µ to experimentally accessible on-shell form factors and scattering amplitudes [27], contrary to what has been previously conjectured [12,24,28,29]. 3 Our dispersive formalism is based on the fundamental principles of unitarity, analyticity, crossing symmetry, and gauge invariance. The derivation of such a systematic framework becomes challenging due to the fact that HLbL scattering is described by a hadronic four-point function whose properties are significantly more complicated than those of the two-point function entering HVP. Besides the opportunity to achieve a datadriven evaluation of a HLbL µ , this approach allows us to unambiguously define and evaluate the various low-energy contributions to HLbL scattering, most notably pion pole and pion loop.
The key result of [26] was a master formula giving the contribution to a HLbL µ from the pion-pole term and intermediate ππ states expressed in terms of the pion transition form factor and helicity partial waves for γ * γ * → ππ. In the present paper we reformulate this approach in a more general context, by employing a different generating set for the Lorentz structures of the HLbL tensor. This new set, which we constructed following the prescription by Bardeen, Tung [31], and Tarrach [32] (BTT), has the property of being explicitly free of kinematic singularities and zeros [33]. Therefore, the scalar coefficient functions associated with each Lorentz structure of this generating set are also free of kinematic singularities and zeros, and thus fulfill a Mandelstam double-dispersive representation. Further motivation for deriving such a BTT set is provided by the fact that the role of gauge invariance and crossing symmetry becomes fully transparent, with constraints from soft-photon zeros incorporated automatically. Moreover, the absence of kinematic singularities makes it much easier to perform the angular averages necessary to calculate the contribution to (g − 2) µ -indeed both the derivation and the final form of the master formula are simplified considerably in the present setting, and allow for a simpler inclusion of partial-wave contributions beyond S-waves. In the derivation of the master formula it is very useful to perform a Wick rotation of the integration variables. This has been already used in several papers in the literature without ever addressing an important subtlety which is relevant in the present case: namely if the Wick rotation goes through without changes even if the integrand is a loop function showing anomalous thresholds. We discuss this point in detail and prove that this is indeed the case. Finally, the BTT formulation facilitates the proof that both pion-pole and pion-loop contributions are uniquely defined in the dispersive approach, and confirms the relation to the pion transition form factor and the pion electromagnetic form factor given in [26]. We explicitly proof that the scalar QED (sQED) pion loop dressed with pion form factors (denoted by FsQED in [26]) is identical to the contribution from ππ intermediate states with a pion-pole left-hand cut (LHC).
In summary, this paper contains three main results: (i) the explicit construction of the BTT decomposition of the HLbL tensor; (ii) the derivation of the master formula expressing the HLbL contribution in terms of the BTT scalar functions; (iii) the proof that in the dispersive approach the FsQED contribution is uniquely defined. Everything which goes beyond the FsQED contribution is amenable to partial-wave expansion and will be discussed in full detail in a forthcoming publication.
The present paper is structured as follows. In Sect. 2, we first review some aspects of the process γ * γ * → ππ, mainly as an illustration of the techniques that we apply afterwards to HLbL. In Sect. 3, we derive the decomposition of the HLbL tensor into a set of scalar functions that are free of kinematics. In Sect. 4, this BTT decomposition is then used to derive a master formula for a HLbL µ that is parametrized by the corresponding scalar functions. We verify the validity of the Wick rotation of the two-loop integral in the calculation of a HLbL µ , even in the case of anomalous thresholds in the HLbL tensor. In Sect. 5, we derive a Mandelstam representation of the BTT scalar functions, and, by comparison to the loop formulation of the sQED pion loop, demonstrate that it indeed coincides with the contribution of ππ intermediate states with a pion-pole LHC. We conclude with an outlook in Sect. 6, while some of the lengthier expressions and derivations are collected in the appendices. 1 For a recently suggested alternative approach to determine the leading hadronic correction to aµ, based on space-like data extracted from Bhabha scattering, see [7]. 2 At this order in the fine-structure constant also two-loop diagrams with HVP insertion appear [8]. Higher-order hadronic contributions have been investigated in [9,10]. 3 For a different approach aiming at a data-driven evaluation of a HLbL µ , based on a dispersive description of the Pauli form factor instead of the HLbL tensor, see [30]. 2 The sub-process γ * γ * → ππ As a prelude, we discuss in this section the process γ * γ * → ππ, which will become important as a sub-process when we write down a dispersion relation for the HLbL tensor. Because the Lorentz structure of γ * γ * → ππ is much simpler than the one of light-by-light scattering, it also allows us to illustrate a technique for the construction of amplitudes that are free of kinematic singularities and zeros, which we will apply afterwards to the more complicated case of HLbL. In order to render this analogy manifest, we slightly change conventions compared to [26]: we identify here γ * γ * → ππ as the s-channel process instead of γ * π → γ * π. This assignment is slightly unnatural when it comes to constructing dispersion relations for a process with these crossing properties, see [34][35][36][37], but makes the inclusion into the unitarity relation in HLbL more straightforward.

Kinematics and matrix element
Consider the process shown in Fig. 1. At O(e 4 ), the amplitude for this process is given by where ξ is an arbitrary gauge parameter for the photon propagators and the tensor W ab µν is defined as the pure QCD matrix element The contraction thereof with appropriate polarization vectors may be understood as an amplitude for the off-shell process where λ 1,2 denote the helicities of the off-shell photons. We define the connected part of such a matrix element by The helicity amplitudes are given by the contraction with polarization vectors: We introduce the following kinematic variables: 4 which satisfy s + t + u = q 2 1 + q 2 2 + 2M 2 π .

Tensor decomposition
The tensor W µν ab can be decomposed based on Lorentz covariance as (we drop isospin indices) where we abbreviate q i = {q 1 , q 2 , p 2 − p 1 } and where double indices are summed. The ten coefficient functions {W 1 , W ij 2 } cannot contain any kinematic but only dynamic singularities. However, they have to fulfill kinematic constraints that are required e.g. by gauge invariance, hence they contain kinematic zeros. Conservation of the electromagnetic current implies the Ward identities hence a priori six relations between the scalar functions {W 1 , W ij 2 }. Only five of them are linearly independent and thus reduce the set of scalar functions to five independent ones.
We now construct a set of scalar functions which are free of both kinematic singularities and zeros. We follow the recipe given by Bardeen,Tung [31], and Tarrach [32]. As was shown in [32], the basis (consisting of five functions) constructed according to the recipe of [31] is not free of kinematic singularities. However, a redundant set of six structures can be constructed, which fulfills the requirement.
In the case of γ * γ * → ππ, there exists in addition to gauge invariance and crossing symmetry of the photons also the crossing (Bose) symmetry of the pions, which is responsible for additional kinematic zeros. In fact, this can be used to circumvent the introduction of a sixth function, reducing the set again to five scalar functions free of both kinematic singularities and zeros. This was first derived in [38] in the context of doubly-virtual nucleon Compton scattering. 5 We define the projector which satisfies i.e. the tensor W µν is invariant under contraction with the projector, but contracting the projector with any Lorentz structure produces a gauge-invariant structure. We apply this projector for both photons: As the functionsB i are a subset of the original scalar functions, they are still free of kinematic singularities, but contain zeros, because the Lorentz structures contain singularities. We have to remove now the single and double poles in q 1 · q 2 from the Lorentz structuresT µν i . This is achieved as follows [31]: • remove as many double poles as possible by adding to the structures linear combinations of other structures with non-singular coefficients, • if no more double poles can be removed in this way, multiply the structures that still contain double poles by q 1 · q 2 , • proceed in the same way with single poles.
It turns out that no double poles in q 1 · q 2 can be removed by adding to the structures multiples of the other structures, henceT µν 2 , . . . ,T µν 5 have to be multiplied by q 1 · q 2 . The resulting simple poles can be removed by adding multiples ofT µν 1 . In the end, we have to multiplyT µν 1 by q 1 · q 2 in order to remove the last pole. We then arrive at the following representation: (2.14) whereT (2.16) The Lorentz structures (2.15) agree with the basis used in [26], upon q 2 → −q 2 , p 1 → −p 1 due to crossing and the identificatioñ T 1,µν = T [26] 1,µν ,T 2,µν = T [26] 3,µν ,T 3,µν = As shown by Tarrach, this basis is not free of kinematic singularities and zeros [32]. The structuresT µν i form a basis for q 1 · q 2 = 0, but are degenerate for q 1 · q 2 = 0. This degeneracy implies that there is a linear combination of the structuresT µν i that is proportional to q 1 · q 2 : Hence, in order to have a generating set even at q 1 · q 2 = 0 this sixth structureT µν 6 has to be added by hand. Projected on the basis, it gives coefficients with poles in q 1 · q 2 . Although the basis is not free of kinematic singularities and zeros, we have found the exact form of the singularities: The functions B i are free of kinematic singularities. The tensor structures can be written in terms of the Mandelstam variables t and u as The derivation of the Lorentz decomposition for γ * γ * → ππ up to this point should illustrate the techniques that will be applied to the much more complicated case of HLbL in Sect. 3. In the case of γ * γ * → ππ, there is in fact a solution to get rid of the additional redundant structure without introducing kinematic singularities [38]. The necessary steps are discussed in the next subsection. Note that it is an open question whether an analogous procedure to get rid of all redundancies in HLbL exists, see the discussion below (3.29).

Kinematic zeros due to crossing antisymmetry
We have now constructed a redundant set of six Lorentz structures. All kinematic constraints from gauge invariance are implemented and the scalar coefficient functions B i are free of kinematic singularities. However, due to the crossing symmetry of the pions, additional kinematic zeros are present in the functions B i . As noted in [38], this can be exploited to eliminate the redundancy and to work again with basis coefficient functions free of kinematic singularities and zeros.
Since the two-photon state is even under charge conjugation, so must be the two-pion state. Therefore, the isospin I = 1 amplitude vanishes. Bose symmetry implies that the amplitude and hence the tensor W µν is invariant under p 1 ↔ p 2 or equivalently q 3 ↔ −q 3 . Under this transformation, the tensor structuresT µν i are even for i = 1, 2, 5, 6 and odd for i = 3, 4. Hence, B 1,2,5,6 are even, while B 3,4 must be odd. They contain a kinematic zero of the form whereB 3,4 are free of kinematic singularities. Crossing symmetry of the photons requires the invariance of W µν under q 1 ↔ q 2 , µ ↔ ν. WhileT µν 1,2,5,6 are invariant under this transformation, we observe the crossing relationT µν 3 ↔T µν 4 . Hence, for fixed Mandelstam variables, B 1,2,5,6 andB 3 −B 4 are even under q 2 We use the notation The scattering angle is then given by . (2.30) We define the polarization vectors by The longitudinal polarization vectors are normalized to one for ξ i = q 2 i . However, as the off-shell photons are not physical states, the choice of ξ i has no influence on any physical observable. The fact that the dependence on the normalization has to drop out in the end provides a useful check on the calculation, so that we will keep ξ i general in the following.
The helicity amplitudes are then given by where λ 12 (s) := λ(s, q 2 1 , q 2 2 ). Since the functions A i are free of kinematic singularities and zeros, we can read off from these equations the soft-photon zeros [39,40]. In the limit q 1 → 0, the Mandelstam variables become We conclude that the helicity amplitudes must vanish at this point apart from terms containing a dynamic singularity, which will be discussed in Sect. 2.6. The second soft-photon limit, q 2 → 0, leads to 34) and the same arguments apply for the helicity amplitudes. Crossing symmetry of the photons implies that under the transformation q 2 1 ↔ q 2 2 (and fixed Mandelstam variables), H ++ , H +− , and H 00 remain invariant, the other two helicity amplitudes transform as H +0 ↔ H 0+ .

Partial-wave expansion
To make the connection to [26], we briefly discuss the consequences of the decomposition (2.32) for the kernel functions that appear in partial-wave dispersion relations for γ * γ * → ππ. Partial waves are most conveniently defined for helicity amplitudes. Using the formalism of [41], we write the partial-wave expansions as where d J m0 is the Wigner d-function, m = |λ 1 − λ 2 |, and the helicity partial waves h J,λ1λ2 depend implicitly on the photon virtualities q 2 1 and q 2 2 . Since the isospin of the two-pion system is I = 0, 2, only even partial waves are allowed. For m = 0, i.e. H ++ and H 00 , the partial-wave expansion starts at J = 0, otherwise at J = 2. Note also that d J 00 (z) = P J (z) are the Legendre polynomials. The partial waves can be obtained by projection: Since the functions A i are free of kinematic singularities and zeros, they are appropriate for a dispersive description. To this end, the Roy-Steiner treatment of [36] can be generalized to the doubly-virtual case. One starts by writing down hyperbolic dispersion relations and K ij JJ are integral kernels. The ellipses in (2.38) stand for the contribution of the partial waves of the crossed channels. The normalization of the helicity amplitudes is reabsorbed in (2.39), so that the partial waves h J,i remain finite in the limit q 2 i → 0. If only S-waves are taken into account, the relation between the scalar functions and the helicity partial waves in the s-channel is given by In this case, the scalar functions depend only on s. If we take into account D-waves as well, the scalar functions A 3 , A 4 , and A 5 no longer vanish, but they depend only on s, while the scalar functions A 1 and A 2 are now second order polynomials in t and u. The explicit expressions are given in App. A.
For the diagonal kernel functions in (2.38), we find: Because we use here a different basis, the kernel functions are slightly different from the ones given in [26]: in order to avoid fictitious poles in q 2 1 − q 2 2 , one has to work with the linear combinations h 2,+0 ± h 2,0+ , which are symmetric/antisymmetric under q 2 1 ↔ q 2 2 . The remaining difference to [26] is a polynomial, which can be reabsorbed into subtraction constants.

Pion-pole contribution
In order to define the contribution to γ * γ * → π + π − due to the exchange of a single pion, we employ a dispersive picture. In this subsection, we demonstrate in detail that the pure pion-pole contribution in the sense of unitarity is exactly the Born contribution in scalar QED (sQED), multiplied by pion vector form factors. The dispersive picture cleanly defines this as the correct dependence on the photon virtualities q 2 i . We also discuss the different meaning of Feynman and unitarity diagrams in this context.
To keep the presentation as transparent as possible, we omit subtleties due to isospin conventions here and simply present the plain sQED expressions. Compared to [26], this results in an overall sign for the amplitude of γ * γ * → π + π − . For HLbL, this convention has no effect, but the choice in [26] is advantageous when transforming to the isospin basis and reconstructing partial waves in accordance with Watson's theorem [42].
We assume that the asymptotic behavior of γ * γ * → π + π − in the crossed Mandelstam variables t and u permits an unsubtracted fixed-s dispersion relation 6 for the scalar functions: whereρ s i;t,u denote the pole residues andD s i;t,u the discontinuities along the t-and u-channel cuts. Both are determined by unitarity. Consider the t-channel unitarity relation: where F V π is the electromagnetic form factor of the pion, defined by Figure 2: Correspondence of the pion-pole contribution to γ * γ * → π + π − in terms of unitarity diagrams and the Born contribution in terms of sQED Feynman diagrams. The dashed lines in the unitarity diagrams indicate a cut line, hence the internal pion is on-shell.
After performing the trivial integral, we obtain t-and u-channel are related by p 1 ↔ p 2 : Note that these expressions are only gauge-invariant due to the presence of the delta function: the contraction of the bracket in (2.46) with q µ 1 or q ν 2 is proportional to t − M 2 π . If we project the imaginary parts of W µν onto the scalar functions, making use of the delta function, we obtain while the one-pion contributions to the imaginary parts of the remaining scalar functions vanish. The pion-pole contribution to the scalar functions is therefore We compare the pion-pole contribution with the Born contribution in sQED: and read off the Born values of the scalar functions: (2.51) We find that the pion-pole contribution corresponds exactly to the sQED Born contribution multiplied by electromagnetic pion form factors for the two off-shell photons. 7 Note that, if we think in terms of unitarity diagrams, we have now considered the pure pole contribution to the scalar functions. However, in terms of Feynman diagrams in sQED this corresponds to a sum of two pole diagrams and the seagull diagram. 8 It is important to be aware of the different meaning of a topology in the sense of unitarity and a Feynman diagram, see Fig. 2. As will be shown in Sect. 5, it is exactly this distinction that makes the sQED pion loop in HLbL coincide with box-type unitarity diagrams representing ππ intermediate states with a pion-pole LHC, although, in terms of Feynman diagrams, it is composed of the sum of box, triangle, and bulb topologies.
3 Lorentz structure of the HLbL tensor

Definitions
In order to study the contribution of HLbL scattering to the anomalous magnetic moment of the muon, we need first of all a description of the HLbL tensor. The object in question is the hadronic Green's function of four electromagnetic currents, evaluated in pure QCD (i.e. with fine-structure constant α = e 2 /(4π) = 0): The electromagnetic current includes only the three lightest quarks: . The contraction of the HLbL tensor with polarization vectors gives the hadronic contribution to the helicity amplitudes for (off-shell) photon-photon scattering: For notational convenience, we define The kinematics is illustrated in Fig. 3. We use the following Lorentz scalars as kinematic variables -these are the usual Mandelstam variables: which fulfill (we will take k 2 = 0 at some later point) Gauge invariance requires the HLbL tensor to satisfy the Ward-Takahashi identities 7 Therefore, the dispersive definition of the pion pole (2.49) coincides with the gauge-invariant pole contribution of the 'softphoton amplitude' in [43]. We thank S. Scherer for pointing this out. 8 The equivalence of the pion pole and the Born term is surprising given the fact that (2.50) contains a term with g µν , while the imaginary parts (2.46) and (2.47) do not. Tracing the above steps backwards, one sees that in the t-or u-channel imaginary parts the coefficient of and hence vanishes due to the delta function.

Tensor decomposition
In general, the HLbL tensor can be decomposed into 138 Lorentz structures [13,44,45]: The 138 scalar functions depend on six independent kinematic variables, e.g. on two Mandelstam variables s and t and the virtualities q 2 1 , q 2 2 , q 2 3 , and q 2 4 . They are free of kinematic singularities but contain kinematic zeros, because they have to fulfill kinematic constraints required by gauge invariance. The Ward identities (3.7) impose 95 linearly independent relations on the scalar functions, reducing the set to 43 functions.
As we did in Sect. 2.2 for the case of γ * γ * → ππ, we will now construct a set of Lorentz structures and scalar functions, such that the scalar functions contain neither kinematic singularities nor zeros. Compared to γ * γ * → ππ, the application of the recipe given by Bardeen, Tung [31], and Tarrach [32] is much more involved. Again, the recipe by Bardeen and Tung does not lead to a kinematic-free minimal basis (which would consist here of 43 scalar functions). 9 Following Tarrach, we will construct a redundant set of 54 structures, which is free of kinematic singularities and zeros.
In a first step, we define the two projectors which have the following properties: i.e. the HLbL tensor is invariant under contraction with the projectors, but the contraction of every Lorentz structure produces a gauge-invariant structure. Hence, we project the tensor Only 43 of the 138 projected structuresL µνλσ i are non-zero, i.e. all constraints imposed by gauge invariance are already manifestly implemented. Since the projected structures are still multiplied by the original scalar functions Ξ i , no kinematic singularities have been introduced into the scalar functions. We now have to remove the kinematic zeros from the scalar functions by removing the single and double poles in q 1 · q 2 and q 3 · q 4 , which are present in the structuresL µνλσ i . We adapt the recipe of [31] (cf. Sect. 2.2): • remove as many (q 1 · q 2 , q 3 · q 4 ) double-double poles as possible by adding to the structures linear combinations of other structures with coefficients containing no poles, • if no more double-double poles can be removed in this way, multiply the structures that still contain double-double poles by either q 1 · q 2 or q 3 · q 4 (the choice is irrelevant in the end), • proceed in the same way with double-single, single-double poles, etc. until no poles at all are left in the structures.
As already mentioned, it is again impossible to avoid introducing kinematic singularities into the scalar functions by applying this procedure [32]. However, the only step where kinematic singularities can be introduced is the multiplication of the structures by q 1 · q 2 or q 3 · q 4 (i.e. the division of the scalar functions by these terms). This means that the only possible singularities are (double or single) poles in q 1 ·q 2 or q 3 ·q 4 . The precise form of these poles can be easily determined: they correspond to degeneracies of the obtained basis of Lorentz structures in the limit q 1 · q 2 → 0 and/or q 3 · q 4 → 0. Therefore, the 43-dimensional basis has to be extended by additional structures, which are found by studying the null-space of the present structures in the mentioned limits. 11 such structures can be found. The extended generating set of 54 structures exhibits all possible crossing symmetries in a manifest way. Explicitly, the resulting representation of the HLbL tensor reads These structures satisfy the following crossing symmetries: where the crossing operators C ij exchange momenta and Lorentz indices of the photons i and j, e.g. 10 All the remaining structures are just crossed versions of the above seven structures, as shown in App. B.1.
Since the HLbL tensor Π µνλσ is totally crossing symmetric, the scalar functions Π i have to fulfill exactly the same crossing properties of the corresponding Lorentz structures as given in (3.15) (note the antisymmetric crossing relations in Π 49 ). Therefore, only seven different scalar functions Π i appear, together with their crossed versions. These scalar functions are free of kinematic singularities and zeros and hence fulfill a Mandelstam representation. 11 They are suitable quantities for a dispersive description. The subset consisting of the following 43 Lorentz structures forms a basis: The corresponding scalar coefficient functionsΠ i , defined by exhibit kinematic singularities in q 1 · q 2 and q 3 · q 4 . The exact form of these kinematic singularities can be determined by projecting (3.13) on this basis. The relation between the basis coefficient functionsΠ i and the scalar functions Π i is given in App. B.2.
While crossing symmetry is manifest if the HLbL tensor is expressed in terms of the redundant set of 54 structures, it is partially obscured in the basis of 43 elements. Still, many basis functionsΠ i are related by crossing symmetry, see (B.3) in App. B.2. The basis coefficient functions are completely specified by crossing symmetry and 9 representative functions, e.g.
We stress again that our decomposition (3.13) is both manifestly crossing symmetric and gauge invariant. In particular, this implies that all soft-photon zeros [39] of the helicity amplitudes (3.3) are already implemented, in contrast to the helicity basis employed in [26]. In this case, only specific combinations of basis functions were free of kinematic singularities and zeros and thus admitted a dispersive representation, while the softphoton constraints needed to be imposed by hand. The BTT set (3.14) provides a Lorentz decomposition that incorporates soft-photon constraints and avoids kinematics by construction. On the other hand, the main motivation for choosing the helicity basis was that the unitarity relations in the form of helicity partial-wave amplitudes become diagonal. Of course, this will not be true in the BTT basis, so that the derivation of an explicit dispersive representation for the BTT amplitudes requires a basis change. In practice, this change between BTT and helicity bases proceeds by means of (3.8) as an intermediate step: the helicity amplitudes can be easily expressed in terms of (3.8), while the remaining basis change can be performed along the lines explained in App. C.

Redundancies in the Lorentz decomposition
Because the set of 54 scalar functions is redundant, there is obviously an ambiguity in the definition of the coefficient functions. However, crossing symmetry and the absence of kinematic singularities restrict this ambiguity in a very specific way. The shifts leaving the HLbL tensor Π µνλσ invariant are easily found by requiring that the basis elementsΠ i remain invariant. We find where the shifts have to satisfy The ellipses in the last line stand for the permutations that are not listed explicitly: ∆ 49 is symmetric (antisymmetric) under all even (odd) crossing permutations. The shifts δΠ i that are not specified follow from the crossing relations (B.1). Apart from the redundancy introduced by extending the basis to a set of 54 elements, there exists another ambiguity in the definition of the scalar function. As it was pointed out in [46], the 138 initial Lorentz structures of the HLbL tensor fulfill two linear relations in d = 4 space-time dimensions. This is not obvious in a covariant notation, but can be easily shown in an explicit reference frame. The two linear relations can be translated into two linear relations between the 43 basis elements (which means that these structures actually do not form a basis in a strict sense). Therefore, in four space-time dimensions there are only 41 independent scalar functions, which matches the number of fully off-shell helicity amplitudes. This is analogous to the case of γ * γ * → ππ, where the number of helicity amplitudes (5) also matches the number of independent scalar functions.
The two linear relations can be further translated into relations amongst the 54 elements of the redundant crossing-symmetric set of Lorentz structures. In fact, it turns out that expressed in this way, the two relations become particularly simple. They are given by: where and where f 1,2 are a priori unspecified functions free of kinematic singularities. The requirement that the scalar functions satisfy crossing relations of the same form as (3.15) and (B.1) imposes crossing relations on the functions f 1,2 : Let us write The first bracket in this expression is given by Π 49 − Π 50 − Π 53 + Π 54 and has exactly the crossing properties of f 1 . Therefore, the ambiguity can be fixed by the condition Very recently, a Lorentz decomposition of the HLbL tensor into 72 redundant structures has been presented in [47]. A comparison with our redundant set of 54 structures, as first given in [33], has been performed, confirming the absence of kinematic singularities. By studying systematically the irreducible representations of the permutation group S 4 , the authors of [47] aim at constructing a minimal basis consisting of 41 totally crossing symmetric elements free of kinematic singularities. Note, however, that [47] neither provides this basis nor a concrete recipe for how to obtain it. Therefore, it remains to be seen if a minimal basis free of kinematic singularities exists at all. For our purpose the question whether such a basis exists is immaterial, since for a dispersive representation of the HLbL tensor all that is required is a decomposition into scalar amplitudes free of kinematic zeros and singularities.
To conclude this subsection, we stress that the presence of kinematic singularities in any known basis requires us to work with a redundant set of 54 Lorentz structures and scalar functions Π i . Although the redundancy of the Lorentz structures introduces an ambiguity in the functions Π i (which is already very restricted by crossing symmetry and will be further restricted by analyticity), this ambiguity will cancel in the calculation of any physical quantity such as (g − 2) µ . In particular, this implies that the additional redundancy in d = 4 dimensions does not affect our calculation in any way, since the minimal set suited for a dispersive representation still involves 54 functions. In the rest of the paper we will continue to work with the BTT amplitudes defined in (3.14).
4 HLbL contribution to (g − 2) µ In Sect. 4.1, we review the definition and calculation of a µ , as well as general techniques for the calculation of a HLbL µ (see e.g. [48]). The well-known general formula requires still a rather long calculation before a number can be finally obtained. For the pion-pole contribution, these steps have been worked out long ago [19]. With our complete set of 54 kinematic-free structures (3.13), this procedure can be repeated for the whole HLbL contribution in full generality, as we explain in Sects. 4.2 and 4.3. In Sect. 4.4, we study the validity of the applied Wick rotation in the presence of anomalous thresholds.

Projector techniques
The T -matrix element of the interaction of a muon with the electromagnetic field is defined by where e = |e|. Diagrammatically Assuming parity conservation, the vertex function Γ µ can be decomposed into form factors as where k = p 1 − p 2 and σ µν := i 2 [γ µ , γ ν ]. F 1 is the electric charge or Dirac form factor, F 2 the magnetic or Pauli form factor. The anomalous magnetic moment of the muon is given by We expand the vertex function Γ µ to first order in powers of k µ : where p := 1 2 (p 1 + p 2 ). Using projector techniques together with angular averages (see [19,48]), the anomalous magnetic moment of the muon can be calculated as where now p 2 = m 2 µ . 12 We are interested in the contribution of the HLbL tensor to a µ , diagrammatically The HLbL tensor is defined in (3.1). Differentiating the fourth Ward identity in (3.7) with respect to k ρ = (q 1 + q 2 + q 3 ) ρ yields It was already argued in [49] that Π µνλσ vanishes linearly with k (i.e. the derivative contains no singularity), and so must Γ HLbL σ . This is easily verified with our tensor decomposition (3.13). Therefore, the HLbL contribution to the anomalous magnetic moment is given by where Γ HLbL We use the Ward identity (4.9) to write Taking the derivative and limit leads to the well-known expression (4.13)

Loop integration
In order to compute the contribution to (g − 2) µ , one has to take the trace in (4.10) and perform the two-loop integral of equation (4.13). Five of the eight integrals can be carried out analytically with the help of Gegenbauer polynomial techniques [50]. To this end, we employ the representation of the HLbL tensor in terms of the 54 14) It turns out that there are only 19 independent linear combinations of the structures T µνλσ i that contribute to (g − 2) µ . It is possible to make a basis change in the 54 structures  15, 18, 23, 24, 27, 28, 29, 30, 32, 35, 36, 37, 38, 41, 44, 45, 48, 49, 52, 53 (4.17) are irrelevant for the calculation of (g − 2) µ . The HLbL contribution to (g − 2) µ can now be written as whereT i (q 1 , q 2 ; p) := and theΠ i are needed for the reduced kinematics s = (q 1 + q 2 ) 2 , t = q 2 2 , u = q 2 1 , q 2 1 , q 2 2 , q 2 3 = (q 1 + q 2 ) 2 , k 2 = q 2 4 = 0. (4.20) The explicit result of the trace calculation and the contraction of the Lorentz indices is given in App. E.1. We can reduce the number of terms contributing to (g − 2) µ further by using the symmetry under the exchange of the momenta q 1 ↔ −q 2 : the loop integration measure and the product of propagators are invariant under this transformation, while the kernelsT i transform under q 1 ↔ −q 2 aŝ For the reduced kinematics (4.20) the exchange q 1 ↔ −q 2 is equivalent to the crossing transformation t ↔ u, With the help of the crossing relations of the scalar functions Π i , it is easy to check that theΠ i transform analogously to the kernelsT i , i.e.
Therefore, it is convenient to write the HLbL contribution to (g − 2) µ as a sum of 12 terms: Note that the first two terms in this sum correspond to the well-known result for the pion-pole contribution [19] (up to conventions: exchange ofT 1 andT 2 , the explicit factor ξ 2 = 2, and symmetrization ofT 1 ). In (4.23), the integrand depends on the five scalar products q 2 1 , q 2 2 , q 1 ·q 2 , p·q 1 , and p·q 2 , where the dependence on the last two is given explicitly (the scalar functions only depend on q 2 1 , q 2 2 , and q 1 · q 2 ). Therefore, five of the eight integrals can be performed without knowledge of the scalar functions. The same integrals as in the case of the pion-pole contribution occur [3,19], which have been solved with the technique of Gegenbauer polynomials [50]. This method has been applied before to the full HLbL contribution in the context of vectormeson-dominance and hidden-local-symmetry models [51,52].
We perform a Wick rotation of the momenta q 1 , q 2 , and p (see Sect. 4.4) and denote the Wick-rotated Euclidean momenta by capital letters Q 1 , Q 2 , and P . Note that Q 2 is a pure number, it does not depend on the direction of the momentum P of the muon, hence we can take the angular average by integrating over the four-dimensional hypersphere: a HLbL µ = dΩ 4 (P ) 2π 2 a HLbL µ . (4.25) The kernelsT i are at most quadratic in p, therefore we need the following angular integrals [3]: where τ = cos θ 4 , defined by Q 1 · Q 2 = |Q 1 ||Q 2 |τ , is the cosine of the angle between the Euclidean four-momenta Q 1 and Q 2 , and further

Master formula
After using the angular integrals (4.26), we can immediately perform five of the eight loop integrals by changing to spherical coordinates in four dimensions. This leads us to a master formula for the HLbL contribution to the anomalous magnetic moment of the muon: where Q 1 := |Q 1 |, Q 2 := |Q 2 |. The hadronic scalar functionsΠ i are just a subset of theΠ i and defined in (D.2). They have to be evaluated for the reduced kinematics The integral kernels T i listed in App. E.2 are fully general for any light-by-light process, while the scalar functions Π i parametrize the hadronic content of the master formula. In particular, (4.28) can be considered a generalization of the three-dimensional integral formula for the pion-pole contribution [3]. It is valid for the whole HLbL contribution and completely generic, i.e. it can be used to compute the HLbL contribution to (g − 2) µ for any representation of the HLbL tensor, irrespective of whether the scalar functions are subsequently specified dispersively or taken from a model calculation. For an arbitrary representation of the HLbL tensor, the scalar functions Π i can be obtained by projection, see App. C and App. F.2. An analogous master formula was derived in [26] in the case of a helicity basis for the HLbL tensor. This step in the calculation is completely equivalent, in both cases the calculation of a HLbL µ proceeds via an evaluation of the trace in (4.10) and subsequent reduction of the two-loop integral with Gegenbauer techniques. From a technical point of view, the BTT approach offers several simplifications, since the D-wave-related amplitudes in [26] require another angular average to define the k → 0 limit as well as more complicated Gegenbauer integrals than the standard ones given in (4.26).
In close analogy to the pion-pole contribution [19], the main benefit of the master formula (4.28) is the fact that it contains only a three-dimensional integral, and thus is well-suited for a direct numerical implementation. In particular, the energy regions generating the bulk of the contribution can be identified by numerically integrating over τ and plotting the integrand as a function of Q 1 and Q 2 [19,[51][52][53].
Before turning to the main part of this paper, the foundations for a model-independent calculation of the scalar functions Π i by making use of dispersion relations, we next consider a subtlety in the application of the Gegenbauer integrals in (4.26). To obtain these integrals a Wick rotation has been performed, whose validity might be questioned if the HLbL tensor possesses non-trivial analytic properties. For instance, already a simple triangle loop function would give rise to anomalous thresholds due to the integration over all photon virtualities, see [54][55][56]. In the following subsection, we prove that the Wick rotation remains valid even in the presence of anomalous thresholds.

Wick rotation and anomalous thresholds
Performing a Wick rotation before calculating a loop integral is a standard procedure. The applicability of this transformation depends on whether the integrand is free of singularities in the region swept by the integration path during the rotation. In our case we have to perform the two-loop integral (4.23), with an integrand containing products of propagators multiplying functions satisfying a Mandelstam representation, so having the analytic properties of one-loop functions. The case of the propagators is standard. The one-loop functions occurring in our integrand, however, need some more discussion. The complication is not due to the presence of cuts instead of poles: cuts of loop functions are usually on the positive real axis above some threshold, and the functions are evaluated on the upper rim of the cut. The Wick rotation can be applied without problem even in these cases. The need for a more detailed discussion is due to the presence of anomalous cuts -cuts which appear first in three-point functions and which, depending on the value of the momenta squared of the external legs, may intrude into the first Riemann sheet away from the real axis. In the present paper we need to deal with three-and four-point one-loop functions, such as C 0 (s; q 2 1 , q 2 2 ) and D 0 (s, t; q 2 1 , q 2 2 , q 2 3 , q 2 4 ) (see App. G). For a discussion of the analytic properties of these loop functions, see [55,[57][58][59][60]. We will illustrate what happens in the presence of such anomalous cuts by discussing the case of the C 0 function: the conclusion is that the Wick rotation can be applied without problem even in these cases.
The function C 0 (s; q 2 1 , q 2 2 ) has an anomalous cut for q 2 1 > 0, q 2 2 > 0, and q 2 1 + q 2 2 > 4m 2 , where m is the mass of the particle running in the loop (assuming it is only one mass). The origin of this anomalous cut is that C 0 admits a dispersive representation having a discontinuity on the positive real axis for s ≥ 4m 2 . The discontinuity, however, has itself a cut between the branch points: While for q 2 1,2 < 0 these branch points lie on the second Riemann sheet, for q 2 1 > 0, q 2 2 > 0, and q 2 1 +q 2 2 > 4m 2 one of them moves into the first Riemann sheet by crossing the real axis above s = 4m 2 . In this case the dispersive representation of the function C 0 needs to be modified by adding a further integral along the anomalous cut from 4m 2 up to s + . The most convenient way to see how the branch points move as a function of q 2 1 and q 2 2 is to introduce the variables [55] The upper half-plane of the complex variables q 2 i are mapped onto the upper unit semicircle of unit radius: the negative real axis of q 2 i is mapped onto the segment 0 < ξ i ≤ 1, the region 0 ≤ q 2 i ≤ 4m 2 onto the semicircle ξ i = e iφi , with 0 ≤ φ i ≤ π, and finally the semiaxis 4m 2 ≤ q 2 i onto the segment −1 ≤ ξ i < 0. The branch points s ± can be expressed as follows in terms of these variables and from this it is easy to check that s − always stays on the second sheet, whereas s + moves into the first one under the conditions given above. It is important to stress that when it enters the first Riemann sheet, s + stays always below the real axis, either on its lower rim (for q 2 1 < 4m 2 and q 2 2 < 4m 2 or for q 2 1 > 4m 2 and q 2 2 > 4m 2 ), or into the lower complex half-plane (for q 2 1 < 4m 2 and q 2 2 > 4m 2 or vice versa), in contrast to what is shown in Fig. 2 of [55].
When we perform the Wick rotation for q 1 and q 2 (to avoid weird trajectories in the complex s-plane we rotate both simultaneously), their squares (in Minkowski space) become transformed into minus their Euclidean squares (−Q 2 i ), and can only take negative values. But also s = (q 1 + q 2 ) 2 is mapped onto s E ≡ −(Q 1 + Q 2 ) 2 , which equally can only be negative. If we start from a positive value of s and perform the Wick rotation, the trajectory of s follows an arc in the upper complex half plane and lands on the negative axis at s E . The potential danger is if we start from s > 4m 2 , just on the upper rim of the cut, but the Wick rotation moves it away both from the normal as well as from the anomalous cut, which, as mentioned above, is always below the Re(s) Im(s) real axis. If we start from negative values of s we also follow an arc in the upper complex half plane and land at s E < s, without ever coming even close to the cuts.
The discussion above is still incomplete, however, because we have so far neglected the fact that as we make the Wick rotation not only s moves, but also the q 2 i and with them s + , which determines where the anomalous cut is. In principle, s + could move into the complex upper half-plane and interfere with the trajectory of s during the Wick rotation. An explicit calculation shows, however, that this does not happen: at the end of the Wick rotation both q 2 i have become −Q 2 i and therefore negative, which implies that the positivity conditions required for the existence of the anomalous cut are not fulfilled anymore. The trajectory of s + during the Wick rotation is shown in Fig. 4. There one sees that in all possible cases s + moves through the unitarity cut back into the second sheet, so that indeed, the anomalous cut disappears. The two figures show the two possible cases, where the anomalous threshold enters the physical sheet. Fig. 4a shows the trajectory of the anomalous threshold for a fixed value of 0 < q 2 1 < 4m 2 and q 2 2 varying from −∞ to ∞ (solid line, light gray indicates the part of the trajectory on the second sheet, whereas the part on the physical sheet is shown in black). Fig. 4b shows the analogous trajectory for a fixed value of q 2 1 > 4m 2 . The dashed lines are intermediate trajectories, when the two momenta are simultaneously Wick rotated by the same phase. The dotted lines indicate the movement of selected points on the trajectory of the anomalous threshold during the Wick rotation. All points of the trajectory that originally lie on the physical sheet move back through the unitarity cut into the second sheet.
As a further check we have calculated numerically a simplified toy example that exhibits all relevant features discussed above. We consider where The function defined by the last two integrals is proportional to the Feynman-parameter representation of the triangle loop function C 0 (s; q 2 1 , q 2 2 ), while the remainder of a toy µ is chosen in as close analogy to the real HLbL integral as possible, including a sufficient number of propagators to make the integral converge. In the standard loop-integral approach the intricacies due to anomalous thresholds are automatically taken into account by the i -prescription in (4.33) and (4.34), leading to the following Feynman-parameter representation In analogy to the HLbL master formula (4.28) the Wick-rotated representation becomes (the Q i have been rescaled by m µ ) We checked numerically that indeed the representations (4.35) and (4.36) are identical. This confirms that the Wick rotation is permitted despite the occurrence of anomalous thresholds, all signs of which disappear in the Feynman parameterization as long as the i -prescription is applied consistently. Similarly, the final result for the Wick rotation displays no remnants of the anomalous thresholds. As shown in Fig. 4, in this case the deeper reason can be traced back to the trajectory of the anomalous threshold during the Wick rotation towards the second sheet. Finally, we note that the above discussion of the anomalous threshold of the triangle diagram is already sufficient for the full HLbL calculation. Box diagrams do not lead to further complications, since, in the master formula (4.23), the limit k → 0 is taken before the Wick rotation of the loop momenta, hence the analytic structure of box diagrams is already reduced to the one of triangle diagrams.

Mandelstam representation
In the previous section, we have obtained a master formula (4.28) for the HLbL contribution to the anomalous magnetic moment of the muon, where the hadronic dynamics is parametrized in terms of the scalar functions Π i . Since these functions are free of kinematic singularities and zeros, they are the quantities that should satisfy a Mandelstam representation [61]. We need to determine seven scalar functions that are not related to each other by crossing symmetry. Due to the complexity of the problem, we cannot obtain an exact solution for the scalar functions but have to rely on approximations. Our strategy is to order the contributions to the scalar functions by the mass of the intermediate states. The lightest intermediate states are expected to give the most important contribution, heavier states are suppressed by the higher threshold and smaller phase space. In this paper, we will consider the two lowest-lying contributions: one-and two-pion intermediate states in all channels, i.e. pion-pole and pion-box topologies. These two contributions can be described without any further approximation. Of course, pion-pole and pion-loop contributions have been considered before in various model calculations. The important point of the dispersive approach is that these contributions become unambiguously defined: the pion pole, with on-shell pion transition form factor, corresponds precisely to one-pion intermediate states, and the sQED pion loop dressed with pion vector form factors to two-pion intermediate states with pionpole LHC. The explicit proof of these identifications will be presented in this section, based on a Mandelstam representation for the BTT scalar functions.

Derivation of the double-spectral representation
For the derivation of a Mandelstam representation of the scalar functions, we follow the discussion in [62]. We assume that the photon virtualities q 2 i are fixed and small enough so that no anomalous thresholds are present. A parameter-free description of the HLbL tensor and therefore an unsubtracted dispersion relation is crucial. The behavior of the imaginary parts, which is determined by the asymptotics of the sub-processes, does indeed suggest that no subtractions are needed. Furthermore, even a quark-loop contribution to the HLbL tensor has an asymptotic behavior that requires no subtractions. 13 Hence, for a generic scalar function Π i , we write a fixed-t dispersion relation without any subtractions:  where c t i is supposed to behave as lim t→0 c t i = 0 and takes into account the t-channel pole. The imaginary parts are understood to be evaluated just above the corresponding cut. The primed variables fulfill If we continue the fixed-t dispersion relation analytically in t, we have to replace the imaginary parts by the discontinuities, defined by Both the discontinuities as well as the pole residues are determined by s-or u-channel unitarity, which also defines their analytic continuation in t. While ρ t i;s,u are due to a one-pion intermediate state, D t i;s,u are due to multi-particle intermediate states, see First, we study the pion-pole contribution by analyzing the unitarity relation: where S n is the symmetry factor of the intermediate state |n . We consider now only the π 0 intermediate state in the sum: After reducing the matrix elements and using the definition of the pion transition form factor we find Im π s Π µνλσ = − 1 2 dp (2π) 4 δ (4) (q 1 + q 2 − p) µναβ λσγδ q 1α q 2β q 3γ q 4δ × F π 0 γ * γ * q 2 1 , q 2 2 F π 0 γ * γ * q 2 3 , q 2 4 = −πδ(s − M 2 π ) µναβ λσγδ q 1α q 2β q 3γ q 4δ F π 0 γ * γ * q 2 1 , q 2 2 F π 0 γ * γ * q 2 3 , q 2 4 . (5.8) By projecting onto the scalar functions Π i , this leads to and, analogously, In order to identify the discontinuities, we project the unitarity relation selecting two-pion intermediate states: Im ππ s e 4 (2π) 4 δ (4) (q 1 + q 2 + q 3 − q 4 )H λ1λ2,λ3λ4 = 1 2 dp 1 dp 2 π + (p 1 )π − (p 2 )|γ * (−q 3 , λ 3 )γ * (q 4 , λ 4 ) * π + (p 1 )π − (p 2 )|γ * (q 1 , λ 1 )γ * (q 2 , λ 2 ) + 1 4 dp 1 dp 2 π 0 (p 1 )π 0 (p 2 )|γ * (−q 3 , λ 3 )γ * (q 4 , λ 4 ) * π 0 (p 1 )π 0 (p 2 )|γ * (q 1 , λ 1 )γ * (q 2 , λ 2 ) , where the subscripts {+−, 00} denote the pion charges. The analytic continuation of the unitarity relation can be obtained if the γ * γ * → ππ matrix element W µν is expressed in terms of the fixed-s dispersion relation (2.42) for its scalar functions: Note that W µν 00 does not contain any pole terms because the photon does not couple to two neutral pions due to angular momentum conservation and Bose symmetry.
If we pick the contribution of the pole terms on both sides of the cut, we single out box topologies: (5.14) where the primed variables belong to the sub-process on the left-hand side and the double-primed variables to the sub-process on the right-hand side of the cut.
We could now apply a tensor reduction to obtain phase-space integrals for the reduced scalar quantities. The projection on the scalar functions Π i would allow us to identify the discontinuities D t i;s due to box structures. The reduced scalar integrals could then be transformed into another dispersive integral. Together with the dispersion integral ds of the primary cut, this produces the double-spectral representation. The case of the simplest scalar phase-space integral is explained in [33].
Unfortunately, the tensor reduction is enormously complex because it involves the inversion of a 138 × 138 square matrix. We will therefore follow a different strategy for the box topologies. We note that the non-zero pole residuesρ s;+− i;t,u contain two electromagnetic pion form factors for the off-shell photons. These form factors can be factored out and multiply then the discontinuity that would be obtained by applying Cutkosky's rules [63] to the sQED pion loop calculation. This becomes clear from the relation between the pole terms and the sQED Born terms as discussed Sect. 2.6. Therefore, the box contribution is nothing else but the sQED contribution multiplied by a vector form factor F V π (q 2 i ) for each of the off-shell photons. As in the case of the sub-process, the difference between unitarity diagrams and Feynman diagrams is absolutely crucial: the sQED contribution consists of boxes, triangles, and bulb Feynman diagrams, but corresponds to the pure box topology in terms of unitarity. In Sect. 5.3, we will prove that this identification is correct and unique.
Finally, there are the contributions with discontinuities either in one or in both of the sub-processes. They contain for example resonance contributions in the sub-process, but also two-pion rescattering effects in the direct channel.

Symmetrization and classification into topologies
In the previous subsection, we have explained how the double-spectral representation can be derived from a fixed-t dispersion relation by taking the analytic continuation in t, which is defined by the unitarity relation. In the s-channel unitarity relation, a fixed-s dispersion relation of the sub-process is inserted (in the unitarity relation for the u-channel contribution, the variable u is kept fixed, which, however, plays again the role of s in the sub-process). Of course, one could have started with a fixed-u or fixed-s dispersion relation in the first place. The requirement that this lead to the same result allows us to identify a symmetric representation, which treats the Mandelstam variables on an equal footing and therefore implements crossing symmetry. In this symmetric representation, we classify the different contributions in terms of topologies. Note that in the case of HLbL, we get the two other possibilities (i.e. taking fixed-u and fixed-s dispersion relations as the starting point) for free, because we consider a totally crossing symmetric process.
If we compare the three different representations of unsubtracted dispersion relations (5.4), we immediately see that the contributions in one representation are either contained explicitly in the other representations, can be understood as part of the respective constant c t i , or correspond to a contribution of neglected higher intermediate states, as we are now going to discuss.

Pion-pole contribution
First, we consider the pion-pole topology (see Fig. 6). The fixed-t dispersion relation contains the poles in the sand u-channel explicitly. Analogously, the fixed-s (fixed-u) dispersion relation contains the poles in the t-and u-channel (s-and t-channel). The t-channel pole contribution, which is not explicit in the fixed-t representation, can be identified with c t i as it vanishes in the limit t → ∞. Hence, the total pion-pole contribution is just given by where the pole residues are products of pion transition form factors: where δ ij is the Kronecker delta.
Since the Lorentz structures of the usual pion-pole contribution, see e.g. [19], coincide with the first three BTT structures, this proves that in a dispersive approach the pion pole is unambiguously defined as given in [26] (and already in [19]), in particular with the transition form factor defined for an on-shell pion.

Box contribution
For the contribution of box topologies (see Fig. 7), Mandelstam diagrams prove very useful for the discussion of double-spectral regions. In Fig. 8, such a diagram is shown for the case q 2 i = 0.5M 2 π . A dashed line indicates a line of fixed t, used for writing the fixed-t dispersion relation. The two cuts are highlighted in gray. The discontinuities along these cuts can be written again as a dispersive integral over double-spectral functions. The three regions of non-vanishing double-spectral functions are labeled in Fig. 8 by ρ st , ρ su , and ρ tu . The s-channel cut receives contributions from the double-spectral regions ρ st and ρ su , according to the unitarity diagrams 7a and 7b, where first the vertical cut, then the horizontal cut is applied. 14 The u-channel cut receives contributions again from ρ su and from ρ tu , according to the unitarity diagrams 7b and 7c. In diagram 7b, the horizontal cut is now applied first.
Hence, the fixed-t dispersion relation leads to a priori four double-spectral integrals: one for each of the regions ρ st and ρ tu and two for the region ρ su . However, it turns out that the sum of the two double-spectral integrals for the region ρ su equals the crossed version of one of the other double-spectral integrals. This is illustrated for the example of a simple scalar box diagram in App. G.3. Therefore, the box contributions constructed from a fixed-t dispersion relation are already crossing symmetric and identical to the box contributions that are obtained from a fixed-s or fixed-u dispersion relation.
This discussion already anticipates the idea of the uniqueness proof in Sect. 5.3: isolating the pion-pole contribution to the LHC in ππ intermediate states, the double-spectral functions become identical to the ones of the box diagrams in the sQED loop calculation. What needs to be shown is that the Lorentz structures in the BTT set not only produce the right Lorentz structures for the sQED box diagrams, but also that the triangleand bulb-type sQED diagrams are correctly reproduced. Moreover, this identification has to be proven for the full set of BTT functions, including those which cannot immediately be identified with basis functions due to the occurrence of the Tarrach amplitudes.

Higher intermediate states
Topologies where we replace either one or both of the pion poles in the sub-process by a multi-particle cut are the last contribution with a two-pion intermediate state in the direct channel. Such contributions could be captured in several ways, either relying on a partial-wave picture or approximating the multi-pion states by an effective resonance description. While the case of S-waves was already discussed in [26], the BTT formulation facilitates the generalization to D-waves, since the technical complications related to angular averages and kinematic singularities, the latter requiring the introduction of off-diagonal kernels, are avoided. The application to rescattering effects will be presented in a forthcoming publication.

Uniqueness of the box topologies
In this subsection, we will prove that the box topologies are equal to the FsQED contribution (sQED multiplied by pion vector form factors) and that this contribution is unique. The proof is based on the Mandelstam representation and consists of two arguments: • the double-spectral densities of box topologies and FsQED agree, • both representations fulfill an unsubtracted Mandelstam representation.
The argument for the first point is the following: in the case of the box topologies, the discontinuities are defined by the unitarity relation as phase-space integrals of products of two pure pole contributions for γ * γ * → ππ. In the case of FsQED, the discontinuities are most easily found by applying Cutkosky's rules [63], i.e. by cutting the loop diagrams. We see immediately that the discontinuities in FsQED are phase-space integrals of products of two Born terms for γ * γ * → ππ. In Sect. 2.6, we have already shown that the Born contribution is the same as the pure pole, hence the discontinuities (and therefore also the double-spectral densities) of the box topologies and FsQED are equal. To complete the proof it remains to be shown that the scalar functions in FsQED fulfill a Mandelstam representation (the scalar functions of the box topologies are defined by the Mandelstam representation). Apart from the first six scalar functions Π 1 , . . . , Π 6 , which are unaffected by the Tarrach amplitudes, the scalar BTT functions Π i are defined only up to the ambiguity (3.22), so that we cannot compare the representations for these functions immediately. Therefore, we will derive the double-spectral representation of the basis functions Π i that is implied by an unsubtracted Mandelstam representation for the scalar functions Π i and show that the FsQED basis functions fulfill a double-spectral representation of this form. This will then complete the proof of the uniqueness and the equality of the box topologies and the FsQED contribution.
First, we calculate the fully off-shell sQED loop contribution. It consists of six box diagrams, twelve triangles, and three bulb diagrams: The three bulb diagrams are given by where the scalar loop function B 0 is defined in (G.1). One of the twelve triangle diagrams can be represented as where the scalar loop function C 0 is defined in (G.4). For the tensor coefficient functions C i , C ij , we use the convention of [64] up to a normalization factor: etc. The remaining eleven triangles follow by crossing the external momenta. Similarly, the fully off-shell formula for the box diagrams can also be expressed in terms of the scalar four-point function D 0 and the tensor coefficient functions D i , D ij , D ijk , D ijkl , but the result is rather long and need not be reproduced here.
In this way, we obtain a representation in terms of scalar loop functions B 0 , C 0 , D 0 , and the pertinent tensor coefficient functions. Next, we project this expression onto the HLbL tensor basis (3.18), see App. C, in order to identify the sQED contribution to the basis functionsΠ i . Using FeynCalc [65], we perform a Passarino-Veltman reduction [66,67] of the tensor coefficient functions, so that the basis functions are given by linear combinations of the scalar loop functions: i D 0 (q 2 1 , q 2 2 , q 2 4 , q 2 3 , s, t, M 2 π , M 2 π , M 2 π , M 2 π ) + d su i D 0 (q 2 1 , q 2 2 , q 2 3 , q 2 4 , s, u, M 2 π , M 2 π , M 2 π , M 2 π ) + d tu i D 0 (q 2 1 , q 2 3 , q 2 2 , q 2 4 , t, u, M 2 π , M 2 π , M 2 π , M 2 π ), (5.20) where the coefficients p i , a i , b j i , c j i , d j i are meromorphic functions of the virtualities q 2 i and the Mandelstam variables s, t, u. The Passarino-Veltman reduction results unfortunately in extremely long expressions for the coefficients, which, in addition, contain many kinematic singularities. However, at the kinematic points of these singularities the scalar loop functions fulfill linear relations which cancel most of the singularities. By evaluating the expressions for the functionsΠ sQED i numerically, we have checked that the only surviving kinematic singularities are exactly the ones required by the projection of the 54 scalar functions Π i onto the basis (B.2).
As already mentioned, the verification of the Mandelstam representation is simplest for the functions Π 1 , . . . , Π 6 , because they are identical to the basis functionsΠ 1 , . . . ,Π 6 , i.e. their projection does not involve kinematic singularities. We have to show that these functions fulfill a Mandelstam representation of the form (suppressing the i ) , (5.21) where the borders of the double-spectral regions t + , u + are defined in (G.12).
As the scalar four-point functions D 0 are the only ones in (5.20) with a double-spectral region, we can immediately identify the double-spectral densities: where ρ 0 is the double-spectral density of the D 0 scalar loop function, defined in (G.15). By evaluating numerically the double-spectral integrals (5.21) at some random kinematic points below the appearance of anomalous thresholds, we have checked that the Mandelstam representation for the functions Π sQED This result implies that the double-spectral densities of the sQED box diagrams (5.22) have the correct form that, when inserted into the Mandelstam representation, the non-residue pieces can be separated in such a way that all triangle and bulb loop functions in (5.20) are reproduced with the correct coefficients. The explicit verification of this property is the crucial point in the argument. Since the double-spectral densities of box topologies and FsQED agree, this proves the claim for i = 1, . . . , 6.
The remaining basis functions exhibit kinematic singularities. In order to work out their double-dispersive representation, we can limit ourselves to the seven remaining representatives of (3.19). All other basis functions are related by crossing symmetry (B.3).
We write the basis functions of interest as where s a := q 2 1 + q 2 2 , s b := q 2 3 + q 2 4 , u a := q 2 2 + q 2 3 , u b := q 2 1 + q 2 4 , t b := q 2 2 + q 2 4 . Next, we introduce subtractions for the scalar functions Π i : each function is subtracted once or twice at the subtraction point that appears in the coefficient of the other function Π i in the above linear combinations. For instance, in the equation forΠ 9 , we subtract Π 9 at u = u b and Π 22 at s = s b . This leads to representations of the basis functionsΠ i where both double-spectral contributions are multiplied by a common prefactor (e.g. (s − s b )(u − u b ) in the case ofΠ 9 ), while only the subtraction terms, which are single-dispersion integrals, have different prefactors. The results for these double-spectral representations of the basis functionsΠ i are given in (F.2) and (F.5) in App. F.1.
The important point is that in these double-spectral representations both the discontinuities of the singledispersion integrals and the combined double-spectral densitiesρ i are unambiguously defined. We do not attempt to split the double-spectral densities into two contributions from the Π i : this splitting involves exactly the ambiguity (3.22). For the sQED contribution, the discontinuities and double-spectral densities can be extracted from the loop representation of the basis functions (5.20) as follows: • The discontinuities are extracted by taking the appropriate limit of the basis functionsΠ i . For instance, the discontinuity of the subtraction term of Π 9 is obtained by considering the limit lim • The combined double-spectral densitiesρ i are (up to a kinematic factor) just the double-spectral densities of the basis functionsΠ i . Therefore, they can be obtained again from the coefficients of the D 0 functions in the loop representation, in analogy to the case ofΠ 1 , . . . ,Π 6 .
After having identified all the discontinuities and double-spectral densities in (F.2) and (F.5), we have checked numerically for random kinematic points (below the appearance of anomalous thresholds) that the dispersive representations of the functionsΠ i agrees with the loop representation. It turns out thatΠ sQED which also fixes the redundancy (3.23). This completes our proof of the uniqueness of the pion-box contribution. Cutkosky's rules tell us that the discontinuities of the FsQED contribution are the same as the ones of the pion-box topologies in the sense of unitarity. The FsQED contribution fulfills the same double-spectral representation as the pure pion-box topologies. Therefore, the two representations are the same. Unitarity and Mandelstam analyticity define the pion-box contribution, i.e. ππ intermediate states with a pion-pole LHC, in a unique way.
Finally, we stress that these calculations also provide a strong test of the Lorentz decomposition (3.13). Apart from the function Π 49 , which does not receive a contribution from the pion loop, all scalar functions have been shown to be free of kinematics, so that the sQED pion-loop amplitude behaves exactly as expected from the general BTT formalism.

Contribution to (g − 2) µ
In this subsection, we insert our dispersive representation of the scalar functions into the master formula (4.28) to obtain the contribution to a µ .

Pion-pole contribution
With (5.15) and using the master formula (4.28), we recover the well-known result for the pion-pole contribution to a µ [19]: where Q 2 3 = Q 2 1 + 2Q 1 Q 2 τ + Q 2 2 and the integral kernels T i are given in App. E.2.

Pion-box contribution
The single-integral discontinuities and the double-spectral densities in the dispersive representations of the basis functions (F.2) and (F.5) are quantities that can be extracted directly from the projected basis functionsΠ i . To the contrary, the separation of the double-spectral densitiesρ i into the two contributions from the different scalar functions Π i is not unambiguously possible, which reflects just the redundancy (3.22). However, such a separation is not necessary: for the calculation of a µ , we need the scalar functions Π i only in the limit k → 0. In this limit, all the scalar functions Π i appearing in the master formula (D.2) can be expressed in terms of singledispersion integrals, where the discontinuities are directly related to the basis functionsΠ i . All the subtracted double-spectral integrals, which are not unambiguously defined, drop out in the limit k → 0.
We stress that the ambiguities in the definition of the BTT functions Π i are related to the fact that these functions are not observables. In the calculation of any physical quantity, such as helicity amplitudes or the HLbL contribution to a µ , any ambiguity due to the redundancy in the BTT set has to drop out.
The pion-box contribution to a µ is therefore given by

Conclusion and outlook
We have derived a decomposition of the HLbL tensor into scalar functions (3.13), following the general recipe by Bardeen, Tung, and Tarrach, which allowed us to derive a master formula (4.28) for the HLbL contribution to the anomalous magnetic moment of the muon (g − 2) µ in terms of scalar functions that by construction are free of kinematic singularities and zeros. In particular, these scalar functions fulfill a Mandelstam double-spectral representation. In the present work, we have considered only the lowest-lying intermediate states, pion-pole and pion-box topologies. These two contributions can be studied without any approximation. The central result is that the pion pole corresponds exactly to the contribution of one-pion intermediate states if interpreted with an on-shell pion transition form factor, and the pion-box, i.e. ππ intermediate states with a pion-pole LHC, to the sQED pion loop dressed with pion vector form factors. These results confirm both the definition of the pion pole and prove explicitly the validity of the separation of Born and non-Born ππ contributions given in [26]. Moreover, while the formulation in terms of helicity amplitudes adopted in [26] required the introduction of off-diagonal kernel functions to account for kinematic singularities, the BTT version avoids these complications and thus facilitates the generalization to partial waves beyond S-waves. In order to exploit unitarity in a partial-wave picture, the partial-wave unitarity relations analyzed in the helicity basis need to be transformed to the BTT basis, the formalism for which we provided in the present paper.
This treatment is based on fundamental principles of particle physics: gauge invariance and crossing symmetry are already implemented in the decomposition of the HLbL tensor into scalar functions. The dispersive description uses analyticity and unitarity to establish a relation between the HLbL contribution and different on-shell quantities. These on-shell quantities are in principle either experimentally accessible or can be reconstructed from data with dispersive methods. For the two cases covered in this article, the corresponding input is parametrized by the pion transition form factor F π 0 γ * γ * in case of the pion pole, and by the pion vector form factor F V π for the pion box, both form factors being needed for negative virtualities of the off-shell photons. For general ππ intermediate states also information on the partial waves for γ * γ * → ππ is required. An overview over which processes can help in the dispersive reconstruction of the pion transition form factor [68][69][70][71] or γ * γ * → ππ [26,36,40,56,72] is provided in [27].
An extension of the presented dispersive treatment is possible within certain limits. It is straightforward to include higher pseudoscalar poles, i.e. the η and η mesons by just adding their contribution in complete analogy to the π 0 pole. The input quantities will be the transition form factors of these heavier pseudoscalars, see [73][74][75] for work towards their dispersive calculation. Although these mesons are unstable in QCD, their decay width is certainly small enough to justify the treatment as a pure pole. A bit more difficult is the inclusion of higher two-particle intermediate states: an extension to KK intermediate states is still straightforward, but the experimental input will be less accurate. In general, the approach applies to the low and intermediate energies 1.5 GeV, where a few channels dominate. All model calculations done so far indeed support the assumption that the lowest-lying singularities govern the HLbL tensor.
A careful numerical analysis of the formalism presented in this paper is in progress. It will reveal the relative importance of the different contributions and should allow us to identify which input quantities will have the largest impact concerning the reduction of the hadronic uncertainty in a HLbL µ . In this way, we are confident that the presented treatment of HLbL scattering shows a path towards a data-driven and thus less model-dependent evaluation of (g − 2) µ .
A Partial-wave representation for γ * γ * → ππ Here, we give the expressions for the scalar basis functions of the sub-process γ * γ * → ππ, defined in (2.26), in terms of the helicity partial waves including D-waves: We note that under q 2 1 ↔ q 2 2 , h 2,4 changes sign, while all the other S-and D-waves are invariant. Therefore, the apparent singularities in

B.2 Basis coefficient functions
The relation between the 43 basis coefficient functionsΠ i and the 54 scalar functions of the redundant set can be obtained by projecting (3.13) on the basis. The result defines the exact form of the kinematic singularities in q 1 · q 2 and q 3 · q 4 in the functionsΠ i : Π 23 = Π 25 − q 1 · q 3 q 2 · q 4 q 1 · q 2 q 3 · q 4 Π 26 ,Π 24 = Π 27 − q 1 · q 4 q 2 · q 3 q 1 · q 2 q 3 · q 4 Π 28 , Crossing symmetry for the functions Π i implies the following relations between the basis functionsΠ i :

C Projection of the scalar functions
Given any representation of the HLbL tensor Π µνλσ , the following procedure allows the identification of the basis coefficientsΠ i in (3.18): • write the tensor in terms of the 138 elementary structures (3.8) and identify the scalar coefficients Ξ i , • take the subset consisting of the following 43 scalar coefficients: • perform a basis change according toΠ where P is an invertible 43 × 43 matrix, defined below.

E.2 Kernels for the master formula
In the master formula (4.28), the following integral kernels appear: − q 1 · q 2 q 3 · q 4 q 1 · q 4 q 2 · q 3Π 21 , Π 31 (Q 1 , Q 2 , τ ) = lim k→0 lim s→s b q 3 · q 4 q 1 · q 4 q 2 · q 3Π 19 , Π 33 (Q 1 , Q 2 , τ ) = lim k→0 C 23 lim s→s b q 3 · q 4 q 1 · q 4 q 2 · q 3Π 19 , Π 39 (Q 1 , Q 2 , τ ) = lim k→0 C 14 lim s→sa q 1 · q 2 q 1 · q 4Π 36 , Π 40 (Q 1 , Q 2 , τ ) = lim k→0 C 12 C 14 lim s→sa q 1 · q 2 q 1 · q 4Π 36 , Π 42 (Q 1 , Q 2 , τ ) = lim k→0 C 12 C 24 lim s→sa q 1 · q 2 q 1 · q 4Π 36 , Π 46 (Q 1 , Q 2 , τ ) = lim k→0 C 14 C 23 lim s→sa q 1 · q 2 q 1 · q 4Π 36 , factors of (t − u) multiply the double-spectral contribution: They have exactly the form of the difference terms in (F.20), hence the difference between the choices of Π 49 or S 49 /Â 49 as the functions to fulfill an unsubtracted Mandelstam representation can be reabsorbed in a shift of D t;s 49 and D u;s 49 . This means that the basis functionsΠ i and the HLbL tensor remain unaffected by the kinematic zero in Π 49 . In particular, the form of the dispersion relation for (g − 2) µ in (F.9) is unchanged: the functions Π 50 , Π 51 , and Π 54 are defined by crossed versions of a fixed-s dispersion relation for Π 49 , where the effect of the kinematic zero is reabsorbed in the discontinuities of the t-and u-channel integrals. Since in the end the discontinuities will be calculated directly from the unitarity relation, the observation that the kinematic zero in Π 49 does not alter the functional form of the dispersion relation completes the proof that there are no practical consequences, neither in the construction of the dispersion relation nor in (g − 2) µ . The scalar two-point function [64] corresponding to the loop diagram in Fig. 9 is given by

G Dispersive representation of loop functions
. (G.1) We define s := p 2 . The two-point function has a normal threshold at s = (m 1 + m 2 ) 2 . According to Cutkosky's rules [63], the discontinuity is given by: Consider the scalar three-point function [64] corresponding to the loop diagram in Fig. 10: . (G.4)