Meson spectrum in the large $N$ limit

We present the result of our computation of the lowest lying meson masses for SU(N) gauge theory in the large $N$ limit (with $N_f/N\longrightarrow 0$). The final values are given in units of the square root of the string tension, and with errors which account for both statistical and systematic errors. By using 4 different values of the lattice spacing we have seen that our results scale properly. We have studied various values of $N$ (169, 289 and 361) to monitor the N-dependence of the most sensitive quantities. Our methodology is based upon a first principles approach (lattice gauge theory) combined with large $N$ volume independence. We employed both Wilson fermions and twisted mass fermions with maximal twist. In addition to masses in the pseudoscalar, vector, scalar and axial vector channels, we also give results on the pseudoscalar decay constant and various remormalization factors.


Introduction
Large N gauge theories [1] sit at the crux of different approaches to quantum field theory. This is one of the motivations for studying these theories. For the most basic case of a pure Yang-Mills theory we face the difficulties of a strong non-perturbative dynamics without the help of supersymmetry. Nonetheless, the theory exhibits many simplifications with respect to theories at finite N . This comes without paying the price of loosing some of the main non-perturbative phenomena, such as confinement or existence of a mass gap, which remain a challenge for a rigorous proof. From a physicist standpoint it is important to extract the values of the main physical parameters of the theory. These values will remain a testing ground for new ideas and methodologies aiming at an analytic or numerical approach to the theory. In particular, as mentioned earlier, the large N theory seems more easily addressable from studies originating from string theory, as the AdS/CFT approach [2][3][4]. The main first principles approach to non-perturbative dynamics is given by Wilson's lattice gauge theories [5]. Thus, it is extremely desirable to attack the problem using this methodology. The importance of this goal has been realized from long time ago by a small fraction of the community, pioneered by the work of Mike Teper and collaborators [6,7] (See Ref. [8] for a review of results). The intrinsic difficulty of this goal comes from the fact that the large N results follow from extrapolation of those obtained at several small values of N . Fortunately for this program, it seems that many of the physical quantities have a relatively mild dependence on N . On the negative side it turns out that the large N gauge theories possess simplifications that are not present in theories at finite N , making the study of the latter harder. To name one, which is of particular relevance for this work, we have the role paid by quark degrees of freedom in the fundamental representation. In the large N limit these quarks are non-dynamical, so that the so-called quenched approximation becomes exact. On the other hand, at finite N fundamental quarks are dynamical. Generating the configurations becomes much more costly, and if the quenched approximation is used one has to be careful in treating the chiral behaviour of the theory. There is certainly an interplay between the small quark mass and large N limits that one has to be careful about [9,10].
The present work follows a completely different approach. The idea is to exploit one of the simplifications that the large N limit produces. In particular, we will employ the property of volume independence, originating from an observation by Eguchi and Kawai [11]. In essence, the statement says that finite volume corrections go to zero in the large N limit. Although the original formulation was proven incorrect in the continuum limit [12], several ways to enforce the property have been found [12][13][14][15]. Here we will make use of an idea originating soon after the work of Eguchi and Kawai by two of the present authors [16,17]. The main point is to use appropriately chosen twisted boundary conditions [18][19][20]. In this way one can preserve the necessary symmetries to guarantee that the original proof of Eguchi and Kawai holds. As a matter of fact one can take volume independence to the extreme and reduce the lattice volume to single point. This gives a matrix model of d = 4 matrices known as the Twisted Eguchi-Kawai model (TEK). The advantage of this formulation is that because of the reduction in the number of space-time degrees of freedom, one can go to much larger values of N . Thus large N extrapolations become unnecessary and some corrections like those associated to quark loops become negligibly small. There are, however, finite N corrections which turn out to be of a different nature than those of the ordinary theory. These corrections are of two types. One takes the form of ordinary finite volume lattice corrections on a lattice of size ( √ N ) 4 . Hence, a fixed value of N determines a scaling window in which the effective size of the box in physical units remains large enough. Going to smaller lattice spacings one would enter the femtoworld like in standard lattice gauge theories. To enlarge the window one has to increase the value of N . The other type of correction is related to an incomplete cancellation of non-planar contributions. The nature of these corrections are related to effects in non-commutative field theory [21][22][23]. Fortunately, the formulation of TEK has a free integer flux parameter that can be tuned to minimize these corrections and used to quantify their effect [24].
In summary, our formulation of the large N lattice theory (the TEK model) has the capacity of producing estimates of the physical observables of the theory with all statistical and systematic errors under control: measurable and improvable. In the last few years we have been running a number of tests [25] to ensure the capacity of our method to produce physical results within the standard computational resources available at present. The first continuum observable that was studied was the string tension. Indeed, the early estimates [26,27] using our method preceded in many years any other large N lattice estimate. More recently [28], a measurement to a few percent precision has been obtained, which is at least of the level of precision as other recent determinations [29,30] The present work aims at producing a calculation of the low-lying meson spectrum in the large N limit with all statistical and systematic errors under control. The work is the result of several years of study. The methodology was developed in Ref. [31]. Several tests and technical improvements, as well as partial results have been presented in other works [32][33][34]. The meson spectrum has also been investigated by other authors using different techniques. Indeed, this was pioneered by extrapolations from quenched studies [35,36]. There are also studies which employed the idea of volume independence to give determinations [37]. More recent determinations based on the extrapolation by other authors [38][39][40] give results with similar precision to ours. The exact compatibility demands that all estimates have well quantified errors. This is our purpose here. It must be said that the existence of different methodologies is very welcome since there are obviously pros and cons in each of them. We have already mentioned some advantages of our method in relation with the quenched approximation and the chiral limit. However, on the negative side our method does not allow the computation of 1/N corrections, which are also important phenomenologically. Combining our results with finite N estimates could be very interesting.
The structure of the paper is as follows. In the next section we review the philosophy and main formulas involved in our methodology, referring to the original references for explicit derivations. Section 3 is of the most technical nature. Readers mostly interested in the results might skip it. However, for us this section is crucial since it lists and analyzes all the steps involved in giving final numbers and the possible sources of errors they might introduce. A lot of effort has been put into it so that our final mass values have trustworthy statistical and systematic errors. Section 4 contains the presentation of our results starting with those associated with the pseudoscalar channel. In this channel we use both Wilson and maximally twisted mass fermions. This turns out to be crucial to obtain a determination of the pseudoscalar decay constant. The vector, scalar, and axial vector meson masses are computed as well. In Section 5 we use the previous information to present our final table of values of the meson masses in the continuum limit. Values are presented with statistical errors and systematic errors mainly arising from the extrapolation to vanishing lattice spacing. Readers interested in the results should address directly to this section. We also compare our results with other determinations and predictions of the meson masses and the pion decay constant. Possible future improvements are discussed. At the end of the paper we give a long list of tables containing the explicit lattice results of our simulations. This might be useful for other researchers who might want to use our bare results for comparison or display. We have ourselves profited from other authors doing the same.

