Hadron-Hadron Interactions from $N_f=2+1+1$ Lattice QCD: isospin-2 $\pi\pi$ scattering length

We present results for the $I=2$ $\pi\pi$ scattering length using $N_f=2+1+1$ twisted mass lattice QCD for three values of the lattice spacing and a range of pion mass values. Due to the use of Laplacian Heaviside smearing our statistical errors are reduced compared to previous lattice studies. A detailed investigation of systematic effects such as discretisation effects, volume effects, and pollution of excited and thermal states is performed. After extrapolation to the physical point using chiral perturbation theory at NLO we obtain $M_\pi a_0=-0.0442(2)_\mathrm{stat}(^{+4}_{-0})_\mathrm{sys}$.


Introduction
Quantum Chromodynamics (QCD) describes, beyond the mass spectrum of stable and unstable hadrons, also their interaction due to the strong force. While originally thought to be not accessible to lattice QCD [1], Martin Lüscher was able to relate the energy spectrum in finite volume to the infinite volume interaction properties in a series of seminal papers [2,3,4,5]. During the last years, extensions of these finite volume techniques have been worked out and they have been applied to a number of interesting systems.
Lüscher's original formalism dealt with two massive bosons in a finite cubic box below the inelastic scattering threshold. It was originally formulated for the center of mass frame. In order to enlarge the applicability of Lüscher's method, extensions have been developed over the years which include: moving frames with non-vanishing center of mass momentum [6,7,8,9], asymmetric boxes [10,11,12], twisted boundary conditions [13,14,15,16], and more. Most of these aim to circumvent the poor resolution in momentum space due to the quantisation of the three-momentum in finite volume. Generalisations beyond the inelastic threshold for the two particle cases are also discussed [17,18,19,20,21,22] and various groups are now working on the more difficult three-particle scenario [23,24,25,26,27,28,29].
In this first paper of a planned series of papers we investigate ππ scattering in the isospin-2 channel. ππ scattering in the I = 2 channel is technically the easiest case to consider in lattice QCD, because so-called fermionic disconnected contributions are absent. Moreover, the predicted quark mass dependence in chiral perturbation theory (χPT) at next-to-leading order (NLO) is governed by only one low energy constant, and so far lattice data showed surprisingly little deviations even from the leading order (LO) χPT predictions, which is parameter free. Isospin-2 ππ scattering is, therefore, an important benchmark system to compare with non-lattice methods and an important test case for future investigations of hadron-hadron interactions, while simultaneously of phenomenological interest.
Consequently, it has been computed in lattice QCD previously [30,31,32,33,34] using N f = 2 or N f = 2 + 1 dynamical fermions, for a recent review see Ref. [35]. In this paper we extend this list by using N f = 2 + 1 + 1 dynamical quark flavours for the first time. The analysis relies on Wilson twisted mass fermions [36] at maximal twist [37] based on gauge configurations provided by the European Twisted Mass Collaboration (ETMC) [38,39]. This allows us to investigate a range of pion masses from 250 to 500 MeV and discretisation effects by using three values of the lattice spacing. For estimating finite volume effects we have several ensembles at our disposal with all parameters fixed but the volume.
We apply different strategies to determine scattering parameters from the data using Lüscher's finite volume method. This gives us an estimate of residual finite volume effects on the scattering length for other ensembles where we do not have multiple volumes at hand.
The main result is an extrapolation of M π a 0 to the physical point utilising chiral perturbation theory at next-to-leading order. We carry systematic effects from different ensemble β aµ aµ σ aµ δ (L/a) 3 × T /a N conf A30. 32 1  Table 1: The gauge ensembles used in this study. For the labelling of the ensembles we adopted the notation in Ref. [38]. In addition to the relevant input parameters we give the lattice volume and the number of evaluated configurations, N conf .
fit ranges through the entire analysis chain. In addition we investigate the stability of the extrapolation by different cuts in the pion mass values. The paper is structured as follows: in section 2 we review the fermion action, the Laplacian-Heaviside smearing technique and the interpolating operators we use to construct the correlation functions. Lüscher's finite volume method is explained in section 3. Our analysis and the treatment of all sources of systematic effects are worked out in section 4. The final results are discussed in section 5 and summarised in section 6.

