Quantum Spectral Curve for a Cusped Wilson Line in N=4 SYM

We show that the Quantum Spectral Curve (QSC) formalism, initially formulated for the spectrum of anomalous dimensions of all local single trace operators in N=4 SYM, can be extended to the generalized cusp anomalous dimension for all values of the parameters. We find that the large spectral parameter asymptotics and some analyticity properties have to be modified, but the functional relations are unchanged. As a demonstration, we find an all-loop analytic expression for the first two nontrivial terms in the small |phi \pm theta| expansion. We also present nonperturbative numerical results at generic angles which match perfectly 4-loop perturbation theory and the classical string prediction. The reformulation of the problem in terms of the QSC opens the possibility to explore many open questions. We attach to this paper several Mathematica notebooks which should facilitate future studies.


Introduction
The exploration of integrable structures in planar N = 4 supersymmetric Yang-Mills theory has led to numerous results which go far beyond the usual restrictions of perturbative QFT calculations [1]. With the help of integrability the spectral problem was reduced to a strikingly simple set of Riemann-Hilbert type equations known as the Quantum Spectral Curve (QSC) [2,3]. They are expected to capture the exact spectrum of single trace operator scaling dimensions and string state energies in the dual AdS 5 × S 5 theory at any value of the 't Hooft coupling λ. As a compact set of equations for only a few functions, the QSC is tremendously more efficient than the preceding infinite system of Thermodynamic Bethe ansatz (TBA) or Y-system equations [4,5,6,7,8].
The efficiency of the Quantum Spectral Curve over any other approach has been already demonstrated in a variety of settings. In [9,10] it was used to reach up to 10 loops in perturbation theory, while the all-loop near-BPS solution in [11] led to new strong coupling predictions. At finite coupling it allows to compute the spectrum numerically with extremely high precision [12]. In addition, the QSC made it finally possible to deeply probe the BFKL regime using integrability. The leading order BFKL predictions were reproduced in [13] and very recently the novel NNLO term was computed in [14]. The QSC has been formulated for the ABJM theory as well [15], leading to weak coupling results [16] and to an exact computation [17] of the interpolating function which enters all integrability-based predictions. Finally, originating in the universal QQ-relations, the QSC is also expected to be helpful in application to 3-point functions, as it should provide the exact wavefunctions in Sklyanin's separated variables.
A natural goal is to extend the Quantum Spectral Curve to various integrable deformations and boundary problems in N = 4 SYM, making possible an in-depth investigation of their properties. The TBA equations/Y-systems for these examples [18,19,20,21,22,23,24,25,26] are quite similar to those in the undeformed theory, suggesting that the QSC could also be adapted with relatively small changes.
In this paper we focus on one of the most intriguing problems of this kind, namely the boundary TBA for the generalized cusp anomalous dimension Γ cusp . This much-studied observable corresponds to the divergence associated with a pair of Wilson lines forming a cusp of arbitrary angle φ, where Λ IR,U V are the infrared and ultraviolet cutoffs. The Wilson lines include an additional coupling to the six scalars of N = 4 SYM, defined by two constant unit vectors n, n θ ∈ R 6 corresponding to the two lines, with angle θ between these vectors (see Fig. 1). We can write the cusped Wilson line explicitly as where Φ is a vector made out of the six scalars Φ i , and x q (t), xq(t) are straight lines which form an angle φ at the cusp. The generalized cusp anomalous dimension can be equivalently understood as the quark-antiquark potential on the three-sphere and is related to many physical quantities, such as radiation power from a moving quark or IR divergences in amplitudes. Recently it was also found to determine the energy levels of a supersymmetric "hydrogen atom" made out of massive W-bosons in N = 4 SYM [27].
The key insight which allowed to derive exact equations for Γ cusp from integrability was to consider the same Wilson lines with L scalars Z = Φ 1 + iΦ 2 inserted at the cusp 1 , This leads to a boundary problem for Γ cusp as the ground state energy in finite volume L with extra reflection phase factors corresponding to the two Wilson lines. The outcome is a set of TBA equations obtained in [28,29] providing the value of Γ cusp at any value of the coupling λ for arbitrary φ, θ and any number of insertions L. The usual definition of Γ cusp corresponds to L = 0. One can also consider more general insertions instead of Z L in (1.3), which are described as excited states in the TBA.
When φ = θ this observable is BPS, and one can study the near-BPS expansion in (φ − θ) which can be written as The leading coefficient in this series is known for L = 0 to all loops from localization [30,31]. For θ = 0 it was reproduced at any coupling in [32] by a direct analytic solution of the TBA, which also gave a new prediction for the case with arbitrary L. This analytic solution was extended to the case with arbitrary θ ∼ φ and any L in [33]. The near-BPS results for L ≥ 1 organize in a curious matrix model type partition function whose classical limit was investigated in [32,34] giving the corresponding classical spectral curve (see also [35,36,37]). In addition, the TBA was solved to two loops at weak coupling for finite φ, θ 1 the scalars inserted at the cusp should be orthogonal to the combinations n · Φ and n θ · Φ which couple to the Wilson lines [38], reproducing direct perturbative predictions which are also known at up to four loops [39,40,41,42].
In this paper we adapt the Quantum Spectral Curve approach to study the generalized cusp anomalous dimension at any values of the parameters. Instead of deriving the QSC from the TBA, we make a proposal based on available data and consistency of the equations, and confirm it by several highly nontrivial tests. We find that all functional equations of the QSC remain unchanged, but the asymptotics at large values of the spectral parameter, as well as some of the analyticity properties, should be modified. In particular some functions acquire exponential asymptotics ∼ e ±φu , e ±θu , as expected by analogy with spin chain Q-functions in the presence of twisted boundary conditions. We also observed that rather subtle cancelations take place resulting in complicated constraints on subleading coefficients in the large u asymptotics of Q-functions. As an application we compute the subleading term (of order (φ − θ) 2 ) in the near-BPS expansion of Γ cusp without scalar insertions, at any coupling and for any φ. Our explicit result (3.79) fully agrees with perturbative predictions.
Our paper is organized as follows. In section 2 we review the original Quantum Spectral Curve and discuss in detail the modifications needed for our problem. We discuss the vacuum state, i.e. with only Z fields inserted at the cusp, but keep L arbitrary. In section 3 we reconstruct the near-BPS solution at any θ and L, and then for L = 0 extend it to the next order in the near-BPS expansion. In section 4 we describe a highly precise numerical method for solving the QSC equations and demonstrate it on several examples. Section 5 contains a discussion of the weak coupling solution at generic angles. We present conclusions in section 6. Several appendices contain various technical details, and we also attach to this paper several Mathematica notebooks.

