Large $N_c$ scaling of meson masses and decay constants
Pilar Hernández Carlos Pena Fernando Romero-López
hep-lat hep-ph















IFIC (CSIC-UVEG), Edificio Institutos Investigación, 
Apt. 22085, E-46071 Valencia, Spain

Departamento de Física Teórica and Instituto de Física Teórica UAM-CSIC,
Universidad Autónoma de Madrid, E-28049 Madrid, Spain

IFIC (CSIC-UVEG), Edificio Institutos Investigación, 
Apt. 22085, E-46071 Valencia, Spain






We perform an  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≃ 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. 

11.15.Pg,12.38.Gc,12.39.Fe
Large N_c scaling of meson masses and decay constants
    F. Romero-López
    December 30, 2023
=====================================================






The 't Hooft limit of QCD <cit.> 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 <cit.>, 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 <cit.>. Several studies have already been performed. In Ref. <cit.> 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. <cit.> 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 <cit.> for a review). In this context, there has been significant progress regarding the properties of mesons <cit.>. 

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  <cit.> and for recent reviews see <cit.>. 

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 <cit.>. 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. Our study  builds on previous lattice determinations of the LECs for N_c=3 <cit.>, whose main results are summarized in <cit.>.
 
 
 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. <cit.>.

This paper is organized as follows. First, we describe chiral perturbation theory predictions and the relation to the large N_c limit in Section <ref>. In Section <ref>,  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 <ref>. In Section <ref> 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 results for theories with different values of N_f, compare with previous literature and discuss systematic uncertainties. We conclude in Section <ref>.





§ CHIRAL PERTURBATION THEORY PREDICTIONS



The light spectrum of QCD is the result of the pattern of spontaneous chiral symmetry breaking, SU(N_f)_L × SU(N_f)_R → SU(N_f)_L+R. 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 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 <cit.>. 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.  <cit.>:

    𝒪(δ) ∼𝒪(p^2) ∼𝒪(m_q)∼𝒪(m_π^2) ∼𝒪(N_c^-1).

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.  



 §.§ SU(N_f) effective theory

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 <cit.>:

    F_π = F[   1   - N_f/2M_π^2/(4 π F_π)^2logM_π^2/μ^2
       +4 M_π^2/F_π^2( L_5^r  +  N_f L_4^r  )  ],

and 

    M^2_π   = 2B m [1  + 1/N_fM_π^2/(4 π F_π)^2logM_π^2/μ^2
        +8M_π^2/F_π^2(  2L_8^r - L_5^r + N_f (2 L_6^r  - L_4^r )   )  ],

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. (<ref>-<ref>) are valid for an arbitrary number of colours, but the LECs scale with N_c as (for a review see <cit.>):

    [ O(N_c): F^2, L_5, L_8;     O(1): B, L_4, L_6. ]

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

    F = √(N_c)( F_0 + F_1/N_c),      B =   B_0 + B_1/N_c,

and

    L_5  + N_f L_4 ≡ L_F = N_c L_F^(0)  + L_F^(1), 
        2 L_8 - L_5 + N_f (2 L_6 - L_4 )≡ L_M = N_c L_M^(0)  + L_M^(1).

Note that according to the scaling of Eq. (<ref>) and the definition of Eq. <ref>:

    L_F^(0) = L_5  N_c + 𝒪(1  N_c), 
         L_M^(0) = 2 L_8 - L_5  N_c + 𝒪(1  N_c).

The NNLO Lagrangian of the SU(N_f) theory is also known <cit.>. At this order we will instead only use the U(N_f), to which we now turn.



 §.§ U(N_f) effective theory

In the U(N_f) ChPT at NLO, i.e., 𝒪(δ^1), the result can be read from Eqs. (<ref>) and (<ref>) and the different N_c scalings of the LECs in Eqs. (<ref>) and (<ref>):

    F_π = √(N_c)(F_0+F_1 N_c)    [1    + 4M_π^2/F_π^2  N_c  L^(0)_F  +𝒪(δ^2)  ],

