A lattice study of $\pi\pi$ scattering at large $N_\text{c}$

We present the first lattice study of pion-pion scattering with varying number of colors, $N_\text{c}$. We use lattice simulations with four degenerate quark flavors, $N_\text{f}=4$, and $N_\text{c}=3-6$. We focus on two scattering channels that do not involve vacuum diagrams. These correspond to two irreducible representations of the SU(4) flavor group: the fully symmetric one, $SS$, and the fully antisymmetric one, $AA$. The former is a repulsive channel equivalent to the isospin-2 channel of SU(2). By contrast, the latter is attractive and only exists for $N_\text{f} \geq 4$. A representative state is $\left( |D_s^+ \pi^+\rangle - |D^+ K^+\rangle\right)/\sqrt{2}$. Using L\"uscher's formalism, we extract the near-threshold scattering amplitude and we match our results to Chiral Perturbation Theory (ChPT) at large $N_\text{c}$. For this, we compute the analytical U$(N_\text{f})$ ChPT prediction for two-pion scattering, and use the lattice results to constrain the $N_\text{c}$ scaling of the relevant low-energy couplings.


Contents
1 Introduction The 't Hooft limit of QCD [1] constitutes a significant simplification of the theory preserving all the interesting non-perturbative phenomena, such as asymptotic freedom, confinement and chiral symmetry breaking. In the limit of large number of colors, N c , and a gauge coupling scaled as g 2 ∝ N −1 c , QCD is believed to reduce to a theory of free, infinitely narrow and colorless resonances [1][2][3] (i.e., glueballs and hadrons).
Some controversy has arisen as regards the existence of exotic states, such as tetraquarks, in this limit. The old standard lore concluded that tetraquarks cannot exist at large N c [3]. However this reasoning has been revised recently and a different conclusion has been reached [4,5]. The old argument simply implies that the mixing of any tetraquark with a two-meson state vanishes in this limit, but says nothing on whether tetraquarks actually exist. This is an interesting and timely question, in view of recent experimental discoveries of exotic states [6].
The description of the interactions and decays of these resonances requires non-vanishing N −1 c corrections. A first-principle method that can quantify such corrections is the lattice formulation, where the number of colors can be varied [7,8]. Indeed, various groups have successfully studied the N c scaling of some observables, such as the hadron or glueball spectrum, as well as matrix elements [9][10][11][12][13][14][15][16][17][18][19]-see Ref. [20] for a recent review.
In this paper we present the first study of the N c scaling of the two-pion scattering amplitude using lattice QCD. We have used dynamical simulations at different number of colors in the range N c = 3 − 6 with N f = 4 degenerate quarks. This setup was used in Ref. [18] to study the N c scaling of the famous ∆I = 1/2 rule in the so-called GIM limit, where the charm is degenerate with the up quark. In this limit, penguin contractions vanish and therefore the enhancement of the isospin-0 over the isospin-2 channel depends solely on long-range QCD effects. The completion of that study requires the incorporation of final state interactions of the two pions, for which the study of ππ scattering is a first step.
We study the scattering in the two simplest channels, which do not involve vacuum contractions. They correspond to two irreducible representations (irreps) of SU (4). First, we have a fully symmetric 84-dimensional irrep, which is equivalent to the isospin-2 channel for N f = 2. Second, a 20-dimensional one which is antisymmetric for both quark and antiquark indices, contains states with four open flavors, and only exists for N f ≥ 4. Interestingly, as we will see, the antisymmetric channel turns out to be an attractive one, in contrast with the symmetric one which is repulsive.
We measure two-pion energy levels in lattice simulations from which we can extract infinite-volume scattering properties, such as the scattering length, using the celebrated Lüscher's method [21,22]. The same observables can also be predicted in Chiral Perturbation Theory (ChPT) [23,24]. These predictions can be found in the literature for N c = 3 and any N f [25]. However in the context of large N c , it is necessary to include the η meson in the effective theory, since it becomes degenerate with the pions. The vacuum manifold is therefore not SU(N f ) but U(N f ) [26][27][28]. Somewhat surprisingly we have not found any computation of pion scattering in U(N f ) ChPT in the literature, so we present it here for the first time.
The comparison of lattice results with ChPT predictions can be used to determine the relevant low-energy couplings (LECs) of the chiral Lagrangian, and study their scaling in N c . These results can be compared with those from large N c inspired phenomenological approaches such as Resonant Chiral Theory [29,30]. On the other hand, unitarized Chiral Theory has been used to predict the fate of resonances at large N c from low-energy dynamics [31,32]. We will briefly comment on this approach for the attractive channel. Other approaches based on dispersion relations can also be found in the literature ??.
The paper is organized as follows. In Sec. 2, we present the results for the scattering amplitudes in ChPT with and without the η . In Sec. 3, we review the finite-volume formalism, and in Sec. 4 we present our lattice setup. The results for two-pion scattering are presented in Sec. 5, while the fits to ChPT and a comparison with other approaches is included in Sec. 6. We conclude in Sec. 7. The paper contains two appendices: in App. A we quote the main results for the scattering amplitudes in SU(N f ) and U(N f ) ChPT, while in App. B we discuss how ChPT can be extended to study discretization effects.

