Pseudoscalar glueballs in the Klebanov-Strassler Theory

In this paper we describe a pseudoscalar subsector of the Klebanov-Strassler model. This subsector completes the holographic reconstruction of the spectrum of the lowest-lying glueball states, which are singlet under the global symmetry group $SU(2)\times SU(2)$. We derive the linearized supergravity equations for the pseudoscalar fluctuations and analyze their spectrum. The system of equation is shown to be compatible with six eigenmodes, as expected from supersymmetry. Our numerical analysis allows to reliably extract four of the corresponding towers. Their values match well the eigenvalues of the $0^{++}$ scalar states known from an earlier work. Assuming the masses of $0^{++}$ as a reference, we compare the lightest states of the holographic spectrum with lattice calculations in the quenched QCD at $N_c=3$ and $N_c=\infty$.


Introduction
In [1] Klebanov and Strassler derived a solution of the type IIB supergravity equations, which describes a holographic dual of a non-conformal N = 1 gauge theory. Contrary to initial expectations [2][3][4][5] this theory did not quite provide a dual of N = 1 supersymmetric Yang-Mills theory (SYM) (the simplest extension of the bosonic gauge Yang-Mills theory) with large number of colors N . Instead it gave an interesting and novel example of the so-called "cascading" theories, with a less conventional RG flow [6].
The difference of the Klebanov-Strassler (KS) cascading theories from the SYM in the IR limit is due to the spontaneous breaking of the baryonic U (1) B symmetry [7][8][9]. The spontaneous breaking is responsible for the presence of the Goldstone modes of U (1) B and for the corresponding tower of light states. Those light states mix with the light states of the SYM sector of the theory.
The KS theory is constructed as a solution of type IIB supergravity equations on AdS 5 × T 1,1 , where T 1,1 is the space S 3 × S 2 with a special choice of the metric compatible with N = 1 supersymmetry. T 1,1 space, and hence the dual gauge theory, has an SU (2) × SU (2) global symmetry. The spectrum of such a theory should be organized in its irreducible representations. The structure of the mutliplets of the conformal version of the KS model on T 1,1 (known as Klebanov-Witten theory [2]) was analyzed in [10]. Meanwhile, in the KS theory the conformal symmetry is broken by the flux of the 3-form field F 3 . The supergravity metric does not have isometries of the AdS 5 : it is a "warped" metric with a non-conformal logarithmic running of the AdS 5 radius, as well as different scaling of S 3 and S 2 in T 1,1 . Hence, in the KS theory superconformal multiplets of [10] break into multiplets of N = 1 supersymmetry. This implies, in particular, that massive bosonic states of this theory appear in pairs.
Fields of the pure SYM sector of the theory are singlets under the global SU (2) × SU (2). One might be interested in the structure of this particular sector as of a cousin of the simplest N = 1 SYM theory, and even as a more distant cousing of the pure glue (quenched) Yang-Mills. The spectrum of glueballs in the latter theory has been computed on a lattice [11,12] and one may wonder if a meaningful comparison with lattice results can be made.
The first observation is that, from the string theory point of view, the classical gravity solution that we use here is only a valid approximation for a very small string tension, which practically means that all states with spin higher or equal than two (with the exception of 2 ++ ) cannot be seen in the gravity description. Therefore one should focus only on the low-spin states, whose mass has the lowest order in the string tension.
The second observation is that classical gravity is only valid in the limit of large number of colors N c . This is less of an obstacle because the glueball masses are expected to weakly depend on N c : The goal of this paper is to complete the study of the SYM sector of the KS theory. Its structure is in general understood from a series of previous works [17][18][19][20][21][22][23]. See also [24] for the summary of this sector, [25][26][27] for some earlier works and [28][29][30][31][32][33] for extensions beyond the singlet regime, and beyond the KS theory. Almost all glueball states from the singlet sector were explicitly constructed as perturbations of the KS background. The masses of the lowest states in the conformal towers were computed. The only missing subsector in the previous analysis is that of the 0 −+ states, although it was explained in [23] that this sector should consist of six modes degenerate with six of the seven of 0 ++ scalars described in [17,19].
Here we will explicitly construct the pseudoscalar modes. In the holographic approach the spectrum is derived from linearized perturbations over a background (vacuum) solution. We will explain the structure of the J P C = 0 −+ perturbations in the KS background and present the corresponding linearized supergravity equations.
The perturbations are described by a complicated set of six coupled second order differential equations. However, it is more practical to write them with a help of an auxiliary field, which makes the equations more compact at the expense of introducing an additional first, or second order equation.
We find no obvious way to decouple the equations. Before analysing their spectrum numerically, we make a few consistency checks. First, the full number of equations obtained in the analysis of the perturbations is eight. We check that one of the eight equations can be derived from the remaining seven in a non-trivial way.
In deriving the equations we worked in explicitly gauge invariant setting. Gauge (diffeomorphism) invariance provides an additional check of the consistency of the derived equations. Finally, after fixing the gauge the system indeed can be reduced to six second order equations, matching the expected number of the 0 −+ modes.
We expect that the spectrum of the pseudoscalar modes matches the spectrum of 0 ++ calculated in [19] 1 . In our numerical simulations we were able to observe what appears as four towers of eigenvalues that indeed match well the masses of 0 ++ . However, our current numerical method does not resolve for the remaining two towers. A different approach with seven unconstrained modes produces two more towers, that do not match the spectrum in [19]. The nature of the latter two towers is not completely clear from our analysis and we expect that the subsequent studies will either confirm or discard them.
We show that a meaningful comparison of the holographic spectrum with lattice calculations can indeed be made. We can compare the lightest states in the spectrum with 0 ++ and 0 −+ of the SU (3) theory [12] and with its the SU (∞) extrapolation [14]. The states in the C-odd sector studied in [21,22] can be compared with the SU (3) values [12]. Some of them can also be compared with SU (∞) [15].
The comparative plot is shown on figure 4. In general the holographic calculation reproduces quite well the structure of the lattice spectrum, which confirms the expected weak dependence of the glueball spectrum on the details of the theory, such as N c and additional matter. Since the KS theory is supersymmetric, one sees more states than there are in the pure glue theory, so the holographic calculation may be expected to provide a reasonable estimate for the structure of the spectrum of a supersymmetric Yang-Mills theory. Further details of the comparison will appear in a separate work [34].
The remainder of this paper is organized as follows. In section 2 we summarize the necessary background material. Section 2.1 contains a brief review of the Klebanov-Strassler solution and section 2.2 comments on the dual theory. In section 3 we review the quantum numbers of particles from both the field theory and supergravity points of view. We also discuss the expected operator content in the pseudoscalar sector by comparing with the structure of superconformal multiplets on AdS 5 × T 1,1 from the analysis of [10]. Section 4 contains the main analytical results of our work. Following the earlier analysis of quantum numbers, in section 4.1 we specify the general ansatz for pseudoscalar fluctuations. In section 4.2 we present the corresponding linearized equations. We analyze their asymptotic behavior in section 4.3. Finally, in section 4.4 we check the consistency of the derived equations by checking their gauge invariance and analysing the number of independent modes. In section 5 we explain the results of the numerical analysis of the spectrum. Concluding remarks are made in section 6, where we also present the results of our comparison with the lattice data. The paper also contains two appendices A and B, in which we describe the asymptotic solutions.

