The $\pi\gamma \to \pi\pi$ transition and the $\rho$ radiative decay width from lattice QCD

We report a lattice QCD determination of the $\pi\gamma \to \pi\pi$ transition amplitude for the $P$-wave, $I=1$ two-pion final state, as a function of the photon virtuality and $\pi\pi$ invariant mass. The calculation was performed with $2+1$ flavors of clover fermions at a pion mass of approximately $320$ MeV, on a $32^3 \times 96$ lattice with $L\approx 3.6$ fm. We construct the necessary correlation functions using a combination of smeared forward, sequential and stochastic propagators, and determine the finite-volume matrix elements for all $\pi\pi$ momenta up to $|\vec{P}|= \sqrt{3} \frac{2\pi}{L}$ and all associated irreducible representations. In the mapping of the finite-volume to infinite-volume matrix elements using the Lellouch-L\"uscher factor, we consider two different parametrizations of the $\pi\pi$ scattering phase shift. We fit the $q^2$ and $s$ dependence of the infinite-volume transition amplitude in a model-independent way using series expansions, and compare multiple different truncations of this series. Through analytic continuation to the $\rho$ resonance pole, we also determine the $\pi\gamma \to \rho$ resonant transition form factor and the $\rho$ meson photocoupling, and obtain $|G_{\rho\pi\gamma}| = 0.0802(32)(20)$.

The first numerical calculations involving the Lellouch-Lüscher formalism were performed for K → ππ, providing an ab-initio Standard-Model prediction of direct * leskovec@email.arizona.edu † smeinel@email.arizona.edu ‡ marcus.petschlies@hiskp.uni-bonn.de CP violation in this process [40][41][42]. More recently, the generalization of the formalism by Briceño, Hansen, and Walker-Loud (BHWL) [31] was applied by the Hadron Spectrum Collaboration to compute the πγ → ππ amplitude, with the ππ system in a P -wave, as a function of photon virtuality and ππ invariant mass [43,44]. This amplitude describes ρ photoproduction and radiative decay [45,46], and also plays an important role in dispersion relations used to calculate the hadronic contributions to the anomalous magnetic moment of the muon [47][48][49][50]. Various theoretical aspects of the πγ → ππ process have also been discussed in Refs. [51][52][53][54][55][56][57]. As far as the finitevolume formalism is concerned, the πγ → ππ amplitude in the ρ resonance region is one of the simplest 1 → 2 processes to study on the lattice, because the ππ scattering is almost completely elastic in the relevant energy region.
In this paper, we report a lattice QCD calculation of the πγ → ππ transition with 2 + 1 flavors of cloverimproved Wilson fermions [58] at a pion mass of approximately 320 MeV, building upon our previous work on ππ scattering [28]. In contrast to the original Lellouch-Lüscher approach to the nonleptonic K → ππ decay, where the lattice parameters need to be tuned such that the final and initial hadronic states have equal energy, the BHWL formalism enables us to obtain the πγ → ππ amplitude for all ππ energy levels and arbitrary momentum transfer.
In Sec. II, we discuss the πγ → ππ amplitude and related quantities in the continuum. The parameters of our lattice calculation are given in Sec. III. We describe the interpolating fields and correlation functions in Sec. IV, and the extraction of the finite-volume matrix elements from these correlation functions in Sec. V. The mapping from finite volume to infinite volume using the Lellouch-Lüscher factor is explained in Sec. VI. We carefully study a model-independent approach for parametrizing the q 2 and s dependence of the πγ → ππ amplitude in Sec. VII, and present our results for the πγ → ππ cross section, the πγ → ρ resonant form factor, the ρ meson photocoupling, and the ρ radiative decay width in Sec. VIII.

