Extensive numerical study of a D-brane, anti-D-brane system in AdS5/CFT4

In this paper the hybrid-NLIE approach of [38] is extended to the ground state of a D-brane anti-D-brane system in AdS/CFT. The hybrid-NLIE equations presented in the paper are finite component alternatives of the previously proposed TBA equations and they admit an appropriate framework for the numerical investigation of the ground state of the problem. Straightforward numerical iterative methods fail to converge, thus new numerical methods are worked out to solve the equations. Our numerical data confirm the previous TBA data. In view of the numerical results the mysterious L = 1 case is also commented in the paper.


Introduction
In this paper in the context of AdS/CFT [1,2,3] we study numerically the ground state of a pair of open strings stretching between two coincident D3-branes with opposite orientations in S 5 of AdS 5 × S 5 . The main motivation for the study is that according to string-theory the ground state of such a configuration is expected to be tachyonic for large values of the 't Hooft coupling [4]. In our work we rely on the perturbatively discovered and later "all loop conjectured" integrability [5] of both the AdS 5 × S 5 super-string and the dual large N gauge theory. For string configurations with D-branes integrability enabled one to describe string configurations ending on different types of D-branes as 1-dimensional integrable scattering theories with boundaries [6,7,8,9,10]. This formulation of the problem makes it possible to go beyond the approaches of perturbative gauge and string theories being valid for small and large values of the 't Hooft coupling respectively, and to determine the exact spectrum of the model at any value of the coupling constant. However, even with the help of the powerful techniques offered by integrability, the exact analytical solution of the problem is not possible. Remarkeble analytical results are available in the small [11,12,13,14,15,32] and large [16,17,19,18] coupling regimes, but the determination of the spectrum at any value of the coupling constant can only be carried out by high precision numerical solution [20,21,22] of the corresponding nonlinear integral equations.
In our paper we consider the case, when the two D3-branes are giant gravitons [23], namely they carry N units of angular momenta in S 5 . If the S 5 of AdS 5 × S 5 is parametrized by three complex coordinates X, Y, Z satisfying the constraint: |X| 2 + |Y | 2 + |Z| 2 = 1, then our D3-brane and anti-D3-brane are given by the conditions Y = 0 andȲ = 0 respectively. They wrap the same S 3 , but with opposite orientation. As a consequence of Gauss law such a system can support On the large N gauge theory side a Y = 0 brane is represented by a determinant operator [24] composed of N copies of the field Y : where a i and b i are color indices and ǫ is a product of two regular epsilon tensors ǫ a 1 ···a N b 1 ···b N = ǫ a 1 ···a N ǫ b 1 ···b N . The local operator corresponding to an open string ending on a Y = 0 giant graviton can be obtained from (1.1) by replacing one Y field with an adjoint valued operator W [25]: (1. 2) The gauge theory description of a pair of open strings stretching between two Dbranes is given by a double determinant operator, such that the string insertions W and V connect the two determinants of the Y fields 1 : Unfortunately, the precise gauge theory dual of the DD-system of our interest is not known. In [4] it was approximated by a double determinant operator similar to (1.3), but in one of the determinants the Y fields are replaced withȲ fields 2 : This observation allows us to apply the boundary Thermodynamic Bethe Ansatz technique (BTBA) [26] to each open string separately. The necessary ingredients of this technique are the boundary reflection factors [8,27,28,29] and the asymptotic Bethe equations of the problem [4]. Unfortunately, apart from some very special cases [30,31,32], it is still unknown how to derive BTBA equations for a general non-diagonal scattering theory in the context of the thermodynamical considerations of [26]. This is why in [4] the Y -system [33,34,35] and the related discontinuity [36] 1 The ground state of such string states is BPS. 2 According to the argument of [4] the correct state might have other structures involving the fields Y andȲ , but should be similar to the double determinant form (1.4) and the mixing with other fields seem to be suppressed at large N equations supplemented by analyticity assumptions compatible with the asymptotic solution [37] were used to derive BTBA equations for the nonperturbative study of the ground state of the DD-system [4].
The BTBA description of the system is an infinite set of nonlinear integral equations. The numerical solution of the equations [4] showed that the ground state energy is a monotonously decreasing function of the coupling constant 3 g. The analytical investigation of the large rapidity and large index behavior of the Y -functions of the BTBA revealed that the usual BTBA description of the system breaks down when the energy of an open string state with angular momentum L gets close to the critical value: E c (L) = 1 − L. This point was interpreted in [4] as a transition point where the ground state becomes tachyonic. Approaching the critical point the contribution of infinitely many Y -functions must be taken into account to get accurate numerical result for the energy 4 . This fact suggests reformulating the finite size problem in terms of finite number of unknown functions. The possible candidates could be the FiNLIE [46], the quantum spectral curve (QSC) [47,48] or the hybrid-NLIE (HNLIE) [38] formulation of the problem. Since at present it is not known (not even for the Konishi problem) how to use the analytically very efficient [19,14] QSC method for numerical purposes, we choose the HNLIE method to reformulate the finite size problem of the DD-system. In this paper we transformed the infinite set of boundary TBA equations [4] into a finite set of hybrid-NLIE type of nonlinear integral equations. We perform the extensive numerical study of these type of equations in order to get as close to the special E BT BA = 1 − L critical point as it is possible.
Our numerical results reproduce the numerical evaluation of the boundary Lüscher formula [27,39] in the linear approximation, and the numerical BTBA results of [4] as well. These numerical comparisons give further numerical checks on the hybrid-NLIE technique of [38]. Unfortunately, as g increases new local singularities enter the HNLIE formulation of the problem. Thus we could not approach very close to the critical point. Nevertheless, in the range of g where physically acceptable numer-3 Throughout the paper the relation between g and the 't Hooft coupling λ is given by: λ = 4π 2 g 2 . 4 This means that the usual truncation procedure for solving the infinite set of TBA equations is not applicable to such a system ical results were obtained, the HNLIE results could give higher numerical precision than that of the BTBA and also some interesting facts could be read off from our numerical data.
During the numerical solution of the HNLIE equations straightforward numerical iterative methods failed to converge, thus new numerical methods were worked out to solve the equations.
The ground state of the L = 1 state is a very special case, since there the critical point is right at g = 0 and so far neither perturbative field theory computations nor the boundary Lüscher formula could provide a finite quantitative answer to the anomalous dimension of this state. On the integrability side the HNLIE approach allows us to get some numerical insight into this problem. The outline of the paper is as follows: Section 2. contains the HNLIE equations. In section 3. the numerical method is described. In section 4. the numerical results and their interpretation is presented. Section 5. contains some comments on the mysterious L = 1 case and finally our conclusion is given in section 6. Various notations, kernels of the integral equations together with the necessary asymptotic solutions are placed in the appendices of the paper.