Glueballs from a holographic model
Holographic approach [35][36][37] provides a powerful tool to analyze a few explicitly known interacting gauge theories in the regime of extremely strong coupling, for reviews see [38][39][40][41][42]. In particular, the spectrum of light states in a theory can be extracted from classical gravity equations. In this section we describe a specific gravity system found by Klebanov and Strassler [1], based on earlier developments in [2][3][4][5], dual to a N = 1 supersymmetric gauge theory with large number of colors at strong coupling.

Brief review of the Klebanov-Strassler theory
The Klebanov-Strassler (KS) model [1] is based on a solution of the equations of motion of type IIB supergravity [43]. The bosonic sector of this theory reduces to the Einstein equation here written in the Einstein frame, and equations for the matter fields (2.5) The following notations are commonly used: Here Φ, B 2 ≡ B M N and G M N are the NS-NS sector fields (respectively, the dilaton, the antisymmetric rank two tensor and the metric, whose Ricci tensor is denoted R M N ) of the type IIB theory in ten dimensions M, N = 0, 1, . . . , 9. C ≡ C 0 , C 2 and C 4 are the R-R fields (the scalar and antisymmetric tensors of rank 2 and 4). For compactness the equations are written in the differential form notations. In particular, * denote the Hodge dual. Additionally, there is a Bianchi identity for the 5-form fieldF 5 , The KS solution of the above equations starts from the metric where h(τ ) is called the warp factor and ds 2 6 is the metric of the deformed conifold, a six dimensional cone with the base S 3 × S 2 [44], Here τ is the radial coordinate on the conifold, measuring the coordinate distance away from its tip τ = 0. is the deformation parameter, controlling the curvature of the conifold at the tip. The base of the conifold can be parameterized by angular coordinates. Equation (2.9) uses the following basis of 1-forms on the base where e 1 ≡ − sin θ 1 dφ 1 , e 2 ≡ dθ 1 , e 3 ≡ cos ψ sin θ 2 dφ 2 − sin ψdθ 2 , e 4 ≡ sin ψ sin θ 2 dφ 2 + cos ψdθ 2 , e 5 ≡ dψ + cos θ 1 dφ 1 + cos θ 2 dφ 2 .
The solution also contains a non-trivial F 3 , the field strength of C 2 , which can be written as and a non-trivial H 3 , conveniently defined through (2.14) In the above equations g s is the string coupling constant, α = M −2 Pl is the string scale parameter and M is an integer explained below. Finally, there is a self-dual 5-form, which is decomposed as where with The above ansatz for the differential form is written in terms of functions F (τ ), f (τ ), k(τ ), l(τ ) and h(τ ). The explicit form of these functions is given by [1] while the warp factor is found from an integral The asymptotic behavior of this integral is with I 0 ≈ 0.71805.

Field theory interpretation and glueball states
Let us also briefly discuss the field theory interpretation of the dual geometry and explain how the spectrum of light particles of the gauge theory can be extracted from it. First, the metric explicitly breaks the conformal group SO(4, 2) of isometries of AdS 5 and hence corresponds to a non-conformal theory. It is only compatible with N = 1 supersymmetry. This geometry can be obtained from N D3-branes and M D5 branes, N M , such that four of the spacetime dimensions of the D5 coincide with those of the D3, and the remaining two are wrapped around the S 2 ∈ T 1,1 . Consequently, the dual theory is a SYM theory with SU (M + N ) × SU (N ) gauge group and global symmetry SU The gauge theory is coupled to two chiral superfields, A 1 , A 2 , in the representation (M + N, N ) and two anti-chiral superfields, B 1 , B 2 , in the representation (M + N , N ). The chiral and antichiral superfields transform as doublets of the respective SU (2) factor, while U (1) B is the baryon symmetry that acts as A i → e iα A i and B i → e −iα B i . The theory has a superpotential of the form (2.26) The gauge couplings of the two factors of the gauge group flow in opposite directions and the theory exhibits a "cascade" of Seiberg dualities whenever one of the couplings becomes infinitely strong: the spectrum of the theory changes and the direction of the flow flips. At the IR end of the cascade the theory becomes a strongly coupled SU (M ) SYM with light excitations corresponding to bound states of gluons and gluinos (the glueballs). The theory also possesses a nonzero condensate spontaneously breaking the baryon symmetry, which comes with an associated tower of light mesons and hybrid states of mesons and glue [7][8][9]. Below we will refer to all the above as glueballs.
The particle spectrum can be found from the poles of two-point correlation functions. In the AdS/CFT dictionary [36,37], the two-point functions can be computed from the solutions of the equations of type IIB supergravity, linearized over the background solution. See [45,46] for early examples of the glueball spectrum calculations. Our purpose here will be to select the perturbations of the type IIB fields corresponding to pseudoscalar modes, singlet under the global symmetry. The particle spectrum is commonly classified by the J P C quantum numbers, where J is the particle spin, P its parity and C the charge conjugation. Hence we will be interested in 0 −+ states.
In the next section we will discuss how the J P C quantum numbers can be determined for the supergravity fluctuations. We will also summarize our expectations about the dimensions of the dual gauge theory operators.
3 Symmetries and quantum numbers of glueballs

