The spectrum of tachyons in AdS/CFT

We analyze the spectrum of open strings stretched between a D-brane and an anti-D-brane in planar AdS/CFT using various tools. We focus on open strings ending on two giant gravitons with different orientation in $AdS_5 \times S^5$ and study the spectrum of string excitations using the following approaches: open spin-chain, boundary asymptotic Bethe ansatz and boundary thermodynamic Bethe ansatz (BTBA). We find agreement between a perturbative high order diagrammatic calculation in ${\cal N}=4$ SYM and the leading finite-size boundary Luscher correction. We study the ground state energy of the system at finite coupling by deriving and numerically solving a set of BTBA equations. While the numerics give reasonable results at small coupling, they break down at finite coupling when the total energy of the string gets close to zero, possibly indicating that the state turns tachyonic. The location of the breakdown is also predicted analytically.

the strings connecting the two [7]. Both are examples of open superstring systems undergoing the wrong GSO projection.
In this paper we study the coincident D-D system within the AdS/CFT correspondence. In flat space, when the two D-branes are coincident, the ground state of the open string connecting them is tachyonic with a masssquared −1/(2α ′ ). If the D-branes are not coincident this mass squared is increased and beyond a distance of order the string length all states become massive.
Unstable D-brane systems have been studied within the AdS/CFT correspondence initially in [8]. In certain cases it is possible to match instabilities in the field theory to those in string theory. Generically, these systems are not amenable to perturbative calculations in either, let alone both, weak and strong coupling. In special circumstances it has been possible to take a scaling limit to get a match between weak and strong coupling [9,10]. Here we study an unstable system beyond such a limit. 1 Our study relies on the integrability of the AdS5 × S 5 superstring, a property which is conjectured to hold beyond the classical string limit. Integrability has led to a great understanding of the spectrum of closed string states in AdS5 × S 5 , as well as certain open-string sectors which are conjectured to be integrable. These correspond to strings ending on different types of D-branes [11,12,13,14,15] as well as macroscopic open strings extending to the boundary of space and representing Wilson loops in the dual 4d gauge theory [16,17,18]. The most studied case is that of a "giant graviton", which is a D3-brane carrying N units of angular momentum on S 5 [19].
Integrability of the giant graviton systems can be seen from both sides of the AdS/CFT correspondence.
From the gauge theory side, the dilatation matrix -calculated perturbatively -coincides with an integrable open spin chain Hamiltonian [12,13]. From the string theory side, integrability is a consequence of the fact that the classical two-dimensional sigma model with boundary admits a Lax pair formulation, which leads to an infinite number of conserved charges [20,21].
Instead of a single D-brane we consider here a pair of coincident D-branes with arbitrary orientation. When the two orientations are identical this is a BPS system, and when opposite, this is a D-D system. Thanks to integrability, we can compute the asymptotic spectrum of open strings on these D-branes by solving the Bethe Ansatz equations with boundaries. To be more specific, let us choose the reference ground state of the Bethe Ansatz as Z L . There are two important orientations of the D-brane, one carrying the angular momentum on S 5 in the same direction ("Z = 0 giant graviton"), or the other in a perpendicular direction ("Y = 0 giant graviton") [13]. The names reflect the fact that the world-volume of the D3-brane is embedded as S 3 ⊂ S 5 ⊂ C 3 , where C 3 is parameterized by complex X, Y, Z coordinates satisfying |X| 2 + |Y | 2 + |Z| 2 = 1. In our problem, one brane satisfies Y = 0 while the other an arbitrary linear equation involving Y ,Ȳ , X andX, which we callŶ = 0 witĥ Y = Y cos θ1 cos θ2 − X cos θ1 sin θ2 +X sin θ1 cos θ2 +Ȳ sin θ1 sin θ2 . (1.1) We will mainly concentrate on the D-D system, which corresponds to θ1 = θ2 = π/2, but many of the calculations can be generalized to arbitrary angles.
In the next section we discuss the gauge theory dual of these operators. The Y = 0 giant graviton is a determinant operator made of N of the Y scalar fields [22]. An open string attached to it is obtained by replacing one of the Y fields with an adjoint-valued word made of other fields (and covariant derivatives). The system studied here should involve two determinants connected by a pair of adjoint valued words with mixed indices. A single adjoint valued word can replace one of the letters in one determinant, but not both, which is why two words are required. The dual statement in string theory is that a compact D-brane cannot support a single open string, due to the Gauss law constraint and must have an even number of strings (with appropriate orientations) attached. In the planar approximation the two open strings should not interact, which we verify in the gauge theory calculation in the next section. Hence we can consider the spectral problem independently for each of the insertions/open strings.
The exact gauge theory description of the D-D system is not known and it requires solving a mixing problem which is quite complicated, because the operator consists of more than N fields. At tree level, an orthogonal basis of gauge-invariant scalar operators is constructed by the Brauer algebra [23] or the restricted Schur polynomial [24] at any N . At loop level, little is known about how to find dilatation eigenstates in the D-D system using these bases [25]. Nevertheless, the mixing problem of our interest seems to simplify at large N . We expect that the D-D system with or without open strings has the gauge theory dual closely resembling the double determinant.
The mass of the (potentially tachyonic) open-string state should correspond to the dimension of the local operator, or more precisely, the contribution to the dimension from the insertion of the word into the determinant operators as discussed in Section 2. In the case of the Y = 0 brane, the insertion of the word Z L corresponding to the ground state of the open string gives a protected operator. The system with Y ,Ȳ and Z L is not protected and we expect that the ground state energy is lifted by 'wrapping type' graphs which involve the interaction between the Y andȲ fields at the two boundaries of the word. We identify a set of such graphs at order λ 2L in perturbation theory, which we conjecture to be the first ones to contribute to the anomalous dimension of these operators. In fact, the leading non-vanishing wrapping correction coming from the integrability formulation derived in Section 3 is exactly of order λ 2L and equal to the UV divergences of the integrals that arise from these graphs. The UV-divergences of these integrals were recently proven to agree with our conjecture [26].
In the integrable description we identify how the string excitations scatter off from the D-branes. Combining in flat space is m 2 = −1/(2α ′ ), which translates to − √ λ/2 in units of the AdS curvature radius. 2 In the case of arbitrary angles θ1, θ2 [29] the expression becomes 3 The dimension of the operator inserted in the determinants is dominated by the charge L at weak coupling, but at strong coupling it should asymptote to m ∝ i 4 √ λ. We therefore expect the dimension to turn imaginary at a finite value of the coupling. To probe this transition we employ boundary thermodynamic Bethe ansatz equations (BTBA).
BTBA was derived first in [30] for models with diagonal S-matrices. If the S-matrix is non-diagonal it is difficult to construct BTBA explicitly by applying the methods of [30]. In specific cases, however, it is possible to overcome the appearing technical problems and derive a BTBA with non-diagonal S-matrices, as was done for example in [17,18]. As our case is more complicated, the approach we take here is to use the Y-system equations together with their analytic properties to derive the BTBA, as was done in [31]. In Section 4 we apply this method (following [32]) to derive a set of BTBA for the ground state of the θ1 = θ2 = π/2 case.
We develop numerical algorithms to solve these equations and evaluate the anomalous dimensions of ground states with different values of L at finite coupling. In all cases we find that the anomalous dimension is a monotonously decreasing function of the coupling. However, when the BTBA energy becomes comparable to 2 If we compare to the mass-squared of the Konishi operator 4 √ λ, [27,28] which matches the first excited closed string state in flat space, a factor of 4 arises from replacing closed strings by open strings, and another factor of 2 from taking the wrong GSO projection. 3 Note that for θ 1 = 0 the rotation in (1.1) mixes only Y and X, so the ground state Z L is still BPS, hence the mass is zero. Likewise for θ 2 = 0.
1 − L, namely when the total energy of the open string gets close to zero, it becomes extremely difficult to obtain the precise value of the energy from BTBA solutions. As a result, the evolution of the energy cannot be traced further toward strong coupling.
Such a pathological behavior can arise for states with negative anomalous dimension. A novel lower bound for the BTBA energy is derived analytically, and the violation of this bound makes the BTBA solution inconsistent.
We expect this breakdown to signal the transition of the states from massive at weak coupling to tachyonic beyond the critical value of the coupling. Beyond this singular point another formalism must be employed to find a continuation of the BTBA equations, whose details are beyond the scope of this paper.
2 The YȲ brane system in gauge theory The Y = 0 giant graviton is described in the gauge theory by a determinant operator [22] where ai and bi are color indices and ǫ is a product of two regular epsilon tensors ǫ a An open string ending on the giant graviton is described by replacing one Y with an adjoint valued local The simplest insertion is the vacuum W = Z L .
One can consider also two giant gravitons by taking the combination OY OY and likewise add an open string attached to one or to the other. But with two giant gravitons we can also consider strings stretched between the two D-branes. Having a single such string is impossible, though. The endpoint of a string serves as a source of charge on the D-brane world-volume, which is compact, and there must be another charge source with the opposite sign. We therefore will consider the case of a pair of open strings with opposite orientation connecting the two D-branes.
The gauge theory description of this system is the double-determinant operator with all fields at a single so one Y was removed from each determinant and then the two words W and V are inserted with the indices crossed.
However, we are not interested here in the case with two identical D-branes. That configuration is BPS and the spectrum of open strings also includes BPS states, which belong to the multiplet of the non-abelian gauge fields on the pair of D-branes. We want to study instead the spectrum of strings stretching between a D-brane and an anti D-brane.
The anti D-brane can be realized by replacing all the Y byȲ , and we shall call it aȲ = 0 giant graviton, as opposed to the original Y = 0 giant graviton. If we parameterize the S 5 part of the target space by three complex coordinates Y , X and Z subject to |X| 2 + |Y | 2 + |Z| 2 = 1, then the Y = 0 giant graviton wraps an S 3 given by Y = 0. TheȲ = 0 brane will wrap the same S 3 but with the opposite orientation.
With this we have constructed a two parameter family of pairs of D-branes interpolating between the pair of identical D-branes and the D-D systems.
We expect the planar dilatation operator to act on this complicated operator as the sum of two independent integrable open spin-chain Hamiltonians, one acting on each of the two words: 4 We will verify this splitting now at the one loop level in perturbation theory, and then assume it holds in general. This allows us to study the spectrum of each of these open strings independently by application of different integrability tools.