Lattice action
The sea quarks are described by the Wilson twisted mass action with N f = 2 + 1 + 1 dynamical quark flavours. The Dirac operator for the light quark doublet reads [36] where D W denotes the standard Wilson Dirac operator and µ the bare light twisted mass parameter. τ 3 and in general τ i , i = 1, 2, 3 represent the Pauli matrices acting in flavour space. D acts on a spinor χ = (u, d) T and, hence, the u (d) quark has twisted mass +µ (−µ ). For the heavy unitary doublet of c and s quarks [40] the Dirac operator is given by The bare Wilson quark mass m 0 has been tuned to its critical value [41,38]. This guarantees automatic order O (a) improvement [37], which is one of the main advantages of the Wilson twisted mass formulation of lattice QCD.
The splitting term in the heavy doublet Eq. 2 introduces flavour mixing between strange and charm quarks which needs to be accounted for in the analysis. However, this is only important for quantities involving valence strange and charm quarks.

Stochastic Laplacian Heavyside Smearing
On the long term we plan to address as many observables as possible. Hence we have decided to use the stochastic Laplacian Heavyside smearing (sLapH), which is quite generally applicable and which is described in detail in Refs. [42,43]. It represents a smearing method based on the covariant 3-dimensional Laplace operator where a, b are colour indices, x, y space-time coordinates, U denote smeared gauge links, and the sum is taken over the three spatial directions. The gauge fields in the Laplace operator, Eq. 3, are smeared with three steps of three dimensional HYP smearing [44] with smearing parameters α 1 = α 2 = 0.62, independently of the volume and lattice spacing. The Laplace operator can be decomposed as follows where V ∆ represents the matrix of all eigenvectors and Λ ∆ is a diagonal matrix containing the eigenvalues. Colour, Dirac and space-time indices are suppressed. The smearing matrix then reads where V s only contains eigenvectors corresponding to eigenvalues smaller than a cut-off σ 2 s . These eigenvectors span what is called the LapH sub-space. The value σ 2 s can be chosen in a way that excited state contaminations are maximally suppressed. In Ref. [43] it has been found that σ 2 s = 0.33 is optimal for excited state suppression independent of the interpolating operator. On a 48 3 × 96 lattice, for example, this amounts to about 960 eigenvectors per time slice. However, our tests show that it is better to have less eigenvectors for these large lattices, namely 660 on 48 3 × 96 and 220 on 32 3 × 64, to suppress excited state contaminations at early Euclidean times.
The building blocks of observables are quark lines, Q, which can be written in the LapH framework as where the middle part, P = V † s Ω −1 V s , is called perambulator with Ω = γ 4 M and M being the Dirac operator. To generate a perambulator it is necessary to invert the Dirac operator for every eigenvector in V s . It follows that for an all-to-all perambulator, since the Laplace operator is diagonal in time and Dirac space, this procedure needs to be repeated for each time slice and Dirac component. Taking the 48 3 × 96 lattice as an example this would result in a total number of 253440 inversions. This number can be reduced significantly by combining the LapH method with a stochastic approach as described in Ref. [43] in detail.
For constructing all-to-all propagators with a stochastic approach, random vectors ρ are introduced which carry Dirac, time, colour and spatial indices. A propagator can be computed by solving for X r [b] . Here the random vectors are not chosen to be random valued in all elements, but diluted, which is used to reduce the stochastic noise, see Ref. [43] for details. Correspondingly, the composed index r[b] counts the total number of random vectors, N R , via r and the total number of dilution vectors, N D , via b. When constructing the all-to-all propagator via the zeros in the diluted random vectors ensure exact zeros in the product ρ r[b] ρ r[b] † , which reduces noise significantly. However, diluting random vectors leads to higher computational costs due to solving Eq. 7 for each of these vectors separately. It is expected that the noise in correlation functions built from diluted stochastic propagators reduces with 1 / √ N R and 1 /N D , which favours more dilution vectors over more random vectors. Aside from this, each quark line needs its own set of random vectors to avoid a bias in the correlation functions.
In the stochastic version of LapH, denoted by sLapH, the random vectors are introduced in time, Dirac and Laph sub-space indices. A quark line can then be estimated stochastically via In our implementation we choose Z 2 random numbers. As mentioned before each quark line as defined in Eq. 9, needs its own random vector to avoid a bias. This means that at least four random vectors are needed to be able to compute the correlation functions relevant for ππ scattering processes. However, we use five random vectors, since the fifth random vector will allow for additional permutations in the four point functions. This improves the signal-to-noise ratio by a factor of 2 by only increasing the computing costs for inversions by 25%.
As dilution scheme we have chosen full dilution in Dirac space combined with a block dilution in time and an interlace dilution in the LapH sub-space. The number of dilution vectors are summarised in table 2.1. They were chosen in such a way that the number of inversions and the noise in our observables are minimised simultaneously.
We remark here that with the sLapH smearing scheme as explained above, only smeared-smeared correlation functions can be computed. Hence, we cannot compute matrix elements of local operators needed for instance for f π without major additional effort.  Table 2: Summary of the number of dilution vectors, N D , used in each index. We use a block scheme in time and an interlace scheme in eigenvector space. The total number of inversions is the number of random vectors, here 5, multiplied by the number of inversions for one quark line.

