Dispersion relations for γ∗γ∗ → ππ: helicity amplitudes, subtractions, and anomalous thresholds

We present a comprehensive analysis of the dispersion relations for the doubly-virtual process γ∗γ∗ → ππ. 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ès 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ès 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 γ∗γ∗ → ππ, which allows us to predict the doubly-virtual responses of the f2(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.

JHEP07(2019)073
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 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 section 2, based on which we then derive all relevant kernel functions that define the full system of Roy-Steiner (RS) equations for the partialwave helicity amplitudes. In section 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 section 4, before we conclude in section 5 and comment on the implications of our results for the application to HLbL scattering. 2 Helicity amplitudes and Roy-Steiner equations

JHEP07(2019)073
and the contraction with polarization vectors finally defines the helicity amplitudes according to Next, the kinematic invariants are taken as 1 satisfying 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 (2.13) 1 We denote γ * γ * → ππ as the s-channel process, which is the canonical choice in the context of HLbL.

JHEP07(2019)073
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

JHEP07(2019)073
where q 3 = p 2 − p 1 and (2. 19) 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 appendix 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 scattering angles, leading to with Born terms A π i given in appendix A. In writing (2.24) we have implicitly assumed that the Born-subtracted 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 appendix 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 (2.26) 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 JHEP07(2019)073 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 appendix A, this leads to the representation From these relations we may read off the kinematic singularities as follows: at threshold the combinationQ 12 (s)) ∼ s for s → 0, the 1/s singularities in σ 2 π (s) cancel and we find 2 2 The signs and factors in hJ,3(s) and hJ,4(s) have been chosen to simplify the form of the final kernel functions.
where the functionsh J,i (s) are regular at the thresholds. s −1/2h J,i (s), i = 1, 2, 5, are finite at s = 0 for J ≥ 2. Further zeros are possible for specific contributions, but should be considered to be of dynamical origin [64]. The functionsh J,i (s) have LHCs, encoded in the Q J (x t ), as well as the right-hand cuts from the direct-channel contribution in A i (s, t ). As a cross check on (2.30), we recover the fixed-s resonance partial waves lim a→∞ h V,HDR J,i (s) if the δ-function imaginary parts of the narrow-width resonance amplitudes are inserted. 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ỹ

JHEP07(2019)073
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 theh J,i , these kernel functions are lengthy and therefore not reproduced here. We will give them in a more convenient basis in section 2.5 and their role in the MO solution will be studied in detail in section 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 section 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ȟ 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:

JHEP07(2019)073
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 section 3.4 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 section 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ȟ 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 are indeed fulfilled, but the two sum rules are violated since the asymptotic behavior of the resonance contribution is worse than what we assume for the rescattering contribution: However, the asymptotic behavior in (2.54) still implies that all S-and D-waves fulfill unsubtracted DRsȟ with ππ phase shifts δ J (s), the solution to the MO problem can be given immediately in terms of the Omnès functions .

JHEP07(2019)073
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 section 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 section 3.4 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.

MO solution: D-waves
With the diagonalization of the D-wave system derived in section 2.5, the MO solutions follow immediately, as the defining DRs are given in decoupled form (2.42) or (2.44). The solution reads 3ȟ where δ 2 is the D-wave ππ-scattering phase shift and Ω 2 the corresponding Omnès function. From the functionsȟ 2,i , we obtain the original helicity partial waves by inverting the basis JHEP07(2019)073 change (2.39): where we have introduced the combinationš where∆ ± 23 and∆ ± 45 are defined in analogy to (3.10). We also remark that for spacelike virtualities, the zeros s ± of λ 12 (s) need to be analytically continued according to JHEP07(2019)073 The complete set of MO D-wave solutions is given by (3.8), together with the basis change (2.39) and its inverse (3.9). In particular, the solution (3.11) amounts to a rewriting of the S-wave solution (3.4), which can indeed be cast into the form (3.11) once expressed in terms of s Ω 0 (s), (3.14) witȟ π ds s 2s + Imȟ 2,1 (s ) + √ 6s + Imȟ 2,3 (s ) + Imȟ 2,4 (s ) .  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.

