Hadron-Hadron Interactions from $N_f=2+1+1$ Lattice QCD: The $\rho$-resonance

We present an investigation of the Rho-meson from Nf=2+1+1 flavour lattice QCD. The calculation is performed based on gauge configuration ensembles produced by the ETM collaboration with three lattice spacing values and pion masses ranging from 230 MeV to 500 MeV. Applying the L\"uscher method phase shift curves are determined for all ensembles separately. Assuming a Breit-Wigner form, the Rho-meson mass and width are determined by a fit to these phase shift curves. Mass and width combined are then extrapolated to the chiral limit, while lattice artefacts are not detectable within our statistical uncertainties. For the Rho-meson mass extrapolated to the physical point we find good agreement with experiment. The corresponding decay width differs by about two standard deviations from the experimental value.

the context of vector meson dominance, for a review see Ref. [3].
Therefore, an investigation of the ρ-meson properties from first principles with lattice QCD is highly desirable. However, unstable particles require special care in lattice QCD: interaction properties can only be computed using the by now famous Lüscher method [4][5][6]. With its help, infinite volume scattering properties can be extracted from finite volume energy shifts. In the meanwhile the Lüscher method has been developed further in many directions, for a review see Ref. [7], in particular also for three particle systems, see for instance Refs. [8][9][10][11][12]. For this paper most relevant is the derivation of the formalism in moving frames [13][14][15], which allows one to map out the phase shift at many different scattering momenta, without the need to study different volumes.
For a long time the Lüscher method was difficult to apply to the ρ in realistic lattice calculations, albeit there are some early attempts [16,17]. By now, there are a number of investigations of the ρ-meson from lattice QCD using the Lüscher method [14,[18][19][20][21][22][23][24][25][26]. The first computation with light dynamical up and down quarks can be found in the pioneering work of Ref. [14]. Subsequent investigations focused on different aspects like large operator bases [27] or asymmetric boxes [24].
Recently, a first investigation involving different lattice spacings and a range of pion masses has been performed [26]. However, in the latter reference chiral and continuum extrapolations could not be taken.
With this paper we fill this gap and present a computation of the ρ-meson applying the Lüscher method using gauge ensembles generated with N f = 2 + 1 + 1 dynamical quark flavours by the ETM collaboration at three different lattice spacing values and a wide range of pion masses [28,29]. This allows us to perform a controlled chiral and continuum extrapolation of the ρ-meson mass and width. Note that in Ref. [30] the mass and width of the ρ-meson has been determined on the same gauge configurations, however, using an inverse Lüscher approach based on the vector current only combined with a parametrisation of the pion form factor. This paper is organised as follows: after presenting the lattice action and its parameters in section II, we discuss our methods in section III. In section IV we present our results, the main result being the continuum extrapolated values of M ρ and Γ ρ at the physical pion mass value reading M ρ = 769 (19) MeV , Γ ρ = 129 (7) MeV .
In section V we discuss our results and put them into perspective, followed by a summary in section VI. More technical details can be found in the appendix.  Table I: The gauge ensembles used in this study. For the labeling of the ensembles we adopted the notation in Ref. [28]. In addition to the relevant input parameters we give the lattice volume and the number of evaluated configurations, N conf .

II. LATTICE ACTION
The lattice details for the investigation presented here are very similar to those we used in our previous studies on hadron-hadron interactions [31][32][33][34]. We use N f = 2 + 1 + 1 flavour lattice QCD ensembles generated by the ETM Collaboration, for which details can be found in Refs. [28,29,35].
The parameters relevant for this paper are compiled in Table I: we give for each ensemble the inverse gauge coupling β = 6/g 2 0 , the bare values for the quark mass parameters µ , µ σ and µ δ , the lattice volume and the number of configurations on which we estimated the relevant quantities.
The ensembles were generated using the N f = 2 + 1 + 1 twisted mass fermion action [36][37][38]. For orientation, the β-values 1.90, 1.95 and 2.10 correspond to lattice spacing values of a ∼ 0.089 fm, 0.082 fm and a ∼ 0.062 fm, respectively, see also Table II. The ensembles were generated at so-called maximal twist, which guarantees automatic order O(a) improvement for almost all physical quantities [36]. The corresponding lattice Dirac operator in the light sector reads  tuning of the Wilson quark mass to its critical value is discussed in Ref. [28], where also the Dirac operator for the strange/charm sector can be found, which is not relevant for the remainder of this paper. In the gauge sector the Iwasaki action is used [39,40].
The biggest disadvantage of Wilson twisted mass fermions at maximal twist is the breaking of isospin symmetry. As a consequence, charged and neutral pions are not mass degenerate, with the splitting in the squared masses vanishing like a 2 towards the continuum limit. This pion mass splitting is also about the only quantity where strong effects of isospin splitting have been observed so far [41].
We are going to study the decay ρ 0 → π + π − in a p-wave. In nature, there is no mixing with two neutral pions possible. Even if there is reduced isospin symmetry (only I z is a good quantum number) in the Wilson twisted mass formulation at maximal twist, such mixing is still not possible due to C-symmetry: ρ 0 is C-odd, while π 0 π 0 is C even. Likewise, non p-wave symmetric combinations of π + π − are C-even, excluding also mixings with I = 2, I z = 0 states. Moreover, also a single π 0 is C-even and cannot mix.
However, due to missing isospin symmetry, there are fermionic disconnected contributions to the ρ 0 lattice interpolating operators. These can be shown, like for the neutral pion, to be purely of O(a 2 ). Thus, we drop them from our calculation, as was also done in Ref. [14]. Note that the neutral to charged ρ-meson splitting was found to be negligible [42].
As a smearing and contraction scheme we employ the stochastic Laplacian-Heaviside (sLapH) approach, described in Ref. [43]. Details of our sLapH parameter choices can be found in Refs. [31,32].

A. Scale Setting
The scale setting for the ensembles used here has been performed in Ref. [44] by extrapolating pseudo-scalar meson masses and decay constants to the chiral and continuum limits and using the physical values of M π and f π as inputs. As an intermediate scale the Sommer parameter r 0 /a has been used. The values for the lattice spacings resulting from this procedure can be found in Table II together with the values of r 0 /a for each β-value. The physical value of the Sommer parameter was determined in Ref. [44] on the same ensembles as the value In this paper we are also going to use the Sommer parameter as intermediate lattice scale. In addition to the physical value for r 0 given above we need the physical pion mass value as input.
Here, we use the value of M π in the isospin symmetric limit [45] (consistent with what was used in Ref. [44]) corrected for QED and strong isospin contributions. The values of r 0 /a were not determined by us on the identical set of gauge configurations. Therefore, we use the values given in Table II with re-sampling (parametric bootstrap). M π and its error are treated in the same way. We use the following augmented χ 2 -function, here for the example of errors on r 0 /a to account for the uncertainties of r 0 /a, where P r 0 /a (β) denote the priors for r 0 /a(β) with error ∆r 0 /a(β) for a given β-value taken from Table II. χ 2 is the usual correlated χ 2 -function. In the same way errors on M π and r 0 are included in the fits.

III. METHODS
In this section we summarise the methodology we applied to extract our results.

A. Scattering in Finite Volume
As is well known, the extraction of scattering properties from lattice QCD in Euclidean spacetime and a finite volume requires the application of the so-called Lüscher method [5,6]. It allows one to relate finite volume induced energy shifts to infinite volume scattering properties of n-particle systems in the continuum. The formalism is based on the following determinant equation where M lm,l m is an analytically known matrix function of the lattice scattering momentum k, see below. δ l is the phase shift of the l-th partial wave and the determinant acts in angular momentum space. In the case of pion-pion scattering the lattice scattering momentum k is related to a given energy value E CM in the centre-of-mass (CM) frame and the pion mass M π via Given the scattering momentum on the lattice, Eq. (4) thus yields δ l . In order to map out the dependence of δ l on E CM , as many values of E CM as possible must be extracted from a lattice calculation.
This is most conveniently done by using several CM momenta, as first proposed in Ref. [13]. For given CM momentum p cm , the relativistic energy reads where p cm is, due to the finite volume, quantised as We classify momentum sectors by |d| 2 and use all allowed lattice momenta in each sector up to d 2 = 4. We denote the set of equivalent momenta as By applying a corresponding Lorentz boost we can compute E CM for given W L and p cm . Adopting the notation of Refs. [15,46], it remains to give details for the matrix M from Eq. (4). Its matrix elements are given by with the convenient notation C lm,js,l m represent coefficients which can be expressed using Wigner 3j-symbols, see Ref. [15].
In a finite volume the symmetry group of rotoflections (rotations and space inversions) is reduced from O(3) to a finite subgroup. 1 Because p cm is an invariant, the group is different for each momentum sector.
In general, an irreducible representation restricted to a subgroup does not remain irreducible.
The prescription to decompose an eigenstate of the l-th partial wave is often referred to as "subduction". To introduce notation, assume the irrep D l decomposes into a direct sum of different irreps Γ i each of which appears n i times such that Let Γ ∈ {Γ i }, and label the basis vectors of Γ by α ∈ {1, . . . , dim(Γ)}. The decomposition can be completely described by a set of "subduction coefficients" denoted by s. Given a basis {|l, m | − l ≤ m ≤ l}, the α-th basis vector of the n-th copy of Γ is given by The derivation of subduction coefficients is discussed in appendix A. Applying the subduction to the matrix M from Eq.
The Lüscher formula Eq. (4) remains formally unchanged except for the space it acts in. In the following, we will neglect all partial waves apart from the p-wave. In this case n i = 1 for all i and Eq. (4) simplifies to The contributions of higher odd partial waves have been analysed and found to be negligible [20,24].
While twisted mass breaks parity and thus even partial waves may enter, the effect is suppressed by O a 2 and also neglected here.
In Table III we list the explicit expressions for M Γ 11,11 used in this work.

B. Extraction of Energy Levels
In order to be able to use Eq. (13), we need to extract interacting energy levels for a given lattice irrep Γ as well as the pion energy E π (p). The latter is, as usual, determined from the Euclidean time dependence of two-point functions with operators O π + (t, p) coupling to the charged pion state with momentum p, see below. Note that in our formulation we have M π + = M π − . The spectral decomposition of C π yields In the limit of large Euclidean times only the ground states survives and allows one to extract E π (p) from its exponential decay.
For irrep Γ we define a list of suitable operators O i Γ (t, p), i = 1, . . . , n, which project to irrep Γ for momentum p. Because the eigenvalues of operators from the same momentum sector and irrep are degenerate up to statistical fluctuations, we compute the correlator matrix by averaging over all moving frames connected by an allowed lattice rotation and rows of the irrep The correlator matrix C Γd 2 (t) is then analysed using the standard variational method [52,53] Here, T is the time extent of the lattice and W i the ith energy level to be extracted. t 0 represents the reference time at which the generalised eigenvalue problem (GEVP) is seeded. The correction to Eq. (17) due to excited states reads at fixed t 0 -value [53] Here, ∆W i is the energy difference of W i to the first state not resolved by the correlation matrix.
For a detailed discussion see Ref. [54].

C. Operator Construction
We start with interpolating operators for pions π ± with definite isospin |1, ±1 I : where u and d denote Dirac spinors for an up and down quark, respectively. α, β denote spin and c colour indices, and Γ π = iγ 5 .
For the ρ-meson, we have to construct operators projected to I = 1. A single ρ 0 can be interpolated by the canonical anti symmetric combination of quarks with isospin |1, 0 I : The projection of a given single particle operator O(t, x) to momentum p is performed via and likewise for two particle operators O(t, x 1 , x 2 ) to momenta p 1 , p 2 , respectively, yielding The projection to a given lattice irrep Γ and basis vector α is performed via the so-called subduction procedure described in appendix A.
Apart from excited state contaminations Eq. (18) there are additional so-called thermal state pollutions, which are relevant with finite time extent T , periodic boundary conditions and in the presence of multi-particle states. This was discussed e.g. in Ref. [14] and subsequently e.g. in Ref. [55].
For the case of pion-pion systems with momenta p 1,2 , the leading thermal pollution to a matrix element of the correlator matrix C Γα reads For p 1 = p 2 this is a constant contribution and the time dependence drops out. The thermal pollution ε t vanishes for T → ∞. However, at finite T it can become relevant for t → T /2.
There are, of course, further pollution terms which are exponentially suppressed compared to the one quoted above. Let us now assume E π (p 2 ) > E π (p 1 ) and concentrate on the corresponding, exponentially decreasing term in Eq. (23). This is sufficient because the signal to noise ratio in the relevant correlator matrices is decreasing exponentially with Euclidean time. Therefore, we will have to extract the signal at relatively small t-values where the second, exponentially increasing term in ε t is not yet relevant.
We can deal with this pollution term by applying the so-called weighting and shifting procedure [55]. It amounts to the following transformation of C: with ∆E = E π (p 2 ) − E π (p 1 ). It is easy to see that this transformation leaves the leading, physical exponential dependence unchanged, while the thermal pollution is removed. As an input for the transformation Eq. (24) we use E π (p) determined from single charged pion two point functions at zero momentum combined with the continuum dispersion relation Eq. (6).
We remark here that we have investigated thermal pollutions in some detail in Ref. [33]. However, the corresponding findings are not applicable here, because the signal does not extend to large enough t-values.

E. Phase Shift Curves
Once the energy levels have been determined for all the irreps mentioned above, the phase shift δ 1 is to be determined from Eq. (13). This requires the evaluation of the Lüscher zeta function Z lm (1, q 2 ) in w lm . Z has poles at q 2 -values corresponding to the free, non-interacting two particle energies. The larger the spatial extent L of the lattice, the closer are the interacting energy levels to these poles.
This structure makes the error estimate for δ 1 difficult in cases where the statistical uncertainty of the interacting energy levels is not small enough: when an energy level is compatible with a pole of the Z-function within errors, a proper estimate of the uncertainty of δ 1 becomes impossible.
However, also when this is not the case, such a situation can still be and actually is triggered in some cases during a bootstrap analysis. Since bootstrap replicates are sampled uniformly random with replacement, it is not unlikely to hit a pole of the Z-function, even if the pole is two or three sigma away from the actual energy level.
To circumvent this problem, we use instead of the bootstrap the jack-knife procedure, which can be understood as a linear approximation to the bootstrap. The standard-deviation over jack-knife replicates is per construction a factor of √ N − 1 smaller than the one over bootstrap replicates, where N is the sample size.
One might wonder whether using the jack-knife procedure introduces additional uncertainties, in particular in the vicinity of a singularity of the Z-function. We have compared the jack-knife and bootstrap procedure for all cases, where bootstrap did not show the aforementioned problem.
For all these cases we found excellent agreement for the error estimate between the two methods.
Thus, we conclude that the systematic error introduced by jack-knife is not significant.
With this procedure we then determine δ 1 as a function of E CM using equation Eq. (13). The next step is to determine the ρ-meson mass M ρ and width Γ ρ from these phase shift points. For this purpose we use a relativistic Breit-Wigner functional form which we fit to our data. Here, g ρππ is the ρ to ππ coupling constant. The width is related to g ρππ Eq. (25) allows to extract the mass and width from the phase shift data at a given pion mass.
We remark that Eq. (25) contains several approximations. The resonance must be isolated and narrow. Additionally tan δ 1 has a pole at E CM = M ρ which was rewritten as a rational function where the denominator is a first-order polynomial in k 2 . For M ρ = 775 MeV the predicted width is improve fit quality, but had no significant effect on the final results. [20,21,25] Since N conf is different on all our ensembles, the jack-knife procedure is not easily applied in such a chain of analyses and we take the jack-knife errors as an input to a parametric bootstrap procedure.
Here we generate the parametric bootstrap replicates such as to have the same correlation between E CM , M π and δ 1 as the jack-knife replicates. Then we fit Eq. (25) to our data for E CM , δ 1 and M π with two free parameters g ρππ and M ρ .

F. Pion Mass Dependence
In Ref. [57] the pion mass dependence of the ρ-meson mass has been computed using effective field theory with infrared regularisation. Up to O(M 3 π ) plus the non-analytic term of order M 4 π , the dependence reads To this order the formula contains four unknown parameters, the ρ mass in the chiral limit M 0 ρ and the parameters c 1 , c 2 and c 3 . Using this mass dependence of M ρ and the KSFR relation [58,59], we can try to relate g ρππ to M π up to order M 3 π using Eq. (27) and the SU(2) chiral perturbation theory formula for f π [60] Here, f π is the pion decay constant, f 0 its value in the chiral limit and the parameters M 0 ρ and c i are the ones from Eq. (27). Note that we follow the convention with f π ≈ 130 MeV [61]. In addition we have used the definitions and the usual low energy constant¯ 4 . Values for f 0 and¯ 4 have been computed on the ensembles used here in Ref. [44] f 0 = 121.1(2) MeV ,¯ 4 = 4.7(1) .
We remark that the KSFR relation [58,59] g ρππ ≈ M ρ /f π is fulfilled in nature to very good approximation. However, it is not clear at all whether it can be extended beyond leading order in the pion mass.
In Ref. [62,63], the pion mass dependence of the ρ-meson mass and width has been calculated with the complex mass renormalisation scheme from an effective field theory with explicit contributions corresponding to the ω-meson. It is based on the assumption of vector meson dominance and, thus, model dependent; see also Ref. [64] for details on the model. However, its advantage is that mass and width can be extrapolated in a combined fit. The squared pole position of the ρ resonance, Z = M ρ − i/2 Γ ρ 2 has the following pion mass dependence where Z χ is the pole position in the chiral limit and c χ , g ωρπ are coupling constants. Higher order corrections in M π are known in principle, which also include logarithmic terms. The non-analytic structure in M ρ is identical to the one of Eq. (27).
In order to apply this formula to our lattice data, we re-express it in units of the Sommer parameter r 0 and add an a 2 term, which represents the leading lattice artefacts for the twisted mass formulation at maximal twist. p a 2 is an unknown complex parameter.

A. Pion Dispersion Relation
In order to extract the energy shift, we need the pion energy not only at rest but also in moving frames. As mentioned before, in order to reduce statistical uncertainties we are going to use the relativistic continuum dispersion relation to compute W π (p) from the zero momentum pion mass value. As a check for the validity of this approach we have also computed W π (p) from two-point correlation functions with momentum.
In Figure 1 we compare the measured W 2 π (p) with the prediction of Eq. (31) with M 2 π at zero momentum as input exemplarily for the A40.32 ensemble. Good agreement within errors is observed up to d 2 = 4. This makes us confident that using the dispersion relation is safe.

B. Energy Levels
One of the major uncertainties in our extraction of energy levels of multi-particle correlation functions is caused by thermal pollutions. For the case of two pions with maximal isospin the onset of thermal pollutions in Euclidean time in the correlators is clearly visible. However, due to the exponential deterioration of the signal-to-noise ratio, this is not the case for the correlation functions investigated here. This manifests itself also in the fact that there is no clear difference visible between principal correlators λ(t, t 0 ) derived from C(t, t 0 ) or their weighted and shifted counterpartsλ(t, t 0 ) derived fromC(t, t 0 ). Therefore, we perform the full analysis with and without weighting and shifting and include the difference as a systematic uncertainty in our error budget.
The other major uncertainty in extracting energy levels from lattice correlation functions stems from the choice of fit range. There have been approaches to make this choice more objective by performing a weighted average over many fit ranges, see for instance Ref. [31], which works well for the case of single pions or two pions with maximal isospin. In contrast, for the case in question here, the ρ channel, the weighted average turns out not to be useful.
Therefore, our procedure is the following: we perform the fitting to the principal correlator λ(t, t 0 ) (andλ) by surveying multiple fit ranges [t min , t max ] and selecting a representative one. We enforce a plateau length of at least four points, which must be compatible within errors and have relative errors below 50%. Additionally we require no significant dependence on t max as this would Finally, all other qualities being equal, we preferred larger t max .
In Figure 2 we show an example for the fit range chosen for ensemble A40.32 where d 2 = 1 and irrep Γ = E without weighting and shifting. In the left panel, we show the ratio of principal correlator λ(t, t 0 ) and the single exponential fit model C th (t, t 0 ) = exp(−W (t − t 0 )). Compared to the effective mass, the ratio is more robust numerically. By definition the central value is 1. In the right panel we show for illustration the result of the correlator fit as a red band along with the effective mass m eff (t) = log C(t) C(t + 1) .
As mentioned above, the effects of thermal states are not visible here. The energy level was determined as aW = 0.4412 (26).
In Figure 3 we show the same plots but this time with weighting and shifting. The size of error bars is increased compared to without weighting and shifting, which can be explained by the reduced correlation of neighbouring time slices. For very large t, points are not depicted because they were compatible with zero. For this reason, t max was chosen smaller compared to before. The fit model was modified as described in Eq. (24) and the calculation of the effective mass in the right panel was changed accordingly. The fit result increased by roughly one standard deviation to aW = 0.4463 (23). Whether this results from the independent choice of a fit range or due to not visible but barely significant thermal states remains hard to decide. By including this difference as a systematic error we are confident that we keep control of both major sources of systematic uncertainties.
In Figure 4    cleaner: in the former case we might be plagued with thermal state pollutions, while in the latter case the fit range might be chosen incorrectly due to noise.
Therefore, we decided to use the weighted mean over results w/o and w/ thermal state removal.
In addition we include the difference ∆Q Y between the weighted mean and w/o and w/ thermal state removal into the error by rescaling the bootstrap distribution with a factor [33] Here, ∆x is the statistical uncertainty of the weighted mean and Y ∈ {w/o, w/}.
All results for M ρ and g ρππ determined by this procedure w/ and w/o thermal state removal are compiled in We have two groups of ensembles with all identical parameters apart from the volume. These are ensembles A40.24 and A40.32 as well as B35.32 and B35.48, which we can use to investigate residual finite volume effects in our results for M ρ and Γ ρ .
In Figure 6 we compare in the left panel the phase shift points for A40. 24   fits happen to result in slightly different values for the resonance parameters, deviations are below the 2σ level and do not show a systematic ordering with volume, see Table IV and Table V.
Thus, the weighted average with error including the systematic uncertainty from thermal state removal should also safely include residual effects from finite volume.
There are a few ensembles where the Breit-Wigner type fits to the phase shift points are problematic. On the one hand this is the case for ensemble with the heaviest pion mass A100.24. The width approaches zero, which leaves the fits little freedom; a fact reflected by the untrustworthy On the other hand, unfortunately the fit on D15. 48 Table V: We give aM π , the finite size correction factor K M π , the ρ width aΓ ρ computed from aM ρ and g ρππ using Eq. (26) w/ and w/o thermal state removal, and the weighted average as explained in the text. In addition we give the reduced χ 2 -values of the corresponding fits to the phase shift data.

D. Chiral extrapolation
We first consider M ρ and g ρππ . In the left panel of Figure 7 we show r 0 M av ρ , in the right one g av ρππ , both as a function of (r 0 M π ) 2 . Note that the error on r 0 /a is not included in the plot, because it is 100% correlated for all data points of the same β-value. Colours and symbols encode the three lattice spacing values. The black diamonds represent the corresponding experimental values. The first observation is that lattice artefacts are not resolvable given our current level of statistical uncertainty. Overall, M ρ appears to show a rather linear dependence on M 2 π , a bit less so g ρππ . The values for aM π can be found in Table V. For the following extrapolations we correct aM π for finite size effects by applying a correction factor K M π computed in Ref. [44], which can also be found in Table V. Next we have tried to fit the pion mass dependence of M av ρ and g av ρππ combining Eqs. (27) and (28) up to the order M 3 π . However, such a fit did not result in convincing results. Even though the chiral log in g av ρππ stemming from f π somewhat compensates the term c 1 M 2 π , a satisfactory description of the data for both the mass and the coupling could not be achieved.  As visible, the extrapolation overestimates both the ρ mass and the coupling at the physical point compared to experiment.
Therefore, we turn now to fits of mass and width using Eq. (30) for the complex valued variable Z. As described in section III, we extrapolate M ρ and Γ ρ to the physical point combined in As we also mentioned already, the error analysis for this fit is performed using the parametric bootstrap procedure maintaining the correlation among M ρ , Γ ρ and M π . We use 1500 bootstrap samples and the values for r 0 /a for the different lattice spacings were resampled from the values compiled in Table II. The actual fit function reads The fit parameters are the following: p 1 and p 2 represent the real and imaginary parts of r 2 0 Z χ and p 3 represents C χ , furthermore p 4 ≡ g 2 ωρπ /(24πr 2 0 ) and p 5 and p 6 parametrise the real and imaginary part of the a 2 lattice artefacts. p r 0 /a is one fit parameter per lattice spacing value for r 0 /a accompanied by a corresponding prior P r 0 /a . Thus, we have in total 6 real-valued free fit parameters.
In the fit we include only the ensembles with the largest volume per pion mass value, i.e. A40.24 and B35.32 are not included in the fit. We do not include ensemble D15.48 in the fit, for reasons Open symbols are not included in the fit. In the right panel g av ρππ is shown also as a function of (r 0 M π ) 2 . The lines with error bands represent independent fits to the data. mentioned above. Moreover, we include only data points with M π ≤ 420 MeV, which excludes ensembles A80.24 and A100.24.
The best fit parameters can be found in Table VI together with the reduced χ 2 -value. We give the best fit parameters for fits with and without lattice artefacts included. Clearly, p 5 and p 6 , which parametrise the a 2 effects in Z are compatible with zero. Also, the remaining parameters do not change significantly with and without a 2 artefact included in the fit. The χ 2 -values for these fits are all a bit too large, indicating a tension in the data in particular Our final result for M ρ and Γ ρ taken from the fit without a a 2 effects included reads corresponding to g ρππ = 5.5(1) .
In addition we find Note that these also correspond to Breit-Wigner parameters determined experimentally from e + e − reactions. The deviation to other reactions can be of the order of 10 MeV. We observe rather good agreement for M ρ , while our value for the width is slightly too low. This is also visible in Figure 9, where we plot the experimental phase shifts from Ref. [2] and compare them to the phase shift curve we obtain using the final values from Eq. (34) and then again assuming the Breit-Wigner form from Eq. (25).
However, this good agreement should be taken with caution. First of all our extrapolation form for M ρ and Γ ρ is not model independent. This is in particular important, because the curvature needed to obtain an M ρ -value close to the experimental one comes from constrains due to Γ ρ . This, as discussed earlier, manifests itself also in a bit too large χ 2 -values in the chiral and continuum fits. Moreover, the ensemble with the lightest pion mass included in the fit is B35.48 with a pion mass of about 300 MeV. Thus, the extrapolation to the physical point is quite long. In addition we have assumed that we can perform a Breit-Wigner type fit to all the phase shift data, which is an approximation. This might also be the reason for the too low value of Γ ρ compared to experiment.
We are currently working on an alternative extrapolation using the inverse amplitude method which might allow us to perform the chiral extrapolation even more reliably [65][66][67][68][69]. Our fitted value for g ωρπ Eq. (36) is in the right ballpark, when compared to the numbers given in Refs. [62,63], where 16 GeV −1 is quoted. From Refs. [3,70] one finds g ωρπ = ±20.7 GeV −1 in very good agreement with our value.
Finally, our determinations of mass and width rest on the assumption that all partial waves apart from = 1 are negligible. This assumption is supported by previous lattice investigations of the ρ meson, but has not been checked by us yet.
On the other hand, our results for M ρ and Γ ρ make a consistent extrapolation to the physical point and to the continuum limit feasible for the first time . This is to a large extent due to merely invisible lattice artefacts in our results. We have different volumes available with otherwise fixed parameters, which allow us to argue that residual finite volume effects are not a dominant source of uncertainty in our results.
In Figure 10 we compare results for M ρ and g ρππ from various lattice collaborations with N f = 2 + 1 or N f = 2 + 1 + 1 dynamical quark flavours. We observe that there are probably lattice artefacts in some of the results for M ρ , in particular in the results from Andersen et al. [26] and from the Hadron Spectrum Collaboration [20]. For g ρππ uncertainties are in general larger and within these large uncertainties the agreement among different lattice collaborations is reasonable.
However, leaving aside lattice artefacts, one could be tempted to conclude from Figure 10 that M ρ is rather linear in M π , very similar to what is observed for the nucleon mass [71]. In fact one finds M ρ = 680 MeV+0.6M π to a good approximation by fitting only the data by Fu and Wang [23] together with our data, which represents yet another version of the "ruler" plot. From an effective field theory point of view this cannot be the correct pion mass dependence and future results will hopefully shed light on this puzzle.
We can also compare to the results of Ref. [30], where M ρ and g ρππ have been determined on the same ETMC ensembles we used, however, using the inverse Lüscher method based only on the , HadSpec [20,22], PACS-CS [19] as well as the experimental value [61].
vector current based on a parametrisation of the pion form factor. Their continuum extrapolated values for M ρ and g ρππ at the physical point are consistent with ours.

VI. SUMMARY
We have presented an investigation of the ρ-meson properties using lattice QCD with N f = 2 + 1 + 1 Wilson twisted mass quarks at maximal twist. With three values of the lattice spacing and a range of pion mass values we could perform chiral and continuum extrapolations of ρ-meson mass M ρ and width Γ ρ with better control than previously possible. The latter two quantities have been determined on our ensembles using a Breit-Wigner type fit to phase shift data assuming that partial waves with ≥ 3 are negligible.
The phase shift curves have been determined applying Lüscher's method using moving frames up to d 2 = 4 and all available lattice irreducible representations. Our final result reads M ρ = 769 (19) MeV , g ρππ = 5.5(1) , Γ ρ = 129(7) MeV , which is determined from a combined continuum and chiral extrapolation of M ρ and Γ ρ . Systematic errors from thermal state pollutions, the chiral and the continuum extrapolation should be covered by the error we quote. M ρ is very close to its experimental value, the width is about two sigma too low. The agreement of our data for M ρ with previously published lattice results is satisfactory.
It is clear that more work is needed to better estimate the width, which likely suffers from e.g.
the use of a Breit-Wigner type fit to the phase shift data. Therefore, we are currently working on using the inverse amplitude method to directly describe the pion mass dependence of the phase shift curves [69]. This should alleviate systematic uncertainties in our current analysis.  defines a projector, where D Γ denote the irreducible representation matrices and α, β ∈ {1, . . . , dim(Γ)} are arbitrary but fixed. We refrain from discussing the modifications for nontrivial multiplicities. We use the Schönflies notation and follow the conventions for D Γ used in crystallography [82], conveniently implemented in Maple [83].
At rest, the symmetry group is the octahedral group O h . Choosing a non-zero CM-momentum p cm further reduces the relevant symmetry group to the "little group" LG(p cm ) ≡ {g ∈ O h ,R g p cm = p cm } , which leaves p cm invariant. The relevant little groups are listed in Table VII.   e i(x 1 ·( 1 2 p cm +R g q)+x 2 ·( 1 2 p cm −R g q)) D Γ (R g ) * αβ O π + (x 1 ) O π − (x 2 ) .

Appendix B: Correlation Coefficients
In the following table we compile the correlation coefficients of the chiral fit without lattice artefacts included in the fit: fit function is thus Eq. (33) without the term proportional to (a/r 0 ) 2 .
The bare data can be made available upon request.