The Quantum Spectral Curve
In this section we first briefly review the Pµ-subsystem of the original Quantum Spectral Curve equations for the spectrum of local operators in N = 4 SYM (full details can be found in [3]). Then we discuss how these equations should be modified to describe the generalised cusp anomalous dimension. We will see that all functional equations and analyticity conditions remain the same, and the only difference is in the asymptotics at large values of the spectral parameter.
The QSC is a system of functional equations for the exact Q-functions of the theory, which are analogs of the Baxter polynomials in spin chains. The analytic properties of these functions play a key role in the construction. A particularly simple basis of Q-functions is given by 4+4 functions P a (u) and P a (u) (a = 1, 2, 3, 4). Their very nice feature is that they have only one branch cut in the complex plane of the spectral parameter u. This is a 'short' cut connecting the branch points at u = ±2g where g = √ λ 4π is related to the 't Hooft coupling λ. Knowing P a (u) and P a (u) one can reconstruct all other Q-functions, which are analytic in the upper half-plane, but typically have infinitely many branch cuts in the lower half-plane. u u Figure 2: The branch cuts of P a and µ ab . While P a have only one cut,P a have infinitely many cuts, shown by dashed lines.
While P a themselves have only one branch cut, by analytically continuing them through the cut we get new functions denoted asP a , which now have infinitely many cuts at u ∈ [−2g + in, 2g + in], n ∈ Z (see Fig. 2). Remarkably, one can get a closed system of equations by introducing an antisymmetric matrix µ ab (u) with unit Pfaffian, which relates P a to the original P-functions:P a = µ ab P b . (2.1) The discontinuity of µ ab on the cut is again expressed in terms of P a and µ ab , The functions µ ab (u) have infinitely many cuts, but are i-periodic if the cuts are chosen to be 'long' (connecting the branch points through infinity). With short cuts we have instead µ ab (u) = µ ab (u + i). There are also similar equations for P-functions with upper indices, where µ ab is the inverse matrix, µ ab µ bc = δ a c . Apart from the branch points, the Pand µ-functions have no singularities in the complex plane. Finally, P's have to satisfy For left-right symmetric states (e.g. in the sl(2) sector) the P's with lower and upper indices are not independent: The relations (2.1), (2.2), (2.3), (2.4) are known as the Pµ system. In the sl(2) sector they can be derived from the TBA equations, and in general stand as a conjecture which has been confirmed by numerous tests. It is important that the Pµ system is a closed set of equations once it is supplemented by asymptotic constraints at u → ∞ which in particular allow to compute the conformal dimension. These constraints are discussed in detail in the next section.

Asymptotics in the original QSC for local operators
To ensure uniqueness of the solution one should supplement the functional equations described above with asymptotic boundary conditions at u → ∞. The asymptotics encode the Noether charges of the state, including the conformal dimension ∆ (in our problem the analog of ∆ is Γ cusp ). Let us recall the derivation of the asymptotic constraints [3] in some detail, as this step is a crucial point in adapting the QSC to the generalized cusp.
First, all Q-functions as well as µ ab have power-like behavior at infinity , and the powers in the asymptotics of P-functions encode the SO(6) R-charges (or the angular momenta of the string on the S 5 part of the background) The charge J 1 also defines the length L = J 1 of the weak-coupling spin chain.
To constrain the leading coefficients A a , A a it is very useful to consider another four Q-functions, denoted as Q i , which are analogous to P a but correspond to dynamics in AdS 5 instead of S 5 . Like P a , these functions can be chosen to have only a single cut, but then it has to be a long cut at u ∈ (−∞, −2g] ∪ [2g, +∞). The Q i can be reconstructed from P a , P a as Q i = −P a Q a|i (2.9) with the functions Q a|i obtained by solving the difference equation where we introduced the notation This is a 4th order difference equation which mixes the Q a|i corresponding to different values of a. The index i then labels the 4 solutions of this equation. Using (2.9) we can also rewrite it as which is in fact one of the canonical QQ-relations. The matrix Q a|i should be normalized such that it preserves the χ ab matrix from (2.6), χQχQ T = 1 (2.13) and should have unit determinant [3], det 1≤a,i≤4 Q a|i = 1 . (2.14) Similarly to (2.7), the powers in the asymptotics of Q i contain the SO(4, 2) charges including the conformal dimension ∆: From consistency of the above equations we can extract the relation between ∆ and the leading coefficients A a , A a in P a , P a appearing in (2.7). To do this we expand Q a|i , P a etc as (asymptotic) series at large u, e.g.
Then from (2.12) we find at leading order Plugging this into (2.12) we see that the coefficients B j cancel and we get These equations fix the values of A a 0 A a 0 for a 0 = 1, 2, 3, 4, leading to the important relations , a 0 = 1, 2, 3, 4 . (2.20) We see that while the SO(6) charges are encoded in powers of the asymptotics of P's, the leading coefficients A a and A a contain the remaining charges S 1 , S 2 , ∆. Equations (2.7), (2.20) supplement the Pµ system relations (2.1)-(2.4) and fix the conformal dimension ∆ as a function of λ and of the other charges.
Let us finally mention that one can close the equations at the level of Q i (together with Q i which also have one long cut), introducing new functions ω ij (u) which are analogs of µ ab . The matrix ω ij is antisymmetric and has unit Pfaffian, but while µ ab are periodic with long cuts, ω ij are periodic with short cuts instead. A related useful property is that while µ ab has powerlike asymptotics at large real u, ω ij tends to a constant matrix. The equations relating ω's with Q's read where ω ij is the inverse matrix to ω ij . These are analogs of the Pµ system relations we described above. When these equations are supplemented with constraints on the asymptotics of Q's similar to (2.20), they form a closed system of equations alternative to the Pµ-system. One can also relate µ's and ω's with the help of Q-functions with four indices Q ab|ij defined as a determinant, Then we have In the next section we will discuss how to change the large u asymptotics to describe the generalized cusp.

