Extraction of isoscalar $\pi\pi$ phase-shifts from lattice QCD

We conduct a two-flavor ($N_f=2$) lattice QCD calculation of the elastic phase-shifts for pion-pion scattering in the scalar, isoscalar channel (the $\sigma$-meson). The calculation is performed for two quark masses corresponding to a pion mass of $315\text{ MeV}$ and $227\text{ MeV}$. The $\sigma$-meson parameters are extracted using various parametrizations of the scattering amplitude. The results obtained from a chiral unitary parametrization are extrapolated to the physical point and read $M_\sigma = (440^{+10}_{-16}(50) - i\,240(20)(25))\text{ MeV}$, where the uncertainties in the parentheses denote the stochastic and systematic ones. The behavior of the $\sigma$-meson parameters with increasing pion mass is discussed as well.


I. INTRODUCTION
The lightest excited state in the spectrum of hadrons is at the same time one of the most controversial. As described in detail in an extensive review [1], its properties (mass and width) and even presence were debated for a long time. Many precise analyses finally led to the currently accepted ranges for mass, 400−550 MeV, and width, 400 − 700 MeV, of the so-called σ or f 0 (500) resonance [2]. This scalar/isocalar excited state has a dominant decay channel to two pions, with a very uncommon shape of the partial wave in this channel.
Lattice QCD is the only method to compute the hadron properties directly in terms of quark-gluon QCD dynamics. In the context of ππ scattering many important results have been reported in the I = L = 1 channel [3][4][5][6][7][8][9][10][11][12][13]. Due to the presence of disconnected diagrams the corresponding calculations in the I = L = 0 channel were not possible for a long time, despite early pioneering works [14][15][16]. The first results have been reported recently by the Hadron Spectrum Collaboration [17], extracting also phase-shifts using the Lüscher framework [18] in combination with moving frames [19]. The data have then been analyzed and extrapolated to the physical point using the Inverse Amplitude Method in Ref. [20]. The isoscalar scattering length at three pion masses has been extracted recently by the European Twisted Mass collaboration [21].
In the present work we report new N f = 2 lattice QCD results for two different pion masses (M π = 227 * guodehua@gwmail.gwu.edu † aalexan@gwu.edu ‡ ramope@if.usp.br § maximmai@gwu.edu ¶ doring@gwu.edu and M π = 315 MeV) and analyze them. To extract a robust energy spectrum in a specific scattering channel, a large and proper interpolating field basis is required.
To have an interpolating field basis that has the correct quantum number and symmetry as the scattering channel and has enough overlap with the relevant eigenstates of the system, both quark-antiquark (qq operators) and two-hadron interpolators are included. We perform a variational analysis [22] with this interpolating field basis to extract several low lying energy states that are in the elastic scattering region. The evaluation of the correlation functions are carried out using the Laplacian-Heaviside smearing method [23,24] for all the possible quark diagrams calculated from the Wick-contraction procedure. To have more phase-shift data points in different kinematic region so as to better describe the energy dependence of the ππ scattering phase-shifts, we implement our calculation in three boxes with different elongation factor and two total-momentum frames, one being the rest frame P = (0, 0, 0), the other one moving along the elongated direction with P = (0, 0, 1), the smallest non-zero momentum allowed by the boundary conditions. Since we are interested in two-particle scattering in the elastic region, the physical observables such as the phase-shift can be obtained from the energy spectrum in finite volume using the Lüscher formula [18].
In the second step of the present work we perform an energy-dependent analysis using directly the set of energy eigenvalues at different elongations and momenta. This allows to take into account the correlations between different energy eigenvalues. To this end, we formulate the scattering amplitude in a similar manner to the K-matrix approach, which allows to access finite-volume energy eigenlevels (as positions of the poles of this amplitude). The free parameters of the parameterizations are fitted to reproduce the lattice data and then used in the infinitevolume formulation to determine the parameters of the arXiv:1803.02897v1 [hep-lat] 7 Mar 2018 isoscalar ππ resonance. We use different parametrizations to estimate the systematic uncertainty from this source. Specifically, we consider a general expansion in an energyvariable conformally mapping the energy plane to the unit disk, similar to the analysis of Refs. [25,26]. Second, we employ a model based on the chiral unitary approach (UChPT), used, e.g., in Refs. [9,27,28]. Subsequently, the UChPT amplitude is extrapolated to the physical point. Our final result, based on all lattice data presented here with and without the isovector channel data [9], reads M phys σ = (440 +10 −16 − i 240 +20 −20 ) MeV and agrees with the result of the most recent analysis of experimental data [1] within the quoted 1σ region. Additionally, the study of the pion mass dependence of resonance mass and coupling to the ππ channel is studied in a broader range of M π , as it was done in Refs. [20,29].
The paper is organized as following: In Sec. II we describe the details of the lattice calculation to obtain the energy eigenvalues and correlation matrices; In Sec. III we describe the scattering amplitudes used for the extrapolation of the lattice results in energy and pion mass; The results of these analyses are discussed in Sec. III C, and the overall summary is given in Sec. IV.