Quantum numbers
In this section we will review how the quantum numbers of the glueball states are determined in the holographic approach. We refer the reader to [45,46] for some original literature. The matter sector of the Klebanov-Strassler theory has SU (2) × SU (2) × U (1) B continuous global symmetry. The particle states should be classified by its irreducible representations. For example, under the SU (2) × SU (2), the states are classified by a pair of half-integer numbers (j 1 , j 2 ). They also carry charge under U (1) B baryon symmetry.
The glueball states of the Yang-Mills sector of the theory are singlets with respect to SU (2) × SU (2), so we will be interested in all states with j 1 = j 2 = 0. They also carry no baryon number. One should keep in mind, however, that this sector mixes with "hybrid" glueballs, containing A and B fields, charged under U (1) B . Due to the presence of these fields there is also a large non-singlet sector, classified by the A and B composition of the hybrids. Invariant combinations of A and B also contribute to the singlet sector.
The axial symmetry U (1) A and the U (1) R symmetry of the SUSY algebra, are anomalous in the KS theory [2]. The U (1) R is broken down to a Z 2M subgroup [47]. The vacuum further breaks the remaining symmetry spontaneously down to Z 2 . 2 Nevertheless, because it is only broken by the anomaly, U (1) R remains a convenient symmetry to classify the supermultiplet structure of the states. Since the superpotential must have R-charge 2, and supercoordinates transform as ϑ → e iα ϑ under the R-symmetry, A and B have charge 1/2, while the charge of the components is defined in such a way that the transformation of the superfield is homogeneous.
Let us now discuss the discrete symmetries. Realization of parity in the Klebanov-Strassler theory is straightforward. It reflects the spatial coordinates of the Minkowski factor, P : Charge conjugation involves complex conjugation of the fields and respectively, of their representations. Since fields A and B belong to (N +M, N ) and (N + M , N ) representations of the gauge group, respectively, one can combine the charge conjugation with an exchange of A and B. Note that the superpotential is odd under this transformation, W → −W. If, together with the exchange of the fields, the supercoordinates are rotated, ϑ → iϑ, the combined transformation will be a symmetry of the action: I : Following [8] we call this I-symmetry. It can also be understood from the embedding to the N = 4 theory, where A andB are combined into the fundamental multiplet of the SU (4) R-symmetry. I-symmetry is an unbroken Z 2 subgroup of that symmetry, mixing the factors of the continuous SU (2) × SU (2). We note that on the pure gauge sector, I-symmetry acts simply as a charge conjugation. So, for the purpose of this paper, C will be the eigenvalue of the I operation.
We now come to the discussion of the realization of the above symmetries in the type IIB SUGRA. First of all, the continuous SU (2) × SU (2) is the isometry of the two S 2 factors of conifold metric (2.9). Since we are interested in the singlet sector, we need to define a basis of differential forms that is invariant under the isometries.
• One forms. The only SU (2)×SU (2)-invariant one-form on T 1,1 is g 5 (2.11), so the full basis is provided by dx µ , dτ, g 5 . • Antisymmetric two-forms. There are four invariant two-forms on T 1,1 and the full basis is provided by where dots stand for external products of invariant one-forms (3.3).
• Symmetric two-forms are needed to construct metric fluctuations. The corresponding basis is provided by where dots denote terms obtained from internal products of invariant oneforms (3.3).
• All higher rank invariant antisymmetric forms can be obtained by evaluating the exterior products of the forms listed above.
The U (1) R symmetry acts by shifts of the coordinate ψ on T 1,1 , ψ → ψ + ζ. The metric and F 3 form have an explicit dependence on ψ, which means that this symmetry is broken by the KS background. The dependence is compatible with the anomaly [47]. In particular, since ψ is a double cover of a circle, there is a remaining Z 2 symmetry ψ → ψ + 2π.
The U (1) B symmetry is not realized as an action on T 1,1 . As U (1) R , it is not an isometry either, because it is spontaneously broken. Consequently, it generates a one-parametric family of deformations away from the KS background [7]. This family is called the baryonic branch [48].
The parity P in gravity theory is realized as inversion of the sign of the purely spatial coordinates x, but also as an action on the internal coordinates (which is a remanence of the higher-dimensional parity of ten-dimensional string theory). This implies that some gravity fields also transform. We assume that parity acts on the angular coordinates as Comparing this with the background solution, parity is a conserved quantity if B 2 (H 3 ) and C 4 (F 5 ) have negative "intrinsic" parity. Besides, the "axion" C is a pseudoscalar. I symmetry is an internal symmetry of the gauge theory. It acts only on the T 1,1 part of the geometry exchanging two S 2 spheres within T 1,1 . In terms of the coordinates, it swaps Besides, the I symmetry flips the sign of the F 3 and H 3 forms. It is useful to classify the invariant differential forms according to their P and I transformations. We summarize the charges of the forms on T 1,1 in table 1. Table 1: Parity, I and R charges of the differential forms.