II. ABOUT THE πγ → ππ PROCESS
The resonance photoproduction process πγ → ρ is obtained from the more general process πγ → ππ, where the final ππ state is in P -wave and couples strongly to the ρ resonance with isospin I = 1, I 3 = 1 and J P C = 1 −− . Throughout this paper (except where stated otherwise), we allow the photon to be virtual, but continue to denote it as just γ. The ππ photoproduction is described by the continuum infinite-volume matrix element ππ|J µ (0)|π , which is constructed from the initial state |π , the insertion of the QED current J µ (defined without the factor of e) and the final state |ππ with I = 1 and fourmomentum P = ( s + P 2 , P ). The latter is projected to the P -wave, so that it couples to the ρ resonance, where the polarization of the system is described by ν (P, m) [59]. Due to the Lorentz symmetry the matrix element decomposes like ππ|J µ (0)|π = 2iV πγ→ππ (q 2 , s) m π νµαβ ν (P, m)(p π ) α P β , where q = p π − P is the photon four-momentum transfer. Above, the current is taken in position space, and the single-pion state is normalized as π, p π |π, p π = 2E pπ π (2π) 3 δ 3 ( p π − p π ).
The P -wave two-pion states with polarization m are given by |ππ, √ s, P , 1, m where |ππ, √ s, P , k cm is a two-pion state with total momentum P , relative momentum direction unit vector k cm in the center-of-momentum frame, and invariant mass √ s. These states are normalized according to ππ, √ s , P , k cm |ππ, √ s, P , k cm where E 1 and E 2 are the individual pion energies, These normalizations of states imply that the matrix element (1) is dimensionless and that V πγ→ππ has units of MeV −1 . Notice that there is no explicit ρ label in the amplitude; this is because the ρ is not a QCD asymptotic state, but rather a resonance in P -wave ππ scattering with I = 1 associated with the pole in the scattering amplitude T ππ→ππ at s P ≈ m 2 R − im R Γ R . The transition amplitude V πγ→ππ depends on both the photon four-momentum transfer q 2 and the ππ invariant mass s. Like T ππ→ππ , this amplitude also has a pole at s = s P ; the residue at the pole gives the ρ resonance photoproduction form factor. For s in the vicinity of s P and at q 2 = 0, the amplitudes T ππ→ππ and V πγ→ππ behave like [32] T ππ→ππ (s) ∼ G 2 ρππ s P − s , where G ρππ and G ρπγ are the couplings of the ρ resonance to ππ and πγ, respectively. The ππ elastic scattering amplitude is related to the scattering phase shift δ(s) via where k is the scattering momentum, defined by √ s = 2 m 2 π + k 2 . Near a narrow resonance, the phase shift is well described by parametrizations of the Breit-Wigner type, where multiple different choices can be used for Γ(s). Inserting Eq. (10) into (9) gives Motivated by Eqs. (7) and (8), we write the photoproduction amplitude V πγ→ππ (q 2 , s) as where the form factor F (q 2 , s) no longer has a pole in s, and becomes equal to the photocoupling G ρπγ for s = m 2 R − im R Γ R and q 2 = 0. More generally, we define the resonant form factor for arbitrary photon virtuality as Note that Eq. (12) explicitly satisfies Watson's theorem. In Ref. [28] we found that our ππ scattering amplitude is well described by the BWI and BWII Breit-Wigner models discussed in Sec. II of that same reference, so we will continue to utilize the Breit-Wigner formulas throughout this work. The nonresonant backgrounds were found to be consistent with zero and are not included in the ππ scattering amplitude here. For convenience, we repeat the definitions of BWI and BWII here: • BW I: where g ρππ is the coupling between the ππ scattering channel and the ρ resonance in the Breit-Wigner model.
• BW II: where k R is the scattering momentum at √ s = m R and r 0 is the radius of the centrifugal barrier [60].
We consider two physically observable quantities we can determine from | ππ|J µ (0)|π |. The first is the π + γ → π + π 0 cross section as a function of π + π 0 invariant mass, which in the center-of-momentum frame is given by [43] σ(π + γ → π + π 0 ; s, The cross section can be measured at q 2 = 0, i.e., with a real photon. In Eq. (16), V πγ→ππ denotes the π + γ → π + π 0 amplitude with the unsymmetrized π + π 0 final state. Because this amplitude is related to the amplitude with the isospin-projected final state through A second physically observable quantity is related to the ρ resonance, which appears in the ππ system. The ρ radiative decay width Γ(ρ → πγ) is determined by the photocoupling Note that the form factor F (q 2 , s), and hence the photocoupling G ρπγ , do not depend on whether one considers the unsymmetrized or the isospin-projected twopion state. In the unsymmetrized case, Eqs. (9) and (12) would read [31] with the same F (q 2 , s).

III. LATTICE PARAMETERS
This calculation is performed on a single ensemble of gauge-field configurations with 2 + 1 flavors of dynamical clover fermions. This is the same ensemble as used in our calculation of ππ scattering [28], and we refer the reader to that reference for further details. The main parameters are summarized in Table I. The strange-quark mass is consistent with its physical value as determined via the "η s " mass [62,63]. The lattice scale was determined from the Υ (1S)-(2S) splitting [62,64], where NRQCD [65] with the physical b-quark mass was used to calculate the masses. The renormalization factor Z V of the local vector current was determined by the LHPC Collaboration as explained in Ref. [66].

IV. INTERPOLATING FIELDS AND CORRELATION FUNCTIONS
To determine the finite-volume matrix elements we are interested in, we need to compute two-point functions for the single-pion system (J P C = 0 −+ , I = 1, I 3 = 1) and for the two-pion system (J P C = 1 −− , I = 1, I 3 = 1), as well as three-point functions with an insertion of the electromagnetic current. The generalized eigenvectors obtained in the spectroscopic analysis of the two-point functions are then used to construct optimized threepoint functions.
A. Two-point functions overlapping with J P C = 0 −+ The projection of the single-pion field to an irreducible representation is trivial, i.e., it resides in the (pseu-  irreducible representation [4] for all momenta. For clarity we suppress the group indices of the single-pion field. We use the following interpolating operator: with momentum p π . The associated correlator C pπ π is The ground-state contribution to the pion correlator, which is obtained in the limit of large t, has the decomposition where the overlap factor is defined as and the finite-volume states are normalized such that π, p π |π, p π F V = 2E pπ π δ pπ, p π .
Because the pion is a stable hadron, its energy is affected only by exponentially suppressed finite-volume effects, which are negligible for our value of m π L. The dispersion relation of the pion was presented in Ref. [28] and follows the relativistic form well.
B. Two-point functions overlapping with The J P C = 1 −− two-point functions with momentum P are constructed using two types of interpolators, the single-hadron and the multi-hadron interpolators: where p 1 + p 2 = P . To project these interpolators to definite irreps Λ of the little group LG( P ), we use the projection formulas with representation matrices [67] O P , Λ, r and where p takes on the values Above, dim(Λ) is the dimension of the irrep, N LG( P ) is the order of the Little Group LG( P ), and Γ Λ rr (R) are suitably chosen representation matrices ofR ∈ LG( P ). In our choice of basis indexing and projecting to finitevolume irreps, we use the x, y, z polarization indices and not the helicity basis like in Refs. [68,69].
In the following, we jointly denote the projected interpolators as where the index i labels the type: i = 1 and i = 2 correspond to the quark-antiquark interpolators with Γ i = γ i and Γ i = γ 0 γ i , respectively, and i = 3, 4 correspond to two-pion interpolators with different values of | p 1 | and | p 2 |.
From the interpolators O P , Λ, r i we calculate a correla- . Its construction in terms of forward, sequential, and stochastic quark propagators is discussed in Ref. [28]. The spectral decomposition of the two-point correlation matrix reads 5 As in Sec. IV A, we define the overlap factors as and the finite-volume states are normalized such that n , P , Λ , r |n, P , Λ, r F V = 2E P ,Λ n δ n,n δ P , P δ Λ,Λ δ r,r .
To extract the energies and overlap factors from the correlation matrix, we use the variational analysis [70][71][72][73] by solving the generalized eigenvalue problem where we fix the normalization to v n, P ,Λ † Throughout this paper, we use the summation convention for repeated indices i or j. The principal correlators asymptotically behave as and we use single-exponential fits to extract E P ,Λ n [28].
The generalized eigenvector v n, P , Λ i (t 0 ) can also be used to construct the optimized interpolator [70][71][72][73] which has a dominant overlap to a single well-defined state labeled with (n, P , Λ, r). Note that, although we perform the analysis independently for different rows r of the irrep Λ, in the infinite-statistics limit the energies and eigenvectors are independent of r. An important quantity related to the energy E P ,Λ n of the ππ system in the moving frame P is the invariant mass which is also used to define the scattering momentum k P ,Λ n via s P ,Λ n = 2 m 2 π + (k P ,Λ n ) 2 .
C. The three-point functions The current insertion that represents the interactions between the photon and the hadrons depends on the photon momentum q, which combined with the initial and final state momenta satisfies momentum conservation: P + q − p π = 0. For the current insertion operator we use with the local current The renormalization coefficient Z V was determined in Ref. [66] and is listed in Table I. The three-point correlation functions are then obtained from the sink/source interpolators and current insertion as where t ππ is the source time, t J is the current insertion time and t π is the sink time. The three-point function is expressed in terms of quark propagators by evaluating Wick contractions. Figure 1 shows the quark-flow diagrams needed to calculate the C pπ, P , Λ, r 3, µ,i three-point functions. The current-disconnected diagrams labeled (a) and (b), i.e. the diagrams where the quark flow goes from the current J µ directly back to the current J µ , are omitted in this study. For the nucleon electromagnetic form factors, the current-disconnected contributions are known to be of order 1% for the quark masses used here [74].
The Wick contractions depicted in Fig. 1 are constructed from point-to-all, sequential and stochastic time-slice propagators. The technique builds upon and extends the scheme used in Ref. [28] for the construction of two-point correlation functions. This combination of propagator types allows for a compromise in flexibility to construct all required diagrams, minimal input of stochastic noise into correlation functions and economy in the cost of producing quark propagators and contractions.
The point-to-all propagators are obtained from the inversion of the Dirac operator D on a fully spin-and colordiluted point-source localized at fixed source location y, The sequential propagator results from performing an additional inversion of the Dirac operator on a pointto-all propagator on sequential source time slice t with insertion of a spin matrix Γ and momentum p according to T (x; t, Γ, p; y) ab αβ = z S(x; t, z) ac αγ Γ γδ e i p· z S(t, z; y) cb δβ .
For the purpose of this work the sequential source time slice always coincides with the source time slice , Finally, the stochastic propagators follow from inverting the Dirac matrix on stochastic time slice sources, whose components on a fixed time slice t y are independently and identically set with Z 2 + i Z 2 noise, such that we have the expectation values As a variant of the stochastic time-slice propagator defined in Eq. (43) the one-end-trick based on spin diluted stochastic time-slice sources is used to construct diagrams (e) and (f) in Fig. 1. The sources in Eq. (44) are thus modified according to The one-end-trick then allows for the representation of a product of quark propagatos by two stochastic propagators through a vertex given again by Γ and p as = S(x; t y , y) ac ακ e i p· y (Γ γ 5 ) κλ S(z; t y , y) bc * β λ (γ 5 ) β β = S(x; t y , y) e i p· y Γ S(t y , y; z) ab αβ .
The quark propagator loops of the connected diagrams (c) and (d) are closed using the stochastic time-slice propagator from current vertex J µ to pion vertexd γ 5 u at sink. Based on the application of point-to-all and stochastic propagator these diagrams are factorized into elementary contractions. For diagram (c), we have where η φ,ξ are contracted, Fourier transformed and stored separately as η φ (t J , q) and η ξ (t π , p π ) for each stochastic sample. Subsequently they are used to recombine the diagram for all required momenta q, p π as well as any vertex Γ i and P at the source. Diagram (d) follows analogously by promoting the point-to-all propagator S(x f ; x i ) in Eq. (50) to a sequential propagator, For diagram (e) in Fig. 1, the one-end-trick setup in Eq. (49) leads to the factorization of the diagram, Finally, diagram (f) is calculated as the product of propagator loop traces using again the one-end-trick, All quark propagators are smeared at their source and sink side in the same way as in Ref. [28], except the end of propagators joining the local current insertion vertex.

D. Optimized three-point functions
The spectral decomposition of the three-point function C pπ, P , Λ, r 3, µ,i (t π , t J , t ππ ), keeping as before only the groundstate contribution for the pion (for large t π − t J ), is For the ππ system we want to project to the n-th state. This will allow us to have a definite invariant mass, s P ,Λ n , and momentum transfer, (q 2 ) pπ n, P ,Λ = (E P ,Λ n − E pπ π ) 2 − q 2 , in our matrix element. To achieve this we utilize the orthogonality between the generalized eigenvectors and overlap factors 1 , and construct the optimized three-point function [75][76][77] Ω pπ, P , Λ, r This gives and we see that the optimized three-point function overlaps only to the single definite state |n; P , Λ, r .

V. DETERMINING THE FINITE-VOLUME MATRIX ELEMENTS
To extract the finite-volume matrix elements π, p π |J µ (0, q)|n, P , Λ, r F V from the correlation functions, we construct the ratio where C pπ π is the pion correlator, λ P ,Λ n is the principal correlator of the variational analysis, ∆t = t π − t ππ is the source-sink separation, and t = t ππ + t π − t J . The t 0 dependence of the optimized three-point function cancels with the t 0 dependence of the principal correlator. Inserting Eq. (57) into Eq. (58) gives (for large time separations) The matrix elements determined from Eq. (58) still contain residual excited-state contamination that decays exponentially for large ∆t, t J − t ππ , and t π − t J . We have data for ∆t/a = 8, 10, 12. There are several ways to proceed from this point on: 1) Set t J −t ππ = ∆t/2 and fit only the ∆t dependence of the matrix element with an excited-state model, as for example in Ref. [78], 2) Fit both the ∆t and t J − t ππ dependence with an excited-state model, 3) Fit constants to the ratios (assuming that only the desired initial and final states contribute), varying the time ranges to assess residual contamination.
We found that that options 1) and 2) did not yield stable fits, because we have too few source-sink separations and the statistical uncertainties are too large. We therefore use option 3), where we investigate whether the various fits are statistically compatible, and estimate a systematic uncertainty associated with the fit choice. In Fig. 2 we present results for the matrix elements | π, p π |J µ (0, q)|n, P , Λ, r F V | at representative kinematic points (plots for the other kinematic points are shown in Appendix A). As explained in the caption of the figure, we perform fits for many different time ranges and then choose one that appears to have plateaued for the further analysis. To estimate the systematic uncertainty associated with the fit range for the ratio, we compute the change in the central value when going from the chosen fit to ∆t/a = 10, as marked with an X in Fig. 2.
As a cross-check, we also tested an alternative method 0.000 : 0. 39 12 ; 3 Examples of results for the finite-volume matrix elements | π, pπ|Jµ(0, q)|n, P , Λ, r F V |. The left three panels show the data as a function of tJ − tππ for the three different source-sink separations. The right panels show the fitted values for multiple different fit ranges, which are indicated at the bottom. There, the first set of numbers are the included source-sink separations, and the second set of numbers are the distances from the mid-point that are included for each of these source-sink separations. The blue bands show the chosen fit result, and the half-crosses mark the fits that are used to estimate systematic uncertainties. The values of χ 2 /dof are also given. The quantity denoted as LD is the kinematic factor appearing next to 2iVπγ→ππ/mπ in Eq. (1).
for extracting the matrix elements, in which we did not use ratios, but fitted the three-point functions (57) after dividing out the time dependence and overlap factors. That method gives results consistent with the ratio method. Because the ratio (59) also depends on the energies E P , Λ n , we additionally include a second systematic uncertainty associated with the choice of fit range used in the spectrum analysis of Ref. [28]. The numerical results for all kinematic points are listed in Tables IV and V in Appendix A. There, both systematic uncertainties have been added in quadrature to the statistical uncertainties.