A. Interpolating basis
As we mention in the introduction, the composition of the σ-meson might include different kinds of components. In order to extract low-lying energy states from the correlation function, we perform a variational analysis. In this study, we are interested in two-pion elastic scattering. To better extract the low lying energy levels in the elastic scattering region, we choose a set of interpolating fields with the same quantum number including the quarkantiquark (qq) and meson-meson interpolating fields in the variational basis. There are several reasons for using a large variational basis. Firstly, it helps resolve energy states that are nearly degenerated. Secondly, it offers large enough overlap with the eigenstates of the Hamiltonian, which can improve the accuracy and stability of the extracted energy spectrum.
The correlation matrix is constructed from two-point functions of all combinations of the interpolating fields in the variational basis. The elements of the correlation matrix are We denote the interpolators as O i with i = 1, ..., N with N being the number of interpolating fields in the basis. The eigenvalues of the correlation matrix can be obtained by solving the generalized eigenvalue problem for a particular initial time t 0 and for each time slice t. The energies of the system are then determined from the long-time behavior of the eigenvalues [30] where the correction term depends on the energy difference ∆E n = E N +1 − E n . According to this behavior, for the low-lying energy states, the larger interpolating basis we use, the faster the correction vanishes. However, the benefit from enlarging the interpolating basis decreases because the energy eigenstates get denser in the higherenergy part of the spectrum. Our goal is to choose a interpolating basis that is good enough to capture the energy eigenstates in the elastic scattering energy region.
In this work, we consider lattices with one spatial direction elongated. The corresponding rotational symmetry group for this elongated box is D 4h which is a subgroup of the full rotational symmetry group SO (3). Therefore, the angular momenta that label the irreducible representations (irreps) of SO(3) are split into multiplets related to irreps of the D 4h group. The resulting split for the lowest angular momentum multiplets is listed in Table I. The σ-meson has angular momentum L = 0 and positive parity. The irrep L = 0 maps to the one-dimensional A + 1 irrep. The lowest state in this channel corresponds to a π-π state for which the two pions decay at rest. Therefore, the lowest state changes very little when varying the elongation of the box. The excited states energy changes as we increase the elongation and we can scan this corresponding energy region.
In order to have more energy values data points in different kinematic region, we implement the boosted frame method. The idea is to boost the whole system with a given momentum P in a certain direction. Due to the relativistic effects, the box along the boosted direction is contracted [19]. A generic boost direction of the box will further reduce the original symmetry group to a subgroup which depends on the direction of the boost. In this study, we implement a boost to an elongated box parallel to the elongated direction. In this case, the boost reduces the rotational symmetry group of the cubic box from O h to D 4h , but it does not change the rotational symmetry group for the elongated box, which is still D 4h .
We note that the states in the A + 1 irrep belong to different irreps of SO (3). The A + 1 irrep couples not only to L = 0, but also to other higher angular momentum channels such as L = 4 and so on in the O h group for the cubic box and L = 2, L = 4 and so on in the D 4h group for the elongated box. However, to study the σ-meson, we are interested in a low energy region where the two-pion states with relatively small scattering momentum. In this case, the effect from the L ≥ 2 channels are small and their contribution can be safely neglected because it is kinematically suppressed through the barrier factor. Indeed, the influence of the D-wave in the extraction of the S-wave from energy eigenvalues in A + 1 has been estimated in Ref. [31] using realistic S and D waves; it was found to be small. As a result, we focus on the states in the A + 1 irrep. For the volumes considered in this study, there are only two or three low-lying energy states in the elastic scattering energy region. As mentioned before, we need a basis that has overlaps both with the resonance state (to take into account possible two-quark components of the σ-resonance) and with the states that have a dominant two-pion content. We include four quark-antiquark interpolators in our basis so that we can see their effects on the energy spectrum. These four quark-antiquark interpolators have the form The u(t) and d(t) denote the up and down quark on the entire t time slice which is a column vector of size Their form in the creation operators is defined as Γ i (p) which can be derived using The details for Γ i (p) and Γ i (p) are listed in the first four row of Table II. The first interpolator is point-like and the next three interpolators involveqq pairs that are separated by several lattice spacing, defined using the covariant derivative The forth interpolator has a different gamma matrix structure.
In previous studies for the ρ-meson resonance in the π-π scattering channel [4,9,12], it was shown that the quark-antiquark interpolators are not sufficient to extract a reliable spectrum in the interacting theory where the actual eigenstates are mixed with quark-antiquark basis states and multi-hadron basis states. The reason is that the quark-antiquark interpolators have little overlap with the multi-hadron states and the overlap is shown to be suppressed by a power of the lattice volume [12]. To solve this problem, we include the pion-pion interpolators in the variational basis. First, we construct the pion-pion interpolators to have isospin I = 0 and I 3 = 0 which correspond to the σ-meson: i where π + , π − and π 0 are given by and As we mention before, in this study, we focus on the A + 1 irrep which mainly couples to S-wave of the pionpion scattering. To construct interpolators transforming according to the A + 1 representation, we project the general ππ interpolators into A + 1 using where R(g) implements the rotation associated with the symmetry transformation g, and χ A + 1 is the character of group element g in the A + 1 irrep. In order to have more energy levels in different kinematic regions, we implement two different total momenta for the system. The ππ operators in the rest frame P 0 = (0, 0, 0) in A + 1 irrep are as follows: In the boost frame P 1 = (0, 0, 1), the lowest scattering momentum p 1 = (0, 0, 1). Therefore, we use the following interpolators ππ (1) where P = {(0, 1, 1), (1, 0, 1), (−1, 0, 1), (0, −1, 1)} are all the possible momenta generated by symmetry transformations R(g)p from p = (0, 1, 1). The value of p 2 in the equations above is imposed by momentum conservation. The reason we only apply the positive total momentum P 1 in the boosted case is that the expectation values for the correlation functions associated with momentum P 1 and −P 1 are the same due to rotational symmetry.
With these six interpolators, we construct a 6 × 6 variational basis. Each entry in the correlation matrix can be calculated through the Wick contraction procedure. There are three types of entries in the correlation matrix, The notation above is defined as where j k+1 is defined to be j 1 and M −1 (t, t ) = u(t)ū(t ) F is the quark propagator between time slices t and t (for more details about the notation see Ref. [9]).
The calculation of the correlation functions involves the evaluation of the all-to-all propagator which is not practical to compute directly from the position space. We evaluate them in another way using the LapH method. The setup details for the LapH method in this study can be founded in our previous study [9]. To calculate the all-to-all propagator projected on the LapH space requires a large number of fermionic matrix inversions. This was done efficiently using our GPU inverters [32].