Operators
The phase shifts are extracted from the finite-volume energies as described in the next section. In order to estimate the finite-volume energy spectrum, we build a matrix of correlators with a set of operators that resemble the ππ isospin-2 system, and then use the variational method [45,4] to extract the spectral information. The SO(3) symmetry in continuum space is reduced to the octahedral group O on the lattice. For non-zero momentum P the symmetry is further reduced to the little group LG( P ) that leaves P invariant. We construct the operators with definite total momentum P in each irreducible representation (irrep) of LG( P ) via: Here, π + ( p, t) = x e i p· xd ( x, t)γ 5 u( x, t), is the single pion operator projected onto momentum p. With periodic boundary conditions, p is quantised as p = 2π L n, where n is a vector of integers. Λ is an irrep of the group LG( P ) and λ is the irrep row. { p 1,2 } * represents the set of vectors {R p 1,2 , R ∈ O}. C( P , Λ, λ; p 1 , p 2 ) are the Clebsch-Gordan coefficient for Λ 1 ⊗ Λ 2 → Λ, where Λ 1 (Λ 2 ) is the irrep that π + ( p 1 )(π + ( p 2 )) resides in. The two pions are separated by one time slice in order to avoid the complications due to Fierz rearrangement [46].
The correlation matrix is computed via: The variational basis contains various combinations of p 1 and p 2 that are allowed by the decomposition Λ 1 ⊗ Λ 2 → Λ. For the case of P = [0, 0, 0] we include for the A 1 irrep the operators with | p 1,2 | = (0, 1, 2, 3, 4), which gives a 5-dimensional correlation matrix. More details about the irreps and other P values can be found in Ref. [34]. The energies are obtained by solving the generalised eigenvalue problem: It can be shown that the eigenvalues λ i (t) behave like where E i is the i-th eigenvalue of the Hamiltonian of the system. However, we focus on zero total momentum where apparently solving the GEVP does not give any advantage over using p 1 = p 2 = 0 only. This is due to a rather weak coupling of different momenta in the matrix of correlators. All results presented in the following are, therefore, obtained directly from the four-point function at p 1 = p 2 = 0.

Removing Thermal States
As discussed in Ref. [31] and references therein, the spectral analysis of the two pion correlation function in the case of total zero momentum deviates from the usual cosh like behaviour. It was shown in Ref. [31] that the diagonal elements of the correlation matrix, Eq. 11, for the two particle system obey the following spectral decomposition in the limit of large Euclidean times with E ππ the two pion energy and M π the single pion mass, respectively, and constants A 0 and c. This spectral decomposition differs from the standard by the term constant in Euclidean time. In the thermodynamic limit T → ∞ this polluting term vanishes, but for finite T it will dominate the correlation function at large Euclidean times. To remove this pollution, it was proposed to take a finite difference first [47] and then build the following ratio [31] with C π (t) the single pion two-point correlation function. One can show that the ratio has the functional form [31] with t = t + 1/2 − T /2 and δE = E ππ − 2M π the energy shift. The generalisation of this procedure to correlation matrices and the variational method persists in shifting the correlation matrix Eq. 11 before using the standard GEVP procedure, see Ref. [34] for details.

Finite Volume Methodology
We are interested in the limit of small scattering momenta for the ππ system with I = 2 below inelastic threshold. Using the finite range expansion, the scattering length a 0 and the effective range r 0 can be related to the energy shift δE by an expansion in 1/L as follows [3] with δE = E ππ − 2M π and coefficients [3,48] More generally, including also non-zero total momentum, Lüscher's method relates the phase shifts δ to the finite volume energy shift via the relation where the matrix elements of the matrix M are given as [6] M Z is Lüscher's generalised Z-function and the lattice scattering momentum. The elements of the tensor C lm,js,l m can be found in Ref. [6]. For the case = 0 with total zero momentum and no mixing with higher partial waves, Eq. 18 reduces to The scattering momentum q is then given as from which its lattice version, Eq. 20, can be computed.
As proposed in Ref. [6] the continuum dispersion relation can be replaced by a lattice modified one. However, we do not see any difference in using one or the other version of the dispersion relation for the case of zero total momentum studied here. Hence, we stick to the continuum dispersion relation for this paper.

