Light- and strange-quark mass dependence of the $\rho(770)$ meson revisited

Recent lattice data on $\pi\pi$-scattering phase shifts in the vector-isovector channel, pseudoscalar meson masses and decay constants for strange-quark masses smaller or equal to the physical value allow us to study the strangeness dependence of these observables for the first time. We perform a global analysis on two kind of lattice trajectories depending on whether the sum of quark masses or the strange-quark mass is kept fixed to the physical point. The quark mass dependence of these observables is extracted from unitarized coupled-channel one-loop Chiral Perturbation Theory. This analysis guides new predictions on the $\rho(770)$ meson properties over trajectories where the strange-quark mass is lighter than the physical mass, as well as on the SU(3) symmetric line. As a result, the light- and strange-quark mass dependence of the $\rho(770)$ meson parameters are discussed and precise values of the Low Energy Constants present in unitarized one-loop Chiral Perturbation Theory are given. Finally, the current discrepancy between two- and three-flavor lattice results for the $\rho(770)$ meson is studied.


I. INTRODUCTION
The ρ(770) meson is the lightest vector meson in the hadron spectrum and one of the most studied hadrons in the literature. It is one of the best examples of a qq resonance well described within the quark model. Its phase shift fits well into a simple Breit-Wigner (BW) parameterization up to small corrections [1,2] and it is usually considered as the prototype of narrow resonance in the light-quark sector. It also dominates the ππ scattering amplitude in the I = J = 1 1 channel below 1 GeV, decaying almost exclusively to two pions [3]. The ρ-meson mass and width are well known from experiment; the Particle Data Group (PDG) quotes for their BW values M = 775.1(3) and Γ = 149.1 (8), respectively [3] 2 From the theory side, the most precise determination of its pole parameters comes from the Roy-equation analysis of ππ scattering [4][5][6][7][8]. The contribution of the ρ(770) is also important for the hadronic total cross section σ(e + e − → hadrons) [9][10][11], which explains applications that go well beyond low-energy meson physics, ranging from the hadronic-vacuum polarization and the light-bylight contributions to the anomalous magnetic moment of the muon (see, for instance, [12][13][14][15]) to the electromagnetic and tensor-nucleon form factors [16][17][18][19]. Fur- * Electronic address: raqumoli@ucm.es † Electronic address: elvira@itp.unibe.ch 1 Where I and J refer to isospin and angular momentum respectively. 2 These values correspond to the average value for the charged meson seen in τ decays and e + e − collisions.
Although, the ρ-meson properties fit well within the naive quark-model picture, its nature in terms of QCD degrees of freedom is still under discussion [31]. Nevertheless, at low energies QCD becomes non perturbative, what hinders the study of hadron composition in terms of the fundamental QCD degrees of freedom. LatticeQCD simulations attempt to tackle this problem, however, several challenges are met when dealing with hadron scattering processes [32,33]. In this regard, the m q and 1/N c expansions [34][35][36] 3 provide model independent predictions to identify different kinds or hadrons. These parameters can be used to study whether the response of resonance properties to a change on N c or m q compares well with the behavior expected for different QCD configurations. For instance, by studying the N c dependence of the ρmeson properties, it was found that it also has a small non-qq component [37][38][39][40][41]. In addition, the analysis of the quark-mass dependence of the ρ(770) parameters by means of the generalization of the Feynman-Hellmann theorem for resonances suggests that it requires nonnegligible corrections beyond the quark model [42]. In this way, the extraction of the light-and strange-quark mass dependence of the ρ-meson parameters from Lat-ticeQCD simulations provides a powerful tool to confront quark-model predictions.
At low energies, Chiral Perturbation Theory (ChPT) [43][44][45] is the Effective Field Theory (EFT) that controls the quark-mass dependence of hadronic observables. ChPT encodes the interactions of the pseudo-Goldstone bosons of the spontaneous chiral symmetry breaking, and hence, it is capable to describe the quark-mass dependence of the light-pseudoscalar meson masses and decay constants at low energies, being this completely inherited from QCD [44,45]. Nevertheless, ChPT is constructed as an expansion in quark masses and momenta and hence it is only valid below a certain scale. Therefore, ChPT does not provide direct information about resonance properties. On the contrary, unitarized Chiral Perturbation Theories (UChPTs) [46][47][48][49][50][51][52] are based on imposing exact unitarity while matching ChPT at low energies. Thus, the region of validity of the chiral expansion is extended, allowing one to generate poles on unphysical Riemann sheets in the complex-energy plane and to access resonance properties. In particular, the Inverse Amplitude Method (IAM) [47][48][49][50] generates the ρ(770) resonance from ππ scattering and provides a tool to study the lightand strange-quark mass dependence of the ρ-meson properties, while reproducing the chiral series at low energies. Thus, in this work we utilize the IAM to investigate the quark mass dependence of the ρ-meson pole parameters, such as its mass, width and couplings to the ππ and KK channels. The analysis of these properties requires the determination of the Low Energy Constants (LECs) involved in the pseudoscalar meson masses, decay constants and meson-meson scattering. A chiral trajectory specifies the way in which the light-and strange-quark masses vary. In UChPT, the behavior of resonance properties over chiral trajectories is controlled by chiral symmetry and unitarity. It is desired to determine LECs which provide a full description of the resonance properties on chiral trajectories where m s and/or m u,d vary. These kind of predictions of an EFT can be tested by lattice QCD simulations.
LatticeQCD (LQCD) is the only known tool to extract non-perturbative information from QCD. It is the instrument to determine the low-energy parameters of the chiral Lagrangian that govern the quark mass dependence of resonance properties, hence, rendering evidence of the EFT predictions. In recent simulations, lattice data on I = J = 1 ππ scattering have been extracted for several pion masses for two light flavors (N f = 2) [53][54][55][56][57][58][59] and including also the strange quark (N f = 2 + 1) [60][61][62][63][64][65][66][67]. See also [68,69] for recent N f = 2 + 1 + 1 simulations. Surprisingly, results for the ρ mass in N f = 2 simulations are at odds with experimental predictions. Namely, the N f = 2 simulation with the lightest pion mass m π 150 MeV by the RQCD Collaboration [58] predicts a ρ-meson mass around 60 MeV below the physical value [58]. Other N f = 2 simulations also show disagreement with the closest pion-mass result for N f = 2 + 1. For example, the N f = 2 GWU simulation at m π 226 MeV [59] gives a ρ-meson mass around 45 MeV lighter than the N f = 2+1 Hadron Spectrum (HadSpec) outcome of the simulation for m π 236 MeV [60]. It has been argued in recent analyses [31,59] with the UChPT model of [46], that this difference can be explained through the effect of KK loops in the ππ → KK → ππ reaction, where the kaon is off-shell. This effect has been shown to be consistent among the N f = 2 simulations [31]. Moreover, while the error ellipses of N f = 2 lattice data analyses do overlap, hence showing consistency among the simulations, the same cannot be stated for the N f = 2 + 1 results, where one finds inconsistencies among lattice simulations [31].
The light-quark mass dependence on decay constants has also been studied in LQCD simulations in [70][71][72][73][74][75][76][77]. However, almost no attention has been paid in the past to their strange-quark-mass dependence, nor of the ρ(770) phase shift. Most of these simulations have been performed with a strange-quark mass kept fixed at the physical point, m s = m 0 s . 4 A reflection of this can be found in the Flag Review [78]; the averaged LEC values from different fits to decay constants with ChPT do not represent a global analysis of data and they do not track other trajectories rather than those with roughly m s = m 0 s . In addition, these LECs do not describe the meson-meson interaction at the energies where the ρ meson begins to resonate. The only exception up to recently was a simulation of the pion decay constant done by MILC over the m s = 0.6 m 0 s trajectory [72]. Thus, in spite of the great advance of lattice simulations, data out of the chiral trajectory m s = m 0 s are still scarce, even though the response of hadron properties to different chiral trajectories could elucidate their strangeness nature, in particular, and dynamical nature, in general. A larger amount of highly precise data on a variety of chiral trajectories are necessary to shed light on the composition of hadrons.
Recently, the CLS Collaboration generated ensembles on different chiral trajectories with TrM = C, (where TrM = m u + m d + m s and C denotes a constant), in large volumes [70,79]. Moreover, this constant varies a little with the inverse gauge coupling, β, of the simulation, which characterizes the set of ensembles generated. These trajectories are of particular interest since the hadron response along the trajectory will manifest as a consequence of both, variations in the light and strange quarks. These recent lattice simulations motivate the present analysis by investigating them in combination with the simulations over m s = m 0 s trajectories. This provides a good ground to study the strange-quark mass dependence of decay constants and the ρ(770) phase shift, which we intend to do here. Of course, new LQCD simulations over trajectories with larger variations on these constants or for different values of the strange-quark mass in the m s = k trajectories would improve the analysis presented here.
The study of hadron properties (ρ meson) we conduct here needs to emphasize the role of pseudoscalar decay constants, which are strongly connected to the coupling of vector mesons to pions. This is supported by the assumption of dominance of vector mesons in the pion-photon coupling, the so-called Vector Meson Dominance (VMD) [80], which connects the size of the pion decay constant (f π ) and the ρ → ππ coupling (g ρππ ) in the EFT [81]. In this context, the large experimentally observed decay width of the ρ(770) meson is directly connected to its coupling to two pions, which explains why the ρ-meson phase shift and ρ-meson properties are tightly related to the size of f π . In this sense, these two observables should always be determined together in lattice simulations. Beyond that, the quark mass dependence of pseudoscalar decay constants fixes the chiral trajectories in the lattice and hence, they can be used to set the lattice scale by letting them go to the physical point.
The analysis we perform here will be useful to further check the KSFR relation [82], which under VMD states that g ρππ = m ρ / √ 2f π in the SU(3) limit where m u = m d = m s . While, it is common that in previous lattice/experimental data analyses of ρ-meson phase shift the ρ-meson mass increases monotonically with m π , so that the KSFR relation is fulfilled [59,83,84], this behavior was not observed in the recent data of [79]. Whether this is a consequence of the lightness of the strange-quark mass used in these simulations or not will also be checked in the present analysis.
Let us make some previous remarks. Experimental phase-shift data on ππ → ππ scattering in the I = J = 1 channel were successfully reproduced using the IAM [83,85]. Here, the LECs are extracted by performing a fit to lattice phase-shift data instead, taking into account the covariance matrix for energy levels, similarly as in [84]. The main differences with the work of [84] are: 1. A global fit to lattice data on two distinct chiral trajectories, TrM = C and m s = k, is done instead of considering trajectories only over m s = m 0 s simulations. 5 As mentioned previously, this includes data from [70,79] for TrM = C and from [60][61][62][71][72][73][74] for m s = k.
2. We perform a simultaneous fit of phase shift and decay constant lattice data. Note that in [84], only phase-shift data were analyzed, while the quark mass behavior of pseudoscalar decay constants was fixed with the LECs obtained in the fit done in [83], which only included lattice data on m s = m 0 s .
3. The theoretical framework used here is the IAM in coupled channels [85] instead of the simplified UChPT model considered in [46]. This is, we include here one-loop diagrams not just in the s channel, but also in the t and u channels, hence, consistently with chiral symmetry at low energies.
4. The systematic error in the lattice spacing is taken into account by using the bootstrap method, assuming that the lattice spacing is normally distributed with the standard deviation associated to the lattice error.
Although in the present analysis we only include N f = 2 + 1 lattice data, this work can be considered as complementary to the previous N f = 2 and N f = 2 + 1 analyses done in [31,59] and [84], respectively. If the strange-quark mass has no effect on the ρ-meson properties extracted from the lattice simulations, then, the disagreement among N f = 2 and N f = 2 + 1 lattice results will be due to the scale setting or other finite volume effects, such as the lattice spacing or the box size. Thus, we study in detail in which particular quark-mass regime the ρ-meson properties in the simulation are sensitive to both, the strange-and light u-, d-quark masses.
This paper is organized as follows. In section II we explain the formalism considered. In section III, we show the results of a global analysis on decay constants over several chiral trajectories. Section IV provides the results of the combined fit both to phase shift and decay constant lattice data. In particular, we first analyze in IV A lattice data over m s = k trajectories, while the same analyses for the TrM = C data are shown in section IV B. Following this, we present our final results on a global fit on both trajectories in section IV C. Finally, the main conclusions are presented in section V.