Chiral Perturbation Theory Predictions
Chiral Perturbation Theory (ChPT) describes QCD at energies below the QCD scale, Λ QCD , in terms of the lightest non-singlet multiplet of pseudoscalar mesons (an octet for N f = 3). These are the pseudo-Nambu-Goldstone bosons (pNGB) arising from spontaneous chiral symmetry breaking, The theory contains a number of unknown couplings, the low-energy couplings (LECs), that need to be determined from experiment or from first principles, i.e., via matching to lattice QCD. Mesonic observables can be computed in terms of these parameters as a perturbative expansion according to the standard power counting [23,24], where m q is the quark mass, and M π and k 2 are the meson mass and momenta, respectively. A subtle point of ChPT in the 't Hooft limit is the treatment of the flavor-singlet pseudoscalar meson-the η . According to the Witten-Veneziano relation [33,34], the mass of this particle receives a large contribution from the U(1) A anomaly: 1 being χ YM the topological susceptibility of the pure Yang-Mills theory, and F π the pion decay constant. 2 Since the χ YM ∼ N 0 c and F 2 π ∼ N c , this contribution is suppressed as N −1 c . Thus, the large N c pattern of chiral symmetry breaking becomes instead and the η is a pNGB that needs to be included in the effective theory at low energies. The required modification of ChPT, sometimes referred to as U(N f ) ChPT, has been studied before [26][27][28][35][36][37]. It requires a modified power counting scheme that includes the scaling of N −1 c . A consistent one has been shown [27,28] to be: While several quantities have been computed in U(N f ) ChPT, meson-meson scattering has not been studied before in this context to the best of our knowledge. We therefore present here the first computation of scattering amplitudes of non-singlet mesons in U(N f ) ChPT, which will be necessary for this work.
In a general theory with N f active flavors, there are seven different two-meson scattering channels 3 of non-singlet mesons [25]. They correspond to different irreducible representations (irreps) of SU(N f ). In this work, we focus on two of them: • The irrep of maximal dimensionality-84 for N f = 4-that is symmetric in both quarks and antiquarks. It is equivalent to the isospin-2 channel of SU(2). We will refer to it as SS channel. A representative state of this channel is the well-known |π ± π ± .
• The irrep containing states with four distinct quark flavors that is antisymmetric in both quarks and antiquarks. This one only exists for N f ≥ 4, and is 20-dimensional for N f = 4. We denote it as the AA channel. A representative state is (|D + s π + − |D + K + )/ √ 2.
Throughout this work, we will generically refer to any of them as R.

SU(N f ) predictions
At leading order (LO) in the chiral expansion, the pion-pion scattering amplitude depends only on the ratio M π /F π . At this order the result is found to be the same for both channels up to a sign where s is the usual Mandelstam variable normalized by the pion mass. At next-to-leading-order (NLO) and for N f ≥ 4 a total of 13 additional LECs become relevant in the chiral Lagrangian, albeit only some linear combinations appear in the observables of interest. The scattering amplitudes for the two channels are known up to next-to-next-to-leading-order (NNLO) for N f degenerate quark flavors [25]. For completeness we reproduce the s-wave projected NLO result in App. A. From these amplitudes, one π π 2 N f where L R ≡ L R (µ) are the following linear combinations of LECs, defined at the renormalization scale µ, Eqs. (2.7) and (2.8) are valid for any number of colors, but the LECs, F 2 π and M 2 π have an implicit N c scaling [38]: The leading N c scaling of L R is thus the same, but subleading effects differ. We therefore expect: (2.12) Another quantity of interest that may be extracted from the scattering amplitude is the effective range, r 0 . At LO both channels have meaning r 0 is positive (negative) in the SS (AA) channel.

U(N f ) predictions
We use the same notation for the LECs in the two chiral theories, even though they will take different values, as we will see when we discuss the matching between the two. At NLO in the U(N f ) power counting of Eq. (2.5), the different observables are given by the tree-level contribution from the LO Lagrangian, and the NLO corrections from the LECs of order O(N c ). For instance, the SS scattering length is given by 14) and similarly for the AA channel. This prediction is however not sufficient [17], and one needs to go to NNLO in the U(N f ) power counting. At this order, loop diagrams-and thus chiral logarithms-start to contribute, and further care has to be taken to include the η meson. Note that, we use the same notation for the LECs in the two chiral theories, even though they will take different values, as we will see when we discuss the matching between the two. The contribution of the η meson is encoded in a few additional diagrams with respect to the SU(N f ) result. These are summarized in Fig. 1, where solid (dashed) lines represent multiplet (singlet) mesons. The NNLO result is given by where ∆Z 2 represents the additional mass renormalization due to η loops-diagram 1(a)and M loop where M η is given by the Witten-Veneziano formula in Eq. (2.3). K R are combinations of products of LECs and are O(N 2 c ). From the results of Ref. [25], we find that the leading N c dependence of these factors is also equal for both channels, i.e., .

(2.18)
It is also possible to get the prediction for the effective range in the U(N f ) theory as was done for SU(N f ) ChPT. The LO result is the same as in Eq. (2.13) while the NNLO result can be reconstructed from the scattering amplitudes in App. A.