B. Finite-volume spectrum
In the I = 0 channel, each entry of the correlation function contains temporally disconnected diagrams in which the trace only has one time slice. Evaluation of these diagrams is difficult. The reason it that this channel has the same quantum numbers as the vacuum and the correlators do not vanish as the time separation is taken to infinity. The constant contribution has to be subtracted in order to get the expected exponential behavior. In our study, we implement two approaches to subtract the vacuum contribution from temporally disconnected diagrams. In the first method we estimate the vacuum contribution by taking the average of the vacuum bubble and subtract this value from the original correlation functions as The second approach to solve this problem is to consider the so called shifted correlators instead of the original correlators in the correlation matrix where d is the time shift between two correlation functions.
In this case, since the vacuum contribution is a constant in the correlator, it is subtracted implicitly when taking the difference of correlation functions. We compare the energy levels extracted from these two approaches at M π = 227 MeV with different elongation factor η in Fig. 1.  The parameters for the ensembles used in this study. The lattice spacing a for each ensemble is listed as well as the number of gauge configurations am N , af π , and af K represent the nucleon mass, pion decay constant and kaon decay constant in lattice units. The two errors for the lattice spacing are stochastic, from the w 0 /a determination, and a systematic one estimated to be 2%.
We note that the results from direct subtraction have smaller error bars, but they seem to be inconsistent with expectation from the UChPT predictions when using the parameters for this model extracted from the ρ study. Furthermore, these results are also inconsistent with the ones extracted using the moving frame correlators that do not require a vacuum subtraction. It turns out that energy levels extracted using this method are very sensitive to the value of the constant used to subtract the correlator, the vacuum expectation value of the one-point functions generated by the interpolators. One possibility is that the values we computed are biased by wrap-around effects in the time direction. The shifted correlator method does not suffer from this problem, but it generates results with larger error-bars. The reason for this is that, for the same fitting window as in the direct subtraction method, the correlators involved are noisier since the shifted values of the correlators correspond to later times. We decided to use the shifted correlator method for the zero momentum states. For the moving states we do not need any subtraction, since the one-point functions vanish in this case.
We also looked at the stability of our results with respect to varying the interpolator basis. Our conclusions were similar to the ones derived in our ρ study [9]. The only noticeable difference is that for the moving states, we found that the energy levels are more sensitive to the presence of the O 4 = γ i ∇ i interpolator in our basis. We believe this is because the otherqq operators have the same γ-matrix structure, all γ = 1, and thus O 4 provides a significantly different overlap with the relevant states.
The energy levels extracted for this channel both for zero-momentum and moving states and the details of the fitting parameters are listed in Appendix A.

