The $I=1$ pion-pion scattering amplitude and timelike pion form factor from $N_{\rm f} = 2+1$ lattice QCD

The elastic $I=1$ $p$-wave $\pi\pi$ scattering amplitude is calculated together with the isovector timelike pion form factor using lattice QCD with $N_{\rm f}=2+1$ dynamical quark flavors. Wilson clover ensembles generated by the Coordinated Lattice Simulations (CLS) initiative are employed at four lattice spacings down to $a = 0.05\,\mathrm{fm}$, several pion masses down to $m_{\pi} = 200\,\mathrm{MeV}$, and spatial volumes of extent $L = 3.1-5.5\,\mathrm{fm}$. The set of measurements on these ensembles, which is publicly available, enables an investigation of systematic errors due to the finite lattice spacing and spatial volume. The $\pi\pi$ scattering amplitude is fit on each ensemble by a Breit-Wigner resonance lineshape, while the form factor is described better by a thrice-subtracted dispersion relation than the Gounaris-Sakurai parametrization.

The improvement in these calculations suggests that the quark-mass dependence of such amplitudes may be investigated quantitatively, providing valuable input to effective theories of low-lying hadron resonances as well as numbers at the physical point relevant for experiment. In order to obtain reliable results however, various systematic errors must be controlled. These include effects due to the finite lattice spacing and spatial volume inherent in lattice QCD simulations, as well as systematics in the calculation of finite-volume two-hadron energies and matrix elements from which the amplitudes are determined.
The relation discussed above between finite-volume energy shifts and real-time two-to-two scattering amplitudes can be written as [35] det[K −1 (E cm ) − B ( ,d) where E cm is the finite-volume two-hadron energy in the center-of-mass frame, K is proportional to the infinite-volume K-matrix, and B is a known matrix encoding the effect of the finite volume. This determinant is block diagonalized so that Eq. (1.1) describes finite-volume energies in a single irreducible representation (irrep) of the little group for a particular class of total momenta d = (L/2π)P tot , where d is a vector of integers. The determinant is taken over total angular momentum (J ), total spin (S), all coupled two-hadron scattering channels, and an index enumerating the possibly multiple occurrences of a partial wave in irrep . The determinant condition in Eq. (1.1) holds up to corrections which are exponentially suppressed in the spatial extent L. However, unlike finite volume corrections to single-hadron observables, the fall-off of these residual exponential finite volume effects may in principle be different than m π . The 'rule of thumb' m π L 4 which is usually applied in single-hadron calculations to ensure that finite-volume effects are at the percent level must be re-investigated in the context of scattering amplitudes. There exists therefore a hierarchy of finite volume effects: those described by Eq. (1.1) are polynomial in L −1 (and constitute the 'signal') while Eq. (1.1) holds only up to (unwanted) terms exponential in L.
As a benchmark amplitude suitable for an investigation of these systematic effects, we consider here the elastic I = 1 p-wave pion-pion scattering amplitude relevant for the ρ(770). In order to extrapolate to the physical point and continuum, a range of pion masses m π = 200-280 MeV and lattice spacings a = 0.050-0.086 fm are employed. Such extrapolations have been performed recently using chiral perturbation theory and its extensions [63][64][65][66][67][68][69][70][71][72][73][74] and treat all scattering data simultaneously. As these extrapolations are somewhat involved and have not yet treated cutoff effects, we present the determination of the amplitudes only in this work and leave extrapolation to the physical quark masses and continuum for the future. Nonetheless, in our results the residual finite volume and cutoff effects are evidently small compared to the statistical errors.
We calculate scattering amplitudes over an energy range E cm ∈ [2m π , E max ]. Since Eq. (1.1) applies below n > 2 hadron thresholds, E max is reduced as m π is lowered. Because of the chiral trajectory employed in this work (which is discussed in Sec. 2.1), for the heaviest quark masses the lowest inelastic threshold is K K, while the lowest n > 2 hadron threshold is 4π . Although levels in the range 2m K < E cm < 4m π could be treated using Eq. (1.1) with a coupled-channel K-matrix, we nonetheless impose a restriction to elastic scattering, namely E max = min(4m π , 2m K ).
In addition to this ππ scattering amplitude, we also calculate the I = 1 timelike pion form factor, which encodes the coupling of an external (timelike) photon to two pions in an isovector configuration. Phenomenologically, it can be extracted from e + e − → hadrons and hadronic τ -decays [75] and is of particular relevance for the hadronic vacuum polarization (HVP), a leading source of theoretical uncertainty in the anomalous magnetic moment of the muon (g − 2) μ [76,77]. Using the optical theorem, the imaginary part of the HVP can be related to where σ (e + e − → hadrons) is the total cross section, α em the electromagnetic coupling, and s = E 2 cm the usual Mandelstam variable. In the elastic region R had is given by the two-pion contribution which contains the timelike pion form factor F π (s). The phase of this form factor is fixed by Watson's theorem, so we are interested in the amplitude only here. Furthermore, we work in the isospin limit, where electromagnetic interactions are ignored and m u = m d . Because of this, the elastic region persists up to either s = 4m 2 K or s = 16m 2 π . Although a precise determination of this form factor is phenomenologically desirable, there exists only the pioneering determination of Ref. [10] which employs a single lattice spacing, heavier quark masses, and a (single) smaller physical volume than this work. It is therefore imperative to also investigate lattice spacing and finite volume effects for this quantity, the former of which may be affected by renormalization and O(a)-improvement of the electromagnetic current.
In addition to its phenomenological impact, the timelike pion form factor is an important stepping stone toward more complicated resonance photoproduction amplitudes. Such amplitudes are relevant for ongoing and future experiments which photoproduce resonances. An additional step in this direction is the πγ → ππ amplitude studied using lattice QCD in Refs. [24,25]. However, the timelike pion form factor calculated here does not require disconnected flavor-singlet Wick contractions, which are ignored in Refs. [24,25].
Preliminary work toward the results reported here is found in Ref. [78]. The remainder of this paper is organized as follows. For completeness we review the gauge field ensembles, methods for calculating finite-volume two-pion energies and matrix elements, and their relation to infinite-volume scattering amplitudes in Sec. 2. Results are given in Sec. 3 and conclusions in Sec. 4.

Lattice QCD methods
The subset of CLS ensembles used in this work is discussed in Sec. 2.1 and application of the stochastic LapH method for all-to-all quark propagation in Sec. 2.2. The analysis strategy used to extract the required finite volume energies and matrix elements from temporal correlation functions is contained in Sec. 2.3, while the relation between finite-volume quantities and infinite-volume scattering amplitudes is given in Sec. 2.4.

Gauge field ensembles
The ensembles of gauge field configurations employed here are from the Coordinated Lattice Simulations (CLS) initiative and are presented in Refs. [79,80]. They employ the treelevel improved Lüscher-Weisz gauge action [81] and non-perturbatively O(a)-improved Wilson fermions [82]. Open boundary conditions [83] are implemented in the temporal direction. Although these boundary conditions were adopted to reduce autocorrelation times of the global topological charge, they also influence finite-temporal-extent effects in temporal correlation functions.
Contributions to two-hadron correlation functions where the hadrons propagate in opposite temporal directions, which for identical particles and zero total momentum yield a constant in time [84][85][86][87], are present with periodic boundary conditions but absent in this setup. Therefore, for large temporal extent T and if both interpolators are far from the boundaries, all two-point correlation functions with open temporal boundary conditions have the form the correlator in the T → ∞ limit, E 0 the lightest state with vacuum quantum numbers and t bnd = min(t 0 , T − t f ) the minimal distance from an interpolator to the temporal boundaries. Since presumably E 0 ≈ 2m π , if m π t bnd 2 then the exponential corrections in Eq. (2.1) are parametrically similar to exponentially suppressed finite-volume effects in singlehadron energies. While correlation functions are affected by temporal boundary conditions, the transfer matrix (and therefore also the spectrum) is unaffected. Although having an interpolating operator near a temporal boundary does not change its quantum numbers, we are after excited states and employ generalized eigenvalue methods requiring hermitian correlation matrices. The source and sink interpolators therefore must not be significantly affected by the temporal boundaries in order to maintain hermiticity. To this end, we always choose a minimum distance to the boundary (t bnd ) of at least m π t bnd 2. As in Ref. [6], our resulting insensitivity to finite-T effects can be demonstrated using the single-pion correlation function. Fitting this correlation function to a single exponential ignores contributions from the temporal boundaries. Fits of this type on a single ensemble are shown in Fig. 1, where the fitted energy is shown to be insensitive to the source interpolator position t 0 . Insensitivity to t 0 in our most precisely determined correlation function suggests that temporal boundary effects may be neglected in subsequent fits.
Another measure of finite-T effects is the ratio which under the asymptotic assumptions of Eq. (2.1) receives corrections to unity of O(e −E 0 t bnd ). R(5a) is also shown in Fig. 1 for various (t 0 , t 0 ) pairs on all ensembles with multiple source times. This ratio shows significant deviations from unity for m π t bnd 3, despite no observable difference in the fitted energies. While the deviation of R t 0 ,t 0 (t) from unity in the single-pion correlator suggests that averaging over source times may affect the hermiticity of correlation matrices, such deviations are not visible in two-pion correlation functions. For the D200 ensemble the ratio with (t 0 /a, t 0 /a) = (32, 52), shown for the pion in Fig. 1 left panel as the left-most point, is consistent with unity for the single-ρ meson correlator at rest and the two-pion correlator (each with a single unit of momentum) in the same channel. While we omit a complete discussion of algorithmic details used in configuration generation, some aspects are relevant for the analysis of correlation functions measured on these ensembles. The CLS ensembles employ twisted-mass reweighting [88] for the degenerate light quark doublet and use the RHMC algorithm [89] for the strange quark determinant. One reweighting factor (W 0 ) is used to change the light quark action to clover Wilson, while another (W 1 ) corrects for the RHMC approximation. Efficient evaluation of these reweighting factors is discussed in Ref. [79]. Measurements of primary observables must be multiplied by the corresponding re-weighting factors on each configuration according to Table 1 Parameters of the CLS ensembles used in this work. The timelike pion form factor is not determined on the coarsest lattice spacing. After the ensemble ID in the first column, we list the gauge coupling, lattice spacing and dimensions, pseudoscalar meson masses, the separation between correlation function measurements in molecular dynamics units (MDU), the number of such measurements, and the minimum distance from an interpolator to a temporal boundary. 3) must also be taken into account in any resampling procedure used to estimate statistical errors or covariances. If both the twisted mass parameter and the range and degree of the rational approximation are chosen appropriately, these reweighting factors are typically close to unity. However, we observe anomalously large fluctuations in the reweighted zero-momentum single-pion correlator for a single source time on each of the C101 and D101 ensembles. These ensembles have the lightest quark mass at the coarsest lattice spacing and large fluctuations may indicate an inefficient choice of the simulated action. Data on these two source times are removed from the final analysis and are not included in Table 1.
There are several possibilities for the trajectory of m s as m l = m u = m d is lowered toward its physical value. Quark masses on the CLS ensembles employed here are tuned to lie on a chiral trajectory with tr M q = const., where M q is the bare quark mass matrix M q = diag(m u , m d , m s ), in order to reduce the quark mass dependence of certain renormalized quantities. An additional chiral trajectory in which m s = const. is presented in Ref. [80]. While it is interesting to investigate the quark-mass dependence of scattering amplitudes on both chiral trajectories, we present results on the tr M q = const. trajectory only.
As discussed in Ref. [90], fixing the trace of the bare mass matrix is not equivalent to fixing the sum of renormalized quark masses. There a Taylor expansion is employed to slightly shift the quark masses in order to satisfy φ 4 = 8t 0 (m 2 K + 1 2 m 2 π ) = const., which is not performed here. At the coarsest lattice spacing, imposing tr M q = const. results in deviations of less than 5% of tr M R /(tr M R ) symm (where the trace in the denominator is evaluated at the symmetric point m u = m d = m s ) from unity at the lightest pion masses considered in Ref. [90]. This small deviation from our desired chiral trajectory presumably has little effect on the observables considered here.
Properties of the CLS ensembles used in this work are given in Table 1, which also contains τ meas , the separation in molecular dynamics units (MDU) between our measurements of hadronic correlation functions, and t bnd , the minimum distance from an interpolator to a temporal boundary. The timelike pion form factor is not determined on the coarsest lattice spacing. A more precise scale determination can be found in Tab. 3 of Ref. [90], while pseudoscalar meson masses and decay constants can be found in Tab. 2 of that work.
As discussed above, open temporal boundary conditions are employed to decrease the integrated autocorrelation time of the global topological charge. However, there is still a significant amount of autocorrelation present in some observables on the CLS ensembles. A method to estimate statistical errors in the presence of large autocorrelations is outlined in Refs. [91,92]. This involves propagating the errors linearly, a method which may not be suitable for our purposes given the non-linear nature of the B-matrix elements given in Eq. (1.1) and discussed further in Sec. 2.3. Therefore, we simply 'bin' our correlator measurements and employ the bootstrap procedure with N B = 800 bootstrap samples. Although no statistically significant autocorrelations are observed in any of our correlation functions, the largest integrated autocorrelation times (τ int ) measured on these ensembles span the range τ int ≈ 30-150 for β = 3.4-3.7, respectively [79].