Large N c dependence
In the U(N f ) theory, the N c scaling of the LECs in Eq. (2.12) holds, with a common leading N c dependence for both channels. It remains to be seen whether it is possible to derive a simple N f dependence for the L  The scattering of two pions in these channels at the quark level involves two topologies for the quark Wick contractions, as depicted in Fig. 2. We denote them as disconnected, D, and connected, C. Their contribution to the meson-meson correlator in the two channels is There are two color loops in the D diagram, and just one in the C diagram. Thus, they are O(N 2 c ) and O(N c ), respectively. In terms of correlation functions, the scattering length is roughly given by: where C π ∼ O(N c ) is the single-pion correlator, which contains a single color loop. One can see that factorizable diagrams do not contribute to the scattering length. For instance, the N c -leading diagrams contributing to D in Fig. 2(a) are planar contributions that involve gluons exchanged within each loop, and get exactly cancelled by C 2 π . Concerning the nonfactorizable contributions in D, they necessarily involve gluon exchange between the two loops, for example Fig. 3(a). This diagram is suppressed as 1/N 2 c relative to the leading scaling of D, and is thus O(N 0 c ). Finally, the scaling C π scales like C. All this can be summarized as follows: where a − e are numerical constants, independent of N c and N f . Combining the scaling of the diagrams, with the linear combinations in Eq. (2.19), one can derive where the correlated and anticorrelated terms, L (1) c and L (1) a , are now independent of N f . Comparing to Eqs. (2.9) and (2.10), one can see that Finally, one may also check that the factors of 1/N f in Eqs. The standard form of SU(N f ) ChPT can be recovered from the U(N f ) one in the limit of large η mass, M η M π . The effect of integrating out this particle remains encoded in the SU(N f ) LECs, and results in corrections that scale as 1/N f and 1/N 2 f . The matching between both theories at NNLO gives the following relations:

Finite-volume formalism
To compute two-pion scattering amplitudes from lattice results, we use the standard approach: Lüscher formalism [21,22]. This technique relates the finite-volume spectrum of QCD, obtained from Euclidean correlation functions, to the two-meson scattering amplitude.
We are interested in the simpler formalism for the s-wave interaction of two identical (pseudo)scalar particles. The relation reads: where δ 0 is the s-wave phase shift of the interaction, k is their momentum in the center-ofmass (CM) frame, and Z is the Lüscher generalized zeta function [21]. The energy levels on the lattice are related to k as, Near threshold, i.e., k/M π 1, one can parametrize the phase shift using the effective range expansion (ERE), A related result is the so-called threshold expansion. It is a 1/L expansion of the energy shift of the ground-or threshold-state of the two particles. This was worked out through O(L −5 ) in [21], while the O(L −6 ) correction was added in Refs. [39][40][41]: Here, I, J and K are known numerical coefficients that can be found, e.g., in Ref. [39]. Note that, up to O(L −5 ), there is a one-to-one relation between the scattering length and the energy shift. While the full formalism in Eq.

Lattice setup
The lattice setup of this work is similar to that of Refs. [17,18]. Our ensembles have been generated with the HiRep [42,43] code. We use the Iwasaki gauge action [44], and O(a)improved Wilson fermions in the sea. The value of c sw is chosen as described in Ref. [17]. That is, for N c = 3 we consider the one-loop result boosted by the plaquette, and we keep it constant across other values of N c . In all cases, we use periodic boundary conditions. A summary of the simulation parameters is given in Table 1.
In this work, we are interested in spectral quantities, so the use of Wilson fermions in the sea is appropriate. Another option, along the lines of our previous work [17,18], is to use a mixed-action setup with twisted-mass fermions on the valence. In order to achieve maximal twist, we tune the PCAC valence mass to a value compatible with zero, and we tune the bare twisted mass, aµ 0 , so that the valence pion mass, M v π , matches the sea value, M s π . There are various advantages for the latter choice. First, automatic O(a) improvement is achieved. 5 Second, the pion decay constant can be estimated from the matrix elements of the pseudoscalar current, P , without the need of renormalization constants: which is relevant for the chiral fits of this work. The combination of both lattice regularizations will be useful to quantify discretization effects.
In Table 2, we summarize our results for the single-pion quantities. We introduce a new quantity, the chiral paremter, ξ, defined as which will be useful when matching our results to ChPT. Note that the leading finitevolume effects from Refs. [47,48] have been taken into account. The typical size of these corrections is found to be of ∼ 1%. In addition, the errors in the tuning of the PCAC mass in the mixed-action setup, measured by |m v PCAC /µ 0 |, are at most at the few-percent level. This implies that the effects in our results are well below the statistical precision.

Scale setting
In order to study the lattice artifacts of the quantities we explore here, it is useful to know the lattice spacing in physical units. For this, we follow the procedure outlined in Ref. [17].  Table 1. Summary of ensembles used in this work. L and T indicate the number of points in the spatial and temporal extent, respectively. β is the gauge coupling, c sw is the O(a)-improvement coefficient, and am s is the bare mass of the Dirac operator in the sea sector. Finally, am v and aµ 0 are the parameters for the mixed-action setup-the bare mass and the bare twisted mass, in this same order.
Ensemble (23) Table 2. Results for the pion mass for the sea and the valence sector (M s π and M v π , respectively), the pion decay constant, the PCAC mass in units of the twisted mass, and the chiral parameter, as defined in Eq. (4.2). The "L" subscript indicates that the leading finite-volume effects estimated from ChPT have been taken into account [47,48].
This requires measuring t 0 /a 2 , for which an improved estimator is computed, t imp 0 /a 2 , and M s π on each ensemble. Our scale setting condition is then where M s ref is the reference mass. We use the physical value: obtained from extrapolating existing lattice results [49][50][51] at N f = 2, 3 to N f = 4. Indeed, the uncertainty in Eq. (4.4) will turn out to be the largest in the determination of the lattice spacing.
For the "A" ensembles detailed in Table 1 the scale setting was carried out in Ref. [17], and we do not repeat it here. The result was a = 0.075 (2) fm for N c = 3, and very close for N c > 3. Here we focus on the newer and finer "B" and "C" ensembles. The necessary values of t imp 0 /a 2 are given in Table 3. Using the ChPT expectation for the flow scale [52], we find for each set of ensembles the value of t imp 0 /a 2 at the reference mass, M ref . The final results for the lattice spacing can be obtained by combining the value of (t imp 0 /a 2 )| M ref and Eq. (4.4). From the last column of Table 3, and assuming O(a 2 ) scaling, the expected reduction of discretization effects from the ensembles "A" to "C" is of about 40%. Ensemble 5.534(12) 6.10(15) 0.059(2) 3C20 6.05(14) Table 3. Summary of the quantities entering in the scale setting of the "B" and "C" ensembles, which are the new ones generated for this work. Following the procedure of Ref. [17], we use the clover definition of the energy density of the field strength tensor and apply the improvement of Ref. [53]. In addition, we interpolate to the reference mass, M ref = 420 MeV, and calculate the value of the lattice spacing in physical units using Eq. (4.4).