Meson masses from reduced models
As explained in the introduction, the goal of this paper is to compute the masses of mesons with small number of flavours in the large N limit of d = 4 Yang-Mills theory. Our method is based on reduced models. In particular we will be using the twisted Eguchi-Kawai model, which is a model involving d SU(N) matrices without any space-time label. Despite its conceptual simplicity this matrix model has observables whose expectation value in the large N limit coincides with that of ordinary lattice gauge theory at infinite volume and infinite N . In particular this refers to Wilson loops. This statement can be proven non-perturbatively under certain assumptions by Schwinger-Dyson equations [11,16], perturbatively to all orders [17,41] and also tested up to 5 decimal places by direct evaluation of Wilson loops in the ordinary and reduced model [25].
Let us briefly revise the probability distribution and action density of the TEK model. The former is given in terms of the latter by means of the partition function where the integration on the SU(N) U µ matrices uses the invariant Haar measure on the group. The action S T EK follows by contracting to a point the action of ordinary lattice gauge theories in a finite box with twisted boundary conditions. Many different types of lattice actions can be used, but most of the numerical results have been obtained using the simple Wilson action where 1/b is the lattice equivalent of 't Hooft coupling λ, and the factors z νµ = z * µν = exp{2πin νµ /N } are Nth roots of unity. The integer-valued antisymmetric twist tensor n νµ is irrelevant in the large N limit, provided it is chosen in the appropriate range. A bad choice can yield very large finite N corrections and even a possible breaking of the symmetry conditions for reduction to hold [42][43][44]. In practice our choice has been the socalled symmetric twist which demands N to be the square of an integer N =L 2 . Then one takes n νµ = kL for µ > ν, where k is an integer coprime withL. The choice is irrelevant provided one satisfies the criteria given in Ref. [41].
Reduction implies that the expectation value of a Wilson loop W (C) for an SU(N) gauge theory at infinite volume and infinite N can be obtained as follows where the right-hand side is an expectation value in the TEK model of the trace of the product of link variables following the sequence of directions specified by C (no space-time labels for TEK links). The factor z(C) is a product of the z νµ factor for a collection of plaquettes with boundary on C.
The previous formula (2.3) can be proven non-perturbatively from Schwinger-Dyson equations assuming center symmetry [11,16], shown to be valid in perturbation theory to all orders [17], and checked directly to very high precision numerically [25]. Furthermore, these results give information about how big has N to be to acquire a given precision. Indeed, part of the finite N corrections assume the form of finite size corrections for a lattice of size ( √ N ) 4 =L 4 . Hence, when computing extended observables it is necessary that they all fit inside such an effective box. Typical loop sizes should then be smaller thanL. Lattice masses should also satisfy ML 1. In practice, due to scaling, this affects the range of values of b to be covered. All these limitations are pretty standard in all lattice gauge theory studies. Here we simply have to look atL = √ N as the effective box size. The whole procedure has been tested quite successfully in computing the string tension [28] in the large N limit and other observables [45]. In those studies we went as far as N = 1369, corresponding to a box size of 37 4 , which is quite satisfactory for lattice standards. Unfortunately, in the present study we are unable to reach those values because of the computational effort involved in computing the quark propagator. Our work is mostly based on results obtained for N = 289 corresponding to a box size of 17 4 . This limits the range of our b values, but still seems sufficient to obtain good estimates of the corresponding masses. We have also obtained some results atL = 13 and 19, allowing us to test the effect of finite N .
The methodology for computing meson masses has been explained in previous publications of two of the present authors [31]. We will assume that the number of flavours of quarks in the fundamental stays finite in the large N limit. Hence, quarks are nondynamical and the configurations are generated with the pure gauge TEK model. Quarks do propagate in the background of these gauge fields. One can allow the quark fields to propagate in a lattice of any size (with restrictions) including infinite size. In practice, since the gauge fields live in an effective box of sizeL 4 , we can restrict the quarks to live in a box of the same size, although we will double it in the time direction for ease in computing correlators.
The meson masses are obtained by measuring the exponential decay in time of correlators among quark bilinear operators projected to vanishing spatial momentum. For example, one can consider objects of the form where O is also a matrix in spinor and colour indices. Correlators of these bilinears, after integration over the fermion fields, become expectation values in the probability distribution provided by the TEK action: where D −1 (x, y) stands for the lattice quark propagator between two space-time points. The lattice propagator is the inverse of the Dirac operator. On the lattice there are several options. Most of our results are obtained for the Wilson-Dirac operator, although we will also use the twisted mass [46][47][48] operator. Other possibilities are feasible but have not been used in this work. The set of meson operators used for our work includes both local and non-local ones and will be explained later. All of the operators are gauge-invariant and projected to vanishing spatial momentum. It is important to use various operators with the same quantum numbers, since the corresponding masses should coincide. This helps in obtaining more precise determination of the masses. Another comment is that in Eq. (2.5) we have omitted disconnected pieces, which are subleading at large N . This is an important simplification, which erases the difference between flavour singlet and non-singlet channels. The description done in the previous paragraph looks completely standard. The main difference in our case is that quarks propagate in the background of volume independent gauge fields (up to a twist). This simplifies considerably the structure of the propagators and the meson correlators. An analogy can be drawn with the situation of electrons propagating in a crystal solid. The background field in that case is the potential created by the ions of the solid, which is periodic with a period of the order of the lattice spacing. The motion of the electron in an infinite or much bigger solid translates into a dependence of the propagator on the Bloch momentum. Our case is not exactly identical since the gauge field is only periodic up to a twist. The gauge field only repeats itself exactly after √ N steps in any direction. However, we only consider situations in which the correlation lengths are much smaller than this scale. For that purpose, obviously N should be large enough.
Following the derivation given in Ref. [31] we arrive at the following formula for the meson correlator in momentum spacê where the trace runs over colour and spinorial indices. In a first approximation the operators O and O are just different matrices of the Dirac Clifford algebra defining the quantum numbers of the meson to be studied. Finally, D −1 (q 0 , q) is the inverse of the Dirac operator in our particular setting. We will first of all focus on Wilson fermions, on which most of our results are based. As explained in Ref. [31], after some algebra one can write the Dirac operator as where U µ are the SU(N) matrices generated by the standard Monte Carlo method for the TEK model and Γ µ are fixed SU(N) matrices satisfying where n µν is the same twist tensor used in the TEK action. With our choice of twist tensor the solution of the previous equation is unique up to similarity transformations and multiplication by elements of the center. The first ambiguity is irrelevant since our observables are traces. The phase ambiguity can be reabsorbed into U µ . In practice, we can pick any particular solution and keep it fixed. A final comment affects the range of values of momenta. This is determined by the size of the lattice in which the quark fields live. If quarks propagate in an infinite lattice, then p µ is an angle. However, since the gluon fields are equivalent to those living in a box of size √ N 4 , it is reasonable to confine the fermions to live in a similar box. In that case, momenta range over integer multiples of 2π/(a √ N ). This choice has a bonus, since then the momentum factors can be reabsorbed in U µ , and one can omit the sum over q in Eq. (2.6). For the total temporal momentum we chose p 0 = πn 0 /(a √ N ) with n 0 the integer ranging from 0 to 2 √ N − 1, allowing propagators to extend longer in the time direction. The correlator in configuration space is then given by It is this observable that has been used to extract masses of Wilson fermions.

