The $\gamma\pi\to\pi\pi$ anomaly from lattice QCD and dispersion relations

We propose a formalism to extract the $\gamma\pi\to\pi\pi$ chiral anomaly $F_{3\pi}$ from calculations in lattice QCD performed at larger-than-physical pion masses. To this end, we start from a dispersive representation of the $\gamma^{(*)}\pi\to\pi\pi$ amplitude, whose main quark-mass dependence arises from the $\pi\pi$ scattering phase shift and can be derived from chiral perturbation theory via the inverse-amplitude method. With parameters constrained by lattice calculations of the $P$-wave phase shift, we use this combination of dispersion relations and effective field theory to extrapolate two recent $\gamma^{(*)}\pi\to\pi\pi$ calculations in lattice QCD to the physical point. Our formalism allows us to extract the radiative coupling of the $\rho(770)$ meson and, for the first time, the chiral anomaly $F_{3\pi}=38(16)(11)\,\text{GeV}^{-3}$. The result is consistent with the chiral prediction albeit within large uncertainties, which will improve in accordance with progress in future lattice-QCD computations.


Introduction
The calculation of hadronic scattering processes constitutes a notoriously challenging task in lattice quantum chromodynamics (QCD), given the complications that ensue once multihadron dynamics are properly taken into account [1]. In the past years, significant improvement has been achieved mainly for ππ → ππ [2][3][4][5][6][7][8][9][10][11][12][13], the simplest hadronic scattering process, with computations even at physical pion mass available [14,15]. Meanwhile, the computation and analysis of more complicated processes remains challenging. Already for the generalization when a pion is replaced by an external electromagnetic current, γ ( * ) π → ππ, up to now only two results are available that correctly describe the resonant nature of the process [16][17][18], both obtained at unphysically high pion masses exceeding 300 MeV. Hence, tools are needed to extrapolate the lattice-QCD results to the physical point and confront them with or even improve upon experimental data.
In addition to defining an ideal test case to extend the lattice-QCD calculation of elastic scattering to more complicated processes, phenomenological interest in γπ → ππ itself motivates a detailed study of the extrapolation to the physical point. At low energies, its form is dictated by the Wess-Zumino-Witten anomaly [19][20][21][22][23], leading to a theoretical prediction that has been tested experimentally at the 10 % level [24][25][26], including the study of higher-order chiral corrections [27][28][29][30][31] and dispersive techniques [26,29,32,33]to be contrasted with the π 0 → γγ anomaly, whose chiral prediction, F πγγ = 1/(4π 2 F π ) = 0.2745(3) GeV −1 , has been confronted with experiment at sub-percent precision, F πγγ = 0.2754(21) GeV −1 [34]. Furthermore, the process γπ → ππ provides input to the data-driven Standard-Model prediction of the anomalous magnetic moment of the muon a µ , constraining both hadronic vacuum polarization and hadronic-light-by-light scattering. In view of the current 4.2σ discrepancy between the resulting prediction   1 and experiment [78][79][80][81][82], further constraints on the hadronic amplitudes would of course be highly welcome. Phenomenologically, γπ → ππ is dominated by the ρ(770) resonance, thereby providing access to the radiative coupling of the ρ to a photon and a pion [33], with the full kinematic dependence required to improve vector-meson-dominance (VMD) approaches. In particular, the amplitudes that enter the hadronic contributions to a µ actually depend on the virtual process γ * π → ππ, an extension that automatically arises in lattice QCD. Equivalently, the process is related to the decay γ * → 3π via crossing symmetry, and thus connected to the vector-meson decays ω(782) → 3π and φ(1020) → 3π when the virtuality of the photon coincides with the respective mass. In this regard, the analysis of lattice-QCD data in the scattering region provides a testing ground for frameworks that aim to analyze three-particle scattering directly, a subject that is currently under intense investigation [83][84][85][86][87]. Accordingly, any model used to describe the γ ( * ) π → ππ lattice data needs to allow for a controlled extrapolation in the pion mass, work in the presence of a resonance, be accurate both in the low-energy region and in the complex plane where the ρ pole is located, and, ideally, respect crossing symmetry and allow for a description of the decay region at the same time. As such, chiral perturbation theory (ChPT) by itself is insufficient, as unitarity is only restored perturbatively and resonances thus cannot be produced without unitarization. Instead, here we use a dispersive approach, based on the fundamental principles of unitarity and analyticity, which are implemented in the so-called Khuri-Treiman (KT) equations [88]. Their solution defines a set of reliable amplitudes for γπ → ππ at the physical point [26,33], including the generalization to non-zero virtualities [89][90][91]. However, with the main input quantity the physical phase shifts for ππ scattering, this representation alone does not constrain the chiral extrapolation of lattice data, thus requiring the combination with ChPT. In particular, we use the inverse-amplitude method (IAM) [92][93][94][95][96][97][98][99] to describe ππ scattering at unphysical pion masses, as input for the solution of the KT equations away from the physical point, based on the implementation from Ref. [100]. In fact, the single-channel SU (2) IAM can again be justified via a dispersion relation, with the only approximation regarding the chiral expansion of the left-hand cut [101]. The resulting expression can therefore not only be used to describe ππ scattering on the real axis, but also to study the pion-mass dependence of resonance trajectories [102,103] or form factors [104,105]. Following the strategy already laid out in Ref. [106], the combined KT + IAM representation for γ ( * ) π → ππ is then fit to lattice-QCD data and afterwards extrapolated to the physical point, where the observables are extracted. At non-zero photon virtualities also the quark-mass dependence of the vector mesons ω and φ starts to play a role, see Ref. [107], which can again be constrained using ChPT arguments [108][109][110].
This paper is organized as follows. The basic form of the γ ( * ) π → ππ amplitude is introduced in Sec. 2. Subsequently, the dispersive framework used to analyze the lattice data and the fit procedure are discussed in Sec. 3 and Sec. 4, respectively. Finally, the fit results are presented in Sec. 5 and conclusions drawn in Sec. 6.
2 The process γπ → ππ The scattering amplitude M of the process γ ( * ) (q)π + (p) → π + (k)π 0 (k ) with four-momenta q, p, k, and k can be expressed in terms of the electromagnetic current J µ = e 2 3 uγ µ u − Here e is the elementary charge, µ is the polarization vector of the photon, s = (k + k ) 2 as well as t = (p − k) 2 denote the Mandelstam variables, and the pion states are normalized in the standard manner [111]. The pseudoscalar nature of the pions allows for decomposing the matrix element in terms of a complex-valued function F as At vanishing energy, the anomalous nature of the process fixes the amplitude in terms of the pion decay constant F π = 92.28(10) MeV [112] as F(0, 0, 0) = eF 3π with [21] In particular, we follow the convention to absorb the class of chiral corrections that corresponds to the quark-mass renormalization of the decay constants of the three external pions into the physical F π , and hence use Eq. (2.3) as the reference point. Only chiral corrections that go beyond the quark-mass renormalization of F π will thus be applied in the matching (3.14), see Refs. [26,27]. Throughout this work, isospin symmetry is assumed to hold. In this case, only odd partial waves contribute, leading to the expansion [113] F s, t, q 2 = ∞ j=0 f 2j+1 s, q 2 P 2j+1 (z). (2.4) Here f J denotes the partial wave of total angular momentum J = 2j + 1, z = cos θ with θ = ∠(k, k ) the scattering angle in the center-of-mass (CM) system, and P J the derivative of the J-th Legendre polynomial. Taking into account ππ intermediate states only, the P -wave fulfills the unitarity relation [89] Im f 1 s, the ππ → ππ P -wave amplitude and the ππ phase space, respectively, where δ = arg[T ] is the P -wave phase shift, M π the pion mass, and the cut of the square-root is chosen along the positive real axis, which leads to σ π (s * ) = −σ π (s) * . The partial wave T in turn obeys a slightly simpler unitarity relation, namely which is a special case of Watson's theorem [114]. Building upon the unitarity relation (2.5), an expression for f II 1 , the P -wave on the second Riemann sheet, can be derived. To that end, we make use of f 1 (s ± i , q 2 ) = f II 1 (s ∓ i , q 2 ) for s ≥ 4M 2 π and → 0, the Schwarz reflection principle, i.e., f 1 (s * ) = f 1 (s) * , as well as the uniqueness of analytic continuation to obtain .
(2.9) Equation (2.9) exhibits a pole at s ρ = (M ρ −iΓ ρ /2) 2 that is associated with the ρ resonance of mass M ρ and width Γ ρ and accompanied by a twin pole at s * ρ . The residue at s ρ factorizes into the coupling of the ρ to the final state, g ρππ , as well as to the radiative coupling g ργπ to the initial state as [33] res f II 1 , s ρ = −2eg ργπ g ρππ . (2.10) Hence to determine the radiative coupling, g ρππ needs to be known. The latter can be fixed via where T II (s) = T (s) 1 + 2iσ π (s)T (s) (2.12) is the ππ P -wave on its second Riemann sheet, whose form is again dictated by unitarity, i.e., Eq. (2.7). Note that to extract this couplings it is necessary to evaluate T in the complex plane, hence a parameterization of T is needed that goes beyond Eq. (2.6), which is valid only along the real axis above threshold. In lattice-QCD computations, the matrix element (2.2) is usually expressed in a different manner. To that end, the ππ final state is expanded into components |ππ, P, J, m of total angular momentum J and magnetic quantum number m. The P -wave contribution to the matrix element is subsequently decomposed as [17] ππ, P, 1, m|J µ (0)|π, p = e 2i M π µναβ p ν * α (m, P )P β A s, q 2 , (2.13) where the total momentum is given as P = k + k , * (m, P ) is the polarization vector of the two outgoing pions, i.e., the standard polarization vector of a spin 1 particle, and A is a complex-valued function. In Ref. [17] it is shown that A(s, q 2 ) ∝ k CM f 1 (s, q 2 ). Re-performing the computation, this time keeping track of all factors, results in with k CM the absolute value of the momentum of a final-state pion in the CM frame, i.e., s = 4(M 2 π + k 2 CM ). For the remainder of this work, we will ignore all partial waves with J ≥ 3. The on-shell cross section is then given as [17,26] 3 Dispersive representation of γπ → ππ In the energy region of interest, where the P -wave dominates, γ ( * ) π → ππ can be accurately described by the KT framework [88]. Building upon dispersion relations and elastic unitarity (2.5), the KT equations take into account not only individual ππ rescattering in the s-, t-, and u-channel, but also mixed rescattering, where pions rescatter, e.g., in the t-channel and subsequently in the s-channel. The starting point is the reconstruction theorem [29], which decomposes F into functions of a single Mandelstam variable only. Here u = 3M 2 π + q 2 − s − t is not a free variable. The KT equations then take the form [26,89,106] B s, q 2 = n−1 k=0 c k q 2 B k s, q 2 , Here n ∈ N is the number of subtractions that are employed in the dispersive integrals, c k are the subtraction functions, the mappings B k are known as basis functions, the Omnès function is given as [115] and t s, q 2 , z = τ s, q 2 + zκ s, q 2 , τ s, q 2 = 3M 2 π + q 2 − s 2 , κ s, q 2 = 1 2 σ π (s) λ(s, q 2 , M 2 π ), The basis functions subsume the ππ rescattering and are fixed as soon as the ππ phase δ is known. While the Omnès function describes ππ scattering in one channel, the integral in Eq. (3.2) incorporates mixed rescattering. That is, the replacement B k (s, q 2 ) → s k Ω(s) amounts to taking into account only ππ rescattering in the individual channels.
In the form of Eq. (3.2) the KT equations are valid only if q 2 < (3M π ) 2 , i.e., as long as the photon cannot decay. By deformation of either one of the integration contours the equations can be analytically continued towards q 2 > (3M π ) 2 [116,117], however, this is not needed for the lattice data of interest. This analytic continuation reveals that the basis functions indeed possess a three-particle cut in q 2 that is associated with pairwise ππ rescattering, but they do not contain any q 2 -dependence arising from genuine three-pion interactions [118]. Such interactions are to be described by the subtraction functions c k , which are not fixed by the KT approach. Thus, to arrive at a complete representation of γ ( * ) π → ππ, we need both a representation of δ as well as a parameterization of the subtraction functions.
For the former, we employ the IAM. In this approach, the ππ P -wave is expanded in SU (2) ChPT, T = T 2 + T 4 + . . . , where T 2 denotes the leading-order (LO) ChPT expression and T 4 the next-to-leading-order (NLO) one. This expansion satisfies Eq. (2.7) only perturbatively, but we can unitarize it to obtain which is precisely the NLO IAM [93][94][95]. Equation (3.5) satisfies Eq. (2.7) exactly, exhibits the correct analytic structure, and is valid in the entire complex plane. Explicit expressions for the ChPT amplitudes in closed analytical form are given, e.g., in Ref. [100], building upon the computations presented in Refs. [119,120]. In addition to the pion mass M π , the amplitudes depend on F , the pion decay constant in the chiral limit, as well a l r = l r 2 −2l r 1 , a single linear combination of the ordinarily renormalized low-energy constants (LECs) l r 1 , l r 2 . It is beneficial to express the amplitudes in terms of F instead of F π , for the former is pionmass independent, see also the discussion in Ref. [100]. In this work we use the N f = 2 + 1 FLAG average of F π /F [121][122][123][124][125][126], which yields F = 86.89 (58) MeV when combined with the PDG value of F π .
Unfortunately, the phase of the NLO IAM does not approach π, instead, lim s→∞ T (s) = −[96πl r + i + 2/(3π)] −1 , which yields lim s→∞ δ(s) < π for all reasonable values of l r . However, for the numerical computation of B k it is beneficial if δ approaches π at a finite value Λ of Mandelstam s, since this provides a natural cutoff of the integral in Eq. (3.2) via T (Λ) = 0. For this reason, and because the NLO IAM is physically reasonable in the elastic region only, we guide δ smoothly to π at energies far above the resonance region. To be precise, we use the IAM below s = 270M 2 π , δ(s) = π for s ≥ Λ = 310M 2 π , and a fourthorder polynomial in between such that the transition is smooth. The suppression of the high-energy region due to the subtractions ensures that the systematic error introduced in this way is negligible compared to the error of the γπ data, this suppression is particularly strong since we use not only one, but n = 2 subtractions, see Sec. 5. The basis functions can be computed numerically via standard methods [117].
The subtraction functions need to be holomorphic in the complex q 2 -plane except for a cut along [9M 2 π , ∞) that is associated with γ * → 3π. Hence we can write down an m-times subtracted dispersion relation of the form denotes the discontinuity along the branch cut. As long as q 2 < 9M 2 π , the Schwarz reflection principle dictates b kj ∈ R and disc[c k ] = 2iIm[c k ]. In the energy region that contributes most to the dispersive integral in Eq. (3.6) the three-pion physics is dominated by the ω(782) and φ(1020) resonances [112], both of which are narrow and (at the physical point) far away from the three-pion threshold. Thus disc[c k ] inside the integral can be reasonably well described by a sum of two Breit-Wigner functions, yielding a dispersively improved variant of a Breit-Wigner parameterization that ensures the correct analytic properties [91,127]. In practice, the lattice data we are going to analyze are obtained at M π > 300 MeV, at which mass the ω becomes a bound state [107]. Accordingly, instead of being incorporated into the dispersive integral, it appears as a pole at q 2 = M 2 ω . This can be taken into account by writing down a dispersion relation in the form of Eq. (3.6) for c k (q 2 )/P(q 2 ) with P(q 2 ) = (1 − q 2 /M 2 ω ) −1 and multiplying the result by the pole factor P.
Since the lattice data are obtained at virtualities significantly below the 3π threshold, Eq. (3.6) can be expanded as a Taylor series, keeping the first m terms yields However, the convergence of the Taylor series is poor as soon as q 2 gets close to the 3π threshold, this drawback goes hand in hand with a wrong asymptotic behavior for large q 2 , i.e., the expression diverges. To improve on Eq. (3.7), a conformal polynomial can be used instead [128]. That is, Eq. (3.6) is approximated by where the conformal variable w reads In this way, the cut along [9M 2 π , ∞) is retained, moreover, the asymptotic behavior is improved, as w is bounded.
At the 3π threshold, Im[c k ] should scale like (q 2 − 9M 2 π ) 4 to be in accordance with the three-particle phase space [129]. This is impossible to obtain with Eq.
for q 2 below threshold the Schwarz reflection principle needs to be fulfilled, hence b kj ∈ R.
Expanding w in powers of x := 9M 2 π − q 2 makes it clear that only odd powers of x contribute to Im[c k ]. This problem is fundamental to the method, for the Riemann mapping theorem implies that each biholomorphic map from the cut complex plane to the interior of the unit disc is of conformal form. It is however easily possible to remove the leading square-root-like scaling via fixing [42] (3.10) Altogether we have six different parameterizations of the subtraction functions: a polynomial, a conformal polynomial, and a conformal polynomial with modified threshold behavior, each either with or without the pole factor P in front. To take into account the pion-mass dependence of M ω appearing in P, we use the result of the analysis in Ref. [107], namely Lastly, the number of terms m needs to be fixed. We use N − k terms with a single global value of N for the k-th subtraction function c k , since it multiplies s k in Eq. (3.2), such that with our choice in the simple polynomial representation (3.7) the highest combined power of Mandelstam s and q 2 has mass dimension 2N . An overview of the different strategies is given in Table 1.
Given the basis functions, the P -wave can be expressed as Hence it requires the computation of B k . Doing so directly from its definition in Eq. (3.2) is in this context inconvenient for two reasons. First, the lattice data contain several data points at q 2 > M 2 π . At these virtualities, Mandelstam t develops a non-vanishing imaginary part due to the square root of the Källén function in Eq. (3.4). Thus B k needs to be evaluated at different lines in the complex plane, one for each q 2 > M 2 π . Second, we aim for an evaluation of f 1 at the resonance pole to extract the radiative coupling by means of Eq. (2.10). Again, this requires the computation of B k at complex values of Mandelstam t, where care needs to be taken to avoid collisions with the branch cut of B k . To circumvent these issues, we resort to the kernel method [33], which allows for evaluation of B k at arbitrary values of Mandelstam s by computing B k along the real axis only. The details are given in App. A.
Altogether, the free parameters of our dispersive representation are l r as well as the variables b kj appearing in Eq. (3.7) or Eq. (3.8). While the former is pion-mass independent, the latter depend on M π . Since there are only lattice data sets at two different pion masses available, only very simple parameterizations of the pion-mass dependence of each b kj can presently be constrained. For that reason, we opt for the simplest ansatz Here the fact that the variables are linear in M 2 π instead of M π is motivated by ChPT, for otherwise the b kj would possess branch points in the quark mass. As soon as more data sets at different pion masses become available, it will be possible to test more refined prescriptions.
With the subtraction functions extrapolated to the physical point, defined by the PDG value of the mass of the charged pion [112], the anomaly (2.3) can be determined via matching the dispersive representation to ChPT. For n = 2 subtractions the matching yields [26] determined via one-loop ChPT and a new LEC has been fixed via resonance saturation [27].

