Rho resonance, timelike pion form factor, and implications for lattice studies of the hadronic vacuum polarization

We study isospin-1 P -wave ππ scattering in lattice QCD with two flavors of O ð a Þ improved Wilson fermions. For pion masses ranging from m π ¼ 265 MeV to m π ¼ 437 MeV, we determine the energy spectrum in the center-of-mass frame and in three moving frames. We obtain the scattering phase shifts using Lüscher ’ s finite-volume quantization condition. Fitting the dependence of the phase shifts on the scattering momentum to a Breit-Wigner form allows us to determine the corresponding ρ mass m ρ and g ρππ coupling. By combining the scattering phase shifts with the decay matrix element of the vector current, we calculate the timelike pion form factor, F π , and compare the results to the Gounaris-Sakurai representation of the form factor in terms of the resonance parameters. In addition, we fit our data for the form factor to the functional form suggested by the Omn`es representation, which allows for the extraction of the charge radius of the pion. As a further application, we discuss the long-distance behavior of the vector correlator, which is dominated by the two-pion channel. We reconstruct the long-distance part in two ways: one based on the finite-volume energies and matrix elements, and the other based on F π . It is shown that this part can be accurately constrained using the reconstructions, which has important consequences for lattice calculations of the hadronic vacuum polarization contribution to the anomalous magnetic moment.


I. INTRODUCTION
The study of hadronic resonances in terms of the underlying theory of QCD necessitates a nonperturbative treatment. Lattice QCD has emerged as a versatile tool enabling ab initio determinations of many hadronic properties [1]. The ρ meson, which is the simplest QCD resonance and decays almost exclusively into two pions [2], is interesting for several reasons: It serves as a benchmark for the finite-volume formalism pioneered by Lüscher [3][4][5], whose practical implementation poses a number of challenging tasks. Furthermore, the relevant correlation function has a rather favorable noise-to-signal ratio compared to those for other resonances, due to the ρ being the lightest isovector resonance.
Beyond its role as a benchmark, the precision study of the ρ resonance has a number of interesting applications. A good understanding of the ρ → ππ channel is a vital component for any study of more complicated resonances, where the ρ is an intermediate decay channel. Thus, the ρ has been subject to many lattice QCD studies already [6][7][8][9][10][11][12][13][14][15][16][17][18][19]. Second, using the approach suggested by Meyer [20] (which is closely related to work by Lellouch and Lüscher [21]), the pion form factor, F π , can be determined in lattice QCD in the timelike region. For first lattice implementations of this method see [16,22].
An interesting and increasingly relevant application of lattice calculations of F π arises in the context of ab initio determinations of the hadronic vacuum polarization (HVP) contribution to the muon's anomalous magnetic moment, a hvp μ . The latter is accessible via the (spatially summed) vector correlator Gðx 0 Þ [23][24][25], which, at large Euclidean times x 0 is dominated by the two-pion channel. Given sufficiently precise data for F π , one can accurately constrain the long-distance regime of Gðx 0 Þ which helps to significantly reduce both statistical and systematic uncertainties in lattice calculations of a hvp μ [26]. The outline of this paper is as follows: In Sec. II we summarize the methods used for determining the isospin-1 scattering phase shift and the timelike pion form factor from our lattice calculations. Section III presents our results for the scattering phase shift, while Sec. IV contains the results for the timelike pion form factor. Implications for the calculation of the leading order HVP contribution to the muon anomalous magnetic moment a μ are discussed in Sec. V. Finally, Sec. VI summarizes our results. Our analysis supersedes previous preliminary results presented in [27,28].

II. METHODOLOGY
A. Determination of the finite volume energy spectrum To study the ρ resonance, we first need to extract a tower of low-lying energy levels. The strategy we use is to build a matrix of correlation functions using interpolating field operators with the quantum numbers of the ρ meson. The lowest states of the spectrum can be extracted using the variational method [29][30][31][32]. We start by forming a correlator matrix from the correlators formed of interpolating operators O i ðtÞ for the ρ and ππ states in a given frame and then solve a generalized eigenvalue problem (GEVP), CðtÞvðt; t 0 Þ ¼ λðt; t 0 ÞCðt 0 Þvðt; t 0 Þ; ð2:2Þ for this matrix. The nth eigenvalue λ n asymptotically decays exponentially with the energy E n of the nth state.
There are different ways of choosing the parameter t 0 in the GEVP; one of them is to keep t 0 constant (the "fixed-t 0 method") and another way is to use the "window method" [32], which keeps the window width t w ¼ t − t 0 constant. For suitable choices, the latter ensures that the leading excited state contamination to λ n ðt; t 0 Þ from the finite correlator basis comes from ΔE n ¼ E Nþ1 − E n [32], where N is the size of the basis.

ð2:4Þ
The momenta p 1 and p 2 of the single pions add up to the frame momentum P, i.e., p 1 þ p 2 ¼ P ≡ 2π=Ld. The single-pion interpolators are defined by π þ ðq; tÞ ¼ 1 2L 3=2 X x e −iq·x ðūγ 5 dÞðx; tÞ; ð2:5Þ π − ðq; tÞ ¼ 1 2L 3=2 X x e −iq·x ðdγ 5 uÞðx; tÞ: ð2:6Þ In a finite hypercubic volume, the rotational symmetry Oð3Þ of the continuum is reduced to that of a discrete subgroup. The operators are therefore classified by the irreducible representations (irreps) of the respective subgroup. The set of irreps depends on the momentum frame used. In this work, we are using a center-of-mass frame (CMF) as well as moving frames with three different lattice momenta with a maximum momentum of P 2 ¼ 3ð2π=LÞ 2 , i.e., d 2 ¼ 3. Correlation functions are computed for all such moving frames that can be realized on a lattice of spatial size L. Frames that share the same absolute momentum are averaged over.
In the rest frame, continuum operators O J with spin J are subduced [34] into the irreps Λ of the octahedral group via where M are the magnetic quantum numbers of J, S J;M Λ;μ are the subduction coefficients, and μ is the row of the finite volume irrep Λ. The J in O ½J Λ;λ is in brackets because, although it was produced only from operators with spin J, the operator can now have an overlap with all other spins which are contained in Λ [35].
In moving frames, there is a further reduction of symmetry, namely into the subgroup of the octahedral group that keeps P invariant [35], which is referred to as the little group [36]. To subduce continuum operators into the lattice irreps of the moving frame, we need helicity operators where λ is the helicity index and D ðJÞÃ Mλ ðRÞ is a Wigner-D matrix [37] for the transformation R that rotates jpjê z into p [38]. This allows a further subduction into little group irreps Λ, forming a so-called subduced helicity operator We construct multiparticle operators from linear combinations of products of single-particle operators with definite momentum. A general ππ creation operator in an irrep Λ can be written [35] ðππÞ where fp 1;2 g Ã is the group orbit of p 1;2 , i.e., the set of momenta that are equivalent under an allowed lattice rotation. C is a Clebsch-Gordan coefficient which couples the irreps Λ 1 and Λ 2 of the single-pion creation operators π † ðpÞ with the irrep Λ of the ðππÞ † operator. These singlepion irreps are either the A − 1 irrep of the cubic group for p ¼ 0 or the A 2 irrep of the little group of p for p ≠ 0. The coefficients relevant for this work are listed in [35,39].
In the isospin limit G-parity allows only contributions from odd partial waves [40]. Taking these reductions of symmetry into account, the relevant irreps of the ρ → ππ channel, where J P ¼ 1 − and where l ¼ 1 is the dominant contributing partial wave, are listed in Table I.
In addition to the correlator matrix CðtÞ, the calculation of the timelike pion form factor in Sec. IV also requires the matrix elements hJ μ ðx ¼ 0; tÞO † i ð0Þi, for both the local (single-site) current, J l μ ðxÞ ¼ where ψðxÞ ¼ ðu; dÞ T and τ 3 ¼ diagð1; −1Þ. In analogy to the single-meson operators, the spatial components of the current operators J μ ðxÞ are projected into the respective irreps Λ, yielding J Λ ðxÞ. In what follows, the superscript Λ will be omitted in all equations where the irreps are treated the same way. We extract the relevant information on the ground state and the first few excited states from the eigenvectors v n ðtÞ of the corresponding eigenvalues λ n ðtÞ determined via the solution of the GEVP of CðtÞ.
The former are used to define operators X n ðtÞ that project on the state with energy E n : The corresponding two-point function is defined as which is the (approximate) projection of the correlation matrix C ij ðtÞ onto the correlator corresponding to the nth state. We investigated the eigenvectors v n on each time slice and have chosen to use the vectors from the earliest time slice after which the absolute value of their components plateaued. At large times, remnant contributions from other states in D nn ðtÞ are expected to be exponentially suppressed such that only the nth state survives: D nn ðtÞ → jZ n j 2 expð−E n tÞ: ð2:15Þ Z n ¼ hΩjX n jni is an overlap factor with state n of the optimized interpolating operator X n . From an exponential fit to D nn ðtÞ we extract jZ n j for our further analysis. The operators X n are then used to form a two-point function with the current insertions at the sink, which again has a large-time behavior dominated just by one state: hJðtÞX † n ð0Þi → hΩjJðtÞjniZ Ã n e −E n t : ð2:17Þ The timelike pion form factor requires the knowledge of the matrix element hΩjJðtÞjni, which can be extracted either by fitting an exponential function to D nn ðtÞ and hJðtÞX † n ð0Þi or by forming the ratio [16]: We also computed two other ratios with the same asymptotic value proposed by the authors of [16]. Similar to that work, we find RðtÞ produces the most precise plateaus of the three and is not reliant on the fit to Eq. (2.15) for the extraction of Z n . Therefore we fit a constant to jR E n ðtÞj 2 ¼ jhΩjJðtÞjnij 2 to extract the plateau value, which we denote jA n j 2 . B. The distillation method The two-pion operators are nontrivial to compute, due to so-called sink-to-sink quark lines which require all-to-all propagators to be computed. To facilitate this task we are using the "distillation" [41] and stochastic Laplacian Heavyside (LapH) smearing [42] methods.
With distillation [41], a smearing matrix S xy ðtÞ is constructed in the following way: We start with the lattice spatial Laplacian, The definition of the actual smearing matrix is v ðkÞ ðx; tÞv †ðkÞ ðy; tÞ ≡ VðtÞV † ðtÞ: ð2:21Þ One main advantage of this approach is that this smearing matrix can be split and used to project propagators into the subspace spanned by the N ev eigenvectors, a much smaller number than the 3N 3 L color fundamental fields on each time slice which are naively needed to save a propagator.
Particularly on larger lattices, because the total computational cost scales with the cube of the physical volume or higher for fixed smearing, distillation is often treated stochastically [42]. In this approach noise partitioning (also referred to as dilution) in the space spanned by the Laplacian eigenmodes [44,45] is used to reduce the variance of the stochastic estimator. With a suitable dilution scheme, using just one noise per quark line typically produces a statistical uncertainty due to the stochastic estimation of the quark propagation that is of the same size or smaller than the one from the Monte Carlo path integral.
A quark line, i.e., a smeared-to-smeared propagator within stochastic Laplacian-Heavyside smearing, can be computed via where v, w are LapH source or sink vectors ϱ,ρ or φ;φ.
Single-meson correlation functions are a product of two such meson functions, for example, ð2:28Þ which uses the Einstein summation convention for the dilution indices b 1 , b 2 . As our correlation functions can contain two-pion operators at both the source and the sink, we evaluate expressions with products of up to four meson functions. Correlation functions with a vector current at the sink require a propagator that is not smeared at the sink. This can be computed via [16,46] where ϕ is a LapH unsmeared sink vector, and likewiseφ ½b ðρÞ ¼ γ 5 ϕ ½b ðρÞ yields an estimator for SD −1 .

C. Gauge field configurations and distillation schemes
We use three gauge field ensembles with two dynamical mass-degenerate light flavors of nonperturbatively improved Wilson quarks [47,48] generated by the Coordinated Lattice Simulations (CLS) consortium using the DDHMC algorithm and software package [49,50]. The ensembles were generated with β ¼ 5.3 corresponding to a lattice spacing of a ¼ 0.0658ð7Þð7Þ fm, and Table II lists key parameters of these ensembles along with the number of configurations used in our study.
We use different dilution schemes for quark lines connected to the source time slice and for sink-to-sink quark lines. Lines connected to the source time slice use full spin dilution and full time dilution. Full Laplacian eigenvector dilution is used on E5, while interlace-12 eigenvector dilution (LI12 in the notation of [42]) is used on F6 and F7. The perambulators for sink-to-sink (sts) quark lines are calculated with full spin dilution and interlace-8 time dilution (TI8). On E5 sink-to-sink lines use LI8, while LI12 is used on F6 and F7.
For the calculation of Laplacian eigenmodes, the PRIMME package [51] is used with a preconditioner built from Chebyshev polynomials [52]. Our code uses the library QDP++ from USQCD [53] and the deflated SAP þ GCR solver from the openQCD package [54]. For cross-checks of the analysis the package TwoHadronsInBox [55] was used.

III. THE ρ RESONANCE
In this section the determination of the energy spectra, the calculation of the phase shift from the energies, and the  resulting resonance mass m ρ and coupling g ρππ are described.

A. Energy spectra
The pion masses on the three ensembles have been extracted using a cosh-fit ansatz. The fit ranges and results are shown in Table III. We also solved the GEVP in the window method for the eight irreps listed in Table I. The extracted energy levels of two selected irreps on the F6 lattice, together with the effective energies E

B. Lüscher formalism
Lüscher's finite volume method [3,4] is used to map the energy levels of the finite-volume lattice box to the continuum phase shift. 1 For the ρ, we are interested in the l ¼ 1 partial wave. In principle, higher partial waves also contribute to the spectrum. While the effect of the l ¼ 3 and l ¼ 5 partial waves has been studied in [11,16], we expect it to be negligible at our level of statistical precision. With this restriction, the quantization condition reads In this equation, k ¼ ð2π=LÞq are the scattering momenta, δ 1 ðkÞ is the l ¼ 1 infinite volume phase shift, and ϕ d Λ ðqÞ is a kinematical function related to modified zeta functions, which can be computed to arbitrary precision. The centerof-mass energy is given by With the spectrum data from the GEVP we can use this relation to map out the infinite volume phase shift in the energy region 2m π < E < 4m π .
The results from this procedure are shown in Fig. 2 for all three ensembles used. The curve in this plot is a fit to a Breit-Wigner parametrization, ð3:2Þ which is motivated in the resonance region by the effectiverange formula. Note that in the limit of a narrow resonance, the parameters m BW ρ and g BW ρππ become the mass and with the covariance matrix lat;i;k −Ē lat;i ÞðE lat;j;k −Ē lat;j Þ; ð3:6Þ calculated from the n jk jackknife samples of the lattice energies and their central valuesĒ lat;i . By minimizing this χ 2 -function on each jackknife sample, we can obtain fit values for the parameters. One advantage of this approach is that we can use any parametrization suitable for the situation. We can compare our form factor results to the Gounaris-Sakurai (GS) parametrization [57] of the phase shift, which is characterized by the parameters m GS ρ and and show the results in Table V. The two fits produce consistent results, although the Gounaris-Sakurai parametrization yields slightly higher values of χ 2 . Figure 3 shows the world data for the coupling g ρππ from various 2 and 2 þ 1 flavor simulations. There is no   [6,7,11,12,[14][15][16]18,19,22], while the lower pane shows results with dynamical light quarks only [8][9][10]13,17]. Where available, the scale-setting uncertainty provided by the authors has been added in quadrature to obtain the errors on the horizontal axis. The value extracted from the physical ρ-meson width is indicated by the magenta star and the black dashed line. The results from this work are the red open squares in the lower pane. With the exception of [22] which quotes result from a modified Gounaris-Sakurai parametrization, the plotted couplings refer to the Breit-Wigner parametrization.  [58]. The difference in jA l j and jA c j is likely due to cutoff effects, which are studied in [58,59].  (22) significant dependence on the pion mass, and the lattice results are generally close to the physical value.

IV. THE TIMELIKE PION FORM FACTOR
To determine the timelike pion form factor we first need to calculate the matrix elements jA l=c j n ¼ jh0jJ l=c jnij. The subscripts l=c refer to the local and the conserved currents, respectively. Our results for the matrix elements jA l=c j are listed in Table VI. There are sizable differences between jA l j and jA c j, likely due to cutoff effects, which are studied in [58,59]. This is a clear indication that an improved version of the currents (defined, e.g., in [60,61]) would be preferable. These differences are supposed to vanish in the continuum limit, but we cannot check this since we are only considering a single lattice spacing.
We now have all the input to compute the timelike pion form factor [20,62], with the Lorentz-boost γ ¼ E E cm . This equation includes derivatives of the infinite volume phase shift δ 1 ðkÞ obtained in the previous section, and of the modified Lüscher zeta functions ϕ d Λ , which were also used in the phase-shift analysis, and which can be obtained to any desired mathematical precision. We want to compare our form factor results to another study [58], which was performed on the same ensembles and which used correlators with one local and one conserved vector current. To conform with that study, we define the local-conserved version of jAj 2 , jA lc j 2 ≡ jA l jjA c j: ð4:3Þ For another literature comparison [63], we use the locallocal version. Equation (4.1) allows us to directly determine F π ðsÞ from lattice data for discrete values of s, using a parametrization of the phase shift as well as the current matrix elements. To get a continuous description of F π ðsÞ, we can use the Gounaris-Sakurai (GS) parametrization [57], given by the resonance parameters m ρ , Γ ρ , ð4:4Þ with the definitions from Eqs. (3.7)-(3.10). The comparison of our lattice-calculated values for F π and the Gounaris-Sakurai curves is shown in Fig. 4. We want to stress that these curves are not fits to the form factor data. The Gounaris-Sakurai curve seems to describe our data reasonably well, but it would be desirable to have a fit to our form factor data extracted from lattice QCD. One way to realize such a fit is an n-subtracted Omnès representation [64,65]  where P n−1 ðsÞ is a polynomial function of degree n − 1. We parametrize the phase shift δ 1 ðs 0 Þ in this equation using the Breit-Wigner form, Eq. (3.2), and our extracted resonance parameters. For the two-subtracted version, the polynomial is a constant, The grey curve is the GS representation of F π , which only takes the fit parameters of the phase-shift fit m ρ ; g ρππ into account-it is not a fit to the data pictured in these plots. The vertical red bars indicate the 4m π threshold for each lattice.
with the square radius hr 2 π i of the pion. The polynomial for the three-subtracted version reads with the curvature c π V of the pion form factor. The integrand of has a pole at s 0 ¼ s, and in order to solve the integral numerically we need to do a subtraction, Z ∞ The integral OðsÞ can now be computed analytically. We divide the lattice data F π ðsÞ by the function O n ðsÞ and fit the result using the function f fit ðsÞ ¼ expðP n−1 ðsÞsÞ. The results of the fit to the three-subtracted version are shown in Fig. 5. In the two-subtracted version, our data were not very well described by the fit function. While the fits describe the F π data much better than the GS representation of the form factor, they tend to have rather large values of χ 2 =d:o:f:, in particular for ensembles F6 and F7. We also found that the off-diagonal elements of the (normalized) correlation matrix are much larger for F6 and F7 compared to E5. Noting that the 500 gauge configurations of ensemble E5 originate from a much longer Markov chain, it is conceivable that the large values of χ 2 =d:o:f: can be explained by the fact that the level of statistics may not be sufficient to estimate the covariance matrix reliably. Nonetheless, we observe that the curves describe the data quite well, which gives us confidence that the method works as intended. This will be further studied in the future using ensembles with N f ¼ 2 þ 1 flavors of Wilson quarks, which generally contain a much higher number of gauge configurations. Results for the square radius hr 2 π i from this fit are shown in Table VII. The results for the two-and three-subtracted versions differ on the level of 2σ, which is another indication that the two-subtracted version is not enough to describe the data accurately. The square radius was previously determined in [63] by fitting the spacelike pion form factor, computed on the same ensembles we are using in our study. This is a completely different approach and provides a very good cross-check of our fit procedure. Because the authors of [63] employ a local current (as opposed to the local-conserved setup used up to this point), we repeated the analysis using jA l j in Eq. (4.1). The results for the square radius from this analysis are shown and compared to the result from [63] in Table VIII. While both results agree very well for ensembles E5 and F6, we obtain a somewhat smaller square radius on ensemble F7. This observation is discussed further in Sec. V. The comparison of this table with Table VII shows again that discretization effects in our currents are sizable.
As a consistency check our results for the square radius and curvature are plotted as a function of the pion mass, along with the values from phenomenological determinations in Fig. 6. For the square radius we compare to the recent determination from Ref. [67], while for the curvature we use the value from Ref. [68], which also provides an overview of various determinations. Note that the pion mass dependence of our results is in good qualitative agreements with the expectations from Ref. [69]. Lattice results for the curvature have previously been obtained in [70].

V. HADRONIC VACUUM POLARIZATION
Recently, it has been realized that the timelike pion form factor has an important application in the context of lattice calculations of the hadronic contributions to the muon g − 2. The hadronic vacuum polarization contribution, a hvp μ , is accessible in lattice QCD via several integral representations involving the vector correlator [26]. A convenient way to evaluate a hvp μ is based on the so-called timemomentum representation (TMR) [23][24][25], with a known kernel functionKðx 0 ; m μ Þ, the muon mass m μ , and the vector-vector correlator,  [63], where hr 2 π i has been computed from a fit to the spacelike form factor. The difference of our results to the corresponding values in Table VII comes from discretization effects, which are also visible in the matrix elements themselves, shown in Table VI [67] and [68], respectively. The inner error bar on the lattice data denotes the statistical uncertainty, while the outer error bar includes the scale setting uncertainty from the conversion to physical units. VII. Square radius and curvature (in units of 10 −2 ) of the pion obtained from the fit to the n-subtracted Omnès representation of the form factor, using a local-conserved current setup. The Sommer scale r 0 is taken from [66]. We want to stress that the curvature can indeed be calculated from the fit parameter we use, but that the result and particularly the error estimate presented here might not be the physical value. 3.59 (7) 4.98 (7) 6.05 (15) Gðx 0 Þδ kl ¼ − Z d 3 xhJ em k ðxÞJ em l ð0Þi: ð5:2Þ Here, J em μ is the electromagnetic current, ðxÞγ μ sðxÞ þ ÁÁ Á :

ð5:3Þ
A definition of the kernel can be found in [58]. This correlator can be decomposed into an isovector (I ¼ 1) and an isoscalar ( It is also commonly decomposed into connected diagrams from each quark flavor and disconnected diagrams. For comparison with Ref. [58], we will focus on the connected light-quark contribution, G ud ðx 0 Þ ¼ 10 9 G I¼1 ðx 0 Þ. While for small x 0 , this correlator can be precisely computed on the lattice, the signal cannot be traced to arbitrarily large values of x 0 , partly due to the deteriorating signal-to-noise ratio, but also due to the finite time extent of the lattice. Getting a good estimate for the long-distance behavior of Gðx 0 Þ, which is needed to perform the integral to infinity, is one of the main challenges. The general idea is therefore to use the direct lattice data up to some cutoff distance x cut 0 and to determine the part above this distance separately. 2 Reference [58] used a simplistic single-exponential model for the large-time part of G ud ðx 0 Þ, where m ρ was a naive estimate for the rho mass, namely the plateau value of a hρðtÞρ † ð0Þi correlator, and c was determined by fitting G ud ðx 0 Þ. We are improving on this method in our work using two different approaches, one using a reconstruction of the finite-volume correlator and one estimating the infinite-volume correlator. The finite-volume approach uses the information we have about the lowest states in the energy spectrum from the GEVP. We can reconstruct the light-quark correlator with the current matrix elements jA l=c j we already used to compute F π , FIG. 7. The light quark contribution to the integrand of Eq. (5.1) for ensembles E5, F6, F7 (top to bottom). The data points computed in [58] are plotted as black filled circles up to x cut 0 and as open circles above the cut. The bands represent the continuation of the correlator above x cut 0 as discussed in Ref. [58]. Colored symbols denote the data from this work using the reconstructed light-quark correlator G ud n max from Eq.  Fig. 7. Blue triangles represent the integrand reconstructed from the isovector correlator G ud n max of Eq. (5.5) for the corresponding value of n max . Data corresponding to the isovector correlator constructed from the GS parametrization and the three-subtracted Omnès representation of F π are shown as red and blue bands, respectively. The different types of extending the correlator above x cut 0 are used to compute the results for a hvp μ presented in Table IX. The difference between the two-and threesubtracted versions of the Omnès representation integral is too small to be seen on this plot. 2 One can also obtain rigorous upper and lower bounds for the long-time contribution [71,72], which can be improved with knowledge of the spectral decomposition of G ud ðx 0 Þ [59,73].
jA lc j 2 n e −E n x 0 : ð5:5Þ This approach has several advantages: Not only do we get a more precise estimate for the large-x 0 behavior of G ud ðx 0 Þ, but we also have a way to determine the number of states required for a reliable estimate. By computing G ud n max ðx 0 Þ for different values of n max , we can see the estimates converging toward each other. In a region where G ud n ðx 0 Þ agrees within errors with G ud nþ1 ðx 0 Þ, we assume that all energy levels n þ 2 and above will not contribute significantly to G ud ðx 0 Þ. The integrand of Eq. (5.1) for different values of n max can be seen in Fig. 7. We compare it to the data obtained by a direct calculation of the vector-vector correlator on the same ensembles, performed in [58]. Even for values lower than x cut 0 , the contribution obtained only from the first level on E5 saturates the contribution from the lowest two levels. On F6 and F7, the contribution from two levels saturates the contribution obtained from three levels, also at comparably low x 0 . This means that the computation of further levels would not contribute significantly to a hvp μ , and it also shows that a one-exponential tail is not well motivated on F6 and F7. Also, on E5 and F6, our reconstructed data saturate the lattice data from [58] around x cut 0 and are much more precise afterwards. On F7, the correlator data lie significantly above the reconstruction, which might be caused by a correlated fluctuation upward that overestimates the vector-vector correlator. Already starting at about 1 fm, the data from the direct lattice calculation on F7 seem to deviate from the expected behavior, leading to this possible overestimation.
In the infinite-volume approach, the long-time part of the correlator is estimated by evaluating the integral Below the 4m π threshold, 3 ρðsÞ can be parametrized by This approach was also used in [58], where the form factor was estimated using the Gounaris-Sakurai [57] parametrization using the naive rho mass m ρ and an estimation of the width Γ ρ based on its experimental value and an assumed scaling Γ ρ ∝ k 3 ρ =m 2 ρ . 4 In this work, we have several parametrizations of F π and can therefore directly evaluate Eq. (5.6). The result of this is shown in Fig. 8, where we compare the vector-vector correlator G ud F π obtained from the Gounaris-Sakurai and Values for a hvp μ obtained using various methods, in units of 10 −8 . The first line shows the accumulated integral over the lattice data up to x cut 0 . The next four lines show the integral over the long-time tail using the following four methods: (1-exp/GS) is the single-exponential (on E5) or the finite-volume GS parametrization (on F6 and F7), which is a reanalysis of the data from [58]. G ud n max is our extension using the reconstruction of the lightquark correlator using Eq. (5.5). G ud F π reconstructs the vector-vector correlator using Eq. (5.6), where the pion form factor F π is parametrized by the n-subtracted Omnès representation for n ¼ 2 and n ¼ 3. We do not show the results of G ud F π reconstructed using the GS parametrization of F π as it does not describe our data well, as can be seen in Fig. 8. The last three lines show the estimate of a correction for finite-volume effects, based on the difference between G ud F π and G ud n max , or based on the GS parametrizations in Ref. [ [58] 0.03 0.07 3 Because the integrand is exponentially suppressed at high energy, we use this parametrization (and the one of F π ) also above the 4m π threshold. 4 We will not compare the infinite-volume GS results from Ref. [58] with ours. In that work, the GS model was also used for a finite-volume extension of the correlator, and we compare those results with ours in Table IX. from the Omnès representations and for comparison show the estimator with the highest n max from Fig. 7 as well as the Mainz HVP data from [58] again. It is obvious that the Gounaris-Sakurai representation with the parameters taken from our phase-shift analysis is not a good parametrization of our data and leads to an integrand that does not saturate the lattice data. Table IX shows our results for the long-time tail computed using the different methods employed in this work and compares them to the naive estimate obtained in [58] without access to the resonance data from this work. Also shown is the full value for a hvp μ , which is the sum of the contribution from the direct lattice calculation and the different long-time tails. One can see readily from Fig. 7 that on F7, our reconstruction of the vector-vector correlator using Eq. (5.5) does not saturate the data from the direct lattice computation of G ud ðx 0 Þ. When comparing our new values for a hvp μ with the ones from [58] and the chiral extrapolation performed in that work (see Fig. 9), one can see that the value for F7 shifts significantly, but that it comes to an overall better agreement with the chiral extrapolation curve. Because the lattice data on F7 seem to show a large correlated fluctuation already at about 1 fm, and because we are using a transition value of x cut 0 ≈ 1.38 fm, the true value for a hvp μ might be even lower. In any case, the published value for F7 lay prominently above the fit curve of the data points sharing the same lattice spacing, and our analysis brought this data point closer to the curve. A similar issue is observed for the pion radius when comparing our results in Table VIII to the published results in [63].
Although asymptotically finite-volume effects in G ud ðx 0 Þ are suppressed exponentially as e −m π L , in practice these effects can be significant. When the large-x 0 region is dominated by a small number of states, the volume dependence is not in the asymptotic regime [23]. Therefore, it is useful to consider the difference between infinite-volume and finite-volume reconstructions, which provides an estimate of a finite-volume correction. This is also shown in Table IX. On ensemble E5, the correction is statistically significant and roughly þ1%. On F6 and F7, where a finite-volume correction was previously estimated in Ref. [58], our results are consistent with zero and also consistent with the previous estimate.

VI. CONCLUSIONS
We have performed an analysis of I ¼ 1 pion-pion scattering on three different N f ¼ 2 ensembles at fixed lattice spacing. Our spectra have been determined using the variational method for a total of eight different irreps in the center-of-mass frame and three different moving frames with lattice momenta up to d 2 ¼ 3. The spectral information was used in a finite-volume analysis to determine the parameters that characterize the resonance. We have studied the consistency of different parametrizations of the phase shift by comparing the results for the mass m ρ and the coupling g ρππ obtained from fits to either the Breit-Wigner or the Gounaris-Sakurai representation. Our results, shown in Table V, indicate that the parameters are only weakly dependent on the parametrization.
We have used our parametrization of the phase shift together with the matrix elements of the local and point-split vector currents to compute the pion form factor, F π , in the timelike region. While the results for F π agree qualitatively with the Gounaris-Sakurai parametrization based on our ρ meson masses and couplings, they are better described by an Omnès representation obtained from a fit to the F π data, taking the parameters m ρ and g ρππ as input quantities.
The thrice-subtracted version provides a particularly good description and also allows for the determination of the (squared) pion charge radius hr 2 π i. Our results compare well to an independent calculation of the charge radius on the same ensembles, obtained from the slope of the pion form factor in terms of the spacelike momentum transfer Q 2 [63]. When lowering the pion mass, our results for hr 2 π i and the curvature approach the phenomenological values [67,68]. The resulting mass dependence of the squared radius is compatible with the results in Ref. [69].
While the characterization of resonances using lattice techniques is interesting in its own right, the gained FIG. 9. Pion mass dependence of a hvp μ at β ¼ 5.3. Red triangles correspond to the data computed in Ref. [58] with the chiral extrapolation at nonzero lattice spacing represented by the band. The blue circles denote the data points determined from the largex 0 tail of our most precise reconstruction of the correlator, which is the G ud n max correlator using the matrix elements jAj as an input. Points are slightly shifted for clarity. The leftmost red triangle corresponds to ensemble G8 (m π ¼ 185 MeV), which was not considered in this work. Our determination of the tail of the isovector correlator allows for a significantly more precise determination of a hvp μ compared to Ref. [58]. Furthermore, we find that the result for F7 moves closer to the curve when our reconstruction of the large-distance tail is used. information can also be put to good use in different contexts. As another important application we have considered the calculation of the hadronic vacuum polarization contribution to the muon anomalous magnetic moment, a hvp μ . The precision of lattice calculations of a hvp μ is typically limited by the long-distance tail of the vector correlator.
By means of a direct comparison with an earlier study [58], we have shown that the precision in a hvp μ can be substantially increased by describing the long-distance tail of the TMR integrand [see Eq. (5.1)] using the spectral information on the first few states in the isovector channel. Alternatively, the tail of the integrand can be much more accurately constrained via the representation of the vector correlator in terms of the pion form factor. These techniques have, in the meantime, been employed in a recent calculation of a hvp μ on CLS gauge ensembles with N f ¼ 2 þ 1 flavors of dynamical quarks [59]. Going beyond that work, we have used the difference between infinite-volume and finite-volume reconstructions to estimate finite-volume effects; the results are consistent with previous estimates using the Gounaris-Sakurai model.