Large $N_c$ scaling of meson masses and decay constants

TL;DR

This study uses lattice QCD simulations to analyze how meson masses and decay constants scale with the number of colors ($N_c$), confirming theoretical expectations and exploring the dependence on the number of flavors ($N_f$).

Contribution

It provides the first ab initio lattice calculation of $N_c$ scaling of low-energy couplings, including subleading corrections and $N_f$ dependence, for $N_c=3$ to 6.

Findings

Good agreement with expected $N_c$ scaling.

Measured coefficients of large $N_c$ expansion.

Supports strong paramagnetic suppression of chiral condensate.

Abstract

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$.

Tables

Table 1. Table 1: Summary of our ensembles: β 𝛽 \beta , sea quark bare mass parameter, m s superscript 𝑚 𝑠 m^{s} , and sea pion mass M π s subscript superscript 𝑀 𝑠 𝜋 M^{s}_{\pi} . We keep c s ​ w = 1.69 subscript 𝑐 𝑠 𝑤 1.69 c_{sw}=1.69 throughout.

Ensemble	$L^3\times T$	$\beta$	$a{m}^s$	$a{M}_{\pi }^s$	$t_0^{\text{imp}}/a^2$
3A10	${20}^3\times 36$	1.778	-0.4040	0.2204(21)	3.263(50)
3A20	${24}^3\times 48$		-0.4060	0.1845(14)	3.491(32)
3A30	${24}^3\times 48$		-0.4070	0.1613(16)	3.740(39)
3A40	${32}^3\times 60$		-0.4080	0.1429(12)	3.855(27)
4A10	${20}^3\times 36$	3.570	-0.3725	0.2035(14)	3.494(45)
4A20	${24}^3\times 48$		-0.3752	0.1805(7)	3.565(26)
4A30	${24}^3\times 48$		-0.3760	0.1714(8)	3.593(29)
4A40	${32}^3\times 60$		-0.3780	0.1397(8)	3.723(23)
5A10	${20}^3\times 36$	5.969	-0.3458	0.2128(9)	3.532(17)
5A20	${24}^3\times 48$		-0.3490	0.1802(6)	3.614(18)
5A30	${24}^3\times 48$		-0.3500	0.1712(6)	3.664(24)
5A40	${32}^3\times 60$		-0.3530	0.1331(7)	3.776(19)
6A10	${20}^3\times 36$	8.974	-0.3260	0.2150(7)	3.619(17)
6A20	${24}^3\times 48$		-0.3300	0.1801(5)	3.696(17)
6A30	${24}^3\times 48$		-0.3311	0.1689(7)	3.721(15)
6A40	${32}^3\times 60$		-0.3340	0.1351(6)	3.820(17)

Table 2. Table 2: Results obtained in the mixed action setup, with Wilson fermions on the sea and twisted mass in the valence sector. We use c s ​ w = 1.69 subscript 𝑐 𝑠 𝑤 1.69 c_{sw}=1.69 , as in the sea sector.

Ensemble	$a{m}_{\mathrm{cr}}$	$a{\mu }_0$	$a{M}_{\pi }^v$	$\left|a{m}_{pcac}^v\right|$	$a{F}_{\pi }$
3A10	-0.4214	0.01107	0.2216(20)	0.0000(3)	0.04405(41)
3A20	-0.4196	0.00781	0.1834(6)	0.0001(2)	0.04023(24)
3A30	-0.4187	0.00632	0.1613(11)	0.0008(2)	0.03678(33)
3A40	-0.4163	0.00513	0.1423(7)	0.0006(3)	0.03554(15)
4A10	-0.3875	0.01030	0.2037(11)	0.0001(2)	0.05131(37)
4A20	-0.3865	0.00844	0.1803(9)	0.0000(4)	0.05037(26)
4A30	-0.3865	0.00778	0.1717(9)	0.0001(4)	0.04913(31)
4A40	-0.3851	0.00546	0.1416(5)	0.0001(2)	0.04608(15)
5A10	-0.3611	0.01225	0.2114(13)	0.0003(4)	0.06125(32)
5A20	-0.3611	0.00906	0.1799(10)	0.0001(4)	0.05767(30)
5A30	-0.3607	0.00824	0.1706(13)	0.0000(4)	0.05647(40)
5A40	-0.3596	0.00509	0.1328(5)	0.0002(2)	0.05278(18)
6A10	-0.3415	0.01298	0.2142(6)	0.0003(2)	0.06813(21)
6A20	-0.3414	0.00956	0.1801(4)	0.0002(2)	0.06435(25)
6A30	-0.3414	0.00803	0.1668(5)	0.0002(2)	0.06278(24)
6A40	-0.3409	0.00542	0.1342(4)	0.0000(1)	0.05929(14)