and 

    M^2_π = 2 (B_0+ B_1 N_c) m [1   +  8M_π^2/F_π^2  N_c L_M^(0) +𝒪(δ^2) ].


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. <cit.>. For degenerate flavours, they simplify to:

    F_π   =√(N_c)( F_0+ F_1 N_c +  F_2  N_c^2) [1    - N_f/2M_π^2/(4 π F_π)^2logM_π^2/μ^2
       +   4 M_π^2/F_π^2( N_c L_F^(0) + L_F^(1)) + N_c^2 K_F^(0)(M_π^2 F_π^2)^2 
       +    𝒪(δ^3)  ] ,

and 

    M^2_π   = 2  m (  B_0+  B_1 N_c + B_2  N_c^2)  [ 
       1   + 1/N_fM_π^2/(4 π F_π)^2logM_π^2/μ^2 - 1/N_fM_η'^2/(4 π F_π)^2logM_η'^2/μ^2
        +   8M_π^2/F_π^2( N_c L_M^(0) + L_M^(1)) + N_c^2 K_M^(0)(M_π^2 F_π^2)^2 +𝒪(δ^3) ]   ,

where  K_F,M^(0) are combinations of  L_F^(0), L_M^(0) and new LECs that appear in the U(N_f) case. For details see <cit.>. 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:

    M_η'^2 = M_π^2 + 2 N_f  F^2χ_t ≡ M_π^2 + M_0^2,

where χ_t is the topological susceptibility in pure Yang-Mills, recently computed in the large N_c limit in Ref. <cit.>. 

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 <cit.> and only affects the coupling B and L_M^(1) of the above <cit.>:

    [B]_SU(N_f)   = [B]_U(N_f)( 1- 1/N_fM_0^2/(4 π F_π)^2λ_0 ),
    [L_M^(1)]_SU(N_f)   = [L_M^(1)]_U(N_f) -1  8 N_f (4 π)^2(λ_0+1),

with λ_0 = logM_0^2/μ^2.











 §.§ N_f versus N_c 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 𝒪(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 𝒪(N_f/N_c) dependence, as can be explicitly seen in the matching of L_M^(1) in Eq. (<ref>).

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_M^(1)]_SU(N_f). 


§ 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 <cit.>. Some of them have been already presented in Ref. <cit.>. We have chosen the Iwasaki gauge action (following previous experience with 2+1+1 simulations <cit.>) and O(a)-improved[For N_c=3, we take the perturbative value  of c_sw = 1 + c^(1)_sw g^2 from Ref. <cit.>, where we use the plaquette-boosted coupling g^2 = 2 N_c/(β P) = O(1/N_c). For other values of N_c, we use the fact that the one loop coefficient is dominated by the tadpole contribution, which is of order N_c (see Eq. 58 in Ref. <cit.>). This way, c_sw is constant up to subleading corrections in N_c, which have an effect of O(a^2/N_c) in physical observables. The full result cannot be easily reconstructed from Ref. <cit.>.   ] Wilson fermions for the sea quarks. Our simulations use the standard Hybrid Montecarlo (HMC) algorithm with Hasenbusch acceleration. We include five layers in each of the fermionic monomials. Interestingly, we observe that the tuning of the integrator at N_c=3 yields similar results at other values of N_c (at similar pion mass) for the acceptance rate, which we keep at 80-90%. The computational cost of each step in Montecarlo time scales as ∼ N_c^2, with the advantage of a more efficient parallelization at large N_c. 
 
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  <cit.>  previously  used in Refs. <cit.>. 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:

    lim_m^ v→ m_ cr m^ v_ pcac≡lim_m^ v→ m_ cr∂_0⟨A^0(x) P(y)|⟩/2 ⟨P(x) P(y)|⟩= 0,

where A^μ(x) ≡Ψ̅(x) γ^μγ^5 Ψ(x) and P(x) = Ψ̅(x) γ^5 Ψ(x). 
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 <cit.>:

    F_π = 2 μ_0  ⟨0 | P | π|_⟩bare/M_π^2.