Integrable spin-chain
To calculate the conformal dimension of an operator, we consider the two point function between two similar operators and find the mixing matrix. We therefore study the pair of operators O W,V Y,Ȳ (0) and OV . It is useful to separate the calculation according to how many of the fields Y andȲ from the determinant operators interact with the W and V insertions. In the case when none do, we perform free field contractions between the Y andȲ of the two operators. Using This is a non-local trace, which is not gauge invariant, but this is due to the fact that 2(N − 1) contractions were already done in a specific gauge. The entire correlator is of course gauge invariant.
In the planar approximation the expectation value in (2.7) factorizes and at tree level we find that W andW ′ are conjugate operators, as are V andV ′ . Each trace gives an extra factor of N . This statement holds as long as the last letter in W and the first one in V are orthogonal to Y and the other ends of the words are orthogonal tō Y . Otherwise there would be extra planar tree-level contractions beyond (2.6) which will mix these states with operators made of a sub-determinant and a single trace operator.
There are also interacting graphs contributing to Tr W(0)W ′ (x) , which by construction do not know about the rest of the determinant operators. These will give the bulk part of the spin-chain Hamiltonian. At one-loop level in the so(6) sector this is the same as the usual Minahan-Zarembo Hamiltonian [34].
We should consider separately the boundary interactions, where the beginning and end of W and V interact with Y andȲ from the rest of the determinant. The first boundary interactions involve just one pair of Y and Y . The trace structure arising from contracting all the determinant fields except for one Y (0) and oneȲ (x) is given by (B.2) in Appendix B.
The one-loop interactions arising from these graphs are very similar to the boundary interaction for the giant graviton open spin-chain [12]. Again the interaction between W and V completely factorizes at large N and for each word one finds the usual one-loop boundary interaction. The only modification is that the first letter is 4 The dilatation operator actually mixes also the structure of the Y andȲ fields, which is an indication that this is not the exact state. We assumed that the correct state without insertions to be protected at large N . Otherwise the first term of (2.5) is a certain function of λ. Figure 1: Graphs for insertion of W = Z (the case of L = 1). The fields are ordered according to the trace structure in (2.9).
projected on states orthogonal toȲ and the last letter should be orthogonal to Y . The one-loop Hamiltonian acting on the word W is Here I, P and K are respectively the usual identity, permutation and trace operators on the spin-chain [34]. Q φ l is a projector whose kernel are all words with the field φ at location l. The Hamiltonian acting on the word V can be obtained by exchanging Y ↔Ȳ .
Although we do not derive here the explicit Hamiltonian at higher loop order, one can still write down the Bethe ansatz, as we do in the next section.