Dual operators
Knowing the parity and I transformations of the forms, it is straightforward to construct an SU (2) × SU (2) invariant ansatz for the 0 −+ modes. We will do this in section 4.1. Before that we can anticipate the spectrum of operators, which will appear in the 0 −+ sector, by looking at the superconformal structure of the modes on AdS 5 × T 1,1 studied in [10].
It is useful to establish the charges of the forms under the R-symmetry. Since only ψ coordinate is transformed, the charge depends on the degree of the trigonometric function, which appears in the form. It is easy to check, that the linear combinations of forms shown in table 1 can be assigned either zero R-charge, or R = ±2.
On should be looking for operators of rational dimension, in particular those, corresponding to trλλ and trF µν F µν and trF µνF µν that combine into short multiplets of supersymmetry. In the superconformal theory the dimensions of the operators in the short multiplets are protected. Indeed, the analysis of the conformal dimensions of the 0 ++ modes studied in [19] showed that all of them have integer dimensions (table 2) and the spectrum contains ∆ = 3 and ∆ = 4 modes. This analysis was originally done in [49].  Table 2: Spectrum of SU (2) × SU (2) singlet 0 ++ scalar operators in the KS theory [49].
From the analysis of the spectrum of Kaluza-Klein modes [10,49,50] one knows that the ∆ = 4 modes come from the fluctuations of B 2 or C 2 proportional to the S 2 volume form ω 2 = g 1 ∧ g 2 + g 3 ∧ g 4 . Besides, fluctuations of the dilaton and the axion are ∆ = 4. Since dilaton Φ and axion C, as well as B 2 and C 2 have opposite parity, we expect to have two pseudoscalar modes of dimension ∆ = 4.
These modes have R = 0.
One pair of dimension three modes comes from the fluctuations of the metric, proportional to (g 1 ) 2 + (g 2 ) 2 − (g 3 ) 2 − (g 4 ) 2 and g 1 · g 4 − g 2 · g 3 . The first one is scalar, and the second -pseudoscalar. Another pair of ∆ = 3 modes comes from the fluctuations of the 3-form potentials proportional to g 1 ∧ g 3 + g 2 ∧ g 4 and g 1 ∧ g 2 − g 3 ∧ g 4 . Again, one of these forms is parity-even and the other parity-odd. Besides these modes have R = ±2. Therefore, we expect two ∆ = 3 modes in the pseudoscalar sector.
One can expect the ∆ = 3 modes to pair up with ∆ = 4 modes and form chiral multiplets of N = 1 supersymmetry. This is compatible with the structure of short superconformal multiplets of [10]. In particular, shortened Vector Multiplets III and IV of [10] contain ∆ = 3 and ∆ = 4 operators. The two multiplets differ by the sign of the R-charge of the lowest component, R = −2 for type III and R = 2 for type IV.
From the analysis of [49] one also observes modes with dimension ∆ = 7, ∆ = 8 and ∆ = 6. The first one comes from the second possible combination of the fluctuations of the 3-form potentials, proportional to g 1 ∧ g 3 + g 2 ∧ g 4 and g 1 ∧ g 2 − g 3 ∧ g 4 , while the other two modes come from linear combinations of fluctuations of traces of the metric on AdS 5 and on T 1,1 . Out of singlet superconformal multiplets of [10] only Vector Multiplet II can accommodate scalars with such a high dimension. This multiplet is not short, but nevertheless has a rational dimension. It corresponds to an unconstrained vector multiplet V of N = 1 symmetry, which accommodates four spin zero fields. Without conformal symmetry, this multiplet decomposes into on-shell massive vector mutliplet and two massive chiral multiplets. The vector component of the vector multiplet of dimension ∆ = 7 was found in [23]. To complete the Vector Multiplet II, one is missing a ∆ = 7 scalar and the longitudinal part of the vector mode.
Below we will explicitly identify the pseudoscalar modes and reproduce their spectrum from the linearized equations. Further details of the multiplet structure and the explicit form of dual operators can be found in [10,24,49].

Singlet Pseudoscalars of the KS theory 4.1 Ansatz for the modes
In this section we construct the ansatz for the SU (2) × SU (2) singlet 0 −+ modes in the KS theory. From the P and I transformation properties the most general form of the ansatz is Here we have considered fluctuations of the metric, R-R scalar, R-R 2-form, NS-NS 2-form potential and R-R 5-form, respectively. * 4 is the Hodge star operator in 4-dimensions with d * 4 dP = 4 P = −m 2 P . G 11 , G 33 and G 55 are the inverse coefficients of the (g 11 ) 2 , (g 33 ) 2 and the (g 55 ) 2 terms in metric (2.9) and h is the warp factor (2.22). The fluctuation of F 5 looks so complicated because it is constructed to satisfy the self-duality condition. The ansatz is constructed in terms of ten modes: 3 , a and A. However, we expect that there are only six physical modes, whose mass spectrum should match six of the seven towers of scalars in [17,19]. In particular, Bianchi identity (2.7) allows to solve explicitly for two modes, So modes φ 1 and φ 2 can be dropped in favor of φ 3 , A and a. Below we will write the linearized equations of type IIB supergravity over the KS background for the remaining eight fluctuations and show that the equations indeed describe six independent modes.

