Two-loop analysis of the pion-mass dependence of the $\rho$ meson

Analyzing the pion-mass dependence of $\pi\pi$ scattering phase shifts beyond the low-energy region requires the unitarization of the amplitudes from chiral perturbation theory. In the two-flavor theory, unitarization via the inverse-amplitude method (IAM) can be justified from dispersion relations, which is therefore expected to provide reliable predictions for the pion-mass dependence of results from lattice QCD calculations. In this work, we provide compact analytic expression for the two-loop partial-wave amplitudes for $J=0,1,2$ required for the IAM at subleading order. To analyze the pion-mass dependence of recent lattice QCD results for the $P$-wave, we develop a fit strategy that for the first time allows us to perform stable two-loop IAM fits and assess the chiral convergence of the IAM approach. While the comparison of subsequent orders suggests a breakdown scale not much below the $\rho$ mass, a detailed understanding of the systematic uncertainties of lattice QCD data is critical to obtain acceptable fits, especially at larger pion masses.

On the technical level, the failure to produce resonant states is related to the fact that unitarity is only restored perturbatively in ChPT, so that any description of resonances requires a unitarization procedure. A widely used approach known as the IAM achieves this unitarization by studying the unitarity relation for the inverse amplitude [37][38][39][40][41][42][43][44][45]. In particular, in the case of SU (2) ChPT the IAM procedure can be derived starting from a dispersion relation in which the discontinuity of the left-hand cut is approximated by its chiral expansion [40,41]. While Adler zeros induce a modification for the S-waves [46], the naive derivation of the IAM survives for the P -wave amplitude: writing the partial wave for ππ scattering t(s) as t(s) = t 2 (s) + t 4 (s) + t 6 (s), (1) with the subscripts indicating the chiral order, the unitarized amplitude at next-to-leading order (NLO) becomes [38][39][40] t NLO (s) = t 2 (s) while at next-to-next-to-leading order (NNLO) [41,44] t NNLO (s) = t 2 (s) 2 t 2 (s) − t 4 (s) + t 4 (s) 2 /t 2 (s) − t 6 (s) .
To assess the chiral expansion of the unitarized amplitude beyond the first term, one thus needs the partial-wave amplitudes at two-loop order [47]. The IAM has been applied to study resonance properties at unphysical pion masses at one-and two-loop order as early as in Refs. [48][49][50][51], with numerous subsequent works confronting the IAM predictions with lattice data [52][53][54][55][56]. However, apart from Refs. [50,51] such studies have been restricted to one-loop order, so that it was not possible to scrutinize the convergence properties of the expansion in the pion mass.
The reason for this situation was twofold: first, while the one-loop amplitudes can be given in analytic form, similarly compact expressions were not available for the two-loop amplitudes, thus complicating their implementation considerably. Second, as shown in Refs. [50,51], the increased number of low-energy constants (LECs) renders the fits more volatile, so that lattice data need to reach a sufficient quality to allow for meaningful two-loop fits. In this paper we address both points: we present compact analytic expressions for the two-loop amplitudes that are straightforward to implement and devise a strategy for stable two-loop fits to current lattice data. While expressions are provided for all partial waves up to J = 2, we concentrate on the application to the ππ P -wave, including the resonance parameters of the ρ meson and its pole residue.

II. PARTIAL WAVES IN CHPT
We express the partial waves t I J (s), where I and J stand for the isospin and angular momentum, respectively, in terms of the pion decay constant in the chiral limit F as well as the pion mass M π (including quarkmass corrections from the LEC l r 3 ). We will follow the conventions of Refs. [3,57] for the one-loop LECs l r i and the two-loop LECs r r i . First, the leading-order (LO) results are [58] t 0 0 (s) 2 At NLO, the partial-wave amplitudes can be written in the form c IJ li (s)l r i , (5) in terms of L(s) = log 1 + σ(s) 1 − σ(s) , σ(s) = 1 − 4M 2 π s , (6) and coefficients as listed in App. E. We find that the NNLO expressions can be brought into a very similar form where σ ± (s) = 2σ(s)/(σ(s)±1) and in addition to powers of L(s) also polylogarithms Li n appear. The contributions from the NNLO LECs are collected in P IJ (s) and the imaginary parts determined by perturbative unitarity Im t 4 (s) = σ(s) t 2 (s) 2 , Im t 6 (s) = 2σ(s) t 2 (s)Re t 4 (s).    [66,67], although the uncertainties especially in r a,b are substantial and difficult to quantify.