Table 3. Table 3: Results for the t 0 / a 2 | M ref evaluated-at subscript 𝑡 0 superscript 𝑎 2 subscript 𝑀 ref \left.t_{0}/a^{2}\right|_{M_{\text{ref}}} and the lattice spacing as a function of N c subscript 𝑁 𝑐 N_{c} . The first error is statistical, the second comes from the uncertainty in t 0 subscript 𝑡 0 t_{0} in physical units, the third stems from the difference in the definitions of E ​ ( t ) 𝐸 𝑡 E(t) after improvement, and the fourth are finite volume effects estimated from Ref. Fodor et al. ( 2012 ) .

$N_c$	${t_0/a^2|}_{M_{\text{ref}}}$	$a$ ($\times {10}^{-2}$ fm)
3	$3.71\left(4\right){\left(7\right)}_{t_0}{\left(12\right)}_a{\left(3\right)}_L$	$7.53\left(4\right){\left(19\right)}_{t_0}{\left(12\right)}_a{\left(3\right)}_L$
4	$3.64\left(1\right){\left(3\right)}_{t_0}{\left(12\right)}_a{\left(3\right)}_L$	$7.60\left(1\right){\left(20\right)}_{t_0}{\left(12\right)}_a{\left(3\right)}_L$
5	$3.69\left(2\right){\left(3\right)}_{t_0}{\left(12\right)}_a{\left(3\right)}_L$	$7.54\left(2\right){\left(20\right)}_{t_0}{\left(12\right)}_a{\left(3\right)}_L$
6	$3.76\left(1\right){\left(2\right)}_{t_0}{\left(12\right)}_a{\left(3\right)}_L$	$7.48\left(1\right){\left(20\right)}_{t_0}{\left(12\right)}_a{\left(3\right)}_L$

Table 4. Table 4: NLO Fits for F π subscript 𝐹 𝜋 F_{\pi} for separate values of N c subscript 𝑁 𝑐 N_{c} .

$N_c$	$aF/\sqrt{N_c}$	$L_F/N_c$	${\chi }^2/dof$
3	0.0088(9)	0.0046(14)	0.7/2
4	0.0155(6)	0.0013(3)	3.9/2
5	0.0175(4)	0.0011(2)	2.2/2
6	0.0188(2)	0.0011(1)	0.4/2

Table 5. Table 5: Fits for M π subscript 𝑀 𝜋 M_{\pi} for separate values of N c subscript 𝑁 𝑐 N_{c}

$N_c$	$aB$	$L_M/N_c$	${\chi }^2/dof$
3	1.564(55)	0.00086(10)	10.2/2
4	1.560(37)	0.00064(7)	1.4/2
5	1.648(30)	0.00031(6)	0.1/2
6	1.610(20)	0.00031(4)	9.5/2

Table 6. Table 6: Different fits for the decay constant as described in the text.

Fit	$F_0$	$F_1$	$F_2$	${\left(F{L}_F\right)}^{\left(0\right)}$	${\left(F{L}_F\right)}^{\left(1\right)}$	$K_F^{\left(0\right)}$	${\chi }^2/dof$
1	0.0255(12)	-0.040(6)	-	$4.7\left(9.5\right)\cdot {10}^{-6}$	$4.8\left(5.1\right)\cdot {10}^{-5}$	-	0.79
2	0.0266(9)	-0.034(8)	-0.033(14)	$-8\left(10\right)\cdot {10}^{-6}$	$5.6\left(4.4\right)\cdot {10}^{-5}$	$7.6\left(6.4\right)\cdot {10}^{-7}$	0.9