Fit to lattice data
The lattice-QCD computations of γ ( * ) π → ππ are based on the formalism presented in Ref. [130], which describes a two-step approach. First, N ππ different ππ CM energy levels E lat k , k = 1, . . . , N ππ are computed in the finite volume, which are related to the phase shift δ via Lüscher's quantization condition [1,131], (4.1) Here Z is a known expression that depends on the kinematics and the characteristics of the lattice. Second, a finite-volume version of the matrix element (2.13) is computed, from which one can extract A FV , the finite-volume analog of A in Eq. (2.13). The latter is related to its finite-volume counterpart by [17,130] A s, where the Lellouch-Lüscher factors are given as 2 These factors are uniquely defined only on the energies that are solutions of Eq. (4.1), that is, they are not defined as functions of arbitrary values of s [133]. Notably, the computation of L requires the computation of the derivative of the phase shift. Since the lattice data points are too sparse, it is not feasible to compute the derivative by an interpolation of the data. Instead, a continuous parameterization is needed, we use the NLO IAM as given in Eq. (3.5). Accordingly, the fit to the data works as follows. First, the NLO IAM is fit to the ππ energy levels by minimizing with respect to the fit parameter l r . Here C ππ is the covariance matrix of the lattice ππ energies E lat k , and E IAM k (l r ) is obtained by substituting δ in Eq. (4.1) with the IAM phase shift and solving the equation for √ s = E IAM k , with the kinematics of the kth lattice energy level, as explained in detail in Ref. [100].
Second, the derivative of the IAM phase for the resulting value of l r is used to compute the Lellouch-Lüscher factors. To be consistent, the factors are evaluated at the energies E IAM k , for these are the solutions of Lüscher's quantization condition with the IAM phase shift. Inserting the factors obtained this way into Eq. (4.2) and combining it with Eq. (2.14) allows us to compute are used to fit the P -wave as computed via KT equations. The lattice-QCD computation yields N γπ values of A FV , corresponding to N γπ different virtualities q 2 lat a , a = 1, . . . , N γπ , at N ππ different energies, all of which have errors. For the data sets at hand, N ππ < N γπ , i.e., several data points are obtained at the same energy. To take into account the errors of the energies and virtualities, we follow the standard approach and introduce an auxiliary fit parameter for each kinematic variable, see, e.g., Ref. [8], leading to with v lat = f lat  Table 2: The outcomes of the IAM fit to the ππ data. The first error arises due to the statistical error of the ππ data, the second due to the error of the lattice spacings, and the third due to the error of the literature value of F . The third and fourth column contain reference values for the LEC from ChPT and lattice QCD [121][122][123][124][125][126], respectively, while the fifth column lists the ρ properties as determined via Roy-like equations.
the auxiliary fit parameters q 2 1 , . . . , q 2 Nγπ , E 1 , . . . , E Nππ , as well as s k = E 2 k , and covariance matrix C γπ . The error of the IAM phase leads to an error of the Lellouch-Lüscher factors. The corresponding covariance matrix is added to the appropriate entries of C γπ . Equation (4.6) is minimized with respect to the auxiliary fit parameters and the variables b kj appearing in the parameterization of the subtraction functions. Since only the absolute value of the partial wave is fit and f 1 is linear in the fit parameters b kj , the latter are fixed by the fit only up to a global phase ±1. To fix this, we impose the upper case of Eq. (2.8), i.e., arg[f 1 (s, q 2 )] = δ(s).
There are different sources of error that need to be taken into account. The ππ energy levels E lat k carry an error due to the statistical nature of the lattice computation, this error is taken into account by jackknife resampling. Furthermore, on the lattice everything is computed in units of the lattice spacing a. The translation into physical units requires the determination of a, the so-called scale setting. The resulting value of a carries a statistical uncertainty, moreover, a systematic error arises for the scale setting is not unique away from the physical point. To keep the impact of the scale setting and the associated error of a as small as possible, we phrase both Eq. (4.4) and Eq. (4.6) in lattice units. However, in one place at the χ 2 -level the lattice spacing enters, namely via the decay constant F , whose literature value is required for the evaluation of the IAM amplitudes and needs to be translated into lattice units. To assess the impact of the statistical error of the lattice spacing on the ππ fit, we perform a parametric bootstrap, see Ref. [100] for further details. We do not attempt to estimate the uncertainty associated with the systematic error of the lattice spacing. To determine the error of the fit parameters of the γπ fit, we simply use the Hessian. In principle, the error of the IAM phase impacts the γπ fit not only via the covariance matrix C γπ , but also via the KT equations. In practice, the error of the phase is negligible compared to the error of |A FV |.    [136]). On the left-hand side, the phase is depicted. Here the markers encode the irreducible representations (irreps) of the residual rotational symmetry on the lattice, while the colors encode the square of the boost momentum d = L 2π P with L 3 the spatial volume of the lattice. On the right-hand side, the comparison is shown on the energy level, with the statistical uncertainty of the lattice data indicated by the light red boxes. All energies are given in lattice units.