Data sample and scale setting
The analysis to be presented is based on several years of work accumulating configurations generated with the TEK probability density based on Wilson action as shown in Eqs. (2.1)-(2.2). The main parameters defining each simulation are the value ofL = √ N , the inverse 't Hooft coupling b, and the value of the twist flux integer k. The total number of configurations accumulated for each set of parameters are given on Table 1. For each value a number of thermalization steps was performed initially. We do not appreciate any Monte Carlo time dependence of our results. The configurations used for the analysis were generated using the overrelaxation method explained in Ref. [49], which gives shorter autocorrelation times that the more traditional modified heath-bathá la Fabricius-Hahn [50]. The number of sweeps performed from one configuration to the next N s , also shown in Table 1, is chosen large to ensure that the configurations are largely independent.
The choice of parameters explained in the previous paragraph was dictated by the standard lattice requirement of having small lattice spacings a(b) in physical units. These values can be extracted from our previous analysis of the string tension [28], which used much larger values of N and many more values of b. Having 4 different values of b in our case will allow us to test the scaling behaviour of our data. Scale-setting is an important step in providing physical values in the continuum limit. The idea is to express all dimensionful magnitudes in terms of an appropriate physical unit. Scaling dictates that choosing one unit or other is irrelevant as we approach the critical point (in our case b = ∞). Apart from the string tension one can choose other physical units. The constancy of the ratio of units as we approach the continuum limit is by itself a check of scaling. Many proposals appear in the literature. A good unit must be easily computable and less affected by corrections such as lattice artifacts or finite volume corrections. One possible choice is given by the scaler introduced in Ref. [28]. Essentially, it is similar in spirit to Sommer scale [51] but defined in terms of square Creutz ratios. Our determination, mainly based in our extensive analysis done for our aforementioned string tension paper, is included in the table. A very popular scale lately is the quantity √ 8t 0 defined in terms of the Wilson flow [52]. We have analyzed this quantity for the TEK model and various values of b and N . Using this information we were able to extrapolate the value of √ 8t 0 to N = ∞. The results are given in Table 1. Errors come from a simultaneous fit to various N and are presumably underestimated. Remarkably a simultaneous fit to the ratio √ 8t 0 /r and b ≥ 0.36 gives 1.042(3) with χ 2 smaller than 1. This is a good check of scaling given that both scales involved are based on quite different observables and procedures. Using our data for b ≥ 0.36 we can give determinations √ 8t 0 σ = 1.078(9) andr √ σ = 1.035 (7). The last number can be compared with the estimate 1.051(5) given in Ref. [28], and involving data up to b = 0.385. In conclusion, the different units used are consistent with scaling and translate into possible errors of 2% in the continuum limit.
The standard lattice condition of having sufficiently large physical volumes now translates into large enough values of N . We recall that √ N plays the role of effective size in lattice units. Using the lattice spacing values described in the previous paragraph one can compute the lattice effective linear size of our data l(b, N ) ≡ √ N a(b). This is given in table 1 in string tension units.

Observables
As explained earlier, our main observables are the meson correlation function in channel γ A and γ B at time distance t Table 2: Quantum numbers of the meson channels analyzed in this work and spin-parity assignment.
For the case of Wilson fermions, the Wilson-Dirac matrix D WD (p 0 ) depends on the value of the hopping parameter κ which is a function of the bare quark mass. D WD (p 0 ) acts on color (U µ ), spatial (Γ µ ) and Dirac (γ µ ) spaces and its explicit form is γ A and γ B assign meson quantum numbers in the continuum limit. The choice of elements and their continuum spin parity correspondence is given in Table 2. Eq. (3.1) corresponds to correlators of ultralocal operators. In order to have a range of operators having the same quantum numbers we make use of the smearing method. Smearing can be implemented by replacing γ A in Eq. (3.1) by the operator [38]: Here, is the smearing level andŪ i is the APE-smeared spatial link variable obtained after iterating n ape times the following transformation with c and f free smearing parameters. Here Proj SU(N) means the projection onto the SU(N) matrix. In this paper, we took n ape = 10, c = 0.4 and f = 0.081. The correlation functions of smeared operators are given by The computation of the traces and the inversion of the Dirac operator is performed by means of a stochastic evaluation based on the use of Z 4 random sources [53,54]. Let z(α, β, γ) be the source vector having color (U µ ) index 1 ≤ α ≤ N , spatial (Γ µ ) index 1 ≤ β ≤ N and Dirac (γ µ ) index 1 ≤ γ ≤ 4. The source vectors take values z(α, β, γ) = 1 √ 2 (±1 ± i) and satisfy z * (α , β , γ )z(α, β, γ) = δ α α δ β β δ γ γ , when averaged over random sources. Now, we can replace the trace in Eq. (3.6) by the following expression: where Q(p 0 ) = D(p 0 )γ 5 . To account for the inversion of the Dirac operator, we solve the following 2 √ N equations for y(p 0 ) and the following equation Averaging over random sources, and taking into account that D s is a Hermitian operator, In the actual average, we use 5 random sources. We are using the Conjugate Gradient inversion, so Q −1 actually means For each value of b, many different values of κ have been used (the list will be given later). By fitting the dependence of the PCAC quark mass (defined later) on κ, one can determine κ c (b, N ), which is the chiral limit at which this mass vanishes. The list of values obtained in our simulation are collected in Table 3.
As explained in the previous section, we have also studied twisted mass fermions with maximal twist [46]. This corresponds to a modified Dirac operator of the form where the Wilson-Dirac operator is computed at the value κ c (b, N ) determined earlier with Wilson fermions. After a chiral rotation of the fermion fields this can be written as the Dirac operator of a fermion field of mass µ plus an irrelevant modified Wilson term. The interesting bilinear operators from which we compute correlators involve two different fields having opposite values of the µ-term. It is equivalent to having non-singlet bilinear operators involving two flavours. This modified Dirac operator has many advantages. For example, the constant µ is proportional to the quark mass and monitors the explicit violation of chiral symmetry. It has also important advantages from the point of view of lattice artefacts [47] and most importantly for our purposes allows the determination of the pseudoscalar decay constant without unknown renormalization factors [55] (for a review of properties and results, plus a list of other relevant references we refer the reader to Ref. [48]). As shown in the last section, the original formula for the correlator in momentum space Eq. (2.6) involves a momentum sum. In principle choosing different ranges only produces finite N corrections. A priori it seems that the formula involving the momentum sum is harder to obtain since one would need to invert the propagator for all values of q. However, there is a trick that can be used that costs no extra effort. The idea is for each gauge field configuration to generate a random value of q within the right range and with a uniform distribution. In particular if we consider that the quark field lives in an infinite lattice, for every gauge field configuration we can generate a random angle α µ and replace the link in the Dirac operator as follows: and then apply the same source and inversion procedure as explained earlier. The method produces results at no extra computational cost and we verified that for the free theory the method gives the right propagator. Since the modified links are in U(N) rather than SU(N), we did not impose the unit determinant condition in the Ape smearing procedure. For Wilson fermions we used both methods giving compatible results within errors. For consistency, all the data and masses presented here for the Wilson fermions were obtained with the simple version with no momentum sum. However, our twisted mass results are obtained with the random angle method. Comparison of the results obtained with both methods provides an explicit test that finite N corrections are under control.
In relation with the previous comment, it is clear that stronger finite N corrections are expected when we approach the chiral limit, since the effective lattice size becomes small in units of the inverse pion mass. To avoid this problem, all our results have been obtained for large enough values of m π l. The results are then extrapolated to the chiral limit. We should point out that because of our large values of N this has less problems that for small ones, since chiral logs are strongly suppressed. This is one advantage of using our method.
Finally, as explained earlier, once we obtain the expectation values of the correlators as a function of momenta we Fourier transform them to a time dependent function. For large enough time-separations these correlation functions would decay exponentially in time, the coefficient of which determines the minimal mass in that channel in lattice units. This is in short the procedure. In practice this becomes much more complicated and requires an efficient treatment of all data. This will be explained in the following subsection.