The HNLIE equations
In this section we transform the previously proposed BTBA equations of [4] for the ground state of our D-brane anti-D-brane system to finite component hybrid-NLIE equations. For presentational purposes we group the equations into 3 types.
There are TBA-type equations, horizontal SU(2) hybrid-NLIE type equations, and vertical SU(4) hybrid-NLIE type equations. They together form a closed set of nonlinear integral-equations, which are solved numerically in this paper. As it is usual, structurally the equations consist of source terms plus convolutions containing coupling dependent kernels and nonlinear combinations of the unknown functions. The objects appearing in the arguments of the source functions are subjected to quantization conditions, but similarly to the boundary TBA description [4], due to the u → −u symmetry of the problem they are tied to the origin of the complex plane, thus extra quantization conditions are unnecessary to be imposed, since they are automatically satisfied by symmetry. Since these source term objects have fixed positions their positions are exactly the same as that of their asymptotic counterparts. This fact saves us from the tedious computation of the source terms, since if we take the difference of the exact equations and their asymptotic counterparts the source terms cancel from the equations. To be pragmatic and save time and space, the equations will be presented in such a difference form. Thus for any combination f of the unknown functions, we introduce the notation δf is the asymptotic counterpart of f . Having introduced this notation, we start the presentation of the equations by the TBA-type part. For the labeling of the Y -functions we use the string-hypothesis [40] based notations of [41]. For the presentation of the equations a few more notations need to be introduced: ).
(2.1) For later numerical purposes we re-parametrize log Y Q by the formula: 2) such that c Q is the constant value of log Y Q at infinity and ε is minus twice the energy 5 : ε = −2 E BT BA . From the TBA equations of the problem [4], it follows that δc Q = c Q −c o Q ≡ δc is Q-independent, and for small g, logȳ Q is a smooth deformation of its asymptotic counterpart, such that δ logȳ Q tends to zero at infinity. 6 Using this decomposition the following notations are need to be introduced: Then the TBA-type equations take the form: The log multiplier of ε in (2.2) is chosen not to modify the constant term in the large u behavior and to satisfy , which is the LHS of an important Ysystem equation divided by its asymptotic counterpart. 6 log Y Q cannot be considered as smooth deformation of log Y o Q , because log Y Q −log Y o Q ∼ ε log |u| diverges for large u at any g. On the other hand logȳ Q − logȳ o Q is small for any u at small g and tends to zero at infinity.
For Y 1 the modified hybrid form [42] of the BTBA equations is used, where p 0 is the index limit starting from which the upper part of the TBA equations is replaced by an SU(4) NLIE of [38] (See figure 1.). For any kernel vector appearing in the TBA equations Ω(K Q ) denotes the residual sum ∞ Q=p 0 −1 L Q ⋆K Q , and following the method of [42] for p 0 ≥ 4 it can be expressed by next to nearest neighbor Yfunctions as follows: where r m = log(1 + Y m|vw ), the kernels s, s 1 2 , σ1 2 are hyperbolic functions [42], while the other TBA kernels can be found in appendix A. As a consequence of the re-parametrization (2.2) the two constants δc and ε also become part of the set of where for any kernel K: CK(u) denotes the constant term in the large v expansion of K(u, v). 8 As we mentioned −ε/2 is the TBA energy, thus (2.12) gives the energy formula in our formulation of the finite size problem. The asymptotic forms of the Y -functions necessary for the formulation of (2.4-2.13) are listed in appendix D. To close the discussion of the TBA-type equations we note that equations (2.7) and (2.8) determine Y ± up to an overall sign factor. The sign factor can be fixed from the asymptotic solution and its value is −1. Thus the fermionic Y-functions can be expressed in terms of the LHS of (2.7) and (2.8) by the formula: 14) The horizontal SU(2) wing of the TBA is resumed by an SU(2)-type NLIE [43,38], which in our case takes the form: where 0 < γ < 1/2 is a contour shift parameter, the kernel G is given by (B.6) and the asymptotic solution for b andb is given in appendix D 9 . The upper SU(4) NLIE 7 Here the ⋆ notation means simply integration from −∞ to ∞. 8 Here we note that only the dressing kernel has logarithmically divergent term in its large v expansion, all the other kernels has either constant term or they simply vanish at infinity. 9 In practice b andb are complex conjugate of each other.
of [38] is attached to the TBA equations at the p 0 -th node. The upper NLIE is for 12 complex unknown functions: b A and d A , A = 1, ..., 6. They are combinations of the T -functions of the upper wing SU(4) Bäcklund-hierarchy [38]. Their relations to the unknowns introduced in [38] are given by (B.3,B.4) in appendix B and their asymptotic forms are given in appendix C. Using the notation B A = 1 + b A and D A = 1 + d A , the equations they satisfy take the form: where the kernels are given in (B.5-B.12). The shifts in the kernels which is equivalent to fixing the lines on which the NLIE variables live, are chosen in a symmetrical way and fixed as follows: 3 , γ 1 , γ 2 , η 3 , η 1 , η 2 , η In practice this reduces to half the number of SU(4) NLIE variables. The vectors E A andĒ A are conjugate to each other and they give the TBA input into the upper NLIE. To give their form we introduce the notations:

25)
(2.28) The last set of equations gives, how the upper NLIE variables couple to the TBA part of the equations.
whereb 2 andd 2 are from the re-parametrization of b 2 and d 2 : (2.37) with η = ±1 being a global sign factor. Similarly to the definition ofȳ Q , also here the benefit of usingb 2 andd 2 is that, for small g, logb 2 and logd 2 are smooth deformations of their asymptotic counterparts, and in addition δlogb 2 and δlogd 2 vanishes at infinity, which is necessary for the convergence of certain integrals. The decompositions (2.36),(2.37) are chosen to be compatible with the functional relation b [38]. Equations (2.4)-(2.35) constitute our complete set of nonlinear integral equations, which governs the finite size dependence of the vacuum of our D-brane anti-D-brane system.

The numerical method
Here we describe our numerical method for solving the hybrid-NLIE equations presented in the previous section. During the iterative numerical solution of the equations we faced with very serious convergence problems, which forced us to work out such a method that overcomes all the difficulties emerged. Our numerical method can be applied to solve other type of nonlinear integral equations as well. The power of the method is shown by the fact that numerical convergence was reached even in such cases, when the solution was physically unacceptable. The numerical method consist of two main steps, namely: • Discretization of the equations • Iterative solution.
The first step involves the discretization of the unknown functions and kernels, furthermore the discrete approximate representation of the convolutions. Having carried out the appropriate discretization of the problem, the equations are considered as large nonlinear algebraic set of equations. Thus eventually instead of integral equations we solve discrete algebraic equations. In this paper we will present two methods to solve them numerically.