Table 7. Table 7: Different fits for the meson mass as described in the text.

Fit	$B_0$	$B_1$	$B_2$	${\left(B{L}_M\right)}^{\left(0\right)}$	${\left(B{L}_M\right)}^{\left(1\right)}$	$K_M^{\left(0\right)}$	${\chi }^2/dof$
1	1.70(11)	-0.5(5)	-	-0.00046(29)	0.0056(15)	-	2.0
2	1.72(7)	-1.8(5)	1.8(1.5)	-0.00017(25)	0.0066(10)	$1.3\left(9\right)\cdot {10}^{-6}$	2.4

Table 8. Table 8: Values for the LECs from the fits in Tables 6 and 7 .

Fit	$L_F^{\left(0\right)}$	$L_F^{\left(1\right)}$	$L_M^{\left(0\right)}$	$L_M^{\left(1\right)}$
1	$1\left(4\right)\cdot {10}^{-4}$	$23\left(13\right)\cdot {10}^{-4}$	-$20\left(15\right)\cdot {10}^{-5}$	$29\left(6\right)\cdot {10}^{-4}$
2	-$3\left(4\right)\cdot {10}^{-4}$	$17\left(18\right)\cdot {10}^{-4}$	-$1\left(1\right)\cdot {10}^{-4}$	$37\left(7\right)\cdot {10}^{-4}$

Equations

$$\mathcal{O}(\delta )\sim \mathcal{O}(p^2)\sim \mathcal{O}(m_q)\sim \mathcal{O}(m_{\pi }^2)\sim \mathcal{O}(N_c^{-1}).$$

$$\displaystyle\begin{split}F_{\pi}={F}\Bigg{[}&1-{\frac{N_{f}}{2}\frac{M_{\pi}^{2}}{(4\pi F_{\pi})^{2}}\log\frac{M_{\pi}^{2}}{\mu^{2}}}\\ &+4\frac{M_{\pi}^{2}}{F_{\pi}^{2}}\Big{(}L_{5}^{r}+N_{f}L_{4}^{r}\Big{)}\Bigg{]},\\ \end{split}$$

$$\displaystyle\begin{split}M^{2}_{\pi}&=2Bm\Bigg{[}1+{\frac{1}{N_{f}}\frac{M_{\pi}^{2}}{(4\pi F_{\pi})^{2}}\log\frac{M_{\pi}^{2}}{\mu^{2}}}\\ &+8\frac{M_{\pi}^{2}}{F_{\pi}^{2}}\Big{(}2L_{8}^{r}-L_{5}^{r}+N_{f}(2L_{6}^{r}-L_{4}^{r})\Big{)}\Bigg{]},\end{split}$$

$$\displaystyle\begin{array}[]{ll}O(N_{c}):F^{2},L_{5},L_{8};&O(1):B,L_{4},L_{6}.\end{array}$$

$$F=\sqrt{N_c}\left(F_0+\frac{F_1}{N_c}\right),B=B_0+\frac{B_1}{N_c},$$

$$L_5+N_f{L}_4\equiv L_F=N_c{L}_F^{(0)}+L_F^{(1)},$$

$$2L_8-L_5+N_f(2L_6-L_4)\equiv L_M=N_c{L}_M^{(0)}+L_M^{(1)}.$$

$$\begin{array}{rl}& L_F^{(0)}=\frac{L_5}{N_c}+\mathcal{O}\left(\frac{1}{N_c}\right),\\ & L_M^{(0)}=\frac{2L_8-L_5}{N_c}+\mathcal{O}\left(\frac{1}{N_c}\right).\end{array}$$

$$\displaystyle\begin{split}F_{\pi}=\sqrt{N_{c}}\left(F_{0}+{F_{1}\over N_{c}}\right)&\Bigg{[}1+4\frac{M_{\pi}^{2}}{F_{\pi}^{2}}N_{c}{L}^{(0)}_{F}+{\mathcal{O}}(\delta^{2})\Bigg{]},\\ \end{split}$$