Adapting the QSC for the cusp anomalous dimension
In this section we will discuss the modifications in the QSC which are needed to describe the generalized quark-antiquark potential. Below we will only discuss the vacuum state, i.e. the Wilson line with L scalar insertions at the cusp (the extension for more general insertions should be straightforward).
The Quantum Spectral Curve equations of [2,3] in N = 4 SYM can be deduced from the TBA equations or the corresponding T-and Y-systems with special analyticity assumptions. In our case the TBA equations for the generalised cusp are almost the same as the original TBA system. The Y-system and T-system equations are exactly the same as for the original problem. Thus it is natural to expect that the QSC equations should also be the same to a large extent. In the TBA there are only two important differences: the extra boundary dressing phase supplementing the BES phase, and the twists which appear as chemical potentials and introduce the angles φ, θ into the problem 2 . We do not derive the QSC from the Thermodynamic Bethe ansatz, rather we will put forward and motivate a proposal which is consistent with several highly nontrivial checks, leaving little doubt as to its correctness.
First, we expect to have the same set of Q-functions and auxiliary functions such as µ ab as in the original problem. All of them will satisfy the same functional relations, for instance the Pµ-system equations (2.1)-(2.4) or the QQ relations are unchanged. However some analyticity properties will change, as we will discuss below, and in particular the P-functions acquire an extra cut going from u = 0 to infinity. In addition, the large u asymptotics clearly need to be modified. Indeed, the twists in the boundary conditions typically correspond to imposing exponential rather than powerlike asymptotics for the Qfunctions (see e.g. [3] and references therein). In our case the angle θ is naturally related to the S 5 part of the geometry, which qualitatively corresponds to the P-functions, so roughly speaking we expect P a ∼ e ±θu at large u. Similarly, the angle φ is associated to AdS 5 leading to Q i ∼ e ±φu . This argument is also supported by the expectation that P's and Q's should be related in the classical limit to the quasimomenta for S 5 and AdS 5 , correspondingly. Similarly, we expect that L should enter the power in the asymptotics of P's, while the power in the asymptotics of Q's should contain ∆.
In the original QSC proposal [3] some guidance to fix the powers in the asymptotics (2.7), (2.15) came from comparison with the Asymptotic Bethe ansatz (ABA) which can be reproduced from the QSC, and also with the classical spectral curve. For our problem the ABA is also available [28,29], and another piece of information is the all-loop solution of the Pµ system to leading order in the near-BPS expansion, based on analytic solutions of the TBA [32,2,33]. In particular these solutions suggest that the large u asymptotics should contain half-integer powers coming from a √ u prefactor which the P's contain. However it turns out that there is an important subtlety -in the near-BPS limit the leading large u coefficient in P 3 , P 4 vanishes, making it not straightforward to guess the correct asymptotics even knowing the all-loop result.
The available data indicates that, first, the boundary dressing phase leads to exponential rather than powerlike asymptotics in µ ab . This was already observed in [2,33]. More precisely, we should have while other components of ω ij become zero at infinity. This translates via (2.23) into e ±2πu asymptotics in some components of the µ ab matrix. It remains to fix the powers in the asymptotics of P's and Q's, and relate their large u expansion coefficients to the charges of the state. To do this we demanded consistency of the equations (2.9), (2.12) expanded at large u. This precisely follows the logic described in the previous section. However, our case turned out to be much more tricky, in particular since some of the twists are the same (e.g. two of the P a functions scale with the same exponent ∼ e θu ) there are many subtle cancellations at the first several orders. It was also convenient at intermediate steps to use (2.14) as well as the 4th order Baxter-type difference equation on Q i with coefficients built from P a , P a -this equation follows from (2.9), (2.12) (see [13] for details on its derivation). Finally, already the near-BPS solution suggests that not all four P a are independent, e.g. P 1 (u) is equal up to a constant to P 2 (−u).
As a result, we found the following large u asymptotics: Here L is the number of scalar insertions at the cusp, while the constant C is unfixed and can be set to 1 by the rescaling symmetry as discussed below (2.31), (2.32). The coefficients should satisfy The relation which includes ∆ ≡ Γ cusp is more involved and we give its full form in Eq.
We see that in contrast to the undeformed case we need to expand P's up to fourth order at large u to extract the conformal dimension! With these asymptotic constraints the Pµ-system becomes a closed a set of equations fixing the cusp anomalous dimension. Notice that the asymptotics of P a contains half-integer powers of u. Thus P a are not as regular as in the case of local operators and should necessarily have extra cuts. Thus we require the regularity on the plane with only Zhukovsky cuts not for P a (or Q i ) but for This is an important additional analyticity condition.
Alternatively to the Pµ-system one can use the Qω system described in (2.21) which is also a closed set of equations provided the proper constraints at large u are imposed. In our case the leading asymptotics of Q i are (2.29) The coefficients in their large u expansion are constrained similarly to (2.25), (2.27), and in particular one can extract from them the R-charge L. We give the corresponding relations in Appendix B. Finally, like in the sl(2) sector of the original QSC we have due to which P a P a = 0 is satisfied automatically. It is useful to note that there is a rescaling symmetry under which while other µ ab are not changed (α is a constant). In particular with this rescaling one can set to 1 the constant C appearing in (2.25). We also have the γ-symmetry transformation [11,3] which reads with constant γ. With this transformation the coefficients in the asymptotics of P's will also change, e.g. for L = 0 The formula (2.27) for ∆ is invariant under this transformation, as it should be. As discussed above, from (2.25) we see that when φ → θ the leading coefficient in P 3 , P 4 is proportional to (φ − θ) 3/2 and thus is not visible at the leading order in the near-BPS expansion. The next coefficients b 1 , b 2 , . . . will scale as 1/(φ − θ) and thus all P a are of order √ φ − θ, as expected from the solution found in [33]. We will reconstruct this solution in the next section.
The asymptotics discussed in this section constitute our main result. They provide the crucial boundary conditions, thus concluding the reduction of the infinite TBA system of [28,29] to the finite set of QSC equations.
In the next sections we will demonstrate the usage of the QSC in several cases. We will compute at all loops the next-to-leading term in the near-BPS expansion, solve the equations numerically and also construct the leading weak coupling solution. All these calculations provide stringent tests of our proposal as well as giving new results.