Correlation function construction
Since we employ two-pion interpolators in which each pion is projected to definite momentum, quark propagators between all space-time points are required. We employ the stochastic LapH method to estimate such all-to-all propagators and efficiently construct correlation functions [2]. Based on Ref. [1], this method endeavors to make all-to-all propagators tractable by considering quark propagation between a low-dimensional subspace defined by the lowest N ev modes of the three-dimensional gauge-covariant Laplace operator, hereafter referred to as the 'LapH subspace'. This projection is a form of quark smearing, with an approximately Gaussian spatial profile and width controlled by the N ev -th eigenvalue. In order to maintain a constant width, N ev must be scaled proportionally to the spatial volume.
This smearing procedure enables more efficient stochastic estimation schemes by employing noisy combinations of Laplacian eigenvectors. It was determined in Ref. [2] that (at least for the range of spatial volumes considered there) with a moderate level of dilution the quality of the stochastic estimator remains constant as the volume is increased while maintaining a fixed number of dilution projectors. This work further demonstrates that the quality of the stochastic LapH estimator does not degrade for even larger volumes. Without significantly increasing the number of dilution projectors, we obtain precise results for scattering amplitudes with stochastic LapH on spatial volumes up to V = (5.5 fm) 3 .
In the stochastic LapH framework, N R stochastic sources {ρ r } are introduced in time, spin, and Laplacian eigenvector indices. These sources are diluted by specifying N dil complete orthogonal dilution projectors {P b } so that an unbiased estimator of the smeared-smeared all-to-all quark propagator is furnished by where rb = V s P b ρ r is the smeared stochastic source, ϕ rb = S Q rb the smeared sink, Q the quark propagator, and S = V s V † s the smearing operator which projects onto the LapH subspace. To date, only schemes where dilution in each of these indices is done independently have been employed. A common strategy is to interlace n dilution projectors uniformly (denoted In) in the index in question. The 'full' dilution limit (denoted 'F') is recovered if n is equal to the total dimension of the index. Full specification of a dilution scheme therefore specifies a prescription in each of time, spin, and Laplacian eigenvector space. For example (TF, SF, LI8) refers to full dilution in time and spin, and eight dilution projectors interlaced uniformly among the Laplacian eigenvectors.
As discussed in Ref. [2], it is typically beneficial to employ different dilution schemes for 'fixed' quark propagators (denoted by the subscript 'F'), where x 0 = y 0 and 'relative' quark propagators (denoted 'R') where x 0 = y 0 . We therefore employ either full or interlace dilution in Table 2 Parameters of the stochastic LapH implementation used in this work. (ρ, n ρ ) are the stout link smearing parameters, N ev the number of Laplacian eigenvectors, N R the number of independent noise sources, N t 0 the number of source times for fixed quark lines, and N D the total number of light quark Dirac matrix inversions per gauge configuration. Notation for the dilution scheme is explained in the text.
time, full dilution in spin, and interlace dilution in eigenvector space. The dilution scheme and other parameters of the stochastic LapH algorithm employed here are given in Table 2. Table 2 also contains information on the LapH subspace and thus the smearing operator S applied to quark fields in our interpolating operators. Before calculating eigenvectors, the gauge link field entering the covariant 3-D Laplace operator is stout smeared [93]. The stout smearing parameters (ρ, n ρ ) together with the number of retained eigenvectors N ev therefore define our smearing scheme. We maintain an approximately constant physical link-smearing radius (r link /a) 2 = ρn ρ by tuning n ρ appropriately. The quark smearing procedure is defined by retaining all eigenvectors with eigenvalue λ (aσ s ) 2 , where σ s = 1 GeV. As the physical volume (V ) is increased the number of eigenvectors must be scaled as N ev ∝ V . The stout smearing parameters and N ev are given in Table 2.
We employ interpolating operators with light quarks only, since we calculate elastic pion-pion scattering amplitudes. The number of required light quark Dirac matrix inversions per configuration, denoted N D , is also given in Table 2. Our treatment of all-to-all propagators enables us to efficiently evaluate all required Wick contractions involving two-pion and single-ρ interpolators, which are enumerated in Ref. [2]. An unbiased estimator results only if each quark line in a diagram employs independent stochastic sources. As discussed in Ref. [94], in each diagram we typically average over some number of multiple noise 'orderings', i.e. different permutations of the N R available quark lines.
The correlation functions used in pion-pion scattering require smeared quark fields only. However, correlation functions for the timelike pion form factor contain the unsmeared vector current operator. These current correlation functions are easily constructed in the stochastic LapH framework although they employ quark fields which are not projected onto the LapH subspace.
As done in Ref. [95], by exploiting γ 5 -hermiticity it can be ensured that quark fields in the vector current bilinear are always unsmeared sinks φ rb = Q rb . This motivates the construction of 'current sinks' defined as where t = x 0 = y 0 and (d, ) denotes projection onto an irreducible representation (irrep) of the little group of total momentum d. J (d, ) (t) has two noise/dilution indices and can therefore be employed in the correlation function construction procedure of Ref. [2] exactly as a smeared ρ-meson interpolator.
As suggested by Eq. (2.4), in order to save disk space the quark sinks are typically projected onto the LapH subspace before they are written to disk. However, the current functions of Eq. (2.5) must be constructed from unprojected sinks. For fixed quark lines the {φ rb } must be kept in memory until calculation of J (d, ) (t) is complete. After construction of the current functions, the quark sinks are smeared and written to disk in the usual way.
The calculation of the Laplacian eigenvectors is performed using a variant of the thick restarted Lanczos method [96], which entails global re-orthogonalizations of the Krylov subspace. These re-orthogonalizations scale poorly with N ev , so that as L is increased calculation of the Laplacian eigenvectors will eventually dominate the computational cost. However, for the L 5.5 fm volumes considered here the Dirac matrix inversions are still most computationally intensive.
We perform these Dirac matrix inversions using the efficient DFL_SAP_GCR solver in the openQCD software suite. 2 In summation, our workflow consists of three main tasks: (1) Dirac matrix inversion, (2) hadron source/sink construction, and (3) formation of correlation functions. Due to their large storage footprint, the Laplacian eigenvectors are computed first in task 1 and not saved to disk. They are then recomputed during task 2, which is implemented entirely in 3-D. Task 3 then no longer requires any lattice-wide objects and is simply tensor contraction of noise-dilution indices.
These different tasks are typically performed on different computer architectures, but a rough breakdown of the relative cost is 70-80% for the Dirac matrix inversions, 20-26% for task 2, and 1-5% for task 3. In total, 1-3% of the total for these three tasks is spent on calculating Laplacian eigenvectors.