A. Chiral Perturbation Theory
At low energies QCD interactions become nonperturbative and EFTs provide the proper framework to perform systematic calculations. The basic premise of EFTs is that the dynamics at low energies (or large distances) do not depend on the details of the dynamics at high energies (or short distances). As a result, lowenergy hadron physics can be described using an effective Lagrangian containing only a few degrees of freedom, hence, ignoring those present at higher energy scales. Chiral perturbation theory is the low-energy EFT of QCD. It is built as the most general expansion in terms of derivatives and quark masses [44,45] compatible with QCD symmetries, which relevant degrees of freedom at low energies are the pseudo Nambu-Goldstone bosons (NGB) of the chiral symmetry spontaneous breakdown, i.e., pion, kaon and eta mesons.
At leading order (LO) in this expansion, the chiral Lagrangian reads where f 0 coincides with the pion decay constant in the chiral limit and χ = 2B 0 M, with B 0 a constant to be related with the quark condensate and M = diag (m ud , m ud , m s ) is the three-flavor quark-mass matrix, where exact isospin symmetry ) collects the contribution of pions, kaons and etas, with By expanding the LO chiral Lagrangian in powers of f 0 , one can identify the mass field terms obtained with the pseudo NGB fields, which yields a relation between meson and quark masses The constant B 0 is related with the quark condensate value in the chiral limit, with q ∈ {u, d, s}, leading to the well known Gell-Mann-Oakes-Renner formula 2m ud Σ 0 = M 2 0π f 2 0 [86], i.e., even though both m ud and Σ 0 are scale dependent quantities, and hence, they are not observables, their product is scale independent.
At higher orders, all terms in the Lagrangian come multiplied by LECs, which contain information about higher energy scales. In addition, they absorb the divergences which appear in the chiral expansion, so that, the theory is renormalizable order by order. Unfortunately, the LECs cannot be determined perturbatively from QCD. While the LECs which multiply energydependent terms can be extracted quite well from dispersion theory [87][88][89][90], Lattice QCD provides in principle a model independent way to determine the values of LECs which fix the quark mass dependence [78,91].
The NLO Lagrangian was first derived in [44] for two flavors. The effect of the strange quark was studied in [45]. Omitting field tensor and vacuum terms, the SU(3) NLO ChPT Lagrangian reads In Eq. (4), L 1 , L 2 and L 3 multiply massless terms and hence they also contribute in the chiral limit. L 4 and L 5 accompany terms depending linearly on the quark masses and they contribute to the renormalization of the NGB wave functions and decay constants. Lastly, L 6 , L 7 and L 8 come together with quadratic terms on the quark mass. These only contribute to the renormalization of the NGB masses and have a minor role in the determination of the ρ(770) meson properties.
One-loop correction to the pion, kaon and eta NGB masses read [45] with The superscript r denotes renormalized LECs, which carry the dependence on the regularization scale µ [45]. This scale dependence cancels exactly in the calculation of any observable. In the following, we will identify the physical NGB masses with the one-loop ChPT prediction above. Nevertheless, note that the quark mass dependence is always expressed in terms of the leading order NGB masses. In addition, while at LO the NGB decay constant f 0 is independent of the quark mass, one-loop corrections in the pseudoscalar decay constants lead to which are also identified with the physical quantities.

B. Meson-meson scattering in ChPT
The scattering of NGB meson is computed in ChPT as an expansion in momenta and meson masses. Denoting as A I (s, t, u) the scattering amplitude of the NGB process a → b with defined isospin I, one has the generic form where s, t and u are the usual Mandelstam variables and A I k = O(p k ), where p means either meson momenta or masses. The LO amplitude A 2 is obtained at tree level from the L 2 Lagrangian. The NLO contribution contains one-loop diagrams from L 2 plus tree-level contribution from L 4 involving LECs.
Using the normalization conventions given in [100,101], the s-channel partial-wave projection of the am- plitude is defined as (13) where N is a normalization factor equal to 2 if all the particles are identical and 1 otherwise. The Mandelstam variables t(s, x) and u(s, x) are defined by the kinematics of the corresponding a → b process and x = cos θ, being θ the scattering angle in the center-of-mass frame.
Being an expansion in momenta and masses, it is clear that ChPT cannot satisfy unitarity, which in the elastic case implies the relation where σ(s) = 2q(s)/s and q is the momentum in the center-of-mass frame. In the following, we only consider the I = J = 1 channel and the superscript index IJ will be suppressed to ease the notation. Nevertheless, ChPT satisfies elastic unitarity perturbatively. For instance, defining as the chiral series of the I = J = 1 ππ partial-wave amplitude, with t 2 (s) and t 4 (s) the tree-level and one-loop ChPT partial-wave amplitudes, in the elastic case one finds the relations Im t 2 (s) = 0, which implies that the unitarity bound in Eq. (14) is increasingly violated in ChPT at larger energy values. In practice, it implies that the chiral series is limited to scattering momenta around 200 MeV above threshold. Furthermore, the ChPT series does not converge equally well in all parts of the low-energy region. This is particularly evident in the scalar-isoscalar channel where strong pion-pion rescattering effects slow the convergence of the chiral series [102]. Finally, at increasingly large momenta, several partial-waves become resonant. Resonances are non-perturbative effects and, as such, they cannot be reproduced within the ChPT power expansion. Furthermore, they usually saturate the unitarity bound in Eq. (14), which implies that elastic unitarity can be violated in the resonance region.