Near-BPS solution
In this section we will describe the solution of the QSC in the near-BPS limit φ → θ. We will first recover the leading order solution at arbitrary θ found in [33], and then extend it to the next order. This calculation is quite similar to the iterative solution of the QSC at small spin studied in [11]. The main outcome is a prediction for the value of Γ cusp at order (φ − θ) 2 to all loops.

Leading order
In the limit φ → θ the generalized cusp anomalous dimension can be written as The first coefficient, also known as the Bremsstrahlung function, was computed at any coupling in [30,31] and later reproduced from integrability in [32,33] by a direct analytic solution of the TBA in this limit. It reads In [33] the leading near-BPS solution was obtained from the TBA and linked to the Pµsystem. Let us rederive this solution using solely the information coming from our asymptotics.
The key simplification is that P a ,P a ∼ √ φ − θ are small. This can be seen from our general asymptotics (2.25), (2.29), (B.1) where we have to send ∼ φ − θ → 0 meaning that in the near-BPS limit we get and (3.4) Notice that the leading coefficient in P 3 and P 4 tends to zero faster than the subleading ones since a 1 − b 1 ∼ 1/ , which modifies the expected behaviour at infinity in this limit (and similarly for Q 3 , Q 4 ). Thus we can write the expansion of P and µ as where the scaling is From (2.2) we see that at leading order the discontinuity of µ ab vanishes so µ (0) ab are periodic entire functions. To fix them we should look in more detail at the functions Q a|i and Q ab|ij , using (2.23) and our prescription (2.24) which states in particular that ω 12 ∼ e 2π|u| , u → ∞. For φ θ the r.h.s. of (2.12) is small so Q a|i are periodic functions. At the same time their large u asymptotics should be consistent with that of Q i and P a from (3.3), (3.4), meaning that Q a|i u N ai e ψ ai u where ψ ai can be equal to ±2θ or to 0 in our limit. From that we conlude that Q a|i must be constants. Moreover the relation (2.12), together with (3.3), (3.4) means that the only nonzero constants are In other words P a and Q i are the same in this limit after a relabeling of their indices (up to a constant factor). This is indeed an expected feature for a BPS configuration where cancellation between S 5 and AdS 5 modes is taking place. Similarly, ω ij and µ ab should coincide after the same relabeling of indices. Together with our requirement (2.24) this means that µ 12 = B 0 + B 1 e 2πu + B 2 e −2πu , µ 13 and µ 24 are constants, while other µ ab are zero. Note that since we should have a u → −u symmetry of the system, of course µ 12 should be either even or odd which further constrains these constants. More formally, from the asymptotics of the P's in (2.25) we see that these functions have the following symmetry under u → −u: where the matrix S reads 3 (3.10) Due to this we have from the first equation in (2.3) together with (2.5) (notice also that S −1 = −S and χSχ = −S). Imposing now the symmetry (3.11) we get µ 12 = A sinh(2πu) and µ 13 = −µ 24 . From Pf(µ) = 1 we also find that µ 13 = ±1. However with µ 13 = −1 we found that the equations on the P's (2.1) at leading order have no solution consistent with the asymptotics (2.25). Thus in summary we get where A is a constant. This also implies that at leading order Therefore the equations on the P's (2.1) to leading order take the form To solve them let us first introduce some notation. We have a very useful expansion where with I n being the modified Bessel function. By x(u) we denote the usual Zhukovsky variable which resolves the cut [−2g, 2g], We also have 3 here one should be careful due to the extra cut in Pa going to infinity, and we understand this equation to hold for Re(u) > 0 and let us introduce In this notation we have e.g.
(notice that applying the tilde amounts to flipping x → 1/x). We see that S + is the part of the Laurent expansion of sinh(2πu)e −2gθ(x−1/x) containing negative powers of x. We can alternatively write it as a contour integral where the contour goes counterclockwise around the unit circle.
Focussing on the case L = 0 we can now write the explicit solution of (3.14): 4 The constant A is arbitrary and is related to the constant C appearing in the asymptotics (2.25), so using the rescaling (2.31), (2.32) one can set either A or C to 1. One can check that this solution is fully consistent with the asymptotics (2.25), noting that, as discussed above, in (2.25) the leading coefficient in P 3 , P 4 vanishes and all b i ∼ 1/(φ − θ). This solution also reproduces via (2.27) the known result for ∆ at the leading order in (φ − θ), We also translated to our conventions the solution for any L constructed in [33] and we present it in appendix C. Remarkably, the result for Γ cusp extracted from this solution via our asymptotic relations (2.25), (2.27) perfectly matches the known predictions from TBA found in [33] (we have checked this explicitly for the first several values of L). This is already a nontrivial check of the proposed large u asymptotics.

Next-to-leading order
Let us now discuss how to solve the Pµ system at the next order in (φ−θ). The calculation is quite similar to the one done in [11] for the small spin expansion, so we will be brief in some cases.