The results for the meson masses and decay constant in the mixed-action setup can be seen in Table <ref>. 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 <ref>). When the tuning to maximal twist is not perfect, we correct the bare quark mass (and thus F_π) as follows (see also <cit.>):

    aμ_0    →    aμ_0 √(1 + ( Z_A am_ pcac/aμ_0)^2) , 
      
    a F_π   →    a F_π√(1 + ( Z_A am_ pcac/aμ_0)^2) .

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:

    μ_0 = Z_A m^ s_ pcac,  for  M_π^ v = M_π^s.








§ SCALE SETTING AT LARGE N_C 


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:

    ⟨ t^2 E(t)|⟩|_t=t_0 = c =0.3.

The leading dependence in N_c is known <cit.> in perturbation theory:

    ⟨ t^2 E(t)|=⟩3/128 π^2N_c^2-1/N_cλ_GF(q),

where λ_GF(q)  is the gradient flow 't Hooft coupling at the scale q=1/√(8t). 








Hence, as in Ref.  <cit.>, we will generalize t_0 to an arbitrary N_c as:

    ⟨ t^2 E(t)|⟩|_t=t_0 = c(N_c)=3 8N_c^2-1/N_c c(3).

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  <cit.> and N_f=3 <cit.> degenerate quarks and at a reference pion mass M_ ref=420 MeV:

    √(t_0)|^N_f=2_M_ ref = 0.1470(14)  fm,
    √(t_0)|^N_f=3_ M_ ref = 0.1460(19)  fm

We can use these to perform a linear extrapolation to N_f=4, motivated by the weak N_f dependence:

    √(t_0)|^N_f=4_ M_ ref = 0.1450(39)  fm.

Our scale setting condition involves therefore the dimensionless quantity

    (M_π√(t_0)) |_M_ ref = 0.3091(83).



In order to reduce discretization errors we have performed a tree level improvement of t_0.  In Ref. <cit.>, lattice perturbation theory is used to improve ⟨t^2 E(t)|$⟩ and thus,t_0. The prescription is: 

    ⟨t^2 E(t)|_⟩a =  ⟨t^2 E(t)|_⟩imp[    1+ ∑_n C_2n(a^2/t)^n   ],

where the coefficientsC_2ndepend on the gauge action, the flow action and the definition ofE(t)(clover or plaquette). 

The coefficients for the Iwasaki gauge action, the plaquette action for the flow and the clover definition ofE(t)are:

    C_2 = -0.262333,  C_4 = 0.0936935, 
      C_6 = -0.048002,  C_8=0.0320211.

The numerical results after the improvement,t_0^imp/a^2, are shown in Table <ref>. 



Finally, Eq. (<ref>) requirest_0atM_ref. The mass dependence oft_0has been studied in chiral perturbation theory in Ref. <cit.>. For degenerate flavours it is given by

    t_0  =   t_0^χ(1 + k  M^2 )  + 𝒪(M^4),
 
wherek ∝1/(F_π)^2 = O(1/N_c)and so the chiral dependence is suppressed inN_c. We have performed  accordingly a linear fit inM^2to extract the reference value. The mass dependence oft^imp_0for the different values ofN_ccan be seen in Fig. <ref>. As expected, the slope is suppressed withN_c. The results of the scale setting can be seen in Table <ref>, where we also include the systematic uncertainties. The leading uncertainty comes from the error on the value oft_0in physical units,  the discretization error is estimated from the difference in two definitions ofE(t)after improvement, and the  finite volume systematic error is estimated from Ref. <cit.>. As it can be seen, the scale setting yields a uniform lattice spacing for all the values ofN_c. From now on, we will quote our results in terms of the lattice spacinga = 0.0754  fm, corresponding toN_c=5.














§ CHIRAL PERTURBATION THEORY FITS


The results forM_πandF_πin the mixed-action setup are presented in Table <ref>.  We want to compare these results to the expectations in ChPT described in Sec. <ref> in order to the extract the LECs and study theirN_cscaling. 