A. Lellouch-Lüscher factors
The mapping between a finite-volume matrix element | π, p π |J µ (0, q)|n, P , Λ, r F V | calculated on the lattice and the corresponding infinite-volume matrix element | π, p π |J µ (0)|s, P , Λ, r IV |, for our normalization of states, is [29,31,32,36] | π, p π |J µ (0)|s, P , Λ, Note that the current in the infinite-volume matrix element is evaluated in position space at x = 0, while the current in the finite-volume matrix element is projected to momentum q. The energy-dependence of the ππ Pwave scattering phase shift δ has to be determined from the Lüscher analysis on the same lattice. We use our Breit-Wigner fits from Ref. [28], as already discussed in Sec. II. The function φ P ,Λ in Eq. (60) appears in the Lüscher quantization condition as where w lm is defined as with the generalized zeta function Z P lm and the Lorentz gamma factor γ.
The quantization conditions for cot φ P ,Λ used are discussed in Sec. VI of Ref [28]; the nonzero factors c P ,Λ lm appearing in elastic P -wave ππ scattering are also listed in Table II. The right-hand side of Eq. (60), known as the Lellouch-Lüscher factor, depends on the ππ system's momentum P , irreducible representation Λ, invariant mass s P ,Λ n , and scattering momentum k P ,Λ n . In Fig. 3 we show the Lellouch-Lüscher factors as a function of invariant mass.
Calculating the derivative ∂φ P ,Λ ∂E in practice means that we must calculate the derivative of w lm (k 2 ): where m 1 , m 2 are the two hadron masses; in the case of ππ scattering m 1 = m 2 = m π . In the rest frame, the derivative of Z lm is again a zeta function: Since this does not hold in moving frames, we compute the derivative numerically. In Fig. 3 we can see that the two different models for the phase shift δ, BW I and BW II, are statistically compatible. Nevertheless, we use both Breit-Wigner models in our analysis to quantitatively assess this.
The fitting systematic uncertainties in E P ,Λ n enter in the Lellouch-Lüscher factors not only via the explicit factor of E P ,Λ n in Eq. (60), but also through the phase-shift parametrization fitted to these energies via the Lüscher quantization condition. In Ref. [28], we estimated the systematic uncertainties in E P ,Λ n by comparing the results of exponential fits with start times t min and t min +a. To correctly propagate these uncertainties to the Breit-Wigner parameters, we then performed the Lüscher analysis and the Breit-Wigner fits for both sets of energies [28]. In the present work, we therefore also repeat the mappings of the πγ → ππ matrix elements (and the subsequent analysis) for both sets of Breit-Wigner parameters.