Extraction of masses
As just mentioned, meson masses are extracted from the exponential decay in time of appropriate correlation functions. Although such correlators receive contribution from a tower of states, the lowest mass in the corresponding channel dominates at large times and can be determined from the coefficient of the late single-exponential decay. Nonetheless, this is not as straightforward as it may seem; the determination of the large time decay in an actual simulation faces two problems: the fact that the time extent of the lattice is finite and the deterioration of the signal to noise ratio at large time separation. Therefore, for a precise determination of the ground-state mass it is advantageous to work with operators for which the projection onto the ground state is maximized and the single-exponential decay sets in early on. Smearing is one step in that direction; it increases the overlap with the ground state wave function by extending the operator range at distances comparable with the inverse wavelength of the particle. Beyond that, further improvement can be obtained by using a variational approach [56,57]. The idea is very simple and is based on constructing a linear combination of operators with the same quantum numbers, with coefficients tuned to optimize the ground state projection. We will describe below the particular implementation of the variational procedure used in this work.
To start with, consider a basis of N op operators with the right quantum numbers, in terms of which a N op × N op correlator matrix is built: The variational strategy looks for solutions to the generalized eigenvalue problem (GEVP): for given t 1 and t 0 , with t 1 > t 0 . This is obviously equivalent to solving the standard eigenvalue problem for the matrix: . The correlation matrix given in eq. (3.15) can be rotated to the v (i) basis, leading to the matrix: which is diagonal at t = t 1 with eigenvalues given by λ (i) (t 1 , t 0 ). With a complete basis of operators and for infinite time extent, one would obtain The objective of the procedure is to select the set of operators and the choice of t 0 and t 1 to optimize the projection onto the lowest states, in particular the ground state. The standard GEVP proceeds by repeating the diagonalization at all values of t = t 1 , looking for a plateaux in the effective masses extracted from λ (i) (t 1 , t 0 ). Given the limited time extent of our lattice, we have preferred instead to fix the value of t 1 and work with the correlator matrix defined byC(t, t 1 , t 0 ). To determine the ground state mass in each channel, we use the diagonal matrix element related to the largest eigenvalue λ max (t 1 , t 0 ) = max i λ (i) (t 1 , t 0 ). Let us denote by v max (t 1 , t 0 ) the corresponding eigenvector in terms of which we introduce an optimal, within the given choice of operators and selected values of t 0 and t 1 , correlation function: The ground-state mass is extracted in the usual way from the exponential decay at large times of this correlator. For the specific implementation in this work, we have taken t 0 = a and t 1 = 2a, with a the lattice spacing. The choice obeys to a compromise to maximize the projection onto the ground state without affecting the signal to noise ratio which becomes worse for larger values of t i . We have tested that other choices in the range of values of t i below 4a give compatible results. The operators in the basis are constructed by using different smearing levels of the quark bilinear operator as described in section 3.2. We have considered up    to 10 operators corresponding to fermion smearing levels: 0, 1, 2, 3, 4, 5, 10, 20, 50, 100. Ground-state meson masses in each channel are extracted by fitting the optimal correlator to a hyperbolic-cosine. The optimization procedure involves, first, a choice of the operators in the basis and, second, the choice of fitting range [t min , t max ]. These choices are the main sources of systematic errors in the mass determination. We will discuss the selection criteria in sec. 3.4. Let us just mention here that they respond to two general considerations: the reduction of excited state contamination and the limitation of finite size effects and signal reduction at large times. While smearing improves the projection onto the ground state, at the same time it leads to larger finite size effects and larger statistical errors at large time, for this reason we vary the selection of operators in the basis to limit the maximal amount of smearing. In general we have used the full basis in the pseudoscalar channel and up to 50 smearing levels for the other meson channels.
In order to illustrate the quality of the fits that can be attained with this procedure, we present in figs. 1, 2 one example in each meson channel. The figure corresponds to a lattice with N = 289 at b = 0.355 and κ = 0.16. For each channel, we display the time dependence of the optimal correlator (left) and the effective mass (right). The red band in the plot shows, for illustration, how the determination extracted from the hyperbolic cosine fit compares with the effective mass plateaux, with the length of the band indicating in each case the fitting range [t min , t max ].
In addition to the ground state mass, we have also attempted to estimate the mass of the first excited state in each meson channel. For that purpose, we have performed double exponential fits keeping the ground state mass, m, fixed to the one extracted from the optimal correlator. To be specific, we look at correlators derived including in the variational basis the 4 lowest smearing levels. This correlator is fitted, in the range t ∈ [(2 √ σ) −1 , t max + a], with a combination of two hyperbolic cosine functions with one mass parameter free and the other fixed to the ground-state mass. This allows to determine the first excited state mass m * . The correlator fits for the pseudoscalar and vector mesons shown in fig. 1 include these two contributions with m and m * fixed to the values determined in this way.

Analysis of errors: statistical and systematic
The goal of this paper is to provide values for the mesons masses and decay constants for the large N gauge theory based on first principles that can be used as a test for other estimates based on other types of ideas and approximations. For that purpose the final numbers have to be supplied with trustworthy errors. Some of these errors are statistical. These are indeed the simplest conceptually. In this paper we have employed the well-known jacknife method, based on splitting the sample into several groups of equal size, computing the quantity in question by averaging over all but one of these groups, and equating the error to the dispersion of these values.
A much more difficult task is the estimate of systematic errors. These have different sources and are not necessarily symmetric around the mean value. Most of the errors that we will be concerned with are typical of all lattice calculations, but our methodology enhances or reduces some of those. A very important characteristic of the lattice approach is that all the errors can be estimated.