Constructing µ's
We will first solve the equation for the correction to µ, which reads where P (0) are given by (3.22)-(3.25) (we understand µ as functions with short cuts). The function we will actually need is not µ itself, but rather is the discontinuity of µ ab on the cut on the real axis. Due to the relation P a P a = 0 one can just as well use µ reg instead of µ in the r.h.s. of the Pµ system equations (2.1). The key point is that µ reg does not have a cut on the real axis and is thus a much nicer function. It should satisfy Our goal is to solve this equation when the r.h.s. is composed from the leading order P's as in (3.28), in particular this means that ∆µ only has one cut at [−2g, 2g]. The formal solution to this equation can then be written as an integral operator acting on ∆µ, where the contour goes clockwise around the cut [−2g, 2g] and the kernel is This expression would work well if ∆µ decayed powerlike at infinity. A novel feature compared to [11] is that we also can have functions with exponential asymptotics e ±2θu in the r.h.s. of (3.31). For them the solution is written in a similar way but with a θ-dependent kernel. Namely, if ∆µ(u) = e ±2θu (c/x + O(1/x 2 )) we have with the kernels Here Φ is the Hurwitz-Lerch transcendent function 5 . Equivalently, One final remark is that if we need to solve (3.31) where ∆µ is not decaying at infinity, we subtract the non-decaying part of ∆µ and then solve the equation for that part separately. For example, if ∆µ = x − 1/x we can write it as The part in brackets is decaying at infinity, and a particular solution of (3.31) for the remaining part is found as As a resut, we can compute in closed form the functions µ reg at next-to-leading order in (φ − θ). Let us denote as in (3.5) Then we find µ reg(1) 12 where we use the notation Σ(h) to denote a particular solution f (u) of the equation It's also important that with these corrections to µ ab the Pfaffian constraint is still satisfied, i.e. reconstructing µ ab via µ ab = µ Having a particular solution of the equation (3.31) for µ reg , one could also add to it zero modes, i.e. i-periodic entire functions. Due to (2.23) and our prescription ω 12 ∼ e ±2πu they can be either constants or exponents e ±2πu . In [11] the zero modes were important to ensure correct asymptotics. Let us see however that in our case all zero modes can be removed by symmetries to leave the unique solution given above. First, as the µ's we construct already satisfy the parity constraint (3.11), the zero modes µ z.m. ab need to satisfy it independently. This immediately shows that the constant zero mode can only appear in µ 13 , µ 24 with µ z.m. 13 = −µ z.m. 24 = c 13 . Moreover, the exponents e ±2πu can originate only from ω 12 and can come therefore either from the O(φ−θ) correction to ω 12 or from the correction to Q ab|12 . In the first case the correction to ω 12 will multiply the leading order Q ab|12 and thus exponents will appear only in µ 12 . They are then restricted by (3.11) to appear only in the combination proportional to sinh(2πu). In the second case the leading order ω 12 will multiply the correction to Q ab|12 so again the exponents will appear as sinh(2πu). By arguments similar to those leading to (3.8), we expect the constant part of Q a|i (in the large u expansion) at order O(φ − θ) to be where for complete generality at order (φ − θ) we introduced constants in all components of Q a|i for which the twists ±φ, ±θ in the asymptotics cancel in the near-BPS limit. From (3.54) we find that out of Q ab|12 we will have a nonzero correction to the constant part only for (a, b) = (1, 2), (1,4), (2,3). Thus a zero mode proportional to sinh(2πu) could appear only in these components of µ ab . In summary, we find the following possible zero modes: where c i are constants of order O(φ − θ). From the Pfaffian constraint we get c 2 = 0. Moreover as µ 12 ∝ sinh(2πu) the constant c 1 can be set to zero by a rescaling transformation (2.32) with α = 1 + const · (φ − θ). Finally we can set c 2 to zero by making a γ transformation (2.34) with γ = const · (φ − θ) as then this zero mode will cancel against the leading order part of µ 12 (notice that µ 34 will change under this γ transformation only at order (φ − θ) 2 which is irrelevant for us). Thus we find that the solution for µ ab at NLO given above is completely general.
For µ 12 one can make yet another test of our prescription (2.24), as its exponential part is nonzero already at leading order in the near-BPS expansion. Let us use (2.23) where as at large u Q 12|12 ∼ u 2∆ we can expand This gives a prediction for the ratio of the coefficients of the log u term and the leading term. In our results for µ 12 the logarithmic part comes only from Γ·(uS + −uS − ) appearing in (3.44), so the coefficient of log u is determined by the 1/u term in (uS + − uS − ) and is straightforward to evaluate. We find where the coefficient γ precisely matches the expression (3.27) for ∆ at leading order in (φ − θ). Thus our result (3.58) agrees with the prediction (3.57).
In the next section we will use the correction to µ ab to construct P a at the next-toleading order.

Constructing P's
Having found the correction to µ ab we can now proceed and compute the correction to P a . It is convenient to parameterize them as The functions E 1 , E 3 , O 1 , O 3 then have powerlike rather than exponential asymptotics at large u. From (2.25) we find that O 1 , O 3 are odd, while E 1 , E 3 are even and at large u Plugging this parameterization into the Pµ system equations (2.1), we find the following equations on these functions: Since we have computed µ reg(1) ab in the previous section, δ a are some explicitly known functions. As in [11] we can write a particular solution of these equations as some integral operator acting on the r.h.s. First, the equatioñ where necessarilyh = h, is solved by To the particular solutions (3.73), (3.71) we can also add zero modes, i.e. solutions of the same equations with zero right-hand side. Then we get

74)
where the constants C 1 , C 2 parameterize the most general zero modes that do not violate the asymptotics (3.64). Now we can compute P (1) 4 and then solve the equations on E 1 , O 1 : Having found E n and O n , we obtain the P-functions from (3.59)-(3.62). Then we fix the constants C 1 , C 2 by imposing the asymptotic constraints (2.25), (2.26).
From the corrected P's that we have now computed we can finally extract the correction to the conformal dimension ∆ at all loops. We will present this result in the next section.

