Lattice study of the chiral properties of large N c QCD

We present a lattice calculation of the low energy constants of QCD with N c = 3, 4 and 5 colors and N f = 2 ﬂavors of degenerate mass fermions. We ﬁt data for the pseudoscalar meson mass, the pseudoscalar decay constant, and the Axial Ward Identity fermion mass to formulas from next to next to leading order chiral perturbation theory. We extract the next to leading order low energy constants and study their behavior as a function of N c . Pre-existing analyses of N c = 3 inform our ﬁtting strategies.


I. INTRODUCTION AND MOTIVATION
The most studied deformation of real-world three-color (N c = 3) QCD is its extension to a large number of colors.In the limit that N c is taken to infinity, it is believed that mesons are narrow quark anti-quark bound states, baryon masses scale as N c , and N c counting rules set the overall scale for hadronic matrix elements [1][2][3].
There is a small lattice literature related to the large N c limit of QCD.Simulations involve the pure gauge theory and fermions included in quenched approximation (mostly at single-digit values of N c ), quenched simulations at very large N c and small volume, and simulations in large volume with a small number of dynamical fermion flavors (so far, N f = 2 and 4).Refs.[4][5][6] are a selection of reviews of the various approaches.The goal of these simulations is to confront predictions of large N c QCD, which are typically (semi) analytic, with nonperturbative lattice-based results.
The qualitative situation with regard to such comparisons is as follows: Choose to fix the lattice spacing in a simulation using some observable taken from a correlation function dominated by gluonic degrees of freedom.Such observables include the string tension, the Sommer parameter [7], or the flow [8,9] parameter t 0 .Perform lattice simulations while tuning the bare parameters (the most sensitive one is the gauge coupling g 2 ) so that the lattice spacing is the same for all values of N c .One will discover that the bare 't Hooft couplings λ = g 2 N c are roughly matched across N c .Alternatively, performing simulations across N c at the same value of λ, one will discover that the lattice spacing had been approximately matched.Once that is done one will see that the masses of flavor nonsinglet mesons, and their dependence on the fermion mass, will be approximately equal across N c .Simple matrix elements (decay constants, the kaon B− parameter) will scale with N c as predicted from the color weight of amplitudes.All of these statements are modified by small corrections, which can be arguably interpreted as 1/N c effects.
We are not aware of any large volume simulations which attempt to take a continuum limit; most of them involve simulations across N c at a single matched value lattice spacing.It seemed to us that it might be possible to extend one of these studies to do that.We make a first attempt to measure at chiral parameters of SU(N c ) gauge theory with N f = 2 flavors of degenerate fermions, in the continuum and large N c limits.Specifically, we studied N c = 3, 4 and 5.We wanted to address the following set of questions: 1. How do the low energy constants of the chiral effective theory scale with N c ? 2. What is their N c → ∞ limit?How do our results compare to those from other approaches (quenched simulations, simulations in small volume)?
3. As N c rises, the low energy effective theory is expected to change from the conventional SU(N f ) L × SU(N f ) R → SU(N f ) L+R pattern of symmetry breaking to one where the flavor singlet meson becomes light [10][11][12][13][14][15].At what value of N c does this "U(N f )" expansion begin to reproduce the data?
There is already an extensive literature about the low energy constants of N c = 3, N f = 2 QCD, summarized in the Flavor Lattice Averaging Group's 2019 summary, Ref. [16].Typical analyses in their summary have much better quality than what we present here; we cannot add much to this already-closed subject.We instead use the results of these pre-existing N c = 3, N f = 2 to inform our analysis, and to calibrate the accuracy we claim for N c > 3.
We should say at the start that while we have made a first attempt to extract continuum to the low energy constants of N c > 3 QCD, we do not regard our results as definitive.The main issues are three: two typical of early QCD studies of any observable, with too large a lattice spacing and too large fermion masses for comfortable extrapolations, the third is that we do not feel that we have enough values of N c .But the outline of the solution is there: once again, N c = 3 is not that different from infinite N c .
Our calculation was strongly influenced by two previous ones.The first is a 2013 computation of the low energy parameters of N f = 2 QCD (N c = 3) in the continuum limit, by the Budapest-Marseille-Wuppertal collaboration [17].Their study has larger volumes, smaller lattice spacings, and smaller fermion masses, and than our exploratory project.We adapted much of their analysis methodology.The other project was a calculation of the low energy parameters of N f = 4 systems by P. Hernández, C. Pena and F. Romero-López [18].This was done at a single lattice spacing, matched across N c .They studied N c = 3, 4, 5, and 6.The addition of N c = 6 allowed for more controlled extrapolation to infinite N c than we could do.This is because large N c predictions are typically involve a power series in N c : and it is useful to have more than j + 1 values of N c to fit data to such predictions at jth order.
Most discussions of the literature of large N c QCD (such as the one we just gave) make the point that comparisons across N c reveal only small differences.But to lattice practitioners QCD's at different N c 's are different.Some of these differences are unfavorable.For example, the up-front cost of a simulation at some N c scales roughly like N 2 c (when matrix times vector multiplication dominates, a situation encountered in inverting the Dirac operator) or N 3 c (for matrix multiplication of gauge links).The simulation autocorrelation time as measured by the topological susceptibility grows with N c , in both quenched simulations and ones with dynamical fermions.However, differences across N c can be favorable for larger N c .That is the case in this project, where we had much greater difficulty generating useful N c = 3 data sets than we did with N c = 4 or 5.This is due to the scaling of the pseudoscalar decay constant F with N c : F ∝ √ N c .This is helpful in two (related) ways.First, the expansion parameter of chiral perturbation theory is (speaking loosely) the ratio of the squared pseudoscalar mass to F 2 , x = M 2 P S /(8π 2 F 2 ).Simulations which make comparisons at identical x across N c involve ever larger values of M 2 P S and become easier to perform as N c grows.Second, finite volume corrections in a box of size L are proportional to ∆(m, L)/F 2 where ∆(m, L) is the propagator of a pseudoscalar"to an image point," which appears along with the "tadpole" or "'snail" graph in finite volume.The 1/F 2 ∝ 1/N c prefactor means that as N c grows there are smaller finite volume corrections at identical pion masses.As a consequence, at bigger N c one can simulate at smaller volume without encountering volume artifacts.
The outline of the paper is as follows: Sec.II sets the conventions we follow for the low energy chiral effective theory, to which we compare our data.Secs.III and IV describes the various parts of the lattice simulation, first the collection of data and its conversion to the dimensionless quantities which are what we fit, and then a discussion of issues involved in the fits themselves.Results of fits which extract the low energy constants of the chiral effective Lagrangian are presented in Sec.V. Some conclusions are found in Sec.VI.