B. Lorentz decomposition of the infinite-volume matrix elements
The infinite-volume matrix elements s, P , Λ, r|J µ (0, q)|π, p π IV obtained from Eq. (60) still carry the finite-volume irrep indices P , Λ, r. The infinite-volume states s, P , Λ, r| are linear combinations of the states labeled by the continuum polarization index m in Eq. (3). The coefficients of these linear combinations are given by the irrep projection formula Eq. (25). We form the same linear combinations of the polarization vectors on the right-hand side of Eq. (1) to obtain the irrep-projected form-factor decompositions. Taking this into account, we can determine the values of the infinite-volume transition amplitude V πγ→ππ . Most kinematic points only have a single possible Lorentzdecomposition factor, but at certain values of (s, q 2 ) there are two, as shown in Fig. 4. We average over the two resulting values of V πγ→ππ , which reduces the full set of 59 matrix elements to 48 distinct kinematic points (s, q 2 ).

A. Parametrization of the infinite-volume transition amplitude
To allow the calculation of observables, the transition amplitude V πγ→ππ (q 2 , s) determined with lattice QCD at 48 discrete values of q 2 and s needs to be fitted to an analytic parametrization. In Sec. II, we factored out the ρ pole in s according to Watson's theorem, What remains is the transition form factor F (q 2 , s), which should not have any additional poles in s in our region of interest. To obtain a modelindependent parametrization of F (q 2 , s), we perform a two-dimensional Taylor expansion in the variables and after dividing out the lowest expected pole in q 2 : The variable S was chosen to be dimensionless and small near the resonance. The definition of z maps the complex q 2 plane, cut along the real axis for q 2 > t + , to  the interior of the unit circle [79][80][81][82][83][84]. The constant t 0 determines which value of q 2 is mapped to z = 0; we choose t 0 = 0. The constant t + should be set to the lowest branch point. For the QED current, the branch cut starts at (3m π ) 2 and the lowest pole is located at m 2 ω . However, because we neglect the disconnected contributions, we use t + = (2m π ) 2 and m P = m ρ 2 In practice, the series (68) needs to be truncated. We organize these truncations into three different families: F2) Order N in z, combined order K: F3) Order N in z, order M is S: The first two families, F1 and F2, cut the series at the combined z and S order, while the third family F3 separately specifies the orders in z and S. In the limit of large K, N , M , all parametrizations become equal.
In the construction of χ 2 , we take into account the uncertainties in all z and s values by promoting these values to nuisance parameters, like we did (for the s values) in Ref. [28]. The covariance matrix, which we estimate using single-elimination jackknife, is therefore a 3N data × 3N data matrix, where N data = 48 is the number of kinematic points. We added the systematic uncertainties associated with the choices of fit ranges in the matrix element fits and spectrum fits in quadrature to the diagonal elements of the covariance matrix. The uncertainties of the best-fit parameters are obtained from the Hessian of χ 2 at the minimum.