Final result
From the next-to-leading order solution of the Pµ system we constructed above, we obtain a new prediction for the generalized cusp anomalous dimension at order (φ − θ) 2 . To compare with the literature we found it convenient to bring our result to the form so that at each order we have a nontrivial function of φ . Our all-loop result reads where both integrals run clockwise around the cut [−2g, 2g] and D + = iS + (y)e −2gφx+2gφ/x+2gφy−2gφ/y We recall that S + was defined in (3.19), while the kernels Γ entering (3.79) are given in (3.33), (3.35), (3.36). It is also interesting to consider a further limit when θ is set to zero and φ is small. This corresponds to an expansion near the straight Wilson line. The first two terms in the series (3.78) scale as where the function f 2 (g) is found by expanding our result (3.79) and is given explicitly in Appendix D. Expanding also the prefactors in (3.78) we can write ∆ as a series in powers of φ 2 , (3.82) Let us now discuss several checks of our main result (3.79) at weak coupling. It is straightforward to expand it for g → 0 simply by expanding the integrand in (3.79) at weak coupling and taking the residue at u x , u y = 0. Then we can make a test against perturbative predictions known up to four loops. In general the structure at weak coupling is expected to be Our all-loop result allows to compute all coeficients γ (2) n (φ) in this expansion. Notice that at each loop order only a finite number of terms in the near-BPS expansion contribute, e.g. the two-loop result is completely determined by the first two terms in (3.78). For arbitrary φ and θ the anomalous dimension was computed directly up to two loops [39,40] giving and in [38] this data was reproduced from the TBA. We found that the weak coupling expansion of our result perfectly matches the prediction (3.87). The cusp anomalous dimension was also computed to four loops in [41,42], giving a prediction for the coefficients γ 4 (φ) which our result should reproduce. Indeed we found a perfect match with the perturbative data. The predictions of [42] are written in terms of harmonic polylogarithms, but match the expansion of our result 6 which does not include more complicated functions than Li n . At three loops our result gives  (7) In fact it is clear that at any loop order our result would generate Li n at most. The reason is that when evaluating the integral (3.79) by residues the most complicated functions that can appear are the Li n (e ±2iφ ) coming from expansion of the kernels (3.35), (3.36). As a further example we computed the novel five-and six-loop coefficients, given in Eq. (E.1) (Appendix E). We attach a Mathematica notebook which allows to reproduce these results and also systematically expand our all-loop result at weak coupling. Thus at weak coupling our result matches known predictions to four loops, which serves as a deep test of the proposed Quantum Spectral Curve equations and of our near-BPS calculation.

Numerical solution
The formulation of the problem in terms of the QSC allows for an efficient numerical analysis of Γ cusp at finite coupling. A highly precise and fast converging numerical method for solving the original QSC for local operators was proposed in [12]. Here we will describe how to modify it in the present case, and then demonstrate several applications. We will focus on the case L = 0, but we expect the discussion in this section should be valid for general L with minor changes.

The numerical algorithm
The first step is to parameterize the P-functions in terms of the Zhukovsky variable x(u). The only difference with [12] is that these functions now have exponential asymptotics, but they still have only one cut. Due to this, after extracting their exponential and leading powerlike asymptotics like in (2.25) they become power series in 1/x convergent everywhere on the main sheet. Thus we can approximate the functions f (u) and g(u) appearing in where M is some large cutoff 7 . Then we build P's from (2.25) (where we set the constant C to 1) in terms of the coefficients c 1,n , c 2,n which are the main parameters in our algorithm. 8 The next step is to close the equations which will give constraints fixing the values of these coefficients. For that we construct the functions Q a|i defined by (2.10) in terms of P-functions, As in [12] we first solve this equation analytically at large u, plugging the asymptotic expansion truncated at some cutoff K, into (4.2) and obtaining the coefficients B ai,n in terms of c 1,n , c 2,n (here θ a = ±θ and φ i = ±φ). It is important here to account for several cancellations taking place due to the asymptotics (2.25). As a result we get a good numerical approximation for Q a|i (u) when u has sufficiently large imaginary part. Then we use the exact equation (4.2) to decrease the imaginary part of u and eventually obtain the functions Q a|i in the vicinity of the real axis, when u ∈ [−2g + i/2, 2g + i/2]. Now we can build the Q-functions on the cut [−2g, 2g] via (2.12), and since Q + a|i do not have a cut on the real axis we also obtainQ i on the cut as The final step is to close the equations in terms of Q i ,Q i and find the free coefficients c 1,n and c 2,n . For that we use the very convenient trick proposed originally in [14]. Let us discuss it in some detail as this is a crucial part of the calculation. We start by noticing that Q i (u) and Q i (−u) should satisfy the same 4th order difference equation following from (2.10), (2.12) with coefficients built from P-functions as the equation is symmetric under u → −u. As we discussed in section 2.2, Eq. (2.28), it is simpler to work with where Ω j i (u) are some i−periodic functions. As Q i have a definite asymptotics with prescribed exponential part (2.29), all Ω j i (u) become constant at large u and furthermore only a few of them are nonzero at infinity, namely Ω 1 2 , Ω 2 1 , Ω 4 3 , Ω 3 4 . We also know thatq i (u) = ω ij (u)χ jk q k (u) where ω ij are i-periodic. Combining these relations we find where α i m = ω mj χ jk Ω i k are i-periodic (being built from periodic functions) and moreover analytic sinceq i (u) and q i (−u) are analytic in the lower half-plane. In addition to this, most of the functions α i m are equal to zero, because according to our prescription (2.24) from section 2.2 the only nonzero components of ω ij at infinity are ω 34 ∼ const · e 2π|u| and ω 13 , ω 24 ∼ const. Using also that most components of Ω j i are zero at large u we get from (4.7) the following equations (it's enough for us to consider only q 1 , q 4 ) where s 1 , s 4 , a, b, c are constants, and moreover a and b are nonzero as Ω 1 2 and the exponential part of ω 34 are nonzero at infinity. Applying tilde to the first equation we also get q 1 (u) = s 1q1 (−u) = (s 1 ) 2 q 1 (u) By comparing to the leading near-BPS solution where ω 12 ∝ sinh(2πu) (see Eq. (3.13)), we see that the first option is the correct one. It remains only to fix the sign of s 1 . For that let us consider the explicit solution (3.22)-(3.25) for P a in the near-BPS limit. We see that for p a = P a / √ u we havep As in the near-BPS limit we expect to identify q 1 and p 3 , comparing this relation with the first equation in (4.8) we see that we should choose s 1 = +1. In summary, we get a remarkably simple set of equations: where A is a constant and we recall that in our notation q i (u) = Q i (u)/ √ u. These are the key equations which are enough to close the system. Let us stress that they are exact and are not restricted to large u or near-BPS limit. In particular, similarly to [14] these equations should be useful for a systematic weak coupling solution. With this approach we can completely avoid computing ω ij as we are able to close the system using various Qfunctions only. Notice also that in [14] the resulting equations were similar but coefficients in the r.h.s. were all constant, while here we also have sinh(2πu). Now, finally, as we know Q i andQ i on the cut, we can evaluate both sides of (4.14), (4.15) at sampling points u k on the cut, and minimize the difference between them. More precisely, we can express the constant A from (4.15) as and we build a function which on the true solution of the QSC should be zero 9 : where Var denotes the variance, i.e. measures the deviation of the function from a constant 10 . Thus we have reduced the problem to minimization of F which is a function of our main parameters c 1,n , c 2,n . It's easy to see that F can be written as the norm of a 2N -dimensional vector where N is the number of sampling points. Therefore to find its minimum we can use the iterative Levenberg-Marquardt algorithm (an improved version of Newton's method) as in [12]. It converges rather fast and robustly, giving the values of coefficients c 1,n , c 2,n . Now we can reconstruct the P's and compute the anomalous dimension from e.g. (2.27).