Extraction of energy levels
In order to determine the energy levels on the lattice, we need to measure the appropriate two-pion correlation functions.
where the two-pion operators for the channels of interest are: Here, the single-meson interpolators are projected to zero momentum, that is, and analogously for K + , D + and D + s , with the corresponding quark flavors. As we have seen in Sec. 2 the two possible Wick contractions are those in Fig. 2 and Eq. (2.19).
In the limit of infinite time extent, the two-pion correlation functions must decay as exponentials in Euclidean time. However, the propagation of particles around the periodic boundary in the time direction gives rise to additional contributions to the correlator. These are the so-called thermal effects, which can be sizeable around t T /2. For the zero-momentum correlators, the leading thermal pollution term is time-independent where A and B are real and positive constants related to the matrix elements of the operator, and E R is the two-pion energy in channel R. A useful observable is the following ratio [54]: which eliminates the constant term, and empirically seems to cancel contributions from excited states. The asymptotic time dependence of this ratio is being t = t − T /2 and ∆E R = E R − 2M π , with R 0 an arbitrary normalization.
We perform correlated fits to Eq. (4.10), and estimate the errors using bootstrapping. In these fits, we fix M π to its value obtained from the one-pion correlator in the corresponding bootstrap sample. In addition, we take into account autocorrelations in the Markov chain by using block lengths that as several times larger than the integrated autocorrelation time. We also consider correlations between time slices using the covariance matrix computed from all initial blocks. The fit range is taken to be [t min , T /2 − 1], where t min is chosen such that the resulting value of ∆E R lies on a plateau. The results for the energies are provided in the second column of Table 4. In addition, two examples of these plateaux are shown in Fig. 4-one corresponding to each channel.  Table 4. Results for the energy shifts ∆E R = E R − 2M π for the two channels, and the corresponding scattering lengths computed using Eq.

Results for ππ scattering
In this section, we use our results for ground state energy levels, together with the 1/L expansion of the ground state-given in Eq. (3.4)-to obtain the s-wave scattering length. Since this is a perturbative expansion valid when a 0 /L 1, we test whether this condition is satisfied a posteriori. We also study discretization effects in the scattering amplitude.
We start by computing the scattering length to O(L −5 ). The results are summarized in the third and fifth columns of Table 4. As can be seen, the SS channel is repulsive, while the AA one is attractive. Qualitatively, the scattering lengths are of the same order for both channels, but with opposite signs. This is expected from ChPT and large N c , as discussed in Sec. 2.
In Fig. 5 we compare the scattering length against the LO prediction of ChPT. We observe good agreement for the SS channel, suggesting that higher order corrections are small in this channel. The situation is different for the AA channel, where deviations from the LO are sizeable. We now study a possible lack of convergence of the threshold expansion.

Convergence of threshold expansion
As explained in Sec. 3, Eq. (3.4) is an expansion in powers of a 0 /L. The truncation of the expression to O(L −5 ) introduces some systematic error, that may not be neglected if a 0 /L ∼ 1. In order to study this effect, a comparison is done between the O(L −5 ) and the O(L −6 ) expansions. Since r 0 enters at O(L −6 ), we will estimate its effect based on its LO ChPT prediction given in Eq. (2.13). Using the value of ∆E for each ensemble, we collate the results for M π a 0 at both orders of the threshold expansion.
In the case of the SS channel, we observe that for all the ensembles, both assumptions produce very similar results, with discrepancies smaller than 1σ. On the left panel of  expansions, respectively, whose width comes from the errors of M π L. The red and green points are the results for the scattering length, which are compatible.  The case of the AA channel is more complicated. For the ensembles with smaller ξ and larger L, we find acceptable convergence, and compatible results for a AA 0 at both orders. However, as ξ is increased we arrive to a point at which convergence fails. An example of this is shown in the right panel of Fig. 6, with the same color code as before. We conclude that some of our ensembles lie in a regime in which the 1/L expansion is no longer valid, and so any results for the scattering lengths may have uncontrolled systematic errors. In view of this caveat, we avoid the expansion and use instead the full Lüscher formalism.