Discretization of the problem
The discretization serves two goals. First it allows us to reduce the numerical problem from solving integral equations to solving algebraic equations. Second choosing the discretization points appropriately it reduces the number of degrees of freedom as much as it is possible to reach the desired numerical accuracy. In our actual numerical computation instead of u of section 2. we used the new rapidity u → u g , because with such a scaling almost all the rapidity difference dependent kernels become g independent. Thus for example Y ± (u) will be defined in [−2g, 2g]. To decrease the number of discretization points the u → −u symmetry of the problem is exploited. This means that the Y -functions are to be discretized only on [0, ∞] or [0, 2g] and as for the NLIE variables it is enough to discretize the b-and d-type variables on [0, ∞]. Since we do not want to introduce any cutoff in the rapidity through the transformation formula: This formula is chosen such that the branch point 2g corresponds to t = 1 for any choice of the parameter a, where a is a global scaling factor which changes from unknown to unknown. We chose the values as follows: for Y 1|w a = 1, for b andb a = 2, for Y Q and Y Q−1|vw a = Q, for η 1 andη 1 a = p 0 , and finally for the b-and d-type NLIE functions a = p 0 . These values are chosen to preserve the smoothness 10 of the transformed functions in the finite interval. After this transformation all of our unknown functions live on a finite interval. To discretize them we used piecewise Chebyshev approximation. This means that we divide the finite interval into subintervals and on each subinterval the functions are approximated by a given order Chebyshev series. The choice of subintervals is not equidistant. The subintervals are placed more densely around the branch points, since the function x(u/g), which governs the decay of the massive Y Q -functions, has the largest change around this point. The advantage of the Chebyshev approximation is that if the function is smooth enough on the subinterval, the coefficients of the Chebyshev series decay rapidly and the order of magnitude of the last coefficient allows us to estimate the numerical errors of the procedure. Now we describe the discretization method in more detail. Our functions are defined on either [0, B(Q)] or on [0, 2g]. This is why two type of subinterval vectors are defined A Q and A ± , such that the endpoints of the subintervals of [0, B(Q)] are put into the vector A Q and the endpoints of the subintervals of [0, 2g] define A ± . Let l k be the order of the Chebyshev approximation, then using the general rules of the Chebyshev approximation, a given function where now the vector A stand for either A Q or A ± , furthermoreT j−1 are a slightly modified Chebyshev polynomials 11 are the Chebyshev coefficients of the function f , which can be computed from the sampling points of the Chebyshev approximation: by the simple formula: where c (i) (l k ) are the zeros of the l k order Chebyshev polynomial: andC k,i is given by: In our method the next step is to formulate the convolutions and the equations themselves in terms of the discrete values of our functions. Here will sketch the basic idea in some typical scenarios appearing in our equations. Then its application to the concrete unknowns and kernels of the problem is straightforward. If one takes the equations at the required discretized points t (k) j the following typical pattern arises: is the symmetrized kernel to exploit left-right symmetry of the problem for reducing to half the number of variables.
) is intended to modelize the variables in the left-hand side of the equations taken at the discretized points of the transformed variable t and L(u) stands for some nonlinear combination of some unknown function of the equations 13 , then the numerical approximation of the right hand side goes as follows; • First the integration variable is changed from v ′ to t, • then on each subinterval L(u(t)) is approximated by its Chebyshev series, • finally the integration is carried out and the convolution is expressed in terms of the discretized values of L(u(t)).
The final approximation formula takes the form: j )), L(A) denotes the dimension of A and K k,j k ′ ,j ′ is the discretized convolution matrix given by the formula: In this manner a convolution is reduced to a discrete matrix-vector multiplication. The other type of typical convolution is when the integration is taken from zero to 2g. In certain cases the function L(u) has square root behavior close to the branch points 14 . For such functions the truncated Chebyshev series does not give accurate approximation. In these cases not the function L(u) is approximated, but that part of it which remains after the elimination of the square root behavior. Namely, we write L(u) = 4g 2 − u 2L (u), thenL(u) is approximated by a truncated Chebyshev series and finally the approximate discretized form of the corresponding convolution is very similar to (3.8): j ) andK k,j k ′ ,j ′ is the square root factor modified version of (3.9); Here depending on the left hand side of the equation u for some a, or it can denote the sampling points on [0, 2g].
Applying our discretization technique to all unknowns and convolutions of our equations, we can reduce the integral equations to a discrete set of nonlinear algebraic equations. However, the transformation from integral equations to algebraic equations is obviously not exact. The typical error comes from the fact that on each subinterval the Chebyshev series is truncated, so the magnitude of the typical errors in our numerical method is governed by the neglected terms of the Chebyshev series, which can be approximated by the magnitude of the last Chebysev coefficient. In our case this is typically somewhere between 10 −5 and 10 −6 .
The last step of our numerical method is the iterative solution starting from the asymptotic solution.

The iterative solution
Here we will describe two methods to solve our integral equations iteratively. Since our actual equations have very complicated form, we will describe our methods using a model example, which has similar structure to our equations. Let the model equations take the form 15 : where G ab are some kernel matrices, f a are some source terms and y a s are the unknown functions of the problem. The solution of (3.12) is expanded around the asymptotic solution and the equations are formulated in terms of the corrections.
To fix the conventions, the correction functions δy a are defined by: As a consequence: The source term is also expanded around its asymptotic counterpart: f a = f o a + δf a . Then equations (3.12) can be reformulated in terms of the δy a functions as follows: (3.14) To define the iterative method, (3.14) are reformulated so that only O(δy 2 a ) terms remain on the right hand side of the equations. Thus the equations are rewritten in the form: Thus at each step of this iterative method a set of linear integral equations must be solved. Using the discretization method of the previous subsection, the problem reduces to solving a set of linear algebraic equations, which is a straightforward task in numerical mathematics. The very first (0th) iteration starts from the asymptotic solution δy a = 0 and it corresponds to the solution of the linearized equations, which in our case gives the Lüscher-formula for the energy.
This (first) method in a certain range of the coupling constant defined a numerically convergent iteration to solve the equations for the ground state of our D-brane anti-D-brane problem, but beyond a certain value of g the method failed to converge anymore. This is why we worked out a second method, which proved to be much more efficient than the first one. This efficiency is manifested in two facts. First it converges much faster than the previous iterative method, second it gives convergent solutions to our equations even when the solution cannot be accepted as physical one 16 .
This second method can be described simply in words. Instead of defining an iteration as above, we simply take the discretized version of (3.14). We consider it as a set of nonlinear algebraic equations. As a first step we solve the linearized discrete equations (i.e. (3.16) with RHS = 0) and starting from the solution of the linearized equations we solve the discrete nonlinear system by Newton-method 17 . 16 Beyond a certain value of the coupling constant the equations in the form presented in section 2. are not the right ones anymore, they should be corrected by some new source terms and quantization conditions, but even for the "wrong" equations the second method shows numerical convergence, giving unacceptable result. 17 In MATHEMATICA language it can be implemented by FindRoot[...,Method→"Newton"].

