Dispersion relation for hadronic light-by-light scattering: two-pion contributions

In this third paper of a series dedicated to a dispersive treatment of the hadronic light-by-light (HLbL) tensor, we derive a partial-wave formulation for two-pion intermediate states in the HLbL contribution to the anomalous magnetic moment of the muon $(g-2)_\mu$, including a detailed discussion of the unitarity relation for arbitrary partial waves. We show that obtaining a final expression free from unphysical helicity partial waves is a subtle issue, which we thoroughly clarify. As a by-product, we obtain a set of sum rules that could be used to constrain future calculations of $\gamma^*\gamma^*\to\pi\pi$. We validate the formalism extensively using the pion-box contribution, defined by two-pion intermediate states with a pion-pole left-hand cut, and demonstrate how the full known result is reproduced when resumming the partial waves. Using dispersive fits to high-statistics data for the pion vector form factor, we provide an evaluation of the full pion box, $a_\mu^{\pi\text{-box}}=-15.9(2)\times 10^{-11}$. As an application of the partial-wave formalism, we present a first calculation of $\pi\pi$-rescattering effects in HLbL scattering, with $\gamma^*\gamma^*\to\pi\pi$ helicity partial waves constructed dispersively using $\pi\pi$ phase shifts derived from the inverse-amplitude method. In this way, the isospin-$0$ part of our calculation can be interpreted as the contribution of the $f_0(500)$ to HLbL scattering in $(g-2)_\mu$. We argue that the contribution due to charged-pion rescattering implements corrections related to the corresponding pion polarizability and show that these are moderate. Our final result for the sum of pion-box contribution and its $S$-wave rescattering corrections reads $a_\mu^{\pi\text{-box}} + a_{\mu,J=0}^{\pi\pi,\pi\text{-pole LHC}}=-24(1)\times 10^{-11}$.


Introduction
The long-standing discrepancy between the standard-model determination and the experimental measurement [1] (updated to the latest muon-proton magnetic moment ratio [2]) a exp µ = 116 592 089(63) × 10 −11 (1.1) of the anomalous magnetic moment of the muon (g−2) µ has triggered substantial interest in the subject on both the theoretical and the experimental side. The ongoing E989 experiment at Fermilab [3] as well as complementary efforts by J-PARC E34 [4] aim at improving the precision by a factor of 4, see [5] for a detailed account of the experimental strategies in both cases. On the theory side, the uncertainty is dominated by hadronic effects [6][7][8], while QED [9] and electroweak [10] contributions are under control at the level of at least 1 × 10 −11 . Currently, the dominant source of hadronic uncertainties is hadronic vacuum polarization (HVP) at O(α 2 ) in the fine-structure constant, closely followed by the O(α 3 ) hadronic light-by-light (HLbL) contribution, depicted in Fig. 1, and with higher-order insertions of the same hadronic amplitudes already under sufficient control [11][12][13][14]. In view of improved data input for the dispersion relation for HVP [15], it is likely that the stumbling block will eventually become the sub-leading HLbL contribution. Current estimates for HLbL scattering in (g − 2) µ are largely based on hadronic models [16][17][18][19][20][21][22][23][24][25][26][27], which despite implementing different limits of QCD, such as large-N c , chiral symmetry, or constraints from perturbative QCD, all involve a certain amount of uncontrollable uncertainties without offering a systematic path forward. In order to improve the determination of the HLbL contribution, we proposed a dispersive framework [28], based on the fundamental principles of analyticity, unitarity, gauge invariance, and crossing symmetry, which opens up a path towards a data-driven evaluation [29]. As the next step [30,31], we presented a comprehensive solution to the task of constructing a basis for the HLbL tensor devoid of kinematic singularities, defining scalar functions that are amenable to a dispersive treatment. In particular, we derived a Lorentz decomposition of the HLbL tensor that manifestly implements crossing symmetry and gauge invariance, with scalar coefficient functions free of kinematic singularities and zeros that fulfill the Mandelstam double-spectral representation. In this framework, we worked out how to define unambiguously and in a model-independent way both the pion-pole and the pion-box contribution. 1 With pion-as well as η-, η -pole contributions determined by their doubly-virtual transition form factors, which by themselves are strongly constrained by unitarity, analyticity, and perturbative QCD in combination with experimental data [38][39][40][41][42][43][44][45][46], we here apply our framework to extend the partial-wave formulation of two-pion rescattering effects for S-waves [28] to arbitrary partial waves. To this end, we identify a special set of (unambiguously defined) scalar functions that fulfill unsubtracted dispersion relations and can be expressed as linear combinations of helicity amplitudes. Their imaginary part, the input required in the dispersion relations, is provided in terms of helicity partial waves for γ * γ * → ππ by means of unitarity. Working out explicitly the basis change to the helicity amplitudes, we generalize the unitarity relation derived in [28] up to D-waves only to arbitrary partial waves. We demonstrate that indeed the summation of the partial waves reproduces the known full result for the pion box, to which the ππ-rescattering contribution is expected to produce the dominant correction. We provide the details of a first numerical analysis [47] of these rescattering effects based on helicity partial waves for γ * γ * → ππ that we construct dispersively from a pion-pole left-hand cut (LHC) and ππ phase shifts from the inverse-amplitude method, an approach that isolates pure ππ contributions and thus, in the isospin-0 channel, provides an estimate for the impact of the f 0 (500) resonance on HLbL scattering. In the same way, our γ * γ * → ππ amplitudes reproduce the phenomenological value for the chargedpion polarizability, thereby clarifying the role of the associated corrections in (g − 2) µ [48][49][50]. These results lay the groundwork for a future global analysis of two-meson intermediate states in the HLbL contribution.
The outline is as follows: Sect. 2 is devoted to a thorough derivation of partial-wave dispersion relations for the HLbL tensor, with tensor decomposition, dispersion relations, sum rules, and partialwave expansion addressed in Sects. 2.1-2.5. A short summary of the strategy is provided at the beginning of Sect. 2, complemented by a summary of the most important results in Sect. 2.6. In Sect. 3, a numerical evaluation of the pion box is provided based on fits of the pion vector form factor to high-statistics time-like and space-like data. The pion box is further used to explicitly verify the general results derived in Sect. 2, in particular to demonstrate the convergence of the partial-wave expansion for its contribution to (g − 2) µ . Rescattering corrections to the pion box are discussed in Sect. 4, including a numerical analysis of the S-wave contribution, before we conclude in Sect. 5. Further details of the formalism are provided in the Appendices.

Helicity formalism for HLbL
In this section, we derive the formalism for the evaluation of the HLbL two-pion contribution to (g − 2) µ . The goal of our treatment is to relate this contribution to helicity partial waves for the subprocess γ * γ * → ππ, which in principle are measurable input quantities or at least can be reconstructed dispersively.
The outline of this derivation is illustrated as a flowchart in Fig. 2. The first step is the decomposition of the HLbL tensor into Lorentz structures and scalar functions that are free of kinematic singularities and zeros. We have solved this problem in [31] and recapitulate the results in Sect. 2.1. This representation, referred to as BTT tensor decomposition [51,52] in Fig. 2, allows us to write the HLbL contribution to (g − 2) µ in full generality as a master formula that involves only three integrals. This master formula (2.25) applies to any conceivable HLbL tensor, as long as it is consistent with general properties that should be fulfilled by any admissible HLbL amplitude: gauge invariance, crossing symmetry, and the principle of maximal analyticity [53], i.e. the principle that the scattering amplitude can be represented by a complex function that exhibits no further singularities except for those required by unitarity and crossing symmetry. Any such singularities are of dynamical origin, and thus have to be contained within the scalar functionsΠ i in the master formula. Phrased differently, if a given amplitude for the HLbL tensor cannot be expressed in the BTT basis, e.g. due to the appearance of kinematic singularities, this automatically implies that this amplitude is at odds with said general properties.
The dynamics of HLbL scattering is thus contained in the scalar functions, which are the objects that we describe dispersively. In [31], we have used the Mandelstam representation for the scalar functions to study the pion-box contribution. In Sect. 2.2, we extend the dispersive treatment and derive from the Mandelstam representation single-variable dispersion relations for general two-pion contributions. Combining these single-variable dispersion relations with unitarity constraints requires a basis change to helicity amplitudes, since the partial-wave unitarity relation becomes diagonal only for definite helicity amplitudes. However, this basis change is complicated by the appearance of redundancies in the representation which, together with the requirement that longitudinal polarizations for on-shell photons not contribute in the final HLbL representation, necessitates a more careful study of the BTT scalar functions and their relation to helicity amplitudes. The solution to this problem is the explicit derivation of a basis that removes all redundancies and apparent contributions from unphysical polarizations, which is presented in Sect. 2.4. As a by-product we obtain a set of physical sum rules to be fulfilled by the scalar functions and thereby the helicity amplitudes.
After the basis change to helicity amplitudes, we can then employ the unitarity relation to determine the imaginary parts in the dispersion integrals in terms of helicity amplitudes for γ * γ * → ππ. In particular, we perform a partial-wave expansion of the helicity amplitudes and generalize the S-wave result of [28] to arbitrary partial waves, which is the main result of Sect. 2.5. In performing this analysis the partial waves for γ * γ * → ππ are treated as known, given quantities, which unfortunately they are not. The lack of experimental information can be partly compensated by theory constraints, in particular by dispersion relations in the form of Roy-Steiner equations [54][55][56][57]. A simplified, S-wave variant of these will be solved in Sect. 4.
A summary of the main results is provided in Sect. 2.6, including a glossary of the notation for the scalar functions. The subtleties in the various basis changes unfortunately require the introduction of different sets of scalar functions, whose dimension, defining equation, and main properties are summarized in Table 1. 2.1 Tensor decomposition and master formula for (g − 2) µ In this subsection, we recapitulate the decomposition of the HLbL tensor into a sum of gauge-invariant Lorentz structures times scalar functions that are free of kinematic singularities. We slightly modify and improve the master formula presented in [30,31] in such a way that crossing symmetry between all three off-shell photons remains manifest. The dynamical input in the master formula is encoded in only six different scalar functions and their crossed versions.

BTT decomposition of the HLbL tensor
The HLbL tensor is defined as the hadronic Green's function of four electromagnetic currents in pure QCD: Π µνλσ (q 1 , q 2 , q 3 ) = −i d 4 x d 4 y d 4 z e −i(q 1 ·x+q 2 ·y+q 3 ·z) 0|T {j µ em (x)j ν em (y)j λ em (z)j σ em (0)}|0 , (2.1) where the electromagnetic current includes the three lightest quarks: j µ em :=qQγ µ q, q = (u, d, s) T , Q = diag Roy-Steiner equations for γ * γ * → ππ Sect. 4 helicity partial waves for γ * γ * → ππ: The hadronic contribution to the helicity amplitudes for (off-shell) photon-photon scattering is given by the contraction of the HLbL tensor with polarization vectors: where q 4 = q 1 + q 2 + q 3 . The usual Mandelstam variables fulfill the linear relation Gauge invariance requires the HLbL tensor to satisfy the Ward-Takahashi identities Based on a recipe by Bardeen, Tung [51], and Tarrach [52] (BTT), we have derived in [30,31] a decomposition of the HLbL tensor with tensor structures reproduced here for completeness (all remaining ones follow from crossing symmetry [31]) The BTT decomposition has the following properties: • all the Lorentz structures fulfill the Ward-Takahashi identities, i.e.
The first two properties make gauge invariance and crossing symmetry manifest, while the third property provides the foundation for writing dispersion relations: in a dispersive treatment, we exploit the analytic structure of the scalar functions dictated by unitarity and we have to make sure that the singularity structure due to the hadronic dynamics is not entangled with kinematic singularities. Since the number of helicity amplitudes for fully off-shell photon-photon scattering is 41, the set of 54 structures {T µνλσ i } does not form a basis, but exhibits a 13-fold redundancy, as we discussed in detail in [31]. While 11 linear relations hold in general, two additional ones are present in four space-time dimensions [58]. Away from four space-time dimensions, a subset of 43 Lorentz structures forms a basis: where the basis-coefficient functionsΠ i are no longer free of kinematic singularities. However, the explicit structure of their kinematic singularities follows from the projection of the BTT decomposition onto this "basis."