Wrapping corrections
One can proceed this way to higher order boundary interactions, but we would like to study the first wrapping corrections, where interactions are communicated between the two boundaries and the energy of the Z L ground state is lifted.
Wrapping graphs come from the interaction of the word W with Y on one side andȲ on the other. The leading wrapping corrections will arise by choosing one Y (0) and oneȲ (0) from each operator and requiring that they all interact with W. We analyze this in Appendix B and find that we should include connected graphs contributing to To be more specific, consider the ground state W = Z L . Since it shares some supercharges with each of the determinants, the interaction with only one boundary will not give rise to an anomalous dimension. This is identical to the state attached to two Y determinants. Only when the interaction involves both a Y andȲ determinant will the ground state energy be lifted. 5 To capture this we can consider the difference between the two cases Figure 2: Graphs for insertion of W = Z 2 . The fields are ordered according to the trace structure in (2.9). In the case of L = 1 there seem to be two types of relevant 2-loop graphs, depicted in Figure 1. The graph (a) exists also in the case where both boundaries are on the Y = 0 brane (i.e., withȲ (0) → Y (0) and Y (x) →Ȳ (x)).
We expect wrapping effects to start at order 2L in perturbation theory. 6 The generalization of the graph in there are many other graphs of the same order as this graph. For example when the box is replaced with a fermionic hexagon. This graph generalizes to arbitrary L and is of order L + 2, so again it should be the same as for the BPS vacuum and cancel against other graphs.
For general L we conjecture therefore that the first wrapping correction comes from the graph in Figure 3.
It gives rise to a UV divergent loop integral depicted on the right of Figure 3, where the ellipses correspond to repeating the structure to generate the total number of 2L loops in the integrals.
Based on explicit data for L = 1, 2, 3 and our conjecture for larger L (from the results of the next section), these integrals were recently shown in [26] to have the same divergence structure as the zig-zag integrals [35,36].
We present the map of the above integrals to the zig-zag integrals in some more detail in appendix C. The results for the overall UV divergences (denoted by calligraphic I) of the integral I2 arising in the L = 1 case from the 6 Note though that in the quark-antiquark case the first interaction happen at "half" wrapping order (L + 1), which can be attributed to the finite density of zero momentum single particle states in the mirror BTBA formulation as encoded by the poles of the Y-functions [17,18]. diagrams in Figure 1 and the integrals I2L depicted on the right of Figure 3, regularized in D = 4 − 2ε (and with the coupling restored) are where our conjectured result I2L was proven recently (cf., Eqn (3) of [26]). Note that for L ≥ 2 the integrals are free of subdivergences, which is indicated by the absence of all higher order poles in ε. Assuming these are indeed the first graphs that contribute, we conclude that the anomalous dimension of the ground state is At L = 1, however, the integral contains a one-loop subdivergence which leads to an inconsistency here: it enters as a simple ε-pole in the anomalous dimension, which has to be finite and independent of ε. In this setup, i.e. for a gauge invariant composite operator in a theory with unrenormalized coupling, such a oneloop subdivergence at two loops can only be cancelled either by further two-loop diagrams or by a one-loop counter term, both of which are associated with the renormalization of this operator. We conjecture that the approximation of the state with L = 1 by (2.4) is inappropriate and that the correct state will be renormalized at one-loop order.
The fact that short states are subtle and can lead to divergences was seen in the past in the integrability-based descriptions [37,38,39,40]. In certain cases (deformed or orbifold setups) a possible resolution of this issue in the field theory was given in [41,42], though the effect observed there has no net effect for the N = 4 SYM theory. We note that the integrability calculation in the next section also gives a divergence for the L = 1 state -so whatever effect lifts this divergence (presumably by a one-loop anomalous dimension) should somehow also alter the integrability description of this state. We leave it to the future to resolve this issue.

General angle
As mentioned above, the two D-branes do not have to be coincident, but can be at arbitrary angles on S 5 , which on the gauge theory side amounts to replacingȲ with an arbitrary linear combination of Y , X,X andȲ which we denoted byŶ in (1.1). It is easy to see that the wrapping graphs we calculated will see only theȲ factor in Y and the result of the first wrapping effect will be multiplied by sin 2 θ1 sin 2 θ2.

Integrable description of the YȲ brane system
In this section we formulate the integrable description of the YȲ brane system. In the integrable formulation we characterize and solve the system in terms of the scattering data of the particle-like excitations of the strings.
The finite-volume energy spectrum of the particles corresponds to the sought-for anomalous dimensions on the gauge theory side.
D-branes provide boundary conditions for open strings, which translate into reflection amplitudes for the elementary particle-like excitations. The bulk scattering matrix of these excitations supplemented by the reflection factors define the theory and enable one to calculate both the asymptotic large-volume energy spectrum and all finite-size effects.
In the integrable description we assume the quantum integrability of the model, and solve the theory in semiinfinite geometry by determining the scattering (reflection) data. Integrability forces the multiparticle reflection process to factorize into pairwise scatterings and individual reflections. Integrable boundary conditions in this point of view can be classified by finding all one-particle reflection matrices, which are compatible with the given bulk S-matrix and residual symmetries.
When two boundaries exist their relative orientation is also important, which is used to break the supersymmetry of the vacuum state of the Y = 0 brane studied in [43]. The new vacuum state acquires a nontrivial anomalous dimension from finite-size effects. Our notation in this section is summarized in Appendix A.

Reflection matrices inŶ
is broken by the presence of the Y = 0 brane. The reflection factor compatible with the unbroken symmetry, which satisfies the boundary Yang-Baxter equation has the following factorized form [13,44] and σ(p, −p) is the BES dressing factor [45]. This reflection factor can be extended for bound-states both in the string/mirror theories belonging to the totally symmetric/anti-symmetric atypical representations of su(2|2) ⊕ su(2|2), respectively [46,47,48]. The totally anti-symmetric representation describing the bound-states of a fundamental particles in the mirror theory has a diagonal reflection factor and its scalar factor is obtained by fusion.
Although the presence of the D-brane breaks the rotational symmetry, this symmetry does not completely disappear from the system. Acting with such a transformation (3.1) will rotate the D-brane itself and acts on the reflection factors in the following way: where O acts as in Eq. (3.1) in the (1, 2) space and as identity in the (3, 4) space. We introduce two rotation angles θ1 and θ2 for dotted and undotted indices. The reflection factor has to satisfy unitarity, boundary crossing unitarity and the boundary Yang-Baxter equations to maintain integrability. The reflection factor (3.5) solves these constraints, as the rotation by O is part of the bulk symmetry which commutes with the S-matrix.
If we choose θ1 = θ2 = π 2 we obtain the reflection factor of theȲ = 0 system. This is the anti-brane counterpart of the Y = 0 brane, and has a reflection factor in which the two labels (1, 2) are exchanged This is nothing but the charge conjugated reflection factor: This picture extends to the reflection factors of the bound states, too: they can be simply obtained by exchanging the labels (1, 2).