Discretization effects
Discretization uncertainties are always a concern in lattice calculations. Previous work in the isospin-2 channel with twisted-mass fermions has shown that discretization errors are not very significant for this channel [55]. However, hints from the two-baryon sector [56] indicate that this may not always be the same in all scattering observables. Indeed, as we will see, we find significant discretization effects in the AA channel. This is, to our knowledge, the first time that such large effects are seen for meson-meson scattering.
We will assess the magnitude of the cutoff effects at N c = 3, since various ensembles at three different lattice spacings are available. First, we compare the results using the pure Wilson and mixed-action setups, which should coincide in the continuum limit. In Fig. 7 we represent the ratio between the energy shifts determined for both regularizations. We observe that for the SS channel all results are close to unity and thus discretization effects are small. However, this is not the case for the AA channel. We observe large discretization effects, which nevertheless seem to decrease as we approach the continuum limit. Since only the AA energy shifts show relevant cutoff effects, we will focus on this channel henceforth.  To understand the absolute impact of discretization effects in both regularizations we consider the continuum extrapolation of the scattering phase shift, k cot δ AA 0 . This implies additional difficulties, since the scattering amplitude is a function of both ξ and the CM momentum k, which differ for each ensemble. To keep a line of constant physics, we first need to extrapolate/interpolate to a reference point (k 2 ref , ξ ref ) for each lattice spacing and fermion regularization.
The procedure is as follows. We first extrapolate each ensemble to k 2 ref = −0.08M 2 π , using Eq. (3.3) with a prior for the effective range, M 2 π r AA 0 a AA 0 ∈ [−5, −1]. This choice is inspired by the LO prediction from ChPT, Eq. (2.13). The width of the interval is chosen based on some recent results that show that this quantity may have sizeable higher order corrections [57]. An example of this extrapolation for one ensemble is shown in Fig. 8. Note that the errors introduced by the width of the prior, depicted by dotted lines, is negligible compared to the statistical one.  Next, at fixed lattice spacing, we interpolate linearly to ξ ref = 0.14, as shown in Fig. 9. Finally, we use the results at the reference point to do a constrained linear continuum extrapolation using both regularizations, and taking into account correlations. The result is shown in Fig. 10. Note that we expect O(a) improvement in both regularizations, since we are using a non-vanishing value of c sw . Our results are consistent with a regularizationindependent continuum limit and O(a 2 ) scaling, as expected.
Henceforth, we stick to the mixed-action setup and explicitly include discretization effects to analytical expressions. To do so we have studied a generalization of ChPT that considers the effects of a Wilson term as well as a twisted mass [58][59][60]. We have not found in ChPT any hint pointing to larger discretization effects in the AA channel compared to the SS one. We have chosen to parameterize the dominant corrections as . Continuum extrapolation of (k/M π ) cot δ 0 . Data from the pure Wilson (blue) and mixed-action (red) setups is fitted simultaneously with a constrained common continuum limit, depicted by an empty diamond.

Fits to Chiral Perturbation Theory
We now match the lattice results for ∆E R to ChPT predictions in Sec. 2 and App. A. This way we extract the relevant LECs and study their large N c scaling.

Fitting procedure
In order to match our lattice results to ChPT, we need the prediction of ∆E R in this effective theory. The starting point is the result for the scattering amplitudes that can be found in Ref. [25] for the SU(N f ) theory, and in Eqs. (A.10) and (A.11) for the U(N f ) case. Following Ref. [17], we choose a convenient renormalization scale, which is related to 4πF π but has no leading N c dependence, Also, in the U(N f ) theory, we take the splitting between M 2 π and M 2 η to be with a 0 ∼ 6.5, as determined using Eq. (2.3) and the topological susceptibility from Ref. [61].
The results for the scattering amplitudes are then projected to an s-wave following Eq. (A.1), and the corresponding phase shift is computed using Eq. (A.2). Discretization effects are also included for the AA channel as explained in Eq. (5.1). Then, the Lüscher equation [Eq. (3.1)] is numerically solved to obtain the finite-volume energies. 6 The energy shifts are finally determined as a function of ξ, M π L and the LECs.
Since we are working only for small values of the momentum, the LECs that appear in the partial amplitudes multiplied by a power of the momentum-L R and L R in Eqs. (A.4) and (A.5) and K R in Eqs. (A.10) and (A.11)-are set to zero. This is justified, as we do not have enough information to constrain them.
Due to non-negligible correlations between ∆E R and ξ, Ordinary Least Squares is not a good option for fitting. We instead use York Regression [64], defining the χ 2 function as where R is the data vector and V is the corresponding covariance matrix. For a simple fit to a single channel, R is composed by a succession of tuples (f (x i + δ i ) − y i , δ i ) corresponding to all the ensembles i that are fitted. Here, x i = ξ i and y i = ∆E R,i are our lattice results, and f the ChPT prediction of the energy shift, computed as explained above. In this case V is a block-diagonal matrix, due to the statistical independence of the different ensembles. By contrast, when we fit discretization effects in the AA channel, additional correlations must be included. This is due inclusion in R of the value of W from Eq. (5.2), and so, V is no longer block-diagonal. Similarly when we consider a simultaneous fit to both channels, the tuples need to be enlarged, Note that, even though the energy shifts depend on both ξ and M π L, only uncertainties in the former are considered. The reason is that relative errors in the latter quantity are much smaller, so we expect them to have a insignificant effect. We have checked this explicitly for some cases and found this to be true.