Master formula for the HLbL contribution to
Based on a projection technique in Dirac space, one can extract the HLbL contribution to a µ := (g − 2) µ /2 from the following expression: There are only 19 independent linear combinations of the structures T µνλσ i that contribute to (g − 2) µ , hence we can make a basis change in the 54 structures in such a way that in the limit q 4 → 0 the derivative of 35 structuresT µνλσ i vanishes. Since the loop integral and the propagators are symmetric under q 1 ↔ −q 2 , in [31] we made sure to preserve crossing symmetry under exchange of q 1 and q 2 , but did not yet exploit the fact that it is even possible to preserve crossing symmetry between all three off-shell photons-the limit q 4 → 0 singles out one of the photons, but the remaining three are completely equivalent. For the sake of simplifying further calculations, we present here new structuresT µνλσ i and the corresponding scalar functionsΠ i , superseding the ones given in [31].
The 19 structuresT µνλσ i contributing to (g − 2) µ can be chosen as follows: do not contribute to (g − 2) µ and are given in App. A. The set of 19 linear combinations of scalar functions that give a contribution to (g − 2) µ is defined by (replacing Eq. (D.1) in [31]) together with the crossed versions thereof Π 8 = C 12 Π 7 ,Π 9 = C 12 C 13 Π 7 ,Π 10 = C 23 Π 7 ,Π 13 = C 13 Π 7 ,Π 14 = C 12 C 23 Π 7 , where the crossing operators C ij exchange momenta and Lorentz indices of the photons i and j, e.g. 2 The following intrinsic crossing symmetries are preserved (we do not list the symmetries involving the fourth photon):Π where the dots denote three more crossing relations that follow from the given ones. Hence, the scalar functionsΠ i contributing to (g − 2) µ fall into only six distinct classes that are closed under crossing symmetry of the off-shell photons 1, 2, and 3. Apart fromΠ 39 , which is fully symmetric, the representatives in (2.15) are picked because they share a common property: their s-channel is special as follows from the observation that the corresponding Lorentz structuresT µνλσ i are (anti-)symmetric under either C 12 or C 34 (or both). This is reflected in the intrinsic crossing symmetries (2.18). 3 The HLbL contribution to (g − 2) µ can now be written as where G := {1, . . . , 11, 13, 14, 16, 17, 39, 50, 51, 54} and . (2.20) As in [31], we perform a Wick rotation, average the result over the direction of the Euclidean fourmomentum of the muon, and use the Gegenbauer polynomial technique [59] to perform five of the eight integrals in full generality, i.e. without prior knowledge of the functionsΠ i . The symmetry properties of the loop integral and the kernelsT i under q 1 ↔ −q 2 allow us to write the master formula for the HLbL contribution to (g − 2) µ containing a sum of only 12 terms: where Q 1 := |Q 1 | and Q 2 := |Q 2 | denote the norm of the Euclidean four-vectors. The 12 scalar functionsΠ i are a subset of the functionsΠ i : They have to be evaluated for the reduced (g − 2) µ kinematics Due to the basis change, the kernel functions T i differ slightly from the ones given in [31]. We provide the explicit expressions in App. B. µ . The border of the integration region is at r = 1 and corresponds to τ = −1 for π/3 < φ < 5π/3 (solid gray line) and τ = 1 otherwise (dashed gray line). The angles φ = π/3, φ = π, and φ = 5π/3 correspond to Q 2 2 = Q 2 3 , Q 2 1 = Q 2 2 , and Q 2 1 = Q 2 3 , respectively. The three points where one of the Q 2 i is zero are singularities of the integration kernels. The height of the equilateral triangle is given byΣ.
In [60] a different parametrization of the (g−2) µ integration region has been proposed, which proved advantageous for the numerical implementation. We perform the following variable transformation in the master formula (note thatΣ = −Σ is the sum of the squared Euclidean virtualities, whereas Σ denotes the sum of the squared Minkowskian virtualities): The range of integration is thenΣ ∈ [0, ∞), r ∈ [0, 1], and φ ∈ [0, 2π]. The integration region in the Mandelstam plane and the meaning of the variables is illustrated in Fig. 3. After the variable transformation, the master formula becomes where Q 1 , Q 2 , and τ are understood as functions ofΣ, r, and φ. The master formula for the HLbL contribution to (g − 2) µ is exact and completely general: given any representation of the HLbL tensor, one can project out the six scalar functionsΠ i in (2.15). Using these and their crossed versions, one can construct the 12 scalar functionsΠ i in (2.22), which encode the entire dynamical content of HLbL scattering relevant for (g − 2) µ . After their insertion into the master formula (2.25), only a three-dimensional integral has to be carried out.
In a next step, we aim at reconstructing the scalar functionsΠ i using dispersive methods, which will be the content of the remainder of this section.

Dispersion relations for the HLbL tensor
In this subsection, we discuss the dispersive framework that we employ for the reconstruction of the scalar functions. In order to calculate the two-pion contributions beyond the pion box, input on the sub-process γ * γ * → ππ is needed. This input will be in the form of helicity partial waves which, in principle, could be measured or, in the absence of data on the doubly-virtual process, have to be reconstructed dispersively [54][55][56][57]. The partial-wave expansion turns, however, the amplitude into a polynomial in the crossed-channel Mandelstam variables, i.e. the cut structure in the crossed channel due to heavier (e.g. multi-pion) intermediate states gets lost. Therefore, with γ * γ * → ππ helicity partial waves as input, one has to use a single-variable dispersion relation. We derive in Sect. 2.2.3 a suitable form for such a dispersion relation that follows from the Mandelstam representation.

Mandelstam representation for HLbL
In [31], we have used Mandelstam's double-spectral representation [61] for the BTT scalar functions Π i in order to split the HLbL contribution to (g − 2) µ into the following sum: a HLbL µ = a π 0 -pole µ + a π-box µ + a ππ µ + . . . (2.26) This sum directly reflects the sum over intermediate states in the unitarity relation in which, by definition, all intermediate states enter on-shell. While unitarity alone defines the imaginary parts, the real parts are obtained from the dispersion integrals. In short, this amounts to the following procedure: • Write down the unitarity relation for the HLbL tensor.
• In the sum over intermediate (on-shell) states, the one-pion state contributes as a δ-function to the imaginary part, which offsets the dispersion integral and defines the π 0 -pole contribution.
• The next-heavier intermediate state in the unitarity relation is a two-pion state. So far, we concentrate on one-and two-pion intermediate states, shown in Fig. 4.
• In the two-pion contribution, write down the crossed-channel unitarity relation for the sub-process γ * γ * → ππ. The one-pion contribution in this unitarity relation defines the π-pole contribution to γ * γ * → ππ. Separating this pole contribution corresponds to further splitting the two-pion contribution to HLbL into different box-type topologies, shown in Fig. 5.
• The two-pion phase-space integral in the HLbL unitarity relation can be converted into a second (crossed-channel) dispersion integral. This nontrivial but essential technical step is described in detail in App. D of [30]. • Finally, the symmetrization over the different channels produces the Mandelstam representation.
The double-spectral representation for the pion box has the following form: , (2.27) where the functions ρ π-box i denote the double-spectral densities, which have been derived (though not given explicitly) in [31]. The borders of the double-spectral regions t + and u + are defined in App. G.3 of [31].
In [31], we have explicitly shown that the Mandelstam representation for the pion box is mathematically equivalent to a scalar QED (sQED) one-loop calculation, multiplied by appropriate pion vector form factors for the off-shell photons. First, the form factors only depend on the virtualities {q 2 i } and can be pulled out of the double-dispersion integral. Second, triangle and bulb diagrams appear in the sQED calculation only in order to ensure gauge invariance: indeed when projected onto our gauge-invariant tensor structures, the analytic structure of sQED is the one of pure box topologies. In order to calculate the pion-box contribution numerically, it is convenient to rather use a Feynman parametrization instead of the dispersive representation. It turns out that in the limit of (g − 2) µ kinematics, the Feynman parametrization of the scalar functionsΠ i defined in (2.15) is very compact. Due to the limit q 4 → 0, only two-dimensional Feynman parameter integrals appear: where F V π is the electromagnetic pion vector form factor and the integrands I i can be found in App. C, written in a way that shows explicitly the absence of kinematic singularities.
The main goal of the present article is to describe two-pion contributions beyond the pion box, i.e. the topologies that involve a crossed-channel intermediate state heavier than one pion in one or both sub-processes.

Two-pion contributions beyond the pion box
Let us examine in more detail the form of the Mandelstam representation as sketched in the previous subsection. The starting point is a fixed-t dispersion relation with a discontinuity given by the two-pion contribution to the unitarity relation for the HLbL tensor: where W µν are the matrix elements for γ * γ * → ππ. The subscripts {+−, 00} denote the charges and p 1,2 the momenta of the intermediate pions. The phase-space factor is In order to analytically continue the unitarity relation, these matrix elements have to be expressed in terms of fixed-s dispersion relations for the scalar functions in a proper tensor decomposition, see [31]: 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: 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. This contribution was the subject of study in [31]. We consider now the contributions with discontinuities either in one or both of the sub-processes: Figure 6: Unitarity diagrams representing the "2disc"-box contributions that are (partially) accessible through a fixed-t dispersion relation.
If the order of phase-space and dispersive integrals are exchanged, the phase-space integrals can be performed by applying a tensor reduction to the quantities The reduced scalar phase-space integrals can then be transformed into another dispersive integral. Together with the dispersion integral ds of the primary cut, this produces a double-dispersion relation. The case of the simplest scalar phase-space integral is explained in [30]. Here, we do not try to calculate explicitly the tensor phase-space integrals, because we are interested just in the analytic structure of the "1disc" and "2disc" contributions, i.e. the boxes with heavier intermediate states in one or both of the sub-processes. In order to obtain the full double-spectral representation, one has to consider not only a fixed-t dispersion relation as a starting point but also the crossed versions, i.e. fixed-s and fixed-u dispersion relations. The symmetrization leads to the Mandelstam representation. For a more detailed discussion in the case of the pion box, see again [31]. We consider now the "1disc" and "2disc" contributions, where the pole in one or both of the sub-processes is replaced by a discontinuity. As the symmetrization procedure is identical in both cases, we only discuss the case of a discontinuity in both sub-processes. Fig. 6 shows the unitarity diagrams corresponding to the double-spectral representations that are generated if we start in our derivation from the fixed-t dispersion relation: the diagrams 6a and 6b generate a cut for s > 4M 2 π , which is the right-hand cut in the fixed-t dispersion relation. The diagrams 6c and 6d are responsible for the left-hand cut for u > 4M 2 π . In all cases the first cut is always the one through the two-pion intermediate state.
As discussed in [30,31], an (st)-box diagram can be represented either by a fixed-s, fixed-t, or fixedu dispersion relation: in the case of a fixed-t representation, there appears only one dispersion integral along the right-hand s-channel cut. Likewise, in a fixed-s representation, only one dispersion integral along the t-channel cut is present. In the case of a fixed-u representation, however, an (st)-box generates two integrals along both the s-and the t-channel cut. This particularity translates directly into the double-spectral representation: the (st)-box can be written as only one double-dispersion integral if one starts from a fixed-s or fixed-t representation. If one starts from the fixed-u representation, one obtains a sum of two double-dispersion integrals, see App. G.3 of [31].
Consider now the Mandelstam diagram in Fig. 7, which shows the double-spectral regions that we generate if we start from a fixed-t dispersion relation. Because we consider in the primary cut only two-pion intermediate states, not all the contributions from the displayed double-spectral regions are generated. We understand from the above discussion of the (st)-box that ρ st and ρ ut are complete, but that the contributions from ρ us and ρ su are not, because only one double-spectral integral for each of these contributions is obtained. However, two double-spectral integrals would be needed to generate the full contribution of these regions: one of the two integrals has a primary cut at the higher threshold 16M 2 π and is neglected in the fixed-t representation. Of course, two more double-spectral regions ρ ts and ρ tu , which correspond to crossed boxes, are completely missing in the fixed-t representation. The complete set of double-spectral regions, which is obtained after symmetrization, is shown in Fig. 8. In the symmetric version, the double-spectral integrals over ρ st and ρ ut are taken from the fixed-t representation, ρ ts and ρ us come from the fixed-s representation, and finally ρ su and ρ tu stem from the fixed-u dispersion relation.
In summary, we can write the contribution of higher intermediate states in the secondary channel as a double-spectral representation (we suppress the explicit dependence on the virtualities): (2.35) The border functions of the double-spectral regions approach asymptotically t + (s) s→∞ −→ 9M 2 π for the "1disc" contribution or 16M 2 π for the "2disc" contribution.