$$\displaystyle\begin{split}M^{2}_{\pi}=2\left(B_{0}+{B_{1}\over N_{c}}\right)m\Bigg{[}1+8\frac{M_{\pi}^{2}}{F_{\pi}^{2}}N_{c}{L}_{M}^{(0)}+{\mathcal{O}}(\delta^{2})\Bigg{]}.\end{split}$$

$$F_{\pi }$$

$$M_{\pi }^2$$

$$M_{{\eta }^{\prime }}^2=M_{\pi }^2+\frac{2N_f}{F^2}{\chi }_t\equiv M_{\pi }^2+M_0^2,$$

$$\begin{array}{rl}{\left[B\right]}_{SU(N_f)}& ={\left[B\right]}_{U(N_f)}\left(1-\frac{1}{N_f}\frac{M_0^2}{(4\pi F_{\pi }{)}^2}{\lambda }_0\right),\\ {\left[L_M^{(1)}\right]}_{SU(N_f)}& ={\left[L_M^{(1)}\right]}_{U(N_f)}-\frac{1}{8N_f(4\pi {)}^2}\left({\lambda }_0+1\right),\end{array}$$

$$\underset{m^{\mathrm{v}}\to m_{\mathrm{cr}}}{\lim }{m}_{\mathrm{pcac}}^{\mathrm{v}}\equiv \underset{m^{\mathrm{v}}\to m_{\mathrm{cr}}}{\lim }\frac{{\partial }_0⟨A^0(x)P(y)⟩}{2⟨P(x)P(y)⟩}=0,$$

$$F_{\pi }=\frac{2{\mu }_0⟨0⟩P{\pi }_{\text{bare}}}{M_{\pi }^2}.$$

$$a{\mu }_0$$

$$a{F}_{\pi }$$

$${\mu }_0=Z_A{m}_{\mathrm{pcac}}^s,\text{ for }{M}_{\pi }^{\mathrm{v}}=M_{\pi }^s.$$

$$\braket{t^{2}E(t)}\Big{|}_{t=t_{0}}=c=0.3.$$

$$⟨t^{2}E(t)⟩=\frac{3}{128{\pi }^2}\frac{N_c^2-1}{N_c}{\lambda }_{\text{GF}}\left(q\right),$$

$$\braket{t^{2}E(t)}\Big{|}_{t=t_{0}}=c(N_{c})={3\over 8}\frac{N_{c}^{2}-1}{N_{c}}c(3).$$

$$\displaystyle\begin{split}\sqrt{t_{0}}\Big{|}^{N_{f}=2}_{M_{\rm ref}}=0.1470(14)\text{ fm,}\\ \sqrt{t_{0}}\Big{|}^{N_{f}=3}_{M_{\rm ref}}=0.1460(19)\text{ fm}\end{split}$$

$$\displaystyle\sqrt{t_{0}}\Big{|}^{N_{f}=4}_{M_{\rm ref}}=0.1450(39)\text{ fm}.$$

$$(M_{\pi}\sqrt{t_{0}})\Big{|}_{M_{\rm ref}}=0.3091(83).$$

