Finite Size Spectrum of SU(N) Principal Chiral Field from Discrete Hirota Dynamics

Using recently proposed method of discrete Hirota dynamics for integrable (1+1)D quantum field theories on a finite space circle of length L, we derive and test numerically a finite system of nonlinear integral equations for the exact spectrum of energies of SU(N)xSU(N) principal chiral field model as functions of m L, where m is the mass scale. We propose a determinant solution of the underlying Y-system, or Hirota equation, in terms of determinants (Wronskians) of NxN matrices parameterized by N-1 functions of the spectral parameter, with the known analytical properties at finite L. Although the method works in principle for any state, the explicit equations are written for states in the U(1) sector only. For N>2, we encounter and clarify a few subtleties in these equations related to the presence of bound states, absent in the previously considered N=2 case. In particular, we solve these equations numerically for a few low-lying states for N=3 in a wide range of m L.


Introduction
Integrable 1 + 1 dimensional quantum field theories on a finite space circle were rather intensively studied in the last 20 years [1,2,3,4,5,6,7,8,9,10]. A great deal of success in the exact treatment of the finite size effects in various integrable QFT's was due to the thermodynamic Bethe ansatz (TBA) approach resulting in a system (in most of the cases infinite) of non-linear integral equations. It was realized that the TBA equations could be rewritten in a functional, Y-system form [2].
For excited states there are certain poles entering the physical strip, and their qualitative structure is guessed from the ABA. The explicit construction is done only for states in the U (1) sector, but we sketch out the generalization to any state. We show how the exact S-matrix of the model (including the CDD factor) naturally emerges from this approach by simple analyticity assumptions.
The presence of additional singularities on the physical strip related to the bound states, absent for N = 2, leads in N > 2 case to significant modifications, already in the expression for the energies of excited states. We find from our NLIE's the finite size (Lüscher) corrections which reveal the presence of the so called µ-terms. We also test our NLIE's analytically, comparing the results with the known analytic data in the ultraviolet (conformal) limit. Finally, we demonstrate the power of our approach by solving the resulting NLIE's numerically, for the vacuum energy and the energies of some low lying excited states as functions of the size mL for N = 3.
One of the principal motivations for our work is the hope to fulfill the same program in the case of recently constructed AdS/CFT Y-system [14,15,16] for the exact spectrum of anomalous dimensions in N=4 supersymmetric Yang-Mills theory. The PCF model, having N − 1 particles (including N − 2 bound states) in its asymptotic spectrum, bears many similarities with the AdS/CFT case where the number of bound states is infinite. The corresponding Wronskian solutions of AdS/CFT Y-system, or Hirota equation with the so called T-hook boundary conditions is already available [17,18].