Finite-volume energies and matrix elements
We consider all elastic energy levels in isovector irreps where the J P G = 1 −+ partial wave is the leading contribution up to total momentum d 2 ≤ 4, which are tabulated in Table 3. To calculate the energies, we follow the procedure of Refs. [6,16], which is outlined below.
Outside the resonance region, interacting finite-volume two-pion energies are close to their non-interacting values, while for levels with E cm near m ρ these gaps are larger. To exploit the small differences outside the resonance region and to treat all energies in a unified manner, we employ the ratio fits described in Ref. [6]. Using this method, we construct ratios where (p 1 , p 2 ) are momenta of the constituent pions in the nearest non-interacting level and C π (p 2 , t) is a single-pion correlation function with momentum p 2 . The vector v n (t 0 , t d ) is a generalized eigenvector of the correlator matrix C(t) solving the generalized eigenvalue problem (GEVP) C(t d )v n = λ n C(t 0 )v n [97,98]. The {v n } are used to define the diagonal correlators Ĉ n (t) between operators with optimal overlap onto the nth level, and are determined for a single (t 0 , t d ) only. The fitted energies vary little as these diagonalization times, as well as the operator basis, are varied. Table 3 Finite-volume irreps (second column) of the little group for various classes of total momenta P tot = (2π/L)d (first column) employed here. The superscripts on the partial waves ( ) contributing to that irrep denote the number of multiple occurrences, while the '+' indicates positive G-parity. The difference E n between an energy and its closest non-interacting ππ counterpart is extracted directly using single-exponential fits to the ratio in Eq. (2.6). Alternatively, the interacting energy may be obtained from single-or two-exponential fits to Ĉ n directly. All of these correlated-χ 2 fits are performed over some time range [t min , t max ], the variation of which should not affect the fitted energies for asymptotically large t . Energies obtained from ratio, single-, and two-exponential fits all typically depend little on t max , while ratio fits typically exhibit a reduced dependence on t min as well. However, the excited state contamination in ratio fits may be non-monotonic leading to 'bumps' in t min plots. As an important consistency check, we check agreement of the energies obtained from these three types of fits, different (t 0 , t d ) combinations, and GEVP operator sets.
Our fit ranges are chosen conservatively so that the systematic errors discussed above due to the GEVP and fit ranges are smaller than the statistical ones. This is demonstrated using extensive comparisons of t min -plots for different fit types, (t 0 , t d ) choices, and GEVP operator bases, similar to Refs. [6,16]. Bootstrap resamples of all reweighted correlation functions are publicly available in HDF5 format, 3 as is a python Juypter notebook 4 which performs the entire analysis chain. This tool not only provides an interface to view systematics related to our choices of fitting procedure, fitting ranges, and GEVP, but also enables direct access to all results at each step. Generally, in physical units we take (t 0 , t d ) ≈ (0.5, 0.9)fm, t min = 0.7-1.3 fm, and t max = 2-2.6 fm.
In addition to determining the energies, on the three finest lattice spacings we calculate matrix elements of the electromagnetic current where the ellipsis denotes contributions from heavier quarks. For the vacuum-to-ππ matrix elements considered in this work, we require insertions of the isovector component and a dimensionfive counterterm required to implement O(a)-improvement V a μ =ψγ μ τ a 2 ψ,∂ ν T a μν = i∂ νψ σ μν τ a 2 ψ, (2.8) where ψ = (u, d) T , τ a the usual Pauli matrices in isospin space, σ μν = i 2 [γ μ , γ ν ], and ∂ μ the symmetrized lattice derivative. The isovector index a is taken to be maximal and henceforth omitted.
The determination of the timelike pion form-factor then employs linear combinations where the coefficients b ( ,d) μ project the current onto (a row of) irrep and spatial momentum d. The vector current bilinear appearing in Eq. (2.9) has been renormalized and O(a)-improved non-perturbatively according to where the renormalization and improvement coefficients Z V , b V , b V , and c V are functions of the gauge coupling only in this mass-independent scheme. We take Z V = Z V 1 + ab V m l + ab V tr M q and c V from the non-perturbative determination of Ref. [99]. An alternative determination of the non-singlet current renormalization constants for this lattice discretization is found in Ref. [100]. Another preliminary non-perturbative determination of Z V can be found in Ref. [101], while non-perturbative determinations of b V and b V are performed in Refs. [102,103].
Operationally, we calculate current correlation functions using Eq. (2.5) for both the dimension four and five operators in Eq. (2.8) projected onto definite momentum d and irrep . These current correlation functions, which are vectors in the GEVP index, are given as where J denotes either of the operators in Eq. (2.8) projected according to Eq. (2.9) and Ô ( ,d) is an interpolator for irrep ( , d).
To extract the finite volume matrix elements 0|V ( ,d) | dn , we calculate the current correlation functions (defined in Eq. (2.11)) containing each the two operators in Eq. (2.8). These are used to form optimized current correlation functions using the GEVP eigenvectorŝ where the inner product is taken over the GEVP index. Using these optimized current correlators, we then construct three ratios which plateau to the desired matrix elements asymptotically for large t (up to GEVP systematics) where A n and E n are determined previously from the ratio fits to Ĉ n (t). The final matrix elements are then obtained from a plateau average of these ratios over a range [t min , t max ].
Each of these ratios possesses different excited state contamination. In analogy with the determination of the energies discussed above, consistency of the matrix elements using different ratios, (t 0 , t d ), and GEVP bases provides a stringent check in their determination. GEVP corrections to A n have a different form than those of the energies [104], and our choices of fit ranges are optimized with the energies in mind. For this reason we take R (1) n (t) in Eq. (2.13) as the best estimate of the matrix elements. Nonetheless, all three ratios are typically consistent. Data illustrating the t min -dependence of these ratios and the comparisons mentioned above may also be found in the HDF5 files and Jupyter notebook. After determining the matrix elements for each current operator, we combine them to form the renormalized combination in Eq. (2.10).