III. FITS TO LATTICE DATA
From here on, we focus on the P -wave of ππ scattering, with both isospin I and angular momentum J equal to one. Its phase δ(s) = arg(t 1 1 (s)) can be computed using lattice QCD via Lüscher's quantization condition [1,68], which allows one to determine the phase shift given ππ energy levels and vice versa. To illustrate the fitting strategy as well as the conclusions regarding the pion-mass dependence of δ and the ρ parameters, we analyze such energy levels as computed on the lattice by two different groups: first, the one from Ref. [17], based on gauge configurations generated by the CLS collaboration, accompanied by a determination of the pion decay constant [69]. There are six data sets (ensembles) at five different pion masses in the range 200 MeV to 284 MeV. Second, we consider the energy levels from the HadSpec collaboration [12,70], using one of their ensembles with M π ≈ 236 MeV and two with M π ≈ 391 MeV. Both lattice calculations involve N f = 2 + 1 flavor simulations, but in either case the changes compared to the physical kaon mass, which determine the corrections to the LECs in two-flavor ChPT [71,72], are negligibly small compared to other sources of uncertainty. In the following, we concentrate mainly on the fit to the CLS data; a detailed description of the fitting procedure as well as an overview over the lattice data is given in App. A, while the fit to the HadSpec data is discussed in App. D. To reduce the impact of scale-setting uncertainties, i.e., the error that arises when determining the lattice spacing in physical units, we work in lattice units wherever possible.
The fit proceeds as follows: at NLO, Eq. (2) is used to compute the phase δ, which is subsequently inserted into Lüscher's quantization condition to determine the energy levels. Their distance to the energies as computed on the lattice is then minimized. Simultaneously, the pion decay constant is fit, using the ChPT expression given in App. E 3, truncated at NLO. In an NNLO fit, the same procedure is applied, with Eq. (3) instead of Eq. (2) and the pion decay constant truncated at NNLO. This means that at NLO only the LECs l r 2 − 2l r 1 and l r 4 appear, while the NNLO expressions depend on l r 1-4 as well as r a,b,c and r r F . The minimization of the χ 2 with respect to the fit parameters-most importantly the LECs-requires a sufficiently powerful algorithm. To find the global minimum, we first employ the differential evolution algorithm [73], whose results are subsequently refined via a modification of Powell's method [74]. The former algorithm allows one to tackle the multi-dimensional, nonlinear optimization problem at hand in a both robust and efficient manner, if its parameters are adjusted carefully. Together with the improved lattice data the choice and tuning of this algorithm are crucial to obtain sound fits that are stable even when ensembles at only a few different pion masses are available, e.g., the two masses used by HadSpec.
There are three sources of error that need to be considered for a reliable uncertainty estimate. First, the statistical error of the lattice data. Second, the error of the lattice spacing, which enters the ChPT expressions indirectly via the renormalization scale µ, see App. B. Third, the error that arises as a result of the truncation of the chiral expansion (1), which we are able to study in detail by a comparison of the IAM at one-and two-loop order. The chiral expansion proceeds in s/M 2 ρ as well as α = M 2 π /M 2 ρ , with the breakdown scale expected to be set by the ρ mass since it is the lowest-lying resonance in the partial wave of interest. The energy dependence is resummed by the unitarization via the IAM, leaving the expansion in the pion mass as the most critical variable. Following Ref. [75], we estimate the truncation error of an observable X as