YȲ system in large volume
In the following we analyze a two-boundary system in finite volume, namely in the strip geometry. We place R − Y on the right boundary but on the left boundary. We are interested in the asymptotic spectrum of multiparticle states and an exact description of the ground state. Both problems can be attacked via double-row transfer matrices and Y-system. The energy of a multiparticle state gives a half of the total anomalous dimension in gauge theory (2.5), and the other half is obtained in an analogous way.
The boundary Bethe-Yang equations (also called boundary asymptotic Bethe Ansatz equations) are valid for large size L, the R-charge of the inserted word. They determine the momenta, {pi}, of a multiparticle state by the periodicity of the wave function as follows. Pick up any particle (with momentum p k say), scatter through the others to the right, reflect back from the right boundary, scatter with momentum −p k on all particles to the left, reflect back from the left boundary and scatter back to its original position; we have to arrive at the same state multiplied by e −2ip k L . During the scattering processes the labels of the multiparticle state are mixed up, and the problem is to find eigenstates of the mixing matrix. This problem is solved by introducing and diagonalizing the double-row transfer matrix. As the diagonalization problem factorizes between the two su(2|2) factors we focus on one copy only, and define this double-row transfer matrix by where we used the non-graded S-matrix S = S0 S ⊗ S, in the su(2) normalization (S 11 11 = 1). On the left boundary R + Y (p) is introduced in (3.9), which is different from the standard definition [49]. Nevertheless the two are equivalent up to an overall normalization [50,43]. Our choice ensures that a "test" particle is brought around the two boundaries in the above sense; if we specify the test particle momenta as p = p k we obtain the k th particle's mixing matrix.
Let us diagonalize the double-row transfer matrix for YȲ states in the su(2) sector. We start by analyzing the ground state of the first level |1, 1, . . . , 1 in algebraic Bethe Ansatz, and denote its eigenvalue by T (p|{pi}) or T for short. It describes an M -particle state in the su(2) sector. Interestingly, similarly to the Y Y system [43], the eigenvalue can be expressed in terms of the diagonal elements: which explicitly read as: The ρi (i = 1, . . . , 4) can be calculated following the considerations in [43]. The result turns out to be Clearly x − ↔ x + exchanges 1 ↔ 2 and 3 ↔ 4, and thus acts as conjugation, under which the YȲ ground state is invariant. The momenta of the state are determined from the following asymptotic boundary Bethe-Yang where we introduced the proper normalization factor which is determined by the boundary Bethe-Yang equations for one-particle states (D.1).
The extension from the diagonal sector to the full sector can be easily done at the level of the generating functional, which we now introduce. The eigenvalues of the double-row transfer matrix, in which mirror boundstate test particles of charge a are scattered and reflected through the multiparticle state, are generated as in the su(2) grading, where D is the shift operator (A.3). Technically it is simpler to renormalize the generating functional and the transfer matrices as The relation between the normalizations is of the fusion type: Explicit calculation gives all the antisymmetric transfer matrix eigenvalues where as in (A.4), B1R3 and R1B3 represent type 1 Bethe roots denoted by yj, and Q2 represents type 2 Bethe roots denoted byμ l . Regularity of the transfer matrix at the roots gives the boundary Bethe-Yang equations.
Type 1 roots are specified as x + (p) = yj , type 2 roots when u(p) =μ l , finally type 3 roots when x − (p) = y −1 j . The corresponding Bethe equations read as The Bethe Ansatz equation which determine the momenta are

Asymptotic Y-system for the vacuum state
From now on we focus only on the unprotected vacuum state. As M = 0, we have R = B = 1, and the expression (3.19) for the eigenvalues of the transfer matrices simplifies considerably: They constitute part of a solution of the T-system on the su(2|2)-hook. For completeness and later applications we provide here the full solution of the su(2|2) T-system. The transfer matrix eigenvalues in the symmetric representations are generated via the inverse of (3.17): which results in (3.26) The T-functions on the boundary of the su(2|2)-hook are The asymptotic Y-functions are defined from the T-functions as for s > 0. For s = 0, (in a similar analysis for the su(2) sector) Y1,0(p k ) = −1 should provide the boundary Bethe-Yang equations (3.14). This allows us to restore the correct normalization: The normalization of the bound state transfer matrix eigenvalues follow from the bootstrap:

Lüscher correction for the vacuum state
In the following we use the asymptotic Y-functions to calculate the leading finite-size -so called Lüschercorrection for the vacuum state. For this we analytically continue Ya,0 in u to the mirror plane: where we denote the analytically-continued variables x [±a] by z ± , which can be parametrized by the mirror momenta q as: We can compare the Ya,0 functions with the integrand of the vacuum Lüscher correction [48,50] calculated directly from the reflection matrices: As charge conjugation exchanges theȲ = 0 boundary with the Y = 0 boundary, we simply square the analytically-continued bound state reflection factor (3.4) and perform the trace. This gives for the matrix part The prefactor was already calculated in [50] R0 It is now easy to evaluate the finite-size correction in the weak coupling limit. At leading order in g 2 we find the following correction for the vacuum: which agrees precisely with the gauge theory result (2.12) for L ≥ 2, and diverges at L = 1.

Generic angle, theŶ = 0 brane
Here we analyze the system with generic angles. We keep R − Y on the right boundary but place R + θ (p) = R − θ (−p) on the left boundary. The reflection factor in the totally antisymmetric representation can be dressed as: 7 (3.37)

Lüscher correction
In order to calculate the Lüscher correction for the ground state energy, we start from the expression in Eq.
(3.33). Only the matrix part is deformed by the angle: which shows that we simply have to include an additional sin 2 θ factor compared to the YȲ system for each su(2|2) wing. The resulting Y-functions are which at leading order leads to the wrapping correction This is precisely what we expect from gauge theory calculations for the Y -Ŷ brane system. It is (3.36) multiplied by the square of the respective angular dependence in (1.1).
The generating functional for the vacuum in case of a generic angle is analyzed in Appendix E.

The YȲ ground state BTBA
In this section we derive the ground state BTBA equations for the YȲ system and analyze them numerically.
TBA equations in the presence of boundaries can be formulated in the same way as in the periodic case, provided that the S-matrix and the boundary reflection amplitudes are diagonal [30]. BTBA follows from the mirror trick, which equates the open string worldsheet partition function in the string region with the closed string transition amplitude between boundary states in the mirror region. In the mirror picture, the boundary state projects the intermediate states to those consisting of an even number of particles with the opposite momentum. As a result, a Y-function in the ground state BTBA is the ratio of the density of particle pairs to that of hole pairs.
When the S-matrix is non-diagonal, it becomes very difficult to compute the source term in BTBA explicitly, which comes from the overlap between a boundary state and the bulk state written in terms of the density of Bethe roots and holes. Thus, a simple alternative approach is called for. Recall that the periodic TBA can also be derived by integrating the Y-system assuming appropriate discontinuity relations and analyticity of Y-functions [51,32]. In this section we apply this method to derive a set of BTBA equations, and solve them numerically.

