Dispersion relations for $\gamma^*\gamma^*\to\pi\pi$: helicity amplitudes, subtractions, and anomalous thresholds

We present a comprehensive analysis of the dispersion relations for the doubly-virtual process $\gamma^*\gamma^*\to\pi\pi$. Starting from the Bardeen-Tung-Tarrach amplitudes, we first derive the kernel functions that define the system of Roy-Steiner equations for the partial-wave helicity amplitudes. We then formulate the solution of these partial-wave dispersion relations in terms of Omn\`es functions, with special attention paid to the role of subtraction constants as critical for the application to hadronic light-by-light scattering. In particular, we explain for the first time why for some amplitudes the standard Muskhelishvili-Omn\`es solution applies, while for others a modified approach based on their left-hand cut is required unless subtractions are introduced. In the doubly-virtual case, the analytic structure of the vector-resonance partial waves then gives rise to anomalous thresholds, even for space-like virtualities. We develop a strategy to account for these effects in the numerical solution, illustrated in terms of the $D$-waves in $\gamma^*\gamma^*\to\pi\pi$, which allows us to predict the doubly-virtual responses of the $f_2(1270)$ resonance. In general, our results form the basis for the incorporation of two-meson intermediate states into hadronic light-by-light scattering beyond the $S$-wave contribution.


Introduction
Apart from the two-photon decay of the neutral pion, the reaction γγ → ππ constitutes the simplest process that gives access to the electromagnetic properties of the pion, most notably its dipole polarizabilities [1]. Experimentally, most information on the scattering process comes from e + e − colliders via the reaction e + e − → e + e − ππ [2][3][4][5][6][7], while the kinematics relevant for the extraction of the polarizabilities is more directly probed in the Primakoff process, where an incident pion scatters of the Coulomb field of a heavy nucleus and produces a final-state photon-pion pair [8][9][10]. The measurement of the polarizabilities, as well as the energy dependence of the γγ → ππ cross section, also provides a key test of chiral perturbation theory (ChPT) [11][12][13][14], not only because it is the simplest electromagnetic scattering process involving hadrons, but also due to the sensitivity to chiral loop corrections, see [15,16] and [17][18][19][20][21] for the one-and two-loop calculation, respectively. While an extraction of the charged-pion polarizability via radiative pion production off the nucleon [22] had been interpreted as a potential tension with ChPT [21]-despite the model-dependence from the extrapolation to the pion pole-the most recent Primakoff measurement [10] confirmed the chiral prediction. In addition to a future update from COMPASS [23], further low-energy measurements that would entail additional information on the charged-pion polarizabilities are planned at Hall D at Jefferson Lab via the Primakoff process with an incident photon [24].
First of all, the ChPT amplitudes for the doubly-(or even singly-) virtual case have only been worked out at one-loop order [45]. Since the one-loop contribution does not display any angular dependence except for the charged-pion Born terms, this implies that chiral predictions for D-and higher partial waves are not available. Second, it was shown in [35] that an adequate description of the D-waves requires the inclusion of vector mesons in the left-hand cut (LHC) of the γγ → ππ amplitudes, most efficiently in terms of the LHC of the partial waves. This strategy extends the standard Muskhelishvili-Omnès (MO) solution [49,50], and, as we will show here, its necessity is related to the high-energy behavior of the vector-meson partial waves and thus potential subtractions in the MO solution. Third, the derivation of partial-wave DRs for the helicity amplitudes has to be based on scalar functions that avoid kinematic singularities and zeros [51,52]. Based on the corresponding set of amplitudes from [41] we explicitly write down the kernel functions that couple the various partial waves and perform a basis change that diagonalizes their MO solution. Finally, we observe that if vector resonances are to be included in the MO solution in terms of the LHC also in the doubly-virtual case, the analytic structure of these amplitude complicates the implementation for sufficiently large (space-like) virtualities. While the occurrence of anomalous thresholds [53] is expected in the time-like  Figure 1: γ * γ * → ππ as a sub-process of e + e − → e + e − ππ [41].
regime [41,44,54], the analytic structure of the resonance LHCs is sufficiently complicated that even in the space-like case a deformation of the integration contour becomes unavoidable.
We first recall the definition of helicity amplitudes, partial waves, and Bardeen-Tung-Tarrach (BTT) invariant functions in Sect. 2, based on which we then derive all relevant kernel functions that define the full system of Roy-Steiner (RS) equations for the partial-wave helicity amplitudes. In Sect. 3 we then write down the MO solution of these equations, and discuss in detail the role of subtraction constants as well as the analytic structure of the resonance partial waves. Some numerical results will be presented in Sect. 4, before we conclude in Sect. 5 and comment on the implications of our results for the application to HLbL scattering.

Helicity amplitudes
We largely follow the conventions of [41], but for completeness repeat the basic definitions. The process γ * γ * → ππ is strictly speaking not observable, but derived from processes with well-defined asymptotic states under certain assumptions. We start from e + (k 1 )e − (k 2 ) → e + (k 3 )e − (k 4 )γ * (q 1 )γ * (q 2 ) → e + (k 3 )e − (k 4 )π a (p 1 )π b (p 2 ), (2.1) shown in Fig. 1, with isospin labels a and b for the pion states and momenta as indicated. At O(e 4 ), the amplitude for this process is given by where ξ is a gauge parameter for the photon propagators and the tensor is defined in pure QCD. The contraction of this tensor with appropriate polarization vectors is then identified as an amplitude for the off-shell process where λ 1,2 denote the helicities of the photons. The connected part is given by and the contraction with polarization vectors finally defines the helicity amplitudes according to Next, the kinematic invariants are taken as 1 For the helicity amplitudes it is also convenient to choose a frame, we construct the helicity amplitudes with the momenta and polarization vectors in the s-channel center-of-mass system. This gives (2.10) and we introduced the notation the s-channel scattering angle , (2.12) as well as the polarization vectors For the particular choice of ξ i = q 2 i the longitudinal states are normalized to 1, but, since the off-shell photons are not physical states, the choice of ξ i cannot affect physical observables, a useful check on the calculation. For convenience, we also define helicity amplitudes that stay finite in the limit q 2 i → 0 and introduce the labelinḡ Finally, the helicity amplitudes are expanded into partial waves according to [59] where m = |λ 1 − λ 2 |.

Tensor decomposition and dispersion relations
DRs should not be derived for the helicity amplitudes directly due to their complicated analytic structure, but instead for scalar functions that are free of kinematic singularities and zeros. Such a basis has been derived in [41] following the general recipe established in [51]. In contrast to the singly-virtual case, however, the doubly-virtual process is sufficiently complicated that an additional limitation first observed in nucleon Compton scattering [52] occurs, i.e. that to cover all kinematic limits a sixth Lorentz structure needs to be provided, in addition to the five expected in correspondence to the five independent helicity amplitudes. Fortunately, the number of required scalar functions can be reduced by using crossing symmetry in the pion system, again in analogy to nucleon Compton scattering [60], finally leading to the scalar functions A i as defined in [41]. Explicitly, we have where q 3 = p 2 − p 1 and In terms of these functions, the helicity amplitudes read as follows where λ 12 (s) = λ(s, q 2 1 , q 2 2 ). In this paper, we will make repeated reference to the expressions that follow from the pion-pole terms as well as the tree-level exchange of vector mesons, with partial waves N J,i (s) and h V J,i (s), see App. A. The form of the DRs for the coefficient functions A i is defined by a second constraint on the Mandelstam variables besides the on-shell condition (2.8). As argued in [55], the optimal choice for a process with crossing properties of γ * γ * → ππ is given by hyperbolic DRs (HDRs), for which the dispersive variables are constrained to lie on hyperbolas of the form which implies the relation for the differentials and for the scattering angles, leading to with Born terms A π i given in App. A. In writing (2.24) we have implicitly assumed that the Bornsubtracted amplitudes fulfill unsubtracted HDRs. We stress that this is not equivalent to assuming unsubtracted HDRs for the full amplitude, since for A 1 the Born term itself goes asymptotically to a constant along the hyperbola. Assuming unsubtracted HDRs for the full amplitude would thus imply cancellations of this constant behavior with a contribution from heavier intermediate states.
A realistic description of γ * γ * → ππ beyond the S-waves requires further contributions to the LHC, most importantly the exchange of vector mesons, see App. A for the explicit expressions. This Lagrangian-based representation suffers from a polynomial ambiguity [35]: choosing a different Lagrangian representation alters the real part of the amplitude, while the residues of the vector-resonance poles are free from such ambiguities. In the narrow-width limit, in which the imaginary parts of the vector-meson exchange collapse to δ-functions, this can be demonstrated by comparing the expression resulting from the HDRs (2.24) with the starting point (A.6). We find the differences The remaining ambiguity maps onto the polynomial obtained when changing the representation of the vector mesons from vector to antisymmetric tensor fields [35,61,62], and due to (2.20) only affects the S-waves. In the following, we will indeed define the resonance LHCs by their a → ∞ limit, which corresponds to a fixed-s DR, because this is the situation encountered in a dispersive approach to HLbL scattering where the different topologies are defined using the Mandelstam representation. Moreover, we will not comment further on the S-wave case-there, the consideration of subtractions is unavoidable to capture model-independently the effect of vector resonances in the LHC-but concentrate on how to extend the dispersive description to D-waves. The basic idea in the derivation of RS equations is then as follows: expand the imaginary parts of (2.24) into partial waves, express the internal angle z in terms of the external angle z by means of the hyperbola conditions (2.21), and project the whole system onto partial waves. In contrast to [36] we will not calculate the kernel functions for the LHC explicitly, but directly work with a narrow-width approximation for the resonances. The resulting system of partial-wave DRs then takes the form with s-channel kernel functions K ij JJ (s, s ). The calculation of these kernel functions is straightforward, but reveals ostensible singularities in 1/s as well as factors involving √ s that originate from the definition of the helicity amplitudes. Before turning to the explicit form of the kernel functions, we therefore first study singularities that may be produced by the partial-wave expansion.

Kinematic singularities and partial-wave expansion
Using the recipe of [63], the kinematic singularities in the helicity amplitudes can be separated according toH and indeed theH i are given by a sum of the A i with coefficient functions that are polynomials in s, t, u. However, we are mainly interested in the kinematic singularities of the partial waves, not the full amplitudes. To derive the corresponding singularities-with critical points s = 0, s = 4M 2 π , and the zeros of λ 12 (s)-let us assume that the scalar functions A i fulfill an unsubtracted fixed-s DR By performing the angular integrals in terms of Legendre functions of the second kind, in analogy to the partial-wave projection in App. A, this leads to the representation . (2.31) From these relations we may read off the kinematic singularities as follows: at threshold the combinatioñ

For the RS system (2.27), the central conclusion of this derivation is that upon the rescaling
of the S-and D-waves, theh J,i (s) do not have further kinematic singularities provided that the full amplitudes A i (s, t) satisfy unsubtracted fixed-s DRs. This situation changes once subtractions are introduced in the fixed-s DR. In this case, both the subtraction polynomial and the dispersion integral display 1/s singularities, whose residues will cancel each other if sum rules exist that reinstate the unsubtracted version. In the derivation of RS equations the amplitudes are expanded into partial waves both at the level of the hyperbolic dispersion integrals as well as in the integrands. This implies that the full amplitudes are approximated by a truncated partial-wave series, which spoils the asymptotic behavior in the crossed channel, so that unsubtracted fixed-s DRs are no longer possible. Therefore, additional 1/s singularities may appear at any finite order in the partial-wave expansion, and these are precisely the singularities observed in the RS kernels in (2.27).

Kernel functions
Motivated by the discussion in the preceding section, we write the RS kernels not for the h J,i , but for theh J,i according to (2.35), replacing (2.27) bỹ The S-wave kernel functions recover the corresponding results from [41] K 11 00 (s, s ) = K 55 00 (s, s ) = as do the diagonal D-wave kernels. The full list of non-vanishing kernel functions reads There is no coupling of D-to S-waves, i.e. K ij 20 (s, s ) = 0, but the S-waves couple to the D-waves through the kernels K ij 02 (s, s ), which depend linearly on the hyperbola parameter a. In the basis of thẽ h J,i , these kernel functions are lengthy and therefore not reproduced here. We will give them in a more convenient basis in Sect. 2.5 and their role in the MO solution will be studied in detail in Sect. 3.2.
We remark that this pattern seems to persist at higher orders in the partial-wave expansion: we have calculated all the kernel functions for J, J ≤ 8 and found that K ij JJ (s, s ) = 0 for J < J . The kernels with J = J do not depend on the hyperbola parameter a, while the J > J kernels are polynomials in a of the order (J − J)/2. This behavior was indeed obtained in [57] and it should be possible to prove the same here with similar methods.

Diagonalization of the kernel functions
The DRs for the helicity partial waves that follow from the RS system are integral equations that relate the helicity partial waves with their imaginary parts. The equations have the form of an inhomogeneous Omnès problem [49,50] and we will discuss the solution in Sect. 3. The MO solution is most easily found by performing another change of basis that diagonalizes the system of equations for a given angular momentum J and decouples the system into a set of independent equations in standard MO form. The S-and D-wave basis changes are given by where s ± = q 2 1 ± q 2 2 2 . In the on-shell or singly-virtual case, the poles in the kinematic prefactors get canceled by the soft-photon zeros at s = s ± = q 2 [65]. By writing the basis change in matrix form we find that the kernels with J = J are diagonalized to Cauchy kernels: hence the new functions fulfill DRs in standard MO form with inhomogeneities∆ J,i that contain the LHCs and the couplings due to the off-diagonal kernels: In the new basis, the kernels that couple the S-to the D-waves turn out to be very compact: where a is the hyperbola parameter and the kernel functions not listed explicitly vanish.
As an alternative to (2.42), the resonance LHC can be written in terms of a dispersion integral over its discontinuity, see Sect. 3.3 for details. This results in a representatioň

Asymptotic behavior and sum rules
The DRs (2.42) are a direct consequence of the HDRs (2.24), following upon partial-wave projection and the basis change (2.39). They can be written without any subtractions provided that the initial HDRs are unsubtracted. The pure rescattering contributionš are functions that contain only the right-hand unitarity cut. If we make the additional assumption that not onlyȟ resc J,i (s) but sȟ resc J,i (s) vanishes for s → ∞, then the DR holds, which in turn implies the sum rules These sum rules are essential to justify unsubtracted Omnès representations in Sect. 3, however they need to be validated. Based on the general consideration of unsubtracted fixed-s DRs for the scalar functions (2.30), we expect an asymptotic behavior of the D-wave rescattering contribution of which implies two additional sum rules For the S-waves, we would expect a behavior s −1 log s, which in general would require a subtraction in the Omnès representation, in line with the discussion in Sect. 2.2. Due to these sum rules, most of the contributions from the off-diagonal kernels in (2.43) in fact vanish and only the simplified kernelš need to be taken into account. In particular, the dependence on the hyperbola parameter drops out, as has to happen to avoid an unphysical dependence on a. In cases wherě vanishes, no couplings of S-to D-waves would survive at all. As a special case we may consider the LHC resonance partial waves, for which the sum rules

MO solution: S-waves
Since the functions defined in (2.39) fulfill Watson's theorem [66] with ππ phase shifts δ J (s), the solution to the MO problem can be given immediately in terms of the Omnès functions .
We start the discussion by considering the restricted S-wave system [42,43]. The MO solution has the formȟ provided that (ȟ 0,i (s) −∆ 0,i (s))/Ω 0 (s) tends to zero for s → ∞. For a phase shift reaching asymptotically δ 0 (s) π, the Omnès function behaves as Ω 0 (s) s −1 , i.e. the sum rules (2.47) are employed to write the MO solution without subtractions. Performing the basis change back to the original helicity amplitudes leads to

(3.4)
In [42,43] this solution was evaluated using the pion Born terms as LHCs and a ππ phase shift that cuts off the f 0 (980) and thus the coupling to the KK channel. Phenomenologically, the pion-pole LHC produces the polarizabilities [42,43] (α 1 − β 1 ) π ± ,π-pole LHC = (5.4 . . . 5.8) × 10 −4 fm 3 , for the charged pion in perfect agreement with the chiral 2-loop prediction 5.7(1.0) × 10 −4 fm 3 [21] as well as the COMPASS measurement 4.0(1.2) stat (1.4) syst × 10 −4 fm 3 [10]. For the neutral pion the chiral prediction −1.9(2) × 10 −4 fm 3 [20] is much smaller, a discrepancy explained by the fact that the neutral channel is much stronger affected by the contribution from vector-meson exchange [67] Γ Therefore, the discrepancy does not necessarily point at a violation of the sum rule (2.47), but rather the approximation of the LHC by the pion Born term only. However, as argued in Sect. 2.2, to get the phenomenology of the neutral-pion dipole polarizabilities right the introduction of subtractions is nevertheless unavoidable, otherwise the vector-meson LHC remains ambiguous. For the quadrupole polarizabilities the vector-meson contribution indeed restores agreement with ChPT even for the neutral pion [43]. For the D-waves, polynomial ambiguities in the vector-meson LHCs do not occur, so that unsubtracted DRs in principle become possible. To this end, a modified MO solution was derived in [35] in which the vector mesons are not included via the inhomogeneities, but directly in terms of their partial waves. This corresponds to the MO solution of the DR (2.44), which in the S-wave case leads to a modification of (3.4) according to where the new integrals extend over the LHC, see Sect. 3.3 for details. Although a subtracted DR was used for the numerical analysis in [35], it was shown that sum rules that would establish an unsubtracted version are nearly fulfilled, indicating that an approximate description should be possible based on an unsubtracted system as well. Surprisingly, the same observation does not hold for the D-wave analog of (3.4), where the vectormeson LHC is treated as part of the inhomogeneity, so before returning to the D-wave MO solutions we first clarify the reason why the two strategies lead to different results. The assumption that an unsubtracted MO solution can be used relies in the standard Omnès representation on the assumption of an asymptotic behavior If the Omnès function behaves as Ω(s) s −1 , these two assumptions are equivalent provided that the vector-meson LHC vanishes asymptotically at least as h V (s) s −2 . According to (2.54), this is not the case for the S-waves and the two D-wavesȟ V,fixed-s 2,i , i = 4, 5. The difference between the two representations is proportional to the Omnès function: Note that the somewhat pathological high-energy behavior of the vector-meson LHC is only present in the real part, while the imaginary part is much better behaved. The MO solution in the standard form keeps the bad high-energy behavior of the resonance LHC, while the modified representation only involves the imaginary part of the vector-meson LHC and imposes a better high-energy behavior on the Born-subtracted partial waves. Indeed, it turns out that the unsubtracted standard MO form does not reproduce the peak of the narrow f 2 (1270) D-wave resonance. This can be understood by noting that the modified MO solution is equivalent to the standard form with a subtraction, where the subtraction constant is effectively calculated in terms of the vector-meson LHC. In this case, the resonance peak is fully described by the subtraction term: let us consider the part without subtraction constant, and let us further consider the simple case in which ∆(s) = α/s. In this case, the dispersive integral can be performed analytically by using the spectral representation of the inverse Omnès function The result is proportional to the Omnès function, as expected, but one finds an additional polynomial whose coefficients are determined by normalization and derivative of the Omnès function at s = 0.
For a narrow resonance with mass M R , as the f 2 (1270) in the D-wave, one has Ω(s) ∼ M 2 R /(M 2 R − s), so that 1 − sΩ(0) vanishes at s = M 2 R . The resonance peak is thus described exclusively by the subtraction term that we dropped in (3.11). Such a situation indeed occurs for some of the D-waves, which is why we develop the formalism to include the vector mesons in the LHC in the modified MO representation as in (3.7) even in the doubly-virtual case, to be able to put forward an unsubtracted DR in the description of the f 2 (1270) resonance in γ * γ * → ππ, see Sect. 3.3.

MO solution: D-waves
With the diagonalization of the D-wave system derived in Sect. 2.5, the MO solutions follow immediately, as the defining DRs are given in decoupled form (2.42)
α + 15 (s) vanishes for s + = s − = 0, so that the onshell process remains unaffected. As argued in Sect. 2.3, the appearance of the 1/s singularities is an artefact of the partial-wave expansion, and accordingly the size ofα ± 15 (s) should be in line with effects expected from higher partial waves to allow for a cancellation in the full amplitude. We checked numerically that this residual coupling between S-and D-waves is indeed small, but this conclusion remains to be verified after integration over the weight functions in the g − 2 integral.

Analytic structure of the resonance partial waves
To include the vector mesons directly in terms of the LHCs of their helicity partial waves, as in (3.7), we need to analyze the analytic structure of their LHCs in more detail. These LHCs are produced by the t-and u-channel exchange of a resonance with mass M V . For the cut structure itself the details of the partial-wave projection, i.e. angular momentum and helicity states, are irrelevant, let us therefore write symbolically where Q(x) is a Legendre function of the second kind with a cut in the complex x-plane between ±1 and x V is given in (A.9). Instead of considering the Legendre function, we can also study the angular integration path in the complex t-or u-plane, with endpoints at In this way, the wrapping of the integration contour around the pole at t, u = M 2 V determines possible singularities. Throughout, we will restrict the analysis to space-like virtualities, q 2 i < 0, as required for HLbL scattering. For time-like virtualities anomalous thresholds are certain to appear in any dispersive representation, even in MO solutions in the standard form (3.4), see [44].
First, possible kinematic square-root singularities at s = 4M 2 π , s = 0, and s = −( −q 2 1 ± −q 2 2 ) 2 = s ± are in fact absent: in the explicit representation, these singularities of the Legendre function are balanced by the kinematic prefactor. Equivalently, in the path-deformation approach, they are lifted by the same factors coming from the Jacobian when switching from z to t or u as integration variables. The only singularities are therefore logarithmic branch points at x = ±1 or, equivalently t ± = M 2 V , given by The other ends of the branch cuts are located at s = 0 and s = −∞, respectively, as can be inferred from the replacement M 2 V → ∞. For q 2 i → 0, the two branch points are at is present, while the other one disappears. By writing the Källén function as (3.28) it follows that the square roots in s ± cut can only produce an imaginary part for 29) i.e. for time-like virtualities. Therefore, for space-like virtualities the cut structure seems to remain simple: one expects just two branch cuts on the negative real axis, one from −∞ to s − cut , the other from s + cut to 0. However, an important subtlety arises that is reminiscent of anomalous thresholds in triangle diagrams, which appear for sufficiently large time-like virtualities: there, the discontinuity itself has singularities that depend on the virtualities and cross the unitarity cut of the triangle diagram. By entering the physical sheet, they require a deformation of the integration contour and add an "anomalous" discontinuity [41,44,54]. Here, the discontinuity itself has the two square-root branch cuts from the kinematic factors in (3.24), i.e. cuts for s ∈ [0, 4M 2 π ] and s ∈ [s + , s − ]. For space-like q 2 i , the second cut in the discontinuity lies between the two LHCs of the partial waves, i.e. s ± ∈ [s − cut , s + cut ], where the points s + and s − cut coincide for This condition can be fulfilled even for space-like virtualities, so that the corresponding points deserve special attention. Let us consider the difference and add a small imaginary part to the virtualities, q 2 i → q 2 i + i . Then, for fixed values of q 2 1 we trace the path of ∆ cut as a function of q 2 2 . We find that : ∆ cut has a positive real part and a small positive imaginary part of order .
: the imaginary part vanishes and the real part is negative (of order 2 ).
: the real part becomes again positive, the imaginary part is negative.
This implies that the square-root singularity of the discontinuity of the partial waves, which for q 2 1 q 2 2 < (M 2 V − M 2 π ) 2 lies on the second sheet of the logarithmic LHCs, moves onto the physical sheet for q 2 1 q 2 2 > (M 2 V − M 2 π ) 2 , see the sketch in Fig. 2. This requires a deformation of the left-hand integration contour. In the case q 2 1 q 2 2 < (M 2 V − M 2 π ) 2 , the left-hand integral consists of two integrals

dispersion integral picks up an anomalous contribution
In the case of the Legendre functions, the normal imaginary part is given by Since the anomalous singularity is a square-root branch cut, the anomalous discontinuity is simply which again can be verified by considering the path deformation in the complex t-plane. The representation (3.33) indeed displays the correct integration regions, including anomalous contributions, but does not yet fully cover the realistic case encountered in γ * γ * → ππ, in which the singularities are stronger than in the schematic example discussed above. In the case of higher partial waves, the discontinuity at the anomalous singularity behaves as (s − s + ) −(J+1)/2 due to Q J (x V ) alone. However, additional kinematic prefactors appear both in the partial waves and in the kernel functions, so that in the realistic cases the anomalous singularity actually scales as (s − s + ) −5/2 for the S-waves, as (s − s + ) −7/2 forȟ 2,1 andȟ ± 45 , and as (s − s + ) −9/2 forȟ ± 23 . Clearly, this is not integrable and the above representation (3.33) has to be modified further.
To resolve this apparent contradiction, the important observation is that the contour integral around the anomalous singularity gives a non-vanishing contribution, so that the full anomalous integral is finite. The total anomalous integral can then be calculated as follows. We write where the first term collects the singular pieces of the integrand and∆(s, s ) vanishes as a square root for s → s + . The coefficients a k (s) can be calculated analytically. The anomalous integral splits into two pieces (3.37) The integral around the singularity cancels exactly the singular pieces of the integral along the real axis, as can be seen by splitting the path into an integral up to s + − and a circular integral around the singularity with radius .
Finally, if the integral over∆(s, s ) is calculated numerically, one faces the problem of numerical instabilities close to s + , given that∆(s, s ) is defined by the difference of two divergent expressions. This numerical issue can be handled by replacing∆(s, s ) close to s + by a fit function that has the same square-root-like behavior, i.e.
with some appropriate power n. The coefficients b k (s) are determined by a fit to∆(s, s ) in the vicinity of s = s + , but outside the region where numerical instabilities occur. The size of this region depends on the values of the virtualities, so that the fit region needs to be adapted accordingly.
We first verified that with this strategy we can indeed recover the original resonance partial waves from a representation such as (3.33), even for large space-like virtualities that exceed the critical point (3.30). The generalization to the unitarized case with Omnès functions as in (3.7) proceeds along the same lines, given that the Omnès functions do not alter the singularity structure, see App. B for more details. In this case, however, the derivatives of the Omnès function need to be provided as well, which in a numerically stable way follow from the spectral representation Im Ω J (s ) (s − s) n+1 , (3.39) or directly by taking derivatives of (3.2). With increasing degree of singularity, numerical stability of the extrapolation becomes more of an issue, but even for the −9/2 singularities ofȟ ± 23 remains under good control as long as the fit region is chosen prudently. However, we stress that for all D-waves besidesȟ ± 45 the standard MO representation still applies, which does not involve integrals over the LHC. We verified that forȟ 2,1 andȟ ± 23 the above recipe for the treatment of the anomalous threshold in the modified MO representation indeed reproduces the same result as the standard MO representation.

Numerics
In this section we present some numerical applications of the formalism developed in Sect. 3, mainly focused on the contribution of the f 2 (1270) resonance to the various helicity amplitudes. Experimentally, there is ample information on the on-shell cross section γγ → π + π − , π 0 π 0 , derived from e + e − → e + e − ππ via suitable cuts on the lepton momenta.

On-shell case
In the on-shell case only the helicity amplitudes H ++ and H +− contribute. Adjusting the flux factor to an actual γγ initial state, one has where the particle-basis amplitudes are related to the isospin ones by the rotation given in (A.2).
To illustrate the behavior of the f 2 (1270), an isospin-0 D-wave resonance, we neglect unitarity corrections in the isospin-2 partial waves and combine our results for the D-waves with the S-waves from [42,43]. The only free parameters are then the photon couplings of the vector resonances C V , which in a narrow-width picture are related to the partial widths by means of (A.7). We find that the physical couplings do not exactly reproduce the cross section. This observation corresponds to the fact  Figure 3: Cross section for γγ → π + π − (left) and γγ → π 0 π 0 (right), in comparison to the data from Belle [5,7], Mark II [3], CELLO [4], and Crystal Ball [2]. The lines indicate the pion Born terms (blue dashed), their unitarization (red dot-dashed), and the full solution (black solid).
that the sum rules for the subtraction constants introduced in [35] are not fulfilled exactly, pointing to a small correction from higher intermediate states not explicitly included in the calculation. 4 To ensure agreement with the measured cross section, we therefore allow the couplings to vary, as a means to include phenomenologically the effect of higher intermediate states.
Note that the experimental cross sections are not integrated over the full angular range, with | cos θ| ≤ 0.6 and | cos θ| ≤ 0.8 for the charged and neutral channels, respectively. The results in Fig. 3 follow this convention. The relevant helicity amplitudes in the on-shell case are In the figure, the blue dashed lines indicate the pion Born terms and the red dot-dashed ones their unitarization. The S-waves are treated as in [42,43], with a phase shift from the inverse-amplitude method as specified in [68]. This phase shift agrees well with dispersive ππ phase shift analyses [69][70][71] at low-energies, but removes the f 0 (980) contribution in a controlled manner, which otherwise would require a coupled-channel treatment of the ππ/KK S-waves. Further, we do not include the S-waves resulting from the vector-meson exchanges, given that these contributions are not relevant for the f 2 (1270) and need to be studied together with the pion polarizabilities to ensure the correct low-energy properties of the γγ → ππ reaction, see [42,43]. These details can be improved most conveniently by introducing subtractions in the S-wave dispersion relations, but instead we focus here on the f 2 (1270) resonance, as it emerges mainly from the unitarization of the vector-meson D-waves, see Fig. 3, using the phase shift from [70]. In the neutral channel, the unitarization of the Born terms alone actually results in a small resonant contribution, while in the charged channel it displays the pathological behavior of the standard MO solution illustrated in (3.13). In both cases the physical couplings need to be reduced by about 13% to match the physical cross section, reflecting the impact of higher LHCs beyond the lightest vector mesons ω, ρ ± , ρ 0 and potentially inelastic effects in the ππ D-wave.

Singly-and doubly-virtual case
Given that the f 2 (1270) resonance in the on-shell process can be largely understood as a unitarization of the vector mesons in the LHC, the only additional information required for the virtual processes concerns the V π transition form factors as introduced in App. A, in analogy to the pion vector form factor for the pion-pole terms. For the ω, this form factor is again available from a detailed dispersive analysis [72][73][74]. In contrast, a rigorous implementation of the ρ should proceed in terms of 2π intermediate states, based on a suitable γ * → 3π amplitude [47,48,[75][76][77]. Here, we illustrate the numerical solution by approximating the dependence on the photon virtuality by a vector-mesondominance (VMD) suppression M 2 V /(M 2 V − q 2 ), which in the case of the pion form factor reproduces the full solution very accurately [42,43]. For the ω transition form factor the deviations from VMD are more sizable, but a refined analysis should address the ρ LHC at the same time.
For the virtual processes the canonical generalization of (4.1) would be but we stress that these cross sections are not actual observables, only the e + e − cross section is, which can be seen from the fact that in the definition of (4.3) we needed to choose a convention for the flux factor and the counting of the polarization states (due to the latter this choice is discontinuous in the limit q 2 i → 0). For these reasons we present our results instead directly in terms of the squared moduli of the helicity partial-wave amplitudes, which are also the most relevant objects for the future   application to HLbL scattering. For convenience, we combine S-and D-waves into Moreover, we focus on the I = 0 amplitudes, where the unitarization effects that produce the f 2 (1270) occur.
The results are shown in Fig. 4 for several singly-virtual cases and in Fig. 5 for doubly-virtual ones. Already for the on-shell case the Born terms appear suppressed compared to the f 2 (1270) peak, due to their enhancement in the cross section by the flux factor 1/s, and that relative size does not change much once the virtualities are increased. In all cases, the H +− helicity amplitude gives the dominant effect, but the other helicity projections become increasingly important for larger virtualities. In addition, the overall size of the contribution decreases rapidly, as expected from the form factor suppression of the vector-meson couplings.

Conclusions and outlook
Dispersion relations for processes involving virtual photons require a careful study of the (helicity) amplitudes to ensure that results are not invalidated by kinematic singularities. A suitable such decomposition for the γ * γ * → ππ amplitudes has been derived before as a precursor to HLbL scattering [41][42][43], with first numerical solutions provided for the S-waves of the process. In this paper, we extended the solution to higher partial waves, introducing a new basis in which that solution takes a simple form. In particular, we studied the role of vector mesons in the left-hand cut of the amplitudes, in terms of which the D-wave resonance f 2 (1270) can be understood as an effect of the ππ final-state rescattering.
Phenomenologically, the D-waves of γ * γ * → ππ are indeed expected to contribute to HLbL scattering mainly via the f 2 (1270) resonance. For this application it is therefore crucial to understand all helicity amplitudes of γ * γ * → ππ including the potential role of subtraction constants. Here, we settled this issue conclusively, detailing how the high-energy behavior of a given partial wave is tied to the necessity of subtractions in particular variants of the Muskhelishvili-Omnès solution, which explains why for the helicity amplitudes most relevant for the f 2 (1270) the standard variant fails, but a description in terms of the left-hand singularities of the vector-meson amplitudes still applies. We then developed a strategy how to cope with the anomalous thresholds that appear in the doubly-virtual case, even for space-like virtualities, and presented some numerical results for the helicity amplitudes that illustrate the role of the f 2 (1270) depending of the photon virtualities.
Our results will be crucial for a model-independent evaluation of the f 2 (1270) contribution to HLbL scattering in the anomalous magnetic moment of the muon, which so far has only been estimated within a Lagrangian-based hadronic model as a narrow resonance [78,79]. To this end, we demonstrated how all helicity amplitudes can be derived numerically from the unitarization of pion-pole and vector-mesonexchange contributions, with parameters determined from the comparison to the measured γγ → ππ cross section. The same intermediate states in the dispersion relation for HLbL scattering should thus allow one to capture effects corresponding to the f 2 (1270) beyond the narrow-width approximation, and without further assumptions on the form factors corresponding to helicity amplitudes that cannot be probed by available data. Even if data for the offshell process γ * γ * → ππ were available, currently under study at BESIII [80] and potentially in the future at Belle II [81], the weighting with respect to energies and virtualities in the g − 2 integral need not resemble the one in the cross section, which makes a detailed understanding of the various helicity amplitudes all the more important. In this way, the f 2 (1270) will be an important test case also for other resonances in the 1-2 GeV region that are hard to describe explicitly in terms of their decay channels, but still need to be reliably estimated to confront the Standard-Model prediction for the muon g − 2 at the level of accuracy anticipated for the E989 Fermilab experiment [82]. In particular, the isospin matrices in (A.2) ensure that the standard form of Watson's theorem [66] holds, i.e. in the elastic regime the phases of the γ * γ * → ππ helicity partial waves agree with the corresponding ππ phase shifts. In the same way, the exchange of a vector meson based on a Lagrangian model [35] leads to where F V π (q 2 ) denotes the V π transition form factor and C V is the coupling in the Lagrangian model, related to the decay width V → πγ by [35] (A.7) The helicity partial waves are h V J,1 (s) = −C 2 V F V π (q 2 1 )F V π (q 2 2 ) 4 σ π (s)λ 1/2 12 (s) s(M 2 V − M 2 π ) 2 + s q 2 1 q 2 2 + M 2 V (s − q 2 1 − q 2 2 ) + M 2 π (q 2 1 − q 2 2 ) 2 − s(q 2 1 + q 2 2 ) λ 12 (s) , x V = s − Σ ππ + 2M 2 V σ π (s)λ

B Anomalous singularities in the modified MO representation
In this appendix, we explain in more detail how the dispersive integrals over the anomalous LHC can be computed in a numerically stable way. The appearance of the Omnès function in the unitarized case leads to additional complications compared to the description in Sect. 3.3. Instead of (3.37), the anomalous integral is given by .
One could proceed as in (3.37) and directly subtract the part of the expanded integrand that diverges for s → s + . However, in this case, the numerical instabilities in the function∆(s, s ) become worse, since the cancellation of the two divergent expressions involves the Omnès function and its derivatives, both of which are calculated only numerically. However, these intricate cancellations can be avoided if not the full integrand including the Omnès function is expanded but only the part involving the resonance LHC, exactly as in (3.36). We define (B.4) The first integral is manifestly finite and the cancellation of divergences in∆(s, s ) is identical to the case without unitarization, in particular no further instabilities are introduced by derivatives of the Omnès function. The second integral can be split into a path up to close to the singularity at s + and an integral circling around the singularity: Im g(s, s ) + 2π −π π dφe iφ 1 Ω(s (φ)) g(s, s (φ)).

(B.5)
For → 0, both integrals contain divergent pieces that cancel in the sum. The result can be obtained by multiple integration by parts, leading to The integral in the first term has to be done numerically but does not introduce any instabilities as the integrand vanishes as a square root for s → s + . The functions b k (s ) depend on the first k + 1 derivatives of the Omnès function at s . The second term denotes the lower boundary term of the integration by parts at s − cut . The divergent upper boundary term at s + − has canceled against the circular integral in (B.5). The above expression is numerically stable as long as s + is not too close to s − cut , i.e. as long as we stay away from the exceptional point q 2 1 q 2 2 = (M 2 V − M 2 π ) 2 . In the vicinity of this singular point, one can further expand the coefficients a k (s) around q 2 1 ∼ (M 2 V −M 2 π ) 2 q 2 2 to obtain a manifestly finite expression.