IV. RESULTS
To fix the LECs it is necessary to control both the s dependence and the mass dependence. Hence, we fit all CLS ensembles from Ref. [17] simultaneously, once working to NLO and once working to NNLO, excluding only the ensemble N401 from the fit, since its pion decay constant has not been determined in Ref. [69]. To render the  including the goodness of the fit, the properties of the ρ resonance at the physical point, as well as the decay constant in the chiral limit. The first error is the statistical one, the second stems from the lattice spacing, the third is the truncation error estimated via Eq. (9). NNLO fit stable, it is necessary to put a constraint on the LEC l r 3 . This parameter governs the relation between the pion mass M π and its value M at LO in ChPT, information on which is not included in our fit. Thus we add a penalty term to the χ 2 that favors values of l r 3 around its reference value 0.8(3.8) × 10 −3 [59]. The LECs obtained at NLO are given in Table I, and the NNLO ones in Table II.
Since the amplitudes as given in Eqs. (2) and (3) have the appropriate analytic structure, they can be continued analytically to the second Riemann sheet, where the pole associated with the ρ resonance is located. Extracting the mass M ρ and width Γ ρ from the pole position s p via s p = (M ρ − iΓ ρ /2) 2 and the coupling g of ρ to ππ from the residue r via g 2 = 48πr/(4M 2 π − s p ) yields the values shown in Table III. Also shown are the goodness of the fit as well as the obtained value of F , the pion decay constant in the chiral limit. The corresponding phase is depicted in Fig. 1. Here and in the following, the physical point is simply defined by the PDG value of the charged pion mass, M π = 139.57 MeV [76], and F is computed using the PDG value of F π as input.
Due to the unitarization via the IAM, the LECs are The pion-mass dependence of the decay constant, the coupling, as well as the real-and imaginary part of the ρ pole as determined via fits to the CLS data, with error bands corresponding to (in order of decreasing color saturation) the data error (statistical plus spacing), the truncation error, and the total one. The dashed lines mark the physical pion mass. The decay constant is given in units of F to reduce the impact of the scale setting. Since the NLO and NNLO fits yield different values of F , their physical points in these units differ. Also shown as black ranges are reference values, the ρ characteristics taken from Ref. [77] and the decay constant from Refs. [60][61][62][63][64][65]76]. expected to deviate to some extent from the ChPT reference values [40][41][42][43]. Accordingly, all LECs agree well with expectations, apart from a large discrepancy in l r 4 both at NLO and NNLO. To understand its origin, we performed an NLO fit to the pion decay constant alone (at NNLO the fit becomes underconstrained), leading to l r 4 = 1.3(1.0) × 10 −3 , in agreement with FLAG, but already in tension with phenomenology. The remainder of the pull displayed in Table I originates from the ππ data. This pull becomes exacerbated at NNLO, but as indicated by the large uncertainties the sensitivity to l r 4 is limited. Indeed, we observe only a moderate increase of the χ 2 if literature values of l r 4 are enforced, as well as a large change to l r 4 = −16 × 10 −3 when employing a different strategy for the scale setting [17]. We conclude that there is a tension between the pion decay constant in the chiral limit and ρ parameters, which at least in part may be related to scale-setting uncertainties. Further details are given in App. C.
In general, we note that the χ 2 /dof improves significantly when going from NLO to NNLO, although a statistically fully acceptable fit would require a more detailed understanding of lattice artifacts. Comparing the obtained ρ characteristics with the ones from Roy-like equations [77] −0.07shows that both the NLO and NNLO results are compatible with these already within statistical errors, with a 1.4σ discrepancy in the width at NNLO and a 2.2σ tension in Im g at NLO. However, only the NLO value of F is compatible with the literature value F = 86.89(58) MeV, which is obtained by combining the PDG value of F π [76] with the FLAG N f = 2 + 1 average of F π /F [60-65].
Our main results are shown in Figs. 1 and 2, for the pion-mass dependence of the phase shift, the decay constant, and the ρ resonance parameters. Most notably, the two-loop analysis allows us to improve the precision con-siderably when going beyond the physical point, once the truncation becomes the dominant source of error. Second, with error bands produced assuming a breakdown scale of M ρ , the NLO and NNLO bands mostly overlap, which indicates that the true breakdown scale of the theory may lie below the ρ mass, but not by much.
Overall, the coupling shows a very mild mass dependence [48,50], as does the ρ mass. Towards the end of the fit range, the central value of the two-loop curve seems to decrease, in disagreement with the phenomenological expectation from both the KSFR relation [78,79] and the expected ordinary qq nature of the ρ meson [80]. This may again, in addition to the χ 2 and the tension in l r 4 , point to the impact of lattice artifacts, which the two-loop IAM becomes flexible enough to mimic.
Similar conclusions can be drawn from the fits to the data by HadSpec, which are described in App. D. They also show a significant improvement of the χ 2 when going from NLO to NNLO and a pion-mass dependence that mimics the one depicted in Fig. 2, with the difference that M ρ does not decrease at high pion masses, providing further evidence that this decrease may arise due to lattice artifacts. Notably, for the HadSpec data the ρ properties at the physical point are closer to the literature values at NLO than at NNLO.