Before addressing the fits, we need to explain some technical issues regarding the finite volume effects, the renormalization scale and the fitting strategy. We then perform fits at a fixed value ofN_cto test the ansätze for theN_cscaling of the LECs in Eqs. <ref> and <ref>. After that, we perform simultaneous chiral andN_cfits. We present a selection of relevant results for the latter, and conclude the section with a discussion on systematic errors.

 


 §.§ Finite volume effects


Our ensembles haveM_πL > 3.8in all cases so we expect finite volume effects to be small and suppressed as1/N_c. Still, we find that for the decay constant they can be ofO(1%)and thus we correct them as <cit.>:

    M_π(L)    = M_π[ 1+ 1/2 N_fξ g̅_1(M_π L) + O(ξ^2)], 
    
    F_π(L)    = F_π[ 1- N_f/2 ξ g̅_1(M_π L) + O(ξ^2)],

withξ≡M_π^2/(4 πF_π)^2, whileg̅_1(x)is given by

    g̅_1(x) 24/x K_1(x) ∼24 √(2)/√(π)e^-x/x^3/2.

We will use the corrected results for the analysis.


 §.§ Renormalization scale

The NLO couplings are usually defined atμ=4 πFor at theρmass,μ=M_ρ. Still, in the context of the largeN_cexpansion 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 πFis instead the scale at which ChPT breaks down, which for large enoughN_cis much higher than the scale at which new resonances appear. In the context of largeN_c, it is therefore sensible to choose a renormalization scale more closely related
to the physical cutoff that does not scale withN_c.  Keeping the scale related to4πF, however, has some advantages for fitting, so we choose:

    μ^2 = 3  N_c (4 π F)^2,
 
which has no leading dependence onN_c. Using this scale, the NLO predictions can be conveniently written as:

    F_π   =    F [1   - 2 ξlog(N_c  3ξ)    + 64 π^2 ξ L_F(μ)  ], 
    M_π^2/m   =    2 B   [1   + 1/4ξlog(N_c  3ξ)    + 128 π^2  ξ L_M(μ)  ],

wherem=μ_0, the bare twisted mass. Note that in this expressionBis bare, since the quark mass is also bare. The value of the non-singlet pseudoscalar normalization constant,Z_P, is thus needed.




 §.§ Fitting strategy


Some care is needed to perform the fits in Eqs. (<ref>) and (<ref>). 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 inxcoordinate. An alternative approach is the York Regression (YR) <cit.>, in which theχ^2function is:

    χ^2 = ∑_i  min_δ_i[ 𝐑^T_i  V^-1𝐑_i ],

where we have defined the two-dimensional vectors:

    𝐑_i(δ_i) ≡( f( x_i + δ_i ) - y_i, δ_i  ),

wherefis the fitting function, andVis thex,y-covariance matrix, estimated using bootstrap samples.  In order to account for autocorrelations, we vary the block-size of the bootstrap samples. We find that blocks of ∼ 20 units of Montecarlo are sufficient, and we do not observe a clear N_c dependence. We also estimate all the errors of the fit parameters via bootstrap resampling.

 §.§ Fit results at fixed N_c


First we consider eachN_cseparately and perform a fit of the data points to extractF, L_F(μ)andB, L_M(μ).  The NLO fit results for these quantities are shown respectively in Tables <ref> and Tables <ref>.  TheN_cdependence of the LECs is shown in Figs. <ref> and  <ref>. It can be seen that the scaling is well described by leading and subleadingN_ccorrections forN_c=4-6, while there seems to be significant1/N_c^2corrections forN_c=3in the case ofFandL_F. In the case ofBandL_Merrors are larger and there is no sign of1/N_c^2. Interestingly, the data suggest that the largeN_climit ofL_M ∼0.




 §.§ Simultaneous chiral and N_c fits


We now consider a global fit including several data points at different values ofN_c. We first perform aSU(N_f)-NLO fit to the subsetN_c=4-6, including leading and subleadingN_ccorrections for all the LO and NLO LECs, as parametrized in Eqs. (<ref>)  and (<ref>). We linearize the fit by considering the following parametrization

    F_π   = √(N_c)(F_0+ F_1 N_c) [1   - 2 ξlog(N_c  3ξ)  ]   
       + 64 π^2 ξ√(N_c)(N_c (F L_F)^(0) + (F L_F)^(1)) , 
    M_π^2/m   = 2 (B_0 + B_1  N_c)   [1   + 1/4ξlog(N_c  3ξ)  ]  
       + 256 π^2  ξ(N_c (B L_M)^(0) + (B L_M)^(1)),