Fits to ππ data
As already stressed, currently there are only two γπ data sets available, one at M π ≈ 317 MeV [18], the other one at M π ≈ 391 MeV [16,17]. Hence, to fix the pion-mass dependence, we need to analyze these two data sets simultaneously. To that end, we perform a combined fit of the NLO IAM to the ππ lattice data of Ref. [8] and Ref. [136], the former being associated with the 317 MeV γπ data and the latter with the 391 MeV ones. Since the two sets were independently generated, the χ 2 is the sum of two terms Ref. [18] Refs. [16,17]   in the form of Eq. (4.4), one for each data set. A graphical comparison of the fit result with the data is shown in Fig. 1, while the goodness of the fit is shown in Table 2 together with the obtained value of l r , the resulting p-value of 20 % being reasonable. There is a 2σ tension with the ChPT value of l r , however, this deviation comes at no surprise, given the unitarization via the IAM [95][96][97][98]. With l r fixed, we continue the ππ P -wave via Eq. (2.12) to the second Riemann sheet to determine the ρ characteristics at the physical point as listed in Table 2. Comparing with the literature values given ibidem, we note a 4σ discrepancy in M ρ and a 2σ tension in Im(g ρππ ), while both the width and the real part of the coupling agree well. This is explained by the fact that the NLO IAM has only a single free parameter, leading to a trade-off between the different ρ properties. To improve on this, the next-to-next-to-leading-order (NNLO) IAM can be employed [100], however, with data at only two different pion masses both exceeding 300 MeV, we find that stable fits are not feasible.