V. CONCLUSIONS
In this work we have presented compact analytic expressions for the two-loop partial-wave amplitudes for ππ scattering up to D-waves, with a first application to an analysis of lattice data for the P -wave amplitude and the ρ parameters. We have shown that two-loop fits do improve the fit quality and, by comparing NLO and NNLO results, found that the breakdown scale of the chiral ex-pansion should not lie much below the expected scale set by the ρ mass. However, we also concluded that the current data sets cannot be described in a statistically satisfactory way, with a more detailed understanding of the lattice data required.
In the future, anticipated improvements in the precision of lattice QCD calculations will increase the need to match that precision in the analysis. In this work, we have demonstrated how to achieve two-loop precision in practice, using the example of the P -wave, but once lattice calculations mature a similar analysis can be performed for other partial waves including the pion-mass dependence of the f 0 (500). Even once data sets at the physical point become available, the IAM will thus provide a tool for a high-precision analysis of lattice data.

Acknowledgments
Finite-volume energy levels taken from Refs. [12,70] were provided by the Hadron Spectrum Collaboration (HadSpec)-no endorsement on their part of the analysis presented in the current paper should be assumed. In addition, we would like to thank John Bulava as well as ETMC for sharing lattice data with us. We thank Carsten Urbach for many useful discussions, and Raúl To understand the precise definition of the χ 2 , first it is expedient to briefly recall Lüscher's method in the context of the P -wave of interest. The procedure is roughly as follows [1]: on the lattice, operators corresponding to discreteworld versions of ππ states with I = J = 1 are constructed. These states are characterized by their irreducible representation (irrep) of the residual rotational symmetry as well as their relative momentum with respect to the rest frame, i.e., the boost momentum d ∈ Z 3 . Subsequently, the energy levels of these states are computed. Such an energy level E lat is related to the scattering phase shift δ via the quantization condition where Z is a known, albeit complicated, expression that can be computed only numerically and depends on the details of the lattice, the pion mass, as well as the irrep and the boost of the state with energy E lat . A given data set (ensemble) contains several energy levels E lat i , i = 1, . . . , N , corresponding to different combinations of boost and irrep as well as different excitations of the states. Moreover, it has a fixed pion mass and, if it has been measured, pion decay constant. Both the CLS collaboration and HadSpec have generated several data sets with different characteristics, the ones of relevance for this work are listed in Table IV. For simplicity, consider first a single ensemble in isolation. To fit the IAM to the data of this ensemble, we proceed as follows: for a fixed choice of LECs, irrep, and boost the phase δ IAM (E) = arg(t(E)), with t given by either Eq. (2) or Eq. (3), is inserted into Eq. (A1), which is subsequently solved numerically to determine an energy E IAM ,  corresponding to the fixed boost and irrep. This is done for all combinations of boost, irrep, and all excitations, yielding one energy level E IAM i for each level E lat i . The goal is then to minimize the distance of these energies via a variation of the LECs. More specifically, we introduce a χ 2 of the form and C the covariance matrix of the lattice data. Here F ChPT π is given by the expressions in App. E 3, L 1 and L 2 are the set of LECs on which the model expressions for the decay constant and the phase shift depend, respectively, and M π as well as F are additional fit parameters. That is, the fit parameters for a single ensemble are F = {F } ∪{M π } ∪L 1 ∪L 2 . M π is introduced as a fit parameter to take into account the error of the pion mass, while treating F as a fit parameter allows us to work almost exclusively in lattice units, see App. B.
We always fit multiple ensembles simultaneously, because in this way the pion-mass dependence of both the IAM and F ChPT π can be controlled, which becomes particularly important at NNLO, where several free parameters are present. Since the data on different ensembles are not correlated, the χ 2 is given as χ 2 = k χ 2 k , with each χ 2 k mimicking χ 2 single . If two ensembles share the same pion mass or decay constant, the corresponding entry is taken into account only once. Moreover, since F in physical units is the same for all ensembles, there is only one fit parameter F for each lattice spacing.
The minimization of the χ 2 requires repeated evaluations of the quantity Z in Eq. (A1), whose computation is numerically demanding. To accelerate it, we slightly modify the approach of Ref. [55]: assuming that the fit is not completely off, for each energy level E lat it will yield an energy level E IAM that is close to the former, i.e., (E lat − E IAM )/E lat ≪ 1. The same holds true for the mass. Hence it is justified to Taylor-expand both sides of the quantization condition δ IAM (E IAM ) = Z(E IAM ) around E lat and M lat π , up to and including linear order. Solving the result for E IAM yields Thus, for a given energy E lat and mass M lat π we need to compute Z and its first derivatives only once, we do so by numerically differentiating the full quantization condition (A1). Then we can compute E IAM for different values of the LECs by merely re-evaluating the IAM. This speeds up the minimization drastically. We explicitly checked that this approximation is justified by performing selected fits twice, once using the full quantization condition and once using Eq. (A3), both yielding the same minimum, i.e., the same values of the LECs. However, the value of the χ 2 obtained using the Taylor expansion tends to be slightly (at most a few percent) bigger than the exact value, hence to obtain the correct χ 2 all results are re-evaluated using the exact quantization condition. Note that this deviation in the χ 2 is expected: the χ 2 depends on the LECs only via the IAM (ignoring F ChPT π for the time being, since it is of no relevance for this argument), so we might write χ 2 ({E IAM i (L 2 )}). A successful minimization will result in values L min 2 for which χ 2 is at a global minimum. Equation (A3) then amounts to an approximation of the energy values for which χ 2 is minimal, therefore its value needs to increase.
Even with this acceleration at hand, the minimization problem provides a challenge. Namely at NNLO, there are eight LECs. In addition, there are the fit parameters M π and F , e.g., in a global fit to the CLS data four fit parameters corresponding to different pion masses and three fit parameters F corresponding to different lattice spacings are required. Altogether, there are up to 15 fit parameters (even more if the error of the lattice spacing is taken into account, see App. B). To still obtain reliable results, a powerful minimization algorithm is needed. The differential evolution algorithm [73] turns out to be a good choice, i.e., it allows for stable fits while still consuming not too much time. Roughly, it works via generating a population of trial solutions, recombining them randomly into a new generation to find solutions of a better fitness (in our case a lower value of the χ 2 ), and iterate this.
To make it work properly, it needs to be configured carefully, a procedure that requires repeated iterations of the same fits while varying the configuration parameters and observing the impact on stability, the resulting χ 2 , and performance. In addition to the population size, there are two core handles: the differential weight and the crossover probability, influencing how strongly the components and how many components of the trial solutions are modified in each iteration, respectively. For the problem at hand, population sizes of 20|F| for fits to the CLS data and 40|F| for fits to the data by HadSpec are sufficient, where F is the set of all fit parameters. With these population sizes, the crossover probability can be set to 1, while the differential weight is randomly changed to a value between 0.5 and 1 before each iteration, a technique known as dithering. The results obtained in this way are refined via Powell's method [74]; for both algorithms we use the implementation in SciPy [81]. It comes at no surprise that the data by HadSpec, which have a lower number of data points and pion masses as compared to the CLS data, require more generous settings of the differential evolution algorithm, i.e., a larger population size relative to the number of fit parameters.
Exemplarily, the energy levels as well as the phase as obtained in the NNLO fit to the CLS data are shown in Fig. 3 for a single ensemble.