where(F L_F)^(0) ≡F_0 L_F^(0), while(F L_F)^(1) ≡F_1 L_F^(0)+ L_F^(1) F_0, and(B L_M)^(0) ≡B_0 L_M^(0)and(B L_M)^(1) ≡B_1 L_M^(0)+ B_0 L_M^(1).

Secondly, we consider theU(N_f)-NNLO  expansion, since we have checked that theU(N_f)-NLO expressions fit the data very poorly. We also linearize the fit by considering 
the following fitting functions:

    F_π   =   √(N_c)(F_0+ F_1 N_c + F_2  N_c^2)  [1   - 2 ξlog(N_c  3ξ) ]
       +    64 π^2 ξ√(N_c)(N_c (FL_F)^(0) + (FL_F)^(1))   
       +     N_c^2 √(N_c)(16 π^2 ξ)^2  K^(0)_F,

    M_π^2/m   =    2 ( B_0 + B_1 N_c + B_2  N_c^2) [  1   +  1/4ξlog(N_c  3ξ)  
       +     -1/4(ξ+a_0 N_c^2) log(N_c  3(ξ+a_0 N_c^2)) ]  
       +   256 π^2 ξ(N_c  (B L_M)^0 + (B L_M)^(1)) 
       -     64  N_c^2 (16 π^2  ξ)^2 K_M^(0),

 where

    a_0 ≡ N_c^2 M_0^2 (4 π F)^2,

andM_0^2is given by the Witten-Veneziano formula for theη'mass valid in the largeN_climit (see Eq. <ref>). 
We use the result for the topological susceptibility from Ref. <cit.>, 

    t_0^2 χ_t = 7.03 (13) · 10^-4.

We convert to lattice units using the value oft_0/a^2in the previous section and substituteF→√(N_c) F_0, as extracted from the globalF_πfit. 
We finda_0 ∼6.5, a value we fix in the fit.


In summary we compare the following fits:

  *  Fit 1: SU(N_f)-NLO fit to Eq. (<ref>) and (<ref>)  including the data subset N_c=4-6. 

  *  Fit 2: U(N_f)-NNLO expansion fit  to Eqs. (<ref>) and (<ref>) including the full data set.









The results for the fitted parameters in the global fits  are shown in Tables <ref> and <ref>, and the quality of the fits is shown in Figs. <ref> and <ref>. We also quote in Table <ref> 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. <ref>. 




 §.§ 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. Usinga= 0.0754(23)fm, we get from our fits at fixedN_f=4:

    Fit 1: F√(N_c)   =  (67(3) - 26(4)N_f N_c)(3%)^a MeV,  
     Fit 2: F √(N_c)   = (70(2)- 22(5)N_f N_c -  86(37) N_c^2)(3%)^a  MeV,

where theN_fdependence assumed is the expected one as discussed in sec. <ref>. Note that noN_fdependence is assumed in the1/N_c^2terms. 
The first error is just the one obtained from the fits in Table <ref> and the second error of3%is the one corresponding to the lattice spacing determination. For two- and three-flavour QCD we get:

    Fit 1:     F^N_c=3, N_f=2 = 86(3)   MeV, 
        F^N_c=3, N_f=3 = 71(3)   MeV, 
     Fit 2:     F^N_c=3, N_f=2 = 81(7)   MeV,  
        F^N_c=3, N_f=3 = 68(7)   MeV,

where we have taken into account the correlations between the different terms in Eq. (<ref>), and we have assumed noN_fdependence on the last term of the Fit 2. These results are in perfect agreement with phenomenological determinations: 

    F^N_f=2   = 86.2(5)  MeV in Ref. <cit.>,   
    
    F^N_f=3   ≃ 71.1   MeV in Ref. <cit.>,

 and also lattice results (see Ref. <cit.>). In addition, we can compare to previous results in the largeN_climit in the quenched approximation:

    F√(N_c)|_N_c →∞ = 56(5)  MeV ,    Ref.  <cit.>.