Numerical results
In this section we summarize our numerical results. We solved numerically the equations for several integer values of the length parameter L. In this section we concentrate on the states with L ≥ 2. The L = 1 special case is discussed in the next section. For the explanation of the numerical data we will mostly use the L = 2 case as an example, because the critical point of this state is the closest one to zero, so it is enough to work with relatively small values of the coupling constant. This is important from the numerical point of view, since by increasing g the numerical method becomes more and more time consuming. First the parameters of the numerical method is discussed. There are three parameters in the nonlinear integral equations (2.4)-(2.35). The most important one is the coupling constant g, then there are two other parameters which allow us to formulate the equations according to our purposes. The two parameters are p 0 and C, where p 0 is a kind of "truncation index", which tells us the node number starting from which the upper TBA equations are replaced by SU(4) NLIE variables (see figure 1.). The parameter C is a free parameter in the asymptotic solution for the upper SU(4) NLIE variables (C.6-C.20) and it enters the equations such that the asymptotic solution around which the equations are formulated contain this parameter. From this discussion it is obvious that g is a physical parameter which means that the energy depends on it, while the other two parameters p 0 and C correspond to different formulations of the same mathematical problem, so the energy does not depend on them. Thus the choice of these parameters is in our hand and we tried to choose such values for them which allows us numerical convergence in the widest range in g. For example the C = 0 choice is the best for numerical purposes since due to (2.22) a u → −u symmetry arises in the SU(4) HNLIE variables minimizing the number of unknowns in the problem. Tuning p 0 might have two advantages. First, numerical experience shows that for large p 0 the Chebyshev coefficients of the unknowns entering the formula (2.10) for Ω, decay faster, which allows for higher numerical precision. Second also from numerics we learn that with p 0 fixed at certain values of g non physical results are obtained from the numerical solution of the problem. This is a consequence of new local singularities entering the problem, but we still did not take them into account in the equations.  We solved numerically our equations for different values of L and with various values of p 0 and C, and in case the numerical result was physically acceptable for all p 0 and C we tried, it was also independent of these parameters within the numerical errors of the method.
So far we discussed the parameters of the continuous integral equations and their role in the numerical solution. Now we turn to discuss the numerical parameters of the equations. The numerical parameters are artifacts of the numerical method, and they arise mostly from the discretization method described in section 3. We note that there is no cutoff parameter in our numerical method, neither in the integration range nor in the index of Y -functions. Everything is treated in an exact manner, the only source of numerical errors is the discretization of the unknowns and the convolutions. Here we give the most used subinterval vectors of our numerical computations.
On each subinterval we used an l k = 10 order Chebyshev approximation. The subinterval vector A ± of [0, 2g] is given by the empirical formula: where the vectors in components take the form: with v having vector components: • The first element of A Q is 1 2 .
• The length of subintervals in the range 3 < t < B(Q) is approximately 1: ∆t 1.
18 These requirements are based on numerical experiences with the choice l k ≥ 10.
In practice the length of the subintervals are slightly "squeezed" with respect to the conditions above to fill the full [0, B(Q)] properly 19 . Finally, we note that for checking the numerical precision, we also did numerical computations with l k = 12, 14, 16 keeping the subintervals fixed and also with keeping l k = 10, but doubling the number of subinterval points. Before turning to present the numerical results we would like to say a few words about the possible tests of the numerical results. Namely, how one can recognize a wrong result. This is also a very important point of the numerical method, since there are a lot of equations with very complicated kernels and it is easy to make mistakes during writing the code of the numerical solution. There are three basic things that we can check from the numerical results.
The first check is dictated by the energy equation (2.12). It is known that the energy starts at the first wrapping order (i.e. e −L ) and this first order correction is exactly given by the Lüscher formula [4]: with Y o Q (u) given explicitly in (D.2). This quantity can be computed numerically with any digits of precision, so its value is known exactly at any values of g and L. The Lüscher-formula (4.5) corresponds to the linearized version of our equations (2.4)-(2.35), this is why solving the linearized set of equations (which is the first step for the iterative solution) we should reproduce the numerical evaluation of (4.5). This is a nontrivial check on the kernels, on the discretization method and on the equations themselves as well. In addition since this test is quantitative it can tell some information also on the numerical precision of the method 20 .
This test can signal problems on solving the linearized problem. The remaining two tests can signal some discrepancies during the solution of the nonlinear problem.
The second testing condition is that from the numerical solution Y p 0 −2|vw must be real. This sound trivial, but it is not trivial at all. If one takes a look at the equation  BT BA from the exact Lüscher result. This quantity gives some information on the numerical accuracy of the method. Finally the column "number of nodes" tells us the cutoff index of the Lüscher formula, which is necessary to get the Lüscher energy with the precision given by ∆E BT BA . This number is not equal to p 0 in our equations. For the L = 2 state, in case of 0 < g < 1.9 we used p 0 = 4, for 1.9 < g < 2.1 we used p 0 = 8, and in the range 2.1 < g < 2.14 we took p 0 = 12. Finally at g = 2.16 we used p 0 = 26 to get acceptable numerical results. Then beyond this point we could not save our equations from the entrance of new singularities by increasing the value of p 0 with a reasonable O(10) amount. Because of this reason we could not get really close to the supposed critical point. There E BT BA ∼ −1, but we could reach only E BT BA ∼ −0.7 at g = 2.16. Apart from this very embarrassing fact, some important features can be read off from the numerical data. First of all it can be seen that in the range g < 2.16 the energy is very slowly varying function of g, so there is no sign of any divergent behavior. What is more interesting is the behavior of the global constant δc. It is negative and it decreases faster and faster as g is increased. From the definition of δc (2.2) it follows that all Y Q -functions are proportional to its exponent: Y Q ∼ ξ = e δc . The fast decrease of δc indicates that though Y Q has worse and worse large u asymptotic by the increase of g, its global magnitude is actually decreasing. This remark can be understood from the TBA formulation of the energy.
Close to the critical point E BT BA is supposed to be finite [4] E BT BA ∼ 1 − L, but naively the sum in the RHS of (4.7) would diverge due to the large Q terms. Since Y Q is small for large Q, in leading order 21 the log(1 + Y Q ) → Y Q replacement can be done: Diverges close to the critical point +..., (4.8) where Q 0 is an arbitrary index cutoff scale and Y Q = ξỸ Q replacement was applied.
Since ξ is Q-independent all the dangerous Q dependence is still inỸ Q . In (4.8) approaching to the critical point the second sum starts to diverge, and the global multiplicative factor ξ must tend to zero in order to ensure the finiteness of both sides of the equation. Our numerical data seems to support this picture. Namely δc → −∞ as going closer and closer to the critical point.
In [4] from Y -system arguments the large Q behavior of Y Q was also estimated by the formula: where δc is defined after (2.2) in section 2.  the Y -functions of the infinite Y -system. This makes it possible to test numerically the correctness of the large Q estimate (4.9). In case (4.9) holds, it implies that δ lnȳ Q = lnȳ Q − lnȳ o Q tends to zero as 1/Q for large Q. In figure 3. the numerical demonstration of this statement can be seen. The plotted functions are defined by the formula: The plots of figure 3. are based on the numerical computation with p 0 = 26 at g = 2.16. Figure 3. nicely demonstrates the expected 1/Q behavior of the functions δF Q .  This fact shows us that the equations we solved numerically are not the right ones anymore. Something is missing from the equations. Either a special object [44,45] or some other local singularities of the T -and Q-functions of the problem, which enter those strips of the complex plane, which are relevant in the derivation of the HNLIE equations. The numerical data for the L = 3 and L = 4 states are given by table 2 and 3. Also in case of these states the appearance of new singularities obstacled us to get close to the critical point in the framework of the HNLIE technique.