Boundary TBA from Y-system and discontinuity relations
The derivation of the equations goes along the lines of ref. [32] relying on the following assumptions: 8 • There exist TBA-type integral equations governing the spectrum of the YȲ system.
• The Y-functions of the BTBA equations satisfy the Y-system functional equations [52] of AdS/CFT.
• The Y-functions satisfy the discontinuity relations of ref. [51], too.
• The Y-functions are real functions. They are meromorphic in the vicinity of the real axis away from cuts prescribed by the discontinuity relations.
• The ground state Y-functions are parity even and left-right symmetric.
• The Y-functions are smooth deformations of their asymptotic limit, so qualitative information on the location of their point-like singularities can be borrowed from the asymptotic solution.
• The massive Y-functions decay at large rapidity at least as 1/u, while the large u behavior of the other Y-functions is the same as that of the asymptotic solution, namely in the u → ∞ limit they tend to stateand coupling-independent constants. In this section the notations of refs. [53,32] are used so that their results could be referred directly. All kernels and source functions of the subsequent BTBA equations can be found in appendix A of ref. [32]. • YQ, Y±, Y m|vw , and Y 2m−1|w have double zero at u = 0 for Q, m = 1, 2, ....
• Y1 and Y 1|w have simple poles at u = ±i/g.
The derivation of the BTBA equations goes along the lines of ref. [32]. Most of the equations can be derived straightforwardly from the Y-system equations by taking into account the residue contributions of the local singularities listed above. The two subtle equations are the discontinuity functions of log Y− and log Y1 denoted by J and ∆, respectively [32]. The derivation of equations for these quantities requires the usage of discontinuity relations of [51]. Since asymptotically e J = 1 for the ground state it is assumed that it has no local singularities on the whole complex plane. So it follows that Y − Y + is given by (5.30) of [32] by taking the set {uj } = ∅ or equivalently R p/m → 1 and B p/m → 1.
The computation of ∆ goes along the lines of section 6. of ref. [32] taking into account the different singularity structure and asymptotic behavior of the Y-functions. Here we introduce the notations: The discontinuity function ∆ satisfies the following equation: where ∆ red is given by 10 and the source term W is given by the integral representation: where [±γ] means that the integration runs along the lines v ± iγ from −∞ to ∞, with γ being a small positive number. Starting from the discontinuity function (4.2) and applying the simplification techniques described in sections 7 and 8 of [32], the hybrid BTBA equations for the massive Y-functions can be derived. Carrying out the whole process the following BTBA equations were derived: The symbols ⋆,⋆ denote the convolutions defined in (A.4) of [32]. Equations (4.10) and (4.11) 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 (4.10) and (4.11) by the formula: For the massive Y-functions we present the hybrid form of the BTBA equations. where the source term fQ is given by: withW The energy of the YȲ ground state after subtraction of the bare dimension is given by the formula  Before closing the subsection we argue that due to the constraints imposed by the Y-system equations, the L-R symmetry and parity, the local singularities of the Y-functions located at the positions {0, ± i g } do not receive any wrapping corrections. As an example we show that the double zero of Y1 located at the origin remains fixed at any value of the coupling constant. As a first step let us invoke the Y-system equations At the level of the asymptotic solution, the double zero of Y1 at the origin is related to a simple zero of 1 − Y−(v) at v = ± i g . Suppose that wrapping corrections change the quantization condition as Since Y−(v) is parity even, this equation should be valid after the parity transformation v → −v, which implies Y1(v) = 0 2 at v = +iδ. Now if we take the asymptotic limit, Y1 possesses a quartic zero instead of a double zero at the origin, which is a contradiction. In other words, the origin v = 0 is a special point where the parity transformation acts trivially, thus a double zero is allowed.

Lower bounds for TBA energy in AdS 5 × S 5
In this subsection we show that the energy which can be computed from the BTBA equations of the YȲ ground state is bounded from below. This means that the BTBA equations can describe the model as long as the energy is real and remains above the lower bound.
To determine a lower bound, the energy formula (4.16) must be studied. In order for the energy to be finite, the individual integrals should stay finite; and having evaluated the integrals, the remaining sum must also converge.
First consider the case of individual integrals. To see their convergence the large u behavior of the integrand must be analyzed. Since dp Q du ∼ g for large u and at any Q, YQ(u) must decay faster than 1/u at infinity. This simple remark constrains the range of the BTBA energy, because the large rapidity behavior of YQ is governed by the exact energy through the formula: 11 log YQ(u) = − (4L + 4EBTBA) log |u| + O(1), (|u| ≫ 1). (4.20) 11 The formula (4.20) can be derived from (4.13). The ∼ L term comes from theẼ Q term, while the ∼ E BTBA terms originate from (2) term by exploiting the following large u expansion of the kernel: This gives the large rapidity lower bound for the energy: The convergence of the sum in (4.16) imposes a stronger constraint on the energy. Since dp Q du = O(1) at large Q, YQ must be sufficiently small for large Q. Let us investigate the large Q behavior of the summand: Since YQ is small for large Q in (4.22) the log can be expanded and at leading order one can write: As it is shown in appendix F the large Q behavior of YQ can be deduced from the BTBA equations, and it can be expressed in terms of the asymptotic solution as follows: where ξ is a real coupling dependent constant andf (u) is the complex conjugate function of f (u). The explicit functional form of f is not important except for its leading large u asymptotics: where the dots mean contributions negligible for large u. Changing variables u → Qs: and expanding all terms for large Q one gets: where the dots stand for subleading corrections in Q. The sum of the energy formula is convergent only if the energy satisfies the inequality as follows:

Numerical results
We solved the BTBA equations for the YȲ ground state numerically for various (g, L) and computed the BTBA energy by using the methods explained in Appendix G. Our numerical results are presented in Appendix G.4 and Figure 4, which we now explain in detail.
The left figure shows that the BTBA energy for the states with L = 3 2 , 2, 5 2 , 3. We examined half-integer values of L to study how the energy depends on L, although they do not correspond to the determinant-like operator (2.4) with W = Z L or V = Z L . One sees that the energy loses precision at some critical value gcr of the coupling, as indicated by huge error bars spreading out toward EBTBA = −∞. The right figure shows a plot of gcr as a function of L. A linear fit is also drawn under the assumption that gcr = 0 for the L = 1 state, because the Lüscher energy at L = 1 diverges logarithmically.
The data points and the error bars of the left figure are computed in the following way. Once we obtained a numerical solution of the BTBA for each (g, L), we calculate the BTBA energy by with Qmax = 6, instead of (4.16). Here we truncate the YQ functions at Q = Qmax to obtain a numerical solution of BTBA. Unfortunately, the truncation can induce large errors particularly around the critical point EBTBA Ecr .
To estimate the order of truncation errors, we extrapolate the energy integrals EQ using the large Q behavior (4.25). The extrapolation function is given by where the extrapolated BTBA energy is given by We solve these two equations simultaneously to determine E (fit) (Qmax + 1) and E (fit) BTBA (Qmax + 1). By repeating this procedure we obtain E

A Notation
A.1 Notation for Section 3 The notation of the paper [43] is used in Section 3 and Appendices D, E, namely The rapidity v (or u) and the momentum p are defined by and the mirror energy is given by x [+Q] x [−Q] = eǫ Q . In addition, the variable q is defined in (3.32). We use the shift operator In the integrable description, generic states are specified by M momenta {p1 , . . . , pM }, m1 type 1 fermionic roots {y1 , . . . , ym 1 } and m2 type 2 bosonic roots {μ1 , . . . ,μm 2 }. The eigenvalue of the double-row transfer matrices can be expressed by the following functions: ,

A.2 Notation for Section 4
In Section 4 and Appendices F, G, we start using another notation commonly used in the TBA equations AdS5 × S 5 (e.g. [54]), which enables the direct comparison with the literature. It should be kept in mind that the finite-size corrections are computed using the mirror region in (A.5) or the anti-mirror region in (A.1) for the respective notations.
The two conventions are related by the ±-flip, x ± → x ∓ , and where we write a bar on the variables of the notation (A.1).
The Y-functions are labeled in different ways between two sections as We may assume the left-right symmetry Ya,s = Ya,−s for the states of our concern.

B Boundary and wrapping interactions
We consider here the index structures arising from the partial contraction of the Y andȲ fields in the determinants leaving behind one or two pairs which give the boundary and wrapping interactions.