Results
Before coming to the actual results, let us first describe our analysis procedure. Statistical errors are always computed using a bootstrap procedure with R = 1500 bootstrap samples. The bootstrap analysis is chained such that also statistical errors for best fit parameters can be determined. The gauge configurations are sufficiently separated HMC trajectories that autocorrelation does not play a role here, as we explicitly checked using the blocked bootstrap procedure. Energy shifts are determined from fully correlated fits of the ratios, Eq. 15 to the data. The single pion energy levels needed for the fit are determined from the corresponding two point functions using a one state exponential (cosh) fit to the data, again fully correlated. The correlation between the single pion energy level and the ratio data is taken into account by using the same bootstrap samples.
In figure 1 we show example plots of the ratio for various ensembles. For a representative choice of the fit range we also show the corresponding best fit curves obtained by fitting Eq. 16 to the data. For the fit ranges shown one observes visually a good agreement between fitted curve and data.
The fits to the ratios and the two point functions are repeated for a large number of fit ranges. We then assign a weight to every of these fits and quantities X = E π , E ππ , where p X is the p-value of the corresponding fit and ∆ X the statistical error of X determined from the bootstrap procedure.
Using the weights w X , we compute the weighted median over all fit results on the original data to obtain our estimate for the expectation value X . The 68.54% confidence interval of the weighted distribution provides an (not necessarily symmetric) estimate for the systematic uncertainty stemming from the different fit ranges. The statistical error on X is computed from the R bootstrap samples of the weighted median.
When derived quantities like q cot δ are being determined, we follow the same procedure, just that the weights are now given by the products of the weights of the different contributing energy levels.
The fit ranges are chosen such that for both the two point function and the ratio at least five time slices are included in the fit. The two point function is always fitted in an interval [t π 1 , T /2] with t π 1 > 6. For the ratio we fitted in the interval [t ππ 1 , t ππ 2 ] with t ππ 1 ∈ {11.5, 13.5, 15.5} and t ππ 2 ∈ {22.5, 21.5, 20.5, 19.5, 18.5} for L = 24, 20 and t ππ 2 ∈ {30.5, 29.5, 28.5, 26.5, 24.5} for L = 32. We note that due to the weighting procedure described above we could have also included smaller values for t π 1 and t ππ 1 without affecting the final result.

Systematic Effects
One of the main issues in this investigation is the question of systematic uncertainties in our analysis procedure. In particular, we have to consider possible contributions by the neutral pion which is much lighter than the charged pion on our ensembles due to twisted mass isospin breaking effects. Possible effects on the I = 2 pion scattering length extraction have been discussed in detail in Ref. [31].
In Ref. [31] the authors were not able to find evidence of neutral pion contributions within their errors. Even though we have significantly smaller statistical errors, we still do not have any evidence for these effects within the statistical uncertainties. In the pseudo-scalar two point function the possible effect is an excited state with mass M π + M π 0 due to the neutral pion having the vacuum quantum numbers in twisted mass lattice QCD [49]. This excited state could, however, never be identified. In the fourpoint function there can be effects from either close-by excited states at small Euclidean times or thermal states at large Euclidean times. Again, we do not see any evidence of those effects in our data.
However, the analysis procedure detailed above is designed such that such effects should be covered by the systematic error we determine from the weighted distribution.  Table 3: δE, q cot δ 0 and M π a 0 computed with total zero momentum.
This systematic error mirrors possible deviations from the theoretically expected curve and contributions by excited states or pollutions at large Euclidean times.
The neutral pion might also contribute to the exponential finite size corrections in our data. As discussed below, for M π and f π we take these twisted mass specific effects into account as determined in Ref. [50] from the data. For q cot(δ 0 ) finite size corrections specific to twisted mass have not been computed in χPT and, hence, we can only include the corrections computed in continuum χPT in Ref. [51]. However, these finite size corrections computed in continuum χPT have little influence on our results. Hence, we do not expect large effects from twisted mass specific finite size effects.
This conclusion is further supported by the fact that we do not observe large discretisation errors in the results for M π a 0 . Any contribution from the neutral pion should show up as a O(a 2 ) lattice artefact. Of course, all these indications are still not enough to finally exclude such systematic effects in our results. But we conclude that within our uncertainties they are negligible.