Amplitudes from finite volume energies and matrix elements
First we determine the elastic I = 1 p-wave pion-pion scattering amplitude using the finitevolume energies discussed in Sec. 2.3. This amplitude is obtained from the determinant condition introduced for general two-to-two scattering in Eq. (1.1). However, considerable simplification occurs for elastic scattering between spinless identical particles. Here the K-matrix is diagonal in orbital angular momentum with trivial structure in the irrep occurrence index n occ . The box matrix B, which encodes the effect of the finite periodic spatial volume, mixes different orbital angular momenta and is dense in n occ .
All the irreps used in this work are given in Table 3 together with the pattern of partial wave mixing induced by the infinite-dimensional B-matrix. Explicit expressions for all B-matrix elements up to ≤ 6 are given in Ref. [35]. Several off-diagonal B-matrix elements vanish for identical particles, preventing partial wave mixing between even and odd . If contributions from ≥ 3 partial waves are neglected, this simplification provides a one-to-one correspondence between energies in the irreps of Table 3 and where q cm is the center of mass momentum and δ 1 the I = 1 ππ phase shift. This approximation is justified by the near-threshold suppression of higher partial waves and to test it we perform global fits including f -wave contributions. Refs. [9,35] perform similar fits and find such contributions negligible.
We turn now to the determination of the timelike pion form factor |F π (E cm )| using the finite-volume matrix elements calculated according to Sec. 2.4. The relations employed here for zero-to-two matrix elements are given in Refs. [10,43,46] which are based on the seminal work of Ref. [47].
We first define the angle φ where B is from Eq. (1.1). This pseudophase is used together with the physical phase shift to relate the finiteand infinite-volume matrix elements where q cm is the magnitude of the center-of-mass (three) momentum, u 2 = L 2 q 2 cm /(2π) 2 , and  (d, ) |d n but also the derivative of δ 1 . This derivative is obtained from a parametrization of the phase shift points described above and covariances between all data are treated explicitly using the bootstrap procedure. Parametrization of δ 1 (E cm ) and |F π (E cm )| is discussed in Sec. 3.