Fits at fixed N c
We first work at fixed N c and fit the energy shifts of both channels to SU(4) and U(4) ChPT, following the procedure described above. In all cases, we determine the values of L R in Eqs. (A.4) and (A.5) from the fits. Additionally, we determine K R , as defined in Eqs. (A.10) and (A.11), when we fit to the U(4) theory. Finally, for the AA-channel, we also fit the parameter W that accounts for discretization effects, as introduced in Eq. (5.1).
The results for the fits are shown in Tables 5 and 6 for the SS and AA channels, respectively. In the former case, ChPT at the order considered cannot describe the behaviour of the lattice results for the heaviest masses. We have thus not included ensembles with ξ 0.14 in the fit. In addition, as we are using the W determination from the continuum extrapolation as an additional input when fitting the AA channel results, the number of degrees of freedom (dof) is increased by one unit in this case.
The N c scaling of the fitted values of L SS and L AA is shown in Fig. 11. It is well described by the expected leading and subleading N c dependence. However, only in the U(4) theory the leading term is found to be consistent in the two channels-see Eq. (2.12).
In Fig. 12 we also study the N c scaling of K R . At the order we are working, only the leading N c dependence should be included, which implies a constant value of the ratio shown in Fig. 12. However, the results found for both channels (solid lines) are not consistent, as would be expected from Eq. (2.18). A common limit can be found from a constrained linear fit that includes both channels (dashed line), with different subleading N c corrections in K R .
N c SU(4) ChPT U(4) ChPT L SS · 10 3 χ 2 /dof L SS · 10 3 K SS · 10 3 χ 2 /dof 3 −1.85 (7 (17) Figure 11. Results for the fits to ChPT at fixed N c for the SS (blue squares, Table 5) and AA (red circles, Table 6) channels. The lines are the best fit to Eq. (2.12), and empty points represent the large N c limit. In the case of the U(4) theory, we have imposed a common limit for the fit.  Table 5) and AA (red circles,

Simultaneous chiral and N c fits
Given the results above, we are confident to perform global fits including all values of N c . We have done this both for the SU(4) and the U(4) theories, parametrizing L R with a leading and subleading N c terms. Regarding K R in the U(4) theory, we have opted to keep two different leading-N c parameters, one for each channel, since we do not have enough information to constrain also subleading corrections, First, we have performed fits to single channels for the two theories. The results are shown in Table 7. The best fit results of L (0) are different for both channels when using the SU(4) theory, but they are compatible in the U(4) case, in agreement with the results from the previous section and theoretical expectations. This led us to perform a simultaneous fit of both channels to the U(4) theory with a constrained common value of L (0) . The results are presented in Table 8 and the quality of the fits is shown in Fig. 13.

Comparison to previous literature
From the global fit to U(4) ChPT, we observe that the leading N c dependence of the LEC combination that enters ∆E R in both channels, L (0) , is anomalously small. The dominant terms in the LECs for the two channels of interest are therefore the ones subleading in N c . If we recall Eq. (2.23) and combine the results for the two channels, we can further obtain  Table 7. Results of fits to the energy shift in both SS and AA channels to SU(4) and U(4) ChPT, as indicated in the first column. (9) 35.6/33 Table 8. Results of a global fit to the energy shift in U(4) ChPT. The goodness of the fit is ilustrated in Fig. 13. Figure 13. Results for the simultaneous chiral and N c fit of the energy shifts for both SS (top) and AA (bottom) channels. Lattice results are depicted as points, with a different marker for different N c as explained in the plot legend. Empty points in the SS channel are not fitted. Horizontal lines correspond to the best fit to U(4) ChPT, summarized in Table 8. The values are multiplied by (M π L) 3 ξ −1 to eliminate leading chiral and N c dependence.
the N c and N f scaling of the following linear combinations of LECs: The combination that is in principle leading, seems to be suppressed (both in the leading and subleading terms) with respect to the subleading one.
Using this scaling, we can compare our results to previous literature at N c = 3, which exists for the SU(N f ) theory and N f < 4. Thus, only the SS channel is available for comparison. We begin from the U(4) results of Eq. (6.6), for which the N f scaling is known, and determine the value of L S S for the desired number of flavors. Then, we use the matching of Eq. (2.25) to translate our result to the SU(N f ) case. Finally, we change the renormalization scale to that used in the literature. 7 Our initial scale is µ = 1.40 (12) GeV, which we have determined averaging Eq. (6.1) over all ensembles. The error of this quantity is a systematic one coming from subleading N c effects, while the statistical error is negligible.
We first compare to N f = 3 results from ChPT fits to experiment. Setting the renormalization scale at the mass of the ρ resonance, M ρ = 0.77 MeV, we obtain 2) · 10 −3 from Ref. [24]. (6.9) We also compare to lattice results for N f = 2 [67]. In this case, they are at a scale µ = √ 2F π ≈ 130.2 MeV and also a different quantity is quoted where the first reported error is statistic and the second is systematic. Finally, we can compare with the results obtained in the context of Resonant Chiral Theory (RChT) [29,30,69]. RChT assumes that the LECs at large N c are saturated by the contribution obtained from the integration of the low-lying resonances. Interestingly, the single-resonance approximation predicts the same value for both channels, Ref. [30] ( Table 1, 2 nd column), (6.14) which agrees with our U(4) result at large N c , but fails to reproduce the results at N c = 3, where subleading effects dominate.