Appendix B: Error computation and impact of lattice spacing
Taking into account the statistical error of the data is straightforward: both CLS and HadSpec provide each energy level E lat as a collection of several hundred values, each corresponding to a different sample of the underlying gauge configurations. Via jackknife resampling of these underlying values and repeating the fit on each jackknife sample the error can be computed taking into account the correlation of the energy levels automatically. The reason to pick a jackknife (i.e., drawing samples via omitting a single underlying value in each run) instead of a bootstrap (i.e., sampling the underlying values randomly with replacement) lies in the nature of the quantization condition (A1) [18]. It has poles at energies that correspond to free ππ states. Sometimes, the values E lat are so close to these poles that resampling the underlying value via a bootstrap yields central values on the other side of the nearby pole, thereby resulting in a completely different value of the phase. The use of jackknife resampling circumvents this issue, for jackknife samples are significantly more narrowly distributed than bootstrap samples. There is one subtlety stemming from the jackknife. It needs the same number of underlying data points for all energy levels. However, the HadSpec ensembles differ in this respect. Accordingly, in a global fit to the HadSpec data we feed the jackknife errors   into a parametric bootstrap, as suggested in Ref. [18]. The other source of error associated with the data is the scale-setting. Namely, on the lattice all quantities are computed in units of the lattice spacing a, the distance between two adjacent sites in one direction. To translate the quantities into physical units, they need to be multiplied by the appropriate power of a. Thus a needs to be determined in physical units, and this so called scale-setting is error-prone. For example, in Ref. [69] two different methods are used to set the scale of the CLS ensembles, strategy 1 (via the Wilson flow) and strategy 2 (via decay constants), both yielding different results that are incompatible within their errors. Hence throughout the entire fit we work in lattice units wherever possible. However, there is one place where the lattice spacing enters the fit: the renormalization scale µ appearing in the ChPT expressions for both the scattering amplitude and the pion decay constant, see App. E. It shows up both explicitly via logarithms and implicitly, since the LECs depend on µ in a way that renders the total amplitude scale-independent. Therefore, to have one single set of LECs that can be used on all ensembles in a global fit and for extrapolations of observables to the physical point, µ needs to be fixed globally, and we make the common choice µ = 770 MeV. When fitting, this requires the translation of µ into lattice units via µ → aµ inside the logarithms. This is the only place where the lattice spacing shows up in our fits.
To compute the impact of the error of a (with a fixed choice of scale-setting) on the fit results, we extend the χ 2 : That is, for each lattice spacing a lat i we introduce a fit parameter a fit i , with C a the covariance matrix of the lattice spacings. The fit is then repeated multiple times with different values of the lattice spacings, obtained via drawing samples from a multivariate normal distribution with means a lat i and covariance matrix C a . To estimate C a , we proceed differently depending on the data set.
The HadSpec data are obtained at two different lattice spacings, corresponding to the two different pion masses in Table IV. To set the scale, the mass of the Ω baryon as determined on the lattice [12,82] is divided by its experimental value, i.e., a lat i = M lat Ω,i /M exp Ω , i = 1, 2. The two spacings are correlated due to the common choice of M exp Ω . We take this into account by resampling M exp Ω using its PDG value and error [76], for each value obtained in this way we sample M lat Ω,i and compute a lat i , i = 1, 2. All samples are drawn parametrically assuming a normal distribution. Repeating this multiple times allows us to estimate C a in the standard way. Given the small error of M exp Ω , the off-diagonal entries of C a are an order of magnitude smaller than the diagonal ones.
The CLS data are obtained at three different lattice spacings (ignoring the ensemble N401, which is excluded from the fits since the corresponding pion decay constant is not available). The procedure in case of the scale-setting via strategy 1 is very similar to the one outlined for the HadSpec ensembles, with the quantity t 0 associated with the Wilson flow replacing the Ω mass and all relevant values given in Ref. [69]. However, since the error of the reference value of t 0 is larger than the one of M exp Ω , in this case the off-diagonal entries of C a are sizable. On the other hand, C a for strategy 2 is assumed to be diagonal.
The results described in Sec. IV were obtained using strategy 1, since it is the one preferred by the authors of Ref. [69]. Since the lattice spacings show up only via the renormalization scale in logarithms, the impact of their   errors is expected to be small. This is indeed what is found, as becomes clear in Tables I, II, and III: the statistical error always dominates the one stemming from the lattice spacing. Note also that the latter increases when going from NLO to NNLO, the reason being the more dominant role the logarithms containing µ play at NNLO, where they appear not only in F π , but also in the scattering amplitude.