Comments on the L = 1 case
The L = 1 ground state is mysterious, since so far the anomalous dimension of this state could not be determined even for small g either from field theory or from integrability considerations [4]. Here we concentrate on the integrability side. There the boundary Lüscher formula [27,39] diverges for this state [4].    takes the form [4]: This small coupling expression diverges for L = 1, since this point sits exactly on the pole of the ζ -function. As for the origin of this divergence; in (5.1) the individual integrals are convergent, but their sum for Q causes the divergence. In [4] it was argued that also for any larger L the TBA energy formula would diverge beyond a certain critical value of the coupling: g c (L). Assuming that the energy is a monotonously decreasing function of g, which is supported by numerical results, this critical point can be expressed clearly in terms of the energy by the criterion: In [4] this point was interpreted as a turning point where the energy becomes imaginary and as a physical consequence the ground state becomes tachyonic. For the L = 1 state the critical point is right at g = 0 assuming that for small g the energy is also small. Now let us turn our attention to the HNLIE description of the problem detailed in section 2. Here there are no infinite sums and even for L = 1 all the convolutions of the integral equations seem to converge 23 . For the first sight there is no sign of any problem in the HNLIE description and it seems that only the TBA description is inappropriate to treat the L = 1 case. But unfortunately this is not the case.
We can write down the discretized integral equations for the L = 1 case as well, and using the Newton-method, we can solve them for small values of the coupling 24 . We always get some numerical solution for the discretized problem, but it turns out that the Chebyshev coefficients of the unknowns, which correspond to the large u subinterval do not form a decaying series. This phenomenon is a typical sign of some weak (probably logarithmic) large u divergence of the unknowns. If one increases the number of subintervals and sampling points the situation remains the same. The 23 If we assume that large u behavior of the unknown functions, which was used to derive the BTBA equations from discontinuity relations and Y-system. 24 Typically g ∼ 10 −1 .
conclusion is that we can solve the discretized problem, but the solution cannot be interpreted as the discretely approximated version of the continuous solution of our integral equations. In other words the continuous HNLIE equations have no solution for L = 1.
In order to get some analytical insight why the solutions become diverging at large u let us consider the TBA formulation of the problem (p 0 → ∞ in HNLIE). It is known [4] that the TBA energy comes from the coefficient of the most divergent log |u| term in the large u expansion of log Y Q : The E BT BA term originates from the RHS of the TBA equations for log Y Q from the convolution term (2) by exploiting the large u expansion of the kernel: (2) has better large Q ′ behavior than that of dp Q ′ dv , since it behaves like 1/Q ′ . As a consequence contrary to the energy formula, the sum of dressing convolutions is convergent indeed. Thus one might think that for L = 1 the problem emerges, because for the derivation of the energy formula we expanded the sum of dressing convolutions term by term for large u. This is why instead of this usual procedure, we consider the sum of dressing convolutions itself, compute it and then at the end of the computation we take the large u expansion. This procedure is carried out in the small coupling limit. We need the leading order small coupling expression of the dressing kernel in the mirror-mirror channel 25 :