Linearized Equations
We will spare the reader the details of the derivation of the linearized equations and simply present the result. We will discuss the consistency of the derived system in section 4.4. Thus, the only independent equations generated by fluctuations (4.1)-(4.6) are the following: Here, we Fourier transformed the modes, substituting e ik·x for the spacetime coordinates with k 2 = m 2 4 . We also renormalized the mass eigenvalue, For the remaining parameters we used the convention α = g s = 1 and M = 2.
As a result, we arrive at a system of eight coupled ODE's (seven second order and one first order) for eight unknown functions. Equation  .11) is the C EoM (2.2) and the remaining three equations come from the Einstein equations (2.1). One can show that not all of the eight equations are independent. For example, one can take τ derivative of equation (4.14) and find a linear combination of the remaining equations, equal to the result of the differentiation.
This system of equations can be simplified if one uses gauge invariance of the gravity equations (see more details in section 4.4). One convenient gauge choice is to set a = 0. Since the eight equations are linearly dependent, one can drop equation (4.13). Equation (4.14) is only algebraic equation for A. It can be solved, so that one remains with six second order ODE's for six unknown functions.
However, the choice a = 0 is not convenient for the numerical analysis, because the coefficient of A in equation (4.14) vanishes for finite τ . The gauge that we will use for numerics below consists of choosing A = 0 and keeping (4.14) as an additional constraint. This constraint eliminates one of the seven modes leaving only six on shell.

Asymptotic behavior of the solutions
Before solving the system numerically we need to carefully analyze the asymptotic behavior of the solutions. In particular, we need to identify fourteen linearly independent solutions of the second order equations, determine which of those are regular, and eliminate those that do not satisfy the first order constraint. The asymptotic analysis will also allow us to extract the dimensions of the dual gauge theory operators and compare them to the analysis made in section 3.2.

UV
We start from the analysis in the UV limit, τ → ∞. The coefficients are not analytic functions in this limit, so that equations can be organized as an expansion in where λ 0 is the leading UV exponent, which will determine the dimension of the dual operator. P (τ ) is some analytic function of 1/τ . We will construct the asymptotic solution only in the leading exponential order, determining a few first expansion coefficients of P (τ ). The detailed results of the leading exponential expansion of the UV solutions can be found in appendix A. The asymptotic form of the seven second order equations can be read from formulas (A.1)-(A.7). The solutions to those equations are summarized in section A.1 (for dominant modes) and section A.2 (for regular modes). In section A.3 we check the compatibility of the regular solutions with the first order constraint (A.22), which is the UV asymptotic form of equation (4.14). Let us give a simplified version of the analysis here.
One first observes that the modes roughly separate into two groups: one containing φ 3 , C + 2 , C and a, and another with C − 2 , B 2 and B. (Specifically, one observes that for most of the modes, the modes within a group contribute the same order to the equations.) One can further decouple the pairs φ 3 and a, C + 2 and C, C − 2 and B 2 and the single mode B. This decoupling is only approximate, as in most cases the modes of the decoupled equations induce the other modes, but it does correctly determine the leading exponents of the modes. In this way one is left with the following system This decoupled system can be solved analytically: As expected there are fourteen independent modes. However, not all of them satisfy the constraint imposed by equation (4.14). In view of the above splitting of the modes one can write the asymptotic form of equation (4.14) as four simple constraints: Out of modes in equations (4.24)-(4.30) only α 1 , α 2 , α 3 , α 4 , α 5 , α 6 , α 8 and α 9 satisfy the above constraints. Besides, the following linear combination also of the two remaining modes also satisfies the constraints: It turns out that more modes satisfy the constraints imposed by equation (4.14), but to see that, one has to go beyond the naive analysis given here. What will be important for the numerical analysis is that the appropriate generalization of modes .  From the exponents associated with these modes we can obtain the spectrum of the dual operators. For this one changes the radial variable according to After an appropriate field redefinition one can cast equations (4.17)-(4.23) for every singular mode (4.37) in the form and use the standard formula, to compute the dimensions. We summarize our dimensions in table 3.
Mode α 2 α 3 α 6 α 8 α 9 α 11 AdS 5 mass, m  We see that the dimensions of the operators, in general, match our expectations. The only subtlety is the dimension ∆ = 5 mode coming from φ 3 . By construction, this mode should correspond to the longitudinal part of the vector appearing in the Vector Multiplet II in the classification of [10]. Consequently, a more natural dimension for it would be ∆ = 7. However, one observes that the φ 3 mode is unusual. It enters in the equations multiplied by the factor of 4 . Consequently, what couples to the operator of dimension 5 is 4 φ 3 , while φ 3 itself should couple to an operator of dimension 7. The remaining modes complete the Vector multiplets II, III and IV in [10].

IR
The analysis for small τ is technically simpler as the solutions are analytic there. We look for the solutions as a power series expansion (4.41) More generally, the coefficients of the expansion can contain log τ terms. We will see that all the relevant solutions (those that also satisfy constraint (4.14)) will not contain logarithmic terms. As in the large τ case, here we will give a simplified version of the analysis of the asymptotic solutions, while the full expansions are collected in appendix B.
For example, equations (B.1)-(B.8) in appendix B correspond to equations (4.7)-(4.14) with the leading asymptotic of the coefficients for τ → 0. In order to find the set of γ i , the exponents in the IR expansion of the modes, as in equation (4.41) we consider the following simplified system: The solution to this simplified system consists of the following modes: 52) C = β 9 τ −1 + β 10 , (4.53) B = 6 1/3 β 11 + 3 1/3 β 12 τ + β 6 τ 3 , (4.54) a = β 13 τ −2 + β 14 τ . (4.55) This solution demonstrates the leading modes present in the system. The complete solution of the asymptotic equations is constructed in appendix B, where the modes are in general, linear combinations of the above modes.
We also have to identify, which of the fourteen modes are regular. This is defined with respect to the action functional for the modes, which appears in equation (B.24) in section B.3 of the appendix. Our analysis shows that all the modes β i with even i are regular, besides all of them, except β 12 satisfy equation (4.14).