Subtractions and the f 2 (1270) resonance
In [35] it was shown that sum rules for the subtraction constants in the modified Omnès representation are nearly fulfilled, making an approximate description possible that is based on an unsubtracted system. Surprisingly, the same observation does not hold for the Dwave analog of (3.4), where the vector-meson LHC is treated as part of the inhomogeneity. Here, we want to 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 an asymptotic behavior If the Omnès function behaves as Ω(s) s −1 , the 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: (3.20) In order for an unsubtracted standard MO solution to work, one would have to assume a cancellation of the bad high-energy behavior of the vector-meson LHC with the highenergy behavior of the rescattering contribution in h std (s), which seems unlikely. Therefore, the standard MO form is expected to work only when subtractions are introduced. 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. While the MO solution in the standard form keeps the bad high-energy behavior of the resonance LHC, 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. According to the general considerations of section 2.6, the asymptotic behavior for the D-wave rescattering (2.48) should make an unsubtracted dispersion relation possible for the Born-subtracted part, hence a priori one would expect the modified MO solution to work even without subtractions.

JHEP07(2019)073
Checking the MO solutions numerically, indeed it turns out that the unsubtracted standard MO form does not reproduce the peak of the narrow f 2 (1270) D-wave resonance. The effect on the resonance peak 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) The resonance peak is thus described exclusively by the subtraction term that we dropped in (3.21). Using an unsubtracted standard MO solution corresponds to fixing the subtraction constant with a sum rule that cannot be expected to hold and therefore leads to an incorrect description of the resonance. Such a situation indeed occurs for some of the D-waves, which is why in the following 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 γ * γ * → ππ. Avoiding the introduction of subtraction constants is advantageous for the generalization to the singly-or doubly-virtual case, because otherwise their q 2 -dependence would need to be addressed as well.

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 -20 -

JHEP07(2019)073
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 (B.1). 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 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 or- : 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 figure 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 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 a k (s) (s + − s ) (2k+1)/2 +∆(s, s ), (3.36) 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 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 -23 -JHEP07(2019)073 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 appendix 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 section 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  . 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, all partial waves), including the I = 0 unitarization of S-and D-waves (red dot-dashed), and the full solution (black solid).
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 Dwaves with the S-waves from [42,43] (as well as the higher partial waves for the pion pole without rescattering). 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 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 figure 3 follow this convention. The relevant helicity amplitudes in the on-shell case are h 0,++ (s) = 1 2ȟ 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 4 A similar observation was made in [40], where the authors argued that the difference between the fit values for the photon couplings and the ones extracted from the radiative widths reflected SU(3) uncertainties. We disagree with that statement: if the deficit were due to SU(3) uncertainties, it should disappear once the known couplings for the individual states, ω, ρ ± , ρ 0 , are used instead of a common SU(3) coupling, but this is not the case.

JHEP07(2019)073
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 vectormeson D-waves, see figure 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.23). 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. 5

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 appendix 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][78]. Here, we illustrate the numerical solution by approximating the dependence on the photon virtuality by a vector-meson-dominance (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,79]. 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 [80,81], 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 figure 4 for several singly-virtual cases and in figure 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 Dwave 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 highenergy 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.

JHEP07(2019)073
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 [82,83]. To this end, we demonstrated how all helicity amplitudes can be derived numerically from the unitarization of pion-pole and vector-meson-exchange 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 [84] and potentially in the future at Belle II [85], 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 [86].

JHEP07(2019)073
for charged (c) and neutral (n) pions and isospin I = 0, 2, respectively, which in practice implies an overall sign in the transition from (A.1) to the helicity amplitudes. In these conventions the Born-term partial-wave projections become [43] 2sσ π (s) 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) -30 -

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 section 3.4. Instead of (3.37), the anomalous integral is given by . (B.1) 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 functioñ ∆(s, s ) become worse, since the cancellation of the two divergent expressions involves (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: 1 Ω(s ) Im g(s, s ) + 2π −π π dφe iφ 1 Ω(s (φ)) g(s, s (φ)). 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.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.