Simulation and inversion errors
Here we gather together the errors that are characteristic of the numerical procedure. We are fairly confident that the generation of configurations is safe. We have previously carried the same procedure for much larger values of N (up to 1369) and much larger values of b (up to 0.385) and found no relevant systematic effects. Thus, this should even be more so in the more restricted range of parameters that we are exploring here. Furthermore we also played with different thermalization steps and updates per measurement and found no deviations.
Concerning the methodological aspects entering in the evaluation of the observables, we made several tests to quantify the errors committed. In particular we used the same programs and methods for the free field theory (coupling b −→ ∞). The advantage is that the propagators and correlators can be computed analytically in this case, so that the comparison with data gives an idea of the numerical errors involved. In this respect we tested the effect of taking a given number of sources in the calculation of the meson correlators. With 200 sources we found complete agreement between the analytic and numerical computation within errors. With 5 sources some systematic deviations are observed for small values of the quark mass. However, the main effect turns out to be in the meson correlators for p 0 = 0 which produces the addition of a small constant to the correlator. In configuration space the effect only shows up at large correlation times. Thus, in our mass extraction algorithms we have set the upper limit of the fitting range to avoid this region, which is also the most affected by other statistical and systematic errors.
In estimating the final meson masses and decay constants systematic errors are expected to arise from finite lattice spacing and finite size effects. In our particular construction, finite size is traded by finite N . There is an antagonistic interplay among these two error sources. For discretization errors to be small, a(b) should be small, while for finite size errors to be small √ N a(b) should be large. This is pretty tight given that most of our data are obtained forL = √ N = 17. The sensitivity to finite size effects is larger when there are large correlation lengths involved. Thus, in our measurements we have stayed away from the worse region, keeping the pion mass times the effective length large enough. Anyhow, to monitor and estimate the size of these errors we have studied the most sensitive quantities with two other values of N (in addition to our four values of b).
It is important to realize that the finite-size effects in the meson correlators have two sources. On one side, the gauge field configurations generated with the TEK model correspond to a finite volume of sizeL 4 , as Wilson loop expectation values show. On the other hand the quarks can be made to propagate in a box of the same size or larger. These two options for the quark box size can be easily implemented, as explained in the previous section. Thus, the effect can be monitored. We tested this first for the case of free fermion fields (b −→ ∞) in which the meson propagator can be computed analytically (by a momentum sum) as well. Indeed, the effects of using different quark box sizes are sizably smaller for the twisted mass Dirac operator than for the Wilson-Dirac one.
For our dynamical configurations at finite b we have also monitored the effect of the fermion propagation volume. The results obtained with the simplest formula, which involves no explicit sum over spatial momenta, corresponds to a fermion box of sizeL 4 . The results corresponding to an infinite quark box were obtained, as explained before, by generating momenta p randomly with a uniform distribution for each configuration. The procedure does not involve any additional cost in computer time and is seen to give the correct results in our free fermion case. We did not find significant differences within errors between the random momenta results and the others.
Finally, we come to the systematic errors which arise from the procedure to extract the masses from the time-dependent meson correlators. This is by far the largest source of systematic errors. In our procedure masses are obtained from a fit to an exponential decay of the time-correlation function of a properly selected combination of mesonic operators with the appropriate quantum numbers. The results then depend on the selection of the operator, the region of fitting and the functional form chosen for the fit. As discussed in sec. 3.3 , the mesonic operator in each channel is extracted from a variational analysis aimed at improving the projection onto the ground state. The basis for this variational analysis is constructed in terms of smeared operators with smearing level 0, 1, 2, 3, 4, 5, 10, 20, 50 and 100. The results that will be presented correspond to a particular choice that we consider optimal and includes all operators up to smearing level 20, 50 or 100, depending on the state under study and the physical volume of the lattice. The selection of fitting range, [t min , t max ] responds to two general considerations: t min is tuned to reduce contamination from excited states, and t max to limit finite size effects and the poor signal to noise ratio at large times. Concrete choices for these quantities depend not only on the channel considered but also on the amount of smearing of the operators in the basis; including operators with larger smearing level improves the projection onto the ground state and allows to reduce the value of t min , but, at the same time, the signal to noise ratio deteriorates at large times, limiting the value of t max . As a rule of thumb, given a meson state and operator choice, t min is fixed in physical units, around the value where the effective mass plateaux sets in, and t max is set to a fraction of the maximal time separation (a √ N ). In order to estimate the systematic effects in the selection of smearing and fitting ranges, we have analyzed the spectrum obtained if only operators with smearing level below 10 are included in the variational basis; as an example we display in fig. 3 the fits to the pseudoscalar correlator, corresponding to the case presented previously in fig. 1. The masses obtained are in good agreement giving m π / √ σ = 1.563(21) and 1.585 (13). In the general case, the results of this exercise lead to masses that are compatible within statistical errors to the ones derived with the optimal operator.
Finally, concerning the functional form used for the fit, the ground state masses are derived from a single hyperbolic cosine fit performed in the range [t min , t max ], however, in the case of the heaviest states (a 0 , a 1 and b 1 ) we also included a constant contribution that improved the χ 2 of the fits at large times. In most cases the constant turned out to be small and the resulting masses compatible within errors with those obtained with no constant added.

Chiral behaviour and the pion mass
We start with the most important channel: that involving correlators among pseudoscalar meson operators. Because of γ 5 -hermiticity the correlation function is positive definite. Our results show a very neat exponential fall-off of the propagator for large enough separations. Using the optimal operator, as explained in the previous section, the exponential behaviour extends to rather low separations. Fig. 1 shows some characteristic correlators illustrating the typical error size. From these correlators a fit allows us to determine the mass of the lowest energy state, which we will call the pion. We show in Fig. 4   1/(2κ). We see that the dependence is linear. Indeed, this is what we expect from the pseudo-Goldstone boson nature of the pion. The quantity 1/(2κ) is a linear function of the quark mass. The vanishing of the pion mass occurs when the quark mass vanishes and we recover the explicit chiral symmetry, which remains spontaneously broken. In Table 3 we list κ π c , the value of κ at which our linear fit predicts the vanishing of the pion mass. We also give the slope and the square root of the Chi-square per degree of freedom of the linear fit (χ ≡ χ 2 /#dof ). The slope is scaled by a(b) in string tension units. This scaling is the expected one if both the lattice pion mass M π = m π a(b) and the bare quark mass M q ≡ 1 2κ − 1 2κc = Z S m q a(b), are scaled with the lattice spacing as shown.  Table 3: Determination of the critical κ from the linear dependence of the pion mass square (κ π c ) or the PCAC mass (κ c ) on 1/(2κ). We also provide the slopes of the linear fits. In the case of the PCAC mass, the slope determines the ratio of renormalization constants Z P /(Z A Z S ). In the case of the pion mass square we express the slope in string tension units. The quantityχ denotes the square root of the reduced Chi-square (χ ≡ χ 2 /#dof ). An alternative way to determine the onset of chiral symmetry is to look directly at the Ward identity and the PCAC relation. This allows the definition of a quark mass, called PCAC mass M PCAC , whose vanishing signals the point where the explicit breaking vanishes. This involves expectation values of the axial vector and pseudoscalar operator between the vacuum and the pion at rest. Hence, the PCAC mass is affected by operator renormalization and depends on the operator that we use. On the lattice we consider the ultralocal version of the axial vector current and of the pseudoscalar operator. To compute it we follow a similar procedure as in Ref. [38] by computing the correlation functions of the two operators in question with the optimized version of the pseudoscalar operator used for the extraction of the pion mass. Hence, the values of the M PCAC masses are determined by fitting to a constant the ratio This new mass is proportional to the quark mass M q with proportionality constant given by the ratio of renormalization factors Z P /(Z A Z S ). Therefore, from the vanishing of M PCAC one can extract an alternative derivation of the value of the critical κ, that should coincide in the continuum, infinite volume limit, with the one determined from the vanishing of the pion mass. The values of κ c extracted in this way, as well as the slopes and values ofχ are also given in Table 3 -an example of the quality of the fits is shown in Fig. 4. This determination of κ c is more precise and will be taken as our best estimate of the chiral limit.
All the previous results concerned Wilson fermions. However, we have also studied maximally twisted mass fermions [46]. For that one has to start with the Dirac operator for massless Wilson fermions, hence at κ = κ c . Then, as explained in Eq. (3.13), one adds a mass term ±i2κ c µγ 5 , where the opposite sign is used for the two propagators involved in a meson correlator (see Eq. (3.13) for the lattice twisted mass Dirac operator). This, after a chiral rotation, is equivalent to having two flavours of equal mass µ and computing the correlation function of non-singlet meson operators. The unfamiliar reader is invited to consult the literature [46][47][48]55]. By construction, the pseudoscalar meson mass vanishes (approximately) at µ = 0. Furthermore, since µ is proportional to the quark mass, the pion mass square should depend linearly on µ. Our results agree with this linear dependence with intercept equal to zero. The fitted slopes are given in Table 4, again scaled with the lattice spacing in string tension units. As an example, fig. 5 Table 4: Slope, in units of the string tension, of the linear dependence of the pion mass square on the twisted mass parameter µ.
Having determined the mass of the pion we can also look at the masses of the excited states with the same quantum numbers. The results of the GEVP for the first excited states turn out not to be precise enough and we have opted for estimating the mass of the excited states by performing double exponential fits to the correlators, keeping the ground state mass fixed. The results for Wilson and twisted mass fermions are presented in tables 12 and 16. An extrapolation to the chiral limit using a linear dependence in the PCAC mass, at given b and N , works well and gives the intercepts presented in table 5, which have rather large statistical errors. For each value of b we have 4 determinations depending on the methodology (Wilson or twisted mass) and the value of N , which are displayed as function of the lattice spacing in fig. 6. This can give an idea of the systematic errors. Given that they are quite large, we cannot perform a reliable continuum extrapolation of the results, a constant fit to all the data gives for instance m    Table 5: Masses of first excited state in the pseudoscalar channel extrapolated to the chiral limit for Wilson (W) and twisted mass (Tw) fermions.