Are the equations correct?
Before moving on to the numerics we would like to discuss the consistency of our analysis. Equations (4.7)-(4.14) look fairly complicated and a reasonable question is how one can know that they are correct. Some non-trivial checks are provided by considerations that the eight equations are not all independent and by supersymmetry must describe only six physical towers of states. This is indeed true. Equation (4.14) is algebraic so that mode A can be eliminated. Some of the remaining modes can be eliminated by a gauge choice, for example a = 0, and the remaining seven equations can be shown to be linearly dependent. The analysis of the asymptotic solutions in the appendix also confirms this: we find that out of seven independent modes of the second order equations (4.7)-(4.13) only six satisfy equation (4.14), which lead to six physical modes.
Gauge invariance provides another non-trivial constraint on the equations. 4 Since we did not fix the gauge, we can check, whether system (4.7)-(4.14) is gauge (diffeomorphism) invariant. Infinitesimal diffeomorphisms are generated by Lie derivatives, which for an arbitrary rank tensor read where ξ is a vector, which for the KS theory we can write as ξ =ã (x, τ )ê ψ . In terms of the modes one finds that (4.56) One can check that the above transformations indeed generate a symmetry of system (4.7)-(4.14). Another check is provided by the analysis of the dimensions of the asymptotic behavior of the fields and by the dimensions of the dual operators. We find a set of integer dimensions, summarized in table 3, which match our expectations.
Finally, the spectrum of the equations should provide one more consistency check of the equations. It is expected to match the eigenvalues of [19] without the tower of [23]. We will explain the results of our numerical analysis in section 5. Studying the spectrum of this particular system of seven coupled equations turns out to be a difficult exercise, and our success is only partial so far. In our numerical studies we are able to see four towers of eigenvalues found in [17] and [19], whose numerical values match quite well our spectrum. As for the two remaining towers, a direct method, which we used, was not able to detect them. Instead, if we add another degree of freedom by removing the constraint, we see two additional towers that do not match the eigenvalues of [19], so their status is unclear. In the meantime, we see indications that our direct method might not have sufficient resolution, which might explain the absence of the remaining modes.

Numerical analysis and spectrum
In this section we describe the results of our numerical analysis of the system of equations (4.7)-(4.14).
The basis of our analysis is the midpoint method, also used in [19], for the analysis of the 0 ++ glueballs. This technique is a generalization of the shooting method. In the regular shooting method, one replaces the boundary problem by the initial value problem and solves equations either going from IR (small τ ) to UV (large τ ) or vice versa. The midpoint method involves shooting from both the UV and IR ends and gluing the solutions at some intermediate point τ mid . Consequently, for n equations for n fields, one sets 2n regular initial conditions in the UV (for the fields and derivatives) and 2n regular initial conditions in the IR. The equations are solved until τ mid , where one builds a square matrix of dimension 2n × 2n, where a i is the set of fields, i = 1, . . . , n, and ∂a i is the set of derivatives. The eigenvalues of the system of equations are values ofm 2 , for which matrix γ has a zero eigenvalue. In particular, the determinant of γ, as a function ofm 2 , will change the sign at the eigenvalue point. In this work we will be solving seven equations for six initial conditions in the UV and the IR, satisfying a first order constraint. One possibility in this case is to remove one of the fields and its derivative from the matrix γ, to form a 12 × 12 matrix. Alternatively, one can calculate a rectangular γ, 12 × 14, keeping all the fields, and compute the singular value decomposition. At the eigenvalue point, the smallest singular value should have a zero. This can be observed as a cusp in singular value as a function of m 2 .
Yet another possibility that we will explore is to solve the system of seven equations with seven unconstrained regular initial conditions and look for zeroes of the determinant of a 14 × 14 matrix. In this case one can expect an additional tower of eigenstates, which could be separated and removed when the spectrum is compared with that of [19].
We first discuss the results obtained by the calculation of singular values of the 12×14 matrix. This approach, based on the analysis of the six physical modes showed the best convergence with the spectrum of 0 ++ for smaller mass values. From now on we renormalize mass parameterm 2 to the conventions of [19].
The results of the calculation can be found in table 4, where the values are compared with the values of [19]. We see that the found masses match rather well the values in [19], however, they do not reproduce the whole spectrum. See figure 1, for example.  Table 4: The m 2 values obtained by the singular value method and the corresponding values of the 0 ++ masses calculated in [19]. The results were obtained using the Mathematica Explicit Runge Kutta method fixed step δτ ≤ 5 × 10 −4 , τ U V = 20, τ mid = τ IR = 0.1, for 1 ≤ m 2 ≤ 23; τ mid = 2, for 23 ≤ m 2 ≤ 38; τ mid = 5, for 38 ≤ m 2 ≤ 46. The value m 2 = 0.543 does not appear as a cusp, but rather as a broad minimum, for some values of τ mid .
By analyzing the periodic pattern of the eigenvalues on figure 1 and similar plots, we conclude that we observe four of the six expected towers. The singular value shows cusp-like dips around some of the values of [19]. For for m 2 ≤ 23 the best convergence so far was achieved by selecting τ mid = τ IR , that is effectively using the shooting method. However, for m 2 > 23 the method starts failing on some of the eigenvalues. The sharp cusps of some of the towers are replaced by smoother minima and eventually disappear. This explains why in principle, the remaining two towers are not seen by our calculation. For the four towers the problem can partially solved by moving τ mid to higher values. In particular, to see the eigenvalues 38 ≤ m 2 ≤ 46, we used τ mid = 5. The general experience tells that whenever we have a sharp dip, the eigenvalue matches very well with the table of [19]. Smoother shallow minima probably indicate the problem with the numerics. Their positions typically do not match the values in [19] so well.
In an attempt to see the remaining modes, we also tried to solve the system, including the non-physical mode, that does not satisfy constraint (4.14). In this case we computed the determinant of the 14 × 14 matrix, and there is no need to through away any information, as would be the case of a 12 × 12 determinant. In this approach we used the τ IR = 0.1, τ U R = 20, τ mid = 1, δτ = 10 −2 , although In the 14 × 14 case, we are able to see more eigenvalues and the periodicity patterns indicates to the presence of six towers. Out of the six towers one can distinguish the four towers discovered in the singular value approach. The two methods give compatible sets of values. However, the two new towers (their eigenvalues appear in bold in table 5) are not the towers seen in [19]. This is clearly seen on figure 3, where a set of zeroes of the determinant insert in the gaps of the heavier part of the 0 ++ spectrum.
Comparison of the two methods allows to separate the spectrum of [19] into groups of four and two towers, in accordance with the separation of values in table 5, although we have to assume that the singular value methods captures all the eigenvalues of the four towers. This is apparently true for higher masses, but can be more subtle for the lightest ones. For example, for τ mid = 5, in the singular value method, we see a broad minimum around m 2 = 0.543, which could be an onset of a lighter mode in the four towers, or a mode from a different tower.     Table 5: The values of m 2 computed from the 14 × 14 matrix which includes an unphysical boundary condition. Apart from the four towers that match those, reported in table 4 the spectrum shows two additional towers (shown in bold) which do not match the eigenvalues in [19]. The calculation is done for δτ = 10 −2 , τ IR = 0.1, τ U V = 20 and τ mid = 1.