Appendix C: Scale setting in CLS fits
As observed in Sec. IV, fitting the CLS data using strategy 1 for the scale-setting (see App. B) yields values of l r 4 both at NLO and NNLO that are in conflict with literature values. To understand this discrepancy, a detail of the scale-setting deserves further attention: as explained in Ref. [69], strategy 1 requires small shifts in the values of the pion masses and decay constants as measured on the CLS ensembles (since the statistical uncertainties in strategy 2 are significantly larger than in strategy 1, these required mass corrections are neglected in strategy 2). However, the ππ energy levels from Ref. [17] are obtained using non-shifted pion masses. Thus, to be consistent, we also use the non-shifted values of the decay constants and masses in the fit, at the price of being inconsistent with strategy 1. To assess if this inconsistency is the reason for the discrepancy in l r 4 , we re-perform both the NLO and NNLO fits, this time using the lattice spacings of strategy 2. In Tables V and VI, the resulting LECs are compared with the ones obtained previously. Clearly, both at NLO and NNLO the central value of l r 4 moves a little closer to its literature values, but the discrepancy remains sizable.
To check if the discrepancy can be further reduced, we again perform the strategy 2 fits, this time putting a constraint on l r 4 to enforce it to be close to its literature value 6.2(1.3) × 10 −3 [59]. At NLO, this leads to an increase of the χ 2 /dof, with the resulting value of l r 4 = 0.4 × 10 −3 , still not perfectly overlapping with the reference value. An increase of the χ 2 /dof can also be observed at NNLO, with the resulting value of l r 4 = 5.4 × 10 −3 compatible with the reference value. However, the values of the other LECs change considerably, most notably l r 2 = −0.21 × 10 −3 now being in conflict with the literature, a shift that can already be anticipated by having a look at the correlation c = −98 % of l r 2 and l r 4 in the strategy 1 fit. In addition, the value of F improves, namely its central value moves from 90.8 MeV to 87.2 MeV. On the other hand, the ρ characteristics at the physical point get worse, the pole is now located at To pin down the source of this behavior, we perform fits to the decay constants only as a function of the pion mass. This forces us to set the scale directly in the beginning, since otherwise the fit would be underconstrained. Even with the scale set, an NNLO fit is impossible, for at this order there are six free parameters. Hence we work at NLO, with F and l r 4 as the sole free parameters. The results are given in Table VII and depicted in Fig. 4. Within the statistical error originating in F lat π and M lat π , strategy 2 is compatible with the literature values of the two LECs, although it is in tension with the one from phenomenology. However, the fit using strategy 1 yields values that are incompatible with the literature, shifting the data does not improve the situation significantly.
These observations point to the following conclusion: independently of the scale setting strategy, the full IAM fits to the CLS ππ data require values of F that are larger than the literature value. This, in turn, pushes the values of l r 4 to unphysically small values. Since a fit to the decay constants only using strategy 2 produces values that are more in agreement with the literature, the discrepancy seems to emanate from the ρ data. However, the shift in l Ref. [59] Refs. [66,67]    X: Results of NLO and NNLO fits to the HadSpec data. The first error is the statistical one, the second the one due to the lattice spacing, the third arises from the error of the literature value of F , the final one is the truncation error.
going from strategy 1 to strategy 2 indicates that scale-setting effects do play an important role. We stress that this analysis is not the final word, an improved fit to the decay constants only (taking into account the error of the lattice spacing) and a better understanding of the overall lattice artifacts would be necessary to completely settle this issue. Ultimately, these observations are likely related to the fact that even at two-loop order the quality of the resulting fits to the lattice data cannot be considered statistically acceptable.