Single-variable dispersion relation for two-pion contributions
When we expand the sub-process γ * γ * → ππ into partial waves, we obtain a polynomial in the crossedchannel Mandelstam variables. This means that we neglect the crossed channel cut of the "1disc" or "2disc" boxes, reducing them effectively to triangle (in the case of "1disc" boxes) and bulb topologies (in the case of "2disc" boxes), as illustrated in Fig. 9. After having applied the approximation, there is no way to distinguish e.g. in Fig. 9g between contributions coming originally from ρ st or ρ su . Therefore, we discuss in the following what kind of single-variable dispersion relation is appropriate in the case of a partial-wave expanded input for the sub-process. Consider again the situation for a fixed-t dispersion relation with the corresponding Mandelstam diagram in Fig. 7. When constructing the Mandelstam representation, we selected from this representation only the contributions from ρ st and ρ ut . After the partial-wave expansion, however, we are no longer able to drop the incomplete contributions from ρ us and ρ su . Instead, let us assume that the neglected contributions from these two double-spectral regions are small: they are only due to the higher thresholds 9M 2 π (in the case of "1disc") or 16M 2 π (in the case of "2disc"). Furthermore, their discontinuities, being generated by multi-particle intermediate states, are phase-space suppressed. Instead of combining the completely reconstructed double-spectral regions from fixed-s, fixed-t, and fixed-u representations, we can simply sum all contributions from all three fixed-(s, t, u) representations. Apart from the neglected higher cuts, each double-spectral contribution appears twice in this sum. The appropriate representation is therefore one half the sum of fixed-(s, t, u) representations: (2.36) In the limit of infinitely heavy intermediate states in the crossed channel this relation is exact. In particular, the dominant ππ-rescattering contributions that we consider in this paper can be understood as a unitarization of the pion pole in the crossed channel on a partial-wave basis. In this case, the dispersion relation (2.36) provides a model-independent representation of the contribution of resonant effects in the ππ spectrum.
µ integration region is selected and defines the external kinematics.
Mandelstam diagram for the selected kinematics of point p. The double-spectral regions for the pion box are shown. Lines of fixed s, t, and u running through the point p with (s, t, u) = (q 2 3 , q 2 2 , q 2 1 ) are shown. They do not intersect any double-spectral region. In the case of (g − 2) µ kinematics, we are interested only in space-like momenta of the virtual photons. The lines of fixed-(s, t, u) therefore never enter the double-spectral regions, see Fig. 10. This implies that a partial-wave expansion is valid without restrictions. This is true even in the case of the pion box, which provides the opportunity to check the partial-wave formalism in a case where we know the full result. However, one has to bear in mind that the double-spectral representation for the pion box differs from the "1disc" and "2disc" boxes: in the case of the pion box, there are only two-pion intermediate states, hence only three box topologies exist and there are only three double-spectral regions. Each fixed-(s, t, u) representation reconstructs already all three double-spectral contributions, so that the full result can be obtained from a fixed-(s, t, u) dispersion relation separately. Hence, in a symmetrized version for the pion-box one has to take one third of the sum of fixed-(s, t, u) representations: and the relation is exact.

Sum rules for the BTT scalar functions
The Lorentz decomposition of the HLbL tensor is only unique up to transformations that do not introduce kinematic singularities, hence there is a fair amount of freedom in choosing a particular representation. One important aspect of such transformations concerns the fact that the different mass dimensions of the Lorentz structures imply different mass dimensions of the scalar functions Π i , which must be reflected in a different asymptotic behavior. Indeed if we assume, as it is natural, a uniform asymptotic behavior of the whole HLbL tensor, i.e. in all Mandelstam variables and for all tensor components, this implies that functions multiplying Lorentz structures of higher mass dimension should fall down even faster for asymptotic values of the Mandelstam variables. In order to have a predictive framework, we require, that all BTT scalar functions satisfy unsubtracted (i.e. parameterfree) dispersion relations, and in particular that those multiplying the Lorentz structures with lowest mass dimensions fall down like the inverse of the Mandelstam variables at infinity. This hypothesis, which will be tested later on, implies that the HLbL tensor behaves asymptotically as and that the BTT scalar functions behave (up to logarithmic corrections) according to: coefficient functions and a second, equally valid set of functions will vanish as a consequence of the sum rules (also known as "superconvergence relations" [62]). Consider for example Π 7 for fixed t = t b = q 2 2 + q 2 4 . At this kinematic point, the Tarrach singularity is absent and Π 7 =Π 7 is unambiguously defined (up to the redundancy in 4 space-time dimensions), see [31]. It fulfills an unsubtracted fixed-t dispersion relation: 4 where s 0 and u 0 denote the threshold in the respective channel. Due to the asymptotic behavior, s Π 7 fulfills an unsubtracted dispersion relation as well: We subtract this equation once using as well as in the u-channel integral to obtain: The comparison with (2.40) gives the following sum rule: In the case of Π 31 , an even higher-degree sum rule is fulfilled. Starting from the unsubtracted fixed-s dispersion relation with s b = q 2 3 + q 2 4 , two subtractions lead to (2.47) Both large brackets have to vanish, producing two independent sum rules for Π 31 . We have verified these sum rules explicitly in the case of sQED, see Sect. 3.2.

Relation to observables
In Sect. 2.2, we have derived the form of the dispersion relation for general two-pion contributions to (g − 2) µ , writing the results (2.36) and (2.37) for a generic BTT function Π i . In a next step, we want to use this dispersion relation for the actual input in the (g − 2) µ master formula (2.25). Our goal is to establish via unitarity a relation between the two-pion contribution to (g − 2) µ and helicity amplitudes for the sub-process γ * γ ( * ) → ππ.
While the BTT decomposition solves the problem of kinematic singularities, the 54 scalar functions Π i have the disadvantage to form a redundant set: there are 11 Tarrach redundancies [31] and two further ambiguities in four space-time dimensions [58], as the number of helicity amplitudes for HLbL scattering is 41. Furthermore, in the on-shell limit of the external photon contributions from its longitudinal polarization must not survive. This reduces the number of helicity amplitudes to 3 3 ×2 2 = 27. In the limit of (g − 2) µ kinematics, the number of independent amplitudes is further reduced to 19, see Sect. 2.1.2. Note, however, that this limit applies to the outer kinematics in the master formula and not to the imaginary parts inside the dispersion integrals, where one Mandelstam variable is integrated over and thus not fixed to (g − 2) µ kinematics (2.23).
While working with a redundant set of functions may, at first sight, seem as a minor nuisance which should in the final result take care of itself and lead to a unique and correct answer, this is not the case in our context. The origin of the problem is that (i) establishing the relation between the physical observables (i.e. the helicity amplitudes) and the BTT functions, (ii) projecting on partial waves, and (iii) writing down dispersion relations, are not necessarily commuting operations. In [31], we constructed single-variable dispersion relations that are free of the Tarrach redundancies. However, for most of the scalar functions, we only found a dispersion relation in one of the three channels, which was sufficient to obtain the dispersive reconstruction of the pion box. For general two-pion contributions (2.36) we need all three fixed-(s, t, u) dispersion relations. Furthermore, the fact that longitudinal polarizations of the external photon do not contribute is not immediately manifest in dispersion relations for the BTT functions Π i . In order to solve all these problems we must construct another basis that is appropriate for the kinematics in the imaginary parts of the dispersion integrals. The scalar functions of this basis, which we will callΠ i , are in one-to-one correspondence with the 27 singly-on-shell helicity amplitudes.
In Sect. 2.4.1, we explain how to derive these singly-on-shell basis functionsΠ i . In the construction, we make use of the sum rules for the BTT scalar functions Π i , derived in Sect. 2.3. Readers who are not interested in the technical details of the derivation may skip the following subsection and jump directly to Sect. 2.4.2, where we present the solution for theΠ i functions. As a by-product in the derivation of the singly-on-shell basis, we find a set of 15 sum rules for fixed-t kinematics, presented in Sect. 2.4.3. These physical sum rules are of relevance for the construction of the input on γ * γ ( * ) → ππ and can be considered a generalization of certain sum rules for forward HLbL scattering from [63].