A resonance pole in the AA channel?
In the case of the AA channel, there is no direct comparison with the literature, since the channel only exists for N f ≥ 4. However, its attractive nature makes it very interesting. An obvious question is whether there could be a resonance in it that might be interpreted as a tetraquark state. This channel has the correct flavor content, and spin-parity quantum numbers to couple to the recently observed X 0 (2900) exotic state, which has been detected in the D − K + channel at LHCb [70,71]. It is well known that ChPT violates unitarity as the CM energy approaches the mass of the lightest resonance. Some approaches have been used to increase the range of validity of ChPT, namely including the resonances in the effective theory [72,73], or using the inverse amplitude method [31,[74][75][76][77][78][79]. The latter restores unitarity by rewritting the s-wave partial amplitude as where T LO 0 and T NLO 0 are the LO and NLO parts of the amplitude, respectively, as explained in App. A. It has been argued that this procedure provides a tool to constrain LECs from the observed resonances, and inversely, to predict the presence of resonances, such as the σ, from the LECs as determined from phenomenological fits [31,[77][78][79].
We can apply this approach to the AA channel looking for the possible presence of a resonance. In Fig. 14, we represent the s-wave phase shift for this channel in the U(4) theory for ξ = 0.1 and N c = 3, as a function of the CM momentum. We use our best fit result for L AA and K AA , while we vary the LECs combinations L AA and L AA -as defined in Eqs. (A.5) and (A.11)-from zero (solid line) to ±L AA (shaded band). K AA is set to zero. We conclude from this analysis that it is plausible that the curve crosses the positive x-axis from above, which would indicate the presence of a resonance.
If such a resonance exits, an interesting question is whether one could consider it a tetraquark state, as opposed to a meson-meson molecule. Weinberg's criterium [80,81] relates the probability of the state being a pure tetraquark, rather than a molecular state, to its field renormalization, that can be obtained in terms of the scattering length, a 0 , and the effective range, r 0 , of the channel as (6.16) Setting ξ = 0.1, N c = 3 and using the LO results from Eq. (2.13) for the effective range, and the NNLO U(4) prediction for the scattering length given in Eq. (2.17), 9 we get a value M π a AA 0 = 0.51 (4), and Z ≈ 0.8. A value fo Z so close to one indicates that this state, if it exists, would more likely be a tetraquark rather than a meson molecule. In future work, we plan to address the possible existence of this resonance, by studying the phase shift at different values of the CM momenta.

Conclusions
In this work we have studied two-pion scattering in the context of the large N c limit of QCD with N f = 4. We have focused on two scattering channels that correspond to irreducible representations of SU(4). The first channel, the fully symmetric SS, is equivalent to the standard QCD isospin-2. The second channel is antisymmetric in both quarks and antiquarks, and is denoted as AA. It is an attractive channel that includes states with four open flavors and is only present in N f ≥ 4.
We have generated lattice simulations with HiRep [42,43] using O(a) improved Wilson fermions, and the Iwasaki gauge action. Most of the ensembles had already been used in previous work [17,18], and we have simulated a finer lattice spacing to gain further insight on discretization effects. In our setup, we use two different valence Dirac operators: Wilson (unitary setup) and twisted-mass (mixed action). We have then determined the two-pion energy levels in all ensembles and regularizations from a fit to the two-pion correlation functions. The technical details of the lattice setup can be found in Sec. 4.
A possibility to analyze the energy levels is the use of the 1/L expansion of the ground state-the so-called threshold expansion. While this works nicely for the SS channel, we have found that the converge of this expansion is poorer for the attractive AA channel. Therefore, we decided to use the full Lüscher formalism that is non-perturbative in 1/L. Since we are working with mesons near threshold, we only need to include s-wave interactions.
Whereas both regularizations in the valence sector should agree in the continuum limit, that may not be the case at finite lattice spacing. Indeed, we found that even at the level of the energy shifts, the two regularizations show significant cutoff effects in the AA channel. Our continuum extrapolation at N c = 3 indicated that cutoff effects account for a ∼ 20% impact in the phase shift at a = 0.075 fm. This result, along with Ref. [56], is one of the few evidences of significant discretization effects in hadron-hadron scattering. We finally analyze our data matching the two-pion amplitude to Chiral Perturbation Theory. In this case, the correct effective description is U(N f ) ChPT, since the flavor singlet meson becomes light at large N c . Since no results for pion-pion scattering were available in the literature for U(N f ) ChPT, we have computed them in Sec. 2-see also App. A. We find that the leading coefficients in the N c scaling of the LECs that enter these observables are parametrically small, and that the value at N c = 3 is dominated by the subleading coefficient in the 1/N c expansion. Interestingly, within the large N c limit, we also obtain the value of the LECs at N c = 3 and N f = 2, 3, that agrees nicely with previous literature. In addition, we comment on the comparison of our large N c result to the LECs predicted in resonant ChPT: we find good agreement with the leading N c result.
The higher-energy behaviour of the AA channel, and the interplay with the large N c limit is a particularly interesting avenue for future work. At threshold, it is an attractive channel, which could lead to the presence of a resonance at higher values of the centreof-mass momentum. Given the presence of four open flavors in this channel, a resonance would be a candidate for a tetraquark state. Indeed, the recently found X 0 (2900) scalar resonance [70,71] has the correct flavor content to overlap with this channel. In addition, a rough estimate based on the Inverse Amplitude Method in Fig. 14 indicates that it is plausible that such a resonance exists.
A SU(N f ) and U(N f ) ChPT scattering amplitudes In order to match the lattice results to ChPT, we need to be able to relate the twopion infinite-volume scattering amplitude to the corresponding finite-volume energy spectra. This is done using Lüscher's formalism, which reduces to Eq. (3.1) in the case of s-wave scattering. The s-wave phase shift required in this expression can be obtained from the s-wave projected amplitude. In general, the partial wave is defined as where s, t and u are Mandelstam variables normalized by the mass of the pion, P are the Legendre polynomials and θ is the scattering angle in the CM frame. Here, = 0 denotes the s-wave result. The corresponding phase shift is then obtained as In the case of SU(N f ) ChPT, the scattering amplitudes have been computed up to NNLO in Ref. [25]. One can project them to s-wave and compute the phase shift. Here, we quote the s-wave projection up to NLO, using x 2 = M π /F π and q 2 = k 2 /M 2 π as shorthand notation, Here F R (q 2 ) are integrals that need to be numerically evaluated, J(x) is a scalar loop integral that can be found in App. C of Ref. [82], and s = 4(1 + q 2 ) and t(x) = −2q 2 (1 − x). Also new linear combinations of LECs have been defined, The imaginary part of these NLO partial amplitudes, originated from J(x) when x > 4 and required to compute the phase shift, can be related to the LO amplitude by the optical theorem, Regarding the U(N f ) theory, there are no results available in the literature to the best of our knowledge. We have computed the scattering amplitudes in the two channels of interest to NNLO in the counting of Eq. (2.5), and used them to predict the lattice energy spectra, as explained in Sec. 6.1 of the main text. The additional contributions needed to compute the scattering amplitude are explained in Eq. (2.15). We obtain where K R , and K R are the leading-N c part of products of O(N c ) LECs, and B 1 (x) and B 2 (x) are new loop integrals appearing in these theory that correspond to Figs. 1(c) and 1(d), This amplitudes can then be projected to s-wave with Eq. (A.1) and the corresponding phase shift computed using Eq. (A.2). Note that, in spite of using the same name, LECs take different values in the two theories, which can be related in the limit of large η masssee Eqs.( 2.25) and (2.26).

B Wilson ChPT with a twisted mass
In this work, we have employed two methods to reduce discretization effects. First, we included a Symanzik term in the action [83][84][85], and second, we employed a twisted-mass in the valence sector. In the continuum, the latter is just a chiral rotation of the quark fields, but since chiral symmetry is broken in the lattice, this transformation cannot be rotated away and may be used to obtain an O(a) improvement of some observables. Discretization errors induced by the chosen setup in mesonic observables can be studied using a modification of ChPT that includes the effects of a finite lattice spacing and a twisted mass (known as twisted-mass Wilson ChPT). The theory is written in terms of a SU(N f )valued matrix, Σ, and two spurion fields, χ and A, all of which transform in the same way under the chiral symmetry, e.g., Σ → LΣR † , with L ∈ SU(N f ) L and R ∈ SU(N f ) R . Σ is related its usual ChPT counterpart via a chiral rotation, while the spurions include the effects of the twisted mass and the lattice spacing, and are set to constant values, with T the twist matrix and B 0 , W 0 two LECs. The chiral Lagrangian then contains one new operator at O(p 2 ), and 8 more at O(p 4 ), all with their corresponding LEC. Finally, the lattice spacing must be included in the counting scheme. One in general uses although O(a 2 ) effects are often large and need to be promoted to lower orders. This theory has been studied in the literature, mainly for N f = 2 [58][59][60]. In this appendix we review the basis of this extension of ChPT and generalize some results for N f = 4 up to O(a 2 ) using T = diag(+1, −1, −1, 1). Note most of the references work in Euclidean signature, while we opt for Minkoski spacetime with the mostly minus choice of the metric.

B.1 Lagrangian of the theory
The O(p 4 ) Lagrangian of this theory reads In order to work with this theory, one first redefines the spurion field, where F is the bare pion decay constant, φ is the usual pion matrix, and ω = ω 0 + ε is the non-perturbative twist angle. At NLO, we can determine ε by requiring the one-point function to vanish, where we have abbreviated s = sin ω and c = cos ω.

B.2 Discretization effects on physical observables
We are now in position to study discretization effect for some mesonic observables of interest. The lattice corrections to the continuum decay constant, F cont π , and field renormalization, Z cont , are, respectively, F π = F cont π + 4âc F π (4W 4 + W 5 ), (B.9) all of which vanish at maximal twist. The effects on the pion mass are a bit more complicated, since distinct mesons get different corrections to its continuum value, M cont π . For charged mesons one finds which also vanishes at maximal twist. However, neutral multiplet mesons get further O(a 2 ) corrections that do not cancel in such limit and even lead to a non-diagonal mass matrix for flavorless ones, These corrections have been discussed in the literature to be as large as ∼ 30% [60]. Next, one can consider NLO corrections to the continuum scattering amplitudes in the SS and AA channels, M cont SS and M cont AA " respectively. These are generated which are consistent with Ref. [60]. Note that, as initially expected, all these corrections vanish at maximal twist, meaning that the usage of a properly tuned twisted mass provides automatic O(a) improvement for these observables. However, this also implies that in order to understand the observed discretization effects for this regularization, we need to go to higher order. For Wilson fermions, on the other hand, O(a) improvement depends on the correct choice of c sw , and can only be empirically tested. At non-zero twist, further corrections arise from new three-leg vertices that lead to the diagram depicted in Fig. 15 also contributing to the scattering amplitudes for two channels.   where t and u are the usual Mandelstam variables normalized by the pion mass. If we focus on the corrections near threshold and consider only the leading N c terms, we get the result explained in Eq. (5.1). Two comments are in place. First, these leading-N c corrections should have the same size for both the SS and AA channels. Since discretization effects are not observed in the former, one could argue that the leading corrections are suppressed. We have tried several parametrizations when matching to ChPT and concluded that the chosen one in Eq. (5.1) provides the best explanation of the lattice results. Second, through this appendix we have worked in a generalization of SU(N f ) ChPT. The study of a similar generalization of the U(N f ) theory is out of the scope of this paper, and is left for future work.