Lüscher formula to O(1/L 5 )
We will discuss several procedures to determine scattering parameters from the data. The first of which is to consider Eq. 17 to the order 1/L 5 and, hence, neglecting the contribution from the effective range, a 0 can be determined from δE, M π and L by numerically solving Eq. 17 for a 0 . The corresponding results for a 0 in units of M π can be found in table 3.
As an illustration of our analysis procedure we show example histograms in figure 2. With plain histograms we mean histograms of the unweighted results for different fit ranges. In the weighted histograms the weights have been applied according to Eq. 23. And with bootstrap data we denote the median of the weighted distribution evaluated on the bootstrap samples. In the histograms we plot the densities of the distribution.
In the upper left panel the plain and the weighted histogram of M π a 0 determined from A ensembles B ensembles D45 Figure 3: We show finite size corrected data for M π a 0 as a function of M π /f π with filled symbols compared to uncorrected data with open symbols. The corrected data include also the systematic uncertainty in order to allow for a comparison.
Eq. 17 on the A40.24 ensemble are compared. The weighting leads to a well defined peak in the histogram, which is representative for the findings on most ensembles. In the upper right panel a comparison of the weighted histogram and the histogram of the bootstrap samples of M π a 0 again for A40.24 is shown. That the distribution of the bootstrap samples is approximately normal can be inferred from the lower left panel, where we show the QQ-plot of the bootstrap sample quantiles versus the theoretical quantiles of a standard normal distribution Nμ ,∆µ . Hereμ is the estimate of M π a 0 and ∆µ its statistical error determined from the standard deviation over the bootstrap samples. This comparison indicates that the systematic and statistical uncertainties are approximately of the same size. This finding is again representative for most of the ensembles.
In the lower right panel of figure 2 we show again the plain versus the weighted histogram of M π a 0 , but this time for ensemble D45. This ensemble shows the largest systematic uncertainty of all the ensembles investigated here, which is also asymmetric. The weighted median and the systematic uncertainty which are indicated by the circle and the horizontal error bar above the histogram, respectively, show this asymmetry clearly, even if it is not easy to identify by eye.
For the successive analysis we also need to consider (exponentially suppressed) finite size corrections to our data. For M π /f π and M π we use the results of Ref. [50] to correct our data. The corresponding data for M π /f π and the finite size correction factors are summarised in table 7 in the appendix. We remark that A40.20 and D45. 32 have not been considered in Ref. [50]. For D45.32 we use the factors computed in Ref. [50] for ensemble D20.48 with almost identical M π L-value. For A40.20 we do not need the finite size corrections in the subsequent analyses.
For a 0 we apply the asymptotic finite size correction formula Eq. (31) from Ref. [13] ∆(q cot δ) = (q cot δ) L − (q cot δ) L=∞ which is valid close to threshold (small q 2 ) only. The c(n) are multiplicities for the n and can be found for instance in Ref. [13]. Note that at zero scattering momentum lim q→0 q cot δ = 1 /a 0 holds which gives us directly the finite size correction for the scattering length. For comparison we show the bare data for M π a 0 as a function of M π /f π with open symbols and the corresponding finite-size corrected data with filled symbols in figure 3. From this plot it is visible that most of the finite-size effects in this analysis stem from the ratio M π /f π . This is because the finite size corrections for M π and f π are opposite in direction. One can also see the effect of including the systematic uncertainties in the errors, which are included for the finite size corrected data points, but not for the uncorrected ones. This allows one to get an impression of their size and for which ensembles they are actually relevant.

Lüscher Formula to O(1/L 6 )
We have three A40 ensembles available, which differ only in their spatial extends, L = 20, 24 and L = 32, respectively. Here, we can apply Eq. 17 in order to estimate the scattering length a 0 and the effective range r 0 from the L-dependence. It amounts to a fit to the data for δE for the three volumes with two fit parameters. For aM π we used the values from the largest volume ensemble A40.32. The result is summarised in the first column of table 4.
According to the χ 2 value (dof = 1) this is not a good fit. The data for δE is shown together with the best fit and an error band in the left panel of figure 4. The result for M π a 0 is lower than the one for the A40.32 ensemble alone (see table 3), but it deviates not more than two σ. The effective range parameter r 0 is determined by the fit only with large statistical uncertainties.
It is likely that in particular A40.20 suffers still from exponential finite size artefacts. Therefore, we repeat the analysis with only A40.32 and A40.24 included. Now the fit basically reproduces the result to O(1/L 5 ) of the A40.32 ensemble, and there is no 5.14 -0.79 -