C. Unitarity and analyticity
Below the first production threshold, located at s = 16m 2 π , ππ scattering is purely elastic and, consequently, it can be described in terms of its phase shift. Above this energy, there are possible intermediate processes such as 2π → 2π n, with n = 2, 3, . . . or ππ →KK, ηη, which, in principle, have to be taken into account. In our case of interest, the P -wave ππ-scattering partial wave, inelasticities are completely negligible below the KK threshold and very small below 1.4 GeV [6,8,[103][104][105][106][107]. Thus, in this work elastic scattering is assumed to occur below the KK threshold and above only the ππ and KK channels are considered.
The unitarity condition for the S-matrix, SS † = 1, implies that, for two-coupled channels, it can be parameterized in terms of only three independent parameters. It is customary to choose them as the ππ → ππ and KK → KK phase shifts, denoted as as δ 1 and δ 2 , respectively, and the inelasticity η. Thus, the S-matrix is expressed as The T -matrix elements t ij of the scattering amplitude are related to S-matrix elements as, with and i, j = 1, 2. The relation between the S-and Tmatrix, Eq. (18), allows one to derive the following unitarity condition for the T -matrix elements or in matrix form, being Eq. (21) implies the coupled-channel unitarity relation is fulfilled. The phase space definition, Eq. (19), ensures that in the elastic case, i.e., below the KK threshold, elastic unitarity is satisfied. In the one channel case, Eq. (23) simplifies to The unitarity conditions in Eqs. (23) and (24) imply that the inverse of the imaginary part of an scattering amplitude in the physical region is completely fixed by unitarity. The strong relation between unitarity and resonances has motivated the development of several ChPT inspired methods based on imposing exact unitarity. Some of them are the so-called K-matrix method [108] and the chiral unitarity approach. The latter was considered first in [46,109] to describe ππ and KK scattering in the scalar-isoscalar channel, leading to fairly precise determinations of the f 0 (500) and f 0 (980) resonance properties. There are also more involved unitarization methods. For example, the Bethe-Salpeter (BS) equations were solved for ππ scattering in Refs. [51,52], both in the on-shell and off-shell schemes, while the N/D method was employed in Ref. [110] providing also results for the rest of lightest scalars, namely the κ(700) and a 0 (980). However, none of them generates the ρ(770) pole in the ππ scattering P wave.
The energy-dependence of an scattering amplitude is also strongly constrained by analyticity. Analyticity is based on the Mandelstam hypothesis [111], i.e., the assumption that an scattering amplitude is represented by a complex function that presents no further singularities than those required by general principles such as unitarity and crossing symmetry. In this way, poles in the real axis are associated with bound states (absent in lowenergy meson-meson scattering) and production thresholds give rise to cuts. Cuts are a consequence of the unitarity condition given in Eq. (21), which, together with the Schwartz-reflection principle, imply that an scattering amplitude must have a cut where unitarity demands its imaginary part to be non-zero. It occurs due to both, direct and crossed channels, leading to a right-(RHC) and left-hand cut (LHC), respectively.
Once analyticity is established, Cauchy's integral formula allows one to construct a representation that relates the amplitude at an arbitrary point in the complex plane to an integral over its imaginary part along the right-and left-hand cuts, the so called dispersion relations. The convergence of the dispersive integral often requires subtractions, which introduce a certain number of a priori undetermined constants. The Froissart-Martin bound [112,113] guarantees that at most two subtractions are needed to ensure the convergence at infinity, but one subtraction is enough for the ππ scattering amplitude in the vector-isovector channel. Thus, a once-subtracted dispersion relation for I = J = 1 ππ scattering reads where the first and second integrals stand for the RHC and LHC contributions, respectively. The subtraction constants involve the evaluation of the amplitude at s = 0, so that, they can be pinned down by matching to ChPT in the regime where the chiral expansion is expected to show better convergence properties. However, while the value of Im t(s) in the physical RHC is constrained from unitarity, the LHC contribution is in principle unknown. On the one hand, most UChPT methods differ in the way the LHC is treated. While the K-matrix and chiral unitarity approach models simply neglect the LHC contribution, the BS and N/D methods approximate it with ChPT. On the other hand, Roy-Steiner equations [114,115] solve this problem exactly using crossing symmetry. They provide a representation involving only the physical region, but which, at the same time, intertwines all partial-waves with different isospin and angular momentum. Although, Roy-Steiner-equation solutions allow for high-precision descriptions of different scattering processes at low energies [4-6, 116, 117], and provide the proper framework to extract resonance pole parameters [7,[118][119][120][121][122], or to evaluate an scattering amplitude in an unphysical region [123][124][125], their analysis requires experimental information for the high-energy contribution and higher partial waves. Thus, they are in principle inappropriate for the analysis of lattice data at different quark masses.
In this article, we follow the IAM, which will be outlined in the next section II D.

D. Elastic Inverse Amplitude Method
The Inverse Amplitude Method exploits the relation between a dispersion relation for the inverse of an scattering amplitude and the ChPT amplitude at a given order. At NLO in the chiral expansion, taking into account that ChPT amplitudes grow as s 2 when s → ∞, one needs three subtractions to ensure the convergence at high energies. Thus, a thrice-subtracted dispersion relation for a elastic ChPT ππ-scattering partial wave reads t 2 (s) =t 2 (0) + t 2 (0)s, where we have used Eq. (16) to fix the absorptive part of t 4 (s) in the physical region. Note that Eq. (26) is strongly related to a thrice subtracted dispersion relation for the function g(s) = t 2 (s) 2 /t(s), so that in an elastic approximation the RHC contribution coincides exactly with that of −t 4 (s). The subtraction constants require the evaluation of the scattering amplitude and its derivatives at s = 0, the kinematic region where ChPT provides a reliable description. Thus, using ChPT at NLO one gets Being the RHC exactly fixed from unitarity, and once the subtraction constants are estimated using ChPT, the only remaining unknown information in Eq. (27) is the LHC. The left-hand cut might indeed play a relevant role below threshold, but it is expected that its contribution should be less important as one moves into the physical region. Thus, for a qualitative description it is sufficient to approximate the left-hand cut using ChPT. At NLO, one finds Im g(s) t 2 (s) 2 Im 1 t 2 (s) + t 4 (s) −Im t 4 (s). (29) Inserting Eqs. (28) and (29) in Eq. 27 one obtains which stands for the well-known equation of the IAM method. The IAM was derived first in [47,48] using only unitarity for ππ scattering. Its derivation from a dispersion relation and application thereafter to πK scattering was investigated in [49,50], whereas the remaining IAM meson-meson scattering processes were studied in [85] to one loop. The two-loop version of the IAM was derived in [126] and its generalization to include the effect of Adler zeros was obtained in [127].
The IAM provides a simple algebraic equation that ensures elastic unitarity while at low energies reproduces the chiral expansion. This fact implies that the IAM can be used to describe the resonance region below 1 GeV, i.e., well beyond the applicability range of ChPT. Furthermore, it is based on a dispersion relation, hence, its use in the complex plane is justified, providing a simple tool to study resonance properties. The main difference between the IAM and the on-shell BS or N/D method, is that, in the IAM only the absorptive part of the left-hand contribution is expanded at low energies. It implies that the left-hand cut energy dependence is still controlled by a dispersion relation instead of being fully given by ChPT. In addition, the IAM generates not only scalar but also vector resonances [46], without involving new additional parameters rather than the ChPT LECs. Hence, it reproduces at low energies the quark mass dependence predicted by ChPT.
Nevertheless, it has also several caveats. While the RHC is solved exactly using elastic unitarity, the LHC is approximated using ChPT. The direct consequence of this fact is that the IAM breaks crossing symmetry. Besides, while the IAM provides higher order ChPT contributions needed to fulfill unitarity, some of the leading order logarithms from higher-order loop graphs appear with the wrong coefficients [128].
In addition, it is worth mentioning that the IAM describes experimental data, including resonance pole parameters, of meson-meson scattering in the region below 1 GeV only within a 10%-15% accuracy [85,129]. This small difference highlights the relevance of the LHC in the physical region below 1 GeV.
Clearly, leaving the LECs as free parameters to be adjusted to data instead of being fixed to the ChPT values improves the description of the experimental data. Indeed, ππ and πK scattering experimental data were described in [83,85] using the IAM with LEC values compatible with pure ChPT determinations. Small LECs changes are indeed expected since the IAM includes contributions that go beyond the pure chiral expansion at a given order. However, it is important to remark that while ChPT is a natural theory in the sense that its predictions are linear in LECs changes, the IAM as well as other UChPT models are strongly dependent on precise LECs determinations. Small changes on the LEC values might produce large effects on the phase-shift and pole parameter predictions.
Finally, let us remark that the dispersive derivation of the IAM only constrains its energy dependence, and hence, it is not clear whether it provides the correct quark-mass dependence. While the IAM reproduces the ChPT series at low energies, thus, ensuring that it provides the quark-mass dependence predicted from QCD close to the chiral limit, it also introduces higher-order contributions that spoil the chiral series at higher energies and for heavier quark masses. Thus, high quality lattice data for different light-and strange-quark masses are key to ensure that the chiral extrapolation performed within the IAM is well consistent with QCD.

E. Coupled channel formalism
The generalization of the inverse amplitude method to coupled channels should be in principle straightforward if one assumes the factorization of the RHC and LHC contribution for the different channels involved. In this case, we can define the matrix version of the function g(s) in Eq. (27) (22)). Similarly as in Eq. (27), a thrice-subtracted dispersion relation for G(s) reads where s th and s L stand for the corresponding right-and left-hand cut branching points, respectively. The numerator of the RHC contribution corresponds to the matrix version of Eq. (16), i.e., and hence, the right-hand cut of G(s) coincides with that of the matrix −T 4 (s). The subtraction constants can be evaluated using ChPT. By means of expanding T −1 as one recovers the equivalent version of Eq. (28) in matrix form. However, the problem now is the evaluation of the left-hand cut. Although the RHC branching point s th = 4m 2 π is common for all the elements of the T-matrix, the LHCs of the various channels do differ. Namely, while the ππ scattering LHC starts at s = 0, the LHC for the KK → KK partial wave opens at s = 4m 2 K − 4m 2 π . In this way, proceeding as we did for the elastic IAM, i.e., taking the perturbative expansion in Eq. (33) for the absorptive part of G(s) along the LHC, one is indeed mixing the LHCs of all T-matrix elements. This translates into a violation of the factorization hypothesis, which produces spurious left-hand cuts breaking unitarity [41,85,[130][131][132]. As a summary, the analogous of Eq. (30) cannot be derived in coupled-channels using a dispersion relation.
Alternatively, one can still exploit unitarity in order to derive a coupled channel version of Eq. (30) valid in the real axis. Taking into account Eq. (23), Im T −1 = −Σ, one can write Now, Re T −1 can be approximated once more with ChPT. Using Eq. (33) one gets which provides the IAM coupled channel unitarization formula. Note that to derive Eq. (35) we have used Eq. (32). Nevertheless, it is important to note that Eq. (35) is only justified in the real axis where the ChPT coupled channel unitarity relation (32) is fulfilled. At this point, it is also important to discuss at which energy the couple-channel formalism should be taken into account. Given the phase-space definition in Eqs. (19) and (22), the unitarity relation in Eq. (21) acquires dimension two only when one crosses the KK production threshold. Thus, Eq. (35) should be used only above the KK threshold, i.e., when its dimension coincides with the number of states accessible and the coupled-channel unitarity relation in Eq. (23) is fulfilled. Below this energy one should consider the one-dimensional IAM equation. Thus, this procedure yields a discontinuity at 4m 2 K , instead of a single continuous function. Alternatively, one can include the KK channel for all energies. This provides a continuous function but it again introduces spurious left-hand cuts, leading to a violation of unitarity. Nevertheless, these violations are in general small, around 2%-5% [85,132]. In this paper we consider the second approach for Eq. (35), but in order to reduce the effect of spurious cuts, we introduce an extra term in the χ 2 of our fit to lattice data, which penalizes unitarity violations of the S-matrix by some factor, as explained in Sect. II H.
Eq. (35) was used in [85] to study all possible amplitudes for meson-meson scattering leading to a fairly good description of all available experimental data below 1.2 GeV with reasonable LEC values. These amplitudes were analytically continued to the complex plane in order to look for poles associated to the lightest scalar and vector resonances [129,133], with determinations compatible with experimental values within uncertainties. This result suggests that the role of spurious LHCs which prevent the dispersive derivation of the coupledchannel IAM formula are also small. Furthermore, we have explicitly checked that by removing the t-and uchannel loop functions ( Fig. 1) in the ππ →KK and KK →KK ChPT amplitudes that generate the spurious cuts, the mass and width of the ρ-meson obtained in the global Fit IV (see Sect. IV C) change less than 1 and 6 MeV, respectively, i.e., within the uncertainties quoted. Nevertheless, the effect of the t and u channels in the ππ amplitude lead to a shift of 6 and 15 MeV for the mass and width of the ρ-meson in Fit IV, respectively (without readjusting the LECs).
To conclude, Eq. (35) is the tool we use to analyze lattice scattering data in the ρ(770) channel. The explicit expressions for the elements of the T 2 and T 4 for ππ → ππ, ππ → KK and KK → KK are given in the appendix of [85].

F. Resonances
Resonances are formally defined as poles lying on unphysical Riemann sheets. An unphysical Riemann sheet is reached when the physical right-hand cut is crossed continuously from the upper-half plane to the lower-half plane above a given production threshold. In the elastic scattering case, there are only two Riemann sheets, the physical and unphysical one, which are called, first and second sheet, respectively. These two Riemann sheets must coincide in the real axis, In addition, the scattering amplitude on the first Riemann-sheet satisfies the Schwartz reflection principle, i.e., S (s + i ) = S * (s−i ), which together with unitarity, SS * = 1, yields the relation The analytic continuation of Eq. (37) into the complex plane implies that a pole on the second Riemann sheet corresponds to a zero in the physical one. By means of Eq. (18) one can translate this relation to the T -matrix elements, leading to where σ(s) = 1 − 4m 2 /s, and its determination is chosen as σ(s * ) = −σ(s) * , to ensure the Schwartz reflection symmetry. When further channels are opened, more unphysical Riemann sheets can be defined by continuing the square momenta of the intermediate states over the different thresholds. Thus, there are 2 n Riemann-sheets for a given number n of opened channels. The generalization of Eq. (38) in a coupled-channel formalism is straightforward where Σ (n) is a diagonal matrix containing the phase space factors of those channels that have been crossed continuously. In particular, for the ππ and KK I = J = 1 coupled-channel case, we will have four different Riemann sheets defined as where σ π = 1 − 4m 2 π /s and σ K = 1 − 4m 2 K are the phase space factors of the ππ an KK channels, respectively.
Therefore, a pole in the T matrix corresponds to a zero of the determinant of the matrix inside the brackets of Eq. (39), which is denoted by where M and Γ stand for the mass and width of the resonance, respectively. In addition, the dynamics of a resonance is strongly related to its coupling to a given channel, which is defined from the pole residue as where p(s) stands for the center-of-mass-system momentum of the corresponding process.

G. Formalism in the finite volume
The Lüscher's approach [134,135] allows one to relate the measured discrete value of the energy in a finite volume to the scattering phase shift at the same energy in the continuum. The volume-dependence of the discrete spectrum of the lattice QCD gives the energy dependence of the scattering phase shift. This method, originally derived for a single scattering process was soon extended to coupled channels for potential scattering [136], nonrelativistic effective theories [137,138] and to relativistic scattering [139][140][141][142]. Extensions of the Lüscher formalism to three-particle systems under certain conditions are also available [143][144][145][146].
The Lüscher's approach is based on the analysis of the dominant power-law volume dependence that enters through the momentum sums in a BS equation, where all quantities are written in terms of non-perturbative correlation functions. In order to extract this dependence one assumes that the BS kernel, which accounts for the LHC and subtraction constant contributions in Eq. (25) and only involves a exponentially suppressed dependence on the volume [134], coincides for large volumes with its infinite-volume form. In this way, the difference between finite-and infinite-volume integrals entering on the BS equations only depends on on-shell values of the twoparticle integrand leading to the the quantization condition where T is the scattering amplitude in the finite volume and F is a matrix that contains sums of the generalized Zeta functions subduced into the relevant finite volume little groups [139,140]. Lüscher's method was subsequently rederived in [147,148] by discretizing the s-channel loop functions which appear in the IAM coupled-channel equation of Eq. (35) and neglecting the t-and u-channel contributions. The consideration of relativistic effects in the Lüscher formulation as in Ref. [149] can lead to a significant difference for small volumen sizes in the phase shifts extracted. The discretization of the t and u channels has been discussed in [150,151]. In the latter, the exponentially suppressed volume dependence of the LHC contribution was explicitly taken into account, concluding that the LHC volume dependence is numerically negligible for lattice sizes L > 2m −1 π while for lattice volumes m −1 π < L < 2 m −1 π , it only affects noticeably the first energy level. Furthermore, note that neglecting the volume dependence of the LHC contribution in the finite volume is by no means equivalent to ignoring the LHC in the continuum; lattice energy levels are non-perturbative quantities and, as such, they include all physical effects, both from the RHC and LHC contributions. The same cannot be stated for the dispersive formalism defined in Sect. II D and II E since one explicitly factorizes the RHC and LHC contributions. However, to extract information from the energy levels and connect them with the T-matrix in the continuum one does need a generalized Lüscher method including all physical effects, which might become particularly difficult, for example, in the case of multi-channel and intermediate states of three or more particles.
In principle, one could use the formalism in [150,151] to evaluate the energy levels and fit them to the lattice data. Nevertheless, in order to avoid the discretization of loops we follow here the method used in [84]. Namely, we fit the phase shift values extracted from the lattice using Lüscher's method, while the eigenenergies are reconstructed by means of a Taylor expansion taking into account the correlation between energy E n and phase shift δ(E n ), as well as the covariance matrix of eigenenergies provided by the lattice. This method is explained in the subsection below.
We analyze lattice simulations on two different chiral trajectories, where either the sum of the three-lightest quarks or the strange-quark mass is fixed to the physical point, i.e., TrM = C or m s = k, respectively. The corresponding tree-level pseudoscalar meson masses relations are for TrM = C and for m s = k, with k = m 0 s or 0.6 m 0 s . As a result from a combined analysis of data on these two kind of trajectories, in Sect. IV C we also show predictions for ρ-meson phase shifts, pseudoscalar meson decay constants and masses in other trajectories where the strange-quark mass is fixed to values smaller than the physical one, m s = k with k < m 0 s , on the SU(3) symmetric trajectory, m s = m ud , i.e., and for trajectories where the light-quark mass is kept fixed at the physical point m ud = m 0 ud , i.e., We employ one-loop ChPT for the analysis of pseudoscalar meson masses and decay constants, see Sect. III, in combination with the coupled-channel IAM discussed in Sect. II E for the ρ-meson phase shifts. The fitting parameters are the LECs entering into our expressions, i.e., L i , with i = {3, . . . , 8}, L 12 = 2 L 1 − L 2 , and the parameters which fix the chiral trajectories in Eqs. (43) and (44), C B 0 and k B 0 . The chiral scale µ is fixed to 770 MeV and the pion decay constant in the chiral limit f 0 is set to 80 MeV. We fixed f 0 because its inclusion as a new fitting parameter did not entail any substantial reduction of the χ 2 . In the following we describe the contributions to the χ 2 .
Meson-meson scattering in the lattice translates into discrete energies which are correlated. In order to take into account those the following function is minimized, where E is the vector of eigenenergies measured on the lattice, C their covariance matrix and E the corresponding energies of the fit function. Nevertheless, we do not fit directly lattice energy levels but phase shifts extracted using the Lüscher formula. In order to take into account the energy correlations, we follow the method considered in [84]. This is, for each energy level , E i , a Taylor expansion of both, the phase shift extracted from the lattice, δ L , and the one evaluated in the IAM, δ IAM , is performed around the energy given by the lattice simulation, E i . If one assumes that both δ L and δ IAM coincide exactly at E i , at leading order, one finds which provides a direct way to evaluate χ 2 E in Eq. (47) in terms of phase shift values.
Regarding pseudoscalar meson masses and decay constants from the lattice, we fit the ratios, h 0 = m K /m π , h 1 = m π /f π , h 2 = m K /f K and h 3 = m K /f π , which are, in principle, more stable against possible discretization effects. Thus, we also minimize where i denotes the different ratios, j = 1, · · · , n are the measurements, and n is the length of lattice data. The superscripts l and p indicate values from lattice simulations and predicted by one-loop ChPT, respectively. Finally, as already discussed in section II E, the coupled-channel version of the IAM generates unphysical LHC contributions arising from the on-shell coupledchannel approximation considered. These contributions produce small violations of unitarity, which translate into undesirable phase shift picks at low energies and in the resonance region, starting below s = 4 m 2 K − 4 m 2 π (this energy corresponds to 880 MeV for the HadSpec lighter pion mass). These small picks are enhanced when there are lattice data around that energy. To eliminate these unphysical artifacts, a term that minimizes S-matrix unitarity violations at a degree controlled by a parameter λ is added to the χ 2 , In summary, the total χ 2 -like minimization function reads as In Fig. 2 we show the value of χ 2 and χ 2 λ in Eq. (50) as a function of λ for the minimization of the Hadron Spectrum Collaboration ρ-meson phase-shift data at m π = 236 MeV [60] together with decay constants from MILC [72]. The LEC values obtained are given in Fig. 3. Clearly, for λ ∼ 40 the LECs become stable while χ 2 λ /λ gets significantly reduced. One could also choose a higher value of λ, however, at the cost of increasing χ 2 . Thus, we set the value of λ to 40 .
There is an additional caveat that one should take into account; ChPT is built as an expansion in meson masses and, as such, the chiral series is only expected to converge for light pions. In order to study the convergence radius of ChPT we perform first individual fits of lattice data sets and discard pion mass results for which the fit does not pass the Pearson's χ 2 test at a 90% upper confidence limit. This restricts the lattice data sets to pion masses below around 430 MeV. Results presented in the following sections beyond that pion mass are merely qualitative.

III. CHPT: DECAY CONSTANT ANALYSIS
In this section, we attempt to perform a global fit of pseudoscalar meson masses and decay constants {m π , m K , f π , f K } from [70][71][72][73][74][75][76]. These data are simu-  A few aspects need to be considered before. First, the role of the renormalization scheme used in the lattice simulations to fix quark masses. Here, we do not adjust quark masses values but pseudoscalar meson masses, which, in principle, should be independent of the renormalization scheme. Still, we checked if the pseudoscalar meson masses in the lattice data sets with different renormalization schemes are compatible. For example, we notice that UKQCD Collaboration uses the MS scheme at 3 GeV [71], while the MILC Collaboration uses the same scheme at 2 GeV [72][73][74]. When we compare both sets of data, we do not observe any substantial inconsistency, but instead, their values do agree quite well.
Second, other important issue is the size of the pion masses used in the simulations. We observe that in general the JL/TWQCD [75] and PACS-CS Collaborations [76] have larger pion and kaon masses. For instance, the JL/TWQCD pion and kaon masses are larger than 300 and 600 MeV, respectively. These values might be too large for the perturbative ChPT expansion and indeed we are not able to fit these data sets in combination with MILC and UKQCD data. Thus, in this fit we only include data from [70][71][72][73][74]. The JL/TWQCD and PACS-CS data are studied in separated analysis in the next section.
Third, we should discuss possible finite volume and lattice spacing effects. In [70][71][72][73][74], the dependence of the decay constant determinations with the lattice spacing was studied carefully and the results were extrapolated to the continuum. These extrapolated data are the input of the fit we show here. Another difficulty that we find to study data from [76] (PACS-CS) is the following. In [76], the chiral trajectory is set in such a way that the physical point of the strange quark is determined and later fixed onto the chiral trajectory of the simulation. Thus, the m K dependence on m π in principle should agree with that from MILC [72][73][74] and UKQCD [71], since these simulations are also performed at the physical strangequark mass. However, we found substantial discrepancies in the behavior of the chiral trajectory in [76] with those from MILC and UKQCD, which may be due to finite volume effects. These discretization effects can be partly absorbed by the free parameters. The result from analyzing PACS-CS data, pseudoscalar meson masses and decay constants [76] together with ρ-meson phase-shift data [63], is shown in the next section.
Lastly, in [70] (CLS collaboration) two different scale settings are provided, called here scale settings A and B. In the first one, scale setting A, the lattice spacing is determined by fixing the chiral extrapolations of f π and f K to the physical point. The second one, scale setting B, uses the Wilson flow (t 0 ) to set the scale by assuming that, for all different ensembles, the data over the TrM = C trajectory intersects the m ud = m s symmetric line at φ 4 = 1.15, with φ 4 = 8 t 0 m 2 K + 1 2 m 2 π . This method requires small corrections in the quark masses from the ones used in the simulations [70], which translates into small shifts for the pseudoscalar masses and decay constants. Nevertheless, the CLS ρ-meson phase-shift data in [79] for the scale setting B were not shifted accordingly, and hence, these corrections could lead to a conflict among the CLS decay constant and phase shift data. Then, we take here the no-shifted values, first rows of Table II in [70]. For each ensemble β, these two scale settings lead to different lattice spacing values a β . Namely, {a 3.4 , a 3.46 , a 3.55 , a 3.7 } = {0.079, 0.071, 0.061, 0.0481} for scale setting A and {a 3.4 , a 3.46 , a 3.55 , a 3.7 } = {0.086, 0.076, 0.064, 0.0498} for B [70]. Nevertheless, we find that scale setting B produces systematically smaller values of f π than A for the same pion masses. For instance, we see a difference of around 4 MeV in f π for pion masses of around 200 − 300 MeV between the two scale settings. This difference is not small, since changes of 80 MeV in m π imply variations on f π of around 4 MeV in these data. Because of these discrepancies, we are only able to find an optimal χ 2 when data with scale setting A are included. Notice that this is the method that fixes the scale using the f π and f K physical quantities. 6 In section IV B we analyze the decay constant data in combination with ρ-meson phase-shift data for both scale settings and discuss the main differences.    In conclusion, it is only possible to do a combined fit of data from [70] (scale setting A) and [71][72][73][74]. 7 In Tables I and II, the values of the fitting parameters obtained from this analysis are presented. This result is called Fit I. We notice that the LECs in this fit are not very sensitive to small variations of the CB 0 and kB 0 parameters, being thus quite stable. Furthermore, they are in line with the compilation of the FLAG Review [78], which only includes results for m s = k data. However, note that we are obtaining much smaller LEC errors compared to the FLAG average. Notice also that since these data include variations of the strange-quark mass, one is able to fix well the strange-quark mass dependence of the pseudoscalar decay constants for the pion masses studied.
The various chiral trajectories studied are shown in Fig. 4 (top-left panel), where one can see that the kaon mass squared data for the TrM = C trajectory [70] differ considerably from the m s = k ones [71][72][73][74]. In addition, two of the ensembles simulated in [70], the ones with β = 3.55 and 3.7, lead to very similar curves and hence to similar values of CB 0 in Table II. Furthermore, the UKQCD [71], MILC [72,73] and Laiho [74] lattice data are in very good agreement. Indeed, ChPT is able to reproduce well the data on these two different trajectories.
The ratios m π /f π , m K /f π and m K /f K are also depicted in Fig. 4. For the ratio m π /f π , it is worth noting that all data, independently of the chiral trajectory, lie almost on the same curve. This suggests that the ratio m π /f π is indeed quite independent on m s . We discuss this further in Sect. IV C. In fact, all lattice data for this ratio fall into the gray error band plotted, which is just an extrapolation of the percentage error of this ratio at the physical point determined by MILC [73]. For this collaboration only m π , m K and f π data are provided, which are shown with dashed black lines. The m s = 0.6 m 0 s trajectory from [72] is denoted by a solid gray line. The data of Laiho [74] is represented by dashed-orange lines. UKQCD data are denoted by black squares, while CLS data [70] are given by dark-green squares (β = 3.4), blue circles (β = 3.55) and yellow pentagons (β = 3.7). Note that the UKQCD Collaboration and Laiho data sets provide very similar values of f K . In addition, we include in Fig. 4 the chiral prediction for the SU(3) m s = m ud trajectory. Both, m s = 0.6 m 0 s and m s =m trajectories, lead to a substantial reduction of the ratios m K /f π and m K /f K . For pion masses larger than 430 MeV, the ChPT prediction begins to differ from the data, which suggests the breakdown of the chiral series. In this section we analyze the ρ-meson phase-shift data from [60][61][62][63][64] and pseudoscalar meson masses and decay constants from [71][72][73][74][75]. All these data are taken from simulations over chiral trajectories where the strangequark mass is kept fixed to the physical value, m s = m 0 s , except for the JL/TWQCD, where k {1.6, 2} m 0 s [75]. In fact, the pion and kaon masses used in the simulations of [75] are larger than in the other simulations. This simulation is studied independently and discussed at the end on in this section.
First of all, we perform individual fits to the pseudoscalar masses and decay constant ratios from UKQCD [71], MILC [72,73] and Laiho [74] together with the ρ-meson phase shift data from the HadSpec Collaboration [60,61] corresponding to m π = {236, 391} MeV. The LECs obtained in these fits are shown in the second, third and fourth columns of Table III, respectively. Al-though some small differences among the individual fits are observed for L 5 and L 6 , they provide in general compatible LEC values within uncertainties. Thus, we conduct a simultaneous analysis of the UKQCD, MILC and Laiho decay constants and HadSpec phase shifts, which is denoted as MUL+HS in the fifth column of Table III. As expected, the fit provides a good description of all data with consistent LECs.
Finally, we include the phase-shift results from [62] (JB) at m π = 233 MeV. This is denoted as Fit II in the sixth column of Table III. Notice that this fit encompasses a large bunch of data on m s = k (k = {1, 0.6} m 0 s ). The LECs obtained in these fits are very similar to the previous ones suggesting consistency among the different data sets. Results for ρ phase shifts together with the fitted lattice data are plotted in Fig. 5. As shown in Fig. 5 (left, top), the extrapolation of Fit II results to the physical point (light-blue solid line) is very close to experimental data, depicted as light-blue squares [103] and orange circles [106].
Regarding decay constant ratios and pseudoscalar meson masses, results from Fit II are very similar to those obtained in Sect. III over m s = k trajectories, and are shown in Fig. 38 in the Appendix VI B.
Unfortunately, we could not obtain additional consistency with the lattice data from Refs. [63,64,75]. Thus, in the following, we analyze the remaining lattice results separately. The simulation of [64] (CA) for ρ-meson phase-shift data does not include decay constant determinations, thus, we analyze this data with the UKQCD meson and decay constant values. If other decay constant data are used instead, as for example, those from MILC, the results are very similar. The resulting LECs, given in the second column of Table IV, are, in general, compatible with the values from Fit II, but we find slightly different values for L 4 and L 5 , and larger discrepancies for L 8 . These differences have a large impact on the phase-shift values. As shown in the right-top panel in Fig. 5, the extrapolation to the physical point provides results incompatible with experimental data.
Concerning the JL/TWQCD collaboration decay constant data [75], we only find good partial fits if we include the three and two lightest pion mass data points for the trajectories m s = 1.6 m 0 s and m s = 2 m 0 s , respectively. This can be due to the breakdown of the ChPT expansion for such large m s values. Since in these simulations decay constant determinations are provided but not ρ phase shifts, we analyze them together with the Hadron Spectrum Collaboration (HS) ρ-meson phase shift results at m π = {236, 390} MeV. The only purpose of this fit is to show the qualitative behavior of the pseudoscalar meson mass and decay constant ratios over trajectories with larger m s values than the physical one. The corresponding LECs obtained in the fit are given in the third column in Table IV. A comparison with the result from Fit II in Table III shows up sizable discrepancies between both fits, which might be due to inconsistencies of the JL/TWQCD data with data included in Fit II apart from the breaking of the chiral series. These phase shift results are also plotted in the top-left panel of Fig. 5 in dashed lines. Nevertheless, the extrapolation to the physical point of this fit turns out to be also very close to the experimental data.
For the PACS-CS collaboration, both ρ-meson phase shift [63] and decay constant [76] data are available and analyzed together. The LECs are given in the fourth column of Table IV, and also the L i 's, i = 4, 6 − 8, differ considerably from the Fit II values. As explained in Sect. III these data have larger kaon masses for the same trajectory m s = m 0 s than data in Fit II. This can be due to sizable finite volume effects in these simulations. As a consequence, these data are in disagreement with the data included in Fit II. In this case, the extrapolation to the physical point, depicted in the bottom-right panel of Fig. 5, fails substantially to describe the experimental data.
Let us note that these ρ-meson phase-shift data were analyzed before in [84] using the UChPT model in [46]. Even though this model neglects the LHC contribution, which now is taken into account, we obtain here similar results to the ones of [84].
The chiral trajectories and decay constant ratios for these fits are shown in Fig. 38 of the Appendix VI B. We find that the pseudoscalar meson mass data on m s = k trajectories fit very well into a linear formula m 2 K = a m 2 π + b with slope a = 0.5, depicted in dotted lines. This behavior is qualitatively similar to the leading order ChPT prediction. For the ratios of decay constants we find similarities with the results of Fit I over the trajectory m s = m 0 s . The ratios m K /f π and m K /f K in other m s = k trajectories as a function of the pion mass are parallel to the ones over m s = m 0 s and take higher values. For the m π /f π ratio, only the JL/TWQCD and PACS-CS data are just a bit out of the error band.
Finally, the ρ(770) pole position on the second Riemann sheet obtained in the different fits are given in Table V. While the values obtained in Fit II and in the JL/TWQCD & HS fits are compatible with the most precise theoretical prediction [7], the results obtained for UKQCD & CA and PACS-CS provide smaller and larger values, respectively. In order to write down results that can be compared to the BW values provided in lattice articles, we perform a refit of the IAM solution to the BW formula in Eq. (60). As we show in Fig. 39 in the Appendix VI B, the data is also well described by a Breit-Wigner (BW) parameterization. The BW mass, coupling and width, normalized to the pion mass, are shown in Table V, where we also provide the result for the extrapolation to the physical point.

B. Chiral trajectories TrM = C
In this section we show the outcome of the analysis of ρ-meson phase-shift [79] and decay constant [70] data of the CLS Collaboration over trajectories where TrM = 2m ud + m s = C. Thus, in these trajectories the kaon becomes lighter as the pion mass increases. Two different scale setting methods were considered in [70]. These two methods lead to differences of around 10 − 20 MeV in m π , 20 − 45 MeV in m K , and 4 − 8 MeV for f π and f K . These differences entail several difficulties. As discussed in section III, we could only find an optimal solution to the minimization problem of Fit I, that also includes m s = k data, when scale setting A was taken for the pseudoscalar meson masses and decay constants over the TrM = C trajectories. When the scale setting B was considered instead, the global χ 2 minimum was found to be around twice larger than with the scale setting A. On the contrary, we observe that, when using scale setting A, the dependence of ρ-phase shift data with the pion mass of [79] cannot be described well within the IAM for all ensembles. While the ensembles D101, J303 and D200 are well described, the ensemble N200 (or N401) cannot be reproduced. This is because the IAM predicts higher values of the ρ meson mass for the pion mass used in this ensemble, see    describing all TrM = C lattice data, i.e., pseudoscalar meson mass and decay constant ratios and ρ-meson phase shift (excluding m s = k data). These phase-shift results are plotted in Fig. 41 (Fit IIIB) in the Appendix VI B.
Nevertheless, it is possible to perform fits of decay constant and ρ-meson phase shift data for ensembles with the same gauge coupling β [79]. Namely, C101, D101 (β = 3.4), N401 (β = 3.46), N200, D200 (β = 3.55) and J303 (β = 3.7). Several of these ensembles use the same pion mass but different volume or lattice spacing. On one side, the ensembles C101 and D101 are simulated with the largest lattice spacing but D101 uses a volume 2.4 times bigger than C101. On the other side, the ensembles N200 and N401 were simulated in the smallest    , g BW and Γ BW obtained from the refit of the IAM solution to the BW formula in Eq. (60) normalized to the pion mass, i.e., E0 stands for E0/mπ. volume but N200 has a lattice spacing 1.13 times smaller than N401. Finally, J303 has the biggest volume and smallest lattice spacing. 8 In this way, possible differences between individual fits in these pairs might highlight finite volume and lattice spacing effects. The resulting LECs are shown in Tables VI and VII for the A and B , respectively. The ensembles C101 and D101 are fitted separately in order to study the finite volume effect. Overall, the values of the LECs L 3 and L 5 are approximately stable, but we find large differences for the others.
We also attempt to perform combined fits including most ensembles for different β in order to check whether these effects can be absorbed in the LECs. Since the D101 and N200 ensembles supersede the C101 and N401 ones, accordingly, we only include the ensembles D101, N200, D200 and J303. These fits are denoted as Fit III A and III B for the A and B scale settings, respectively, and they also include the corresponding pseudoscalar meson mass and decay constant ratio data. The LECs obtained are given in the last columns of Tables VI and VII.
In general, the LECs of Fit IIIA agree better with those obtained for the m s = k trajectories. The chiral trajectories and decay constant ratios of these fits are depicted in Fig. 38 (right) in the Appendix VI B, where we also show the result of the m s = k fits for comparison (left panel). Results for the fits IIIA and B are plotted in like blue-solid and green-double-dot-dash lines, respectively. The kaon mass dependence on the pion mass for the trajectories TrM = C also fit well into straight lines, m 2 K = a m 2 π + b, but now with a slope close to a = −0.5 instead of 0.5, as we found for the m s = k ones. In this way, the IAM is able to reproduce very well the Tr M = C trajectories, which appear as three close decreasing curves intersecting the symmetric line, m s = m ud . At pion masses of 300 MeV, the kaon mass is around 60 MeV lower than for the m s = k trajectory.
The ratios of decay constants are also well reproduced. For the scale setting A, the ratio m π /f π agrees well with the m s = k data, emphasizing that this ratio is almost independent of m s . In the case of scale setting B, it falls a bit out of the m s = m 0 s error band, depicted in a lightbrown color. Note that this behavior is different from the m s = 0.6 m 0 s trajectory (MILC, dotted-gray), which lies inside the error band and does not show any substantial difference with the m s = m 0 s curve. This suggests that there could be small dependencies with the strange-quark mass. We comment more on this issue in the next section.
Results for ρ-phase shifts are provided in Figs. 40 and 41 of the Appendix VI D. As commented before, except for the ensemble N200 in scale setting A, all the other phase shifts can be described qualitatively well in these fits. In Figs. 6 and 7 we show the ρ-meson phase shifts obtained for the different gauge coupling fits. The IAM allows one to describe the ρ-meson phase-shift data in TrM = C trajectories for every ensemble. Nevertheless, note that one can not observe a trend of the overall data indicating that the ρ-meson mass 9 increases monotonically with the pion mass. For scale setting A, the N200 and N401 ensembles give rise to a lighter ρ-meson mass than the ensembles D101 and C101 even when they are simulated with heavier pion masses. In addition, the ρmeson mass takes about the same value for the J303 and D101 ensembles, although the pion mass used in J303 is around 20 MeV larger.
At low energies, phase shifts decrease as the pion mass grows, as expected from the p-wave centrifugal barrier and the chiral expansion. For scale setting B one observes that the trend of the ρ-meson mass dependence on the pion mass is flatter. Noticeably, the ρ-meson becomes lighter for pion masses around 300 MeV in both scale settings. In both cases, systematic effects due to a finite volume and lattice spacing are reflected in around 8 MeV difference in the ρ-meson mass between the C101 and D101 and 14 MeV between the N200 and N401 ensembles, respectively.
The corresponding pole positions and couplings for both scale setting are given in Table VIII. In addition, given the large discrepancies observed between the scale settings we perform a new fit with their average for each gauge coupling β, denoted as Fit C in Table VIII. For comparison, the result of the global fit including both data on m s = k and TrM = C trajectories, discussed in the next section (Fit IV) is also shown. The values are normalized to the pion mass, so that the dependence of the ratio m ρ /m π with the pion mass and scale setting used is visible. Overall, we see that the results for different lattice spacings are quite similar. For the ensemble J303 the dependence on the scale setting considered is negligible, while for other ensembles it produces shifts of less than 1% for the normalized ρ(770) mass and less than 5% for the couplings. Regarding finite volume and lattice effects, the systematic differences between C101 and D101 are of around 1% in the normalized ρ mass, and 1.5% between the N200 and N401, while these are of less than 2% in the couplings in both cases. Finally, the comparison between the individual fit solutions obtained using A and B scale settings is given in Fig. 8, where it can be seen that the differences produced in phase shifts as a function of E/m π are in general reasonably small, and negligible for the J303 ensemble. The largest difference is coming from the size of the lattice spacing used in the simulation, i.e., the difference observed between the N200 and N401 ensembles.
Finally, we can compare the result of Fit C in Table VIII with the result of Fit IV which includes also m s = k data. There are small differences between these two fits of less than 3% in the normalized ρ-meson mass and less than 6% in the couplings. We discuss this further below. In this section we perform a simultaneous analysis of lattice data over both m s = k and TrM = c trajectories. This final study will be denoted as Fit IV and it analyzes lattice ρ-meson phase shift data in N f = 2 + 1 of [60][61][62]79] in combination with pseudoscalar meson masses and decay constants from [70][71][72][73][74]. Thus, this analysis takes into account all data included in the fits II and III of Sects. IV A and IV B.
As discussed in Sects. III and IV B, we were not able to find a solution with the IAM describing data on TrM = C trajectories using neither of the scale settings in [70,79] in combination with data over m s = k. Hence, in order to attempt a global fit some remarks are necessary. First of all, the ensemble C101 of [79] will be discarded since the simulation for the ensemble D101 is performed in a larger volume. Note that even when the N401 ensemble has larger lattice spacing than the N200, the former has more data points and its uncertainties are smaller, therefore, we include both ensembles in the present analysis. Secondly, it is important to highlight that, according to Tables V and VIII, the CLS result for the ratio m ρ /m π of the D101 ensemble is very close to the one from the HadSpec (HS) collaboration at m π = 236 MeV, Fit II; the difference is only of around 2%. This fact points out that the pion masses used in these simulations should also be very similar. Nevertheless, only the average between the scale setting A and B has a similar pion mass (m π = 233 MeV). This facts motivates us to consider that the average between both scale settings provides a reasonable estimate to be used in order to perform a global fit of data. Hence, we perform a bootstrap of the lattice spacing for the TrM = C ensembles assuming that for every β, it is normally distributed around the average of scale setting A and B and the standard deviation being half the difference between them. Not only the lattice spacings, a β 's, but also decay constant ratios and energy levels (normalized respect to the pion mass) for each ensemble are generated from a normal distribution accordingly to their lattice data errors. Regarding the lattice energy levels, the resampling is performed assuming a multivariate normal distribution with the original covariance matrix. 10 Remarkably, following this strategy we could reproduce decay constant and phase-shift data simultaneously on both trajectories. The resampling is performed 300 times and the error is evaluated from that sample. This number of fits turns out to be enough, since the average, median and fit solution (taking the average of the lattice spacing) are indeed very close to each other. Namely, they produce differences in the ρ-meson mass of less than 1 − 2 MeV. Since we are interested on interpreting the results in terms of probability and confidence intervals, our central results are represented by the median or first quartile and the uncertainties will be described by the 68% and 95% confidence intervals (CI), which are represented as darker and lighter error bands, respectively.      In Figs. 9, 10 and 11, the chiral trajectories and pseudoscalar meson mass and decay constant ratios studied are plotted. The lattice data fitted correspond to the extrapolation to the continuum limit with finite volume effects corrected. In more detail, for the m s = m 0 s trajectory we include the UKQCD [71] (purple diamonds), MILC [72,73] (black dashed curves with light-brown error bands 11 ) and Laiho [74] (orange dashed curves and error bands) lattice data. For other m s = k trajectories there is not much data except for the ratio m π /f π extracted by MILC [72] for m s = 0.6 m 0 s (gray dotted line). The TrM = C data from the CLS Collaboration are given for the different lattice gauge cou- 11 The data was sent to us by C. Bernard and the error band is extrapolated from the physical point as suggested by him. plings β = 3.4 (green squares), 3.46 (red circles), 3.55 (blue triangles) and 3.7 (yellow pentagons). The error in the pion mass (x-axis) corresponds to half the difference between the pion mass using the two A and B scale settings. Although in principle chiral trajectories for the several gauge couplings β are different, in practice, we obtain very similar curves when the error in the lattice spacing is considered, which only start to separate more clearly when these cross the symmetric line. This is, we get c (β=3.55) c (β=3.7) and only a small difference for β = 3.4. In addition, we include in Figs π . These three trajectories start at a small value of m s (m s B 0 = 2 MeV 2 ), then, they cross the symmetric line and end up at the m s = m 0 s curve. All ratios m K /m π , m π /f π , m K /f π and m K /f K are reproduced well inside the 95 % CI till m π 400 MeV, when the ChPT predictions start to deviate. Therefore, the predictions for pion masses between m π = 400 − 500 MeV are merely qualitative.
From Fig. 10 one sees that the ratio m π /f π does not depend much on m s . However, this does not necessarily mean that f π is independent of the strange-quark mass. In fact, both m π and f π depend on m s . This dependence is shown explicitly in Figs. 12 and 13. In Fig. 12, the squared leading order mass, M 2 0π , is depicted as a function of the pion mass for different strange-quark masses. Indeed, one can see that m ud kept constant does not imply that the pion mass is constant as well. In fact, our analysis at one-loop level predicts that the pion mass grows with m s for a constant value of m ud ; while for m s = 0 one obtains m K 1/ √ 2m π (see the red line in Fig. 9) and M 0π m π (notice the almost quadratic 12 In this case, it is understood that m d = mu. behavior of the red curve in Fig. 12), consistently with the leading order ChPT prediction, effects coming from the kaon and eta particles in f π become more relevant as m s increases. Although this effect is invisible for very light pion masses and small for physical pion masses (for a constant value of m ud , the difference between the physical pion mass and the m π value at m s = 0 in Fig. 12 is around 14%.), it becomes larger for heavy pions. On the contrary, in Fig. 13, where the dependence of f π on m π for different strange-quark masses is shown, one sees that this dependence is more noticeable for light pion masses and smaller for heavier pion masses, when the uncertainties for f π increase. At the physical point, one indeed finds a difference of 5 MeV in f π between its value at m s = 0 and m s = m 0 s . This is intrinsically connected with the contribution of the terms which involve kaons and etas in Eqs. (5) and (9),  with which are These terms, which involve kaon and eta meson loops (tadpoles) and kaon mass contact terms, called t K,η from now on, contribute slightly for m s = 0, around 1 MeV in the f π value, but they account for about 6 − 7 MeV at m s = m 0 s . Later, we show that this relatively small variation in f π is translated into a visible difference in the ρ mass.
Beyond that, one might wonder what would happen in a world where kaons and etas are not present. In order to answer that question, we could explicitly set to zero the t K,η contribution in Eqs. (5) and (9), and define M π and F π in Eqs. (53) and (54) as the pion mass and decay constant in this world. In such a way, one obtains a dependence of the new pion decay constant, F π , with the pion mass, as the orange-dashed line in Fig. 13. Effectively, this limit is equivalent to set the coupling of pions to kaons and etas to zero, which can be achieved when, first, f K,η in Eqs. (10) and (11) are sent to infinity, 13 and, 13 Note that it can only achieved if one breaks the SU(3) symmetry second, when m K,η are set to zero. This is discussed in Sect. VI A. Furthermore, note that this is different from the so-called SU(2) formalism (corresponding to take the in the partially conserved axial current (PCAC), i.e., one has to differentiate between the pion and the kaon and eta decay constants in the chiral limit. See explanation in Sect. VI A.    Tables VI and VII and for each scale setting. The IAM pole position is denoted by E0, while the BW parameters m BW ρ , g BW and Γ BW are obtained by refitting the IAM solution to the Breit-Wigner formula. The fit C stands for the average of both scale settings, while IV denotes the global fit discussed in section IV C, included for comparison, which is obtained performing a resampling of the lattice spacing and lattice data. The quantities with tilde are normalized to the pion mass, i. e. , E0 stands for E0/mπ. limit m s → ∞), where the effect of pions interacting with kaons and etas is absorbed in the SU(2) LECs, the constant B 0 and f 0 , which are fixed, together with f π , to get observables in the physical world with more than two flavors.
When t K,η = 0, f π reduces its value around 6 − 7 MeV at the physical point, while it approaches the value of the chiral trajectories m s = 0 and m s = m ud as the pion mass decreases, reaching f 0 in the chiral limit. These in- teracting terms involving kaons and etas have a similar effect to the one caused by a change in m s for light pion masses. In fact, varying m s from zero rises the contribution of these terms, increasing the value of f π . Later, we show the contribution of these terms also in the ρ meson mass (using the one-channel ππ IAM equation), but before, we continue discussing results from the analysis with the (ππ − KK) coupled-channel IAM.
The extrapolation to the physical point for the mass and decay constant ratios is given in Table IX The values of the LECs and remaining fit parameters are given in Table X, where errors also represent 68% and 95% CI. The quark condensate Σ 0 can be estimated for a given strange-quark mass from Table X We can also compare the ratios obtained here with the ones of the N f = 2 simulation of [59]. For the pion masses used in that simulation, m π = 225 and 315 MeV, the ratios {m K /m π , m π /f π , m K /f K } are {2.25, 2.31, 4.508}  Table III of [59] are compatible with our error bands. This indicates that the setup of the simulation of [59] is in line with the result of this analysis for the m s = m 0 s chiral trajectory. Thus, possible deviations in the ρ-meson parameters with the ones obtained here might be caused by a different reason. Notice also that the method used to determine these ratios in [59] is to take m K /f K to the physical point in a strange-quark quenched approximation. As a consequence, the values of f π obtained are consistent with the ones from Fit IV at m s = m 0 s and its extrapolation to the physical value. However, since the real world has more than two flavors, that approach misses the m s dependence of f π as discussed before.
Phase shift lattice data and solutions from Fit IV at the corresponding pion masses are depicted in Figs. 14, 15, and 16, where one can see that lattice data are very well described also inside the 95 % CI. The only exceptions are few data points from the CLS data for the N200 and N401 ensembles (right-bottom panel in Fig. 15), which lie outside the error band and are also far from the bulk of data. Beyond that, most of data for these ensembles are inside of the 95% CI error bands. This is, the N200 and N401 ensembles are compatible within uncertainties, as concluded also in [79]. Note that the error bands are larger for the TrM = C data since they include the variation of the lattice spacing a β between A and B scale settings.
The extrapolation to the physical point in comparison 14 Where we have divided f K from [59] by a factor √ 2 according to the normalization of f K used here.
Tr M(β = 3.4) 268 (20) Tr M(β = 3.55) 254 +11(7) −18 (18) Tr M(β = 3.7) 257 +12(7) −17 (19) msB0 224 +14(10) −18 (20)  with experimental data is plotted in Fig. 17, where one can see that it indeed provides an excellent description of the experiment. In Fig. 18, the CLS D101 ensemble and HS data for m π = 236 MeV are plotted together. We can see that indeed both results are compatible. The phase-shift solution for the pion mass used in the D101 ensemble (dark green) is shown in comparison with the result from an individual fit of the C101 data (light green) in Fig. 19. 15 Note that even when the C101 ensemble was not included in the global Fit IV, both solutions, and the bulk of data itself, lie well inside the 95% CI. The difference between both fits for the energy at which the phase shift crosses 90 • (E (δ = 90 o )/m π ) is of around 3%. This indicates that the deviations due to the volume size are not large.
To show the trend of the TrM = C data, the mean solution of Fit IV together with the lattice data are represented in Fig. 20. 16 As can be seen, these data are now well described. Phase shift data corresponding to higher pion masses fall more to the right and the ρ-meson mass increases monotonically with the pion mass. The dependence of the ρ-meson mass 17 with the pion mass is depicted in Figs. 21 and 22 for the m s = k and TrM = C trajectories, respectively, where we also show the values of the ρ-meson mass given in the corresponding lattice papers. 18 The m ρ /m π ratios are also represented in the 17 Defined as the value of the energy for which δ = 90 o . 18 For the TrM = C trajectories, the error due to the use of scale right panels. For clearness, we present again the lattice data and the resulting curves separately in Fig. 23. In both trajectories m ρ increases with m π . Furthermore, for the trajectories m s = m 0 s and TrM = TrM 0 , we find almost identical results till pion masses of around 400 MeV, when these start to separate. The reason for this behavior is well understood. On one side, the ρ(770) meson becomes a bound state at pion masses around m π = 450 MeV in the m s = m 0 s trajectory (above this value, the settings A and B is also depicted. On the other side, it starts to decay into KK in the TrM = TrM 0 trajectory when the ρ-meson pole crosses this threshold and the kaon gets lighter than the pion. Indeed, it becomes a pole in the IV Riemann sheet as defined in Eqs. (39) and (40). Conversely, other TrM = C trajectories tend to be flatter than the m s = k ones. This is actually in line with the trend of lattice data. It is relevant to note that close to the physical point we do not observe any relevant change in the ρ-meson mass. This suggests that the ρ-meson properties are quite stable against small variations of the strangeness around the physical point, however, its coupling to the KK channel is still large, around 60% of its coupling to ππ. 19 . Nev- 19 This can be seen in Figs  strange quarks, the effect in the real part of the pole, see Fig. 26, becomes significant.
This behavior of the ρ-meson mass is also visible in the corresponding m ρ /m π plots (right panel of Figs. 21 and 22), where the errors in the y-axis are reduced. These plots also show that the error due to the lattice spacing (or scaling setting used) is smaller than the reduction of the ρ-meson mass around the m s = 0 limit. The behavior of the ρ-meson mass and width respect to the kaon mass is depicted in Fig. 26 for the m u = c and m π = m 0 π trajectories. When m K decreases, both mass and width decrease, as commented before. In Figs. 24 and 25, the continuation of the m ρ /m π ratios for smaller pion masses are also depicted. This difference with respect to the physical point is more abrupt as the quark masses get smaller. The symmetric line is also plotted in dot-dashed lines.
In Figs. 27, 28, 29 and 30 we provide the real and imaginary parts of ρ-meson pole position in a 3D plot respect both, the pion and kaon mass. To render some references, we also give in Table XI Table XI, consistently with previous analyses [42]. Nevertheless, as m π increases, the ρ-meson mass moves slower than the ππ threshold, so that, eventually, the ρ(770) meson becomes a ππ bound state with a mass around 908 MeV for a pion mass of around 450 MeV. Note that, in the case of the m s = m 0 s trajectory, the ρ-meson pole is always below the KK threshold for the pion masses analyzed here. We do not start to appreciate significant changes in the behavior of the ρ-meson mass till the strange-quark mass is reduced in half its physical value. For instance, for m s = 0.6 m 0 s , we still get a pole at 732 − i 81 MeV in the chiral limit, which transforms into a bound state with a mass of 905 MeV also for pion masses of around 450 MeV. For lighter strange-quark mass trajectories relevant changes are observed. Both, ρ-meson mass and width, decrease consistently, so that we obtain E 0 = 678 − i 38 MeV for m s = 0 in the chiral limit. Furthermore, for m s ≤ 0.4 m 0 s , both the ππ and KK thresholds get closer to each other as m π increases, in such a way the kaon becomes lighter at a given point. In this regime, the ρ meson becomes a pole in the fourth Riemann sheet 20 when its mass gets below the ππ threshold. In this case, the ρ-meson decays only into KK and its width starts increasing again until it gets a maximum, after which the ρ eventually becomes a KK bound state as m π increases.
This behavior is even more noticeable for the TrM = C trajectories, depicted in Figs. 29 and 30. In this case, the strange-quark mass decreases as m π grows reaching the symmetric m s = m ud line for pion masses of around 450 MeV. Once the symmetric line is crossed, the KK channels opens below the two-pion threshold and the ρ(770) meson becomes again a pole on the fourth Riemann sheet. Nevertheless, the kaon mass in this trajectory decreases till it ends up in the m s = 0 line (red-solid curve). Hence, the ρ-meson mass (width) starts decreasing (increasing) at a given point (when the kaon gets lighter than the pion after crossing the symmetric line) till ending at the zero strangeness line. This behavior suggests that strangeness plays an important role in the ρ(770) meson near the SU(3) flavor limit. This can also be inferred from the increase of its coupling to KK, as discussed below. 21 The pion mass dependence of the ρ-meson couplings to the ππ and KK channels, g ππ and g KK as defined in Eq. (41), are shown in Figs. 31 and 32 for the m s = k and TrM = C trajectories, respectively. On one hand, g ππ varies smoothly with m π before the the ρ-meson transition into a bound state, decreasing as it approaches the KK threshold and increasing with m s . Once the ρ meson becomes bound, its coupling rises sharply till m π 480 MeV. Overall, it takes values g ππ 5.5 − 6.3 for the range of pion masses studied. On the other hand, and contrary to the g ππ behavior, g KK decreases significantly with the mass of the strange quark. In addition, while the pion-mass dependence of g KK flattens for lighter strange quarks, it becomes larger as m s reaches the physical value. All in all, it takes values within the range g KK 3.2 to 4.6.
More information can be extracted when the ratio of both couplings is depicted, we refer to Fig. 33. Remarkably, for those regions where m π ≤ m K one observes the ratio g KK /g ππ ≤ 1/ √ 2. On the contrary, g KK /g ππ > 1/ √ 2 when m π > m K . In the symmetric line we obtain exactly g KK /g ππ = 1/ √ 2. This is not a coincidence. In the SU(3) limit, the decomposition of a I = 1, I 3 = 0 state of the antisymmetric octet representation into two-Goldstone-Boson states with well defined isospin reads where, Y stands for the hypercharge and I is the isospin. Thus, taking into account the kaon degeneracy due to strangeness, the ρ-meson coupling to pions should be a factor √ 2 times larger than for kaons. Notably, the IAM analysis presented here reproduces exactly the SU(3) limit prediction.
In Fig. 34, the ratio √ 2g ππ f π /m ρ is depicted. 22 This ratio lies within the interval [0.95, 1.1], i.e., close to 1, for the quark masses studied in this work, as predicted by the KSFR relation [82]. Thus, we find that KSFR is qualitatively valid, being more accurate near the chiral limit, specially for the m s = {0, m ud } curves, and well applicable also around the physical point, with deviations from KSFR of less than 4%. The largest deviations (still smaller than 8%) are found for pion masses between 200 and 300 MeV.
To compare the LECs obtained in the different analyses done here with the Flag average [78], we depict them in Fig. 35, where we show the results from Fit I (pseudoscalar meson mass and decay constant ratios), Fit II, (analysis of data over the m s = k trajectories), Fit III, (TrM = C trajectories), and Fit IV (mean and standard deviation as a result of the study including both m s = k and TrM = C trajectories) together with the Flag average (pink color). Indeed, we see that LECs from fits I and IV are very close, being also consistent with the 22 In this figure, mρ means the real part of the pole position, ReE 0 , so that we are able to plot the ratio for a larger range of pion masses and after the transition. Since pole positions are slightly lower than the energy corresponding to δ = 90 0 , this ratio gives values lower in around 1.7%, than if E (δ = 90 • ) is taken.
FLAG average, which has larger errors. In general, fits I, II, IIIA and IV give closer results, while the LECs L 6 , L 7 and L 8 , from the analyses of PACS-CS and JL/TWQCD data, strongly disagree with other analyses. Notice the precise values of the LECs provided by Fit IV. As we did before when discussing f π , we can study the ρ-meson properties in a world where there are no kaons or etas by setting to zero their corresponding interacting terms, i.e., contact mass terms and diagrams involving loops with kaons and etas in Fig. 1 and pions in the initial and final state, (what we called t K,η ). In practice, this means that we solve the one-loop ππ IAM equation, see Eq. (30), taking the limits f K,η → ∞, and m K,η → 0. We refer the reader to a more extended explanation in Sect. VI A. Then, the pion mass and decay constant are given now by M π and F π in Eqs. (53) and (54). The result is shown in Fig. 36. Taking into account that we did not include N f = 2 data in our analysis but this result comes out as a prediction from our SU(3) IAM analysis, the agreement with the N f = 2 data is astonishing.
In fact, starting from a three-flavor formulation, the N f = 2 formalism should be in principle obtained when one decouples the strange-quark contribution by sending m s to infinity [45]. In this case the effect of kaons and etas is encoded into the bare pion decay constant, f 0 , the constant B 0 , and in terms proportional to ν K,η = log(m 2 K,η /µ 2 ) + 1 /32π 2 (m K,η refer to the limit where m ud → 0, which can be absorbed in a redefinition of the LECs. Here, though, we are simply studying the world where there are no kaons or etas relying on the fact that their interaction with pions comes from terms where their masses and decay constants appear explicitly. Lattice N f = 2 simulations typically use either f K [59,153] or the nucleon mass (via the QCD static potential) [54][55][56]58] to fix the scale, which in general requires to perform a chiral extrapolation. Nevertheless, given the observed dependence of m π and f π on the strange quark, and the fact that f π is correlated with f K and the nucleon mass, these methods seems to neglect this dependence. In fact, the agreement between N f = 2 lattice simulations and the IAM prediction without kaons and etas on the ρ-meson properties suggests that N f = 2 lattice simulations leave out the contributions coming from the strange quark, and hence, they describe a world where the strange-quark is missing.
In Fig. 37 we compare this result with Fit IV over the chiral trajectories m s = {m 0 s , m ud , 0} and m u = m 0 u . Moreover, we also show the result of solving the onechannel (ππ) IAM equation with m s = {0, m 0 s }, this is, keeping the t k,η terms in the ππ channel. Remarkably, the one-channel IAM result for t K,η = 0 (orange line), m s = 0 (dashed-red), and the two-coupled channel solution over m s = m ud (dotted-red), provide very close results for the ρ mass, which are also consistent with the N f = 2 lattice data. This is explained because the contribution of t K,η terms for light strange-quark masses is small. As explained before, these terms contribute in around 1 − 1.5 MeV of the f π value when m s = 0, and It is also interesting to see what happens if one keeps the t K,η terms in f π and in the ππ scattering amplitude when solving the one-channel ππ IAM equation in the m s = m 0 s trajectory (dot-dashed blue). We see that this trajectory is consistent with the coupled-channel IAM solution for m s = m 0 s , telling that the effect of the offdiagonal elements t 12 in Eq. (20) is very small for physical strange-quark masses. Nevertheless, the coupled-channel effect becomes appreciable for lighter m s , as as one can see by comparing the difference between the dashed and continuous red lines. The small contribution of the offdiagonal elements at the physical point is in contradiction with the results in [31,59], where the absence of these elements are found to be responsible for the dropping of the ρ mass in the N f = 2 case. However, this is natural since in these works the same value of f π was used in the N f = 2 and N f = 2 + 1 predictions, and then, the effect of the kaon and eta contributions were absorbed in the off-diagonal elements. Indeed, we obtain here similar predictions for the ρ meson mass in N f = 2 and N f = 2 + 1 simulations over m s = m 0 s than in [31,59].
Nonetheless, we have gone through a deeper analysis in this work; by studying the m s dependence of f π , we have obtained that m s regulates the contribution of the kaon and eta interacting terms and that the effect of these loops are absorbed in f π instead. Consistently with the predictions done in the works of [31,59], we also obtain that the ρ mass is reduced when these terms are omitted. In Ref. [59], the pion decay constant was determined and these values were used to make predictions for the ρ mass using the UChPT model in [46] with two and three flavors. Nevertheless, the pseudoscalar meson decay constants in Ref. [59] were determined following the method of Ref. [153] where the kaon is introduced later in the quenched approximation and m K /f K is fixed to the physical point. This leads to extrapolated values of f π in N f = 2 simulations consistent with the experimental value, where more than two flavors do exist. However, by doing this, one is missing the m s dependence and the effect of the kaon and eta loops in the pion decay constant, that we have shown here. These missing effects can lead to discrepancies between observables in N f = 2, such as the values of the ρ mass and pion decay constant determinations. We have shown that a lower value of the ρ mass than the physical one should be reflected also in lower values of f π . In summary, assuming that the pion decay constant is the same in the two and three flavor simulations totally misses the effect of the strange quark and loops containing kaon and eta particles in pion observables.

V. CONCLUSIONS
For the first time, we have studied simultaneously both the light-and strange-quark mass dependence of pseudoscalar meson masses, decay constants and ρ-meson properties, such as its mass, width and couplings to the pion and kaon channels. Our analysis is based on recent lattice data of these observables on the chiral trajectories m s = k and TrM = C. In the analysis we resample pseudoscalar meson observables, energy levels (taking into account covariance matrices), and lattice spacings, providing a satisfactory solution at the 95% confidence level. The IAM proves itself to be able to explain the pseudoscalar meson masses, decay constants and ρ-meson properties over different chiral trajectories. Therefore, the LECs obtained here are the most precise and the only ones up to now that are able to describe the strangeness dependence of these observables. The chiral extrapolation of ρ-meson phase shift data is also in remarkable agreement with experiment.
The dependence of the pion decay constant, f π , with  the strange-quark mass, m s , is also studied for the first time. We have shown that, although to assume that the ratio m π /f π is independent of m s can be a good approximation, the variation of f π with m s is abrupt for light pion masses. Furthermore, this dependence is acting as a regulator of the size of the contribution of loops and contact terms involving kaons and etas. For instance, these terms contribute slightly to f π for m s = 0 but account for around 6 − 7 MeV at m s = m 0 s . This contribution to f π is sufficiently large, so that, their absence is able to explain successfully the lower values of the ρ-meson mass obtained in N f = 2 simulations. Even when we did not analyze here N f = 2 lattice data, but only N f = 2 + 1 simulations, the IAM has demonstrated to be able to describe simultaneously both, the ρ-meson mass over chiral trajectories in two and three flavor lattice simulations. Regarding this last aspect, the results obtained here are consistent with the ones of Refs. [31,59]. However, we obtain here that the ρ-meson mass reduction in two-flavor calculations is due to the absence of the strange-quark mass and the contribution containing strange particles on the pion decay constant. Fixing the pseudoscalar decay constants to the physical point in two-flavor lattice simulations misses this dependence.
Some other interesting effects observed involve the KK channel. First, as m s decreases the ρ-meson mass reduces; when m s approaches zero it drops around 70 MeV respect to its values at the physical point. This is effectively more visible in the m u = m 0 u trajectory. Second, in the m s = m 0 s trajectory, as m π increases and the ρ-meson mass gets closer to the KK threshold, its coupling to KK increases, becoming eventually a bound state. Around m π = 450 MeV, it starts to decay into KK in the TrM = TrM 0 trajectory. For other trajectories, these transitions occur at different pion masses when the kaon becomes lighter than the pion. Third, the coupling ratio g ππ /g KK = √ 2 at the symmetric line, factor which comes from a SU(3) Clebsch-Gordan coefficient. Thus, SU(3) flavor symmetry is recovered in the symmetric line.
Our analysis also shows the operators that could be relevant in the energy region and for the light-and strangequark masses considered in the lattice simulation . We hope that the results obtained here motivate the lattice community to investigate more on the hadron properties over different chiral trajectories, which indeed provide useful information to understand their dynamical nature. In this section we analyze kaon and eta contributions into the ρ-meson properties, as well as we discuss how it is possible to disconnect their effect. In the IAM coupledchannel formalism, kaons and etas contribute to pionpion scattering through terms of the kind: Regarding (1), the amplitudes t 12 and t 22 are proportional to 1/f K and 1/f 2 K , respectively, since they involve diagrams with two and four external kaon legs. 23 It is clear that by sending f K → ∞ these contributions disappear. Note, though, that f K is related to f 0 through Eq. (10). Thus, in practice, taking this limit entails breaking the SU(3) symmetry in PCAC [154][155][156], i.e., to assume that the pion and kaon bare decay constants f 0,K and f 0,π do differ. For the same reason, the terms in (2), i.e., kaon and eta tadpoles and one-loop diagrams, are all proportional to 1/f 2 K,η and they also vanish when f K,η → ∞.
Finally, concerning (3), while kaon and eta tadpoles entering in the pion mass and decay constant, Eqs. (5) and (9), vanish when f K,η → ∞, there are still kaon mass contact terms, which can be removed only when one takes the limit m K → 0. Thus, taking the limits m 2 K = 0, f K , f η → ∞, one obtains the pion mass and decay constant in a world where kaons are etas are both decoupled. Namely, and F π =f 0π 1 − 2µ π + 4M 2 0 π f 2 0π (L r 4 + L r 5 ) .
The above relations are the same given in Eqs. (53) and (54), which are depicted in Figs. 12 and 13, and discussed in the paragraphs around these figures.
Then, once all the contributions in (1)-(3), which are called t K,η in the text, are removed, one can solve the one-channel ππ IAM, Eq. (30), t(s) IAM = t 2 (s) 2 t 2 (s) − t 4 (s) , using as input the pion mass and decay constant given in Eqs. (58), (59), to obtain pion-pion scattering in a world where kaons and etas are absent. This is plotted in the orange line in Fig. (37) referred to as t Kη = 0, and discussed in the Sect. II C.

B. Chiral trajectories and decay constant ratios obtained in fits II and III
In Fig. 38, the chiral trajectories and decay constant ratios obtained in fits II and III are depicted.

C. Breit-Wigner reanalyses of IAM solutions
In Fig. 39 we show the refit of the IAM solution for the m s = m 0 s lattice data analyzed in section IV A. The Breit-Wigner parameterization used in these fits for the phase shift is D. Global fits of decay constant and ρ-meson phase-shift data of section IV B In Figs. 40 and 41 we provide the solutions of the TrM = C global fit to ρ-meson phase shift and decay constant lattice data performed in Sect. IV B using separately the scaling methods A or B. As it can be seen, the ensemble N200 (or N401) is not well reproduced with method A.

E. Covariance matrix
In Eq. 61 we provide the correlation matrix of our fitting parameters.