Towards an improved understanding of $\eta \to \gamma^* \gamma^*$

We argue that high-quality data on the reaction $e^+e^-\to \pi^+\pi^-\eta$ will allow one to determine the doubly-virtual form factor $\eta\to \gamma^*\gamma^*$ in a model-independent way with controlled accuracy. This is an important step towards a reliable evaluation of the hadronic light-by-light scattering contribution to the anomalous magnetic moment of the muon. When analyzing the existing data for $e^+e^-\to\pi^+\pi^-\eta$ for total energies squared $k^2>1\text{GeV}^2$, we demonstrate that the effect of the $a_2$ meson provides a natural breaking mechanism for the commonly employed factorization ansatz in the doubly-virtual form factor $F_{\eta\gamma^*\gamma^*}(q^2,k^2)$. However, better data are needed to draw firm conclusions.


Introduction
Transition form factors contain important information about the properties of the decaying particles. Additional interest in meson decays with one or two virtual photons in the final state comes from the fact that the theoretical precision of the Standard Model calculations for the anomalous magnetic moment of the muon is significantly affected by the one of hadronic light-by-light (HLbL) scattering, where the latter appears as a sub-amplitude; see Fig. 1. A few years ago, the Bern group reported on important progress towards a model-independent determination of the HLbL contribution based on dispersion theory [1][2][3]. In principle this allows for an analysis of similar rigor as commonly applied for the hadronic vacuum polarization [4]. 1 The updated experimental value for (g − 2) µ , combining the BNL measurement [6] with the first Fermilab results [7,8], shows a 4.2σ tension as compared to the current theory consensus on the Standard Model value [4,, and the upcoming improvement a e-mail: holz@hiskp.uni-bonn.de b e-mail: kubis@hiskp.uni-bonn.de 1 See also the alternative approach to analyze the muon's Pauli form factor dispersively [5]. in the experimental statistics [34] demands further efforts to reduce the theoretical uncertainty. Reference [4] identifies the η and η ′ pole terms as major outstanding contributions in the dispersive approach to HLbL, next to axial vectors [35,36] and the matching to short-distance constraints [21,[29][30][31][37][38][39][40]. A refined understanding of the doubly-virtual η transition form factor is the object of this article; see also the review Ref. [41]. In Ref. [42] the isovector contribution of the singlyvirtual form factor F ηγ * γ (q 2 ) ≡ F ηγ * γ * (q 2 , 0) was calculated for virtualities q 2 ≪ 1 GeV 2 from data on η → π + π − γ and the pion vector form factor via a dispersion integral. This was done with the help of a convenient parameterization of the corresponding ππ invariant-mass distribution derived in Ref. [43], see also Ref. [44]. In particular, it was demonstrated that when the high-statistics data of Ref. [45] were used to fix the η → π + π − γ input, this procedure leads to a determination of F ηγ * γ (q 2 ) with an accuracy higher than that of the most recent direct measurements [46][47][48][49]. Especially this is the case since the isoscalar contribution is negligibly small.
In this article, the program to pin down the η transition form factor with high accuracy will be extended to the doubly-virtual form factor F ηγ * γ * (q 2 , k 2 ) in the kinematic regime q 2 ≪ 1 GeV 2 and k 2 > 1 GeV 2 . The analysis is based Fig. 2 The reaction e + e − → ηπ + π − . The wiggly, dashed, solid, and double lines denote photons, pions, η, and a 2 meson, respectively. The gray blob stands for the ππ final-state interactions.
on input from the reaction e + e − → ηπ + π − , see Fig. 2, which plays a very similar role to data on e + e − → 3π for the π 0 transition form factor [26,27,50]. The method provides access to the doubly-virtual form factor in a kinematic regime where it is yet unknown. In particular it allows for a test whether the factorization ansatz F ηγ * γ * q 2 , k 2 =F ηγ * γ q 2 F ηγγ * k 2 (1) for the normalized form factorF ηγ * γ * (q 2 , k 2 ) = F ηγ * γ * (q 2 , k 2 )/F ηγ * γ * (0, 0) is still valid in the kinematic regime specified above. 2 It should be stressed that although a direct measurement of the doubly-virtual form factor via e + e − → ηe + e − is in principle possible, we still expect our method to lead to higher accuracy, simply because the hadronic rates for γ * → ηπ + π − are a factor 1/α 2 em larger than those for γ * → ηe + e − .
In Ref. [51] it was argued that the a 2 (1320), a tensor resonance with I G (J PC ) = 1 − (2 ++ ), should provide the leading left-hand-cut contribution to the decay amplitude for η → π + π − γ, and accordingly distort the spectra significantly, however, only beyond the kinematic region accessible in the direct measurement of the decay. Interestingly all necessary parameters can be fixed from data directly. The claim was corroborated by an analysis of the data for η ′ → π + π − γ [52,53]. In this work we therefore also include the a 2 contribution. Diagrammatically this amounts to including a 2 crossed-channel exchange as shown in Fig. 2(b), in addition to the structureless vertex of Fig. 2(a). Unfortunately, the differential data presently available for the reaction e + e − → ηπ + π − from the BaBar [54,55], SND [56], and CMD-3 [57] collaborations are rather limited: only BaBar and CMD-3 provide ππ spectra, and both with the shortcoming that they are not given at a fixed value of k 2 , but for √ k 2 integrated in ranges from 1.0 − 4.5 GeV [54], 1.4 − 1.8 GeV [55], or 1.3 − 1.8 GeV [57]. As we argue below, analyzing the two BaBar data sets, this integration limits the extraction accuracy of the form factor, since changes in the ππ spectrum cannot only be induced by changes in the explicit dependence of the amplitude on the ππ invariant mass, but also by changes in the parameterization of the line shape of the ρ ′ . The latter is the first excited resonance of the ρ meson and dominates the total cross section.
The paper is structured as follows. In Sect. 2 we generalize the formalism of Ref. [42] to the doubly-virtual η transition form factor, based on amplitudes for η → π + π − γ * . Subsequently, we show how the dispersive representations for η → π + π − γ discussed in Refs. [43,51] can be extended to the virtual-photon case of interest here. In Sect. 3 we discuss the parameterization used for virtualities beyond 1 GeV 2 , which are beyond the range of applicability of the dispersive approach employed in this work. In Sect. 4 we present and discuss the results of our analysis of e + e − → ηπ + π − total and differential cross section data, without and with the inclusion of the a 2 contribution. Section 5 reflects on the impact of our findings for the η transition form factor. We close with a summary and discussion; some technicalities are relegated to an appendix. The η → γ * γ * transition form factor (TFF) F ηγ * γ * (q 2 , k 2 ) is defined by the matrix element where j µ = (2ūγ µ u −dγ µ d −sγ µ s)/3 denotes the electromagnetic current, q and k are the photon momenta, and the convention ε 0123 = +1 is used. The decay into two real photons is driven by the chiral anomaly [58,59], which fixes the normalization F ηγγ ≡ F ηγ * γ * (0, 0) of the TFF, where α em = e 2 /(4π) is the fine-structure constant. Isospin symmetry allows us to decompose F ηγ * γ * (q 2 , k 2 ) into an isovector-isovector (I = 1) and an isoscalarisoscalar (I = 0) part, which, for the normalized doublyvirtual TFFF ηγ * γ * (q 2 , k 2 ) ≡ F ηγ * γ * (q 2 , k 2 )/F ηγγ , can be written according tō As vector-meson-dominance estimates demonstrate, the isoscalar part is strongly suppressed [41,42], such that for the TFF's low-energy properties, its contribution is within the isovector part's uncertainty. A single-variable dispersion relation for the isovector part is at low energies dominated by the two-pion intermediate state, see Fig. 3, which yields where F V π (t) is the pion vector form factor, and F ηππγ * (t, k 2 ) denotes the projection of the η → π + π − γ * amplitude onto the pion-pion P-wave, to be defined formally below. We expect the dispersion relation Eq. (5) to hold up to q 2 ∼ 1 GeV 2 , even for a much wider range in the second argument k 2 . We will briefly comment on the extent to which two pions alone saturate this isovector dispersion relation below in Sect. 5.

Dipion amplitudes
The decay amplitude for η(p η ) → π + (p + )π − (p − )γ ( * ) (k) can be written as with the Mandelstam variables given as s = (p η − p + ) 2 , An expansion of F (s,t, u, k 2 ) in pion-pion partial waves proceeds in odd waves only, and is totally dominated by the P-wave [51], denoted by F ηππγ * (t, k 2 ) in the following, where and the Källén function is defined by In the real-photon case k 2 = 0, the η → π + π − γ differential decay rate with respect to t can be written as [43] dΓ η On the other hand, for k 2 > (M η + 2M π ) 2 , the matrix element (6) can be assessed in the reaction e + e − → ηπ + π − . We write the differential cross section as where the normalization is kept unspecified due to the BaBar differential cross sections being provided in arbitrary units [54,55]. The total cross section data sets used, however, are provided in meaningful units. Hence, for the total cross section we use A central element of the amplitude F ηππγ * (t, k 2 ) is the universality of the elastic pion-pion P-wave final-state interactions as encoded in the Omnès function [60] where δ 1 1 (t) denotes the pion-pion P-wave phase shift. In particular, Ω (t) contains the physics of the ρ(770) resonance. One of the most obvious applications to make use of the Omnès function is the pion vector form factor, which can be extracted directly from data on e + e − → π + π − or τ − → π − π 0 ν τ . It can be written as where R(t) is a function free of cuts up to inelastic effects, which can in practice be approximated by a linear function up to t ∼ 1 GeV 2 [42,52]. Utilizing the phase shift parameterization specified in Appendix A.1 and fitting to τ decay data [61], this linear function is found to be R(t) = 1 + 0.126(2) GeV −2 t. Comparable parameterizations [62,63] were used in previous studies of the η TFF [42,43]. Similarly, is a function free of right-hand cuts in t up to inelastic thresholds. In Refs. [43,45] it was shown that for the real-photon case, the function P(t, 0) is sufficiently well approximated as a linear polynomial where A η ππγ is a normalization constant, to describe the highaccuracy η → π + π − γ decay data obtained by the KLOE collaboration [45]. The slope parameter α Ω is found to be [51] α Ω = 1.52 (6) with only the statistical uncertainty included. 3 As P(t, 0) is not expected to grow for large t on fundamental grounds, the representation Eq. (16) is unlikely to hold beyond the small range in t accessible directly in the decay η → π + π − γ. Indeed, high-statistics data for the analogous decay η ′ → π + π − γ require a polynomial of second order [52,53] instead, motivating an alternative parameterization that we will also test below. A possible physical motivation for such a curvature term was given in Ref. [51], where it was demonstrated that the left-hand cut induced by the a 2 (1320) distorts the linear behavior of P(t, 0) significantly beyond t ≤ M 2 η . In order to include the left-hand-cut contribution, we need to generalize the dispersive representation of F ηππγ * (t, 0) according to [51] with where F a 2 (s,t, u, 0) comprises the a 2 s-and u-channel exchange amplitudes (for real photons), see Ref. [51] for details. The subtraction constant α Ω [a 2 ] that takes over the role of the slope parameter α Ω is shifted only marginally: extracted from a fit to η → π + π − γ (where t < 0.25 GeV 2 ), it reads [51] α Ω [a 2 ] = 1.42(6) GeV −2 .
3 The analyses [43][44][45] employ Eq. (15) with Ω (t) replaced by F V π (t), due to the observation that the linear slope in R(t) relating both quantities is small. The analogously defined parameter α was then determined to be α = 1.32(13) GeV −2 , combining statistical and systematic errors, which is slightly shifted compared to Eq. (17).
We define the function as a straightforward generalization of Eq. (15) (for k 2 = 0). This function is still free of right-hand cuts, but now contains the left-hand cut due to a 2 exchange. P (a 2 ) (t, 0) can be approximated by in the range 4M 2 π ≤ t ≤ 1 GeV 2 . It turns out that with α a 2 = 0.28 GeV −2 and β a 2 = −0.66 GeV −4 , this approximation works to better than 1% accuracy. The constant A η ππγ [a 2 ] is chosen such that F ηππγ * (t, 0) is normalized to the experimental rate for η → π + π − γ. Note that if the a 2 left-handcut contribution is turned off, the constants A η ππγ and A, see Eq. (20), coincide. The overall strength of the a 2 contribution is phenomenologically known from a 2 branching fractions to better than 10% accuracy [51]. Given the size of the errors in the data analyzed in the present study, we will neglect this source of uncertainty in what follows.
The same method that allowed us to connect η → π + π − γ to the isovector component of η → γ * γ permits to connect γ * → ηπ + π − to the isovector component γ * → ηγ * and thus to the doubly-virtual η transition form factor. To this end, we need to generalize the dispersion relation Eq. (5), and hence the description of F ηππγ * (t, k 2 ), to an off-shell photon with invariant mass squared k 2 . We begin by noting that the left-hand-cut contribution of Eq. (19) changes naturally for k 2 = 0 due to the changed kinematics. We define The angular integral over the a 2 exchange amplitudes can be done analytically and yields the hat function (cf. Appendix A.2 for details on the analytic continuation) where the coupling constant g T (c T ) has been determined from the branching ratio of a 2 → πη (a 2 → πγ) [51] and F π = 92.28(10) MeV [64] is the pion decay constant. The dispersion integral contained in can be carried out numerically via path deformation into the complex plane in order to circumvent the singularities of the integrand, inspired by the methods of Ref. [65]. Due to the need to deform the integration path into the complex plane, we utilize a ππ phase shift δ 1 1 calculated via the inverseamplitude method (we employ an approximation to the twoloop representation of Ref. [66]), which provides us with an analytic expression; see Appendix A.1. Equations (25) and (26) encode a nontrivial entanglement of the dependences on t and k 2 that clearly does not factorize in a simple manner.
We now make the following ansatz for the complete γ * → ηπ + π − form factor: i.e., any additional k 2 dependence is assumed to be multiplicative (or factorizing). Note thatF ηγγ * (k 2 ) cannot be treated within the formalism of Ref. [42] for an application to e + e − → ηπ + π − data, since with k 2 > 1 GeV 2 it is to be evaluated in a kinematic regime where the original method, relying on elastic unitarity, is no longer applicable. However, as will be shown in Sect. 3, for the analysis aimed at here we only need a convenient parameterization ofF ηγγ * (k 2 ). Alternatively, we will revert to the simpler approach to replace F (a 2 ) ηππγ * (t, k 2 ) by F ηππγ * (t, 0), given by either the linear or the quadratic approximations P (1,2) (t, 0). In these limits, F (a 2 ) ηππγ * (t, k 2 ) and, via the dispersion integral (5), the (isovector part of the) doubly-virtual TFF, become factorizing in their dependences on the two respective variables, and we can check the consistency of the assumption by testing to which extent the linear subtraction constants of Eqs. (20) and (26) can (on average) be approximated by a constant consistent with the value for the real-photon case. In this limit (only),F ηγγ * (k 2 ) also corresponds to the normalized singly-virtual η TFFF ηγγ * (k 2 ).
In order to fix α * Ω [a 2 ](k 2 ) andF ηγγ * (k 2 ), which provide the input for the calculation of the doubly-virtual transition form factor F ηγ * γ * (q 2 , k 2 ), data on e + e − → ηπ + π − need to be analyzed, with particular attention to the distributions with the dipion invariant masses. However, data for the ππ spectrum are available only with a simultaneous integration over the initial energy k 2 . Therefore it is not possible to extract values for α * Ω [a 2 ](k 2 ) as a function of k 2 , but only to extract some averaged valueᾱ * Ω [a 2 ]. Moreover, since the value ofᾱ * Ω [a 2 ] also influences the shape of the total cross section for e + e − → ηπ + π − , a combined fit to the ππ spectrum and the total cross section is mandatory, which calls for a parameterization of the data in kinematic regimes that cannot be captured by the dispersion integrals employed in this work. We describe this parameterization in the following section.
3 Form factor parameterization above 1 GeV

Inclusion of excited isovector vector resonances
For the fit to the total cross section we need to parameterize the functionF ηγγ * (k 2 ) for k 2 > 1 GeV 2 . As one can see from Fig. 4, the total cross section for e + e − → ηπ + π − shows a very prominent contribution from an excited ρ resonance, the ρ ′ or ρ(1450). It turns out that the low-energy side of the spectrum calls for an inclusion of the ρ(770) resonance, which enters through its upper tail. Moreover, the addition of a higher resonance, the ρ ′′ or ρ(1700), improves the data fit significantly. Note that most experimental or phenomenological analyses of the reaction e + e − → ηπ + π − employ two or three vector resonances [67][68][69][70]. It should be stressed that what is crucial for the analysis is not that we have the correct physics optimally built into the parameterization of the functionF ηγγ * (k 2 ), but all we need is a convenient repre-sentation that describes the data for the quantity of interest for this analysis: the k 2 -integrated ππ spectrum. Sincẽ F ηγγ * (k 2 ) enters here as a weight factor-cf. Eq. (11)-it is sufficient to parameterizeF ηγγ * (k 2 ) by a sum of Breit-Wigner functions. Clearly in doing so we should not expect the resonance parameters to agree with those of the Particle Data Group (PDG), since Breit-Wigner parameters are reaction dependent and a sum of Breit-Wigner functions violates unitarity. More explicitly, the Breit-Wigner (BW) functions we employ are given by with j ∈ {ρ, ρ ′ , ρ ′′ } and the energy-dependent widths written as [71] Γ where the center-of-mass momentum is given by Additionally, we employ another variant, in which the ρ ′ and ρ ′′ widths are parameterized by the decay into the ωπ final state: We parameterize the functionF ηγγ * (k 2 ) ∝ ∑ j c j M 2 j BW j (k 2 ), where j ∈ {ρ, ρ ′ , ρ ′′ }. Through the normalization constraintF ηγγ * (0) = 1 we eliminate the coupling c ρ : The remaining couplings are used as fit parameters. While the (large) Breit-Wigner coupling to the ρ ′ is taken to be real as required by unitarity, we see a substantial improvement in the fit quality if the (small) coupling to the ρ ′′ is allowed to take complex values, which likely subsumes the effects of yet higher resonances and the more complicated coupled-channel structure in that energy region. Additionally, the masses M ρ ′ and M ρ ′′ as well as the widths Γ ρ ′ (M 2 ρ ′ ) and Γ ρ ′′ (M 2 ρ ′′ ) of the ρ ′ and ρ ′′ are fitted to the data. For the ρ we use the PDG values for mass and width.
Since we do not expect the form factor ansatz in Eq. (24) to hold up to arbitrarily high energies, we employ the following prescription for a certain cutoff value s c : and the dispersion integral .
This continuation ensures that the form factor drops off like 1/t and 1/k 2 for high energies and is continuous everywhere when going from one region to another. Through treating the cutoffs in t and in k 2 in the same way, it reflects a symmetry that is important for a later application in the η transition form factor in the HLbL contribution to the muon g − 2. We test the independence of the choice of continuation point s c by varying it in the fits according to s c ∈ {3, 4, 5} GeV 2 .

Finite-width effects in the a 2 exchange
In the tensor model used to describe the a 2 -exchange diagrams, the a 2 meson is assumed to have no width. However, in reality the a 2 (1320) has a width of about 100 MeV, and in the reaction e + e − → ηπ + π − , where the a 2 can be produced quasi on-shell, this effect is surely not negligible. To cure this caveat, we define a dispersively improved Breit-Wigner function BW (s), using the imaginary part of BW (s) as in Eq. (28) as the spectral function above a certain threshold s thr [72,73], where in the last step we have written the Cauchy kernel as a Breit-Wigner function with mass parameter √ x and width zero. We therefore generalize both dispersion integral D(t, k 2 ) and hat functionF a 2 (t, k 2 ) with the a 2 mass dispersed according to this spectral function, where C ∈ {D,F a 2 }. We employ the parameterization [74,75] Γ (s) = Γ a 2 ∑ i∈{η,ρ} where M ρ (M η ) is the mass of the ρ (η) meson, Γ a 2 the total a 2 width, and the barrier factor is given by R = 5.2 GeV −1 [75]. This energy-dependent width explicitly takes into account the a 2 decays into πη and πρ with branching ratios p η = 0.17 and p ρ = 0.83. We evaluate C(t, k 2 ; M 2 a 2 ) for 20 different values of M 2 a 2 in the range s thr , 2.3 GeV 2 , where s thr = (M η + M π ) 2 . For D we use a grid of 107×150 values for t, k 2 ∈ [4M 2 π , 3 GeV 2 ]. The integral of Eq. (37) was then carried out using the trapezoidal rule for each point of the (t, k 2 ) grid. In the fitting routine we make the replacementŝ F a 2 (t, k 2 ) →F a 2 (t, k 2 ) and D(t, k 2 ) → D(t, k 2 ) .
4 Data analysis of e + e − → ηπ + π − In the following, we analyze the e + e − → ηπ + π − data by the BaBar collaboration [54,55] in three successively refined models for the ππ invariant mass distribution: assuming a linear or a quadratic polynomial multiplying the final-state interactions as parameterized by the Omnès function, and using the full-fledged description of left-hand cuts due to the a 2 , which automatically generates breaking of the factorization assumption. We fit to data from BaBar published in 2007 [54] and 2018 [55]. The differential cross section data from 2007 are integrated over √ k 2 in [1.0, 4.5] GeV, the ones from 2018 in the range √ k 2 ∈ [1.4, 1.8] GeV. The ππ spectrum is dominated by the ρ, we choose to fit the differential cross section data up to √ t = 1 GeV. The tail of the Breit-Wigner parameterization ofF ηγγ * (k 2 ) is rather uncontrolled for larger k 2 , thus we fit the total cross section data up to √ k 2 = 1.9 GeV only.

Linear model
As in the first dispersive attempts that dealt with the η TFF [42,43] we employ a linear function as part of F ηππγ * (t, k 2 ), which does not include any effects of the a 2exchange left-hand-cut contribution: For values of the ππ invariant mass squared larger than a certain cutoff parameter s p we let the form factor fall off like 1/t by setting the linear polynomial to a constant. The constant term in the polynomial part of Eq. (40) is fixed by fitting the corresponding linear representation to η → π + π − γ data measured by KLOE [45] to A = 18.0(4)GeV −3 . The results of these fits with cutoff parameter √ s p ∈ {0.9, 1.0, 1.2} GeV are shown in Table 1. The central finding is that the extracted values for α * Ω are consistently and significantly smaller than those found for α Ω , see Eq. (17).

Quadratic model
As noted in Ref. [51] for the real-photon case, the inclusion of the a 2 -induced left-hand-cut contribution adds a curvature part to the form factor F ηππγ * (t, k 2 ). In the quadratic model, we do not use the explicit form of the a 2 exchange amplitudes, but restrict ourselves to the approximation in Eq. (23). The fit of this representation to η → π + π − γ decay data yields A = 17.4(4)GeV −3 . Note that this form omits all the k 2 -dependence in F ηππγ * (t, k 2 ) other than through the multiplicative factorF ηγγ * (k 2 ), cf. Eq. (32): where the linear coefficientᾱ * Ω [a 2 ] is used as a fit parameter and the parameters α a 2 and β a 2 are held fixed at their respective values given below Eq. (23). For values of the ππ invariant mass squared larger than a certain cutoff parameter s p we, again, let the form factor fall off like 1/t by the prescription The fit results for cutoff parameter √ s p ∈ {0.9, 1.0, 1.2} GeV are displayed in Table 2. Compared to the results of the linear model, the linear parameterᾱ * Ω [a 2 ] comes out substantially closer to the value of the linear subtraction constant α Ω [a 2 ] for all tested values of s p , cf. Eq. (21), which has been obtained from η → π + π − γ decay data, and compatible within uncertainties throughout. The sizable shift of the slope between the linear and quadratic model can be understood in the following manner: in the Table 1 Results for the simultaneous fit to differential and total cross section data from BaBar (2007) [54] and BaBar (2018) [55], employing the linear model of Eq. (40). The differential cross section data have been fitted up to 1 GeV, the total cross section data up to 1.9 GeV. (χ 2 /dof) y c labels the reduced χ 2 for the respective data sets, c =dif/tot for the differential/total cross section of BaBar in year y. The parameter N y labels the unphysical normalization of the differential cross section data of BaBar in year y. Columns 2-4 refer to fits employing the BW parameterization of ρ ′ and ρ ′′ widths according to ππ decays, columns 5-7 display those relying on the ωπ running width model.   (20) real-photon case, the main impact of the a 2 contribution on the ππ spectrum is to provide a downward curvature, with the quadratic coefficient β a 2 being sizable and negative. Accordingly, a fit with a linear term only must find a reduced strength in order to reproduce the data in the region of the ρ resonance.

Full model
Here, the full representation of F ηππγ * (t, k 2 ) in Eq. (27) is used for the e + e − → ηπ + π − fits, i.e., including the k 2 -dependent a 2 -induced left-hand-cut contribution. The parameter A = 17.2(4)GeV −3 in Eq. (26) is fixed by fitting the representation in Eq. (19) to η → π + π − γ decay data measured by KLOE [45]. The fit results for cutoff parameters s c ∈ {3, 4, 5}GeV 2 according to the high-energy prescription in Eq. (33) are summarized in Table 3. We observe that the ωπ parameterizations of the ρ ′ and ρ ′′ widths result in a slightly better fit quality. The mass and width parameters of the ρ ′ and ρ ′′ , which must not be confused with pole parameters, vary strongly between both descriptions of the energydependent width, most prominently for the ρ ′ . A compar-  = 400(60)MeV) is therefore not particularly meaningful. However, the linear subtraction constant, the central result for our investigation of potential factorization breaking, only varies within the error margins between the two fits. Additionally, the resulting spectra for cutoff parameter s c = 3 GeV 2 are shown in Figs. 4 and 5, where the error bands are generated by the fit uncertainties taking their correlations into account. The description of differential spectra turns out to be slightly worse than for the linear/quadratic models. Those, however, do not capture the nontrivial k 2dependence of the amplitude illustrated in Fig. 6: the change of the left-hand-cut contribution with k 2 is seen to be quite significant. In Fig. 5 the fit result of the full model can be compared to the central value of the quadratic one. Differences between the two curves can only be seen at the lower flank of the ρ resonance. Moreover, if the uncertainty band were added, it would overlap with the curve of the full model almost everywhere. The large values of the reduced χ 2 of the differential cross section in the full model fit are primarily caused by the flanks of the ρ resonance, as one can see in Fig. 5. Especially below the ρ peak, the χ 2 receives a large contribution. The tension specifically for these ππ invariant masses is surprising, as the dispersive representation should be perfectly reliable there. Doubly-differential data (in the variables t and k 2 ) would be highly desirable to better understand the origin of this discrepancy. The Breit-Wigner coupling of the ρ ′ resonance comes out consistently in all fits; in particular, in contrast to mass and width, it only changes mildly between the two different parameterizations of the energy-dependent width. Moreover, to good accuracy it fulfills the expectation (where in our convention Eq. (32), c ρ = 1), derived from a superconvergence relation for V → P transition form factors that are expected to drop ∝ k −4 asymptotically [76][77][78][79]. The coupling to the ρ ′′ resonance comes out complex, however it is two orders of magnitude suppressed compared to the ρ and more than one order of magnitude with respect to the ρ ′ coupling. Also the mass and width parameters of the ρ ′′ deviate somewhat from the PDG values, M PDG ρ ′′ = 1720 (20)MeV and Γ PDG ρ ′′ = 250(100)MeV, which may be explained by the occurrence of higher resonances in the spectrum that are not included in our parameterization.
The linear subtraction constantᾱ * Ω [a 2 ] comes out consistently for different cutoff parameters s c , but deviates from  15), extrapolated from e + e − → ηπ + π − using both the full and the quadratic model, and the one directly obtained from η → π + π − γ. The band for the full model includes the uncertainty due toᾱ * Ω [a 2 ] as well as the different values for s c ∈ {3, 4, 5} GeV 2 ; see Table 3. In the same way the band for the quadratic model comprises the uncertainty due tō α * Ω [a 2 ] as well as √ s p ∈ {0.9, 1.0, 1.2} GeV; see Table 2. The onset of the high-energy prescription can be seen in the quadratic band around t = 0.81 GeV 2 .
the linear constant found in the polynomial part fitted to η → π + π − γ decay data of 1.42(6)GeV −2 . This is illustrated in Fig. 7, where the e + e − → ηπ + π − form factor in the extrapolation to k 2 = 0 is compared to the representation in Eq. (19) fitted to η → π + π − γ decay data. This is in contrast to the quadratic model, which allows for a smooth connection between the extrapolated e + e − → ηπ + π − and the η decay data. We therefore have to conclude that the dominant left-hand-cut contribution of the a 2 exchange provides a factorization breaking mechanism that is too strong compared to the available data, and likely calls for further effects that mitigate this to some extent.