Fits to γπ data
Next, we fit the γπ data. Since the pion-mass dependence of each fit parameter b kj is described by two free parameters, compare Eq. (3.13), we can perform the fits to the two γπ data sets independently, the fit parameters being the values of b kj at the two different pion masses. Hence at this stage we can work in lattice units, with different units for each data set. We use n = 2 subtractions in Eq. (3.12), for once subtracted KT equations fail to describe the energy dependence of the data correctly and thus do not allow for statistically acceptable fits. Note that increasing the number of subtractions to n = 3 does not provide additional flexibility, since the reconstruction theorem (3.1) is invariant under the shift B(s, q 2 ) → B(s, q 2 ) + λ(q 2 )(3s − 3M 2 π − q 2 ) with λ an arbitrary function, hence one subtraction function can be eliminated. This shift is forbidden for n = 2 due to the high-energy behavior of B, but becomes possible for n = 3. In addition, we pick N = 2, that is, we have three fit parameters b kj in c 0 and two in c 1 . If instead N = 1 is used, the fit quality becomes poor, while at N = 3 the fit stability deteriorates. The exception are the strategies III and IIIP, where we pick N = 3, which again amounts to five fit parameters due to Eq. (3.10).
To obtain statistically acceptable fits to the data at M π ≈ 391 MeV, we need to exclude the six data points at the highest energy, E lat ≈ 1096 MeV. These points lie far above the resonance region, for although at this pion mass the ρ is heavy, i.e., M ρ = 846.1(3.1)(3.2)(0.1) MeV (errors as in Table 2), its width Γ ρ = 10.8(8)(9)(1) MeV is tiny. Moreover, several of the data points with the smallest absolute errors of |A FV | are located at this energy. Hence the six excluded data points provide rather strong constraints on the asymptotic high-energy behavior of the KT equations instead of the resonance physics in which we are primarily interested.
We carry out fits for each strategy enumerated in Table 1, with an overview of the fit qualities given in Table 3. While the goodness of the fit at the lower pion mass is rather insensitive to the parameterization of the subtraction functions, the data at the higher pion mass is more selective, because the relative error of |A FV | at the higher pion mass is smaller than the error at lower mass. Notably, we observe improvement when including a pole factor, this is true for all strategies, with the overall p-value improving by at least an order of magnitude in each case, and even by two orders of magnitude when going from strategy I to IP. As soon as a pole factor is included, it does not matter much if the remaining q 2 -dependence is parameterized by a plain polynomial, a conformal one, or a conformal one with modified threshold behavior, the overall p-values of strategy IP, IIP, and IIIP are similar, with a slight improvement when using conformal parameterizations. Hence in the following we group the results of the three parameterizations including a pole together. If no pole is used, at higher pion mass the fit clearly disfavors a plain polynomial and instead prefers a conformal one, with only a very slight further improvement when modifying the threshold scaling. Thus we exclude strategy I and combine strategy II and III. As a representative of each group, we pick strategy II and IIP. The corresponding partial waves are compared with the two lattice data sets in Fig. 2 and Fig. 3. As can be observed, independently of the presence of a pole factor, the magnitude of f 1 increases with growing q 2 , in accordance with phenomenology [91].
To check if we are sensitive to the mixed rescattering effects included in the KT equations, we re-perform the fits with the replacement B k (s, q 2 ) → s k Ω(s). At M π ≈ 317 MeV we obtain a p-value of 4.97 × 10 −2 with strategy II and 6.26 × 10 −2 with strategy IIP, while at M π ≈ 391 MeV we obtain 8.77 × 10 −3 and 8.49 × 10 −2 , respectively. Comparing with the corresponding entries of Table 3, the observed difference is insignificant, thus we conclude that mixed rescattering does not need to be taken into account to describe the data at the present level of precision. The results of two fit strategies in comparison with the γπ lattice data of Ref. [18] at M π ≈ 317 MeV. Shown are slices of constant energy. For convenience, the results are displayed in physical units, but the fit is carried out in lattice units, thus the error bands represent the statistical error only.