C. Phase-shift formulas
As mentioned in the introduction, in this study, we only consider two-pion scattering below the inelastic threshold. To connect the two-hadron state energies determined from lattice QCD with the physical observables, i.e., the phaseshifts in the continuum, we use Lüscher's formula [18] and its extensions to the elongated box [33] and to states boosted along the elongated direction [34]. Note that the extraction of resonance parameters from the phase-shifts requires a more delicate analysis in case of σ-resonance. This is due to the fact that the σ cannot be described by an ordinary Breit-Wigner shape. Careful implementation of analyticity and unitarity in the scattering amplitude is necessary and will be discussed in the next section.
We illustrate Lüscher's formula for studying σresonance phase-shifts in the I(J P C ) = 0(0 ++ ) scattering channel. The corresponding irrep for the S-wave scattering channel is A + 1 in both the O h and D 4h groups. The A + 1 irrep in D 4h couples to angular momenta J = L = 0, 2, 4.... Using the argument for the angular momentum cutoff we discussed in section II A, we assume that the contribution from J > 0 is negligible. As a result, we consider the A + 1 irrep under the condition that it is dominated by J = 0. Lüscher's formula in A + 1 is then given by The boost along the elongated direction does not change the symmetry group of the elongated box. Therefore, the boosted version of Lüscher's formula has the same form as Eq. (17) in A + 1 irrep of the D 4h group but with the modification that comes from the boost factor. The detailed derivation can be found in Ref. [9].