Consequences for the η transition form factor
The observations made on the different models for F ηππγ * (t, k 2 ) in the previous section immediately translate into the η transition form factor. Figure 8 shows the results for the normalized η → γγ * transition form factor based on a once-subtracted variant of Eq. (5) in the time-like regime (thus enforcing the correct normalization by hand), compared to experimental data [46][47][48][49]. Here, the dispersion integral has been truncated [42], with the upper limit of integration varied between M 2 η ′ and 2 GeV 2 . The error bands in  Fig. 8 Results for the normalized singly-virtual η transition form factor. The very small isoscalar contribution has been omitted. We show η → e + e − γ data from the A2 (squares [47], triangles [49]) and η → µ + µ − γ data from the NA60 (circles [46], diamonds [48]) collaborations. The insert magnifies the low-q 2 region. The gray band shows the result from the full model for s c = 3 GeV 2 . The red (yellow) band represents the outcome of the quadratic (linear) model for s p = 1 GeV 2 . Additionally, for comparison, the result derived from η → π + π − γ is shown as the blue band. The central results are represented by dashed lines colored accordingly throughout. Fig. 8 are generated from the error of the linear subtraction constantᾱ * Ω [a 2 ], the variation of the integration upper limit, the pion vector form factor, and the errors of the branching ratios entering the result. The central results shown as the dashed lines stem from evaluating the integral with cutoff 1 GeV 2 . Not surprisingly, the transition form factor based on the extrapolated quadratic model agrees well with the direct calculation based on η → π + π − γ, while the extrapolated linear and full models yield a form factor significantly below the others. This is also illustrated by the χ 2 /dof, where χ 2 as usual refers to the weighted sum of residuals between data and the central values of the predictions and dof stands for the number of data points (as there are no fit parameters to be adjusted). The full model yields χ 2 full /dof = 2.09, when compared to the linear model with χ 2 lin. /dof = 1.54 and the quadratic model with χ 2 quad. /dof = 0.48. This is to be compared to χ 2 ππγ /dof = 0.78 for the TFF calculated directly based on η → π + π − γ decay data. In particular the data from Refs. [46,48] disfavor the TFF based on the extrapolated full model rather strongly.
These differences are also reflected in the slope parameter of the TFF defined bȳ  ηγ * γ * (q 2 , k 2 ) in the singly-virtual limits k 2 = 0 and q 2 = 0. Results are shown from the full model for s c = 3 GeV 2 . The error band for the k 2 = 0 case is generated as described before, whereas the error band of the q 2 = 0 curve additionally takes the fit errors of the parameters iñ F ηγγ * (k 2 ) and their correlations into account.
According to Eq. (5), its isovector part fulfills a sum rule The slope parameters of the respective models are listed in Table 4. The isoscalar contribution b [42] is inside the quoted uncertainties throughout.
Bose symmetry requires the TFF as well as its isovector part to be symmetric under the exchange of the two arguments q 2 and k 2 . This symmetry is by no means manifest in our model, which in contrast is built onto information in two disjoint kinematic ranges: the dispersion relation (5) is designed for the region q 2 ≤ 1 GeV 2 , while the dependence on the second variable is fitted in the range 1 GeV 2 ≤ k 2 ≤ (1.9 GeV) 2 . Indeed, Fig. 9 demonstrates that the two differ-ent singly-virtual limitsF ηγ * γ * (q 2 = 0, k 2 ) do not coincide in the region of the ρ(770) peak, although the disagreement hardly exceeds the combined width of the two uncertainty bands. Obviously, in order to calculate and predict the limit q 2 = 0, we are required to use the unsubtracted dispersion relation (5), without the integral being cut off. Instead, the high-energy prescription in Eq. (33) has been used above a certain cutoff varied between 1 and 2 GeV 2 . Note that the TFF shown in Fig. 9 is normalized using F ηγγ based on the physical value of the two-photon partial width.F (I=1) ηγ * γ * (0, 0) is hence a test of the hypothesis that two-pion intermediate states saturate the isovector dispersion relation to a large extent: indeed, the normalization is reproduced to 90(8)%. Figure 9 suggests that the remaining deficit, which is due to omitted heavier intermediate states (such as 4π, πω, KK, etc.), likely increases somewhat around the ρ(770), but not dramatically so.