Construction of the singly-on-shell basis
The most efficient way to obtain a representation for the two-pion HLbL contribution to (g − 2) µ involving only physical helicity amplitudes is the construction of a basis {Π i } for singly-on-shell kinematics that can be used together with unsubtracted single-variable dispersion relations. In such a basis, contributions from longitudinal polarizations of the external photon are manifestly absent. As we will see, this construction is possible due to the presence of the sum rules for the BTT scalar functions derived in Sect. 2.3. The rather surprising fact that contributions from unphysical polarizations are not trivially absent in a representation involving redundancies is explained in App. E.1.
Let us define the transformation from the BTT functions toΠ i as the 54 × 54 matrix t: In terms of the Lorentz structuresT µνλσ i , both the 11 Tarrach redundancies (r) and the two ambiguities in four space-time dimensions (a) can be written as linear relations: (2.49) Next, we study the unphysical polarizations: they multiply structures that are proportional to q 2 4 or q σ 4 . Hence, we determine all linear dependencies of the tensor structures in the limit q 2 4 , q σ 4 → 0, which leads to a matrix u of rank 25: If we join u with the two 4d ambiguities, the rank is 27: Since 54 − 27 = 27, this is consistent with the fact that in the singly-on-shell limit there are 27 independent helicity amplitudes. In this limit, the 11 Tarrach redundancies r are linearly dependent on a and u. Moreover, in the singly-on-shell limit the transformations u and a can be interpreted as an ambiguity in the scalar functions:Π where we denote byū the 54 × 27 matrix (u, a).
We consider now the limit q 2 4 → 0 and t → q 2 2 , which is relevant for the fixed-t dispersion relation. For a suitable choice of u and a, we still have rank(ū) = 27 in this kinematic limit. The goal is now to find all linear combinations of scalar functionsΠ i that are invariant under the transformation (2.52) and satisfy an unsubtracted dispersion relation. Hence, we have to determine the matrixp, such that for arbitrary ∆ j , which corresponds to the null-space ofū: However, the requirement thatp kiΠi satisfy an unsubtracted dispersion relation does not allow arbitraryp. Here, the sum rules for the BTT scalar functions Π i , derived in Sect. 2.3, are employed in an essential way: the linear combinations must only involve coefficients p kj that depend linearly on s for j ≥ 7 or at most quadratically for j ∈ {31, . . . , 36}, because the scalar functions Π j satisfy the linear sum rule for j ≥ 7 and the quadratic sum rule for j ∈ {31, . . . , 36}. Meanwhile, the coefficients p kj can have an arbitrary dependence on the virtualities q 2 i , which in the dispersion relation are fixed external quantities. Hence, we write and bear in mind the mentioned restrictions for p kj1 and p kj2 . Solving this linear algebra exercise is the major problem of the calculation. With the help of computer algebra, we obtain a 42 × 162 matrix (p kj0 , p kj1 , p kj2 ), whose contraction p kj (s) has again rank 27 and is in one-to-one correspondence with the 27 singly-on-shell helicity amplitudes. In a last step, we consider the limit s → q 2 3 (which is now equivalent to q 4 → 0) and search for linear relationŝ 11, 13, 14, 16, 17, 39, 50, 51, 54, (2.57) for all the functions contributing to (g − 2) µ , where the coefficients b ik are functions of the virtualities q 2 1 , q 2 2 , and q 2 3 . The solution of this system is not unique: p kj is a 42 × 54 matrix of rank 27, hence there exist 15 null relations again with coefficients n ik depending only on q 2 1 , q 2 2 , and q 2 3 . With the constructed solution for p kj , we can build a singly-on-shell basis by selecting a convenient set of 27 independent linear combinations. We choose the basis functionsΠ i in such a way that only the first 19 contribute to (g − 2) µ : 1. They are linear combinations of the BTT functions Π i for fixed-t with coefficients depending on s only in such a way that the sum rules in Sect. 2.3 allow an unsubtracted dispersion relation for (2.60) 2. In the limit q 4 → 0, a subset of 19 functions reproduces the input for the master formula (2.25) 3. They are free from Tarrach redundancies [31] and the ambiguity in four space-time dimensions [58]. 4. A basis change relates them to the 27 singly-on-shell helicity amplitudes, hence the imaginary parts in the dispersion integrals (2.60) can be expressed in terms of physical helicity amplitudes for γ * γ * → ππ and γ * γ → ππ.
The first point reflects the need to obtain a parameter-free prediction for the two-pion contribution to (g − 2) µ . The second point implies that we can construct a dispersive representation forΠ i of the form (2.36) (or (2.37) for the pion box), by summing fixed-(s, t, u) representations. The fixed-t representation is given directly by (2.60), while fixed-s and fixed-u representations follow from the crossing relations (2.16). The last two properties mean that we can relate the two-pion contribution to (g − 2) µ to observable quantities. In particular, longitudinal polarizations for the external photon must drop out in the limit q 2 4 → 0. The 19 functions contributing to (g − 2) µ can be written as (for q 2 4 = 0 and t = q 2 2 ) are given explicitly in App. D. The coefficientsd ij andd ij depend only on q 2 1 , q 2 2 , and q 2 3 . To verify that the functionsΠ i fulfill unsubtracted fixed-t dispersion relations, we observe where we have used the sum rules for the BTT functions: and written both channels schematically as one integral. This proves that the dispersion relation foř Π i is indeed fulfilled. In particular, the limit s → q 2 3 provides a fixed-t representation forΠ g i , the input for the (g − 2) µ master formula. The solutions for fixed-s and fixed-u follow immediately from the crossing relations (2.16) and (2.18).
Unfortunately, it turns out that it is not possible to find a representation for the functionsΠ i with coefficientsd ij andd ij in (2.63) free of all kinematic singularities. This is a final relic of the redundancy in the tensor decomposition which is, however, not a real problem at all. Indeed the contribution of ∆ i and∆ i in the dispersion relation forΠ g i vanishes due to the sum rules, and the same is true for the residue of any kind of kinematic singularity in the coefficientsd ij andd ij . The residue is defined in terms of physical quantities only and can thus be subtracted explicitly, to obtain a representation that is manifestly free of kinematic singularities.
Using the above sum rules, we can optimize the representation to a certain degree. We have chosen our preferred representation in App. D according to the following criteria: • We have avoided for scalar functions Π i that receive S-wave contributions to mix into other functions in (2.62).
• We have made the singularity structure of the coefficientsd ij andd ij as simple as possible.
• We have optimized the convergence of the partial-wave representation of (g − 2) µ for the pion box.
The minimal singularity structure for the coefficientsd ij andd ij consists of simple poles in 1/(q 2 1 + q 2 3 ) and singularities of the type 1/λ(q 2 is the Källén triangle function. The first singularity lies on a straight line outside the (g − 2) µ integration region, see Fig. 3. Writing we see that the second type of singularity lies on the border of the (g − 2) µ integration region. In the (g − 2) µ master formula (2.25), we subtract the residue of this singularity at r = 1, which vanishes due to the sum rules, to obtain a representation without any kinematic singularities in the (g − 2) µ integration region.

Physical sum rules
In the derivation of the singly-on-shell basis functionsΠ i , we have encountered the 15 null relations (2.58), which lead to sum rules involving only physical (singly-on-shell) quantities. We build the 15 functions By using the null relations (2.58), we subtract zero on the right-hand side and obtain where the second equality follows from the fact that p kj (s) = p kj (q 2 3 ) is constant for j < 7. For j ≥ 7, p kj (s) is linear in s or quadratic for j ∈ {31, . . . , 36}. Hence, we can write These 15 sum rules are special: they are free of any ambiguity and only involve physical helicity amplitudes, i.e. amplitudes with a transversely polarized external photon. They can be used to modify the fixed-t representations (2.62) of the 19Π i functions contributing to (g − 2) µ . The 15 sum rules can be written in very compact form in terms of the singly-on-shell basis functionsΠ i , defined in App. D: 8,9,10,12,13,16,20,21,22,23,24, where fixed-t kinematics is implicit. These sum rules are related to certain sum rules for forward HLbL scattering derived in [63], although we have derived them for a different kinematic situation (non-forward scattering but q 2 4 = 0). A detailed comparison is provided in App. E.2.

Helicity amplitudes and partial-wave expansion
In order to determine the two-pion contribution to the scalar functions in the master formula (2.25), we write fixed-(s, t, u) dispersion relations of the form (2.36), where we take only the contribution of the two-pion intermediate state to the imaginary parts into account. The scalar functions that fulfill single-variable dispersion relations and reproduce the scalar functions in the master formula are given in (2.62). The last missing piece in the formalism for two-pion contributions to (g − 2) µ is thus the link with helicity amplitudes and partial waves for γ * γ ( * ) → ππ. Unitarity determines the imaginary part of the scalar functions, which is the input in the dispersion relations, and is most conveniently expressed in the basis of helicity amplitudes, expanded into partial waves: for helicity partial waves the unitarity relation is diagonal. Furthermore, the input on γ * γ ( * ) → ππ is available in the form of helicity partial waves: these are in principle observable quantities, even though given the absence of double-virtual data they will have to be reconstructed dispersively by means of the solution of a system of Roy-Steiner equations [28,31,55]. In Sect. 4, we will provide a first estimate of the two-pion rescattering contribution by solving the Roy-Steiner equations for S-waves, using a pion-pole LHC and ππ phase shifts based on the inverse-amplitude method [64][65][66][67][68][69].
The step from the singly-on-shell basis to the basis of helicity amplitudes for HLbL is again rather tedious. The helicity amplitudes can be easily expressed in terms of BTT scalar functions or the singly-on-shell basis by contracting the HLbL with appropriate polarization vectors, but expressing the scalar functions in terms of helicity amplitudes requires the analytic inversion of a 27 × 27 matrix, which is a formidable task. Here, we present the solution to this problem and discuss the subtleties of the partial-wave expansion in connection with (g − 2) µ . In Sect. 2.5.1, we recall the definitions for the helicity amplitudes from [31]. In Sect. 2.5.2, we comment on the implication of the sum rules for the partial waves. In Sect. 2.5.3, we discuss the result for the dispersion relation in terms of helicity partial waves, generalizing the S-wave result of [28] to arbitrary partial waves. Some technical parts of the calculation are relegated to App. F.

Unitarity relation in the partial-wave picture
Although the direct inversion of the 27 × 27 matrix is feasible, see App. F.2 for a summary of how to achieve this task, there is a more elegant way to derive the partial-wave unitarity relation without the need for a full inversion. We checked that both derivations lead to identical results, but pursue the latter, more physical approach in the main part of the paper.
The strategy that avoids the inversion of the matrix describing the basis change relies on the following idea: by expanding the sub-process γ * γ * → ππ into helicity partial waves, we can explicitly calculate the phase-space integral in the unitarity relation and determine the imaginary part as a sum of products of helicity partial waves. The phase-space integrals become more and more complicated for higher partial waves, but due to the fact that unitarity is diagonal for helicity partial waves, the contribution of arbitrary partial waves is determined as soon as the S-, D-, and G-wave discontinuities are calculated.
In phenomenological applications, we expect the contribution of partial waves beyond D-waves to be negligible. However, the calculation of higher partial waves allows us to check the convergence of the partial-wave series to the full result in the test case of the pion box and provides a very strong test of the formalism for the single-variable partial-wave dispersion relations. The numerical checks of the convergence will be discussed in Sect. 3.3.
In the following, we define the helicity amplitudes for HLbL and the sub-process γ * γ ( * ) → ππ. The definitions of angles and polarization vectors can be found in [31].
The helicity amplitudes of γ * γ * → ππ are defined as For two off-shell photons, there are in principle 3 2 = 9 helicity combinations. However, due to parity conservation and with our convention for the polarization vectors, we have the relation which implies that only 3 2 −1 2 + 1 = 5 amplitudes are independent: Similarly, for the HLbL helicity amplitudes, defined by there are 3 4 helicity amplitudes, but only 3 4 −1 2 + 1 = 41 independent ones. We introduce rescaled helicity amplitudes that remain finite in the limit q 2 i → 0: where and ξ i refers to the normalization of the longitudinal polarization vectors. Since only theH amplitudes appear in the final results, this procedure avoids any confusion that might originate from a particular choice of normalization. We define the helicity partial-wave expansion for γ * γ * → ππ bȳ where m = |λ 1 −λ 2 |, z is the cosine of the scattering angle, and d J m 1 m 2 (z) denotes the Wigner d-function. For HLbL, we expand the helicity amplitudes as follows into partial waves: where Unitarity is diagonal for helicity partial waves, i.e.
where S is the symmetry factor of the two pions and account for the sign convention in (2.78). We find the relation where the ratio of η f factors compensates the sign (−1) λ 3 +λ 4 from (2.73). The HLbL tensor is written in terms of the redundant BTT Lorentz decomposition as For fixed t = q 2 2 and q 2 4 = 0, we have defined the singly-on-shell basis functionsΠ i . The helicity amplitudes form a basis of the HLbL tensor, hence (2.84) The coefficients c ij contain 13 redundancies, thec ij still two (in four space-time dimensions). In the relation forΠ i , fixed-t kinematics is implicit and the coefficientsč ij are free from redundancies. We define the following "canonical" ordering of j:  The meaning of these subsets is the following: the subset {l j } j corresponds to helicity amplitudes with Hj = ±H j , wherej := {λ 1 , λ 2 , −λ 3 , −λ 4 }. For the subset {k j } j , the Wigner d-functions for j andj are identical up to a sign, while for the subset {n j } j this is not the case.
The imaginary parts of the scalar functions are given by where the signs The explicit Wigner d-functions are where the signs are due to the use of relation (2.89).
In order to identify the coefficientsč il j and (č ik j + ζ jč ik j ), it is sufficient to know the contribution to the unitarity relation from the lowest partial waves h J l j and h J k j (which are either S-or D-waves). However, as the Wigner d-functions d J n j are different from d J n j , we need to know the contribution from the two lowest partial waves h J n j in order to identify the coefficientsč in j andč in j separately. Therefore, the generalization to arbitrary partial waves is possible as soon as the contributions from S-, D-, and G-waves are determined.
The explicit calculation of the partial-wave unitarity relation involves rather complicated phasespace integrals, see App. F.1. By calculating the fully-off-shell unitarity relation, projecting onto BTT, and working out the imaginary parts of the functionsΠ i , we have verified explicitly that the coefficientš c ij for the longitudinal polarization λ 4 = 0 vanish. Therefore,č ij is effectively an invertible 27 × 27 matrix. As mentioned above, we have also computed the matrixč ij by direct inversion of the basis change from helicity amplitudes to the scalar functions, see App. F.2. The fact that the result agrees with the one from the phase-space calculation provides a very strong cross check, and in addition the full inversion allows one to separate theč ik j andč ik j coefficients.

Approximate partial-wave sum rules
Before returning to the final result, we comment on the role of the sum rules in the context of a partial-wave expansion. In Sect. 2.4.3, we have derived a set of 15 sum rules for theΠ i functions, which, after a basis change, can be written in terms of the 27 singly-on-shell helicity amplitudes for HLbL scattering. By construction, these sum rules only hold true for the full helicity amplitudes.
In particular, when expanding the imaginary part of the helicity amplitudes into partial waves and truncating the partial-wave series, there is no reason why the sum rules should still be satisfied exactly: sum-rule violations of a size consistent with higher partial waves are expected, so that the sum rules are fulfilled only approximately. This has some important consequences.
Due to the presence of the sum rules, the formal relation between the master formula inputΠ i at q 4 = 0 and the singly-on-shell basis functionsΠ i is not unique, but can be modified by linear combinations of the sum rules. If the sum rules hold exactly, all these representations are equivalent. Violating the sum rules by a truncation of the partial-wave series implies that a dependence on the precise representation of theΠ i functions is introduced. Our preferred representation of theΠ i functions, discussed in Sect. 2.4.2 and App. D, leads to a fast convergence of the partial-wave expansion in the test case of the pion box, see Sect. 3, but we also checked other variants and convinced ourselves in each case that indeed the sum-rule violations are consistent with a meaningful partial-wave expansion and only slight losses in the rate of convergence.
The dependence on the representation of theΠ i functions or the violation of the sum rules concerns only the truncated higher partial waves. Hence, we can reverse the argument: with the assumption that sufficiently high partial waves are negligible, the included partial waves have to fulfill the sum rules, which also removes the dependence on the representation. This can be used as a check of or even a constraint on the input for the γ * γ * → ππ helicity partial waves, in a similar way as sum rules for forward HLbL scattering have been used to derive constraints on transition form factors of higher resonances [63,70,71].
Out of the 15 sum rules, only a single one involves helicity amplitudes starting with S-waves. If we truncate the partial-wave expansion after S-waves, this sum rule reads Im h 0 00,++ (s ) + higher waves, (2.91) where λ 12 (s) := λ(s, q 2 1 , q 2 2 ) denotes the Källén triangle function. Verifying that the corresponding sum rule is approximately fulfilled for the γ * γ * → ππ amplitudes constructed in Sect. 4 provides an important check on the calculation. In fact, it is precisely this sum rule that proves that the S-wave result derived here based on the BTT formalism and the one from [28] are equivalent. We note that in the limit of forward kinematics the sum rule (2.91) reduces to the S-wave approximation of the sum rule (27b) in [63].

Result for arbitrary partial waves
The calculations of the previous sections allow one to reconstruct the full result for the dispersion relation for HLbL two-pion contributions to (g − 2) µ . The imaginary part of the functionsΠ i , which have to be inserted into the dispersion integrals, are provided by (2.87). Evaluated at s = q 2 3 , the dispersion relations give the s-channel contribution for the fixed-t representation of all 19Π i functions that contribute to (g − 2) µ . Using the crossing relations (2.16) and (2.18), we obtain the five other contributions: the u-channel contribution for fixed-t as well as both channels in the fixed-s and fixed-u representations. Hence, all six integrals in a dispersion relation for the functionsΠ i of the form (2.36) or (2.37) can be calculated.
The crucial ingredient in this calculation is the basis changeč ij from scalar functions to helicity amplitudes, which enables the generalization of the S-wave result of [28] to arbitrary partial waves. The matrixč ij contains two types of ostensible kinematic singularities: 1. The kinematic singularities of the singly-on-shell basisΠ i are present, as explained in Sect. 2.4.2.
In the dispersion relation, their residues vanish due to the sum rules, hence they can be subtracted explicitly in the master formula for (g − 2) µ .
2. Additional kinematic singularities (−q 2 2 ) −n/2 , n = 1, . . . , 4, show up in the coefficientsč ij . They are introduced by the basis change to helicity amplitudes, i.e. they cancel against kinematic zeros in the helicity amplitudes, present in (2.87) in the Wigner-d functions for fixed-t kinematics.
Unfortunately, the matrixč ij is too lengthy to be shown here in full, but is provided as supplemental material in the form of a Mathematica notebook. 5 In contrast, the explicit results for the two-pion dispersion relation in the S-wave approximation are very compact: where the dependence of the helicity amplitudes on the virtualities is not written explicitly. This result agrees with [30]. It slightly differs from the S-wave result presented in [28], but, as explained in the previous section, this difference is precisely of the form of the sum rule (2.91) and thus simply related to a different choice of basis. The above result is given in a form that corresponds to the dispersion relation (2.36). In order to apply it to the pion box, one has to use (2.37), hence the dispersion integrals in (2.92) need to be multiplied by a factor 2/3. For the proper evaluation of the ππ-rescattering corrections, the contribution of the pion box to the partial waves has to be subtracted: we define the operator S, which takes care of the symmetry factor and the subtraction of the pole × pole term [28]. The imaginary part for the ππ-rescattering contribution is then given by The superscripts refer to charged (c) and neutral (n) pions, respectively, and N J,λ i λ j denotes the partial-wave projection of the pure pion-pole term, explicitly given in App. G.

Summary of the formalism
Arguably the most important result of this paper, especially in view of future applications and generalizations, concerns the derivation of theΠ i functions, which allows us to establish a direct correspondence between singly-on-shell helicity amplitudes and the scalar functionsΠ i that enter the master  formula (2.25) for the HLbL contribution to (g − 2) µ . The key quantities in this construction are the various scalar amplitudes, a glossary of which is provided in Table 1, including a reference to the equation where they are defined and a short definition and explanation. They can be roughly divided into four classes: first, the Π i andΠ i are related to the general BTT decomposition of the HLbL tensor, irrespective of any application to (g − 2) µ or dispersion relations. Second, theΠ i andΠ i isolate the functions actually relevant for (g − 2) µ , by forming suitable subsets and taking the appropriate kinematic limit, but are otherwise still completely general. Third and fourth, theΠ i are constructed as the crucial intermediate step in the derivation of single-variable partial-wave dispersion relations, by eliminating redundancies in the representation and thereby allowing a well-defined transition to helicity amplitudesH j . In combination with partial-wave unitarity, this last step completes the derivation of the dispersion relation for two-pion intermediate states in the HLbL contribution to (g − 2) µ .

The pion box: test case and numerical evaluation
The interest in the pion box is twofold. On the one hand, it gives a unique meaning to the notion of a pion loop, by virtue of its dispersive definition as two-pion intermediates with a pion-pole LHC, and is expected to provide the most important contribution to HLbL scattering beyond the pseudoscalar poles. Phenomenologically, the pion box is fully determined by the pion vector form factor, which allows us to pin down its numerical value to very high precision, as we will show in Sect. 3.1 including an error analysis for the form factor input. On the other hand, the pion box constitutes an invaluable test case for the partial-wave formalism that we have developed in Sect. 2. Given a certain representation of the pion vector form factor, the full pion box is known exactly, see App. C. Since the partial-wave expansion and the single-variable dispersion relations are valid not only for the rescattering contribution but also for the pion box, provided the correct prefactor in (2.37) according to the counting of double-spectral regions is taken into account, we can use the pion box to check whether the partial-wave representation converges to the full result upon resummation of the partial waves, and we can study the details of the convergence behavior numerically.
In a similar way, the pion box provides a test case for the sum rules for the HLbL scalar functions. In Sect. 3.2 we demonstrate that they are indeed fulfilled, which is a prerequisite for the unsubtracted single-variable dispersion relations derived in Sect. 2. In Sect. 3.3, we investigate the convergence behavior of its partial-wave representation and discuss the implications for applications beyond the pion box, such as the ππ-rescattering contribution discussed in Sect. 4.

Evaluation of the full pion box
For the numerical evaluation of the pion box, the representation in terms of Feynman-parameter integrals given in App. C proves most efficient. This representation is based on the equivalence of the pion box with the FsQED amplitude [28], which we proved in [31]. It requires the numerical evaluation of two-dimensional Feynman integrals with the pion vector form factor as the only input. For a reliable evaluation of the pion-box contribution to (g − 2) µ , we therefore need a precise representation of the pion vector form factor in the space-like region.
Since about 95% of the final pion-box (g − 2) µ integral originate from virtualities below 1 GeV, it is most critical that the low-energy properties be correctly reproduced. Experimentally, the available constraints derive from e + e − → π + π − data, which determine the time-like form factor [74][75][76][77][78][79], and space-like measurements by scattering pions off an electron target [80,81]. We have also checked that our representation is consistent with extractions of the space-like form factor from e − p → e − π + n data [82][83][84][85], although due to the remaining model dependence of extrapolating to the pion pole we do not use these data in our fits. To obtain a representation that allows us to simultaneously fit spaceand time-like data, and thereby profit from the high-statistics form factor measurements motivated mainly by the two-pion contribution to HVP, we adopt the formalism suggested in [86,87] (similar representations have been used in [88][89][90][91][92][93]), whose essential ingredients will be briefly reviewed in the following.
The form factor is decomposed according to The Omnès factor would provide the exact answer if only the elastic ππ channel contributed to the unitarity relation of the form factor. It is fully determined by the P -wave phase shift δ 1 1 . Next, G ρω describes the isospinviolating coupling to the 3π system, which becomes relevant in the vicinity of the ρ peak as reflected by ρ-ω mixing. In practice, a one-parameter ansatz proves indistinguishable from a dispersively improved version that eliminates the imaginary part below the 3π threshold [86,87]. Finally, G inel parameterizes the effect of higher inelastic channels. We use a conformal mapping where s ωπ = (M ω +M π ) 2 denotes the threshold where phenomenologically 4π inelasticities first start to set in and the second parameter is fixed at s 1 = −1 GeV 2 . The ππ phase shift is taken from the extended Roy-equation analysis of [94], which determines δ 1 1 up to s m = (1.15 GeV) 2 in terms of its values at s m and s A = (0.8 GeV) 2 . Our representation thus involves 3 + p free parameters: the ρ-ω mixing parameter ρω , the two values of the phase shift at s m and s A , and p parameters from the conformal expansion of G inel . This representation ensures that the form factor behaves as 1/s asymptotically as long as the phase shift approaches π, up to logarithms in agreement with the expectation from perturbative QCD [95][96][97][98][99]. We impose this asymptotic behavior by smoothly extrapolating δ 1 1 to π from the boundary s m of the applicability of the Roy solution, but checked that introducing effects from ρ , ρ excitations as suggested in [40] does not impact the space-like form factor. The form of G inel can be further constrained by requiring that the imaginary part exhibit the expected P -wave behavior and respect the Eidelman-Łukaszuk bound [100], but again the impact on the space-like form factor proves to be small.
We fit this representation simultaneously to the space-like data from [81] as well as one of the timelike data sets [74][75][76][77][78][79] (restricted to data points below 1 GeV). Moreover, we varied s 1 , p = 1, 2, and constructed an error band for the uncertainties in δ 1 1 apart from the phase shifts at s m and s A . We find that the results for the space-like form factor are extremely stable to all these variations, the largest effect being produced by the differences between the time-like data sets. For the accuracy required in HLbL scattering we can therefore simply take the largest variation among them as an uncertainty estimate, without having to perform a careful investigation of the statistical and systematic errors that are crucial when combining the different data sets for HVP. The result for the space-like form factor is shown in Fig. 11, leading to a numerical evaluation 6 for the pion box of a π-box µ = −15.9(2) × 10 −11 . (3.5)

Verification of sum rules
In Sect. 2.3, we have presented sum rules for the BTT scalar functions that follow from a uniform asymptotic behavior of the HLbL tensor and ensure the independence from the choice of the tensor basis. These sum rules prove essential for the derivation of single-variable dispersion relations that can be used with input on the γ * γ * → ππ helicity partial waves. Furthermore, an important consequence of the BTT sum rules are the physical sum rules in Sect. 2.4.3, which can be expressed in terms of helicity amplitudes. An important test case for our partial-wave formalism is the pion box: the fact that we know the full result allows us to test the convergence behavior of the partial-wave approximation. Before turning to the tests of the full formalism in Sect. 3.3, here we check that the sum rules as a necessary prerequisite for the single-variable dispersion relations are indeed fulfilled in the case of the pion box. Due to the equivalence of the pion box with the FsQED amplitude [28,31], these tests can be directly performed with sQED.
The BTT sum rules can therefore be formulated as the requirement that the functions in (3.6) fulfill an unsubtracted Mandelstam representation. In contrast, in [31] we only verified that subtracted Mandelstam representations which follow from unsubtracted ones for the BTT functions Π i are actually fulfilled. In analogy to [31], we extract the sQED double-spectral densities of these functions from the explicit expression of the loop calculation in terms of Passarino-Veltman amplitudes [102,103]: in such a decomposition into scalar loop functions the double-spectral densities are given by the coefficients of the D 0 functions times the D 0 spectral densities. By inserting the double-spectral densities into an unsubtracted Mandelstam representation of the form we have verified numerically that the functions (3.6) are reproduced. Surprisingly, even q 1 · q 2 q 3 · q 4Π21 fulfills an unsubtracted Mandelstam representation, which is not expected from the mass dimension and implies an even higher sum rule in the case of sQED. Single-variable dispersion relations then follow from the Mandelstam representation in the appropriate limit, including the explicit cases discussed in Sect. 2.3. While the imaginary parts for the single-variable dispersion relations extracted in this way need to be calculated numerically, it is also possible to obtain analytic expressions starting from a Feynman-parameter representation of the BTT functions. The results again confirm the validity of the sum rules, in agreement with the more general approach via the Mandelstam representation.
Although the sum rules for the BTT scalar functions are crucial ingredients in the derivation of the single-variable dispersion relations, the physical sum rules (2.71) have a more direct significance as they are formulated in terms of physical quantities for the kinematics of the (g−2) µ single-variable dispersion integrals, and thus allow one to impose constraints on the γ * γ * → ππ helicity amplitudes used as input for a numerical evaluation. We have verified that these sum rules are fulfilled in the case of the pion box, by extracting the imaginary parts of theΠ i functions from the sQED calculation and calculating the integrals numerically. These sQED tests thereby allow one to establish the validity of nearly all sum rules-except for the last one involvingΠ 17 −Π 18 +Π 19 =Π 40 −Π 41 , which vanishes identically in sQED. It should be stressed that the underlying assumptions follow solely from demanding a uniform asymptotic behavior of the HLbL tensor, but as the discussion in App. E.2 shows, similar conclusions can be drawn from Regge models as well. Together with the explicit checks in the case of sQED there is therefore compelling evidence for our assumptions regarding the asymptotic behavior of the HLbL tensor.

Convergence of the partial-wave representation
In the following, we perform tests of the helicity partial-wave dispersion relations developed in Sect. 2 by applying the formalism to the pion box. In this case, a dispersion relation of the form (2.37) has to be used in order to account for the fact that only three different double-spectral regions are present. We emphasize that in this test case each single-variable dispersion relation reconstructs the full pion box. Therefore, we can test the three channels separately-each must converge to the full result upon resummation of the partial-wave series.
The input for the γ * γ * → ππ helicity partial waves in the case of the pion box is given by the partial-wave projection of the pure pion-pole terms, see App. G. In order to simplify the convergence checks, we use a simple vector-meson dominance representation for the pion vector form factor: (3.8) Such a form factor leads to a π-box, VMD µ = −16.4 × 10 −11 , which is very close to the full result obtained with the dispersive representation of the form factor discussed in Sect. 3.1. The convergence behavior of the partial-wave expansion is not affected by the details of the form factor implementation.
Since our formalism for single-variable dispersion relations is valid for arbitrary partial waves, we can extend these tests in principle to an arbitrary angular momentum J. In practice, our numerical implementation becomes less reliable for large values of J, so that we performed the numerical tests up to J = 20 and estimated the truncation error by extrapolation.
The HLbL contribution to (g − 2) µ is given as a sum of 12 terms in the master formula (2.25), which, in principle, are completely independent. However, in the case of the pion box it turns out that especially for the lower partial waves a numerical cancellation occurs that leads to a faster convergence of a µ than for the individual terms. Therefore, we define the following vector in the 12-dimensional space of the contributions to the master formula:   In order to quantify the convergence behavior, we define the following two quantities: the relative deviation between the full pion-box contribution to (g − 2) µ and its partial-wave approximation where | · | denotes the 12-dimensional Euclidean norm. Due to cancellations between the 12 terms in the master formula, δ Jmax will indicate a faster convergence than ∆ Jmax , which is more robust against cancellations. Table 2 shows the results of a detailed study of the convergence behavior of the partial-wave representation for the test case of the pion box. Both measures δ Jmax and ∆ Jmax for the deviation from the full pion box result are displayed for fixed-s, fixed-t, and fixed-u dispersion relations, as well as for the average of the three single-variable dispersion relations (2.37). In this context we use the notion of fixed-(s, t, u) as follows: it defines the dispersion relation for each of the six representatives in (2.15), while the remaining scalar functions are obtained via the crossing relation (2.16). In particular, this implies that in the so-called fixed-s evaluation, we do use a fixed-s dispersion relation forΠ 1 , but a fixed-t dispersion relation forΠ 2 and a fixed-u dispersion relation forΠ 3 .
Next, we comment on the following two observations: 1. The S-wave approximation shows a particular pattern: the fixed-s representation vanishes, while fixed-t and -u agree.
2. The fixed-s representation exhibits a slower convergence than the other two dispersion relations.
In order to understand the first point, consider the explicit S-wave representation for theΠ i functions, (2.92). We note that S-waves contribute only to the s-channel discontinuities inΠ 4 andΠ 17 , while the t-and u-channel discontinuities for these two functions start with D-waves (the situation for the other functions follows from crossing symmetry). A discontinuity in the s-channel contributes to a fixed-t and fixed-u dispersion relation, while in a fixed-s dispersion relation the integral runs only over t-and u-channel discontinuities. This means that in the fixed-s representation in Table 2, no S-wave discontinuity is encountered at all, hence in this representation the first non-vanishing contribution is obtained from D-waves. Furthermore, because the S-wave s-channel discontinuity has no angular dependence, it contributes identically to a fixed-t and fixed-u dispersion relation, which makes the fixed-t and fixed-u results in Table 2 agree at J max = 0.
The second point can be understood as follows. For each of the six representatives in (2.15) the s-channel is special with respect to the other two. This is due to the fact that the associated Lorentz structure exhibits an s-channel symmetry, either C 12 , C 34 , or both, the only special case beingΠ 39 , which is totally crossing symmetric in all three channels. For instance, the Lorentz structureT µνλσ 1 is the one that belongs to the s-channel (pseudo-scalar) π 0 -pole contribution to HLbL scattering, whilê T µνλσ 4 can be related to an s-channel scalar amplitude, which manifests itself as the S-wave s-channel ππ contribution. It is therefore not surprising that even in the case of the pion box, the s-channel discontinuity for the functions (2.15) is more important than the other two discontinuities. Since this is the discontinuity that evades the fixed-s dispersion relation, we observe a slower convergence pattern in this case.
We have performed these convergence tests not only with our preferred representation for theΠ i functions, but also with different versions that are modified by terms that vanish due to the sum rules (2.71). While the exact numbers do differ-as expected given the fact that the sum rules only hold for the full amplitudes but not the individual partial waves-the sum rule violations in the case of the pion box due to the partial-wave approximation are reasonably small and the overall picture remains the same.
To fully understand the partial-wave convergence of the pion box we also studied the remaining deviation from the full result at J = 20. Empirically, we observe that the size of the individual terms for given J is well described by a fit function 7 a π-box, PW µ,J ∼ cJ x , (3.13) which in a double-log plot produces the straight lines in Fig. 12. The fact that the first few terms do not fall on this line indicates that the form (3.13) is only asymptotic, and might also be related to the abovementioned cancellations for low J (the fit therefore excludes the points for J ≤ 6). The figure shows that the rate of convergence is actually similar for fixed-s and fixed-u, both of which yield an exponent x ≈ −3, while the fixed-t representation converges with x ≈ −4. The slower convergence of the fixed-s results seen in Table 2 is therefore a remnant of the missed S-wave contribution that leads to larger deviations for small J, not the overall rate of convergence. The resummation of the terms with J > 20 based on the fit function then removes all remaining discrepancies, providing a strong check of the partial-wave formalism developed in Sect. 2. Finally, we discuss the consequences for the application of the formalism to the case of two-pion contributions beyond the pion box, most importantly the unitarity (or rescattering) correction. The most important difference is related to the fact that for these applications, instead of (2.37), the dispersion relation (2.36) applies, where due to the different double-spectral regions an overall factor 1/2 instead of 1/3 is required. However, this means that for the rescattering contribution the slower convergence of the fixed-s dispersion relation is of no significance: let us assume that an important resonant contribution shows up in a partial wave in the s-channel. This resonance will be captured by the fixed-t and fixed-u dispersion relation (though not by the fixed-s dispersion relation). Since the full result is given by the sum of the three dispersion relations weighted by 1/2, this behavior is actually expected and the absence of the resonance in the fixed-s representation does not impact the average in the symmetric representation (2.36) (in contrast to the pion box, where the missed S-wave contribution needs to be recovered by the higher partial waves). Therefore, the average in the case of the pion box in Table 2 should rather be regarded as a worst-case scenario-for the convergence behavior in the case of rescattering contributions the fixed-t and fixed-u dispersion relations are more representative. The second important difference concerns the presence of resonances in the rescattering contribution, a feature that does not occur in the pion box. We expect the rescattering contribution to be dominated by resonant effects, whereas the convergence behavior established for the pion box can be understood as a weighting of the partial waves. In the truncated partial-wave series, the resonances in the included partial waves are fully reproduced. The approximate fulfillment of the sum rules indicates then whether neglected higher partial waves still play an important role, to the effect that the size of the sum-rule violations allows one to estimate the accuracy of the calculation.

Application: two-pion rescattering
The natural application of the partial-wave formalism developed in the main part of this paper concerns ππ rescattering effects, which can be considered a unitarization of the pure pion-pole LHC that defines the pion box. To isolate this contribution, it suffices to subtract the pure pion-pole piece in the partialwave unitarity relation, and insert for the remainder phenomenological input for the γ * γ * → ππ partial waves. The construction of such input is by itself challenging, given that direct experimental results, at least for the doubly-virtual case, are not expected in the near future.
In the on-shell case, available data on γγ → ππ [104][105][106][107][108][109] (in combination with γγ →KK [110][111][112][113][114][115][116]) are now sufficient to perform a partial-wave analysis [117], but such an approach appears unrealistic to control the dependence on the photon virtualities. However, approaches that exploit more comprehensively the analytic properties of the amplitude, see [54,55,118] for on-shell photons, can be extended towards the off-shell case with limited data input required to determine parameters, as demonstrated for the singly-virtual process in [56]. The essential features of the generalization towards the doubly-virtual case, i.e. the appearance of anomalous thresholds for time-like kinematics [57] and the modifications to tensor basis and kernel functions [28,31], have already been laid out in previous work, but the practical implementation involves a number of challenges: due to the strong coupling between the ππ/KK channels in the isospin-0 S-wave a single-channel analysis is limited to rather low energies [54,56,117,118], assumptions for the LHC and number of subtractions need to be carefully studied to reliably assess the sensitivity to the high-energy input in the dispersive integrals [55], a full analysis of the generalized Roy-Steiner equations [28,31,55] involves solving coupled S-and D-wave systems of various helicity projections, and last but not least constraints on the γ * γ * → ππ amplitudes from asymptotic behavior and the sum rules derived in Sect. 2.4.3 need to be incorporated. A full analysis along these lines will be left for future work.
To obtain a first estimate of rescattering effects, we concentrate on S-waves and consider a further simplified system: first, we use ππ phase shifts from the inverse-amplitude method, which reproduces the phenomenological phase shifts as well as the f 0 (500) properties at low energies, and in addition allows one to separate the ππ rescattering from theKK channel in a well-defined manner, all while providing a reasonable extrapolation for high energies. In addition, we restrict ourselves to a pionpole LHC in the solution of the Roy-Steiner equations, which has the advantage that the off-shell behavior is still described by the pion vector form factor. In the following, we lay out the details of this calculation, and discuss the consequences for rescattering effects in (g − 2) µ .
4.1 γ * γ * → ππ helicity partial waves from the inverse-amplitude method Unitarization within the inverse-amplitude method (IAM) [64][65][66][67][68][69] is based on the observation that elastic unitarity Im t(s) = σ π (s)|t(s)| 2 (4.1) for a ππ partial-wave amplitude t(s) implies which together with the chiral expansion t(s) = t 2 (s) + t 4 (s) + O(p 6 ) and perturbative unitarity Im t 2 (s) = 0, Im t 4 (s) = σ π (s)|t 2 (s)| 2 , (4.3) already concludes the naive derivation of the IAM prescription However, in the single-channel case the IAM approach can be justified much more rigorously based on dispersion relations, where the only approximation involves replacing the LHC by its chiral expansion [119]. In this way, one can also remedy the fact that the standard IAM fails to correctly reproduce the Adler zero [120,121], and is thus not fully consistent with chiral symmetry. The modified form of the IAM (mIAM) becomes [119] t mIAM (s) = t 2 (s) where the additional term 8 This form of the IAM thus correctly describes the low-energy phase shifts as well as resonance properties, and has indeed been used in recent years to determine the quark-mass dependence of σ and ρ resonances [122,123]. For our purposes, the single-channel IAM for ππ scattering conveniently separates the ππ channel from its mixing toKK in the vicinity of the f 0 (980) and defines a reasonable continuation to high energies, without compromising the low-energy physics. We use the 1-loop IAM with low-energy constants as specified in [123], which produces the phase shifts shown in Fig. 13. As expected, there is good agreement throughout, apart from the fact that the IAM I = 0 phase shift avoids the rise related to the f 0 (980) and the coupling to theKK channel. We also checked that the σ properties [126] are reproduced: for the pole position we find √ s σ = (0.443 + i0.217) GeV, to be compared to √ s σ = (0.441 + i0.272) GeV [127] and similar numbers from other recent dispersive extractions [118,128]. Accordingly, the width comes out a bit too low, as does the residue at the pole g σππ . This deviation is consistent with earlier IAM analyses, see e.g. [122] for the analogous calculation including the mIAM correction, and can certainly be tolerated to obtain an estimate for the HLbL rescattering contribution, which, after all, only requires the amplitude on the real axis, not the analytic continuation into the complex plane where the slight discrepancy in the width would matter most. Similarly, one can check the coupling to two photons |g σγγ /g σππ | ∼ 0.014, well in line with |g σγγ /g σππ | = 0.014 and 0.015 from [55] and [118], respectively. With the input for the ππ phase shifts specified, the γ * γ * → ππ amplitudes follow by solving the generalized Roy-Steiner equations derived in [28,31] for doubly-virtual kinematics. For the S-waves, these dispersion relations take the form (isospin indices are suppressed for the time being) Im h 0,++ (s ) , (4.8) with LHC singularities represented by the inhomogeneities ∆ 0,++ (s) and ∆ 0,00 (s). These equations can be rewritten as h 0,++ (s) ± q 2 1 q 2 2 h 0,00 (s) = ∆ 0,++ (s) ± q 2 1 q 2 2 ∆ 0,00 (s) (4.9) The new combinations still fulfill Watson's theorem [129] Im h 0,++ (s) ± q 2 1 q 2 2 h 0,00 (s) = sin δ 0 (s)e −iδ 0 (s) h 0,++ (s) ± q 2 1 q 2 2 h 0,00 (s) θ s − 4M 2 π , (4.10) so that the dispersion relation reduces to a standard Muskhelishvili-Omnès (MO) problem [130,131], whose solution reads h 0,++ (s) ± q 2 1 q 2 2 h 0,00 (s) = ∆ 0,++ (s) ± q 2 1 q 2 2 ∆ 0,00 (s) (4.11) with the Omnès function .
For convenience, we finally rewrite the result in terms of the original helicity amplitudes according to h 0,++ (s) = ∆ 0,++ (s) h 0,00 (s) = ∆ 0,00 (s) For a pion-pole LHC ∆ 0,++ (s) and ∆ 0,00 (s) simply correspond to the partial-wave projection of the Born terms, given in App. G, which shows that the dependence on the virtualities, apart from the modified kernel functions in the MO solution, is still governed by the pion vector form factor. In particular, the corresponding factor F V π (q 2 1 )F V π (q 2 2 ) can be moved out of the integrals in (4.13), so that one can simply calculate a reduced amplitude, with the dependence on the pion form factors fully factorized. Further, in the solution of Roy-Steiner equations, a MO representation similar to (4.13) is often required for the low-energy region only, in order to match to some known high-energy input, and to this end a finite matching point is introduced [55,[132][133][134][135]. In case the amplitudes are assumed to vanish above the matching point, it effectively acts as a cutoff both in (4.13) and in the Omnès function. We will use this variant of the MO solution to estimate the sensitivity to the high-energy extrapolation of the phase shifts, referring for more details of its implementation to [55,132].
Finally, the justification why an unsubtracted representation such as (4.13) is still expected to provide a decent description is two-fold: first, by removing theKK intermediate states the Omnès functions are smoothened considerably around the nominal f 0 (980) position, which eliminates most of the need for subtractions necessary otherwise in a single-channel description to suppress the corresponding peak in the Omnès function. Second, while in general a precision description does require subtractions [54,55], we observe in the on-shell case that the results particularly for the charged channel are reasonably close to the twice-subtracted variants studied in [55], see Fig. 14 for a cutoff Λ = 1 GeV. The upper panel shows the modulus |h I 0,++ | for isospin I = 0 and I = 2, which for the unsubtracted IAM emerges remarkably close to the twice-subtracted variant in both cases. However, this agreement is largely driven by the projection of the Born term, while a more realistic picture can be obtained by considering the rotated amplitudes and subtracting the Born term in the charged channel. In this way, we find that the agreement is still very good for the charged combination, while the neutral channel is less well reproduced based on the pion-pole LHC alone, see lower panel in Fig. 14. To improve the quantitative agreement, the introduction of subtraction constants becomes unavoidable. These subtraction constants can be identified with pion polarizabilities and were taken from 2-loop ChPT [136,137] in [55]. The agreement in the charged channel implies that the corresponding sum rules for the subtraction constants, just based on the pion-pole LHC, are reasonably well fulfilled, while significant corrections are expected in the neutral channel. This interplay with the pion polarizabilities will be discussed in more detail in Sect  Table 3: Results for the S-wave rescattering contribution to (g − 2) µ in units of 10 −11 . The cutoff refers to the finite-matching-point analog of (4.13).  the HLbL integral become Im h 0,I ++,++ s; q 2 1 , q 2 2 , q 2 3 , 0 = σ π (s) 32π h I 0,++ (s; q 2 1 , q 2 2 )h I 0,++ (s; q 2 3 , 0) − c I N 0,++ (s; q 2 1 , q 2 2 )N 0,++ (s; q 2 3 , 0) , Im h 0,I 00,++ s; q 2 1 , q 2 2 , q 2 3 , 0 = σ π (s) 32π h I 0,00 (s; q 2 1 , q 2 2 )h I 0,++ (s; q 2 3 , 0) − c I N 0,00 (s; q 2 1 , q 2 2 )N 0,++ (s; q 2 3 , 0) , (4.15) with isospin factors c 0 = 4/3, c 2 = 2/3. The numerical results for the S-wave contribution then follow from (2.36) together with the dispersive representation for the scalar functions derived in Sect. 2. Since the full integration becomes numerically costly-with the dispersion integral in (4.13), the (g − 2) µ dispersion integral, and three integrals in the master formula (2.25) this would amount to a delicate 5-dimensional integral, wherein in addition the Omnès factor requires the numerical evaluation of yet another integral-we calculate the γ * γ * → ππ amplitudes on a three-dimensional grid in (s, q 2 1 , q 2 2 ) and then interpolate in the remaining 4-dimensional (g − 2) µ integration. Using up to 50 grid points in each variable the results become insensitive to the interpolation uncertainty, and we obtain the values listed in Table 3. As expected based on the size of the phase shifts, the I = 2 contribution is much smaller than its I = 0 counterpart, while in both cases the variation with respect to the cutoff amounts to about one unit. Accordingly, this estimate can be interpreted as evidence for a rescattering contribution corresponding to f 0 (500) degrees of freedom of about −9 × 10 −11 in the HLbL contribution to (g − 2) µ .
Another check on our input for γ * γ * → ππ follows from the sum rule (2.91). In fact, it is precisely this sum rule that ensures that the S-wave rescattering contribution as formulated in [28] and the one from Sect. 2.5 are strictly equivalent. Furthermore, this observation immediately suggests a way how to condense the full sum rule into a single number: the difference between the two representations amounts to a shift inΠ 4 of the size   [136,137]. The two numbers in the case of the charged-pion quadrupole polarizability refer to two different sets of low-energy constants.
and accordingly inΠ 5 andΠ 6 from crossing, so that the convolution in the (g − 2) µ integral should be done with the corresponding kernel function. Still subtracting the pion-pole terms since the validity of the sum rule in sQED is already known, we find the results for the separate contribution from h 0 ++,++ and h 0 00,++ as listed in Table 4. The expected cancellation already works at the level of 10% with S-waves only, and even better for the larger values of the cutoff. Such a 10% error on the actual rescattering contributions from Table 3 would yield a very similar uncertainty estimate as the variation observed from the cutoff dependence before. In total, these results lead us to quote a ππ,π-pole LHC µ,J=0 = −8(1) × 10 −11 (4.17) for the S-wave rescattering corrections to the pion-pole LHC.

Role of the pion polarizabilities
The low-energy behavior of the on-shell γγ → ππ amplitudes is strongly constrained by the pion polarizabilities, which therefore encode valuable information on the two-pion rescattering contributions to HLbL. The precise relation can be expressed in terms of the expansion for the Born-term-subtracted on-shell amplitudesĥ 0,++ = h 0,++ −N 0,++ . Here, α 1 −β 1 and α 2 −β 2 refer to dipole and quadrupole polarizabilities, respectively. The soft-photon zero required as a consequence of Low's theorem [139] ensures thatĥ 0,++ indeed vanishes for s → 0. Accordingly, the representation (4.13) implies the following sum rules for the pion polarizabilities whereΩ 0 (0) denotes the derivative of the Omnès factor at s = 0 and the first term in each line disappears for a pion-pole LHC. The numerical evaluation for ∆ 0,++ = N 0,++ , see Table 5, confirms the observation from Sect. 4.1 that the charged-pion amplitude is better reproduced than its neutral-pion analog. In fact, the chargedpion dipole polarizability comes out in perfect agreement with ChPT [137], as well as with the recent measurement by COMPASS (α 1 − β 1 ) π ± = 4.0(1.2) stat (1.4) syst × 10 −4 fm 3 [140]. The quadrupole polarizability is more sensitive to poorly-determined low-energy constants, but the sum-rule value lies within the range quoted in [137] and is also close to (α 2 − β 2 ) π ± = 15.3(3.7) × 10 −4 fm 5 obtained in [55] by combining the more stable chiral prediction for the neutral-pion quadrupole polarizability with a finite-matching-point sum rule for I = 2.
In contrast, both neutral-pion polarizabilities differ by about 10 units each from the full result, a deficiency that signals the impact of higher contributions to the LHC, as we will demonstrate in the following. The next such contribution is generated by the exchange of vector-meson resonances V = ρ, ω, whose impact can be roughly estimated within a narrow-width approximation. Starting from a vector-pion-photon coupling of the form 20) with coupling constant related to the partial width according to we obtain [54,141] Unfortunately, the polynomial piece ∝ s is ambiguous and would even appear with a different sign in an antisymmetric-tensor description of the vector-meson fields [54,142]. It is for this reason that in a full Roy-Steiner approach only the imaginary parts are employed, while the low-energy parameters enter via subtraction constants. However, in order to predict the numerical values of the polarizabilities in terms of the lowest contributions to the LHC in γγ → ππ we do need the full amplitude in (4.22). Parameterizing the ambiguity according to s → ξ V s, we find (4.23) Adding ρ, ω contributions using masses and partial widths from [143], the quadrupole polarizabilities are shifted by (α 2 − β 2 ) π ± V = 0.9 × 10 −4 fm 5 and (α 2 − β 2 ) π 0 V = 10.3 × 10 −4 fm 5 , which explains how vector-meson contributions can restore agreement with ChPT for the neutral pion without spoiling the charged channel. In fact, the hierarchy can be attributed almost exclusively to the large ω → π 0 γ branching fraction 24) which ensures that the same mechanism applies for the dipole polarizability as well.
In any case, such corrections are not contained in our estimate (4.17), but at least at the on-shell point the impact is expected to be moderate due to the fact that the charged-pion intermediate states are most important. In particular, the physics related to the low-energy constantsl 6 −l 5 , which appear at two-loop level in the chiral expansion for the HLbL tensor [48], only contribute to the charged-pion polarizability (a more detailed comparison to ChPT is provided in App. H). Our calculation therefore demonstrates in a model-independent way that such next-to-leading-order corrections are moderate in size, in agreement with [50], but in contradiction to the large corrections suggested in [49]. This conclusively settles the role of the charged-pion dipole polarizability in the HLbL contribution to (g − 2) µ .