The principal chiral field model in the large volume
In this section we will give the definition of the PCF model, remind the reader the scattering theory for the physical particles and the ABA equations, and describe the finite size equations in terms of the Y-system. where the lowest mass scales with the bare charge e 0 according to the asymptotic freedom Ne 2 0 ( Λ is a cut-off). Its wave function transforms in the fundamental representation under each of the SU (N ) subgroups. The exact S-matrix for bi-fundamental particles, found from the conditions of factorizability, crossing, unitarity, analyticity and the bound state structure [19], reads 2 [20,21]: where we introduced the standard SU (N ) R-matrixR L,R (θ) = θ + iP L,R andP is the permutation operator exchanging the left/right spins of the scattering particles. In particular, crossing and unitarity lead to the following identity on the scalar (dressing) factor. We can use this S-matrix to study the spectrum of N particles on a periodic space circle of a sufficiently big circumference L ≫ m −1 imposing the wave function periodicity which quantizes the momenta of the physical particles. The asymptotic spectrum is then given by where θ j are solutions to the system of nested Bethe equations following from the diagonalization of (2.6). This diagonalization can be performed by means of the algebraic Bethe ansatz and leads to the ABA equations (4.20) and (4.9) [22,23]. 3 We will derive these equations as a large L limit of the Y-system of the model -a system of equations valid at any finite volume L and presented in the next subsection. The eq. (4.20) represents the diagonalized version of the periodicity condition (2.6). (4.9) is the set of 2(N − 1) nested Bethe equations for the auxiliary right and left magnon roots u (k) j and v (k) j following from a regularity condition, as we will see in section 4.1. Note that the ABA equations (4.9) remind the Bethe ansatz equations for two inhomogeneous SU (N ) spin chains with the inhomogeneity parameters θ j given by the rapidities of physical particles. Their dynamics is defined by the periodicity equation (4.20). So the large L limit can be also called the "spin chain limit". 2 In the N = 2 case, these definitions give χ CDD = −1, which can also be expressed by multiplying S0 by i. 3 In what follows we will measure all dimensional quantities in the units of the mass m, so that we put everywhere m = 1. The only continuous parameter of the problem is now the volume L. The generalization of these ABA equations to any length L is achieved by the TBA trick [2]: the system is put on the space time torus, with a finite space period L and a big euclidean "time" period R → ∞. Then, using the relativistic invariance, we exchange the roles of time and space and rather solve the problem for the same system but with infinite space extent R placed into a periodic time L interpreted as the inverse temperature [24]. The full energy spectrum of such an infinite system can be found from the nested BAE (2.6) and from (2.7) by means of the so called string hypothesis. The resulting equations for the densities of bound states are presented in [25], following the direct solution of the PCF given in [26] (see also [27]). The free energy calculation at a finite temperature for such an infinite volume system can be done thermodynamically, using the saddle point approximation due to [28]. Then the resulting integral TBA equations can be rearranged into the where, by definition, Y 0,s = Y N,s = ∞ and we have the following boundary conditions at θ → ±∞: This Y-system describing PCF is an infinite set of functional equations (2.8) with the functions Y a,s (θ) of the spectral parameter θ defined in the nodes marked by black and white bullets in the interior of the infinite strip in a, s lattice represented in fig. 2.1.
A direct but rather tedious derivation of this Y-system was performed in the Appendix A of [11] for N = 2. The generalization of this calculation to N > 2 is rather straightforward, but the Y-system (2.8) is known from other considerations [29] and is a very universal system of equations describing the integrable Hirota dynamics. 4 To make many formulas less bulky, the shifts of the spectral parameter will be often denoted as follows As we will see later, the expression for the momentum p a (θ) is the only one compatible with the Y -system and relativistic invariance, up to a normalization that can be absorbed into the definition of the size L of the spin chain 5 . As a result of (2.9) we see that the middle node Y-functions, Y a,0 , a = 1, 2, · · · , N − 1, are exponentially suppressed at large L or at large |θ|.
Obviously, the Y-system (2.8) has many solutions and to specify a solution we have to describe its analytic properties. To have a qualitative idea of the analyticity we have to consider a certain limit for the solution where we know the corresponding Y-functions entirely as analytic functions of the spectral parameter θ. The most convenient limit is L → ∞ where we can solve the Y-system directly, with the appropriate physically natural analyticity assumptions, to obtain explicitly all Y-functions and make a link with the exact scattering matrix and the resulting ABA equations, as it was done for the N = 2 case in [11]. We will give this asymptotic solution in the next section. Then the Y-system, in the form of TBA equations, can be in principle solved numerically by iterations, starting at large L and then adiabatically approaching finite, and even very small L's corresponding to the ultraviolet CFT behavior. The method was successfully used for various integrable sigma models, including the SU (2) PCF [11,30,31,7]. It will be also the main method of this paper devoted to the SU (N ) PCF for N > 2. For the vacuum state, the information from ABA is trivial: it suggests that we don't have any singularities in the physical strip −iN/4 < Im(θ) < iN/4, at least for not too small L's, since there are no Bethe roots. 6 The TBA procedure described above leads to the following expression for the vacuum energy With certain modifications in the analytic properties of Y-functions, described in the next section, the equations (2.8-2.10) appear to be appropriate not only for the vacuum state, as it was originally derived from the string hypothesis, but also for the excited states [3,4]. Y -functions for various excited states differ by their analytic properties which can be qualitatively inferred, as it was mentioned above, from the same states in the ABA. A naive heuristic proposal which worked well for N = 2 case is that the excited states correspond to the appearance of logarithmic poles in the integrand of (2.10) at the points If the contour is deformed so that it encircles these singularities, the pole calculation will give a contribution N j=1 m cosh 2πθ j N which fits well the prediction of the ABA formula (2.7). However, the situation appears to be more complicated at N ≥ 3, already because of the fact that unlike the N = 2 case of [11], the solutions θ j of (2.11) are not necessarily real and this naive prescription should be slightly modified in order to get a real energy. This will be explained in detail in the section 5. One should admit that the right formulas for energies of excited states in the integrable sigma-models are still rather a matter of a natural guess then of a reliable derivation. More insight is needed into this issue.
To solve the Y -system equation (2.8) we will often use it in the form of the Hirota equation

Central node equations
The central node Y-functions Y a,0 related to the black, momentum carrying nodes on the fig. 2.1 play a special role in the Y-system. It will be useful for the future to solve the corresponding Y-system equations for these functions entering the l.h.s. of (2.8) in terms of the r.h.s. Let us rewrite the Y -system (2.8) in the form At s = 0 it can be rewritten using (2.15) as follows: where we introduced a discrete D'Alembert operator ∆ on the interval a ∈ [1, N −1] defined by the formula 7 for any function F a (θ), and where δ is the Kronecker symbol, used to add the counter terms necessary to satisfy (3.1) even at a = 1 and a = N − 1. By a subscript (L) in T in the right hand side of (3.1) are gauge invariant, and we are allowed to write each Y -function in terms of T 's each taken in a different gauge. The meaning and the notation of the gauge (L) will be explained later. We can act by ∆ −1 on both sides of (3.2), to get where p a = cosh( 2θπ N ) sin( aπ N ). The first factor is the only zero mode of ∆ fitting the relativistic invariance and it gives the right large L asymptotics (2.9). Furthermore, the action of the operator ∆ −1 can be easily identified by the fourier transform in θ, a variables, so that the final expression is 8 where the "fusion" operator Π s is defined as the following product Back in the θ-space it takes the form 4. The large L, "spin chain" limit of Y-system and its relation to ABA We will derive in this section the large L, ABA equations (4.20),(4.9) directly from the Y-system (2.8). Following the logic of [11] we use the fact that the Y-functions of the momentum carrying (black) nodes are exponentially small in this limit: This implies that the two wings, left (for s < 0) and right (for s > 0), of the Y-system (2.8) are almost decoupled and can be treated separately.