III. ANALYSIS OF THE ππ SCATTERING AMPLITUDE
The extracted phase-shifts and their covariances carry the full information about the interaction of (unphysically heavy) pions at discrete values of energy (or momentum) above and even below the corresponding ππ-threshold. From the point of view of scattering theory these points are interconnected by a function of energy (partial wave  projection is assumed) which carries specific analytic properties. These properties are the guiding principle for the construction of fit functions to extract its parameters from lattice data (as described in Section II and collected at the end of this paper in Appendix A.) Given the S-matrix and scattering amplitude T for two-to-two scattering with S = 1 − i T , the unitarity constraints the imaginary part of the scattering amplitude projected to definite isospin (I) and angular momentum (L) to be 16π Im T −1 IL (s) = 1 − 4M 2 π /s, where s denotes the square of the total four-momentum of the system. The above equation does not fix the amplitude entirely, but only up to a real-valued function in the physical region, where K IL (s) is a real-valued function, and G(s) is the two pion loop-function. In dimensional regularization it reads where p cm is the modulus of the three-momentum in the center-of-mass system. The regularization scale and the subtraction constant are fixed throughout this study and in accordance with the discussion of Refs. [9,28] to µ = 1 GeV and a(µ) = −1.28, respectively. This value results from a fit to experimental data but it can be varied in a range of ±0.5 without changing the result noticeably, i.e., the change is well absorbed into the values of the low-energy constants even if the amplitude is not explicitly scale invariant [28]. Note that since 16π Im G(s) = − 1 − 4M 2 π /s, the above formulation of the T -matrix differs from the usual K-matrix formulation only by a re-shuffling of the Re G(s) part. The form of the function K IL (s) is not fixed by unitarity. In this work we will use four versions of two different types of the K-matrix to gauge the systematic uncertainty tied to a particular choice.
For the direct use of lattice data, i.e., energy eigenvalues and their covariance matrices, we formulate the finite volume version of the scattering amplitude, i.e., of Eq. rest frame, q = 2π/L (n x , n y , n z /η). For boosts with momentum P, one performs a Lorentz transformation, see Ref. [31]. Note also thatG(η, s) is independent of q max but does depend on the subtraction constant a(µ) via G(s) as given in Eq. (19). In other words, the infinite-volume extrapolation is cut-off independent and equivalent to the Lüscher formalism up to exponentially suppressed contributions. At the same time, the functionG contains a dispersive real part relevant for the parametrization of the infinite-volume amplitude itself. See Ref. [31] for further details. The discrete energy spectrum obtained from the finite volume scattering amplitude (see Eqs. (18) and (20)), depends explicitly on the form of the K-matrix. In the following we will specify different parametrizations for it, which in every case depend on some free parameters. These will be adjusted by minimizing where s i are the fit parameter-dependent solutions of Eq. (21), ordered in a vector in Eq. (22), √ s 0 indicates the vector of eigenenergies measured on the lattice, and C is the covariance matrix of the correlated data. Eq. (22) implicitly contains a summation over contributions to the χ 2 from different elongation factors η and boosts P. The covariances between data from different moving frames are taken into account. We analyze data from different pion masses and channels individually and also simultaneously. In cases the data of the isoscalar channel (from this work) and of the I = L = 1 ρ-channel [9] are simultaneously fitted, Eq. (21) changes accordingly for the data from the ρ channel, K 00 → K 11 , and Eq. (22) becomes a sum over the two channels. Similarly, it becomes a sum over contributions from different pion masses in the respective simultaneous fits. Further statistical tests, such as Pearson's χ 2 test will be discussed below.