(5.4)
Then the formula, the large u expansion of which accounts for the small coupling expanded energy, is given by: is the leading small coupling expression of (D.2) at L = 1. The second derivative of O(u, Q) can be computed explicitly by simple Fourier space technique. We take the Fourier form of each functions under integration, the convolution is the product of the individual Fourier transforms, the sum for Q ′ can be easily done in Fourier space and at the end of the process everything is transformed back to the u space. In such a manner one gets a bulky, but explicit expression for d 2 du 2 O(u, Q), which we do not present here, only its large u expansion: Integrating twice the large u expansion at small coupling becomes: From (5.7) it is obvious why the naive Lüscher energy formula diverged. Because the leading order large u term is not the expected ∼ log |u|, but ∼ (log u) 2 . This is the key point of the problem, since in this case after this first iteration Y Q acquires an unwanted type of large u term, which makes Y Q divergent for large u: This large u divergence contradicts to what was assumed about the large u behavior of Y Q at the derivation of the integral equations, since it was supposed to decay. In this example we have shown in the small coupling limit, that during the iterative solution of the BTBA equations, log Y Q acquires an extra ∼ (log |u|) 2 behavior at infinity, which made Y Q an exploding function at infinity. This means that the iterative solution of the TBA equations leaves the class of physically acceptable solutions.
One might ask the question, whether it is possible to keep somehow the qualitative large u behaviors that we assumed at the derivation of the equations? Here we sketch a possible idea for small coupling to the L = 1 case. Let us assume that we managed to modify the TBA equations, such that all Y -functions have the large u behavior we want. Since most of the TBA equations reflect the structure of the Y -system functional equations we expect to modify only those equations which are affected by also the discontinuity relations. It follows, that for large Q, the formula for the estimate for Y Q (4.9) remains the same. Now, we assume that for small g the energy is also small and take the simultaneous small g and small energy expansion of the RHS of the TBA energy formula (4.7). In leading order the large Q terms will dominate: whereŶ Q denotes the large Q estimate (4.9) of Y Q ,ξ = ξ g 4E BT BA as a consequence of the u → u/g change of variables and the pole term in E BT BA comes from the pole of the ζ-function. In our HNLIE approach the energy E BT BA and the constant δc are parts of the equations which means that they are not simply expressed by explicit formulas based on the solution of the equations, but must me obtained by solving the set of non-trivially entangled equations. In this sense (5.9) defines an equation for E BT BA for small g. Its leading order solution is: Ifξ > 0 then E BT BA becomes imaginary as it would be expected from string-theory expectations [4]. To decide the sign ofξ, the equation (2.13) has to be analyzed in the context of the small g and E BT BA expansion. It turns out thatξ is positive and O(1) for small g, so according to (5.10) E BT BA is imaginary. Another remarkable fact is that according to (5.10) E BT BA starts at O(g 2 ) instead of the O(g 4 ) prediction of the boundary Lüscher formula (5.1). This might be another explanation why the coefficient of g 4 diverges in the Lüscher formula for the L = 1 case. Finally, we note that in the small g and E BT BA expansion of the L = 1 state, the energy is pure imaginary only at leading order in g, but in higher orders it acquires real part as well.
For the first sight, it might seem that without modifying the equations one immediately gets imaginary energy when going through the critical point. But, the situation is a bit more subtle. There is a hidden tacit modification of the equations. This is realized in (5.9) by the replacement: For the L = 1 case it is an identity for Re(E BT BA ) > 0, but for Re(E BT BA ) < 0 it is not an identity anymore, but a nontrivial analytical continuation in E BT BA .
Such an analytical continuation would require the exact determination of complicated sums of convolutions of the TBA equations as functions of the energy. Since this does not seem to be feasible in practice, we give such an alternative modification of the TBA equations which preserves the infinite sum structure of the equations, but the sums will converge everywhere for Re(E BT BA ) > −L except at the critical value E cr = 1 − L.
The basic idea of the modification comes from the sum representations of the ζ-function. The usual one converges for Re(s) > 1: but there is another representation which converges for Re(s) > 0: Then the original TBA equations are modified through their infinite sums by the replacements: where s E = 4(L + E BT BA ) − 3. Taking into account the large Q behavior of all Y Q functions and all the kernels of the infinite sums of the TBA equations, the new representation will converge for Re(E BT BA ) > −L. This slight modification of the TBA equations might make it possible to go beyond the critical point and get solution of the TBA equations with large u asymptotics being in accordance with the ones used for the derivation of the equations.
The conclusion of this heuristic argument is that to keep the expected 26 qualitative large u behavior a nontrivial modification of the TBA equations must be carried out, which might lead to complex energies.