Conclusions
In this paper we presented an in-depth derivation of the general formalism required for the analysis of two-pion-intermediate-state contributions to HLbL scattering in (g − 2) µ . As a first step we gained a detailed understanding of the properties of the HLbL tensor, including its decomposition into scalar functions, projection onto helicity amplitudes, and the relation between the different sets we needed to introduce in the course of our derivation, see Table 1. Some of the more subtle issues that arose in this derivation are related to the fact that, in order to write down dispersion relations for the HLbL tensor, we had to start with a redundant set of functions. At first sight, the relation between the latter and the physically observable helicity amplitudes seems to suffer from ambiguities. To show that this arbitrariness is only apparent we invoked a set of sum rules, which follow from a simple assumption on the asymptotic behavior of the HLbL tensor. These sum rules allowed us to construct a basis for kinematics with one single on-shell photon (singly-on-shell) that satisfies unsubtracted dispersion relations. In addition they lead to physically relevant sum rules that constrain the helicity amplitudes for γ * γ * → ππ. After working out the basis change from the singly-on-shell basis to helicity amplitudes, we combined this general formalism with a partial-wave expansion to address two-pion-rescattering contributions.
In a second step we thoroughly tested our formalism using the example of the pion box, whose full result is known thanks to an exact relation to the scalar QED pion loop we established earlier. In particular, we demonstrated that the sum rules that follow from our assumptions on the asymptotic behavior of the HLbL tensor are fulfilled. Moreover we studied whether the partial-wave expansion of the pion box converges to the full answer after resummation, and demonstrated that it does so sufficiently quickly. Given that the pion-box contribution can be expressed exactly in terms of the pion vector form factor-much as the HVP contribution of two pion intermediate states is completely determined by this form factor-we showed that by fitting a dispersive representation of the pion vector form factor to a combination of space-and time-like data, the space-like form factor required for the HLbL application can be constrained to a very high precision, leading to a π-box µ = −15.9(2) × 10 −11 for the pion-box contribution.
The main motivation for developing a partial-wave framework is to be able to calculate rescattering corrections, since only in a partial-wave basis for helicity amplitudes do unitarity relations become diagonal. Accordingly, as a first application of the formalism developed here we studied the unitarization of the pion box, a correction whose evaluation requires the use of partial-wave amplitudes. Concentrating on S-wave ππ-rescattering effects, we presented a first numerical estimate, which, together with the pion-box evaluation, combines to a π-box µ + a ππ,π-pole LHC µ,J=0 = −24(1) × 10 −11 (5.1) for the leading two-pion contributions to (g − 2) µ . The improvement in accuracy with respect to previous model-dependent analyses is striking. It derives: (i) from our model-independent approach based on dispersion relations that allows us to express this contribution, in a rigorous way, in terms of hadronic observables, and (ii) from the fact that all quantities needed in this calculation (the pion vector form factor and the ππ S-wave phase shifts) are very well known. Remaining two-pion contributions that have not been addressed yet are likely to lead to larger uncertainties, but given that the error quoted in (5.1) lies an order of magnitude below the experimental accuracy goal, we are confident that the final estimate for the total HLbL contribution should be sufficiently accurate to make these measurements of (g − 2) µ a sensitive test of the Standard Model. Many of the technical advances described here are not specific to the two-pion intermediate state but completely general and thus lay the groundwork for a full phenomenological analysis of HLbL scattering. Armed with these, we are now poised to study other contributions and apply further refinements to the numerical analysis of the two-pion channel and beyond. ber PIEF-GA-2013-622527 and P.S. by a grant of the Swiss National Science Foundation (Project No. P300P2_167751).
A Transformed tensor decomposition for the contribution to (g − 2) µ For the calculation of (g − 2) µ , we make a linear transformation of the BTT tensor decomposition: Only 19 of the new structuresT µνλσ i contribute to (g−2) µ , which is the minimal number of independent contributions in the (g − 2) µ kinematic limit. The symmetry under q 1 ↔ −q 2 reduces this to 12 terms in the master formula.