Results
Let us now present the numerical results we obtained. First, we have evaluated Γ cusp for a wide range of the coupling from g = 0 up to g = 0.85 at fixed values of the angles φ = π/4, θ = 4π/10. The results are given in Table 1. A fit of our data at weak coupling gives Γ cusp φ = π 4 , θ = 4π 10 , g 0.8843331608401797458041129816 g 2 (4.18) which agrees with the analytical perturbative result of [39,40,41,42] with 10 −29 g 2 + 10 −25 g 4 + 10 −21 g 6 + 10 −18 g 8 error. The terms g 10 and g 12 above also give a numerical prediction for the five-and six-loop coefficients. One could try to get an analytic prediction for them by fitting the numerical data as a combination of some basis harmonic 9 As in (4.17) we have sinh(2πu k ) in denominator we should make sure the sampling points do not include u k = 0. We choose N sampling points as u k = 2gz k where z k are zeros of the N -th Chebyshev polynomial TN (z).
10 Var [f k ] = k |f k −f | 2 wheref is the average of all elements f k .  Table 1: Numerical data used for the plot in Fig. 3. We give the values of Γ cusp at finite coupling for φ = π/4, θ = 4π/10. Precision is decreased to fit the page. The full data set is available as attachment to the arXiv submission.
polylogarithms. This would require higher precision of course but should be possible to do (e.g. in [14] more than 60 digits of precision were reached). At strong coupling only the leading classical result is known in explicit form at generic angles. It can be extracted from [40,32] which gives the ∼ g coefficient. For φ = π 4 and θ = 4π 10 it gives Γ classical cusp 0.3122881g. Fitting our data we get which agrees nicely with the AdS/CFT prediction. Let us mention that at strong coupling it requires some effort to get high precision since we need to keep many terms in the expansion (4.1). It would be interesting to compare our result for the g 0 term with the 1-loop prediction of [40] which is written in an implicit form. One should also be able to derive the one-loop correction in a simpler and more general way by using the algebraic curve as in [43]. On Fig. 3 one can see that our data clearly interpolates between gauge and string theory results. In addition, on Fig. 4 we show our numerical data for the generalized cusp anomalous dimension at φ = π/4 for various values of θ and of the coupling. One can clearly see in particular the straight lines corresponding to the BPS regime φ = θ when Γ cusp is zero. We covered the full range of θ from −π to π, and on the plot one can see that as expected Γ cusp is a smooth and 2π-periodic function of this angle, invariant under θ → −θ.  [39,40,41,42]. Dashed lines indicate the leading strong coupling prediction for the slope of the function at g → ∞.

Weak coupling solution
In section 3 we constructed the solution of the QSC in the near-BPS limit φ − θ → 0. In this section we will describe the solution for arbitrary angles, at leading order in g. We will discuss the case L = 0.
At weak coupling the cuts degenerate into poles, but the singular part is typically suppressed by the coupling so one could expect P a to be regular at leading order. However the asymptotics (2.25) mean that we have to allow a 1/ √ u singularity in P 1 , P 2 . This leads to the ansatz Then all the coefficients are completely fixed by asymptotics (up to a rescaling (2.31)), giving where b = 2 cos θ cos φ + cos 2θ − 3 2 sin θ(cos θ − cos φ) (5.3) Figure 4: A 3d plot of Γ cusp at fixed φ = π/4 in a range of values of the coupling g and the angle θ, generated from ∼ 800 data points. We also added a semi-transparent purple plane located at Γ cusp = 0, and two red lines corresponding to the BPS configuration θ = ±φ for which Γ cusp = 0 (i.e. θ = ±π/4 in our case). and is defined in (2.26).
Let us now discuss µ ab . At leading order in the weak coupling expansion we expect that in the general expression only ω 12 will contribute, in analogy with the undeformed QSC [2,3,9] as this also what happens in the asymptotic large L regime. Based on our large u prescription ω 12 ∼ e 2π|u| and the near-BPS solution (3.12), it is natural to take In fact, for computing higher orders in the weak coupling expansion it should be better to completely avoid calculating ω ij and apply instead the equations (4.14), (4.15) we used in the numerics. For the functions Q ab|12 we can make an ansatz as polynomials whose degree is determined by the asymptotics of Q ab|12 , times e ±2θu in accordance with asymptotics again. Also, we expect that those of the functions Q ab|12 which do not have exponential asymptotics should be either even or odd. Thus we use the following ansatz: To fix the constants D i appearing here we use the difference equation on µ ab following from the Pµ-system equations (2.1), (2.2): where P a are related to P a by (2.30). This equation fixes all the constants except one, and we get Going to higher orders in g (see below) we also found that the constant R and ω 12 scale as ∼ 1/g 2 . The Q-functions can be found from the 4th order Baxter equation on Q i with coefficients built from P a (see [13] for its derivation). They turn out to be written in terms of generalized η-functions defined as η z 1 ,...,z k s 1 ,...,s k (u) ≡ n 1 >n 2 >···>n k ≥0 For the case when all twists z i are equal to 1 such functions already appeared in the weak coupling calculations of [44,9]. Importantly, all operations needed for the iterative procedure of [14] (e.g. expressing the product as a linear combination or solving equations of the kind f (u + i) − f (u) = η z 1 ,...,z k s 1 ,...,s k (u)) can be carried out for these functions as we describe in Appendix F. For future applications, we attach to this submission a Mathematica notebook implementing some of these operations on the generalized η-functions.
In terms of η-functions we found the following four linearly independent solutions of the fourth order Baxter equation: where z = e 2iφ ,z = e −2iφ . The true Q-functions should be identified with appropriate linear combinations of these four solutions.
To fix the anomalous dimension Γ cusp one needs to go to higher orders in g. This can be done using the iterative algorithm of [14] for which P a and Q i we have found serve as a starting point. Notice that the weak coupling algorithm of [9] is not directly applicable in our situation, as all P a are of the same order ∼ g 0 and none of them are small at weak coupling. In particular, none of the five independent equations among (5.7) decouple from the rest at leading order. However the universal iterative method of [14] works well, and we used it to compute the Pand Q-functions at higher orders 11 . In particular we reproduced the one-loop prediction directly at any φ and θ. The details of this calculation will be presented elsewhere. Using this method it is certainly possible to also reach much higher loops.