B.1 Boundary interactions
The first boundary interactions involve just one pair of Y andȲ . To evaluate it we need to know the contraction of only N − 2 indices from each determinant If we take Y (0) andȲ (x) we get a combinatorial factor of (N − 1) 2 from choosing the two fields and the resulting expression is These terms come with different powers of N . In the large N limit the interactions factorize to the individual traces in the product. When contracting Y (0) andȲ (x) with free propagators, the first term in brackets scales like N L+L ′ +4 , the second and third like N L+L ′ +3 , and the last one like N L+L ′ +2 . The first seems to dominate, but the combinatorics assumed that Tr[Y (0)Ȳ (x)] interacts with one of the other traces, otherwise it was accounted for already in (2.7) (and gets multiplied by the factor (N − 1) −1 ). We should therefore consider only graphs with interactions that involve the distinguished Y (0) andȲ (x) and some fields in the adjoint words. In this case, the first term will involve connected graphs, which do not give additional powers of N , similar to the last term. In the second and third term, however, these interactions generate planar contributions at leading order in N .
The second and third terms on the r.h.s of (B.

B.2 Wrapping corrections
Wrapping graphs come from the interaction of one of the words, like W with Y on one side andȲ on the other.
The leading wrapping corrections will arise by choosing one Y (0) and oneȲ (0) from each operator and requiring that they all interact with either W or V. With two copies of (B.1) we get the index soup 12 There are 16 possible contractions arising from this expression. Focusing on the wrapping corrections to W we take the free contractions of V andV giving The piece with the maximum number of traces is of the form The planar contractions of this expression will give another factor of N L+5 , but the combinatorics assume that there are interactions between all Y andȲ (otherwise we have to multiply the whole expression by (N − 1) −2 and recover (2.7) again). The connected correlator Tr[Y (0)Ȳ (x)] Tr[Ȳ (0)Y (x)] is independent of both W and 12 Other graphs of the same order (or lower) will involve interactions of two Y s taken from one operator and twoȲ s from the other.
These will not lead to wrapping effects, as they will all sit on one side of W or V and will not lift the energy of the vacuum.
V. Such contractions arise also in the absence of the insertions and lead to a mixing of the determinant operators themselves through the action of the full non-planar dilatation operator [55] starting at 1-loop. As mentioned in the main text, the mixing problem for such determinant-like operators with a total of 2N fields, half Y and half Y has not been solved and we will ignore this interaction term, which is not directly related to the insertions W and V.
The leading wrapping correction comes from the single trace term in (B.4), which leads to (2.9).

C Solution of the integrals
The loop integrals for the wrapping corrections are obtained by unifying in Figures 1 and 3 all fields Y (0),Ȳ (0) and Z(0) at the space-time point 0 and regarding the fields Y (x), Y (x) andZ(x) as external. This means, one removes the composite operator at x and the propagators that are connected to it and thus obtains the integral I2L, which for generic L ≥ 2 is shown on the right in Figure 3. Since I2L has an overall UV divergence, it contributes to the renormalization of the composite operator. The divergence has to be absorbed into the renormalization constant Z, which contains the negative of the sum of the UV divergencies of the diagrams. The anomalous dimension is determined by this constant, which is a matrix if mixing with other operators has to be taken into account. Here, i.e. in a CFT and for gauge invariant operators, the anomalous dimension is extracted as δ∆ = 2ελ∂ λ ln Z. Consistency of renormalization requires that ln Z is free of higher-order ε-poles.
For L = 1 the two-loop integral I2 that is found from Figure 1 and its overall UV divergence I2 read where in D = 4 − 2ε dimensions the operator K extracts the poles in ε, while R subtracts the one-loop subdivergence. The subdivergence is given by the simple one-loop integral I1 built from two propagators, which has an overall UV divergence K[I1] = λ (4π) 2 ε . In the expression for the overall UV divergence I2 in (C.1) the presence of the subdivergence is indicated by the occurrence of a quadratic ε-pole. If Z contains no one-loop contribution, then this subdivergence is in contradiction with the consistency requirement that ln Z is free of higher order ε-poles. Hence, the L = 1 state must either have a one-loop divergence, e.g. by mixing with other states, or at two-loops there must be further diagrams which cancel the quadratic ε-pole or even the entire contribution from the diagram in Figure 1. Such a one-loop mixing or additional two-loop contributions could e.g. be related to the fact that the considered gauge theory state admits interactions of the Y andȲ fields at large N . As mentioned in the main text, this state is not known and it is possible that for this state ln Z is free of higher order ε-poles.
For generic L ≥ 2, the integrals are given by the second of the figures 3. They are free of subdivergences and a unifying expression can be found for its pole parts, giving them as a function of L. 13 Our conjecture together with analytical and numerical data and the discussion with one of the authors led to [26], in which a map also of these integrals to the zig-zag integrals IZ n with n loops was presented. A conjecture for the pole parts of IZ n was made almost 20 years ago [35] and it was proven recently in [36]. In our conventions, the result of [35,36] reads In the following, we will summarize in brief the argument presented in [26], i.e. the mapping of the integrals in Figure 3 to the zig-zag integrals. We set n = L − 1 ≥ 2 such that the respective loop integral contains 13 For other examples of such integrals see e.g. [56] and [57]. 2L = 2n + 2 loops and consider the dual graph 14 Note that in the following we understand equal signs as equalities only of the overall divergencies on both sides.
For n = 1 the pole part is directly given by the one of the 4-loop zig-zag integral I4 = IZ 4 . For n = 2 we see that the dual diagram is the zig-zag integral, and hence we find I6 = IZ 6 .
For n ≥ 3, following [26], we add a vertex at ∞, and connect it with propagators to the three-valent vertices such that the integral becomes conformally invariant. This involves adding a line of negative weight between the two (n + 3)-valent vertices at zero and infinity. Then, we apply the twist identity of [58] as follows: first, we identify four vertices subject to the condition that the integral decomposes into two disconnected pieces when these points are erased. Moreover, these vertices should not be connected each by one propagator only to a common further vertex. The selected vertices will be depicted in blue. Then, we add auxiliary lines connecting the four chosen vertices in a particular way [26]. They form two sets of paired lines as will be indicated by using different colors for them. In a next step, all lines which belong to the left part of the diagram and enter the chosen points are rearranged: their endpoints are permuted among the upper and lower two pairs of selected points. This may lead to vertices that are no longer four-valent and hence conformal invariance appears to be broken. It is, however, restored in the next step where the paired auxiliary lines come into play: propagators running parallel to the auxiliary lines can be shifted to run along the respective paired auxiliary lines, thereby ensuring that the vertices become again four-valent. One step of applying this procedure, i.e. the twist identity 14 Taking the dual means going from the coordinate-space to the momentum-space representation of the diagram.
of [58], can be visualized as follows where the subgraph z 2k is obtained from z 2k−2 as follows i.e. z 2k is given by a zig-zag line with 2k + 1 triangles. The above procedure can be applied n − 2 times. After the last step, we obtain where in the second step we have removed the vertex at ∞ in the conformally invariant integral. Then, the remaining four triangles extend the zig-zag line z2n−4 with 2n − 3 triangles to a zig-zag line z2n with 2n + 1 triangles. The upper horizontal propagator connects the two two-valent vertices of this zig-zag line such that one obtains the zig-zag integral IZ 2n+2 with the known overall UV-divergence given in (C.2). Note that the result also extends to n = 1.

D One-particle Lüscher correction
In order to provide further support for the correctness of the asymptotic Y-functions obtained from the generating functional, we compute in this Appendix the Lüscher correction for the 11 particle reflecting between the YȲ boundaries in two different ways: first from the Y-functions obtained by using the generating functional, and second by appropriately modifying the "direct" Lüscher computation done for the Y Y case in [50]. These two computations should give identical results, and it is shown below that this is indeed the case.
We start by solving the boundary Bethe-Yang equations determining the momentum of the reflecting particles is obtained. Therefore, in the weak coupling limit pn = π L + 1 n for bosons , pn = π L + 1 n + 1 2 for fermions , n = 1, . . . , L .

(D.4)
Note that this is consistent with supersymmetry being broken.
For the computation using the Ya,0 functions, one can repeat the same steps taken for the Y Y case. Since S0(p, p1) is the same as for the Y Y case, eventually the complete normalization becomes (see Eqs. (4.9-4.11) of [43]) (D.5) From the weak coupling limit of (3.19) we find for the eigenvalue of the double-row transfer matrix (3.19) Ta,1 = (−1) a aqP (q, a, u(p)) where P (q, a, u) = (−3 + q 4 + 2a 2 − 8u 2 + 2q 2 (−1 + a 2 − 4u 2 ) + (a 2 + 4u 2 ) 2 ) . (D.7) Thus Ya,0 ≃ fa,1 T 2 a,1 e −2ǫaL = g 4 a 2 q 2 P 2 (q, a, u(p)) and we can compute the leading weak coupling correction to the energy of the 11 particle as Res Ya,0. (D.9) As for the "direct" Lüscher computation, recall that the approach in [50] works directly in the mirror theory where -using the boundary state formalism -this contribution can be depicted as Here c refers to the particle type whose energy correction we are calculating (which, in the present case, is 11), and the other indices run over (4a) 2 components of the a-th atypical representation of su(2|2) ⊕ su(2|2).
The boundary state amplitudes, Kl i for the Y = 0 andK jk for theȲ = 0, are related to the reflection factors (Eq. (3.4)) by analytical continuation K ij (z(q)) = C iī R j i ( ω 2 − z(q)), where z(q) is the uniformization parameter on the rapidity torus [59]. These boundary state amplitudes are the only ones in the whole computation we have to change in the YȲ case. Since both the bulk S-matrix and the boundary state amplitudes factorize as S = S0 S ⊗ S and K = K0 K ⊗ K, the energy correction can be written as In this expression we have to change only the matrix part compared to [50] when making the summation over the bound state polarizations. Decomposing the 4a dimensional bound state representation 4a = (a+1)+(a−1)+a+a as in [50], the non-vanishing components of the boundary state amplitudes that correspond to the YȲ boundaries are K 11 j,a−j =K 11 j,a−j = (−1) j , Substituting these into the sums over the bound state polarizations given in [50] leads to Tr K (q)S(q, p)K(q)S(−q, p) = − a P (q, a, u(p)) (q + 2u(p) + i(a − 1))(i(a − 1) − q + 2u(p))(2u(p) + i) 2 . (D.13) Using this together with the fact thatK0(q)K0(q) ≡ R0(z(q) − ω 2 2 )R0(−z(q) − ω 2 2 ) in place of Eq. (3.35) and the explicit form of S0(q, p) given in [50] S0(q, p) = (2u(p) we find that the expression (D.11) exactly reproduces the integrand, Ya,0, of the previous computation (D.8).
For completeness we list here the first few cases of leading Lüscher corrections computed from Eq. (D.9) (5)

E Generating function in the rotated case
In this Appendix we analyze the generating functional for the vacuum state with generic angle (E.5).

E.1 Asymptotic solution of the T-system
We now calculate the solution of the T-system (3.24) in the rotated case from the already explicitly calculated anti-symmetric transfer matrix eigenvalues: Choosing the same boundary condition we had before T0,s = Ta,0 = 1, we can easily find The expression for the symmetric transfer matrix eigenvalues has a complicated form (cf. Eq. (3.26)): where a0,s = Remarkably, these transfer matrices can be generated from the following generating functional where F (z) is given by Evidently F (z) can also be written as which is a simple deformation of 1 − z. Indeed, for θ = 0, the vacuum generating functional (E.5) reduces to W = 1 (i.e., T1,s = 0 for all s > 0, which is consistent with the supersymmetry of the Y Y system); while for θ = π/2, F (z) = 1 + z and therefore (E.5) reduces to the YȲ generating functional (3.25). We now verify that the inverse of this generating functional reproduces our previous result (E.1) for the transfer matrix eigenvalues for anti-symmetric representations: To this end, we expand F (z) in powers of z: By expanding the inverse operators, we can write (E.10) Comparing with (E.8), we can read off the following transfer matrix eigenvalues which coincides with our previous result (E.1). In passing to the last line we used that where P P denotes the principal part of the Laurent series, i.e. terms with negative powers of x.
The generating functional can be rewritten in the following form.
We conjecture that this quantity is equal to the 'dual' generating functional in the sl(2) grading, W −1 = W −1 sl(2) , just as in the Y Y case (see Appendix C of [43]).
We calculate the middle gg term as We now use the relation between the generic θ case and the θ = π/2 case: Ta,1(θ) = sin 2 θ Ta,1(π/2) , a > 0 , (E. 15) which implies that Substituting this result for W −1 (θ) into (E.14), and recalling that W −1 ( π 2 ) is given by the inverse of (3.25), we obtain (E.17) Clearly the square roots have disappeared. This function can be written as where A is the solution of the following equation Although we have not managed to find its solution for general angle, explicit solutions can be obtained at special values of θ, such as where ψ is the digamma function and Φ(−1, 1, x) ≡ ∞ k=0 (−1) k /(k + x) is the Lerch transcendent. We expect the generating functional for states in the sl(2) sector to be of the form: In this appendix the large Q behavior of the YQ functions is determined. As a first step the large m behavior of the Y m|vw functions is investigated with the help of (4.5). For large m the factor (1 + Ym) becomes unity and Y m±1|vw ≃ Y m|vw substitution can be done. Thus for lager m at leading order Y m|vw satisfies the equation: such that the large u asymptotics is given by: The solution of this equation modulo 1 m corrections is given by the asymptotic Y • m|vw functions. Thus for large index Y m|vw tend to its asymptotic counterpart. Now we can turn to investigate the large Q behavior of the "massive" YQ functions. The relevant equation to be studied is (4.9). For our considerations it is worth to convert it into its Y -system form: In the large Q limit it becomes: Now, that solution of (F.3) should be found which has the properties as follows: • The structure and positions of the local singularities within the fundamental strip 15 are the same as those of the asymptotic counterpart Y • Q . • The large u behavior of the solution is given by (4.20).
To satisfy the first requirement and solve the nontrivial part of the Y -system equation, we search the solution of (F.3) in the form: YQ = σQ Y • Q , where σQ is introduced to connect the different large u behavior of the exact and asymptotic YQs. Then it follows that σQ must be a zero mode of the LHS of (F.3): with the properties as follows: • The large u asymptotics is governed by the energy: σQ(u) = ξ u −4 E BTBA (1 + . . . ).
• σQ is real and even function.
• σQ has no zeroes or poles in the fundamental strip.
Thus σQ can be represented as a product of left and right mover modes: where ξ is a real and coupling dependent constant, the function f has the large u expansion: f (u) = u −2 E BTBA (1 + . . . ) with dots denoting negligible terms for large u andf is the complex conjugate of f . Putting everything together we get that the large Q behavior of YQ is given by formula (4.23).

G Solving the YȲ BTBA
The details of numerical computation which yielded the results in Figure 4 will be given below.
We want to solve the BTBA equations, which consist of a set of nonlinear integral equations. It is convenient to divide the whole problem into two subproblems, nonlinear root-finding and numerical integration. The nonlinear part of the problem is solved by relaxed iteration, and the integration part by interpolation and extrapolation of the integrand. Our algorithms are implemented as Mathematica scripts which are executed by CPU clusters.
For notational simplicity, we write the BTBA equations (4.5-4.13) as where we introduce the symbol (a = Q), In the first line we take the positive sign for bosonic Y's and the negative sign for fermionic Y's. 16 The normalized variables L (1 + Y ) have better analytic and numerical behavior than log(1 + Y ).

G.1 Algorithm for nonlinear problems
Iteration is one of the simplest methods for nonlinear root-finding problems. We solve the equation (G.1) by iteration of one-dimensional integrals as is close to the exact solution provided that the initial conditions Y Slower iteration algorithms are generally more stable. If one wants to get a reasonable solution from inexact initial data, a large number of iteration steps are needed. Relaxation is an example of slower algorithms, and used to solve nonlinear integral equations in the literature [60,61]. In the relaxed iteration, instead of (G.3) we update the solution by where µ (n) > 0 is a relaxation parameter. For simplicity we choose the same relaxation parameter for all Y-functions.
The updating rule (G.4) says Y (n+1) is related to Y (n) , Y (n) to Y (n−1) and so on, which is repeated until one reaches Y (0) . However, the computation using recursively-defined variables demands large memory. Thus We used ρ = 4 or 5.
The Y-functions should be updated carefully in the beginning because they change a lot after one step of iteration. This means that the relaxation parameter should be small. In fact, when we plot the energy at each step of iteration, we find that the energy typically increases for the first few steps, and then decreases monotonically. We used 0.1 ≤ µ ≤ 0.3 at the beginning of iteration, 0.3 ≤ µ ≤ 0.75 in later steps, depending on (g, L).

G.2 Algorithm for integration
Our next problem is to evaluate the set of one-dimensional integrals (G.5). In numerical analysis, one needs to approximate integrals by finite sums by choosing an appropriate distribution of sampling points {t1 , t2 , . . . tN p }.
If one wants to achieve the best precision at a fixed number of sampling points Np , one should look for the best distribution of sampling points {ti} for each integral. However, since (G.5) consists of numerous one-dimensional integrals, it is impractical to construct different sampling points {ti} for different integrals.
One solution is to choose different sampling points for different Y-functions. With this method, we construct a fitting function L (Y (t)) based on {ti} and use it to compute the integrals (G.5). This method is similar to the one used in [28], and has the advantage that we can easily keep track of the explicit shape of the Y-functions.
Let us explain our numerical integration scheme in detail. For each Y b at fixed (g, L) we introduce a rapidity cutoff M b . Then we construct a piecewise continuous interpolation function for t ≤ M b and an extrapolation function for t ≥ M b , and combine them together as and compute the integral .

(G.9)
Recall that all Y-functions are even under the parity transformation t → −t. As for Y∓ functions we only need the interpolation function for t ∈ [0, 2]. We also construct fitting functions for some kernels in the hybrid BTBA equation (4.13), namely the dressing phase kernel K Σ Q ′ Q (t, v) and s ⋆ K Q ′ −1, Q vwx (t, v) to accelerate computation. We used the third-order spline to obtain L (Y  If there are multiple solution to these equations, we use the largest one as the width. Then we choose the rapidity cutoff within the range from The distribution of sampling points for a bosonic Y-function is determined from its width. 18 Concretely, we 17 M Q > W Q means that the numerical coefficient in (G.11) is smaller than 1 20 , and similarly for other W b . 18 As for Y ∓ (t) the uniform distribution over t ∈ [0, 2] is used. used the following semi-uniform distribution for t ∈ [0, W b ]: where ∆t k = t k+1 − t k is the distance between the two adjacent sampling points. There are (Np − np) points in total for 0 ≤ t < W b . The outermost np points are used to construct the extrapolation, For small g we used t k = W b × 2 (k−Np+np) for YQ .

G.3 Numerical parameters
For clarity the numerical value of various cutoff parameters is given below.
We must truncate the number of Y-functions appearing in the BTBA equations. Let us denote the index cutoff for YQ , Y M |vw , Y M |w by Qmax , M vw|max , M w|max , respectively. We used Qmax = 6, M vw|max = 14, M w|max = 10. (G.14) The Y-functions beyond the index cutoffs are fixed at the asymptotic values.
Recall that Np is the number of total sampling points used to construct the fit of a Y-function, and np is the number of sampling points greater than the width as in (G.13). For the first argument of the kernels K Σ Q ′ Q (t, v) and s ⋆ K Q ′ −1, Q vwx (t, v), we constructed suitable distributions of sampling points in a manner similar to Appendix G. For small g we start iteration from the asymptotic Y-functions. This part of the computation is completely parallelizable. However, as g increases the finite-size corrections get larger, and one needs to start iteration from the solution at the previous step, i.e. smaller g. 20 The total number of iteration steps is typically 20 to 40 depending on (g, L).
With this parameter choice and using a node of CPU clusters with 48 cores, it took around 90 minutes to generate the dressing kernel data, and 6 hours to finish 20 steps of iteration.
Error estimates. There are three sources of errors: (i) truncation of BTBA by a finite number of Y-functions, (ii) finite number of iterations, (iii) discretization of the integrals.
The first source of errors is significant around g = gcr , as represented by the huge error bars in Figure 4.
The second source of errors makes our results unreliable at the third or fourth digits. The larger number of iterations does not always indicate the more precise results, because errors may accumulate during iterations.
The third source of errors is negligible compared to the first two. We carefully choose the distribution of sampling points for each Y-function at each (g, L), and set PrecisionGoal no less than 6 in computing various integrals in Mathematica. 19 Usually the outermost 3 or 4 points are not used, because the optimal order of extrapolation is np = 3 or 4. 20 For example, our iteration started from the asymptotic solution for g ≤ 2.4 and L = 2.

G.4 Table of numerical results
In Table 1 we present the numerical results drawn in Figure 4, in which the top end of error bars corresponds to raw data and the bottom end of error bars to the fitted data.
The raw data (Qmax = 6) are not sensitive to the large Q singularity (4.25), and thus contain the points g > gcr . When the raw data hits the large u singularity (4.21), we cannot compute the finite BTBA energy further in a reliable way. Moreover, when the energy is close to 1 4 − L, we always find that E(Q) defined in (4.22) for different Q's have comparable order of magnitude.
Note that our data also include some points with EBTBA < 1 4 −L, because we determined the large u behavior of YQ(u) by fitting the numerical data instead of using (4.20).