The $\rho$-resonance with physical pion mass from $N_f=2$ lattice QCD

We present the first-ever lattice computation of pi pi-scattering in the I=1 channel with Nf=2 dynamical quark flavours obtained including an ensemble with physical value of the pion mass. Employing a global fit to data at three values of the pion mass, we determine the universal parameters of the rho-resonance. We carefully investigate systematic uncertainties by determining energy eigenvalues using different methods and by comparing inverse amplitude method and Breit-Wigner type parametrizations. Overall, we find mass 786(20) MeV and width 180(6) MeV, including statistical and systematic uncertainties. In stark disagreement with the previous Nf=2 extrapolations from higher than physical pion mass results, our mass value is in good agreement with experiment, while the width is slightly too high.

Introduction.-Quantum Chromodynamics (QCD)the theory of strong interactions -gives rise to a fascinating plethora of hadronic states: mesons and baryons. A theoretical understanding of these states from first principles requires a non-perturbative method, as provided by lattice QCD. As most of the hadrons are not stable under the strong interaction and decay, such an investigation must include resonance and interaction parameters.
One very prominent such state is the so-called ρresonance (ρ(770)), which decays predominantly in a pwave into two pions with isospin I = 1. It is experimentally observed as a peak in cross-sections at an energy of M ρ ∼ 775 MeV with width Γ ρ ∼ 150 MeV [1]. The corresponding p-wave phase-shift curve is a prime example for a resonance phase-shift, as can be seen from Fig. 1, where the experimentally measured phase-shift δ 1 [2,3] is depicted as a function of the center-of-mass energy. Via vector meson dominance, the ρ plays a fundamental role in our theoretical understanding of many processes [4] and, since it is well investigated experimentally, it represents a benchmark resonance state for lattice QCD simulations.
The ρ-resonance has been investigated in lattice QCD previously [5][6][7][8][9][10][11][12][13][14], most recently in Ref. [15] including a continuum and chiral extrapolation. The latter was needed because the lattice simulations were performed at unphysically large values of the pion mass. One of the interesting conclusions from Ref. [15] is that the chiral extrapolation is difficult, even though there is guidance from effective field theories. The main reason for this was the lack of ensembles with sufficiently light pion mass, the lightest being at 230 MeV. In addition, there is an ongo- We compare experimental data [2,3] to our prediction from three global IAM fits, see text. The blue band visualises the statistical and fitting uncertainties, the gray band includes in addition the estimated lattice artefacts added in quadrature.
With this letter we fill the gap to realistic pion mass values by including an ensemble directly at the physical point for the first time. We present ρ-resonance parameters determined on three N f = 2 lattice QCD ensembles generated by the Extended Twisted Mass Collaboration (ETMC) [20] with pion masses between 132 and 340 MeV. These allow us to interpolate rather than extrapolate to the physical pion mass and, therefore, to directly compare to experiment. Moreover, by comparing to other lattice investigations we can shed new light on the question of the importance of theKK threshold for the ρresonance phase-shift.
The main result of this letter is summarised in Fig. 1, where, in addition to the experimental data for the phaseshift, the blue band shows our interpolation for the phaseshift at the physical pion mass value. The width of the band represents fitting and statistical uncertainties. The gray band includes in addition an estimate of lattice artefacts.
Lattice Computation.-The results for the ρ-resonance properties presented in this letter are based on gauge configurations generated by the ETMC at a single value of the lattice spacing, a = 0.0914 (15) fm [20]. These employ the Iwasaki gauge action [21] and two dynamical mass-degenerate flavours of Wilson twisted mass clover fermions at maximal twist [22,23]. With this action, physical quantities are O(a) improved [24] such that discretisation effects only appear at order a 2 in the lattice spacing. The three ensembles considered for this letter are compiled in Table I together with the lattice volume (L/a) 3 × T /a, the number of configurations N conf and the pion mass value in physical units. For more details we refer to Ref. [20].
On these ensembles, we compute center-of-mass energy levels E Γ cms for irreducible representations (irreps) of the lattice rotational symmetry group Γ. We follow the procedure detailed in Ref. [15] and compute Euclidean correlation matrices C Γ,p 2 averaging over all equivalent momenta p. The operators . .) t are chosen to project to irrep Γ for total squared momentum p 2 . The list of irreps Γ considered is T 1u , A 1 , E, B 1 , B 2 up to d 2 = 4 with p = 2πd/L.
The basis operators used to construct the operators O Γ,P are two pion and single vector meson operators with Γ ρ ∈ {iγ i , γ 0 γ i }.
We apply the generalised eigenvalue method (GEVM), i.e., solve the generalised eigenvalue problem [25,26] for eigenvalues λ (n) (t) and eigenvectors η (n) , where n labels the contributing states. Energy levels of these can be determined from the exponential fall-off of λ (n) (t) at large t. In addition, we apply the so-called Prony generalised eigenvalue method (PGEVM) in form of a matrix pencil on top of the GEVM [27] to reduce excited state contaminations.
The confidence in our energy eigenvalue extractions is increased by employing the following three methods: A1: direct fit to each λ (n) (t) using a fit range chosen by eye.
A2: direct fit to each λ (n) (t) using the fit range which yields the fit with the best p-value.
A3: fit to the principal correlator of the PGEVM obtained from λ (n) (t) [27] using a fit range chosen by eye.
Energy eigenvalues for which the three methods do not yield consistent results are discarded. To account for residual deviations, we perform the resonance parameter determinations based on energy eigenvalues obtained using each method and take the maximal difference between the resulting parameters as a systematic uncertainty. When considering multi-particle operators with periodic boundary conditions, the corresponding correlation functions are polluted by contributions from the so-called thermal states and there exist several methods to reduce or remove these. However, in Ref. [15] we have shown that in the I = 1 channel, the extracted energies agree within errors with or without thermal state subtraction if the fit-range is chosen carefully. We have checked that this is the case also here and thus use energy levels extracted without thermal state subtraction. Like in Ref. [15], we use so-called stochastic Laplacian Heaviside smearing [28] with algorithmic parameters identical to Ref. [29]. For the determination of the pion decay constant on the same gauge configurations, we also employ local time slice sources and the so-called one-end-trick, for details see Ref. [30].   Phase-shift determination.-The discrete and real valued lattice energy levels E Γ cms are mapped to the infinite volume scattering quantities using Lüscher's method [26,32,33]. In case of the ρ-meson and under the assumption that higher partial waves can be neglected, the p-wave phase-shift δ 1 is related to the energy levels via where M Γ is an algebraically known matrix function [15,34,35] of the lattice scattering momentum k 2 (E 2 cms ) = E 2 cms /4 − M 2 π and the pion mass, M π . Note that Eq. (3) is valid below inelastic threshold (4M π ) only. This represents a limitation in particular for the ensemble cA2.09.48, where only five energy levels lie below this threshold for our L-value. In a more general sense this also implies that independently of the number of points below threshold, the resonance region of the ρ-meson can never be mapped out using Lüscher's method only, because 4M phys π < M ρ . Given only discrete values of E cms , one needs to parameterize the scattering amplitude as a function of a continuous E cms . One example for such a parametrization is a simple Breit-Wigner (BW) form with M ρ the ρ-resonance mass, g ρππ the ρ − ππ coupling and s the center-of-mass energy squared.
Supplementary to the experimental measurements, additional information about the dynamics of the ππ system resides in the pion-mass dependence, which can be explored with lattice calculations. Being in the unique position of having data at the physical, as well as heavier than physical pion mass values, we use the IAM parametrization of the scattering amplitude [36][37][38]. This approach preserves unitarity exactly, has the correct pion mass dependence up to next-to-leading order (NLO) in chiral perturbation theory [39,40] and fulfills further nonperturbative constraints on the chiral trajectory [41].
In IAM, the phase-shift δ 1 is parameterized as (for more details see Ref. [18]) where T 2 denotes the leading chiral order amplitude and T 4 the NLO one without s-channel loop diagrams. The two-meson loop in dimensional regularization is denoted by J(s). The corresponding amplitude is regularization scale independent and depends on one combination of low-energy constants (LECs) [39]l 12 :=l 1 −l 2 as well as the pion decay-constant in the chiral limit (f 0 ). Note that both T 2 and T 4 are expressed in terms of M 2 π /(4πf 0 ) 2 . The expressions (4) and (5) are fitted directly to the energy eigenvalues using Eq. (3) without computing the phase-shifts as an intermediate quantity, or performing any scale setting. Fit parameters are the BW parameters or aforementioned LECs, depending on the considered fit form. In both cases M π is fitted to the lattice values including finite size corrections as described in Ref. [42]. Additionally, in the case of the IAM we also fit f π with respect to a further constant, Λ 4 , related to the NLO LECl 4 [39]. For details see Ref. [31]. In all fits we take full account of correlations and compute statistical uncertainties using the bootstrap. All bare data is publicly available in a data repository [43].
Results.-Here we present and discuss mainly the results obtained with method A1 to estimate energy levels if not mentioned otherwise. For methods A2-3 see [31]. In Fig. 2 we show k cot(δ 1 ) as a function of the centerof-mass energy E cms both in units of the pion mass for the three ensembles separately. The solid red lines with 1σ error band correspond to the best fits of Eq. (5) Table II. The complex ρ-resonance pole position E ρ = M ρ + iΓ ρ /2 is given in Table III, [20]. The results of the single (IAM, BW) and global (IAM) fits are depicted in Fig. 3 together with the physical result [1] and previous N f = 2 lattice determinations [5,6,9,44]. The complex pole positions of our global fit are visualised in Fig. 4  Finally, in Fig. 1 we show the experimental phase-shift data [2,3] as a function of E cms and compare to our global IAM fit prediction for the physical pion mass value. For the latter we plot in blue the envelope area of all the error bands of the three global IAM analyses A1-3, thus, visualising statistical and IAM uncertainties. The grey band includes our estimate of the lattice artefacts added in quadrature.
Discussion.-First, we observe very good agreement between BW and IAM fits on the three ensembles separately, as can be see in Tables II and III, with smaller errors in the IAM fits. This is even the case on the physical point ensemble, however, with much too low pole mass and width. The latter can be attributed to only five points being included in the fit, with much of the curvature triggered by a single data point, because all the points are in a region where δ 1 is close to zero.
This emphasises the importance of including ensembles with larger than physical pion mass value in the analysis: only on those ensembles can the Lüscher method cover the resonance region fully.
Interestingly, if one were to ignore the inelastic thresh-old on the physical point ensemble and use Eq. (3) to obtain values for δ 1 , the such obtained phase-shift points compare reasonably well with the experimental data [31]. This is due to the fact that ρ → ππ is almost elastic up to 1 GeV, which is actually also the assumption to obtain the experimental phase-shift points (see also Ref. [45]). The three global IAM analyses based on energy determinations A1-3 agree very well with each other (see Tab. 1 in [31]) and we conclude that they lead to consistent chiral extrapolations. Though the physical point ensemble cA2.09.48 might not, due to the few energy levels, add much information to the global fit, it is very important, because it anchors the fit at slightly below the physical pion mass value. This turns our determination of the ρ-resonance properties into an interpolation, making it much more reliable.
Below, we take the maximum of the ± statistical errors as our statistical uncertainty, the maximal deviation between the three analyses A1-3 as a systematic uncertainty. In addition we assign a generic 2.5% uncertainty to (undetermined) discretisation artefacts, which are generically of order a 2 Λ 2 QCD . Therefore, we quote as our final result the results from the global IAM fit with analysis method A1, reading Note that our parametrically estimated lattice artefacts are the largest source of uncertainty in our results.
Compared to other analyses of N f = 2 lattice QCD data for the ρ-meson, our results are much closer to the values quoted in the PDG. In particular, we obtain a slightly larger value for the pole mass, which is in contrast to the claim made in Ref. [16] that the missing KK channel pushes this value downwards in N f = 2 QCD. From our results we can only conclude that the influence of the missing strange quark in our ensembles is marginal. There are other effects, most likely lattice artefacts, which are able to explain the findings in Ref. [16], see also recent analyses [17][18][19]. Note also that a different scale setting procedure (see Ref. [20]) would not change our results sufficiently.
Compared to Ref. [15] with N f = 2 + 1 + 1 dynamical quark flavours, we have presented in this letter a significantly better controlled chiral extrapolation thanks to the included physical point ensemble. While the pole mass is similarly close to the experimental value, our width is larger than the physical one, whereas in Ref. [15] a lower value was found. A final estimate will require a continuum extrapolation at the physical pion mass value.
Conclusion.-We have presented a lattice QCD analysis of the ρ-resonance including an ensemble with slightly lower than physical pion mass value for the first time. This allows us to interpolate to the physical point using the inverse amplitude method including the pion mass dependence up to NLO in the chiral expansion.
With all our uncertainties added in quadrature, our results for the ρ-meson mass and width read While M ρ agrees well with the PDG value, the width is too large by 20%. The low energy constants are in very good agreement with the corresponding FLAG lattice averages [46]. This is not necessarily expected, since the IAM resums higher order effects due to unitarisation.
Eventually, this computation needs to be repeated with N f = 2 + 1(+1) dynamical quark flavours, several values of the lattice spacing and physical point ensembles included. Particular emphasis should also be on different spatial volumes at the physical point.
Acknowledgments.-We thank all members of the ETMC for the most enjoyable collaboration. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer JUQUEEN [47] and the John von Neumann Institute for Computing (NIC) for computing time provided on the supercomputers JURECA [48] and JUWELS [49] at Jülich Supercomputing Centre (JSC). We thank Ulf-G. Meißner for useful comments on the manuscript. This project was funded in part by the DFG as a project in the Sino-German CRC110. FP acknowledges financial support from the Cyprus Research and Innovation Foundation under project "NextQCD", contract no. EXCELLENCE/0918/0129. The open source software packages tmLQCD [50][51][52], Lemon [53], QUDA [54][55][56], the hadron package [57] and R [58] have been used. MM thanks R. Brett and A. Alexandru for assistance in implementation of FinVol routines.