A.2 Scalar functions
In terms of the BTT functions Π i , the transformed scalar functionsΠ i that contribute to (g − 2) µ are defined in (2.15) and (2.16). The ones that do not contribute to (g − 2) µ are given by: C Feynman-parameter representation of the pion box In the limit q 4 → 0, the pion-box contribution to the scalar functions that appear in the master formula can be written as a two-dimensional Feynman parameter integral: The remaining functions entering the master formula can be obtained with the crossing relations (2.16).

D Scalar functions for the two-pion dispersion relations
Here, we give the explicit solution for the scalar functionsΠ i , which fulfill unsubtracted single-variable dispersion relations and only depend on physical helicity amplitudes. First, we define the following linear combinations of BTT functions: as well as Π c i := C 13 Π i . The 19 functions that contribute to (g − 2) µ can be written in the form (for q 2 4 = 0 and t = q 2 2 ) All other∆ i and∆ i are zero. We use the abbreviations q ijk := q 2 i + q 2 j − q 2 k , Σ = q 2 1 + q 2 2 + q 2 3 , and λ 123 := λ(q 2 1 , q 2 2 , q 2 3 ) for the Källén function. We define five additional scalar functionsΠ i that appear in sum rules: The singly-on-shell basis consists of 27 elements. The three functionsΠ 25 ,Π 26 , andΠ 27 are not given explicitly as they have no significance in the connection with (g − 2) µ .