Chiral extrapolation
Equipped with the KT fit results, we can determine the pion-mass dependence of the fit parameters via Eq. (3.13). To that end, we need to translate the fit parameters associated with the two different data sets to a common set of units, hence the errors of the lattice strategy II strategy IIP Figure 3: As Fig. 2, but for the lattice data of Refs. [16,17] at M π ≈ 391 MeV.
spacings enter the picture. Since the scale of the two data sets is set in different ways (via the Υ(1S)-Υ(2S) splitting at M π ≈ 317 MeV and via the Ω baryon mass at M π ≈ 391 MeV), an additional systematic error arises, which is difficult to quantify (this also applies to the fit of the ππ data). However, compared to the sizable statistical uncertainty of the data and the systematic error of the chiral extrapolation to be discussed in Sec. 5.4, the systematic error associated with the scale setting is likely irrelevant at present. Therefore, the uncertainty associated with the lattice spacing given in the remainder of this work will always refer to its statistical error only. 3 The pion-mass dependence of the fit parameters is depicted in Fig. 4. While the leading parameters in the series expansion (3.8), b 00 and b 10 , are constrained more strongly by the data at lower pion mass than the one at higher pion mass, the opposite is true for the highest-order term associated with b 02 . The latter comes at no surprise, for the data at higher mass contain much larger virtualities in the spacelike region, exceeding in absolute value the timelike virtualities of both data sets significantly, giving thus more weight to the b 02 term. With decreasing pion mass, the ω pole moves from the real axis below the 3π threshold on the first Riemann sheet into the complex plane on the second sheet. Since strategy II strategy IIP Figure 4: The pion-mass dependence of the fit parameters for two different fit strategies. The dashed gray lines mark the physical pion mass and the ones of the two lattice data sets, the dashed error bands correspond to the error of the lattice spacings, while the filled ones are associated with the statistical error.
the pole factor P that is present in strategies IP, IIP, and IIIP describes a bound state, naturally the question arises if the change in the nature of the pole needs to be reflected in the extrapolation in the pion mass for these strategies. A resonant ω could be implemented via a dispersively-improved Breit-Wigner parameterization, which, in practice, is almost indistinguishable from a pole ansatz unless very close to the singularities. Given the large uncertainties of b jk at the physical point, this change is thus immaterial, especially, since for the extraction of the observables the subtraction functions need to be evaluated at vanishing virtuality only, and for every strategy c k (0) = b k0 holds. Taking care of the pion-mass dependence of the fit parameters, the IAM, and the KT equations, we can extrapolate the partial wave to the physical point. Computing the cross section via Eq. (2.15) yields the line shape shown in Fig. 5 exhibiting the characteristic resonance peak. In both fit strategies, the error increases when moving beyond the resonance, which reflects the fact that most data points lie around the resonance region. In principle, the omitted data points at the highest energy could provide further constraints, but since no acceptable fits could be found when including these points, we conclude that with the currently available lattice data the asymptotic form of the cross section remains largely unconstrained. In this regard, we remark that the KT basis functions with n = 2 subtractions increase too fast asymptotically compared to expectations from the Froissart bound [137], so that a proper high-energy completion needs to be imposed [33]. However, these considerations become relevant only well beyond 1 GeV and thus do not affect the current fit, for which the n = 2 subtraction scheme provides the adequate number of free parameters to be able to describe both the chiral anomaly and the ρ-meson properties [33].