Finite Range Expansion for q cot δ 0
Instead of using the expansion in 1/L as defined in Eq. 17, one may also determine q cot δ 0 close to threshold directly and perform the effective range expansion as follows: where δ 0 is the s-wave phase shift, a 0 the corresponding scattering length, r 0 the finite range parameter and q the scattering momentum. Using the three ensembles A40.32, A40.24 and A40.20 again, we are able to obtain q cot δ 0 for three values of the squared scattering momentum q 2 . We then apply finite size corrections using Eq. 24.
The corresponding data are shown in the left panel of figure 5 together with the best fit to expression Eq. 25. The statistical error of the fit is indicated by the grey band. In addition, the value extrapolated to q 2 = 0 is shown.
The best fit result is summarised in the third column of table 4. In contrast to the fit of the 1/L expansion, the finite range expansion provides a good description of q cot δ 0 . The χ 2 -value indicates a good fit, bearing in mind that there is little freedom left. The statistical uncertainties are smaller than the ones obtained from the 1/L expansion. In particular, a statistically significant value for the effective range parameter r 0 is obtained.
Like in the previous case, where we used Lüscher's direct method, we also performed an extraction of a 0 and r 0 with only the largest two volumes L = 24, 32 with the effective range formula. The result is summarised in the last column of table 4 and the data are shown in the right panel of figure 5. While the scattering length and effective range are in good agreement with previous extractions, we lose statistical significance for r 0 again. However, the error on a 0 is only slightly bigger than found for the other three fits. Note that the finite size corrections to q cot δ 0 do not change the fit result significantly.

Chiral Extrapolation
Because we have determined the scattering parameters for bare quark masses corresponding to larger than physical quark masses, we need to extrapolate to the physical point. For the I = 2 ππ scattering length χPT is particularly well suited, because it depends on only one low energy constant at NLO. As suggested in Refs. [52,53], it is convenient to write the χPT expression for M π a 0 as a function of M π /f π , because all quantities are dimensionless and no scale input is needed. The NLO χPT expression for M π a 0 written as a function of (M π /f π ) at a χPT renormalisation scale µ R = f π,phys reads [52,53] with ππ related to the Gasser-Leutwyler coefficients¯ i as follows [54] ππ (µ R ) = Moreover, one can show in Twisted mass χPT that the leading lattice artefacts to M π a 0 are of order O(a 2 M 2 π ) [51]. At NLO we, hence, consistently describe our data with the continuum χPT formula provided above. Since we are not (yet) able to determine the pion decay constant f π with the sLapH approach, we use data presented in Ref. [55] for the ratio M π /f π . The values for M π /f π for all ensembles are compiled in table 7 in the appendix. In this table we also give the values for the finite size correction factors as determined in Ref. [50]. The finite size corrections to a 0 are analytically computed using Eq. 24 by setting q cot δ 0 = 1 /a 0 . For the finite size corrections of M π in M π a 0 we also use the factors compiled in table 7 as discussed above.
The ratio M π /f π can be determined with significantly smaller uncertainty than M π a 0 . Therefore, we do not expect that the missing statistical correlation between M π /f π and M π a 0 plays any role in the following analysis.
For the fit we propagate the errors on M π /f π and the finite size corrections on M π and f π using resampling. For the error on M π a 0 we add the statistical and systematic uncertainty of M π a 0 and the statistical uncertainty of the finite size correction factor K Mπ in quadrature. In the χ 2 minimisation we take errors both on M π a 0 and M π /f π into account. We only include the large volume A40.32 ensemble (and not A40. 24 and A40.20) in the fit.
In the left panel of figure 6 we show the fit to all the data (i.e. with no cut in M π /f π ). The solid line represents the best NLO χPT fit to our data, the dashed line the LO