Pion decay constant
Another very interesting quantity to study is the pion decay constant f π . This is given in terms of the matrix element of the axial current between the vacuum and the pion states.
In contrast with what happens with masses, this quantity is sensitive to normalization and renormalization of the operator involved. The naive extension of the standard QCD definition diverges in the large N limit as √ N . Thus, as previous authors, we define the large N decay constant as follows where A 0 is the temporal component of the axial vector current. The definition is cooked to coincide with the standard one for N = 3. As it stands, the definition is symbolic if we do not specify how the |π; p = 0 state is normalized. The formula assumes the relativistically invariant normalization p|q = (2π) 3 (2p 0 )δ 3 ( p − q).
In practice the matrix element appearing in Eq. (4.2) can be extracted from the correlation function of the temporal component of the axial current with any operator that acts as a pion interpolating field. Dividing out by the square root of the correlation function of the interpolating field with itself, we eliminate the dependence on the choice of interpolating field. The mechanism is true both on the continuum as on the lattice. As a lattice interpolating field we can use the optimal operator selected by the GEVP method, reducing the possible contamination with excited states. In any case, this can be monitored by verifying that the correlator decreases in time as exp(−tm π ). This is essentially the same method used by other authors as for example Ref. [38]. A final word of warning comes from the aforementioned necessity of normalizing and renormalizing the lattice Axial current operator to match with the continuum definition. The normalization is rather standard. With our choice of the Wilson Dirac operator the lattice quark field is given by Ψ L = a 3 /(2κ)Ψ. In addition the Axial current is a composite operator and gets a renormalization factor Z A . This depends on the choice of lattice discretized operator. Here we restrict ourselves to the ultralocal version A L 0 (x) =Ψ L (x)γ 0 γ 5 Ψ L (x). Finally, we can express the lattice spacing and the pion mass in string tension units to determine F π /( √ σZ A ) for each of our datasets. The results are collected in the Table 15. Do our data scale? Examining the results for N=289, one can see that if we compare data points having the same value of the pion mass, the tendency is that F π /( √ σZ A ) tends to decrease with increasing b. For example the data at b = 0.36, κ = 0.1577 has very similar pion mass to that of b = 0.365 and κ = 0.1562, but its value of F π /( √ σZ A ) is 10% higher. However, one cannot deduce from this a sizable violation of scaling, since the renormalization factor Z A introduces a explicit b-dependence of the result which goes in the same direction. Hence, to verify this one would need to know the value of Z A for each b.
Driven by this problem we decided to try a different method for determining the value of F π 1 . This is actually the main reason why we decided to include twisted mass fermions.
One of the advantages of this method is that it allows a determination of F π without a renormalization factor [55]. Furthermore, the procedure is somewhat different and serves as an additional test of the robustness of our results. The value of f π is extracted from the expression Deriving the appropriate reduced lattice formula poses no problem. Again, the matrix element is extracted from the correlation function of the ultralocal pseudoscalar operator with the optimal one. The propagator is now the inverse of the twisted mass one, given in the previous section. This includes the Wilson-Dirac operator set at the value of κ c determined earlier for Wilson fermions. Namely the value at which the extrapolated PCAC mass vanishes. The results then depend on the additional parameter µ, for which we choose 4 values for each data set. As mentioned earlier, the pion mass square obtained in each case vanishes linearly in µ. The corresponding results for F π are given in Table 17.
We still observe that some of the results having similar values of m π have noncompatible values of F π . However, after analyzing the results in detail, we discovered that this phenomenon is due to a sizable dependence on the effective volume. In Fig. 7 we display the values of F π as a function of the pion mass square arranged into sets having similar values of the effective box size a √ N = l √ σ. The data scale well within errors. However, the values of F π tend to decrease sizably for small volumes. Having results for three different values of N has been crucial in discovering the origin of this disagreement.
In order to get maximal information from our data we have tried to parameterize the volume dependence of the results by fitting the data to an exponentially decreasing function of the effective size: F π (m π ) + B(m π ) exp{−C √ σl}. The coefficient C turns out to be 1.36 (11) which explains why our results obtained for relatively small volumes are still sizably dependent on the volume. Curiously we did not find any significant dependence in the determination of the masses. The fitted function F π (m π ) gives our estimate of the pion decay constant at infinite volume and infinite N . In the region covered by our data it can be described well by a second degree polynomial in m 2 π . The same goes for the coefficient B(m π ) of the exponentially decreasing term. Fig. 8 shows all our data for N = 289 and N = 361 together with the solid lines giving the result of the fit. An important quantity is the value of F π in the chiral limit. We have done various fits with different parameterizations and in all cases it seems that a value of F π (0) = 0.22 (1)(2) is a safe estimate. The first error being statistical and the second systematic.
Once we have obtained F π from the twisted mass data we can come back to analyze the results obtained for Wilson fermions. Knowing the value of F π the Wilson data allows a non-perturbative determination of Z A . Unfortunately, this cannot be done point by point since the twisted mass and Wilson data do not correspond to exactly the same values of m 2 π . However, it is highly non-trivial that all the Wilson fermion data can be well described by the same functional dependence on m 2 π up to a multiplicative factor 1/Z A . Remarkably it can, as shown in Fig. 8, where the twisted mass data are displayed together with the Wilson Dirac data after a suitable rescaling.
As a spin-off we have obtained a non-perturbative determination of Z A . The corresponding values are listed in Table 6. This can be compared with the value computed in perturbation theory to two loop order in Ref. [58] giving It is well-known that truncated perturbative results on the lattice give poor estimates when computed in terms of the bare coupling constant λ = 1/b. Much better results follow when using improved couplings λ I of which there are various examples. In Table 6 we use various customary definitions to re-express the perturbative formula Eq. 4.4 in terms of those improved couplings. In general, the perturbative estimates tend to give higher values than our non-perturbative result. Remarkably the same phenomenon takes place for SU (3) as commented in Ref. [38].