B. The fit results
For each of the different families of parametrizations F1-F3 we investigate several fits while keeping the power of z below 3, and power of S below 4. We find that when the z-expansion goes to order n = 3 or higher, the additional parameters are consistent with zero and 2 Because m 2 P > t + , it is not actually necessary to factor out the pole, but there is no harm in doing so. no longer contribute to the description of the data; similarly, for the S expansion, at order m = 4 the parameters become statistically consistent with zero. We drop all parametrizations yielding fit parameters with uncertainties larger than 100 times their central values. We also remove parametrizations that lead to χ 2 dof > 1.1, which includes those that are of 0-th order in the z-expansion. The list of models that we keep in our analysis, and their corresponding values of χ 2 dof , are given in Table III. In Fig. 5 we present the fitted V πγ→ππ combined with the data points in a three-dimensional plot as a function of √ s and q 2 . ues of √ s are plotted as a function of q 2 in Fig. 7, where the upper panel shows the slices with √ s ≥ m R while the lower panel shows the slices with √ s < m R . We can see that the parametrization describes both the √ s and q 2 dependence of the data well. Qualitatively, we can see two main features in V πγ→ππ : the amplitude is falling off as q 2 decreases, and shows the expected enhancement in √ s attributed to the ρ resonance. The amplitude vanishes at the threshold 2m π , then rises and falls steeply as the resonance region is crossed. This can also be seen in Fig. 8