Appendix D: Fits to HadSpec data
Compared to the CLS data, the data by HadSpec have fewer energy levels at only two different pion masses, see Table IV. Nevertheless, stable fits are possible if more generous settings of the minimization algorithm are used, as explained in App. A. On top of this, it is necessary to set l r 3 to a fixed value, a constraint alone is insufficient, for the fit is insensitive to the literature value within its errors. Hence we fix l r 3 = 8.2 × 10 −4 . Another crucial difference is that the pion decay constant has not been measured on the ensemble with the lower pion mass. Contrary to what has been done in the case of the N401 CLS ensemble, excluding the ensemble from the fit is not an option, because this would leave only a single pion mass. Consequently, instead of fitting the pion decay constant, we set F to its literature value F = 86.89 (58) MeV [60][61][62][63][64][65]76]. Since F is needed in lattice units, the lattice spacing appears now both in the translation µ → aµ and F → aF , increasing its relevance for an error computation. In addition, the error of the literature value itself needs to be taken into account.
The numerical results of a global fit at NLO and NNLO are shown in Tables VIII, IX, and X, while the phase at the physical point as well as the pion-mass dependence of the ρ characteristics are depicted in Fig. 5. The conclusions drawn from the fits to the CLS ensembles carry over in large parts, with a few differences: due to the more important role of the lattice spacing, its error is sometimes dominant. In addition, the pole position of the ρ at the physical point at NNLO deviates significantly from its literature value, with both the mass and width being too low. On the other hand, the mass does not start to decrease at high pion masses, as opposed to the CLS fit and in agreement with phenomenology [78][79][80].
In Ref. [52], the NLO IAM was fit to the ensemble at M π ≈ 236 MeV only, yielding χ 2 /dof = 1.26 as well as M ρ = 755(2)(1) +20 −2 MeV and Γ ρ = 129(3)(1) +7 −1 MeV at the physical point. Comparison with our NLO fit (see Table X) shows that our χ 2 /dof is significantly larger, but at the same time Γ ρ comes out much closer to its literature value Γ ρ = 146.4 +2.0 −2.2 MeV [77]. The origin of these differences is twofold: first, the NLO fit presented in this work is not restricted to one ensemble, but instead three ensembles are fit simultaneously. Second, we express the ChPT partial waves in terms of F , while in Ref. [52] the amplitudes are expressed in terms of F π , which contrary to F depends on the pion mass, see App. E 3. In consequence, the NLO expression for the ππ amplitude in Ref. [52] depends on both l r 2 − 2l r 1 and l r 4 , but, as we argued here, the dependence on l r 4 represents a spurious higher-order effect (even the two-loop IAM does not depend on l r 4 once expressed in terms of F ). Since both LECs are fit only to ππ data, this additional freedom improves the fit, but at the expense of adjusting l r 4 not from the pion-mass dependence of F π , but from higher-order effects in ππ scattering. The resulting value of l r 4 = −28(6)(3) +1 −11 × 10 −3 indeed deviates considerably from its literature value, albeit in the same direction as we observed in the CLS fits, see App. C. We write the NLO partial-wave amplitudes as with The coefficient functions are In particular, there are no contributions involving l r 4 when the amplitudes are written in terms of F instead of F π . The P -wave amplitude depends on a single combination of LECs, 2l r 1 − l r 2 , which is scale independent.