Expressions for T -functions in the spin chain limit
Eq.(4.1) suggests that either T a,1 ∼ e −Lpa(θ) or T a,−1 ∼ e −Lpa(θ) . Which one does so (whereas another one is finite) is a matter of choice of a gauge for T-functions.
We will work with two different gauges (R) and (L), such that in the large L limit we have T (R) will be called the "right-wing-gauge" and (L) the "left-wing-gauge", and when this superscript will be omitted it will be implicitly assumed that we are working in the (R) gauge.
In the large L limit, the T-functions of the left (L) and right (R) gauge both describe the same Y functions but (up to exponential corrections) they satisfy Hirota equation restricted to the wings s ≥ 0 (resp s ≤ 0). Moreover, these T -functions are in this limit analytic on the whole complex plane, and therefore polynomial.
Such a solution of Hirota equation is well known in applications to the fusion procedure in similar spin chain systems, bosonic [13] or even supersymmetric [32,33] . . , N by means of the following generating functional These functions can be further expressed as follows in terms of some Q (W ) -functions, where the superscript W = R, L indicates the wing that we study (either right or left). These Q functions in the corresponding gauge are polynomials characterizing different solutions of Hirota equation in the large L limit, their roots are the Bethe roots describing various excited states -solutions of the Y -system 9 : 9 We assume here that there exists a gauge such that the large L limit is described polynomial Q functions. Although it needs a better understanding from the point of view of Y-system, it is the case if we start treating the large L limit from the S-matrix by the ABA approach.
In particular, we have from (4.3) T 1,1 should be free of poles, i.e. polynomial. But for each Bethe root w j = u By requiring their cancellation in the sum (4.8), we get a constraint on the position of w j , which we will call the auxiliary Bethe equation: The rest of the T-functions in the right wing can be expressed through the Bazhanov- and they are also automatically polynomial in virtue of (4.9). Among these Q-functions, the polynomial function Q N = ϕ, encoding, as its roots, the rapidities of all physical particles, will be of a particular importance, and the vanishing of T (R) a,−1 or T (L) a,1 due to (4.1) implies the following asymptotics 10 These relations translate all the zeroes θ j of Y 1,0 (θ + i N 4 ) + 1, according to the Bethe equation (2.11), into the zeroes of T -functions.

Asymptotic Bethe ansatz (ABA)
Now we will reproduce from the Y-system in large volume limit L → ∞ the ABA equations (4.20),(4.9) for the spectrum of energies. In this spin chain limit, (3.5) can be employed to compute Y a,0 to the leading order, by using the asymptotic behaviors (4.2).
At this point, it is interesting to notice that the crossing relation implies that up to a zero mode of K N which gives for instance As a consequence, the large L limit of equation (3.5) is where the factor χ CDD was added as another zero mode, necessary to transform the double poles and double zeroes of S 2 into the simple ones [21]. We will also see in section 7.3.2 that this factor arises in our Y-system formalism in a natural way.
In conclusion, we have shown here that the Y system implies the familiar ABA equations [25,21]. In the next sections we will see how these ABA equations for the spectrum of PCF can be generalized to any finite size L.

