Large $N_c$ scaling of meson masses and decay constants

We perform an $\textit{ab initio}$ calculation of the $N_c$ scaling of the low-energy couplings of the chiral Lagrangian of low-energy strong interactions, extracted from the mass dependence of meson masses and decay constants. We compute these observables on the lattice with four degenerate fermions, $N_f=4$, and varying number of colours, $N_c=3-6$, at a lattice spacing of $a\simeq 0.075$ fm. We find good agreement with the expected $N_c$ scaling and measure the coefficients of the leading and subleading terms in the large $N_c$ expansion. From the subleading $N_c$ corrections, we can also infer the $N_f$ dependence, that we use to extract the value of the low-energy couplings for different values of $N_f$. We find agreement with previous determinations at $N_c=3$ and $N_f=2, 3$ and also, our results support a strong paramagnetic suppression of the chiral condensate in moving from $N_f=2$ to $N_f=3$.

We perform an ab initio calculation of the Nc scaling of the low-energy couplings of the chiral Lagrangian of low-energy strong interactions, extracted from the mass dependence of meson masses and decay constants. We compute these observables on the lattice with four degenerate fermions, N f = 4, and varying number of colours, Nc = 3 − 6, at a lattice spacing of a 0.075 fm. We find good agreement with the expected Nc scaling and measure the coefficients of the leading and subleading terms in the large Nc expansion. From the subleading Nc corrections, we can also infer the N f dependence, that we use to extract the value of the low-energy couplings for different values of N f . We find agreement with previous determinations at Nc = 3 and N f = 2, 3 and also, our results support a strong paramagnetic suppression of the chiral condensate in moving from N f = 2 to N f = 3. The 't Hooft limit of QCD [1] is well known to capture correctly most of its non-perturbative features, such as confinement and chiral symmetry breaking. Large N c inspired approximations are often employed in phenomenological approaches to hadron physics [2][3][4][5][6][7][8][9][10][11], but systematic errors from subleading N c corrections are only naively estimated.
Lattice Field Theory offers the possibility of ab initio explorations of the large N c limit of QCD, by simulating at different values of N c [12,13]. Several studies have already been performed. In Ref. [13] a thorough study of mesonic two-point functions was carried out in the quenched approximation, a limit that captures correctly the leading order terms in N c , but modifies subleading corrections in an uncontrolled way. Furthermore, in Ref. [14] a similar study was performed for N c = 2 − 5 using N f = 2 dynamical fermions at rather high pion masses.
In addition to the standard approach, the study of QCD in the large N c limit can also be achieved using reduced models (see [15] for a review). In this context, there has been significant progress regarding the properties of mesons [16][17][18][19][20].
Besides, lattice simulations have been used to perform studies of various observables in theories with different number of colours, flavours or fermion representations in the context of Beyond-the-Standard-Model theories. Some recent results can be found in [21][22][23][24][25][26] and for recent reviews see [27,28].
In this work, we use previously generated lattice configurations with N c = 3 − 6 and four dynamical fermions. Our particular choice of N f has also advantages for weak matrix elements [29]. On these ensembles, we compute meson masses and decay constants as a function of the quark mass at the different values of N c . We fit these to chiral perturbation theory (ChPT) in order to extract the leading order and next-to-leading order low-energy chiral couplings (LECs). We then study their N c scaling and extract the first two terms in the 't Hooft series.
Interestingly, within the large N c expansion, the 1/N c corrections have a well-defined linear dependence on N f , while the 't Hooft limit is independent on N f . Using this fact, we can predict the low-energy couplings at different values of N f up to higher orders in N c . This allows us to compare with previous determinations, and check the prediction of paramagnetic suppression at large N f of Refs. [46,47]. This paper is organized as follows. First, we describe chiral perturbation theory predictions and the relation to the large N c limit in Section I. In Section II, we present the lattice setup that involves a mixed-action formulation. Next, we explain our scale setting procedure at different N c consistent with 't Hooft scaling in Section III. In Section IV we present the results of our chiral fits to the meson mass and decay constant, first at fixed N c and then combined with the large N c expansion. We also present predictions for theories with different values of N f , compare with previous literature and discuss systematic uncertainties. We conclude in Section V.