NLO-LO ChiPT
A ensembles B ensembles D45 extrapolated Figure 6: Left: M π a 0 as a function of M π /f π determined from Eq. 17 to O(L −5 ) with FS corrections for M π /f π and M π a 0 . The A-ensembles do not include A40. 24 and A40.20. We add the LO χPT prediction as the dashed line. The best fit to the data by the NLO χPT expression is shown as the solid line with error band. M π a 0 | phys is plotted as the triangle. Right: the same but with LO χPT (M π a 0 ) LO subtracted.
parameter-free χPT prediction for M π a 0 . One observes that the LO χPT prediction already describes the data surprisingly well, and the NLO fit makes only a small correction to the LO curve. We also show the value of M π a 0 extrapolated to the physical point. In the right panel of figure 6 we show the same, but with the LO χPT prediction (M π a 0 ) LO subtracted. The data for the A-and B-ensembles fall to a good approximation on a single curve, whereas the only D-ensemble deviates slightly. While not inconsistent with a statistical fluctuation, this might be due to the rather small physical volume of the D45.32 ensemble. Another possible reason might be a 2 M 2 π lattice artefacts together with higher order terms in continuum χPT.
We can explore this by restricting the fit range by applying an upper cut in the M π /f π values. For the case of the fit including only data points with M π /f π < 2.4 this is shown in figure 7, which is otherwise identical to figure 6. The cut value is indicated with brackets in the plots.
The results of the fits to our data with different fit ranges are summarised in table 5. We observe that different fit ranges do have little effect on the extrapolated value of M π a 0 | phys . However, the fit quality improves with decreasing upper fit range, with the best fit for the cuts M π /f π < 2.4 and M π /f π < 2.5. The only fit parameter ππ shows a ππ M π a 0 | phys χ 2 dof p-value M π /f π -cut 5. 13

NLO-LO ChiPT
A ensembles B ensembles D45 extrapolated Figure 7: Like figure 6, but for the fit including only data points with M π /f π < 2.4 as indicated by the brackets.
large variation, though the fitted values are mostly consistent within errors. The reason for this large variation is visualised in the right panels in figures 6 and 7. The curvature is mainly driven by the points around M π /f π = 2. It is also visible that the fitted curves with cut still describe the data reasonably well within errors also for M π /f π values larger than the applied cut.
As a final result we quote the weighted average of the two fits with cuts M π /f π < 2.4 and M π /f π < 2.5, respectively with the systematic error estimated from the maximal deviation of the single results in table 5 to the finally quoted one. The uncertainty stemming from the different fit-ranges of the fit to the ratio data turns out to be sub-leading.