E.1 Unphysical polarizations
In the following, we explain why unphysical polarizations are not trivially absent in any representation. In short, although unphysical polarizations cannot contribute to any observable, the absence of such unphysical contributions is manifest only if the basis is well chosen. Otherwise, their apparent contribution vanishes only due to the presence of sum rules for the scalar functions. Suppose we have a decomposition of the HLbL tensor into a "physical" and an "unphysical" piece, where the scalar functions Π phys i are linear combinations of helicity amplitudes with only transverse polarizations of the external photon. The scalar functions Π unph i contain also contributions from the longitudinal polarization. Because these scalar functions cannot contribute to an observable, the unphysical tensor structures have to fulfill Such structures do not contribute to (g − 2) µ , because the derivative with respect to q ρ 4 either vanishes for q 4 → 0 or is symmetric in ρ ↔ σ.
Next, we apply the following transformation, which mixes the physical and unphysical part: Because not all tensor structures have the same mass dimension, the coefficient α can be dimensionful, e.g. α = q 3 · q 4 if the mass dimension of T µνλσ b,unph is larger by two units than the one of T µνλσ a,phys , all while avoiding kinematic singularities. The new structure T µνλσ b,unph − αT µνλσ a,phys still cannot contribute to (g − 2) µ if α ∝ q 4 . However, we have introduced a new combination of unphysical and physical helicity amplitudes into the scalar coefficient functions of T µνλσ a,phys . If we make such a transformation in the discontinuity appearing in an s-channel dispersion integral, the factor α = q 3 · q 4 becomes in the (g − 2) µ limit where we have replaced the Mandelstam variable s by the integration variable of the dispersion integral s . This factor cancels with the Cauchy kernel 1/(s − q 2 3 ), producing an apparent polynomial contribution that depends on both physical and unphysical helicity amplitudes. As shown in Sect. 2.3 this polynomial contribution actually vanishes due to sum rules, but in practice it can be tedious to identify the combination of physical and unphysical helicity amplitudes that corresponds to this vanishing polynomial, and, worse, in a partial-wave expansion these sum rules are only fulfilled after resumming all partial waves. Since the above example implies that setting by hand only the unphysical polarizations to zero leads to a wrong result, a practical implementation requires a basis where this contribution is manifestly absent from the beginning. The construction of this basis is performed in Sect. 2.4.1.