II. THEORETICAL BACKGROUND AND CONVENTIONS FOR CHIRAL OB-SERVABLES
The goal of this project is to compare the relationship between quantities relevant to the low energy properties of SU(N c ) gauge theory with two flavors of degenerate mass fundamental representation fermions, across N c and in the continuum limit.We focus on the fermion mass, the (squared) pseudoscalar mass and pseudoscalar decay constant as measured in simulations.We label these quantities as m q , m 2 P S , and f P S .The comparisons are done in the context of the parameters of a low energy effective chiral Lagrangian.There are two choices for an effective theory, which differ in the number of light degrees of freedom they represent.
The effective field theory for N c = 3 QCD is based on the spontaneous breaking of chiral symmetry SU(N f ) L × SU(N f ) R → SU(N f ) L+R (plus its explicit brealing by quark masses).Its small expansion parameter is where M 2 is the squared mass of the fields in the Lagrangian (taken to be massless Goldstone bosons in the m q → 0 limit).We will refer to this effective theory as "SU(N f ) effective theory."(Of course, for us, N f = 2.)Only flavor nonsinglet mesons appear in the Lagrangian of SU(2) effective theory.In particular, the flavor singlet meson, the eta-prime (in N f = 3 language) gets a large mass from the anomaly and does not appear in the low energy theory of the N c = 3 world.
However, the situation at large N c is different.The size of anomaly term falls with N c and the squared eta prime mass decreases as 1/N c , so at some point the eta prime must be included in the low energy effective theory.The resulting theory is called U(N f ) (here U(2)) effective theory and has a slightly different chiral expansion.Which chiral expansion scheme works better for any particular value of N c is an open question, which we hope to explore.
We begin by discussing SU(2) chiral perturbation theory.We follow the usual convention for lattice simulations and define the pseudoscalar decay constant f P S in terms of the matrix element 0|ūγ 0 γ 5 d|π = m P S f P S . ( This leads to the identification f π ∼ 132 MeV in QCD.Note that the continuum chiral PT literature uses a "93 MeV" definition: the extra √ 2 is an isospin raising factor.See Ref. [19] for a compilation of conventions for this quantity. There are two commonly used expansion schemes for SU(N f ) chiral perturbation theory.The two leading order constants are B and F .We will mostly use the "x-expansion" which where Our analysis convention is to fix a M = −1/2 and a F = 1 in any fit.Then at NLO there are four LEC's which can be measured, B, F , l 3 and l 4 .NNLO adds an additional three LEC's, (l 12 , k M , k F ) to the collection of quantities to be fit.The LEC's are expected to show the following N c scaling (the quick list can be found in in Ref. [18]): This means that the l i 's for i = 3, 4 scale as i .
Our analysis of data using SU(N f ) chiral perturbation theory is done with separate fits for each N c .We can then ask whether the LEC's scale with N c as given by Eqs.12-13.
The "ξ-convention" rewrites the chiral expansion in terms of a ratio of observables (in the "132 MeV" definition for the decay constant).We have only implemented the NLO expressions, which are and The ξ parameterization is often used because the coefficients of higher order terms (here order( ξ 2 )) are smaller than in the x parameterization.It is more awkward to use in a fit since the expressions have to be inverted (to express measured quantities such as F P S ±∆F P S in terms of fit parameters).At NLO, the relevant expressions are and The U(N f ) chiral Lagrangian includes the eta prime in its degrees of freedom.The appropriate effective field theory was developed in Refs.[10][11][12][13][14][15].The power counting is (Notice that 1/F 2 is O(δ) and the SU(N f ) chiral expansion parameter x of Eq. 4 is O(δ 2 ).For f P S and m 2 P S the many LEC's of the complete theory reduce to two, each with its own expansion in N c , as given by Eq. 13 as in the SU(N f ) case.The factors of N c in the various l i 's of the SU(N f ) expressions give O(1/δ) scaling factors in the U(N f ) perturbative expansion.Thus, all that survives from the x 2 terms T M and T F in Eq. 11 are constants proportional to N 2 c .The eta prime enters as an additional Nambu-Goldstone boson of the non-anomalous (in the large N c limit) singlet U(1) A symmetry.It appears in the perturbative expansion as an extra tadpole contribution.The perturbative calculations for F P S and M 2 P S can be found in Ref. [21], are succinctly presented for the degenerate-flavor case in Ref. [18], and can be re-assembled (if desired) by slightly editing the calculation of Ref. [22].
The formulas for the chiral expansion need the mass of the eta prime and its decay constant as a function of N c .The calculation of the eta prime mass for N c = 3 is already a difficult problem (the statistics of a recent study [23] dwarf those of our little project) and so we have recourse to theory.Witten and Veneziano [24,25] related the eta prime mass to the quenched topological susceptibility χ T and the pseudoscalar decay constant F (recall the 132 MeV convention), where the zero-quark mass eta prime mass is Note that M 2 0 is also O(δ).There is indirect lattice evidence in favor of this relation from the study of the N f = 2 topological susceptibility across N c from Ref. [26].A phenomenological formula which relates the quenched topological susceptibility to the topological susceptibility at small quark mass due to Refs.[10], [27] and [28] agreed qualitatively with lattice results.
We thus proceed, taking Eq.20 for the eta prime mass, assuming that the eta prime and the pions share a common decay constant, and defining and We use Eqs.23-25 for fits to the data at an individual N c value.Notice that, in that case, the NLO fits give l (0) i , but the NNLO fits can only give the full i .We can present results for the LEC's as we do for the SU(N f ) case, a plot of, for example, F versus 1/N c .
We can also imagine doing combined fits to the data.In that case we would make an ansatz, that We would insert these expressions in Eqs.24-25, combine them systematically along with the other LEC's, and fit a range of N c values.In order to do this, of course, the U(N f ) chiral expansion has to be well behaved for all the individual N c values included in the fit.We return to this point when we present fit results in Sec.V, but for now we remark that our N c = 3 data is incompatible with the U(2) expansion leaving only two N c 's to work with; too few values of N c to think about NNLO fits.So we did not pursue this line.We note that the LEC's of SU(N f ) chiral perturbation theory and U(N f ) chiral perturbation theory are not identical [15,29].In particular, there are shifts The shift in the B ′ s is formally order 1/N 2 c but the prefactor happens to be large.Note also that, in contrast to the SU(N f ) case, in the U(N f ) expansion a shift in the regularization point µ affects the value of B. FLAG [16] quotes values for the LEC's using µ 2 set to the physical squared pion mass.This choice causes a large shift between the NLO and NNLO determinations of B at low N c .To avoid this large logarithm, we choose a regularization point closer to the eta prime mass; we take µ 2 to be an N c − independent quantity which we choose to be µ 2 = 8π 2 f 2 π where f π = 132 MeV in all our comparison of lattice data to the U(N f ) expansion.When we compare the U(2) LEC's to the SU(2) ones, there will also be a shift in the SU(2) l ′ i s due to the different choices for µ, see Eq. 8.

III. DETAILS OF THE NUMERICAL SIMULATIONS
The project involves simulations at three values of N c .Each N c value requires its own set of simulations at four bare gauge couplings and 4-6 bare quark masses per gauge coupling.At each bare parameter value six quantities are measured: the Axial Ward Identity fermion mass, the pseudoscalar mass, the pseudoscalar decay constant, two lattice to continuum renormalization factors (for the fermion mass and decay constant), as well as a parameter to set the lattice spacing.All these quantities must be checked for correlations.The lattice regulated quantities are converted into dimensionless continuum regulated quantities which are then fit to the formulas described in the last section, plus lattice artifacts.Doing this is a long and somewhat recursive task.We believe that the methodology we use to generate, collect, and analyze our lattice data sets is completely standard.However, all lattice calculations involve many choices.This section describes what we did in perhaps too much detail.
A. Simulation details: lattice action, updating algorithm, data sets, observables Our simulations are done with the usual Wilson plaquette action, with the bare gauge coupling g 0 labelled by β = 2N c /g 2 0 .Two flavors of degenerate mass fermions (discretized as Wilson-clover fermions) are included.Configurations are generated using the Hybrid Monte Carlo (HMC) algorithm [30][31][32] with a multi-level Omelyan integrator [33] and multiple integration time steps [34] with one level of mass preconditioning for the fermions [35].
The fermion action uses gauge connections defined as normalized hypercubic (nHYP) smeared links [36][37][38].Simulations use the arbitrary N c implementation of Ref. [39].The bare quark mass m q 0 appears in the lattice action along with the lattice spacing a via the hopping parameter κ = (2m q 0 a + 8) −1 .The clover coefficient is fixed to its tree level value, c SW = 1.The gauge fields obey periodic boundary conditions; the fermions are periodic in space and antiperiodic in time.
Lattice volumes are a mix of 16 3 × 32 and 24 3 × 32 sites.The lattice sizes were chosen so as to minimize finite volume effects (which we describe in Sec.III C below).Briefly, we began simulation runs at nearly all our bare parameter values taking 16 3 spatial volumes and collected enough data to determine the mass and decay constant of a pseudoscalar meson, and then to estimate the size of the finite volume correction.When this appeared to be larger than a few per cent, we replaced the 16 3 data sets by 24 3 ones.A glance at Tables I-III shows that there are more 24 3 N c = 3 data sets than there are for N c = 4 or 5.This is because of the scaling of the pseudoscalar decay constant with N c , as described in the Introduction.
Lattices used for analysis are spaced a minimum of 10 HMC time units apart, so individual bare parameter sets contain 30 to 200 stored lattices.All data sets at individual bare parameters (β, κ) are based on a single stream.We check for thermalization at the start of each data stream (which typically begins from an equilibrated configuration at a nearby bare parameter set) and typically discard the first 100 trajectories (more than this for simulations at strong coupling).The data sets are extensions of ones presented in Refs.[39,40].
Data for the pseudoscalar mass, the pseudoscalar decay constant, and the fermion mass come from hadron correlators using propagators constructed in Coulomb gauge, with Gaussian sources and p = 0 point sinks.We tune the width of the source to minimize the dependence of the effective mass (defined, as usual, to be m ef f (t) = log C(t)/C(t + 1) in the case of open boundary conditions for the correlator C(t)) on the distance t between source and sink.All results are based on a standard full correlated analysis involving fits to a wide range of t's.Best fits are chosen with the "model averaging" ansatz of Jay and Neil [41].In one sentence, this gives each particular fit in a suite (with a chi-squared value χ 2 and N DOF degrees of freedom) a weight in the average which is proportional to exp(−(χ 2 /2 − N DOF )).
The pseudoscalar decay constant is defined in terms of the matrix element of the timelike component of the axial current A 0 (x, t) = ū(x, t)γ 0 γ 5 d(x, t) between the vacuum and a pseudoscalar state as in Eq. 3. It is determined via a correlated fit to two propagators with a common source, The fermion mass am q is taken to be the Axial Ward Identity (AWI) fermion mass, defined through the relation of the axial current matrix element to the pseudoscalar again taken from a simultaneous fit to two two-point functions.O can be any convenient source, of course, and for us it is a Gaussian.Conversion factors from lattice regulated quantities will be needed for f P S and m q and their calculation will be described in Sec.III D below.
The pseudoscalar mass am P S is determined from fits involving correlators of pseudoscalar currents.
We give is a brief summary of the data sets.We attempted to collect data at similar values of ξ = m 2 P S /(8π 2 f 2 P S ) for our different values of N c , since that is the expansion parameter for chiral perturbation theory.
For N c = 3 we have four beta values and 24 bare parameter sets in total.The ξ range is from 0.06 to 0.19.t 0 m 2 P S ranges from 0.05 to 0.35, for pion masses in the range 300 -776 MeV (taking a nominal value for the flow parameter √ t 0 = 0.15 fm from Ref. [42]).t 0 /a 2 ranges from 1.0 to 2.4 for lattice spacings between 0.10 -0.15 fm.
The N c = 4 data sets come from four beta values and there are 26 bare parameter sets.ξ ranges from 0.04 to 0.2 and t 0 m 2 P S is in the range 0.06-0.67.t 0 /a 2 is in the range 1.1-3.34 for a = 0.08 − 0.14 fm.
And there are four beta values and 21 bare parameter sets in the N c = 5 sector.ξ ranges from 0.05 to 0.16, t 0 m 2 P S is in the range 0.10-0.63,while t 0 lies between 1.5 and 3.4, or a = 0.08 − 0.12 fm.
A useful marker point is ξ = 0.1, for which t 0 m 2 P S is about 0.10, 0.22, and 0.35 for N c = 3, 4, or 5.This quantity will re-appear in Sec.V.
The data are presented in Tables I-III.

B. Setting a scale
The conversion of simulations quantities whose scale is set by the lattice spacing into dimensionless parameters which can be used in continuum extrapolations is done using the Wilson flow parameter t 0 [8,9].The determination of t 0 from each data set (with its own β and κ) is done in the standard way, from the observable E(t 0 ) extracted from the field strength tensor, C(N c ) is chosen to match what most other large N c simulations take, with C = 0.3 the usual value used in SU(3).Our procedure for computing t 0 from the data is identical to what was done in Ref. [26] and details may be found there.
In contrast to our earlier work, where conversion from dimensionful to dimensionless parameter ((am P S ) 2 to t 0 m 2 P S , for example) was done at individual bare parameter values with a t 0 (β, κ), we need to use a mass-independent definition of t 0 .This is because the chiral expansion for t 0 has its own set of NLO analytic contributions which combine with the ones for our desired observables (f P S and m 2 P S ) [43] and would contaminate the chiral LEC's.We also need a scale choice which is superficially identical across N c .
We investigated several possibilities • t 0 at a fixed value (N c /3)ξ -This is a bit awkward to produce, since it needs Z A , the axial current matching factor.
• t 0 at fixed m 2 P S /m 2 V (where m V is the mass of the vector meson) -This was noisier than the first choice.
• t 0 at fixed t 0 m 2 P S -This was also noisier than the first choice (and with greater correlations).
Then we followed the following procedure to produce a value of t 0 at each β value.Each data set gives it own ξ 0 ± ∆ξ 0 and t 0 ± ∆t 0 .For each β value, we performed a polynomial fit t 0 (i) = n C n (ξ(i) − ξ 0 ) n to produce a value of t 0 (ξ 0 , β).The fit included the errors ∆t 0 and ∆ξ.(See Sec.IV for a discussion about fit methodology.)We checked subsets of the data to see if correlations in t 0 and ξ affected the fit.They did not (differences in correlated and uncorrelated fits for t 0 (ξ 0 ) were well under a standard deviation) so we did not include them in our analysis.We typically fit to a cubic, bracketing ξ 0 so that the fit was essentially an interpolation.
Our choice of ξ 0 is purely heuristic.We have to make sure ξ 0 lies inside the range of measured (N c /3)ξ values recorded for each bare parameter set, and then we simply pick a ξ 0 for which ∆t 0 is minimized.Like many other decisions made in the data analysis, we looked at many possible choices but we saw little effect on results.In the end, we took ξ 0 = 0.12.
An example of the extraction of t 0 from data, for N c = 5, is shown in Fig. 1, and the data are collected in Table IV.

C. Correcting for finite volume
Our smallest quark mass data could potentially be affected by our simulation volumes.For chiral observables, the cause is well known: tadpole graphs, where a pseudoscalar particle is emitted from some location and returns to the same point, are replaced by a set of contributions connecting the location to its image points.(We found the discussion in Ref. [44] to be quite helpful.)If we write the pseudoscalar correlator for a particle of mass m in a box of length L µ in direction µ as the infinite volume propagator, call it ∆(m, x), is the n = 0 term in the sum.The finite volume tadpole is where Ī1 (m, L) is the sum over images.If a typical infinite volume observable has a chiral expansion then the finite volume correction is We replace the simulation observable in our data set O(L) by Ī1 (m, L) can be found from by summing over all image points until the expression saturates.Eq. 36 also needs f P S , the decay constant in the chiral limit.This is, of course a parameter in the fit.However, if we just make an estimate for it and evaluate Eq. 36, we find that for all our data sets, the correction is at most on the order of a few per cent for a few of the lightest-mass points, so we only need to input an approximate value of f P S given all the other uncertainties in the calculation.We use √ t 0 f P S = 0.078, 0.107, 0.128 for N c = 3, 4, 5.

D. Lattice to continuum regularization conversion
We will use the "regularization independent" or RI scheme [45] for computing matching factors (labelled as Z i for current i).The method is standard and so we will only briefly describe our particular implementation of it.
One computes quark and gluon Green's functions in Landau gauge, regulated by giving all external lines a large Euclidean squared momentum p 2 = µ 2 , and uses combinations of them to determine the Z's.Specifically, we define where p i = 2n i π/L and p t = (2n i + 1)π/T since the temporal boundary conditions are anti-periodic.(Recall that the lattices are L 3 × T sites.) The matching factor is defined by computing the ratio of the matrix elements of the operator between off-shell single particle momentum eigenstates in the full and free theories or equivalently where the factor of 1/4N c counts Dirac spins and colors.We use the RI ′ scheme to define Z Q from a projection against the free propagator S 0 (p): We are using clover fermions, for which there is an overall 2κ factor in the field definition as opposed to a continuum-like normalization.We define our Z Γ factors for observables built from clover fermions by where Z Γ ∼ 1 + Cg 2 ... in a perturbative expansion.The 2κ factor gives clover Z's equal to unity for free field theory.
One needs the scale µ to be large enough not to be affected by confinement effects but not so large that is is affected by lattice artifacts.Given our lattice spacings and volumes, this means in practice a value of µa around unity.
We need a renormalization factor for the axial current, Z A , to define F and one for the fermion mass Z m to carry the quark mass to a continuum regularization.The computations are straightforward.We generate lattice data for Z A (aµ) at many values of (aµ) 2 = j a 2 p 2 j for lattice momenta ap j .The data is averaged under single elimination jackknife and sorted in aµ.Fits to a linear function are performed over "windows" of data around a central value and checked for systematic dependence on the chosen central value.We typically use about 30 lattices per bare parameter value to do this.The resulting uncertainties are completely dominated by systematics, to be described below, so it is not worthwhile to work with larger data sets.
1. Z A Fig. 2 illustrates our results for Z A .It shows large scale views of Z A for a 16 3 and a 24 3 spatial volume lattice.The 16 3 volume shows clear discretization artifacts above aµ ∼ 1.4.This represents an upper limit on the range of aµ values which can be analyzed.The finer grained 24 3 data do not have this issue.
Fits about some central value of aµ typically show a value close to unity with a statistical uncertainty of a few parts per mille.That Z A is close to unity for the fermion action used here (nHYP) fermions has been well documented in the literature (cf.Ref. [46]).
However, there is another issue.Z A should be a constant across aµ because the operator has no anomalous dimension.Our data sets show a small but noticeable drift of Z A (aµ) versus aµ.We account for this drift by comparing linear fits to Z A (aµ) versus aµ over windows of aµ.We pick (somewhat arbitrarily) a central value aµ = 1 and a range aµ = 0.8 − 1.2.(For comparison, µ = 2 GeV corresponds to aµ ∼ 1.4 at our coarsest lattice spacings down to 0.8 at our finest spacings.)Our systematic uncertainty is taken to be |Z A (aµ = 1.2) − Z A (aµ = 0.8)|/2.This gives a systematic uncertainty in the range of 0.005-0.008as compared to statistical uncertainties from Z A (aµ = 1) of 0.002-0.003.We use the first quantity as the uncertainty for Z A .

Z S and Z m
Z m is a bit more complicated to compute.We take Z m from the renormalization factor for the scalar current Z RI S as Z m = 1/Z S , defined at some value of µa.We must then convert it to an MS value at a nominal scale µ = 2 GeV.(This is done to make contact with pre-existing SU(3) results.)This involves two steps after Z RI S (aµ) is determined; the analysis to determine this quantity is identical to what is done for Z A .First, there is the conversion to MS at the scale at which Z RI S is determined.An explicit three loop formula (for general N c and N f fundamental representation fermions) for this conversion is given by Chetyrkin and Retey [47].(See also Ref. [48].)To use it we need the MS coupling at scale µ = (µa)/a.The lattice spacing is given by the flow parameter in lattice units (t L 0 = t 0 /a 2 where t 0 = 0.15 fm from the compilation in Ref. [42]).The MS coupling comes from the plaquette, where in the "alpha-V" scheme [49,50] where q * = 3.41/a for the Wilson plaquette gauge action and C F = (N 2 c − 1)/(2N c ) is the usual invariant.Then the conversion to MS is given by and run to α M S (2 GeV) with the usual two-loop formula with In all these equations β 0 and β 1 are the two lowest coefficients of the beta function.Recalling that we note the appearance of the 't Hooft coupling N c α s in all these formulas.(The RI to MS formula, which we do not quote, involves α s C F as well.)We did not find Eq.44 for N c = 3 in the literature but it can be reconstructed from expressions in Ref. [51].
Finally, the MS Z m is run to the final scale µ = 2 GeV; this is done with the usual two loop formula m q (q) = m q (µ) α s (q) α s (µ) This is always a small effect since µa ∼ 1 puts µ in close vicinity to 2 GeV already.The γ i 's are the one and two loop terms for the mass anomalous dimension.
The need to match and run makes the determination of Z m much more fraught that the case of Z A .The issue can be seen in the RI to MS conversion factor.(Here we describe the situation for N c = 3.) Working at aµ ∼ 1, the scale µ ranges from about 1.4 to 2 GeV and α(µ) is in the range 0.18-0.2.The conversion factor can be written as The fraction of R from the highest order term, C 3 α 3 s /R, is about two per cent (slightly bigger at smaller lattice spacing, slightly smaller at bigger lattice spacing, of course).This issue has already been noticed and described by Chetyrkin and Retey [47].The suggested cure, to match at bigger aµ, is not possible for our data sets due to the discretization effects we have already described.We must include a systematic effect in the error budget.We elect to do this by taking C 3 α 3 s /R as a fractional systematic uncertainty on m M S q .(Compare the discussion in Ref. [17], from a study at much smaller lattice spacing.)We again bracket the match between aµ = 0.8 and 1.2.
We conclude with some figures showing the analysis for Z m , in Fig. 3.The errors on Z m (µ) are C 3 α 3 s term.One positive remark about Z M S m worth mentioning: the points shown at different aµ involve slightly different amounts of correction (since the conversion from RI ′ to MS is done at different aµ values with different running distances to µ = 2 GeV).Nevertheless, the results are reasonably independent of aµ even discounting the overall systematic uncertainty.

IV. FITTING STRATEGIES
We believe that little in our fitting strategy is novel.Fits involve minimizing a chi-squared function where the sum runs over all quantities determined by simulation and the entries are X(i) = E(i) − T (i) where E(i) is the measured data value, ∆E(i) is its standard deviation of the mean, T (i) is the model (containing parameters to be fit) and C(i, j) is the correlation between the different experimental values.If the data points are assumed to be uncorrelated, C(i, j) −1 = δ(i, j) in our convention.We include correlations in Eq. 49 by jackknife.Recall that under a jackknife the average a quantity x and its uncertainty ∆x are given in terms of N individual jackknife averages The correlation function is computed similarly, with a convention that the diagonal entries are the identity: Recall all the ingredients we need: am q , am P S , af L , t 0 /a 2 , Z A , Z m .We combine these into dimensionless continuum-regulated quantities, which are then used as input data to fits which determine the LO and NLO chiral parameters.The dimensionless quantities we need, compulsively retaining the lattice spacing, are In practice, only am q , am P S and af L show measurable correlations and we retain only these correlations in the chiral fits.The uncertainties on mq , fP S , and m2 P S come from the uncertainties in their ingredients, combined in quadrature.We also need t 0 /a 2 (actually, its inverse) as a (squared) lattice spacing to add lattice dependent terms in the chiral fits.
In Sec.III B we mentioned the issue that we have to perform a fit determining a set of parameters {a} from y(i) = f (x i , {a}), where both the y i 's and x ′ i s have uncertaintities.We deal with that issue by promoting the x i 's in the fitting function to additional terms in the χ 2 function.For example, in Sec.III B the y i 's were values of t 0 ± ∆t 0 at a set of N bare parameter values, and the x i 's were a set ξ i ± ∆ξ i .The set {a} were a set of J coefficients on a polynomial fit.We expand the χ 2 function to (we suppress the obvious correlation term), so that we now have 2N terms in χ 2 and J + N quantities to be determined.If there were no errors on the x i the counting of degrees of freedom would be N − J, and with the procedure we have outlined it is still 2N − The ξ fits present one final annoyance, since the fitting function must be written entirely in terms of quantities to be fit.This means that in Eq. 18 we must remove the error-bearing m 2 P S from the right hand side of these expressions.We do this by adding a (y(i) − T (i)) = (m P S (i) 2 − M P S (i) 2 ) term to the chi-squared formula.For N d bare parameter values there will be 3N d data points (F P S , m q and M P S ) to be fit, with 4 + 3N d fit parameters, for N d − 4 degrees of freedom.
Finally, we must discuss how we deal with lattice artifacts.Our lattice action has order a 2 artifacts, and this means that the LEC's which will come out of fits will also carry O(a 2 ) corrections.In principle, for x fits, that will be the case for all four LO and NLO parameters and the three additional NNLO ones.To include all these correction terms in the fitting function would be very unwieldy (and, given our uncertainties, quite unstable).Instead, we follow the lead of Ref. [17] and proceed empirically, adding terms and seeing how they affect the χ 2 of a fit.We discover that we are sensitive to lattice artifacts in two places, basically a modification of the right hand sides of Eqs. 10, 17 and 18.To be explicit, we fit the following functional forms.For x fits, Note again that all input variables (m q , M 2 P S , F P S ) have been rescaled by the appropriate power of t 0 .This means that an NLO fit involves six parameters and NNLO, nine.Of these two nuisance parameters, c B is much more important.(It was the only term kept by Ref. [17].)We comment on their values in the next section.
Priors are included in the standard way, as additional terms in the χ 2 function, χ 2 → χ 2 + χ 2 P where We discover that while both x and ξ NLO fits are uniformly stable over a wide range of choices of fit parameters, the NNLO fits are unstable without inclusion of priors.We made a set of fits with a broad range of priors for l 12 , k F , and k M , and we found that values of the four fitted LO and NLO quantities are reasonably independent of the choices we made, but the fitted values and (especially) their uncertainties of l 12 , k F , and k M are completely determined by the priors.This result occurs for all N c 's.We present one illustration, for N c = 5.This is a fit to the data sets with the smallest 15 ξ values.Fig. 4 shows the LO and NLO quantities for ten choices of priors, set identical for l 12 , k F , and k M .The values, from left to right, are (P, ∆P ) = (0.0, 1.0), (0.0,2.0), (0.0,3.0), (1.0,1.0),(1.0,2.0),(1.0,3.0),(2.0,1.0),(2.0,2.0),(2.0,3.0),(3.0,3.0).The fitted values for l 12 , k F , and k M are shown in Fig. 5.In our NLO fits in Sec.V, we will keep the (1.0,2.0)prior.We remark in passing that one lattice artifact which we looked for but did not observe was an extra term in the relation m 2 P S = 2Bm q due to nonzero lattice spacing.This is an extra term on the right hand side of the first equation in the set of Eq. 54.We checked this most thoroughly for NNLO x fits.The addition of terms m 2 P S = 2Bm q (. . . ) + C a (a 2 /t 0 ) or m 2 P S = 2Bm q (. . . ) + C a (a/ √ t 0 ) make a tiny change to the χ 2 (comparing the same fit ranges) compared to leaving them out.All that happens is that the uncertainties on the other fit parameters (mostly F ) grow and the central values drift by a σ or so.A broad prior on C a is needed to stabilize the fit.The authors of Ref. [17] use an action with smeared gauge links which is similar to ours and do not include this term either.Our extraction of LEC's from our data sets is done by performing a series of fits varying the range of quark masses, typically varying the maximum value of ξ in the data set and using model averaging to produce a set of results.The lattice data for f P S and m 2 P S are typically very smooth functions of the quark mass.The fitting functions are also very smooth.It is usually possible to get a good quality fit for any ξ range.The issue is then whether the fit parameters are stable under variation of fit ranges.
The values of SU(N f ) LEC's were already a closed subject before we began our project.They are summarized by FLAG [16] and come from data sets of much higher quality than ours.The study of N c > 3 is our goal.We therefore use N c = 3 as a fiducial: can we reproduce FLAG results?
Before doing any fits, we can just compare our data to NLO curves with these values for the LEC's.Fig. 6 shows two plots: √ t 0 f P S versus √ t 0 m q in panel (a) and √ t 0 m 2 P S /m q versus √ t 0 m q in panel (b).The different plotting symbols correspond to different values of the bare gauge coupling.
The figure shows one loop results with the FLAG parameters for the LEC's.The black curves show the NLO x curve while the dashed red curves the NLO ξ curve.The data, especially for F P S , are much straighter than either NLO fit would favor.In an NLO fit to the data, the fit parameter for l 4 will drift upward as higher mass points are kept, to try to straighten the curve.The data sets we collected mostly lie outside the range where NLO chiral perturbation theory is applicable.To fit all our data sets it is necessary to do NNLO fits.
As a second preliminary picture we display √ t 0 f P S / N c /3 and √ t 0 m 2 P S /m q versus ξ, in Fig. 7.The quantities plotted on the x and y axes of these plots are quite degenerate, but the figure does display the landscape of the data.The different plotting symbols label the different beta values (and hence different lattice spacings) of the data sets.The catalog is: To determine the LEC's, we organized a set of NNLO fits by sorting each N c 's data into a file with increasing ξ and performed a series of fits beginning with a subset of points with the smallest ξ values and extending up to some maximum ξ.We monitored the χ 2 per degree of freedom and confidence level of the fits.Fig. 9 shows results from N c = 3 for B, F , l 3 and l 4 , along with the chi-squared per degree of freedom.The flatness of the fit quantities versus ξ max indicates that we can perform model averaging over our suite of fits to determine the LEC's.
Figs. 10 and 11 show the same information, but for N c = 4 and 5.We show one representative example of our determination of the cutoff dependent terms c B and c F and the NNLO parameters k F , k M and l 12 , for N c = 4, in Fig. 12. Recall that the NNLO fitted parameters are strongly controlled by priors.
Results are listed in listed in Table V and displayed in Fig. 13.The uncertainties on l 3 are much larger than those of l 4 .This is similar to what is seen in the FLAG averages for N c = 3 [16].
As a final check of the consistency of our results, we break apart the NNLO predictions for f P S and m 2 P S /m q into their LO+NLO and NNLO components and plot them.This is shown in Fig. 8.The NNLO piece remains small over the range of √ t 0 m q values used in fits.

B. U (2) chiral perturbation theory
The analysis path for U(2) chiral fits exactly parallels that for SU (2).We performed fits to the NLO and NNLO formulas, studying individual N c data sets.We included additional nuisance parameters to account for finite lattice spacing effects.Specifically, where the NLO and NNLO terms are given in Eqs.24-25.The NLO fits involve six parameters B, F , l 4 , c B and c F ).The NNLO fits add two more: T F and T M , for eight parameters, and replaces the l We need one more input parameter, the quenched topological susceptibility χ T .Ref. [52] measured it to be (See also [53] for earlier determinations of χ T .) We experimented with priors for T F and T M .We discovered that if the fitting range in ξ or √ t 0 m q was large, no priors were needed, while when the range became small, setting a prior 0.2 ± 0.2 for T F and −0.2 ± 0.2 for T M stabilized the fit.Notice that the NLO fits are "fits to a straight line" (plus lattice artifacts) of our data.[16], and fancy diamonds, [17].There are two large N c limits of quenched data, the large volume results of [55] as squares and the small volume ones of [54] as fancy crosses.Purple fancy squares show N f = 4 results from [18].
As in the case of the SU(2), we can generally find fits which have low chi-squares to any data set; the issue is whether the fit parameters are stable under data set size.A second consideration is whether the fit makes sense in that the NNLO contribution is small compare to the NLO one.
N c = 3 is a special case: U(2) chiral fits fail.This can be seen in Fig. 14, where we break up a typical NNLO fit into its component parts.For N c = 3 the actual NNLO component is huge compared to the NLO piece.The fitted B also has a very large uncertainty as it We note that the authors of Ref. [55] extrapolate all their data in 1/N 2 c .Next, we attempt a comparison with the results of Hernandez et al [18].Their N f = 4 data sets are fit to U(4) chiral perturbation theory and they perform interpolations in N f assuming that N f appears in the combination N f /N c (so that tuning N c is a proxy for varying N f ).
They have tables of aF/ √ N c which we convert to √ t 0 3/N c F in our conventions using √ t 0 = 0.15 fm and their a = 0.075 fm.Their N f = 4 results are lower than our N f = 2 ones.
The coefficient of N f /N c in their fitting function is negative and their N f = 2 extrapolations are 0.093(3) fo N c = 3, 0.101 (3) for N c = 4 and 0.105 (3) for N c = 5.Their extrapolation to infinite N c is 0.124(3) -0.130(3) for the two functional forms they use.
We present √ t 0 B = Z s √ t 0 /a(aB) from their tables, converted to the SU(N f ) B, in Fig. 13b.There seems to be little N f dependence on this quantity -which their N f /N c parameterizations also show.We translated their data for L F into l 4 (µ = m 2 π ) by Fig. 13c shows that their N f = 4 results for N c > 3 are consistent with ours.Their prediction of N c = 3, N f = 2 numbers (l 4 = 5.1(3) or 4.1 (11) for two fits) are consistent with ours.Finally, l 3 .We again convert their U(4) quantity to an SU(4) one and plot it in panel d of Fig. 13.This quantity is apparently strongly N f dependent, in addition to being very noisy.They quote an N c = 3, N f = 2 value of l 4 = 0.4 ± 1.6 (at µ = m π ) to be compared with our 3.2 ± 1.5 and FLAG's 3.41(82).
It's a bit dangerous to extrapolate our data to N c → ∞, we think.But (just to save readers the work of doing it themselves) a linear fit to the SU(2) values for √ t 0 B and √ t 0 3/N c F of the form c 1 +c 2 /N c gives large N c values of √ t 0 B = 1.72 (22) and √ t 0 3/N c F = 0.136 (10).
A glance at Fig. 13 shows that these predictions seem to be unexpected.Of course, a convincing plot would benefit from continuum predictions from all the possible approaches (listed at the beginning of the paper) to the large N c limit.Presumably all approaches would converge to the same large N c limit, but would do so in different ways, which might illustrate how varying numbers of fermion flavors affect the LEC's.

VI. CONCLUSIONS
The results shown in Fig. 13 show that we have not answered the list of questions which motivated this project, though we may have made a start: 1. How do the low energy constants of the chiral effective theory scale with N c ?What is their limit?
The expected leading scaling B ∼ N 0 c , F ∼ √ N c , l 3 , l 4 ∼ N c seem to hold.We can see subleading N c behavior in F , l 3 , and l 4 .With three values of N c we cannot say anything about corrections beyond 1/N c .

Is there a crossover to U(2) chiral behavior?
This we cannot tell.U(2) fits fail for N c = 3.We see that in our N c = 4 and 5 data sets both chiral expansions give comparable results for the LEC's, when they are converted to a common scheme.We are aware of no direct calculations of the mass of the flavor singlet meson for N c > 3, so the U(2) fitting functions depend on the validity of the Witten -Veneziano relation.
What could we have done better (and why)?The list is obvious: smaller fermion masses would have been a big help, to be able to do chiral fits deeper in the NLO regime.This then implies a need for bigger lattice volumes, so as not to be compromised by finite size effects.Fortunately, the larger the value of N c , the less this is an issue.Our large lattice spacings led to big uncertainties in the value of lattice to continuum matching factors, and the SU(3) literature makes the obvious point that smaller lattice spacing allows more controlled matching at larger momentum scales.And finally, more values of N c would certainly have been helpful.Trying to extrapolate three values of N c leaves little room to deal with the possibility that quantities scale nonlinearly in 1/N c .We suspect, though, that data sets of the size of the ones reported in FLAG for SU(3) could complete the story of the chiral and continuum limit of large N c QCD.
Of course, to really give high precision answers to the questions we have asked will probably involve some kind of super -analysis involving a variety of approaches including simulations at several fixed N f values, as well as quenched simulations in small and large volumes.But in the meantime, we think that the low-statistics take away continues to be that N c = 3 QCD is not that different from its large N c limit.

FIG. 3 .
FIG. 3. Examples of data sets and results for Z S and Z m .Panels (a) and (c) are from β = 5.35, κ = 0.12775, L = 16 set and (b) and (d) from β = 5.35, κ = 0.1282, L = 24.Panels (a) and (b) show Z S from lattice data versus a µ around µ = 1 while panels (c) and (d) show Z M S m from linear fits to Z S (aµ) centered at the shown aµ values, then converted into Z M S m (µ = 2 GeV).The error bar includes the factor C 3 α 3 s /R, as a systematic.

FIG. 4 .
FIG. 4. Fit values of LO and NLO chiral parameters from an NNLO x fit to N c = 5 over a set of priors for l 12 , k F , and k M .The diamond shows the fit result without priors.

FIG. 5 .
FIG. 5. Fit values of NNLO parameters l 12 , k F , and k M from NNLO x fits to N c = 5 over a set of priors for them.The diamond shows the fit result without priors.
shown in black, squares for β = 5.25, diamonds show β = 5.3, octagons show β = 5.35, crosses show β = 5.4.and √ t 0 f P S vs √ t 0 m q .The different plotting symbols label the different beta values (and hence different lattice spacings) of the data sets as we have listed.The red points are the fit values associated with each data point.The solid lines are the continuum result.The difference between the line and the data points is due to the nuisance parameters C F and C B in Eq. 54.

FIG. 8 .
FIG. 8. NNLO x fits showing √ t 0 f P S and √ t 0 m 2 P S /m q vs √ t 0 m q .(a) and (b) show N c = 3, (c) and (d) show N c = 4, and (e) and (f) show N c = 5 data sets.The fits involve the seven parameters of the NNLO expression, with the NNLO LEC's stabilized by priors as described in the text and two O(a 2 ) correction terms as described in Eq. 54.The different plotting symbols label the different beta values of the data sets and are given in the text.Red points are the fitted values associated with the black data points.The lines show the continuum limit of the fitting functions.The decomposition of the SU (N f ) fitting functions for f P S and m 2 P S /m q (shown as solid lines) is split into their NLO (dashed lines) and NNLO (dash-dotted lines) components.

FIG. 14 . 4 / 3 / 4 / 3 /
FIG. 14. Decomposition of the U (N f ) fitting functions for f P S and m 2 P S /m q (shown as solid lines) into their NLO (dashed lines) and NNLO (dash-dotted lines) components.(a) and (b) N c = 3; (c) and (d) N c = 4; (e) and (f) N c = 5.The data are shown with the same plotting symbols as in Fig. 8.

TABLE III .
Lattice data for N c = 5.The entries C 1 , C 2 and C 3 correspond to the correlation coefficients C mq,m P S , C mq,f P S , and C m P S ,f P S .

TABLE V .
LO and NLO LEC's of SU (2) chiral perturbation theory from our N c = 3, 4, and 5 data sets.

TABLE VI .
LO and NLO LEC's of U (2) chiral perturbation theory from our N c = 4 and 5 data sets.The "converted" columns show conversions to SU (2) quantities according to Eq. 27.In the column we have listed l 3 and l 4 ; note the rescaling.