Chiral anomaly and radiative coupling
Finally, we can determine the anomaly F 3π and the radiative coupling at the physical point via Eq. (3.14) and Eq. (2.10), respectively. The values are listed for the different fit strategies in Table 4. Since the outcomes of the different fit variants are highly correlated with only minor differences in fit quality and very similar statistical errors, we do not compute weighted averages, but instead only perform plain averages to determine the central values. Doing so for the acceptable fits without a pole factor, i.e., averaging over strategy II and F 3π × GeV 3 Re(g ργπ ) × GeV Im(g ργπ ) × GeV   with errors as in Table 4, while the strategies including an ω pole, i.e., IP, IIP, and IIIP, yield

(5.2)
Both values of F 3π are compatible with the prediction (2.3), albeit only due to their large errors. Fits including the pole ansatz do display a better fit quality, but not at a level that would conclusively demonstrate the necessity of the pole. Since, further, both fit variants agree within statistical uncertainties, we conclude that the current lattice data cannot discriminate between Eq. (5.1) and Eq. (5.2) and quote the resulting spread as an additional systematic error. This error also arises due to the absence of lattice data at several different pion masses by one collaboration, forcing us to fit our representation to two data sets by two different collaborations at only two different pion masses, which makes it impossible to fix the pion-mass dependence of the subtraction functions beyond the simple ansatz (3.13). Averaging over all fit results except for strategy I, we finally quote where the last error is our estimate of the systematic uncertainty associated with the parameterization of the subtraction functions. The resulting value of F 3π is perfectly consistent with the chiral prediction (2.3), but carries a large uncertainty. This is the first extraction of this low-energy parameter from lattice-QCD calculations, and will improve accordingly once better data become available. The residue g ργπ is currently not known better than from an SU (3) VMD estimate [138], which suggests |g ργπ | = 0.79(8)GeV −1 [33], again compatible with Eq. (5.3) (within 1.2σ). 4 The difference to the VMD estimate increases to 2.3σ for Eq. (5.1), while there is full agreement with Eq. (5.2). This provides a-posteriori evidence for the presence of an ω pole in the subtraction functions, as does the final result for the cross section shown in Fig. 5 when compared to the expected peak cross section around 20 µb [144]. The radiative coupling has also been extracted in Ref. [18] under the assumption that the pion-mass dependence of |G ργπ | = |g ργπ |M π /2 is weak, leading to |g ργπ | [18] = 1.15(5)(3) GeV −1 . This value differs from the VMD estimate by 3.6σ, a discrepancy that went unnoticed in Ref. [18] because it is mitigated by a missing factor 2 in Eq. (17) for Γ(ρ → πγ) therein [145]. Moreover, our analysis shows that the uncertainties especially from the chiral extrapolations are substantially larger. In particular, a pion-mass independent |G ργπ | renders the residue divergent in the chiral limit, while at M π = 317 MeV one has |g ργπ | Mπ=317 MeV [18] = 0.507 (20)(13) GeV −1 as well as |g ργπ | Mπ=317 MeV = 0.552 (18)(18)(0), the latter being the average (5.3) at this pion mass. We conclude that |g ργπ | instead of |G ργπ | is approximately pion-mass independent, thus avoiding the divergence in the chiral limit.