This value is2σaway from the results in Eq. (<ref>). This discrepancy may be explained however with the lack of non-perturbative normalization constant and discretization effects, which in their case are ofO(a).


Regarding the coupling,B≡Σ/F^2, we do not have a non-perturbative value ofZ_P, up to this factor we get:

    Fit 1: Σ F^2   =    Z_P (  1.70(11)- 0.12(12) N_f N_c),
     Fit 2: Σ F^2   =    Z_P (1.72(7) - 0.45(37)N_f N_c -  1.8(1.5)  N_c^2).











From Ref. <cit.>, we can obtain the the 1-loop perturbative result for the normalization constant: 

    Z_P(N_c=3) = 0.555,

 which at the order we are working is independent ofN_f. With this, we obtain forN_c=3:

    Fit 1:      N_f=4    ⟶    Σ/F^2= 2.26(11) (7)^a  GeV,  
         N_f=3    ⟶    Σ/F^2= 2.31(5) (7)^a  GeV, 
         N_f=2    ⟶    Σ/F^2= 2.35(3) (7)^a  GeV,

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. (<ref>) and (<ref>), we obtain:

    Σ^1/3(N_f=2)    = 257(2)(9)^a  MeV, 
    Σ^1/3(N_f=3)    = 223(4)(8)^a  MeV,

which is compatible within1σwith the numbers quoted in Ref. <cit.>.
We can also consider the ratio of condensates forN_f=2andN_f=3, where theZ_Pfactor drops (up to subleadingN_fdependence):

    Σ(N_f=2)/Σ(N_f=3) = 1.49(10),

which shows good agreement with the prediction  

    Σ(N_f=2)/Σ(N_f=3) = 1.51(11)  in Ref. <cit.>.


Regarding the NLO LEC for the decay constant, we get from Fit 1:

    L_F(μ)/N_c· 10^3 =     0.1(4) +  0.6(3) N_f/N_c + O(N_c^-2),

while for the NLO LEC for the mass, we can only give theN_cscaling atN_f=4:

    L_M^N_f=4(μ) /N_c· 10^3 
      =     -0.2(2) +  2.9(6) /N_c+ O(N_c^-2).

In the case of Fit 2, we can provide both results together:

    L_F(μ)/N_c· 10^3 =  -0.3(4) +  0.4(4) N_f/N_c + O(N_c^-2), 
       L_M(μ)/N_c· 10^3 =  -0.1(1) +  0.9(2) N_f/N_c+ O(N_c^-2).

From Eqs. <ref> and <ref>, we can infer theN_c=3, N_f=3results:

    Fit 1: L_F(μ)    =  2.1(3) · 10^-3,   
    
     Fit 2: L_F(μ)    =  0.4(2.1) · 10^-3,  
     L_M(μ)    =  2.4(8) · 10^-3.



ForN_c=3, N_f=2, it is more common to quoteℓ̅_3andℓ̅_4:

    ℓ̅_3 = 2 log(  4 π F_π M_π^phys) - 16 (4 π)^2 L_M^N_f=2, N_c=3, 
    ℓ̅_4 = 2 log(  4 π F_π M_π^phys) + 4 (4 π)^2 L_F^N_f=2, N_c=3.

This way, we obtain:

    Fit 1:    ℓ̅_4 = 5.1(3), 
    Fit 2:    ℓ̅_3 = 0.4(1.6),    ℓ̅_4 = 4.1(1.1).

We stress thatU(N_f)ℓ̅_3in fit 2 is not  the same as the standardℓ̅_3inSU(N_f).ℓ̅_4agrees instead at1-2σwith the results quoted in Ref. <cit.>.


























































 §.§ 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 forF_πdepend upon the mixed pion mass <cit.>:

    M_π^mixed 2B(m^v + m^s),