E.2 Comparison to forward-scattering sum rules
In [63], sum rules have been derived for the case of forward HLbL scattering. In the following, we compare them to our fixed-t sum rules derived in Sect. 2.4.3. To this end, we consider the case of general forward kinematics, i.e. q 3 = −q 1 , q 4 = q 2 , (E. 5) which implies for the Lorentz invariants t = 0, u = 2q 2 1 + 2q 2 2 − s, q 2 3 = q 2 1 , q 2 4 = q 2 2 .
(E.6) external momenta with tensor integrals that depend only on a single momentum and can be solved by a standard tensor decomposition. Unfortunately, with this method the expressions become very large, which makes the computation already at the level of D-waves extremely inefficient. In order to avoid the excessive amount of contractions in the calculation of the phase-space integral, here we present an alternative way to calculate the tensor integrals directly by taking derivatives of scalar integrals. This allows us to calculate even the G-wave unitarity relation, as required for the present application in Sect. 2.5. In the case of D-waves, we have checked that both methods give the same result.

F.2 Direct matrix inversion
The expressions for the helicity amplitudes in terms of the scalar coefficient functions in the tensor decomposition are easily obtained by contracting the HLbL tensor with the polarization vectors. Expressing the scalar functions in terms of the helicity amplitudes requires the inversion of these relations. If we consider the singly-on-shell case, this amounts to the inversion of a 27 × 27 matrix. The direct analytic inversion of a general matrix of this size is not possible, but in this case it can be reconstructed along the following lines. Let us define the basis change from the singly-on-shell helicity amplitudes to scalar functions as where η is a 27 × 27 matrix. Its inverse is effectively the matrixč in (2.84) (restricted to λ 4 = 0) that we need to determine in order to obtain the imaginary parts of the scalar functions through unitarity. First, we note that the basis of helicity amplitudes suffers from the presence of kinematic singularities, which makes the expressions for η more involved. These singularities can be removed by applying the general recipe of [62] for the construction of amplitudes free of kinematic singularities: first, the singularities at the boundary of the physical region can be removed bŷ where m 1 = λ 1 − λ 2 , m 2 = λ 3 − λ 4 , and z is the cosine of the scattering angle. In our case of fixed-t singly-on-shell kinematics, we have z = . Next, the parity-conserving amplitudeŝ H j ±Ĥj (F. 16) are formed (see Sect. 2.5.1 for the notation). Finally, the remaining singularities can be removed by multiplying with the appropriate powers of √ s, λ 1/2 12 (s), and λ 1/2 34 (s), see [62]. We note that the Martin-Spearman amplitudes constructed in this way are free of kinematic singularities, but have an asymptotic behavior that is much worse than the one of the BTT scalar functions. Since all square-root singularities have been removed, the basis change from the scalar functionš Π i to the Martin-Spearman amplitudes is now meromorphic in s, q 2 1 , q 2 2 , and q 2 3 . We determine all matrix entries with partly numerical methods as follows.
Numerically, the inversion of the 27 × 27 matrix is straightforward. The denominators of the meromorphic matrix entries can be guessed from the pole structure of the numerical inversion: they are products of simple polynomials such as λ 123 , λ 12 (s), (q 2 1 − q 2 2 + q 2 3 ) etc. We calculate numerically the matrix inversion as a function of each of the Lorentz invariants in turn, keeping the other three invariants fixed. A plot of the matrix entries as a function of the varying variable reveals the poles and therefore the exact form of the denominators. This simple but tedious task has to be performed for all 27 × 27 entries. The remaining numerators are then polynomials of the form i+j+k+l=n a ijkl s i (q 2 1 ) j (q 2 2 ) k (q 2 3 ) l , (F. 17) where the mass dimension of the numerator is 2n and known beforehand. In most cases, n is a small number, although for very few entries we encounter a maximal value of n = 9, which results in a polynomial with 220 terms. We perform the numerical inversion on a grid consisting of 9 4 points in the four-dimensional space of s, q 2 1 , q 2 2 , and q 2 3 and determine the integer coefficients a ijkl for each of the numerators of the 27 × 27 matrix entries by a fit. In contrast to the determination of the denominators by hand, this fit of the numerators can be easily automatized.
Combining the results with the (simple) basis change from helicity to Martin-Spearman amplitudes then leads to the full analytic expression for the inverted basis changeč. In particular, it is straightforward to check analytically that the matrixč determined partly with numerical methods is indeed the exact inverse of η. The result is provided as supplementary material in the form of a Mathematica notebook.
G Partial-wave expansion of the γ * γ * → ππ pion-pole contribution In order to test the partial-wave formalism, we expand the pion-pole contribution to γ * γ * → ππ into partial waves. The scalar functions are given by [31] (with isospin conventions from [28]): The helicity amplitudes become: (G.2) We calculate the partial-wave expansion thereof: H Pion polarizability and γγ → ππ in ChPT that exploits a sum rule for the relevant low-energy parameters, largely saturating the phenomenological value. In particular, in this framework the impact of higher contributions such as the a 1 to the LHC on the polarizability itself is expected to be small-in fact, the exchange of vector mesons would contribute first-so that significant corrections would require a weighting in the g − 2 integral that emphasizes kinematics away from s = 0, where the polarizabilities are defined. Such contributions cannot a priori be excluded, all the more since the comparison with the physical polarizability only determines the on-shell properties, but not the dependence on the photon virtualities.