Discussion
The results presented in the previous section require some discussion. The first important point is that contributions from thermal and excited states as well as the contributions from the neutral pion can be excluded within the statistical uncertainty. With the use of the ratio, Eq. 16, thermal states constant in time are minimised if not cancelled completely. In section 4.2, e.g. figure 2, we show histograms created from fits to a large number of fit ranges. Even for D45 where we see the largest systematic uncertainties stemming from the different fit ranges, the systematic uncertainty is not much bigger than the statistical uncertainty. This lets us conclude that we cannot resolve any of the aforementioned systematic uncertainties within our current statistical precision. Thanks to the three ensembles A40.32, A40.24 and A40.20 we were able to extract scattering parameters using three different methods: i) 1/L expansion, Eq. 17, to the order 1/L 5 for every A40 ensemble separately, ii) the same expansion but to the order 1/L 6 and iii) a direct fit of the effective range expansion to the values of q cot δ 0 (q 2 ). Comparing the results in the sections 4.2, 4.3 and 4.4 shows that all three methods give compatible results for M π a 0 , apart from method i) for A40.20. The reason for the deviation of method i) for A40.20 can be twofold: on the one hand M π L = 2.99 might be too small and exponentially suppressed finite volume effects cannot be ignored. On the other hand the physical volume might be too small so higher orders in 1/L in Lüscher's formula and higher orders in q 2 in the effective range expansion are needed. Excluding A40.20 in methods ii) and iii) leads to smaller M π a 0 values, though still compatible within errors to the case when including A40.20. These smaller values, however, are in better agreement to method i) for both A40.32 and A40.24. From this fact alone it is hard to deduce which kind of volume effect we are dealing with here. That the physical volume of A40.20 might be too small is supported by the fact that the D45 ensemble, which is the other ensemble with small physical volume, gives a lower M π a 0 than expected although M π L = 3.87 is in the range of the other ensembles. The important consequence of this finding is that very likely all other ensembles do not suffer from significant finite volume effects. On every other ensemble than A40.20 and D45 the M π L-value and the physical volume is at least as big as on A40.24 which seems to be free from volume effects within our current precision.
For methods ii) and iii) also the effective range r 0 can be extracted from fits to the data. The best fit values vary significantly and have mostly very large uncertainties. Only for one fit (see table 4) the r 0 -value is significant. We conclude that at the current level of available data and statistical precision we are not able to reliably extract the effective range parameter.
The chiral extrapolation we present in section 4.5 includes three values of the lattice spacing. The A-and B-ensembles agree quite well, only at larger pion mass values small deviations show up. Unfortunately, we have currently only one ensemble for the smallest lattice spacing value available, namely D45.32. D45.32 has a quite small physical volume, while M π L ≈ 3.8 which is compatible to the other ensembles. Currently, we cannot finally investigate the reason for this discrepancy. It could be a statistical fluctuation, which is supported by the exceptionally large systematic uncertainty for this ensemble,   [30], ETM (2013) [31], this work denoted as ETM (2015), Yagi et al. [58], Fu [59] and PACS-CS [60].
or a finite volume effect, as discussed above. However, we can also not exclude a lattice artifact. We are working on additional D-ensembles, which will allow us to investigate this further. In table 6 we provide a compilation of theoretical determinations for M π a 0 and ππ from the literature: this includes the LO χPT prediction as well as the value determined using χPT and Roy equations from Ref. [56] denoted as CGL. For the lattice results we have decided to include only direct determinations for which the chiral extrapolation has been performed: CP-PACS [57], NPLQCD (2006) [52], NPLQCD (2008) [30], ETM (2013) [31], this work denoted as ETM (2015), Yagi et al. [58], Fu [59] and PACS-CS [60]. We quote statistical and -where available -systematic uncertainties separately. For NPLQCD (2008) there is only the combined statistical and systematic uncertainty.
The CP-PACS and the here presented ETM (2015) results are based on three values of the lattice spacing, two have been use by the authors of ETM (2013) and by Fu. The others have used only one value of the lattice spacing, however, NPLQCD (2008) and PACS-CS have used χPT to estimate discretisation effects. Finite size effects are estimated in the works of NPLQCD (2008), ETM (2013), Yagi et al. and Fu using chiral perturbation theory. In addition to the estimate from chiral perturbation theory we have studied in this work several volumes to obtain an estimate of residual finite volume effects.
As can be observed visually in figure 8, the values for M π a 0 show an overall good agreement. In particular, there is no definite dependence on N f . The exception are the high values of NPLQCD (2006) and PACS-CS. PACS-CS quotes in addition a small uncertainty, which makes the deviation to our result statistically significant. If the systematic uncertainty is added to the statistical error, the results are compatible again. PACS-CS has studied the smallest pion mass value around 170 MeV of all the works considered in this comparison. However, they claim they cannot fit this point with NLO χPT. Therefore, they had to include higher order terms. More results at close to physical pion mass values might, therefore, clarify whether this is a statistical or systematic effect in the data, or whether the chiral extrapolation with NLO χPT is misleading. The values for ππ show a large variation, which is, however, covered by the large errors on this LEC.

Summary
In this paper, low-energy ππ scattering in the isospin I = 2 channel is studied within Lüscher's finite volume formalism in lattice QCD. We use for the first time N f = 2+1+1 dynamical quark flavours based on gauge field configurations provided by ETMC. The list of ensembles covers three values of the lattice spacing, several volumes and a large range of pion mass values, see table 1. We determine energy shifts δE = E ππ − 2M π using the stochastic LapH method and convert them in scattering length values applying Lüscher's formalism.
We apply different, though closely related, methods to determine scattering parameters and find compatible results for M π a 0 . The finite range parameter r 0 cannot be determined with sufficient certainty.
Due to the use of ensembles with a broad range of pion masses we are able to extrapolate M π a 0 towards the physical pion mass point. After correcting for finite volume errors of the scattering length a 0 we are able to make a rather smooth extrapolation towards the physical pion mass value, the result of which is shown in figure 6. We do not observe significant discretisation effects between A and B-ensembles. The only D-ensemble D45 shows a deviation which can be equally well explained by discretisation effects, finite volume effects or by a statistical fluctuation. To clarify this point, we are generating data on more D ensembles with smaller pion masses.
Our final result for a I=2 0 is, M π a 0 = −0.0442(2) stat ( +4 −0 ) sys . We have compared our result to a list of other lattice determinations and χPT combined with Roy equations and found mostly remarkably good agreement.
The currently ongoing extension of the study presented in this paper is ππ scattering with I = 1, where the ρ resonance is present. With the perambulators ready from the distillation process, this is straightforward to do and first results are available. For the pion-pion scattering, the channel I = 0 is more complicated, particularly for the twisted mass formulation due to isospin breaking at finite lattice spacing values. Another currently ongoing extension is to address low-energy scattering of other mesons: e.g. πK and KK scattering and more than two mesons. Such techniques can also be extended to charmed meson scattering processes relevant for the recently discovered XYZ states.