Results
We first present results for the elastic ππ scattering amplitude. As discussed in Sec. 2.4, if ≥ 3 contributions to Eq. (1.1) are neglected there is a one-to-one correspondence between finite-volume energies and K 11 (E cm ) defined in Eq. (2.14). This energy dependence is parametrized by a Breit-Wigner shapẽ involving two free parameters g 2 ρππ and m 2 ρ /m 2 π . A correlated-χ 2 fit of all points is performed according to the 'determinant residual' method of Ref. [35] with μ = 10, although without ≥ 3 contributions the K −1 -and B-matrices are one-dimensional so that the determinant is trivial. Results for these fit parameters, which are both constrained to be positive, as well as the χ 2 per degree of freedom, χ 2 = χ 2 /N d.o.f , are given in Table 4 for each of the ensembles employed here.
The influence of ≥ 3 partial waves is assessed by enlarging the determinant condition of Eq. (1.1) to include the = 3 contributions noted in Table 3. For this fit the f -wave is parametrized by an unconstrained constant K −1 33 (E cm ) = −(m 7 π a 3 ) −1 yielding results which are also shown in Table 4. We see therefore that there is little dependence on including ≥ 3. Interested readers may perform further fits using App. A, where energies and phase shift points for all ensembles (neglecting ≥ 3) are tabulated and plotted, or the Jupyter notebook described in Sec. 2.3 where bootstrap samples of all energy levels are available. However, it is worth comparing some of the ensembles in this data set here. An explicit check of finite volume effects using the C101 and D101 ensembles, which have the same parameters but different volumes, is shown in Fig. 2. That figure also shows a comparison between the N401 and N200 ensembles, which have (approximately) the same quark masses but different lattice spacing. It is thus evident on these ensembles that both effects are not visible within our statistical errors. Finally, all results for m ρ and g ρππ in shown in Fig. 6, where they are converted to physical units using the scale determined in Ref. [90].
We turn now to results for the I = 1 timelike pion form factor, which are determined according to Sec. 2.3. Results for the form factor from Eq. (2.15) are also tabulated in App. A. As discussed in Sec. 2.3, we employ a non-perturbative determination [99] of c V multiplying the dimension-five counter term in Eq. (2.10). Apart from the one from Ref. [99], there is the preliminary determination of Ref. [101] which obtains values larger in magnitude using a different improvement condition. If c V is non-negligible, the relative magnitude of the leading order matrix elements to this counterterm is of interest. Their ratio is shown in Fig. 3. Given its 5-15% size, a larger c V could indicate larger cutoff effects in the form factor than we observe using c V (g 0 ) from Ref. [99], which is at the few-percent level. We now turn to parametrization of the form factor. Ref. [10] employs the Gounaris-Sakurai parametrization [105] F GS π ( where the notation is from Ref. [106] and q ρ is the center-of-mass momentum at the resonance energy. This parametrization depends only on m ρ and g ρππ , and therefore describes the form factor with no additional free parameters. Additional parametrizations are suggested by unitarity constraints. In the elastic approximation, the form factor satisfies the n-subtracted dispersion relation This dispersion relation has the Omnès-Muskhelishvili solution [107,108] of  of F π . Given the Breit-Wigner parametrization for the phase shift, the twice-subtracted dispersion relation (n = 2) has a single additional parameter (p 1 ) and the thrice-subtracted (n = 3) two additional parameters, p 1 and p 2 . These parameters appear in ln Q n (s), while n [δ 1 ](s) depends only on the Breit-Wigner parametrization of δ 1 . For these fits, we isolate ln Q n (s) by constructing ln(|F π |/ n [δ 1 ]) and fit it to the appropriate polynomial. An example of a thrice-subtracted fit is shown in Fig. 3, while the twice-and thrice subtracted fits are compared in Fig. 5 for the D200 ensemble. Finally, the Gounaris-Sakurai parametrization is compared to the thrice-subtracted fit on the J303 and D200 ensembles in Fig. 4.
The results from these fits are compared in Table 5. The large χ 2 is due to the significant correlation between the horizontal and vertical errors, which is visible in Fig. 3. Fits with four subtractions (n = 4) do not significantly reduce χ 2 . Results from the N401 and N200 ensembles, which have similar quark masses but different lattice spacings, are shown in Fig. 5 together with thrice-subtracted fits. Agreement between these two ensembles indicates that cutoff effects are also under control in the form factor. Table 5 Results from twice-and thrice-subtracted dispersive fits (see Eq. (3.4)) to the form factor on each ensemble with the three finest lattice spacings. The large χ 2 is due to the significant correlations between E cm and ln Q n .

Conclusions
This work presents an N f = 2 + 1 calculation of the I = 1 elastic p-wave ππ scattering phase shift and timelike pion form factor which addresses systematic errors due to the finite lattice spacing, mixing of higher partial waves, and residual (exponential) finite volume effects.
While we do not perform continuum and chiral extrapolations here, our data can be used for future such extrapolations. Chiral extrapolations of lattice scattering data using unitarized extensions of chiral effective theory [65,[67][68][69][70][71][72][73][74] have been performed, although to date cutoff effects have not been considered. Nonetheless, it is evident in our data that cutoff effects in both the scattering amplitude (shown in Fig. 2) and the timelike pion form factor (Fig. 5) are small with respect to our statistical errors. The coarser lattice spacing in this comparison is a = 0.075 fm. Furthermore, for the scattering amplitude we also check explicitly (in Fig. 2) that finite volume effects are also insignificant at our coarsest lattice spacing. The two volumes used here have m π L = 4.6 and 6.1.
As discussed in Sec. 1, complete extrapolations of the energy dependence of amplitudes are left for future work. However as a necessary ingredient to determine the form factor, we use Breit-Wigner fits to model the energy dependence of δ 1 (E cm ). A summary of the fit results for the resonance mass and coupling are shown in Fig. 6. The CLS ensembles employed here adjust m s as m l is lowered to its physical value such that tr M q = const. is fixed. This is to be contrasted with the more common strategy of fixing m s to its physical value for all m l such as the recent N f = 2 + 1 results in Ref. [7] which employs rooted staggered fermions. Compared with the more standard trajectory, it appears that the slope of m ρ is somewhat flatter here. Although not shown in this work, recent summaries of existing results for the ρ-resonance parameters are found in Refs. [5,6].
In addition to future fits and extrapolations, our results for the timelike pion form factor can be used to extend the vector-vector correlator as described in Ref. [76] and implemented in Ref. [109] without reliance on experiment. This may significantly improve lattice determinations of hadronic vacuum polarization contribution to anomalous magnetic moment of the muon, a HVP μ . Finally, the computational effort expended on the CLS lattices for this work can be largely re-used for other two-to-two amplitude calculations. First work in this direction for Nπ scattering has already appeared in Ref. [26]. The set of ensembles used here will also be augmented by several others at lighter quark masses, include one with L = 6.5 fm at the physical point, which is presented in Ref. [110]. and Nordrhein-Westfalen (MIWF). The USQCD QDP++ library [111] was used in designing the software used for the calculations reported here.

Appendix A. Finite volume energies and scattering amplitudes
In this appendix, we tabulate and plot the finite-volume energies, scattering amplitude, and timelike pion form factor for all ensembles in Table 1. The scattering amplitude and form factor tabulated here employ the truncation to ≤ 1. The C101, D101, N401, N200, D200, J303 ensembles are tabulated in Tables 6,7,8,9,10,11 and plotted in Figs. 7,8,9,10,11,12, respectively.   Fig. 7 for the D101 ensemble. One state in each of the A + 1 (3) and B + 2 (2) irreps which have very large errors have been removed from the plot. Table 8 Results from the N401 ensemble. The form factor is omitted for a single level in the (d 2 ) = A + 1 (4) irrep for which a plateau could not be identified.  Fig. 7 for the N401 ensemble. Bottom row: the timelike pion form factor (left) and the ratio employed in the thrice-subtracted subtracted dispersive fit (right), which is also shown. Table 9 Results from the N200 ensemble.  10. Same as Fig. 9 for the N200 ensemble.  Fig. 11. Same as Fig. 9 for the D200 ensemble.  Fig. 9 for the J303 ensemble.