Conclusions
In this work we analyzed state-of-the-art γπ → ππ lattice-QCD data using a combination of dispersion relations and ChPT, to be able to describe both the momentum and the pionmass dependence in a reliable manner. Extrapolating to the physical point, we determined the cross section, extracted the radiative coupling of the ρ meson, and, for the first time, the chiral anomaly F 3π , see Eq. (5.3) for the final results, and Eq. (5.2) for a variant that imposes the ω pole in the subtraction functions. These results agree with expectations from ChPT and phenomenology, albeit within large uncertainties. By combining KT equations with the pion-mass dependence as described by the IAM, we could thus confront predictions from the KT framework with lattice QCD, emphasizing the role of the reaction γπ → ππ to develop methods that could subsequently be applied to more complicated processes. While the current lattice data are not yet sensitive to the mixed rescattering effects included in KT equations, this work demonstrates that the framework is not only a valuable tool for the description of experimental data, but that it also applies to the analysis of lattice calculations.
Future lattice-QCD computations are expected to reduce the statistical uncertainties, hence the data will be able to differentiate between parameterizations of the pion-mass dependence of the subtraction functions, which currently represents the largest source of systematic uncertainties, e.g., with the presence of an ω pole in the subtraction functions preferred by comparison to phenomenology, but not resolved within the lattice calcula- 4 The branching fractions cited in Ref. [112] imply |gργπ| = 0.72(4) GeV −1 for the charged channel and |gργπ| = 0.73(6) GeV −1 for the neutral one. However, these values derive from high-energy Primakoff measurements [139][140][141] and VMD fits to e + e − → π 0 γ data [142,143], respectively, and thus involve a substantial model dependence.
tions. Moreover, once data at more different pion masses become available, more refined parameterizations can be employed, and if such data are analyzed with a single scale-setting strategy, this will remove the systematic error stemming from the simultaneous use of different ways to set the scale. Such refined lattice-QCD calculations, in combination with the analysis tools developed in this work, have the potential to improve several important lowenergy parameters of QCD, most notably the chiral anomaly F 3π , complementary to future experimental determinations, e.g., from the COMPASS Primakoff program [144]. We also expect that the strategy pursued here, combining dispersion relations with effective-fieldtheory methods for the pion-mass dependence, will have broad applications in particular to other processes dominated by ππ dynamics.