The vector channel
We have also analyzed the lowest lying spectrum in the vector channel using in this case Wilson fermions. The analysis, following the strategy presented in sec. 3.3, has allowed a determination of the ρ mass, as well as the mass of the first excited state in this channel, ρ * . The final results for all our lattices are compiled in     Slope m   The value of the ρ mass shows a linear dependence in the PCAC quark mass, as expected at leading order in chiral perturbation theory in the large N limit [59]. This dependence is displayed in fig. 9 for all our lattices. By performing a simultaneous fit to all the points with physical size l(b, n) ≥ 3.4, we obtain an intercept in the chiral limit given by m  (26), slightly worsening the value ofχ from 0.44 to 1.0. We have also performed independent fits to each data set, keeping b and N fixed. The corresponding intercepts and slopes are given in table 7. In fig. 10, the values of the intercepts are displayed as a function of the physical volume, l(b, N ); with the exception of the smallest lattice volume, there is no significant finite volume effect in the results. The scaling with the lattice spacing is very good and we see no clear lattice spacing dependence within the estimated errors.
A common way of displaying the dependence of the ρ mass on the quark mass is to use  √ σ > 2.9, are plotted in this way on fig. 11. The dependence on the pion mass square is linear with a slope given by 0.307(5) and with χ = 0.45 .
We have as well determined the masses of the first excited state in the vector channel, following the same procedure discussed for the pseudoscalar. The results are presented in table 13. The intercepts of the linear extrapolation to the chiral limit are given in table 13 and displayed as a function of the lattice spacing in fig. 6. Assuming that the lattice spacing dependence is within errors, and fitting all the results to a constant gives a mass for the excited state of 5.0(2) in units of the string tension with a value ofχ = 1.2.

Other quantum numbers
The ground-state masses in the a 0 , a 1 and b 1 meson channels, for N = 289 and various lattice spacings, are presented in table 14 in Appendix A.1. They have been obtained using Wilson fermions. As mentioned in section 3.3, we extract the masses from a fit to the optimal correlator, including in this case operators with smearing up to 50 in the basis. This allows to determine the ground state masses with good accuracy, although our results are not precise enough to extract the masses of the excited states in each channel.
As in the case of the ρ meson, the extracted masses depend linearly on the PCAC mass. Figures 12 -14 exhibit this linear dependence for the three different channels. A     √ σ. The corresponding slopes are given in table 9, where we also collect the parameters resulting from fits done at a fixed value of the inverse lattice coupling b. As shown in fig. 15, the results show good scaling towards the continuum limit. The dependence on the lattice spacing is negligible in the case of the a 1 and b 1 states, while the scalar meson a 0 shows a tendency to decrease towards the continuum limit. The bands in the plot are obtained from a constant fit to all the results, which works well in all cases except for the a 0 mass. This issue will be further discussed in section 5, where we will present our final estimates for the continuum extrapolated masses.