Conclusions and discussion
In this work we have analysed the 0 −+ pseudoscalar SU (2) × SU (2) singlet fluctuations of the Klebanov-Strassler theory. We constructed the ansatz for the perturbations of the bosonic type IIB supergravity fields over the KS background, with 0 −+ quantum numbers, and derived the linearized equations.
The linearized equations that we have found can be reduced to a system of six coupled second order equations, so they describe six independent towers of 0 −+ states, however, this system has a more compact presentation in terms of seven second order equations, with a first order constraint. We analyzed the asymptotic behavior of the compact system and recovered asymptotic expansions of the solutions at the two ends of the KS geometry (UV for τ → ∞ and IR for τ → 0). This allowed to extract the spectrum of the 0 −+ operators in the dual gauge theory and complete their classification in terms of the multiplets of a N = 1 superconformal theory dual to strings on AdS 5 × T 1,1 [10].
We discussed different self-consistency checks of our complicated system, which included a non-trivial linear dependence of eight differential equations, gauge invariance and consistency of the operator spectrum.
As an ultimate consistency check we have been analyzing the numerical mass spectrum of the six 0 −+ modes produced by the linearized equations and compar-ing them with the known spectrum of 0 ++ modes [19]. We were able to show that our equations reproduce four of the six towers of the 0 ++ modes with a reasonable accuracy. However, our main method of the numerical analysis was not able to resolve for the remaining two towers. We expect that a more powerful numerical approach should be able to solve this problem. We left this for a future work [34].
As an alternative method we have been analyzing the spectrum of seven equations without the constraint. In this case we were able to capture six independent towers of eigenvalues: four compatible with the previous analysis and with the 0 ++ spectrum, and two other towers quite different from the remaining eigenvalues of [19]. We have to interpret the presence of a new pair of modes. One of those modes could be an unphysical mode, which does not satisfy constraint (4.14). But there are two modes, so there is a possibility that the spectrum of our equations is different from that of [19]. So far, our numerical analysis does not allow to definitely conclude this. One problem is that we have not seen all the seven modes. Moreover, analyzing the different initial conditions, we observe that even if one starts from an unphysical initial condition at one end, the system will evolve into the physical subspace, satisfying the constraint, in the other end. So, the solution is projected onto a smaller subspace in the evolution. This is demonstrated by the plots on figures 5 and 6 shown in appendices A and B.
It is interesting to compare the holographic prediction of the spectrum with the spectrum of the pure glue SU (3) theory calculated on the lattice [11] and its extrapolation to SU (∞) [13][14][15]. One can observe a reasonable match of the pattern of all of the six lowest glueballs of the lattice spectrum, which can in principle be captured by the classical gravity approximation, including 0 ++ , 0 −+ , 2 ++ , 1 +− , 1 −− and 0 −− . Our comparison is shown on figure 4.
To compare, we normalized all the glueball masses to the value of the 2 ++ state mass. In the C = −1 sector, the match is quite good (right panel of figure 4), as previously discussed in [22,24]. There seem to be somewhat less lattice data in this sector, so we only compare with the SU (3) results. See however [15].
In the C = +1 sector, the lightest scalar on the lattice, corresponding to the trF µν F µν operator, matches the second lightest state in the 0 ++ spectrum of [19]. One has to remember that the KS theory is supersymmetric, so the lightest expected 0 ++ state in the supersymmetric theory is dual to the gluino bilinear trλλ, not present in the pure glue theory. The same is true for the 0 −− state that appears in the C = −1 sector, and some heavier states that we do not show here. Moreover, from our analysis, we find it plausible, that the state that roughly matches the first excited state 0 * ++ on the lattice is the first excited state of the lowest bosonic holographic 0 ++ . In other words, the matching can be extended to the excited states as well. There is also a holographic state that matches the lightest pseudoscalar on the lattice. This state is dual to the operator TrF µνF µν . We can also speculate on the possibility that the two towers of pseudoscalar states seen in our analysis, which do not match the spectrum 0 ++ computed -- We only compare states of low spin, which can be obtained from the holographic analysis (0 P C , 1 P C and 2 ++ ). All masses are normalized to the mass of the 2 ++ state. Lattice data are taken from [14] in the positive C sector and from [12] in the opposite one. The column BHM refers to the result of [19], MR compares the results of the present work, DMS and BDKS are the results of [22] and [21] respectively.
in [19], are, in fact, the correct values. The lightest states of this spectrum are shown in the last column of the left panel of figure 4. The values would give a comparable match for the lightest states in the pure glue SU (∞) sector. Besides the 14 × 14 spectrum turns out to be more fittable with quadratic fits in comparison to the spectrum of [19]. This in particular supports the expectation of the third and the fifth states being the first excited states, of trλλ and trF µν F µν respectively. In the following work, with a better numerical resolution, we expect to validate or discard our conjectures.