$$\displaystyle\braket{t^{2}E(t)}_{a}=\braket{t^{2}E(t)}_{\text{imp}}\Bigg{[}$$

$$\begin{array}{r}{C}_2=-0.262333,C_4=0.0936935,\\ C_6=-0.048002,C_8=\mathrm{0.0320211.}\end{array}$$

$$t_0=t_0^{\chi }\left(1+k{M}^2\right)+\mathcal{O}(M^4),$$

$$M_{\pi }(L)$$

$$F_{\pi }(L)$$

$${\bar{g}}_1(x)\stackrel{x\gg 1}{\to }\frac{24}{x}{K}_1(x)\sim \frac{24\sqrt{2}}{\sqrt{\pi }}\frac{e^{-x}}{x^{3/2}}.$$

$${\mu }^2=\frac{3}{N_c}(4\pi F{)}^2,$$

Full text

The ’t Hooft limit of QCD ’t Hooft (1974) 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 Bardeen et al. (1986, 1987); Chivukula et al. (1986); Sharpe (1987); Pich and de Rafael (1996); Hambye et al. (2003); Buras et al. (2014); Gisbert and Pich (2018); Peláez et al. (2010); Nebreda et al. (2011), 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}$ Lucini and Teper (2001); Bali et al. (2013). Several studies have already been performed. In Ref. Bali et al. (2013) 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. DeGrand and Liu (2016) 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 Lucini and Panero (2013) for a review). In this context, there has been significant progress regarding the properties of mesons García Pérez et al. (2018, 2016); González-Arroyo and Okawa (2016); Hietanen et al. (2009); Narayanan and Neuberger (2005).

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 Appelquist et al. (2014); DeGrand et al. (2015); Aoki et al. (2017); Appelquist et al. (2016); Hansen et al. (2018); Nogradi and Szikszai (2019) and for recent reviews see Brower et al. (2019); Del Debbio (2018).

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 Donini et al. (2016). 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$ Dowdall et al. (2013); Frezzotti et al. (2009); Baron et al. (2010a); Brandt et al. (2013); Gülpers et al. (2015); Bazavov et al. (2010); Beane et al. (2012); Borsanyi et al. (2013); Dürr et al. (2014); Boyle et al. (2016); Baron et al. (2010b); Aoki et al. (2009); Noaki et al. (2008); Baron et al. (2010c); Drach et al. (2018), whose main results are summarized in Aoki et al. (2019).

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. Descotes-Genon et al. (2000); Bernard et al. (2012).

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 results 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, $SU(N_{f})_{L}\times SU(N_{f})_{R}\rightarrow 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, $\lambda=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 $\eta^{\prime}$. 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 $\eta^{\prime}$ 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 Di Vecchia and Veneziano (1980); Rosenzweig et al. (1980); Witten (1980); Kawarabayashi and Ohta (1980); Herrera-Siklody et al. (1997); Kaiser and Leutwyler (2000). 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. Herrera-Siklody et al. (1997); Kaiser and Leutwyler (2000):

[TABLE]

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.

I.1 $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 Bijnens and Lu (2009, 2011); Bijnens et al. (2013):

[TABLE]

and

[TABLE]

in terms of the LO couplings, $B,F$, and the NLO Gasser-Leutwyler coefficients, $L^{r}_{4,5,6,8}(\mu)$, defined at the renormalization scale $\mu$.

Eqs. (2-3) are valid for an arbitrary number of colours, but the LECs scale with $N_{c}$ as (for a review see Manohar (1998)):

[TABLE]

Loop corrections are suppressed in $1/F_{\pi}^{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

[TABLE]

and

[TABLE]

Note that according to the scaling of Eq. (5) and the definition of Eq. 8:

[TABLE]

The NNLO Lagrangian of the $SU(N_{f})$ theory is also known Bijnens and Lu (2009, 2011); Bijnens et al. (2013). At this order we will instead only use the $U(N_{f})$, to which we now turn.

I.2 $U(N_{f})$ effective theory

In the $U(N_{f})$ ChPT at NLO, i.e., ${\mathcal{O}}(\delta^{1})$, the result can be read from Eqs. (2) and (3) and the different $N_{c}$ scalings of the LECs in Eqs. (6) and (8):

[TABLE]

and

[TABLE]

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. Guo et al. (2015). For degenerate flavours, they simplify to:

[TABLE]

and

[TABLE]

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 Guo et al. (2015). Note that for degenerate quarks, there is no $\eta$-$\eta^{\prime}$ mixing.

The $\eta^{\prime}$ mass in this expression can be taken in the large $N_{c}$ limit, where it is given by the Witten-Veneziano formula:

[TABLE]

where $\chi_{t}$ is the topological susceptibility in pure Yang-Mills, recently computed in the large $N_{c}$ limit in Ref. Cè et al. (2016).

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 $\eta^{\prime}$. The matching of the two theories starts at NNLO Kaiser and Leutwyler (2000); Herrera-Siklody (1998) and only affects the coupling $B$ and $L_{M}^{(1)}$ of the above Herrera-Siklody (1998); Kaiser and Leutwyler (2000):

[TABLE]

with $\lambda_{0}=\log\frac{M_{0}^{2}}{\mu^{2}}$.

I.3 $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}\rightarrow\infty$ limit is $N_{f}$ independent and the NLO is ${\mathcal{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_{\pi}$. 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_{\pi}\ll M_{\eta^{\prime}}$, that is when the chiral limit is taken first, anomalous $1/N_{f}$ terms appear coming from an expansion in $M_{\pi}/M_{\eta^{\prime}}$. 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 ${\mathcal{O}}(N_{f}/N_{c})$ dependence, as can be explicitly seen in the matching of $L_{M}^{(1)}$ in Eq. (15).

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})}$.

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 Del Debbio et al. (2010). Some of them have been already presented in Ref. Romero-López et al. (2018). We have chosen the Iwasaki gauge action (following previous experience with 2+1+1 simulations Finkenrath et al. (2018); Alexandrou et al. (2018)) and $O(a)$-improved111For $N_{c}=3$, we take the perturbative value of $c_{sw}=1+c^{(1)}_{sw}g^{2}$ from Ref. Aoki and Kuramashi (2003), where we use the plaquette-boosted coupling $g^{2}=2N_{c}/(\beta 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. Aoki and Kuramashi (2003)). 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. Aoki and Kuramashi (2003). 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 $\sim 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 Bar et al. (2003) previously used in Refs. Herdoíza et al. (2018); Bussone et al. (2018a, b). Maximal twist is ensured by tuning the untwisted bare valence mass, $m^{\rm v}$ to the critical value for which the valence PCAC mass is zero:

[TABLE]

where $A^{\mu}(x)\equiv\bar{\Psi}(x)\gamma^{\mu}\gamma^{5}\Psi(x)$ and $P(x)=\bar{\Psi}(x)\gamma^{5}\Psi(x)$. The bare twisted mass parameter $\mu_{0}$ is tuned such that the pion mass in the sea and valence sectors coincide, $M_{\pi}^{\rm v}=M_{\pi}^{s}$. The normalized meson decay constant $F_{\pi}$ can then be obtained from the bare combination Shindler (2008):

[TABLE]

The results for the meson masses and decay constant in the mixed-action setup can be seen in Table 2. We have achieved a good tuning of $m_{\rm pcac}$ and the pseudoscalar masses are compatible within one or two sigma with their pure Wilson value (see Table 1). When the tuning to maximal twist is not perfect, we correct the bare quark mass (and thus $F_{\pi}$) as follows (see also Shindler (2008)):

[TABLE]

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:

[TABLE]

III Scale setting at Large $N_{c}$

The scale setting for different values of $N_{c}$ is performed using the gradient flow scale $\sqrt{8t_{0}}$, via the determination of $t_{0}/a^{2}$. In QCD, with $N_{c}=3$, the standard definition of $t_{0}$ is:

[TABLE]

The leading dependence in $N_{c}$ is known Lüscher (2010) in perturbation theory:

[TABLE]

where $\lambda_{\text{GF}}(q)$ is the gradient flow ’t Hooft coupling at the scale $q=1/{\sqrt{8t}}$. Hence, as in Ref. Cè et al. (2016), we will generalize $t_{0}$ to an arbitrary $N_{c}$ as:

[TABLE]