Summary of Results and Discussion
In this section we will summarize our results and analyze them in relation with other estimates and expectations of meson masses at large N . We have provided a lot of information in tables which could help other researchers to draw their own conclusions. Our methodology has the advantage of having negligible standard large N corrections. However, other type of large N corrections arise, taking the form of finite volume corrections. This can be understood and avoided by working at large enough effective volumes in physical units. Indeed, after monitoring for this effect we do not see appreciable dependence in the masses. On the contrary, the effect is quite sizable in the pion decay constant f π , and extrapolation is needed to obtain sensible results. Fortunately, in this case we have two independent methods (Wilson and twisted mass fermions) to provide additional support to our analysis. To obtain the physical masses from lattice computations one has to extrapolate to the continuum limit. We have expressed the masses obtained at all lattice spacing values in units of the string tension, so that in the continuum limit they should approach a constant value, which is the final continuum result. All the meson masses, with the exception of the pion mass, depend linearly on the quark mass. On the other hand spontaneous chiral symmetry breaking implies that the pion mass square also depends linearly on the quark mass. For a more physical comparison we preferred to combine the two previous results and express the meson masses at each lattice spacing as follows where X specifies the meson in question. Indeed, this is precisely what we see. However, the value of the masses in the chiral limit m (0) X √ σ at each lattice spacing value, were obtained by extrapolating to vanishing PCAC mass, since that quantity is more precise and less affected by finite effective volume corrections than the pion mass itself. Once that is fixed, one can determine the slope coefficient Y X from the remaining data. For the case of F π our results can be well fitted with the following simple parameterization: In this case a large N extrapolation is also needed as explained in the previous section. The final stage is the extrapolation of the chiral limit meson masses and slopes Y X to the continuum limit. Our lattice results, obtained for four different values of the lattice spacing, have to be extrapolated to vanishing spacing. In general, to the leading approximation, the continuum value is approached linearly or quadratically in the lattice spacing depending on the quantity, the system under study and particular discretization employed. For Wilson fermions we expect a linear approach. However, the majority of our physical observables are consistent with a negligible small slope compared to the statistical errors. This can be seen from Figs. 10, 15. Indeed, fitting the values obtained for masses and slopes for all lattice spacings to a unique value gives values ofχ which are well below 1, signalling a good fit. In accordance with that, we have given as our best estimates of the continuum values the ones obtained by this simultaneous fit. These values are collected in Table 10. The described situation applies for all the slopes Y X as well as the masses of the ρ, a 1 and b 1 mesons in the chiral limit. The only exception is the value of the a 0 mass. In this case the fit is not statistically favourable giving a value ofχ = 1.7. More importantly, the data shows a systematic trend towards a decrease with the lattice spacing. Thus, our best estimate for the mass of a 0 will follow from the procedure explained below.
In any case, it is not possible to exclude a possible linear or quadratic dependence on the lattice spacing even for the cases in which the data is consistent with a constant value. The extrapolation of the data using linear or quadratic dependencies on the lattice spacing also depends on which observable ( √ σ,r or t 0 ) is used to fix the scale, despite the good agreement seen in table 1. For the case of the ρ, a 1 and b 1 meson masses, the different extrapolations do not ameliorate theχ of the constant fit. However, we made use of the span defined by all the different extrapolations to give an estimate of the systematic errors, since they exceed those coming from other sources. This shows the typical problem involved in extrapolations. The boundaries of the range covered by the different continuum limit extrapolations relative to the best statistical average define a positive and a negative  Table 10: Masses of mesons in the chiral and continuum limit and slopes corresponding to Eq. (5.1). The central value is our best statistically significant estimate, corresponding to a constant fit for the ρ, a 1 and b 1 mesons and a linear fit in a(b) for the a 0 mass (see text). The column vector give the maximum positive and negative shifts of our best estimate needed to cover the range of all different continuum limit extrapolations. They can be interpreted as the systematic errors of our continuum limit results (see text).
value, shifting up or down the best value. These shifts are shown in Table 10 as a column vector of values to be added or subtracted to the last significant digits of the best value. The situation for the a 0 mass is quite different. In this case, all the extrapolations based on the different choices of scale-setting predict a decrease with the lattice spacing a, accompanied with a better statistical significance, signalled by values ofχ of order 0.75-0.9. The continuum limit values obtained by fitting a linear dependence in a/r and a/ √ 8t 0 are quite consistent with each other. After converting them back to string tension units usinḡ r √ σ = 1.035 and √ 8t 0 σ = 1.078 we get m (0) a 0 / √ σ = 1.829 and 1.816 respectively. Hence, we take the weighted average as our best estimate shown on Table 10. The systematic errors are depicted in the same fashion ranging from the constant fit to the linear fit using a √ σ.
In the case of the slopes Y X (shown in Table 10) the systematic errors associated to the continuum limit extrapolations lie well within the range of statistical errors. Our best estimates for the pion decay constant and its corresponding slope are given, followed by the statistical and systematic error. The latter also depends on the large N extrapolation presented in the previous section. Now let us compare our results with other predictions and calculations. The first comparison can be made with QCD and its meson spectrum. Setting the string tension to √ σ = 440M eV , the values of the physical I=1 meson masses in string tension units are given by 1.76, 2.82, 2.82 and 2.23 for ρ, a 1 , b 1 and a 0 respectively, which are not too far from the 1.71, 2.99, 3.175 and 1.855 of our large N best estimates for the physical value of m 2 π /σ ∼ 0.1. Hence, it seems that, in what respects to meson masses, large N provides a fairly good approximation to the real world.
Our next comparison is with other lattice determinations of the spectrum. As mentioned in the introduction, results using the quenched approximation extrapolated to large N started some years ago [35,36,60]. Results based on ideas of volume independence similar to ours appeared even earlier [37,61]. Very interestingly, some results have also been obtained with two flavours of dynamical quarks, with values reasonably close [39]. The most complete published lattice work is that of Ref. [38]. It involves a very detailed and thorough work producing results on large N spectroscopy by extrapolation of the re-sults obtained at various N , including some at N = 17. These results are obtained at a unique value of the lattice spacing, so that an analysis of the lattice spacing dependence is not possible. However, since the authors use Wilson fermions and a Wilson action with a value of the coupling which seems to correspond to b=0.36 at large N , we can compare their results with what we obtain at that same coupling. The chiral limit masses for a 1 and b 1 and F π are perfectly consistent within errors. Our result for the ρ in the chiral limit (1.63(8)) is slightly larger than their result (1.54(1)), while for the a 0 it is the other way round (2.25(7) versus 2.40 (3)). In any case, there is a very good qualitative agreement, which is remarkable given the very different methodology used by the two determinations. Incidentally our results on the rho mass and F π do not agree with those of Refs. [37,61], despite they being closer methodologically to ours.
An effort to obtain results in the continuum limit was done by some members of the same collaboration. In particular a detailed study with 4 lattice spacing (like ours) for SU (7) has appeared in proceedings of conferences [62]. A more detailed chart including an estimate of the continuum extrapolation of meson masses at large N was presented in the thesis dissertation of Luca Castagnini [63] 2 . The study covers 4 values of the lattice spacing with the finer lattices similar to ours and one point coarser than our data. The final results involve a double extrapolation to N = ∞ and a = 0, which is achieved by a 4 parameter fit to the data of each observable. Hence, the results have to be taken cum granum salis, as recognized by the author. Nonetheless, the final table is strikingly similar to our best values giving chiral limit masses of 1.687(24), 2.93(11), 2.97(13) for ρ, a 1 and b 1 respectively and F π = 0.197 (20). The slopes in a obtained by the fit for a 1 and b 1 are small, in numerical agreement with our results which are consistent with vanishing slopes within errors. There is certain disagreement in their predicted slope in a for the ρ. Mysteriously we end up having the same estimate for the rho mass in the chiral limit and large N .
The case for the scalar meson a 0 deserves special attention. Both our results as well as those of Refs. [62,63] point towards a sizable positive slope with a √ σ. Our best estimate for the a 0 mass in the chiral limit 1.83 (15) is remarkably close to the value 1.81 (17) given in Ref. [63]. This agreement gives a stronger evidence that our observed a-dependence is not a statistical fluctuation but rather a genuine effect. It seems that the scalar mesons are rather special, also having a larger quark mass dependence Y a 0 than the remaining mesons. Indeed, the same happens in QCD and a long lived discussion has been centered about them (see Ref. [64] for a recent review). Actually, it has been proposed that the study of the behaviour of the masses at large N could help in settling some points [65]. A good deal of the controversy has to do with the isoscalar f 0 (500), which in the large N limit is degenerate with the a 0 . Hence, studying the leading 1/N 2 corrections is presumably crucial. Furthermore there are certain predictions of quite different nature [66,67] that suggest that the large N a 0 mass might coincide with that of the rho meson, a possibility which is not inconsistent with our results.
Altogether, it is fair to say that our continuum results are consistent with the ones by Castagnini and, as emphasized earlier, this is more remarkable given the differences in methodology of our approaches. Even the scale-setting is different, since we use our own measurements of the string tension and two other methods to fix the scale.
Our next goal is to comment on other works and results concerning the meson spectrum at large number of colours obtained with other methods. A new paradigm has arisen from the AdS/CFT correspondence [2][3][4], mapping field theory problems into others involving string theory, supergravity or just gravity. This is, no doubt, an interesting breakthrough introducing a new perspective in describing some field theoretical phenomena and allowing the computation of some observables. The large N limit seems to be crucial in making these new methods feasible. The original ideas concern theories which are very different from QCD, being conformal, supersymmetric, and having fermions and scalars in the adjoint representation. As we move away from this situation the level of rigour in the connection decreases. Nonetheless, theories with quarks in the fundamental [68], with reduced or nosupersymmetry and with running couplings, have been studied. Interesting calculations including meson masses have been performed [69][70][71][72]. We address the reader to the review in Ref. [73] for a very nice account and a more complete list of references. It is worth mentioning that even in theories which are not exactly QCD some observables become rather close to our results. For example, the slope Y ρ obtained in Ref. [70] is indeed consistent with our result within errors. We should also emphasize that more phenomenological methods inspired by holography have been proposed [74][75][76]. As a general rule, some of these papers can use our results to fix some of the parameters of their models, but they have the potential of predicting higher excited states which are more difficult to obtain by lattice methods. Along these lines it also worth mentioning the modified string proposal of Ref. [77] which gives meson spectrum results which in some cases are quite compatible with our results.
A final comment concerns possible future improvement of our work. Given the effort involved, the large uncertainties in taking the continuum limit are somewhat discouraging. To obtain a better control it is sometimes good to go to coarser lattices where the effect is more pronounced. However, for coarser lattices using the simplest parameterization becomes more doubtful and adding more parameters spoils the advantage. From our point of view it is better to go to finer lattices and to reduce the errors. For that the most important limitation is the value of N , whose square root translates into an effective lattice size. The present limit is only computational and enters in the cpu resources needed to compute the propagator. Reaching larger values allows a longer time extent of the correlators and hence longer plateaus. Furthermore, one could reach larger values of b and, hence, smaller values of the lattice spacing without running into finite effective size problems. In relation to this, it should be mentioned that there is no need to take volume reduction to the extreme and simulate the 1 point lattice TEK model. The results could be achieved by running in a small lattice of size L 3 × L 0 . However, we advise researchers trying to follow this road to use appropriately chosen twisted boundary conditions. In particular, using symmetric twist the effective lattice size would become L √ N and good results could be obtained with much smaller values of N at very small volumes.