Summary and discussion
In this article, we have investigated the consistency of pionpion spectra in the closely related reactions η → π + π − γ and e + e − → ηπ + π − , which serve as central input to dispersiontheoretical analyses of the singly-and doubly-virtual η transition form factor, respectively. This comparison has been quantified in terms of the residual dependence on the ππ invariant mass squared t that multiplies universal elastic rescattering, as encapsulated in the Omnès function. While the small phase space available in the decay η → π + π − γ is sufficiently accurately described by a residual linear tdependence [43,45], both a study of the dominant left-hand cuts, induced by crossed-channel exchanges of the a 2 tensor meson [51], and data on the closely related η ′ → π + π − γ decay [52,53] demonstrate that at least a quadratic polynomial is required to describe the ππ spectra up to 1 GeV 2 .
The present analysis of the pion-pion distributions in e + e − → ηπ + π − confirms that, indeed, a linear parameterization is insufficient to describe these with a universal slope parameter independently of the dilepton squared invariant mass k 2 . On the other hand, a quadratic polynomial, with the quadratic term in t fixed from the size of curvature induced by the left-hand cut at k 2 = 0, seems well compatible with a naïve factorization assumption for the t and k 2 dependences, and as a result, for the doubly-virtual η transition form factor varying with its two arguments according to Eq. (1). The somewhat puzzling observation is that the natural breaking of factorization, induced by the nontrivial k 2 dependence of the a 2 exchange contribution, is too strong to be consistent with the data: the amplitude we have built on this physically well-motivated picture does not link the available e + e − → ηπ + π − data to those for the η decay in a consistent manner; factorization breaking has to be miti-gated significantly by other mechanisms not investigated so far.
One way how such a damping might be achieved is through different k 2 -dependences of the two considered isobar components in the γ * → ηπ + π − amplitude: we have implicitly assumed that the dominant "ηρ" and the supplementary a 2 π contribution responsible for the left-hand cut evolve with k 2 in parallel. The only way to understand deviations from this simple picture is to analyze differential distributions in ηπ + π − in different bins of k 2 , as opposed to the k 2 -integrated ππ spectra available to date [54,55,57]. Such information, combined with a doubly-dispersive construction based on a η → 2(π + π − ) amplitude at low energies (cf. the corresponding η ′ → 2(π + π − ) decays [80,81] and the discussion in Ref. [41]), would then pave the way towards a more complete understanding of the doubly-virtual η transition form factor. where and σ = σ π (s). In those expressions F represents the pion decay constant in the chiral limit, which can be extracted from the ratio F π /F = 1.062(7) [85][86][87][88][89][90]. The combination of low-energy constantsl 2 −l 1 is later used as a free parameter. Utilizing the IAM, we write the unitarized amplitude as [91][92][93] t IAM (s) = (t 2 (s)) 2 t 2 (s) − t 4 (s) .
Now the ππ P-wave phase shift can be extracted from this expression as δ 1 1 = arg t IAM (s). However, in order to enforce the phase to converge to π for s → ∞, we add terms, inspired by the next-to-next-to leading order amplitude, to t 4 by hand: t 4 (s) → t 4 (s) + t 2 (s) 48π 2 F 2 l s s 2 +l π M 4 π , (A.5) where two additional parametersl s andl π have been introduced. In order to fix the three parameters, we fit the phase parameterization from the IAM amplitude to the solution of the Roy-equation analysis of ππ scattering [82]. The fit up to √ s = 1.3 GeV is displayed in Fig. 10. The parameters resulting from this fit arel 2 −l 1 = 4.47(3),l s = 1.74(3)GeV −2 , andl π = 560(20)GeV −2 .

Appendix A.2: Hat function
The analytical hat function obtained through the angular integral over the a 2 exchange amplitudes in Eq. (25) is implemented in a smooth way by choosing the correct sheets of the logarithm appearing in Q(y), explicitly given by for Im y+1 y−1 < 0 and