Notice that the choice here is not unique. In particular, one could choose another coupling in a different scheme (such as $\overline{\text{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$ Bruno and Sommer (2014); Sommer (2014) and $N_{f}=3$ Bruno et al. (2017) degenerate quarks and at a reference pion mass $M_{\rm ref}=420\text{ MeV}$:

[TABLE]

We can use these to perform a linear extrapolation to $N_{f}=4$, motivated by the weak $N_{f}$ dependence:

[TABLE]

Our scale setting condition involves therefore the dimensionless quantity

[TABLE]

In order to reduce discretization errors we have performed a tree level improvement of $t_{0}$. In Ref. Fodor et al. (2014), lattice perturbation theory is used to improve $\braket{t^{2}E(t)}$ and thus, $t_{0}$. The prescription is:

[TABLE]

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:

[TABLE]

The numerical results after the improvement, $t_{0}^{\rm imp}/a^{2}$, are shown in Table 1.

Finally, Eq. (26) requires $t_{0}$ at $M_{\rm ref}$. The mass dependence of $t_{0}$ has been studied in chiral perturbation theory in Ref. Bar and Golterman (2014). For degenerate flavours it is given by

[TABLE]

where $k\propto 1/(F_{\pi})^{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^{\text{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 3, 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. Fodor et al. (2012). 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\text{ fm}$, corresponding to $N_{c}=5$.

IV Chiral Perturbation Theory Fits

The results for $M_{\pi}$ and $F_{\pi}$ in the mixed-action setup are presented in Table 2. 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 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. 6 and 8. 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.

IV.1 Finite volume effects

Our ensembles have $M_{\pi}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 Gasser and Leutwyler (1987); Colangelo et al. (2005):

[TABLE]

with $\xi\equiv\frac{M_{\pi}^{2}}{(4\pi F_{\pi})^{2}}$, while $\bar{g}_{1}(x)$ is given by

[TABLE]

We will use the corrected results for the analysis.

IV.2 Renormalization scale

The NLO couplings are usually defined at $\mu=4\pi F$ or at the $\rho$ mass, $\mu=M_{\rho}$. Still, in the context of the large $N_{c}$ expansion these are two very different choices, since the former scales with $\sqrt{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 $\rho$. The scale $\mu=4\pi 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\pi F$, however, has some advantages for fitting, so we choose:

[TABLE]

which has no leading dependence on $N_{c}$. Using this scale, the NLO predictions can be conveniently written as:

[TABLE]

where $m=\mu_{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.

IV.3 Fitting strategy

Some care is needed to perform the fits in Eqs. (34) and (35). The complication comes from the fact that both coordinates, $(x,y)=(\xi,F_{\pi})$ or $(x,y)=(\xi,M^{2}_{\pi}/\mu_{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) York et al. (2004), in which the $\chi^{2}$ function is:

[TABLE]

where we have defined the two-dimensional vectors:

[TABLE]

where $f$ is the fitting function, and $V$ is the $x,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 $\sim 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.

IV.4 Fit results at fixed $N_{c}$

First we consider each $N_{c}$ separately and perform a fit of the data points to extract $F,L_{F}(\mu)$ and $B,L_{M}(\mu)$. The NLO fit results for these quantities are shown respectively in Tables 4 and Tables 5. The $N_{c}$ dependence of the LECs is shown in Figs. 2 and  3. It can be seen that the scaling is well described by leading and subleading $N_{c}$ corrections for $N_{c}=4-6$, while there seems to be significant $1/N_{c}^{2}$ corrections for $N_{c}=3$ in the case of $F$ and $L_{F}$. In the case of $B$ and $L_{M}$ errors are larger and there is no sign of $1/N_{c}^{2}$. Interestingly, the data suggest that the large $N_{c}$ limit of $L_{M}\sim 0$.

IV.5 Simultaneous chiral and $N_{c}$ fits

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. (6) and (8). We linearize the fit by considering the following parametrization

[TABLE]

where $(FL_{F})^{(0)}\equiv F_{0}L_{F}^{(0)}$, while $(FL_{F})^{(1)}\equiv F_{1}L_{F}^{(0)}+L_{F}^{(1)}F_{0}$, and $(BL_{M})^{(0)}\equiv B_{0}L_{M}^{(0)}$ and $(BL_{M})^{(1)}\equiv B_{1}L_{M}^{(0)}+B_{0}L_{M}^{(1)}$.

Secondly, we consider the $U(N_{f})$-NNLO expansion, since we have checked that the $U(N_{f})$-NLO expressions fit the data very poorly. We also linearize the fit by considering the following fitting functions:

[TABLE]

[TABLE]

where

[TABLE]

and $M_{0}^{2}$ is given by the Witten-Veneziano formula for the $\eta^{\prime}$ mass valid in the large $N_{c}$ limit (see Eq. 14). We use the result for the topological susceptibility from Ref. Cè et al. (2016),

[TABLE]

We convert to lattice units using the value of $t_{0}/a^{2}$ in the previous section and substitute $F\rightarrow\sqrt{N_{c}}F_{0}$, as extracted from the global $F_{\pi}$ fit. We find $a_{0}\sim 6.5$, a value we fix in the fit.

In summary we compare the following fits:

1. i)

Fit 1: $SU(N_{f})$-NLO fit to Eq. (38) and (39) including the data subset $N_{c}=4$-6. 2. ii)

Fit 2: $U(N_{f})$-NNLO expansion fit to Eqs. (IV.5) and (IV.5) including the full data set.

The results for the fitted parameters in the global fits are shown in Tables 6 and 7, and the quality of the fits is shown in Figs. 4(a) and 4(b). We also quote in Table 8 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.

IV.6 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$:

[TABLE]

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_{c}^{2}$ terms. The first error is just the one obtained from the fits in Table 6 and the second error of $3\%$ is the one corresponding to the lattice spacing determination. For two- and three-flavour QCD we get:

[TABLE]

where we have taken into account the correlations between the different terms in Eq. (44), 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:

[TABLE]

and also lattice results (see Ref. Aoki et al. (2019)). In addition, we can compare to previous results in the large $N_{c}$ limit in the quenched approximation:

[TABLE]

This value is $2\sigma$ away from the results in Eq. (44). 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\equiv\Sigma/F^{2}$, we do not have a non-perturbative value of $Z_{P}$, up to this factor we get:

[TABLE]

From Ref. Alexandrou et al. (2012), we can obtain the the 1-loop perturbative result for the normalization constant:

[TABLE]

which at the order we are working is independent of $N_{f}$. With this, we obtain for $N_{c}=3$:

[TABLE]

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

[TABLE]

which is compatible within $1\sigma$ with the numbers quoted in Ref. Aoki et al. (2019). 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):

[TABLE]

which shows good agreement with the prediction

[TABLE]

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

[TABLE]

while for the NLO LEC for the mass, we can only give the $N_{c}$ scaling at $N_{f}=4$:

[TABLE]

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

[TABLE]

From Eqs. 58 and 60, we can infer the $N_{c}=3,N_{f}=3$ results:

[TABLE]

For $N_{c}=3,N_{f}=2$, it is more common to quote ${\bar{\ell}}_{3}$ and ${\bar{\ell}}_{4}$:

[TABLE]

This way, we obtain:

[TABLE]

We stress that $U(N_{f})$ $\bar{\ell}_{3}$ in fit 2 is not the same as the standard ${\bar{\ell}}_{3}$ in $SU(N_{f})$. ${\bar{\ell}}_{4}$ agrees instead at $1-2\sigma$ with the results quoted in Ref. Aoki et al. (2019).

IV.7 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_{\pi}$ depend upon the mixed pion mass Chen et al. (2007):

[TABLE]

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:

[TABLE]

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_{\pi}]_{c_{sw}=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 $\sim 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. in $N_{c}$

We 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.

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. (44), (49) and Table 8). 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}=2$ and $N_{f}=3$ at $N_{c}=3$. We find that our results nicely agree with those in the literature regarding $L_{F},L_{M}$ and $F$ (see for example Ref. Aoki et al. (2019) for a summary of results). For $B$ we need to improve our determination, including a non-perturbatively 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. Descotes-Genon et al. (2000); Bernard et al. (2012).

We would like to stress that the results presented in this paper are complementary to similar studies that can be performed in reduced models García Pérez et al. (2018, 2016); González-Arroyo and Okawa (2016); Hietanen et al. (2009); Narayanan and Neuberger (2005) or the quenched approximations at large $N_{c}$ Bali et al. (2013), since both of these approaches must yield the leading order result as $N_{c}\to\infty$. Given the strong correlations presents in our results (see Fig. 5), a precise determination of the dominant $N_{c}$ term would significantly improve the determination of the subleading $N_{c}$ corrections, and hence the determination of the physical values at $N_{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 the $K\to\pi$ matrix elements (see Donini et al. (2016); Romero-López et al. (2018) 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$ $\pi\pi$ 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 Peláez et al. (2010); Nebreda et al. (2011); Bernard et al. (2011); Ruiz de Elvira et al. (2017) .

Acknowledgements.

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.