Summary and conclusions
In this paper we studied the ground state energy of a pair of open strings stretching between a coincident D3-brane anti-D3-brane pair in S 5 of AdS 5 × S 5 . The main motivation for the study is that string-theory predicts that the ground state of such a configuration becomes tachyonic for large values of the 't Hooft coupling [4].
In [4] it was shown that the usual integrability based BTBA approach always give real energies for the ground state and it breaks down at latest when the energy gets close to the critical value: During the numerical solution of the HNLIE equations the usual iterative methods failed to converge, this is why we worked out two numerical methods to reach convergence. The most effective one is, if one transforms the integral equations into discrete nonlinear algebraic equations and solves them by Newton-method. The power of this method is demonstrated by the fact that it gives convergent results even if the numerical solution is not physically acceptable.
Unfortunately, in our numerical studies we could not get very close to the critical point, because new singularities entered the HNLIE equations taking into account of which would have required an enormous amount of additional work. Nevertheless, in the range where we could get physically acceptable results, the precision of the HNLIE data were higher than those of BTBA and the HNLIE approach could give a deeper understanding of the problem. For the ground state of the L = 1 state the critical point is right at g = 0 and neither perturbative field theory computations nor the boundary Lüscher formula could provide a finite quantitative answer to the anomalous dimension. Even in this special case the numerical solution of the HNLIE equations was possible. The results showed that without an appropriate modification of the equations, they cannot give physically acceptable results. In this case, it means that the solution of the dicretized problem cannot be considered as a discretized solution of the continuous nonlinear integral equations. Moreover the large rapidity behavior of the numerical solution is incompatible with the one assumed for the derivation of the equations. This phenomenon is analytically analyzed in the framework of BTBA and an idea is sketched to preserve the expected large rapidity behavior of the unknowns. This method is based on an appropriate modification of the TBA equations which would lead to complex energies beyond the critical point.
Hopefully the L = 1 case at g = 0 could be treated analytically in the framework of the quantum spectral curve method [47,48,14], solving the mystery of this state in the context of integrability.

A Notations, kinematical variables, kernels
Throughout the paper we use the basic notations and TBA kernels of ref. [41], which we summarize below. For any function f , we denote f ± (u) = f (u ± i g ) and in general f [±a] (u) = f (u ± i g a), where the relation between g and the 't Hooft coupling λ is given by λ = 4π 2 g 2 . Most of the kernels and also the asymptotic solutions of the HNLIE-system are expressed in terms of the function x(u): which maps the u-plane with cuts [−∞, −2] ∪ [2, ∞] onto the physical region of the mirror theory, and in terms of the function x s (u) which maps the u-plane with the cut [−2, 2] onto the physical region of the string theory. Both functions satisfy the identity x(u) + 1 x(u) = u and they are related by the x(u) = x s (u), and x(u) = 1/x s (u) relations on the lower and upper half planes of the complex plane respectively.
The momentump Q and the energyẼ Q of a mirror Q-particle are expressed in terms of x(u) as follows: Two different types of convolutions appear in the HNLIE equations. These are: The kernels and kernel vectors entering the HNLIE equations can be grouped into two sets. The kernels from the first group are functions of only the difference of the rapidities, thus actually they depend on a single variable. The other group of kernels composed of those, which are not of difference type.
We start with listing the kernels depending on a single variable: The fundamental building block of kernels which are not of difference type is: Using the kernels K(u, v) and K Q (u − v) it is possible to define a series of kernels which are connected to the fermionic Y ± -functions. They are: and Further important kernels entering the Y ± related TBA-type equations are defined as follows: . (A.10) The kernels entering the right hand sides of the equation (2.9) for Y 1 are , (A.11) and the dressing-phase related kernel K QM sl(2) (u, v), which is built from the sl(2) Smatrix of the model [49]. It is of the form where Σ QM is the improved dressing factor [50]. The corresponding sl(2) and dressing kernels are defined in the usual way Explicit expressions for the improved dressing factors Σ QM (u, v) can be found in section 6 of ref. [50]. Here for our numerical computations we used the single integral representation given in [21].
Finally we mention that along the lines of [42] in the derivation of the formula (2.10) for Ω(K Q ), it was exploited that all the necessary kernels: K Q , K Qy , K Q1 xv , s ⋆ K Q−1,1 vwx , K y1 , K Q1 sl(2) satisfy the identity: (A.14)

B Kernel matrices of the vertical HNLIE part
In this appendix the kernel matrices appearing in the upper HNLIE part of our equations (2.18,2.19) are presented. Here the kernel matrices are different compared to those published in [38]. The difference comes simply from a reformulation the equations in the language of new unknown functions. In [38] the unknowns are 6 b-type functions:

C Asymptotic solutions of the vertical HNLIE
In this section along the lines of [38] the asymptotic solutions of the upper SU(4) NLIE variables are presented . In the asymptotic limit the T-hook of AdS/CFT splits into two SU(2|2) fat-hooks. The basic building blocks of the asymptotic solution are the nine Q-functions corresponding to the left and right SU(2|2) fathooks. Due to the left-right symmetry of the Y -system it is enough to give the right Q-functions. They can be derived from the asymptotic solution of the Y-functions given in [4]. They take the form: w o (u) = − i Λ 2 e −π g u 4 ((g u) 2 + w c ), y o (u) = i Λ 2 e −π g u 4 ((g u) 2 + w c − i C), (C.5) 28 Here for correspondence we use the same letters for the names of different unknowns as in [38].
where w c and C are arbitrary constants. Using the building blocks listed above, the system can be obtained from the Bäcklund functions above by appropriately shifting their arguments:  where * denotes complex conjugation. In our numerical studies we mostly use the C = 0 asymptotic solution to setup the equations to solve. In this case the exact equations guarantee the fulfillment of (2.22), which reduces to 6 the number of independent complex functions of the upper NLIE part.