large
√ s is slower than what would be expected for purely resonant behavior, indicating that the πγ → ππ transition probability remains sizable even when the invariant mass is far above the resonance position. This is also reflected in Figs. 9 and 10, where we plot the function F (q 2 , s) that does not contain the Breit-Wigner factor. The slow falloff of V πγ→ππ as a function of √ s corresponds to growing F . The other parametrizations show the same behavior, confirming a nontrivial s-dependence of the function F (q 2 , s).

VIII. OBSERVABLES
As discussed in Sec. II we consider two main observable quantities, both with a real photon (q 2 = 0): the π + γ → π + π 0 cross section and the ρ radiative decay width. The π + γ → π + π 0 cross section (16) evaluated with our nominal parametrization "BWII F1 K2" of V πγ→ππ (s, q 2 = 0) is shown in Fig. 11. Note that we evaluated Eq. (16) using the heavier-than-physical pion mass of this ensemble, m π ≈ 320 MeV. Because the ρ resonance is narrower than in nature, the peak value of the cross section is higher [44].
To determine the ρ radiative decay width, Γ(ρ → πγ), we must first determine the photocoupling G ρπγ , which requires us to analytically continue the transition amplitude V πγ→ππ to the pole position. The resulting resonant form factor F πγ→ρ (q 2 ), defined in Eq. (13), is presented in Fig. 12. We find that the imaginary part of the resonant form factor is consistent with zero, and the real part slowly rises as a function of q 2 . The resonant form factor at q 2 = 0 is equal to the photocoupling, G ρπγ = F πγ→ρ (0). Our results for G ρπγ , now for all The π + γ → π + π 0 photoproduction cross section as a function of π + π 0 invariant mass, computed with the nominal parametrization "BWII F1 K2" of the amplitude, for our pion mass of mπ ≈ 320 MeV. The inner (darker) shaded region indicates the statistical and systematic uncertainties, and the outer (lighter) shaded region also includes the parametrization uncertainty, estimated as explained in the caption of Fig. 10.
fourteen amplitude parametrizations that gave good fits, are shown in Fig. 13. The real and imaginary parts of the resonant form factor Fπγ→ρ(q 2 ) obtained by analytically continuing the nominal parametrization "BWII F1 K2" of the πγ → ππ amplitude to the ρ resonance pole. The inner (darker) shaded region indicates the statistical and systematic uncertainties, and the outer (lighter) shaded region also includes the parametrization uncertainty, estimated as explained in the caption of Fig. 10  We find that the photocouplings extracted from the different parametrizations are consistent with each other. Nevertheless, we estimate a systematic uncertainty associated with the choice of parametrization as where x i is the photocoupling determined from the ith parametrizations, N = 14 is the number of different parametrizations, and x chosen is the value obtained from the nominal parametrization, "BWII F1 K2". Our final result for the photocoupling is where the first uncertainty includes the statistical uncertainty and the systematic uncertainty from the twopoint and three-point function fits, while the second uncertainty is our estimate (72) of the parametrization dependence. The kinematic factors in Eq. (17) lead to a strong pionmass dependence of the ρ radiative decay width. We can calculate the decay width for the physical pion mass under the assumption that the pion-mass dependence of the photocoupling is negligible. This gives Γ(ρ → πγ) = 168 (13) where we used m ρ = 775 MeV and m π = 140 MeV. For comparison, the experimental value of the ρ ± radiative decay width is 68 (7) keV [61].