Fit strategy and results
We employ two distinct parametrizations of cot δ 1 as discussed in the main part of the article. The first one, of Breit-Wigner type, contains two parameters M ρ and g ρππ , which are determined in a fit to energy eigenvalues using Lüscher's method [26,32]. The width Γ ρ is related to M ρ and g ρππ via The fit is performed by minimizing the fully correlated χ 2 with respect to 3 free parameters, i.e., {M π , M ρ , g ρππ }.
Since this parametrization does not include the pion mass dependence, only single fits to each ensemble are performed.
The results for the dynamical parameters are given in Table II for method A1 and in Table IV for all methods A1-3.
A second analysis is based on the so-called inverse amplitude method, see the main text for the explicit parametrization. As discussed there, this method currently allows for the most reliable extrapolation of ππ dynamics along the energy and pion mass directions. However, to avoid biases, a careful discussion of its implementation with regard to scale setting is in order. At the level of our fits, we completely avoid the need for scale setting by expressing leading and next-to-leading chiral order amplitudes (T 2 andT 4 ) in terms of M π /f 0 . f 0 is then related to f π -for which we have data -using the corresponding NLO ChPT relation involving Λ 4 . In addition to f 0 and Λ 4 , we are left with only one combination of free parameters, namely (l 1 −l 2 ). In individual fits the set of free parameters is given by {aM π , af 0 , (l 1 −l 2 ), aΛ 4 }. This allows to fit the energy eigenvalues together with the finite volume pion mass (M π (L)) and decay constant (f π (L)) by relating for with K 1 denoting the Bessel function of the second kind. Multiplicities m(n) and further definitions can be found in Ref. [42]. The fit is performed by minimizing the fully correlated χ 2 , expressing everything directly in lattice units. The global fit is performed in a very similar way by simply extending the set of free parameters to include two further pion masses, i.e., {aM 1 π , aM 2 π , aM 3 π , af 0 , (l 1 −l 2 ), aΛ 4 }. Note that all other parameters are pion mass independent and, thus, are common to all three ensembles. The results of all fits are given in Table IV for all methods A1-3 (see also

Extrapolation to the physical point
The IAM-based global fit strategy for extracting complex pole positions allows for an extrapolation to the physical point. We define the latter consistently using the current FLAG [46] values as M phys π = 135 MeV and f phys π = 130.41 MeV. Note that since cot δ IAM 1 is parameterized byl 1 −l 2 and M π /f 0 , the only unknown quantity required to study the extrapolation to the physical point is M phys π in lattice units. We found that the most consistent way to fix the latter is to solve for aM phys π with fitted af 0 and aΛ 4 as input. a is then obtained by setting M phys π = 135 MeV. Note that this relation is valid up to the next-to-leading chiral order. We perform this for each of the bootstrap samples and methods A1-3 separately. Besides phase-shifts and pole positions at the physical point, this also allows one to extract f 0 in physical units as well as the low-energy constant of interestl 4 = 2 log(aΛ 4 /(aM phys π ).

Fits above inelastic thresholds
Depending on the analysis method (A1-A3), the spectrum obtained from the cA2.09.48 ensemble contains up to 30 energy eigenvalues. The value of having data close to the physical pion mass is clear, however, only ∼ 5 energy eigenvalues lie below the inelastic (4M π ) threshold. To avoid additional bias only these data points are used in the main part of this letter.
For completeness, we also attempted to fit all the existing data using the two-body Lüscher formalism, simply to see whether or not 4π and 6π interactions contribute a sizable amount to the isovector ππ interaction. It is known from ChPT that four-particle effects are highly suppressed [45], which might justify such an approach.
In addition, while approximate, this analysis can provide an estimate of further systematic effects, beyond what is discussed in the main text.
To avoid mixing with other systematic effects, we perform such exploratory fits separately for each of the used parametrizations (BW and IAM) for each of the methods A1-A3. To this end, we either include all energy eigenvalues or restrict to regions below 6M π or 4M π , respectively. The results of these fits are depicted in Fig. 5. In reassuring agreement between all cases we find that fits of data above the 6π threshold lead to very poor correlated χ 2 values, possibly hinting at sizable effects due to the 6π channels, either on the lattice (e.g., operator basis size), or the analysis part (e.g. 6-body quantization condition). On the other hand, fits up to the 6π threshold converge to values of correlated χ 2 dof 2 for both IAM and BW parametrizations. This suggests that the effects due to 4π interactions might be indeed suppressed in this channel.
While further interpretation of these observations needs a more careful study, we point out that all fits lead to similar complex pole positions of the ρ-meson of E ρ ≈ (827 − 86i, 856 − 125i, 821 − 92i) MeV, with the three values