Acknowledgments
The lattice data taken from Refs. [16,17,136] were provided by the Hadron Spectrum Collaboration (HadSpec)-no endorsement on their part of the analysis presented in the current paper should be assumed. In addition, we would like to thank the authors of Refs. [8,18] for sharing lattice data with us. We thank Marcus Petschlies and Raúl Briceño for valuable discussions, and Christopher Thomas for comments on the manuscript. Financial support by the Bonn-Cologne Graduate School of Physics and Astronomy (BCGS), the DFG (CRC 110, "Symmetries and the Emergence of Structure in QCD"), and the Swiss National Science Foundation, under Project No. PCEFP2_181117, is gratefully acknowledged.

A Kernel method
Here we describe a way to facilitate the computation of the hat function B k as defined in Eq. (3.2) and needed in the evaluation of the partial wave via Eq. (3.12). To that end, we extend the discussion given in Refs. [33,146] to non-vanishing virtualities. The starting point is the simple dispersive representation B k s, q 2 = P k (s) + s d π ∞ 4M 2 π Im B k x, q 2 (x − s)x d dx, (A.1) with d the number of subtractions and P k = d−1 j=0 h kj s j the subtraction polynomial whose real coefficients h kj can be determined via matching to Eq. For low d the kernel W d simplifies to [33,146] W 1 s, q 2 , x = W s, q 2 , x − 2 x , W 2 s, q 2 , x = W 1 s, q 2 , x − 2τ s, q 2 x 2 , (A.6) with W s, q 2 , x = 3 2κ(s, q 2 ) 1 −1 1 − y 2 ξ(s, q 2 , x) − y dy = 3 κ(s, q 2 ) 1 − ξ s, q 2 , x 2 Q 0 ξ s, q 2 , x + ξ s, q 2 , x , where Q 0 (ξ) = 1 2 1 −1 1 ξ − y dy = 1 2 log 1 + ξ 1 − ξ − iπ sgn(Im(ξ)) (A.8) is the lowest Legendre function of the second kind. Here we use the principal branch of the logarithm with a cut along the negative real axis. If s = 4M 2 π or, in case of timelike virtualities, s = ( q 2 + M π ) 2 , κ(s, q 2 ) vanishes. Equation (A.7) implies that W has singularities whenever this happens. To prove that these are removable and to derive a representation that is suitable for numerical implementation, we note that κ → 0 implies ξ → ∞, hence we expand Q 0 around ξ −1 = 0, arriving at W s, q 2 , x = 6 x − τ (s) ∞ j=0 1 (2j + 1)(2j + 3) 1 ξ(s, q 2 , x) 2j . (A.9) According to the ratio test, Eq. (A.9) converges absolutely as long as |ξ| > 1. In particular, we see that lim ξ→∞ W is manifestly finite. In passing, we note that the terms Because the high-energy region in the integral in Eq. (A.10) is strongly suppressed by the asymptotic behavior of W 2 , the integral can be cut off at high energies.