A UV asymptotics
In this section we analyze the asymptotic behavior of the solutions to the system of linearized equations (4.7)-(4.14) for large τ (UV regime). We will construct the solutions as expansions in 1/τ , working in the leading exponential approximation, as in equation (4.16). In such an approximation equations (4.7)-(4.13) take the following asymptotic form.
In section 4.3 we presented a simplified version of this system, see equations (4.17)-(4.23), which was solved analytically and the modes were classified according to equations (4.24)-(4.30). Here we will construct solutions to the full system.
Below, in section A.1 we will construct the singular modes of the seven above equations, while in section A.2 we will do the same for the regular modes. The modes will be constructed in the leading exponential and up to the next-tonext-to-leading order (NNLO). In section A.3 we justify the classification of the modes as singular and regular and analyze their compatibility with the constraint imposed by equation (4.14).

A.1 Singular UV solutions
The UV asymptotic equations (A.1)-(A.7) possess seven singular and seven regular modes. One can choose the following basis of the singular modes, which roughly follows the classification (4.24)-(4.30): Note that this mode has linear τ behavior of C + 2 and is "locked" with the constant C mode.
• α 9 mode This is a subtle mode, because the expansion of C + 2 occurs in the subleading exponential order, while the leading constant solution drops from most of the equations, which only depend on the derivative.

A.3 Analysis of the UV modes
Now we would like to check the compatibility of the regular UV solutions with equation (4.14). The leading exponential asymptotic of this equation reads As expected the physical system contains only six independent modes. One can make some numerical tests of the exact constraint (4.14). For each regular α i we plotted the (logarithm of the) ratio of the left hand side of equation (4.14) (which is supposed to vanish) over the sum of absolute values of the terms of the same equation. The results are shown on figure 5. As can be seen from that figure, all but two modes (α 10 and α 14 ) satisfy equation (4.14) with good accuracy. In fact, in the limit τ → 0, the numerical solution evolves α 10 and α 14 to modes that satisfy the constraint in the IR. Similarly, one can show that the above linear combination of α 10 and α 14 always satisfies the constraint.
To conclude the analysis of the UV asymptotic we need to comment on the regular versus singular classification of the modes. For this we evaluate the action functional for the modes substituting the fields in the type IIB action. This gives: For this expression to be finite, the following conditions must be satisfied: We conclude that the above separation of the modes into singular and regular lists was correct.

B IR asymptotics
In this section we will summarize the asymptotic solution of linearized equations (4.7)-(4.14) in the regime τ → 0 (IR).
First, we summarize the form of the equations when only the leading τ → 0 terms in the coefficients are taken into account. We also neglect the terms which are subleading inhomogeneous (e.g. a ∼ a/τ 2 a/τ ).
Note that we diagonalized the pair of the equations containing the second derivatives of C ± 2 . The pair of the resulting equations imply a very simple equations for the sum P = C + 2 + C − 2 : We have not attempted to completely solve the reduced system, but instead considered a simplified system (4.42)-(4.48), which gave the same spectrum of exponents β i , according to the notation in (4.41). Below we list the expansion of the solutions to the full system.
In the IR limit the expressions are typically more bulky, so we restrict ourselves to either the leading order (LO) or NLO. The general principle of organization of the expressions below is that only those orders are listed that can be fit in a single line. Indeed, more complex form of the coefficients typically mean that they correspond to higher orders and are likely to be neglected. In the numerical calculations we use fuller expressions if necessary, to control the precision. The expressions here are listed in order for the results to be reproducible.
We label the modes in accordance to the classification of the simplified system (4.49)-(4.55). We will first list the modes with singular behavior at τ → 0 in section B.1 and then the modes with regular behavior in section B.2. We will justify the classification in section B.3.

B.1 Singular IR solutions
The system of seven second order differential equations has seven modes with singular behavior at τ → 0. These modes are analogous to β i with odd i of the simplified system (4.49)-(4.55). Here we list the modes following their simplified classification • the mode similar to β 1 in (4.49) is • β 7 mode (this mode in fact is a combination of β 7 and β 13 , while the proper • β 13 mode (this mode mixes with β 7 and hence has logarithmic coefficients)

B.2 Regular IR solutions
In this section we list the seven regular asymptotic solutions of system (4.7)-(4.14) in the IR limit. In terms of the modes (4.49)-(4.55) these contain  • β 8 mode • β 12 mode (note that this mode contains logarithmic terms in the expansion)

B.3 Analysis of the IR modes
In the previous sections we separated the IR modes in the singular and the regular ones. This separation follows from the analysis of the asymptotic form of the action functional for these modes, which we deduced from the type IIB action.
The relevant terms in the functional read Here index µ labels only the Minkowski directions, while index M can also take value M = τ . The τ integral of this function is convergent for small τ if These, however, are sufficient, but not necessary conditions. For, example, one can substitute the conditions for C ± 2 by a weaker condition P = C + 2 + C − 2 < τ 3/2 .

(B.32)
It is important that the most singular term of C ± 2 in (B.24) only depends on the sum P .
One can check that all the modes β i with even i satisfy the above conditions with a subtlety, that the modes β 4 and β 12 only satisfy the weak condition.
We should also find out which of the modes are physical, that is satisfy the constraint imposed by equation (4.14), whose leading asymptotic form is equation (B.8). We checked that up to the NNLO all the modes β i with even i, except β 12 do satisfy equation (4.14), which confirms the fact that the system we are studying contains only six physical modes.
Finally, one can make a numerical check of equation (4.14), similar to the one made in appendix B. As demonstrated by figure 6, six of the seven regular mode satisfy the constraint with the numerical accuracy. Again, one sees that the mode that does not satify the constraint in the IR, tends to do so at the opposite end.