IX. CONCLUSIONS
We have presented a (2 + 1)-flavor lattice QCD calculation of the πγ → ππ process, where the ππ system has I = 1 and J P C = 1 −− . The ensemble used has light-quark masses that correspond to a pion mass of approximately 320 MeV, while the strange-quark mass is approximately at its physical value. For the ππ system, we utilized the same moving frames and irreducible representations as in our previous study of ππ scattering [28]. We determined the transition amplitude V πγ→ππ (q 2 , s) with few-percent uncertainty in a broad kinematic region around the ρ pole in invariant mass s and around zero momentum transfer q 2 , using model-independent parametrizations based on a series expansion in the variables z and S, defined in Eqs. (67) and (66). The results obtained from several different truncations of the series are consistent with each other. We observe the expected enhancement of the amplitude associated with the ρ resonance, but find that for large √ s the amplitude falls off slower than expected for purely resonant behavior. In our analysis, we compared two different Breit-Wigner parametrizations of the ππ scattering phase shift (with and without a Blatt-Weisskopf barrier factor). These parametrizations yield consistent results for V πγ→ππ (q 2 , s) in most of the kinematic range, but differ for large √ s.
By analytically continuing V πγ→ππ (q 2 , s) to the ρ pole, we also determined the πγ → ρ resonant form factor and the ρ photocoupling. All truncations of the series used for V πγ→ππ (q 2 , s), and both Breit-Wigner functions, lead to consistent results for the photocoupling, as can be seen in Fig. 13. Our final result for this coupling is |G ρπγ | = 0.0802 (32) (20), which is approximately 1.6 times the value extracted from the measured ρ ± radiative decay width [61] using Eq. (17), |G ρπγ | exp = 0.0508 (26). A significant difference between the lattice result and the experimental value can be expected, given that we did not perform an extrapolation to zero lattice spacing and physical pion mass, and we did not estimate the resulting systematic errors. Most of the past lattice studies of this quantity [77,[86][87][88] were performed in the singlehadron approach, in which the coupling of the ρ to the ππ system is not taken into account. The authors of Refs. [43,44] used the multi-hadron approach at a pion mass of approximately 400 MeV and obtained a value of |G ρπγ | around 0.12, as can be seen in Fig. 12 of [44].
Future calculations at lower pion masses, larger volumes, and additional values of the lattice spacing are needed to extrapolate to the physical point. One aspect that also requires more attention is the residual contamination from higher excited states in the ratios used to determine the matrix elements from the correlation functions. Better control over this contamination can be achieved by using more than three source-sink separations and employing more advanced analysis methods [89].
The lattice methods used here to compute a 1 → 2 transition are also applicable to many other processes of interest in nuclear and high-energy physics. An important example is the rare decay B → K * (→ Kπ) + − [90,91]; new lattice calculations of the B → K * form factors that take into account the strong decay of the K * are needed.

ACKNOWLEDGMENTS
We are grateful to Kostas Orginos for providing the gauge-field ensemble, which was generated using resources provided by XSEDE (supported by National Sci- FIG. 14. As in Fig. 2, for additional kinematic points.
FIG. 19. As in Fig. 2, for additional kinematic points. . The systematic uncertainties from the fits to the ratios (cf. Sec. V) and from the spectrum fits (cf. Ref. [28]) have been added to the statistical uncertainties in quadrature.