A. Conformal mapping
The first type of the parametrization of the scattering amplitude relies on a general form of K-matrix as an analytic function of energy. As discussed in Refs. [25,26] the convergence of a power series in energy is limited but can be improved, mapping it onto the interior of a disk limited by the right-and left-hand cuts lying on the boundary circle. We use two versions of such a mapping slightly adapted to our approach. We refer to them as [cm1] and [cm2] with the respective expansion variable where α, s thres , and c are parameters of the mapping. In particular, s thres in parametrization [cm1] is the position of the next threshold opening above the ππ threshold. Here, we do not have a KK channel as in Refs. [25,26], but we can interpret the parameter to take account of the opening of the four-pion threshold. The [cm2] parametrization is obtained in the limit s thres → ∞. The quantity √ c is the expansion point in the √ s plane connected to α through c = s thres α 2 /(1 + α 2 ). For the [cm1] parametrization, √ c = 778 MeV (for α = 1 and √ s thres = 550 MeV); for [cm2], √ c = 1 GeV. We have checked that the results do not depend on these choices.
For the isoscalar channel (I = L = 0) the chiral symmetry dictates that T 00 (s) must vanish for s = s A ∼ M 2 π /2. To account for this fact the K-matrix in this channel takes for both versions of the mapping variable ω [..] the following form where throughout the further calculations s A is set to its leading chiral order value, i.e. 2s A = M 2 π . The number of polynomial terms in the latter equations (∼ B i ) is not restricted a priori. For the present case two polynomial terms (B 0 , B 1 ) turn out to give sufficient flexibility in the energy variable ω(s) to fit the data. The lattice data consist of the energy eigenvalues and covariance matrices of the energy levels at two different pion masses (M π = 227 and M π = 315 MeV) in the I = L = 0 channel which we fit individually.
Here and in section III B the lattice energy levels for the σ-meson are fitted up to √ s ≈ 970 MeV and √ s ≈ 1070 MeV for the light and heavy pion masses, respectively. This corresponds in both cases to p 2 ≈ 0.187 GeV 2 for the magnitude of the center of mass three-momentum p. The data at higher p 2 cannot be fitted by the considered parametrizations. Fitting with more flexible parametrizations would, however, require more data points in the high momentum range.
The results of the fits are collected in the first four entries of Table IV, which all pass the Pearson's test with the total χ 2 ≈ 9 lying inside of the 80% confidence interval of our two-tailed test, i.e (6,19) and (7,20), for the light and heavy pion masses respectively. The error bands depicted in Fig. 2 are obtained from the error ellipses of the fitting parameters. Such a procedure for the propagation of statistical uncertainties is used in the remainder of the paper. There is a clear overlap of the fits with the data. Only for the heavy pion mass and energies deep below threshold there is some discrepancy that could be significant. Note again that the phase-shifts are not fitted directly, but rather the energies extracted from lattice QCD which makes the fit results and the shown phase-shift data difficult to compare; for example, the error bars in x and y-direction of the phase-shift data are perfectly correlated (one may think of inclined error bars), and the correlations between different phase-shift data can obviously not be visualized.
Finally, having fixed parameters of both parametriza-tions we perform an analytical continuation to the complex energy plane. On the second Riemann sheet of this plane we determine the position (z 0 ) and residuum (g 2 ) of the σ-resonance pole. The results are collected in the last three columns of Table V. They are discussed in Section III C together with the results from the chiral unitary approach discussed in the next section. The second option for the form of the K-matrix explored in this work is inspired by chiral perturbation theory in two-flavor formulation. As such it contains symmetries of QCD, while it also relates different interaction channels as a full fledged quantum field theory. These two facts allow for a chiral extrapolation between different pion masses as well as for a simultaneous description of the isoscalar and isovector channels. However, the price to pay is that the full K-matrix would contain infinitely many terms. To make a practically feasible approach, a truncation of a chiral series is required. Different versions are in use; their theoretical properties are discussed in detail in Ref. [41]. In the following we will use a version of chiral unitary approaches (UChPT), which, on the one hand, allows to address both (σ and ρ) channels of the ππ scattering simultaneously, see, e.g., Refs. [9,27,28,42]. On the other hand it contains only local terms (of the leading and next-to-leading chiral order) which makes the analysis of the discrete finite volume spectrum feasible in the same way as in the case of conformal mapping, see Section III A. The K-matrix reads in both considered channels where f π is the pion decay constant fixed to its value at the given pion mass (using the data from Table III). The low-energy constants (LECs) of the next-to-leading chiral order read [43] L a = −36l 1 + 44l 2 + 20(5L 2 + 6L 6 + 3L 8 ) , The model used here is identical to the one of Ref. [9], where using the potential of Eq. (26), the energy levels of the ρ-meson were analyzed. The same data set for the ρ-meson is considered here as well. Therefore, we refer to Ref. [9] for a detailed discussion of the results concerning the ρ-meson. Here, we perform a combined fit of the energy levels at the two given pion masses (M π = 227 and 315 MeV) in the σ-resonance channel, which depends on the three combinations of LECs (fitting parameters) of Eq. (25), L a , L b and L c , and combined fits of the σ and ρ channels (at one and two pion masses), being the fitting parameters in this casel 1 ,l 2 , L 2 and L 68 := L 8 + 2L 6 . Note that in Eq. (27) L i are LECs of three-flavor ChPT [43], which, however, appear here only in four linear independent combinations corresponding to the LECs of two-flavor ChPT [44].
In the following, we refer to [chm1] and [chm2] being the fits in the isoscalar or both isoscalar and isovector channel, respectively. The lattice data for the isovector channel includes only energy eigenvalues, corresponding to the center of mass energies √ s ∈ {m ρ − 2Γ ρ , m ρ + 2Γ ρ } where m ρ , Γ ρ are the mass and width of the ρ-resonance. For the detailed discussion see Ref. [9] The best fit parameters are collected in Table IV, see the entries "chm1 σ 227,315 " and "chm2 σ 227,315 ρ 227,315 ". These fits pass the two-tailed Pearson's test with χ 2 ≈ 29 for [chm1], and χ 2 ≈ 44 for [chm2], which both are inside of the corresponding 80% intervals of (17,36) and (29,52), respectively. The corresponding phase-shifts are depicted in Fig. 3.
For both best fits (fitting in both cases M π = 227 and 315 MeV simultaneously) we perform an extrapolation to the physical point, predicting the phase-shifts and comparing them with the experimental data from Refs. [35][36][37][38][39][40]. The result is depicted in Fig. 4, and shows a good agreement with the experimental data below energies of 950 MeV. Some deviation from the data becomes evident at higher energies, associated with the presence of f 0 (980) not captured by the two-flavor parametrization of the K-matrix. Furthermore, we perform an analytic continuation to the second Riemann sheet of the complex energy plane, determining the pole position and the coupling to the ππ channel. The results including a prediction at the physical point are collected in the rows 5-8 of Table V. In principle, the chirally inspired approach used here allows for a prediction of the full chiral trajectory of the pole positions and their residua, which will be discussed in the next section.
In tables IV and V also other fits are quoted. In particular, the chiral unitary approach has been fitted to the lighter pion mass alone ("σ 227 , ρ 227 ") and to the heavier pion mass alone ("σ 315 , ρ 315 "). For the low-energy constants we observe agreement between these cases, and also between these cases and the combined fit to both masses. A similar agreement is observed for the pole positions and residues which shows that the data are consistent under the fit hypothesis and, independently of the pion masses in the lattice calculation, lead to similar predictions for the σ properties at physical pion masses. As demonstrated in Figs. 2 and 3 all considered parametrizations lead to a good agreement with the fitted lattice data. Also the experimental data from Refs. [35,[37][38][39][40]45] lies in the 1σ error band around the extrapolated phase-shift from both parameterizations [chm1] and [chm2].
The analytical continuation of the scattering amplitude to the second Riemann sheet reveals the presence of a pole corresponding to the σ-resonance. The real and imaginary parts of the pole position as well as the residuum give the information about the mass, width and coupling to the ππ channel of this exited state, respectively. The results are collected in Table V. For comparison, we also quote there the result of the dispersion analysis of the experimental data from Ref. [ For the higher pion mass, the pole position in the individual [chm2] fit for this mass is in good agreement with the result of the conformal parameterizations. In the chiral fits to both pion masses combined, "[chm1], σ 227 , σ 315 ", the data from the lighter pion mass pushes the state towards the real axis making it, for some value in the uncertainty area, a virtual bound state. As a result, the pole position and residue have large uncertainties (see Fig.  5 and explanation below). Overall, the systematic uncertainty tied to the use of one or another parametrization appears to be smaller than the statistical one.
With this in mind we make a prediction of the σ pole position and the corresponding coupling to the ππ channel as a continuous function of the pion mass based on [chm2] fitted to both sets of lattice data, σ and ρ, and both pion masses simultaneously. To remind the reader these sets are obtained from calculations at M π ≈ 1.65 M phys π and M π ≈ 2.3 M phys π . The result of the extrapolation is depicted in Fig. 5 and exhibits three major regions: 1. With pion mass increasing from the physical one, the σ-resonance becomes lighter in units of pion mass, but couples more strongly to the ππ channel. The same happens with the reflected pole in the positive half-plane of the second Riemann sheet. At M π ≈ 2.5 M phys π both poles meet at the real energy axis below threshold on the second Riemann sheet (i.e. w.r.t. the cms-momentum -positive imaginary halfaxis), becoming virtual bound states.
2. With higher pion mass, both poles evolve on the real axis towards and away from the ππ threshold (p = 0), respectively. At M π ≈ 3M phys π one pole reaches the two-pion threshold, where the coupling g vanishes. 3. At higher pion masses the σ-resonance becomes a bound state, thus appearing as a pole on the first Riemann sheet below threshold or negative imaginary half-axis in the complex cms-momentum plane. The coupling to two pions increases in this region monotonically.
The behavior described above has been pointed out first in Ref. [29] using experimental data only and later in Ref. [20] with the input from lattice data of Refs. [12,13,17]. The present study favors a bound σ state at a pion mass of around M π ≈ 3M phys π . However, we note that the uncertainties grow rapidly with increasing pion mass as can be seen in the size of green 1σ error area on pole positions and error bars on g at heavy pion mass (see Fig. 5).
From Table V we note that the un-extrapolated results of conformal parametrizations ([cm1] and [cm2]) indicate that the systematical uncertainty tied to the choice of the (K-matrix) parametrization and approximation of the left hand cut is well under control, at least at the present level of statistical uncertainty of the lattice data. We have estimated that systematical uncertainty due to fitting window and lattice spacing in the calculation of the lattice data may have higher importance and lead to several percent uncertainty for the mass and width of the σ-resonance.

IV. SUMMARY
We performed a calculation of the phase-shifts in the isoscalar/scalar ππ channel in the elastic region for two quark masses corresponding to M π = 227 MeV and 315 MeV. For each quark mass we used ensembles with three different volumes to help us determine the phaseshifts in different kinematic regions.
To extract the parameters of the σ-resonance, we have to use a parametrization that satisfies physical constraints, in particular unitarity, analyticity, and proper chiral behavior. To gauge the systematics associated with the choice of such parametrizations, we used two types of approaches (each in several variants): a generic one that makes no assumption about the underlying dynamics, and a chiral perturbation theory inspired one that allows to extrapolate the resonance parameters to different (i.e., physical) pion mass.
The systematic errors associated with the choice of parametrization are about 10% for the pole position. The other sources of systematic errors that we assessed are the discretization errors and the fit-window for the extraction of the lattice QCD energies. The discretization errors are estimated based on the spread of lattice spacing results determined using various methods to calibrate it; we estimated this to be around 2-3% (see [9] for more details). For the fit window systematics, we computed the pole position based on energies extracted using slightly shifted fit windows; the shift on the pole position was at the level of 1%.
One of the strengths of the chiral parametrization is that it allows us to fit simultaneously both the σ and ρ channel, for both pion masses. We find that the model describes the data well and that the results extracted from the simultaneous fit to both channels agree well with the σ-channel fit results. We use the combined fit to extrapolate to the physical point and, based on the position of the pole in the complex energy plane, we find that M σ = (440 +10 −16 (50) − i 240 (20)(25)) MeV. Here the first error is the stochastic error and the second one is the combined systematic error discussed above.
The extrapolation to the physical point agrees with the experimental phase-shifts and the pole mass and width of the σ is compatible with the result of recent analyses based on experimental data.  corresponding to m π ≈ 315 MeV and on the right E 4,5,6 corresponding to m π ≈ 227 MeV. The order of the levels in each matrix corresponds to the order they appear in Table VI.