Expressions for the energy of excited states
There exists no full proof procedure to generalize the formula (2.10) to the excited states. The analytic continuation of [4] with respect to the mass is difficult, if possible at all for a general state in our model. The procedure of [4,8] claims that for the excited states a set of logarithmic poles (different for each state) appears under the integral in (2.10). From (4.11) and the ABA we know that at very large L, is a polynomial encoding all real roots 11 For finite L the roots θ j will be shifted and in general become complex. These exact Bethe roots, as opposed to the approximate 11 We consider here for simplicity only the situation when the Bethe roots θj are real in the asymptotic limit. The case when they occur in complex conjugated pairs should not be very different but at the moment we did not try do to it.
) and manipulations with the contours when N = 3 ones given by (4.20), are defined by the exact Bethe equations T a,0 (θ (a) There is a whole family of such roots when a ∈ [0, N ], because even though the two functions T a,0 (θ) and T a+1,0 (θ + i/2) have the same limit at large L, they do not necessarily have the same roots at finite size. Each of these roots also gives rise to two zeroes and two poles in the Y -functions, namely, as we see from (2.15), 1+Y a,0 (θ (a) Among these families of finite size Bethe roots, we will actually restrict ourselves to the roots θ for odd N . We will argue that only those ones will contribute as poles caught by an integration contour.
Separating the logarithmic poles (where 1 + Y a cancels) in the contour integral (2.10) should give a familiar contribution j cosh 2π N θ j to the energy of a finite L state. This appears to be the right, though not completely well understood and justified, answer for some models, including the PCF at N = 2 [11].
For PCF at N > 2 this procedure encounters another difficulty: the zeros under the logarithm in (2.10) appear to correspond in general to complex Bethe roots θ j . We have to decide what is the right integration contour in (2.10) when this formula is applied to an excited state. We are not aware of any well justified procedure for fixing the contour but we shall try to guess it on the basis of our numerical observations and the symmetry considerations.
In the rest of this section, we consider the formula for excited states of the U (1) sector -the one which corresponds to the wave function |Ψ having the the maximal value of total spin S L = S R = N /2 w.r.t. the SU (N ) R and SU (N ) L symmetries. In this case J . In what follows, we shall distinguish even and odd N 's.

Energy of state in the U (1) sector at odd N 's
It is believed that the energy of an excited state can be obtained from (2.10) by an analytic continuation in the parameter L. This continuation has the effect of appearance of new singularities of the integrand in the physical strip in (2.10) and a certain choice of the integration contour, enclosing some singularities of the integrand [3,4]. How it happens in each particular model or state is usually a rather complicated question. It implies the analysis of positions of these singularities at a finite L but the large L asymptotics often serves as an important guiding principle.
Here we propose a formula for the energies of excited states in the U (1) sector which seems to work well for any odd N . It is based on our numerical and analytic observations, in particular for the N = 3 case. It reads as follows so that we have the straight integration contours parallel to the real axis and shifted by −N i 4 + a i 2 . 12 Let us explain the reason for such a choice of contours. First, let us note that in order to have a real energy from (5.1) we should impose the following property of Y-functions under complex-conjugation: Y a,s (θ) = Y N −a,s (θ). We will restrict ourselves to the gauges where this property is a consequence of the relation For finite L, we will focus on the roots θ j defined 13 by T a+1,0 T a−1,0 , each θ j gives rise to two zeros and poles. In particular, 1 + Y N−1 2 ,0 (θ) has a zero and a pole 14 at respective positions θ j − i/4 and θ j − i/4 because T + N−1 2 ,0 (θ j − i/4) = 0 and T N−1 2 ,0 (θ j − i/4) = 0. In the large L limit these zero and pole almost coincide since θ j is almost real. By complex conjugation, we can also say that 1 + Y N+1 2 ,0 (θ) has a zero and a pole at respective positions θ j + i/4 and θ j + i/4. This structure is illustrated for N = 3 in figure 2. From the Lüscher corrections 15 , we can say that the pole occurs below the zero for 1 + Y 1,0 and vice versa for 1 + Y 2,0 , at least for roots with even momentum numbers 16 . This is important to ensure the right answer if we want the contours to be straight.
In (5.1) we chose the integration contour to pass, for the N −1 2 and N +1 2 'th term in the sum, between those zero and pole. Deforming the contour to the real axis and computing 12 One of the advantages of this straight contour is that it can be easily implemented in numerics. We will see indeed that the Y functions can be most easily computed on exactly these lines. We will also see that the statement holds only for roots with even momentum number, but for odd momentum number, the (slightly modified) contour stays very close to this straight line. 13 In this section, we will denote θj for θ , because the other types of finite size roots don't contribute. 14 In addition to this zero and pole, 1 + Y N −1 2 ,0 (θ) has another zero at θj + 3i/4 and a pole at each root of T N −3 2 ,0 , but this will not have any consequence in the contour argument. 15 In section 8.1, we detail how this is proved in the asymptotic limit. Our numerics suggests that it is still true at finite size, and even in the conformal limit. 16 So that the contour will actually have to be slightly modified for roots having odd momentum number.
This will be done in such a manner that (5.3) will stay true. the contributions of the logarithmic poles enclosed by the contour during that deformation, one gets the following formula In the thermodynamic limit (L ≫ 1), the Bethe roots θ j become real and the second line of (5.3) reduces to the asymptotic result (2.7), whereas the term in the first line appears to be O(e −mL ).

Energy of state in the U (1) sector at even N 's
When N is even, the corresponding contour cannot be chosen as a straight line. We will conjecture here the analogue of (5.3) to be simply where the roots θ j are defined by T N 2 ,0 (θ j ) = 0, so that the second term in (5.4) is real due to the reality of T N 2 ,0 . The corresponding contour is shown in figure 3. We should admit here that this formula for the masses at even N has a status of a natural conjecture. We have not enough of numerical, or analytic evidence to be 100% sure in it. It would be good to verify it at least for the mass gap at N = 4, numerically and by means of the Lüscher corrections at large L.

Hirota form of Y-system and wronskian solution
For the principal chiral field, T a,s is defined for a = 0, 1, . . . , N , while Y a,s is defined for a = 1, 2, . . . , N − 1. We can solve the Hirota finite difference equation (2.12) (and the corresponding Y-system) with the appropriate boundary conditions using its integrability. Any solution of (2.12) is gauge equivalent to a solution where T 0,s (θ) = T 0,0 (θ − s i 2 ) and T N,s (θ) = T N,0 (θ + s i 2 ). The most general solution can thus be expressed [13] as an N × N determinant, in terms of 2N unknown functions q j and q j 17 : At this point, q j is not necessarily the complex-conjugate of q j and the gauge freedom reduces to two independent functions g and g As an example of this determinant solution, the large L, spin chain limit solution corresponding to the U (1) sector state with roots θ i is easily identified, in the (R) (or (L)) gauge, by plugging the following definitions for q j into (6.1) To see that this is the correct parameterization of the U (1) solution, first, we can convince ourselves that T a,−1 = 0, T a,0 = ϕ(θ − i N −2a 4 ), and second, that it reproduces the T 1,s generated by (4.3-4.5) where all Q j (θ)| j<N are set to 1. For the vacuum state P ∞ (θ) = Now we will explain how to generalize this large L solution to any finite L.

Solution of the Y -system for PCF at a finite size
This section describes how to solve the Y -system by reducing it to a finite number of non-linear integral equations (NLIEs), that can be solved in its turn by iterative numerical methods. We will focus on U (1) sector states, although the method is in principle applicable to any excited state (see the discussion in subsection 9).

Definition of the jump densities
We propose here an ansatz for the finite size L solution by adding to the large L polynomial expressions (6.4-6.6) for q's certain finite L terms decreasing for θ → ±∞ and exponentially small for L → ∞ or θ → ∞ : The finite L q j 's take thus the form When j < N and Im(θ) ≤ 0 (7.1) When Im(θ) > 0 (7.6) and the polynomial P has the same degree as P ∞ = lim L→∞ P given by (6.6) 18 . As a consequence of these definitions, we have so that f j is actually the jump between the functions q j and q j on the real axis. Eqs.(7.1-7.4) define q j only below the real axis and q j above the real axis, so that the determinant only allows to compute T a,s inside the strip Im(θ) ∈ [− N −2a . We can already see that these strips are the minimal strips to compute the Y functions on the integration contour of equation (5.1), and we will see that it enables to compute the exact 19 energy of states in the U (1) sector, at any length L.
The jump densities f j (η) are well defined on the real axis where they take only imaginary values, which follows from (7.7) and they are exponentially suppressed at large L or large cosh( 2π N η), as can be inferred from (2.9).

Relation to the analyticity of T functions
From the ABA (4.19) and the finite size equation (3.5), we can see that T a+1,0 T a−1,0 has a proper behavior in this strip, being a meromorphic function regular at infinity. On the other hand, when |Im(θ)| = N 4 , 1 + Y a,0 18 The way we fix this polynomial will be explained in section 7.4, where the finite size Bethe equations are discussed. 19 up to the precision of our numerical procedure solving these NLIE's.
oscillates at Re(θ) → ∞, and it diverges when, e.g., |Im(θ)| ∈ [ N 4 , 3N 4 ]. By that reason we conclude that the analyticity strip of 1 + Y a,0 = T + a,0 T − a,0 T a+1,0 T a−1,0 is {θ, |Im(θ)| < N 4 }. From this we can identify the strips where the asymptotics (4.11) hold for T a,s : These conditions ensure the proper analyticity of 1+Y a,0 = T + a,0 T − a,0 T a+1,0 T a−1,0 , and the boundaries of the analyticity strips of each 1 + Y a,0 are given by the boundaries of the analyticity strips of the corresponding T functions 20 . Now, since we know that the T functions are described by wronskian determinants, these analyticity strips suggest that q j is analytic when Im(θ) < 1/2 (7.12) q j is analytic when Im(θ) > −1/2 (7.13) The fact that the analyticity strip for q j ends up at Im(θ) = 1/2 is reflected for instance in the fact that T N −1,0 is not analytic when Im(θ) > N/4 + 1/2, which explains that Y N −1,0 isn't analytic when Im(θ) > N/4. The equations (7.12,7.13) teach us that the analyticity domain is a bit bigger than what is necessary for (7.1-7.4). It tells us that in the definitions (7.5,7.6), the contour can be shifted up to ±i/2. In other words the functions f j (η) are analytic on the strip Im(η) < 1/2.
It is noteworthy that even with these contour deformations, the determinant expressions (7.1-7.4) describe the function T a,s (θ) inside the strip Im(θ) ∈]− N −2a 4 − s 2 −1, − N −2a 4 + s 2 + 1[, which is narrower than in equations (7.9,7.11). But we will show that the relatively narrow strips given by this ansatz are sufficient to solve the Y -system and compute the energies.

Closed system of NLIEs
The gauge freedom under the transformation (6.2,6.3) can be used to impose F 1 (θ) = F 1 (θ) = 0, which leaves only N − 1 independent densities to compute. N − 1 equations on this densities can be obtained by asking that the state is symmetric, i.e. that Y a,−s = Y a,+s and inserting this requirement into the middle node Y -system equation (3.1) 21 , in the same manner as [11]. 20 The analyticity strips for T0,0 and TN,0 can be chosen as half plane thanks to an appropriate gauge. 21 As a consequence, we are solving the Y -system under the two following constraints : Ya,s ∼ e −Lpa(θ)δ s,0 × consta,s on the one hand, and T (L) a,−1 = Ta,1 one the other hand. This second constraint is specific to symmetric states (which includes the states in the U (1) sector), such that Ya,−s = Ya,s.

Equation (3.5) is then reduced to:
On the other hand, we can write the determinant expression (6.1) for T a,−1 θ − i N −2a These c k,l are the coefficients of the determinant (6.1) defining T a,−1 θ − i N −2a

4
, and finally equations (7.14,7.15) can be recast into This is a closed system of equations on {f j (θ)} θ∈R because all coefficients d a,j , and all T 's can be computed out of f j 's through several convolutions. The solution of the Y -system is therefore achieved by solving this system of N − 1 equations on N − 1 densities. The simple inversion of the linear system (7.15) brings (7.17) into the form This H j is a contraction mapping in some vicinity of f j = 0 when L is sufficiently large since it leads to an exponentially small f j (θ). This implies that in some vicinity of L = ∞, the mapping H j has a fixed point that can be numerically found through repeated iterations of H. The way we solve Y -system is therefore simply the iteration of (7.18) and a good news is that, at least for N = 3, even at very small L, this procedure seams to converge to a fix point of H j , giving a complete solution of (7.18) and thus of the Y -system. 22 The two terms in (7.16)  (which corresponds to a principal value in (7.5,7.6)).

Numerically workable form for the NLIE's
One difficulty of this process is that in order to compute the right hand side of (7.17) it is necessary to compute But we have already seen that the description it terms of q j 's having cuts on the real axis allows to compute T 0,0 (θ) only when Im(θ) < − N 4 + 1 2 . So, for instance, T 0,0 (θ + i N 4 ) cannot be computed when the spectral parameter θ is real. The denominator can nonetheless be computed if in the convolutions we shift appropriately both the argument of the kernel and of the T-functions The same idea, applied to the numerator, would give T 0,0 (θ + N i . But the equality fails because K N has a pole at i N −1 2 . Instead, one can use the following relation which simply uses the fact that for ∀f regular, (Π N [f ]) * K N = f . This is true only up to a zero mode of K N , discussed in the section 7.3.2. Finally, the last factor in the eq.(7.17) can be put into a numerically workable form by rewriting This will help us to transforms the eq.(7.17), once the appropriate zero mode is added, into a really closed system of NLIEs where the right hand side can indeed be computed by knowing the functions f j only on the real axis.

χ CDD factor
It is clear from the derivation of eq.(7.17), as well as of the eqs. (7.20, 7.22) that they are fixed only up to a zero mode of the operator K N . A zero mode Z therefore has to be added to (7.17) Such zero mode can include for instance the factors e −L cosh( 2πθ N ) and χ CDD from equation (4.19).
In the asymptotic limit (L → ∞), the zero modes in the equation (7.14) can be obtained by comparison with (4.19). In the same manner, the zero modes implicitly present in equation (7.22) can be computed in the asymptotic limit, by replacing T a,0 's by their asymptotic values in terms of ϕ, see eq.(4.11), so that we can directly compute the zero mode Z when L → ∞. We notice that, although both (7.17) and (7.22) are true up to a non-trivial zero mode, the zero mode in (7.23) happens to be 23 Z = 1, (7.24) at least in the asymptotic limit. In the numerical solution of the Y -system we therefore assume that Z = 1 holds even at a finite size, i.e. that the analyticity structure of the zero modes is the same at finite L as at L → ∞. We explicitly see that at L → ∞, χ CDD defined in (2.4) is taken into account in (7.23). In addition we can check that at finite size, we obtain Y functions having simple poles only.

Finite size Bethe equations
Bethe equations emerge in this procedure as a regularity requirement on the jump densities f j 's. Let us illustrate it for a general U (1) state in the SU (3) case, and also show why these finite L analogues of Bethe equations are equivalent, at large L, to the ABA Bethe equations on the roots of ϕ.
For such a state, the linear system (7.15) can be written as 23 Actually the right hand side of (7.23) is itself defined up to a factor of e i 2kπ N , because any f ⋆K = e K * log f is defined up to a e 2iπ K corresponding to the choice of the branch of the log. As a consequence, a more exact statement for (7.24) is Z = e i 2kπ N . The choice of k is done to reproduce (4.19), where the phase in (2.3) is chosen in such a manner that one particle at rest (θ1 = 0) is a solution of the Bethe equation (4.20).
For states with zero momentum (like the mass gap and the vacuum), this extra phase can also be obtained by requiring a θ → −θ symmetry, which these states should exhibit.
Inverting the matrix A B −A −B , some singularity could occur at the zeroes of the function AB − AB, i.e. when the determinant is zero. If we want f j 's to be regular, we need the numerator to vanish at the same θ to cancel this pole. This gives the following finite size Bethe equation: One can notice that at suchθ j the two conditions in the r.h.s. are equivalent. The claim thatθ j are a finite size analogue of the Bethe roots is supported by the fact that at large L, the roots of AB−AB are precisely the Bethe roots. Indeed, at large L, B ≃ i and A ≃ P ++ − P , giving AB − AB ≃ i(P −− − P + P ++ − P ) = iϕ. Moreover, the second relation in the r.h.s. of (7.26) reduces then to the reality condition Using the leading-order large L expression of Y a,0 in terms of S, eq.(4.19), we get at large L Using the fact that ϕ(θ i ) = 0 at all Bethe roots θ i , and dividing by the complex conjugate, the large L regularity requirement becomes Using the crossing relation, the left hand side becomes simply S(θ j )e iL sinh( 2π 3 θ j ) , so that the finite size regularity condition stated above is equivalent at large L to the asymptotic Bethe equations (4.20).
As a consequence, the iterative solution of the closed, finite size equations (7.23), should start from the expression (7.1,7.4) where P = P ∞ is given in terms of the asymptotic Bethe roots by (6.6), and then at each iteration, this polynomial is updated in order to incorporate this regularity condition. This procedure for finite-size Bethe equations was described here for odd N and for states having real Bethe roots in the asymptotic limit. The subtlety which arises when N is even, or for the states having, in the asymptotic limit, complex-conjugated complexes of Bethe roots, is that the zeroes of the determinant do not lie on the real axis but approximately on R ± i/2. The above procedure can in principle be applied anyway, but its interpretation is less clear because the regularity condition is imposed at the very boundary of the analyticity strip.

Numerical results
As seen in the figure 4, this method allows to compute the energies of U (1) sector excited states for a wide range of lengths L. We can see that in the IR, L → ∞ limit, the energies of individual states behave mainly according to the number of "particles" forming the state -the number of the Bethe roots θ j : The vacuum energy tends to 0, while the energies of the states θ 0 , θ 1 and θ 2 tend to 1, and the energy of θ 00 tends to 2). In the conformal limit L → 0 we will see that the behavior is defined by the "particle's" mode numbers: The energy goes to 2π L (− N 2 −1 12 + n) where n is the total momentum mode number.
Numerical restrictions As the length L is decreasing, the algorithm looks worse and worse converging, and the densities become more and more peaked around the endpoints of the distribution. By choosing a small enough interpolation step (the densities f j are numerically defined by polynomial interpolation from a finite number of values), it was nevertheless possible to make the algorithm reasonably convergent for the considered states and lengths L, when N = 3. Decreasing further the interpolation step means increasing the computation time, and the necessary amount of memory, which puts a limit to our precision and the minimal length 24 . 24 We are aware of the fact that our precision for the energy levels is not sufficiently good. Our numerical results are still quite preliminary, we will try to improve them in the nearest future by using a better interpolation procedure and using more powerful computers. We definitely believe that our procedure can give a much better precision than these first results in the N = 3 case, and we hope also to tackle the N ≥ 4 case.   Unfortunately, at N ≥ 4 the calculations become heavier and with the size of interpolation steps we can afford our algorithm becomes instable already for L of order ∼ 1 (which means we cannot really check the conformal limit for instance). At the moment we cannot say whether this has a physical meaning (like some symmetry breaking down, or some new type of singularity appearing) or whether it is just a numerical artifact, due to a poor numeric accuracy, or to the choice of the equations. For instance, it could be that the function we iterate over stops to be a contraction mapping but still has a fixed point, and maybe even that by rewriting slightly the functions it could become a contraction again, and extracting its fixed point would be possible by iterations.

Leading order results at large L
The approach of this paper allows to compute the first exponential Lüscher type correction, the so called Lüscher correction, to the energy at large L, as we will show now on a few examples.

Vacuum
The large L behavior of vacuum is given by the condition that Y a,0 | L→∞ Lüscher = (T a,1 ) 2 e −Lpa (8.1) where T a,1 is equal to the binomial coefficient N a . This is obtained from (4.19) where ϕ = 1, and can be plugged directly into (2.10) to get the energy to the leading order. By construction, the resulting energy fits well our numeric results when L is large enough 25 .

Mass Gap at N = 3
When N = 3 it suffices to compute Y 1,0 to get the energy, because Y 2,0 = Y 1,0 . Moreover the previous analysis shows that that enables to compute at large L the leading order value of the integral term in (5.3). The analogous analysis for the second term in (5.3) is a bit more tricky, as it involves the position of the Bethe root. This position can be estimated by computing the densities to the leading order, to deduce the first correction to T 1,0 in order to solve the equation For the mass gap, this root should be at the origin, up to exponential corrections in L. Moreover one can show 26 that T 1, That gives which is in very good agreement the numerical results, as can be seen in fig. 5 Moreover, this expression (8.4) coincides exactly with the so-called µ-term [35,36], which is known to dominate the finite-size corrections in the presence of the bound states.

Conformal limit at L → 0
Let us start from the vacuum. At very small L, the effective coupling constant becomes very small e 2 0 (L) ≃ 2π | log L| and we can linearize the field on the group manifold in the vicinity of g(σ, τ ) = I as g −1 ∂ µ g ≃ ∂ µ A, where A(σ, τ ) is a Hermitian N × N traceless matrix field. The SU (N ) PCF model should become a 2d CFT of N 2 − 1 massless bosons: R(L) is very big, the action will be 25 The precision is over 10% when L ≥ 4.  In the ground state, the Casimir effect will define the limiting energy: 12 . The energies of excited states are where n k = (n k , · · · , n (N 2 −1) k ) are the momentum numbers of particles constituting the state. We see that the small L asymptotics of our plots are well described by this formula: The vacuum, and the states θ 0 , θ 0,0 , . . . have total momentum zero, and their energy satisfies L 2π E(L) ≃ − π(N 2 −1)

12
. On the other hand, a state like θ 1 has momentum number equal to 1 and L 2π E(L) ≃ − π(N 2 −1) 12 + 1, etc. The qualitative behavior of the states θ 0 , θ 00 , etc, at very small L's can be explained by the fact that the quantum fields are dominated by their zero modes. Since the momentum modes are not excited the field Y (σ, τ ) does not depend on σ. The action and the hamiltonian become: where the g(τ ) represents the coordinate of a material point (a top) on the group manifold, andĴ is the corresponding angular momentum on the group. The quantum mechanical spectrum of this system is well known: the quantum states are classified according to the irreducible representations of su(N ) characterized by highest weight with components (m 1 ≥ m 2 ≥, · · · , ≥ m N ) usually represented by a Young tableaux λ with N rows with the lengths m j , j = 1, · · · , N . The operator trĴ 2 is nothing but the second Casimir operator with the well known eigenvalues, so that For instance, for a state with only M real roots in the asymptotic limit (and without self-conjugated complexes of roots), we have m 1 = M , m k≥2 = 0, and L 2π This formula explains well the fact that the corresponding plots on fig.4 converge slowly, as inverse logarithm of L, to −(N 2 − 1)/12.
The perturbative calculation of the mass gap [E θ 0 (L) − E 0 (L)] for L ≪ 1 was done in [37] and was compared with the numerical results following from the TBA approach in [6]. We only recall that for the mass gap, in the logarithmic approximation, which is in the perfect agreement with (8.8), as well as with our numerics. The formula (8.8) also gives a prediction for state θ 00 with two zero-momentumparticles, namely 3 . This result is consistent with our numerics when L is between approximately 10 −2 and 1 (see the fig. 5). For smaller L, although the algorithm seams to converge, the numeric precision is to poor because the integrand of (5.1) is very peaked, and the compatibility with L → 0 predictions cannot be checked very accurately in this regime.
Although our approach is based on adding some terms to the finite size solution, it reproduces correctly these conformal expressions, which shows that this description is not only accurate in some vicinity of L = ∞, but even in the conformal limit where L is very small. This confirms that these terms where added in a sufficiently general manner to describe the relevant solutions of the Y -system.

Discussion
We have presented here, on the example of the SU (N ) principal chiral field model, a potentially powerful and rather general approach to the study of the finite volume spectrum of various integrable 1+1 dimensional sigma-models. The approach continues the ideas of [11] where the method was proposed on the example of the SU (2) PCF, but for N > 2 the method has to be seriously reconsidered due to many new physical features w.r.t. the N = 2 case. In particular, the presence of the bound state particles and the non-reality of the Bethe roots at finite L show a few qualitatively new features within our approach.
For virtually all integrable sigma models, the TBA-like approach of Al.Zamolodchikov can be summarized in a very universal system of functional equations, the Y -system. The Y-system is equivalent to the famous Hirota equation -the Master equation of integrability describing in this case the integrable discrete dynamics with respect to a pair of "representational" variables, a, s and the spectral parameter (rapidity) θ. The boundary conditions for a, s are defined by the symmetry algebra of the model, whether as the analytic structure w.r.t. the θ variable is in general the most complicated issue, really defining the dynamics of the model. It is plausible that in fact even the possible analyticity structures are greatly constrained by Hirota dynamics and by the symmetry algebra. It would be interesting to classify possible types of analyticity stemming from Hirota dynamics and some simple physical arguments (relativistic invariance, crossing, absence of certain singularities, etc.) related to the finite volume sigma models, similarly to the S-matrix bootstrap theory of Al.&A.Zamolodchikov valid only at infinite volume. This could lead to an interesting classification of sigma models themselves and eventually to the discovery of new integrable models.
In this paper, we managed to transform one such relativistic σ-model, the SU (N ) × SU (N ) principal chiral field into a finite set of NLIEs, by solving its finite L Y-system in terms of wronskian determinants of a finite number of Q-functions and parameterizing these Q-functions by N − 1 densities correcting their large L asymptotics to any finite L.
Generalizes the analytic and numerical results of [11] to N ≥ 2, and we could numerically check, at least when N = 2, 3, that this procedures solves the Y -system, and enables to compute energies for a wide range of lengths L, compatible with the UV conformal limit and the IR finite size (Lüscher) corrections. On the way, we conjectured the natural generalization of the energy formula for excited states for the U (1) at finite L, to N ≥ 2. This generalization appears to be unexpectedly non-trivial and looks different for even and odd N . The question of definition of the energy formula for excited states deserves a better understanding and hopefully the eventual derivation.
The analysis was done for U (1) sector states, and it certainly can and should be generalized to any excited state, as was done in [11] for N = 2. To do this, one will have to understand the asymptotic terms, and the structure of zeroes. In particular some extra zeroes would appear, which might affect the way the energy is computed by contour manipulation. Apart from that, the main difference with U (1) sector would be that for non-symmetric states (i.e. when Y a,s = Y a,−s ), it would be necessary to introduce N − 1 densities for the right wing and N − 1 densities for the left wing. One would have to write . But a similar approach based on the wronskian solution of Y-system should a priori still enable us to compute the energies of these states.
An interesting problem which our approach might help to solve is the planar N → ∞ limit in PCF. This model has a rich history of its comparison to QCD and it might provide an important example of the exactly solvable 2+1 dimensional bosonic string theory, similarly to the matrix quantum mechanical model of the 1+1 dimensional, c = 1 non-critical string theory proposed and solved in [40]. An exact and explicit solution for it was given in the case of infinite volume L but in the presence of specific "magnetic" fields [41].
As concerns the numerics, our algorithm converges very well for any length when N ≤ 3, but for N ≥ 4 it is very unstable for small enough length L, already at L 1 (which means for instance that we cannot really check the convergence to the conformal UV limit). Hopefully this instability has no direct physical meaning and is just a numerical artifact, due to a poor numeric accuracy, or to a bad choice of the iteration procedure for our NLIEs.
We believe that this method of deriving a finite system of NLIEs for integrable sigma models is general and powerful enough to work for much more complicated cases of AdS/CFT correspondence, such as the superstring on the AdS 5 × S 5 background dual to N=4 SYM theory, and the so called ABJM model where the Y -system was already discovered [14,38,39]. The wronskian quasiclassical character [17] and even the full quantum solution of the Hirota dynamics for AdS 5 /CFT 4 [18] are already available. It is left there to understand the very rich and complicated analyticity structure of Q-functions for short operators to complete the derivation of the AdS/CFT NLIE.