Conclusions
In this paper we present the modifications needed in the Quantum Spectral Curve to describe the generalized cusp anomalous dimension. We show that the main new ingredient of the boundary TBA formulation -the boundary reflection phase [28,29] -is mapped to a simple modification of the ω 12 asymptotics. Our proposal is consistent with the known near-BPS solution, and we also computed the subleading term in the near-BPS expansion at any coupling. The result matches perfectly the known perturbative predictions, providing a deep test of the QSC for this model. As supplementary material to this arXiv submission we added a Mathematica notebook with the perturbative expansion of our all-loop result. We also attach a notebook which should be useful for a perturbative solution of the QSC at generic angles, and a file with numerical data at finite coupling having ∼ 20 digits precision. Curiously, our modification of the asymptotics for the component ω 12 of the periodic anti-symmetric matrix ω ij is very similar to that needed for the analytic continuation in Lorentz spin for the twist-2 local operators where the ω 13 asymptotics was relaxed to be exponentially large [11,13,12]. It seems to be a common feature of non-local operators. It would be interesting to classify all consistent asymptotics of this kind and find the corresponding integrable observables.
Our results together with [46] elucidate the structure of the QSC for models with twisted boundary conditions. Extension to the q-deformation seems also to be straightforward. Generalization of the QSC for other boundary problems such as [22] should help to understand some of their still mysterious features.
The drastic simplification of the TBA we have achieved calls for a systematic exploration of Γ cusp in various regimes, with the hope of revealing new structures. One should now be able to reach much higher loop orders in the perturbative expansion with arbitrary φ, θ using the methods of [14,9,10], study analytically various special cases such as the "ladders" limit [41,47], and try to extend the link to matrix models observed in [32,33,34].
It would be also interesting to explore the connection to the supersymmetric hydrogenlike bound states of massive W-bosons in N = 4 SYM [27].
While a numerical solution of the TBA is additionally complicated by the infinite sums which diverge for real φ and θ [38], the simple high-precision numerical method of [12] for the QSC is applicable almost directly. Computing Γ cusp numerically in a wide range of the coupling we found perfect interpolation between gauge theory and string theory predictions. It is of course also interesting to develop a systematic analytical method for the strong coupling expansion and extend the celebrated string theory predictions [48,40].

B. Asymptotics of Q-functions
Similarly to the asymptotics of P a given in (2.25) in the main text, we found that the asymptotics of Q i have the form (with C an arbitrary constant) while Q's with upper and lower indices are related as in (2.30), The coefficients are constrained by While ∆ enters the powers in the asymptotics of Q i , the remaining conserved charge L is encoded in the large u expansion coefficients as 1 csc 2 θ csc 2 φ(cos φ − cos θ) 2 (cos θ cos φ − 1) where we denote csc θ ≡ 1/ sin θ, sec θ ≡ 1/ cos θ C. The leading near-BPS solution at any L Let us present explicitly the leading order near-BPS solution of the Pµ system at any L. It was constructed in [33] and below we write it in our conventions. Most importantly, imposing the asymptotics (2.25) and (2.26) we recovered from (2.27) the all-loop results of [33] for the near-BPS cusp anomalous dimension at nonzero L, providing a stringent test of the asymptotics we propose in this paper 12 . The solution has the following form. First, the components of µ ab are Second, the P-functions read Here A is a constant which can be set to 1 via a rescaling (2.31), (2.32) while the constant K ∼ √ θ − φ can be fixed from asymptotics (2.25), (2.26). The function F (x) is a power series which satisfies and is fixed as where [f ] + denotes the part of the Laurent expansion of f (x) with positive powers of x. Finally, the Laurent polynomial P L (x) reads 12 We checked the matching explicitly for the first several L's where Notice also that From this solution using (2.27) we recover the result of [33] for the cusp anomalous dimension,

E. Weak coupling predictions at five and six loops
From our all-loop result (3.79) it is straightforward to obtain a prediction for a part of the full anomalous dimension at five and six loops, namely for the coefficients γ

F. Generalized η-functions
We found that the solution of the QSC for arbitrary angles at weak coupling involves the following generalized η functions η z 1 ,...,z k s 1 ,...,s k (u) ≡ For the case when all twists z i are set to 1, the η-functions were encountered in the weak coupling computations of [44,9]. In our calculation of Γ cusp we had to deal with the case where twists are present. Below we summarize some useful relations analogous to those found in [44,9].
Let us denote a solution of the equation A useful property is η Z,z A,a = Zz(η Z,z A,a ) [2] + Z (η Z A ) [2] u a (F.5) where A is a set of indices A i and Z in the superscript is a set of twists Z i , while z is a single complex number. The prefactor Z in the r.h.s. denotes the product i Z i . Using this relation we find In these expressions a, s = 1, 2, 3, . . . while n is arbitrary. Finally we have the 'stuffle' relations which express a product of two η functions as a linear combination of some other η's. They are obtained by splitting the region of summation in the product of η functions and are directly analogous to those for polylogarithms or mutiple zeta values (see e.g. the pedagogical review [49] and references therein): where in case two of the s indices are combined in the r.h.s. the corresponding twists are mutiplied, exactly as in the stuffle relations for polylogarithms. For example, η w 2 η z 3 = η wz 5 + η w,z 2,3 + η z,w 3,2 (F.10) The operations described above are essential for the iterative procedure of [14] and should allow to run it to very high orders in the weak coupling expansion with any φ, θ.