I. CHIRAL PERTURBATION THEORY PREDICTIONS
The light spectrum of QCD is the result of the pattern of spontaneous chiral symmetry breaking, ChPT represents accurately the dynamics of the expected pseudo-Nambu-Goldstone bosons (pNGB), i.e., the lightest non-singlet multiplet of pseudoscalar mesons (the octet for N f = 3), at sufficiently small quark masses. The increase in the number of colours while keeping the 't Hooft coupling constant, λ = g 2 N c , is not expected to modify these features. On the other hand, in the large N c limit, QCD reduces to arXiv:1907.11511v1 [hep-lat] 26 Jul 2019 a theory of narrow and non-interacting resonances and, as a result, the interactions of pNGB within the effective theory decrease with N c , improving the convergence of the perturbative series. One complication of the large N c expansion is the role of the singlet pseudoscalar meson, i.e., the η . Its mass originates in the explicit U (1) A breaking by the anomaly. In QCD this contribution to the mass is at the cutoff scale of the chiral effective theory and it is therefore integrated out. However, the anomalous contribution to the singlet mass decreases with N c and in the large N c limit the η becomes degenerate with the remaining pNGBs. The effective theory should consequently include an additional singlet pseudoscalar meson in the spectrum. The corresponding effective theory has been studied long ago [48][49][50][51][52][53]. A new power-counting is needed which involves a simultaneous expansion in 1/N c and the usual chiral expansion in the quark mass and momenta. A consistent power counting was implemented in refs. [52,53]: In the following we will concentrate on the non-singlet multiplet masses and decay constants. We now compare the usual SU (N f ) ChPT to the U (N f ) ChPT for these observables.
At Leading Order (LO) in the standard SU (N f ) chiral expansion there are only two couplings for any number of degenerate flavours, related to the chiral condensate and the meson decay constant. At Next-to-Leading Order (NLO), and for an arbitrary number of degenerate flavours (N f > 3), 13 more LECs are needed, but only two combinations enter in the observables of interest. For N f degenerate flavours, ChPT predicts at NLO [54][55][56]: and in terms of the LO couplings, B, F , and the NLO Gasser-Leutwyler coefficients, L r 4,5,6,8 (µ), defined at the renormalization scale µ.
Eqs. (2)(3) are valid for an arbitrary number of colours, but the LECs scale with N c as (for a review see [57] Loop corrections are suppressed in 1/F 2 π = O(1/N c ), and hence the loop expansion is expected to converge better at larger N c .
Keeping only leading and subleading dependence on N c a convenient parametrization is and Note that according to the scaling of Eq. (4) and the definition of Eq. 7: The NNLO Lagrangian of the SU (N f ) theory is also known [54][55][56]. At this order we will instead only use the U (N f ), to which we now turn. and The NLO corrections are not enough to explain the data in this case, therefore going to NNLO is essential. At NNLO new features appear, because the singlet contributes to the mass loop corrections. The necessary results can be found in Ref. [58]. For degenerate flavours, they simplify to: where M and new LECs that appear in the U (N f ) case. For details see [58]. Note that for degenerate quarks, there is no η-η mixing.
The η mass in this expression can be taken in the large N c limit, where it is given by the Witten-Veneziano formula: where χ t is the topological susceptibility in pure Yang-Mills, recently computed in the large N c limit in Ref. [59].
Note that even though we use the same notation for the LECs in both chiral expansions, they are different: in the SU (N f ) ChPT the LECs encode the effects of integrating out the η . The matching of the two theories starts at NNLO [53,60] and only affects the coupling B and L (1) M of the above [53,60]: with λ 0 = log C. N f versus Nc dependence A diagrammatic analysis of fermion bilinear two point functions shows that within the large N c expansion, the leading order N c → ∞ limit is N f independent and the NLO is O(N f /N c ). We should confirm this expectation also in ChPT formulae above, in particular given the explicit dependence on N f . It turns out that within the U (N f ) expansion, the large N c expansion yields the expected behaviour: the terms in 1/N f exactly cancel when the large N c expansion is taken at fixed M π . We expect therefore that the LECs should also satisfy this same scaling.
On the other hand within the SU (N f ) expansion or in the U (N f ) when M π M η , that is when the chiral limit is taken first, anomalous 1/N f terms appear coming from an expansion in M π /M η . In the U (N f ) expansion such dependence is explicit, but in the SU (N f ) it permeates to the LECs which can no longer be assumed to have the expected O(N f /N c ) dependence, as can be explicitly seen in the matching of L (1) M in Eq. (14). This way, at the order we are working, we can assume the expected scaling in N f of the U (N f ) and SU (N f ) couplings except in the case of [L

II. LATTICE SETUP
We have generated ensembles for SU (N c ) gauge theory with N f = 4 degenerate dynamical fermions, varying N c = 3-6, using the HiRep code [61]. Some of them have been already presented in Ref. [62]. We have chosen the Iwasaki gauge action (following previous experience with 2+1+1 simulations [63,64]) and O(a)-improved Wilson fermions for the sea quarks. In order to achieve automatic O(a) improvement and avoid the need of a non-perturbative determination of normalization factors, we employ maximally twisted valence quarks, i.e., the mixed-action setup [65] previously used in Refs. [66][67][68]. Maximal twist is ensured by tuning the untwisted bare valence mass, m v to the critical value for which the valence PCAC mass is zero: where The bare twisted mass parameter µ 0 is tuned such that the pion mass in the sea and valence sectors coincide, M v π = M s π . The normalized meson decay constant F π can then be obtained from the bare combination [69]: The results for the meson masses and decay constant in the mixed-action setup can be seen in Table II. We have achieved a good tuning of m pcac and the pseudoscalar masses are compatible within one or two sigma with their pure Wilson value (see Table I). When the tuning to maximal twist is not perfect, we correct the bare quark mass (and thus F π ) as follows (see also [69]): where the axial normalization constant, Z A , can be obtained non-perturbatively by matching the valence bare twisted mass with the PCAC mass measured in the sea sector:

III. SCALE SETTING AT LARGE Nc
The scale setting for different values of N c is performed using the gradient flow scale √ 8t 0 , via the determination  of t 0 /a 2 . In QCD, with N c = 3, the standard definition of t 0 is: The leading dependence in N c is known [70] in perturbation theory: where λ GF (q) is the gradient flow 't Hooft coupling at the scale q = 1/ √ 8t. Hence, as in Ref. [59], we will generalize t 0 to an arbitrary N c as: Notice that the choice here is not unique. In particular, one could choose another coupling in a different scheme (such as MS ), and this would induce corrections at order O(N f /N c ) in dimensionful quantities. We also need the value of t 0 in physical units. This is known from lattice simulations for N f = 2 [71,72] and N f = 3 [73] degenerate quarks and at a reference pion mass M ref = 420 MeV: We can use these to perform a linear extrapolation to N f = 4, motivated by the weak N f dependence: Our scale setting condition involves therefore the dimensionless quantity In order to reduce discretization errors we have performed a tree level improvement of t 0 . In Ref. [74], lattice perturbation theory is used to improve t 2 E(t) and thus, t 0 . The prescription is: where the coefficients C 2n depend on the gauge action, the flow action and the definition of E(t) (clover or plaquette). The coefficients for the Iwasaki gauge action, the plaquette action for the flow and the clover definition of E(t) are: The numerical results after the improvement, t imp 0 /a 2 , are shown in Table I.  and the lattice spacing as a function of Nc. The first error is statistical, the second comes from the uncertainty in t0 in physical units, the third stems from the difference in the definitions of E(t) after improvement, and the fourth are finite volume effects estimated from Ref. [75].
Finally, Eq. (25) requires t 0 at M ref . The mass dependence of t 0 has been studied in chiral perturbation theory in Ref. [76]. For degenerate flavours it is given by where k ∝ 1/(F π ) 2 = O(1/N c ) and so the chiral dependence is suppressed in N c . We have performed accordingly a linear fit in M 2 to extract the reference value. The mass dependence of t imp 0 for the different values of N c can be seen in Fig. 1. As expected, the slope is suppressed with N c . The results of the scale setting can be seen in Table III, where we also include the systematic uncertainties. The leading uncertainty comes from the error on the value of t 0 in physical units, the discretization error is estimated from the difference in two definitions of E(t) after improvement, and the finite volume systematic error is estimated from Ref. [75]. As it can be seen, the scale setting yields a uniform lattice spacing for all the values of N c . From now on, we will quote our results in terms of the lattice spacing a = 0.0754 fm, corresponding to N c = 5.

IV. CHIRAL PERTURBATION THEORY FITS
The results for M π and F π in the mixed-action setup are presented in Table II. We want to compare these results to the expectations in ChPT described in Sec. I in order to the extract the LECs and study their N c scaling.
Before addressing the fits, we need to explain the some technical issues regarding the finite volume effects, the renormalization scale and the fitting strategy. We then perform fits at a fixed value of N c to test the ansätze for the N c scaling of the LECs in Eqs. 5 and 7. After that, we perform simultaneous chiral and N c fits. We present a selection of relevant results for the latter, and conclude the section with a discussion on systematic errors.

A. Finite volume effects
Our ensembles have M π L > 3.8 in all cases so we expect finite volume effects to be small and suppressed as 1/N c . Still, we find that for the decay constant they can be of O(1%) and thus we correct them as [77,78]: with ξ ≡ M 2 π (4πFπ) 2 , whileḡ 1 (x) is given bȳ We will use the corrected results for the analysis.

B. Renormalization scale
The NLO couplings are usually defined at µ = 4πF or at the ρ mass, µ = M ρ . Still, in the context of the large N c expansion these are two very different choices, since the former scales with √ N c , deviating from the physical cutoff of the chiral effective theory, which is expected to be set by the lighter resonances, such as the ρ. The scale µ = 4πF is instead the scale at which ChPT breaks down, which for large enough N c is much higher than the scale at which new resonances appear. In the context of large N c , it is therefore sensible to choose a renormalization scale more closely related to the physical cutoff that does not scale with N c . Keeping the scale related to 4πF , however, has some advantages for fitting, so we choose: which has no leading dependence on N c . Using this scale, the NLO predictions can be conveniently written as: where m = µ 0 , the bare twisted mass. Note that in this expression B is bare, since the quark mass is also bare. The value of the non-singlet pseudoscalar normalization constant, Z P , is thus needed.

C. Fitting strategy
Some care is needed to perform the fits in Eqs. (33) and (34). The complication comes from the fact that both coordinates, (x, y) = (ξ, F π ) or (x, y) = (ξ, M 2 π /µ 0 ) have correlated errors. In particular the Ordinary Least Square (OLS) method is not appropriate, since it assumes no errors in x coordinate. An alternative approach is the York Regression (YR) [79], in which the χ 2 function is: where we have defined the two-dimensional vectors: where f is the fitting function, and V is the x, ycovariance matrix, estimated using bootstrap samples. Moreover, we estimate all the errors of the fit parameters via bootstrap resampling.

D. Fit results at fixed Nc
First we consider each N c separately and perform a fit of the data points to extract F, L F (µ) and B, L M (µ). The NLO fit results for these quantities are shown respectively in Tables IV and Tables    We now consider a global fit including several data points at different values of N c . We first perform a SU (N f )-NLO fit to the subset N c = 4 − 6, including leading and subleading N c corrections for all the LO and NLO LECs, as parametrized in Eqs. (5) and (7). We lin- earize the fit by considering the following parametrization where M . Secondly, we consider the U (N f )-NNLO expansion, since we have checked that the U (N f )-NLO expressions fit data very poorly. We also linearize the fit by consid-ering the following fitting functions: and M 2 0 is given by the Witten-Veneziano formula for the η mass valid in the large N c limit (see Eq. 13). We use the result for the topological susceptibility from Ref. [59], We convert to lattice units using the value of t 0 /a 2 in the previous section and substitute F → √ N c F 0 , as extracted from the global F π fit. We find a 0 ∼ 6.5, a value we fix in the fit.
In summary we compare the following fits: i) Fit 1: SU (N f )-NLO fit to Eq. (37) and (38) including the data subset N c = 4-6.
The results for the fitted parameters in the global fits are shown in Tables VI and VII, and the quality of the fits is shown in Figs. 4(a) and 4(b). We also quote in Table  VIII the results for the NLO LECs from these fits. Errors are large, but there are significant correlations between the parameters as can be seen in Fig. 5.

F. Selected results
We will now quote some results that can be inferred from our fits. We first focus on the decay constant in the chiral limit. Using a = 0.0754(23) fm, we get from our fits at fixed N f = 4:

Fit2
: where the N f dependence assumed is the expected one as discussed in sec. I. Note that no N f dependence is assumed in the 1/N 2 c terms. The first error is just the one obtained from the fits in Table VI and the second error of 3% is the one corresponding to the lattice spacing determination. For two-and three-flavour QCD we get: Fit2 :F Nc=3,N f =2 = 81 (7) MeV, where we have taken into account the correlations between the different terms in Eq. (43), and we have assumed no N f dependence on the last term of the Fit 2.
These results are in perfect agreement with phenomenological determinations: [80], and also lattice results (see Ref. [45]). In addition, we can compare to previous results in the large N c limit in the quenched approximation: This value is 2σ away from the results in Eq. (43). This discrepancy may be explained however with the lack of non-perturbative normalization constant and discretization effects, which in their case are of O(a).
Regarding the coupling, B ≡ Σ/F 2 , we do not have a non-perturbative value of Z P , up to this factor we get: Fit2 : Σ F 2 = Z P 1.72(7) − 0. 45(37) From Ref. [82], we can obtain the the 1-loop perturbative result for the normalization constant: which at the order we are working is independent of N f . With this, we obtain for N c = 3: where the first error is systematic, the second comes from the scale setting, and we omit any systematic errors regarding the normalization constant. Combining these results with the ones in Eqs. (44) and (45), we obtain: which is compatible within 1σ with the numbers quoted in Ref. [45]. We can also consider the ratio of condensates for N f = 2 and N f = 3, where the Z P factor drops (up to subleading N f dependence): which shows good agreement with the prediction (11) in Ref. [47].
Regarding the NLO LEC for the decay constant, we get from Fit 1: while for the NLO LEC for the mass, we can only give the N c scaling at N f = 4: In the case of Fit 2, we can provide both results together: For N c = 3, N f = 2, it is more common to quote¯ 3 and This way, we obtain: We stress that U (N f )¯ 3 in fit 2 is not the same as the standard¯ 3 in SU (N f ).¯ 4 agrees instead at 1 − 2σ with the results quoted in Ref. [45].

G. Comments on systematics
The most important systematic uncertainty comes from the finite lattice spacing. Even though a continuum extrapolation would be needed to quantify this error properly, we can get an estimate by comparing the pion mass made of different combination of sea and valence quarks. In particular, Chiral Perturbation Theory in the mixed-action setup predicts that the chiral logs for F π depend upon the mixed pion mass [83]: where m v is the renormalized quark mass in the valence sector and m s in the sea action. We have measured this mixed pion in one ensemble: obtaining a result which is compatible within errors with both, the sea and valence quark pions. A different estimate comes from the dependence on c sw in the valence sector. We have recomputed the decay constant for c sw = 0 in the 3A10 ensemble, obtaining [F π ] csw=0 = 0.04303 (40), within 2% of the value at the nominal c sw . The effects of a change in c sw are in principle O(a 2 ), which can be estimated at ∼ 2% for this observable. This concerns however only the charged meson sector, since the neutral pion is known to have higher discretization effects with twisted mass. That last point is out of the scope of this work, and it will be addressed in future publications.

V. CONCLUSION AND OUTLOOK
In this work we presented the first lattice determination using dynamical fermions of the N c scaling of the couplings in the chiral Lagrangian that contribute to the meson masses and decay constants (see Eqs. (43), (48) and Table VIII). We have been able to disentangle the leading and subleading terms and we found that the subleading contributions are typically non negligible. In fact, within statistical precision, the large N c limit of L M and L F is compatible with zero.
From our chiral fits and theoretical expectations, we have been able to infer the values of the couplings for theories with different numbers of flavours, N f = 2 and N f = 3 at N c = 3. We find that our predictions agree with those in the literature regarding L F , L M and F (see for example Ref. [45] for a summary of results). For B we need to improve our determination, including a nonperturbatively determined renormalization factor. On the other hand, as long as this factor has a small N f dependence, we can estimate the ratio of B and the chiral condensate for N f = 2 and N f = 3. We find excellent agreement with the prediction of paramagnetic suppressions of Refs. [46,47].
We would like to stress that the results presented in this paper are complementary to similar studies that can be performed in reduced models [16][17][18][19][20] or the quenched approximations at large N c [13], since both of these approaches must yield the leading order result as N c → ∞. Given the strong correlations presents in our results (see   As for the future, we would like to mention that our ensembles have a big potential to study other physical observables. We plan to use them to analyse the scaling of other quantities, such as the K → π matrix elements (see [29,62] for previous results). We also believe that the study of scattering amplitudes is a relevant quantity of study at large N c : on one hand quantities such as the I = 2 ππ scattering length give access to LECs of the chiral Lagrangian; on the other hand the study of the behaviour resonances at large N c is interesting, as it may shed light about their nature [10,11,84,85] .