wherem^vis the renormalized quark mass in the valence sector andm^sin the sea action. We have measured this mixed pion in one ensemble:

    Ensemble 3A10  → aM_π^mixed = 0.2201(26),

obtaining a result which is compatible within errors with both, the sea and valence quark pions.

A different estimate comes from the dependence onc_swin the valence sector. We have recomputed the decay constant  forc_sw =0in the 3A10 ensemble, obtaining[F_π]_c_sw =0 = 0.04303(40), within2%of the value at the nominalc_sw. 
























The effects of a change inc_sware in principleO(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 issue is out of the scope of this work, and it will be addressed in future publications. 
inN_cWe end this section with a last word on the chiral fits. We find that our data is well described by ChPT at the order we worked. Still, we cannot exclude that higher order corrections might be relevant in the range of masses we are considering. A robust study on the convergence of ChPT would require simulations at lighter quark masses and a proper continuum extrapolation.


§ CONCLUSION AND OUTLOOK


In this work we presented the first lattice determination using dynamical fermions of theN_cscaling of the couplings in the chiral Lagrangian that contribute to the meson masses and decay constants (see Eqs. (<ref>), (<ref>) and Table <ref>). We have been able to disentangle the leading and subleading terms and we found that the subleading contributions are typically non negligible. In fact, we find that the value for L_M at N_c=3 seems to be dominated by the subleading corrections, and the fit result suggests an accidental cancellation of 2L_8 -L_5 in the large N_c limit.



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=2andN_f=3atN_c=3. We find that our results nicely agree with those in the literature regardingL_F, L_MandF(see for example Ref. <cit.> for a summary of results). ForBwe need to improve our determination, including a non-perturbatively determined renormalization factor. On the other hand, as long as this factor has a smallN_fdependence, we can estimate the ratio of B and the chiral condensate forN_f=2andN_f=3. We find excellent agreement with the prediction of paramagnetic suppressions of Refs. <cit.>. 

We would like to stress that the results presented in this paper are complementary to similar studies that can be performed in reduced models <cit.> or the quenched approximations at largeN_c<cit.>, since both of these approaches must yield the leading order result asN_c →∞. Given the strong correlations presents in our results (see Fig. <ref>), a precise determination of the dominantN_cterm would significantly improve the determination of the subleadingN_ccorrections, and hence the determination of the physical values atN_c=3. We are willing to provide the bootstrap samples if requested.

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 theK →πmatrix elements (see <cit.> for previous results). We also believe that the study of scattering amplitudes is a relevant quantity of study at largeN_c: on one hand quantities such as theI=2ππscattering length give access to LECs of the chiral Lagrangian; on the other hand the study of the behaviour resonances at largeN_cis interesting, as it may shed light about their nature <cit.> .











We thank Andrea Donini for very useful discussions and previous collaboration on related work, as well as  M. García Pérez, A. González-Arroyo, G. Herdoíza,  A. Ramos, A. Rusetsky, S. Sharpe, C. Urbach and A. Walker-Loud for useful comments and suggestions. 
We are particularly grateful to Claudio Pica and Martin Hansen for providing us with a SU(N_c) lattice code. This work was partially supported by grant FPA2017-85985-P, MINECO's “Centro de Excelencia Severo Ochoa” Programme under grant SEV-2014-0398, 
and the European projects 
H2020-MSCA-ITN-2015/674896-ELUSIVES and H2020-MSCA-RISE-2015/690575-InvisiblesPlus. The work of FRL has also received funding from the European Union Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 713673 and "La Caixa" Foundation (ID 100010434, LCF/BQ/IN17/11620044). Furthermore, CP thankfully acknowledges support through the Spanish projects FPA2015-68541-P (MINECO/FEDER) and PGC2018-094857-B-I00, the Centro de Excelencia Severo Ochoa Programme SEV-2016-0597, and the EU H2020-MSCA-ITN-2018-813942 (EuroPLEx)

We thank Mare Nostrum 4 (BSC), Finis Terrae II (CESGA), Tirant 3 (UV) and Lluis Vives (Servei d'Informàtica UV) for the computing time provided.