Fast QSC Solver: tool for systematic study of N=4 Super-Yang-Mills spectrum

Integrability methods give us access to a number of observables in the planar N=4 SYM. Among them, the Quantum Spectral Curve (QSC) governs the spectrum of anomalous dimensions. Low lying states were successfully studied in the past using the QSC. However, with the increased demand for a systematic study of a large number of states for various applications, there is a clear need for a fast QSC solver which can easily access a large number of excited states. Here, we fill this gap by developing a new algorithm and applied it to study all 219 states with the bare dimension $\Delta_0 \leq 6$ in a wide range of couplings. The new algorithm has an improved performance at weak coupling and allows to glue numerics smoothly the available perturbative data, resolving the previous obstruction. Further ~ 8-fold efficiency gain comes from C++ implementation over the best available Mathematica implementation. We have made the code and the data to be available via a GitHub repository. The method is generalisable for non-local observables as well as for other theories such as deformations of N=4 SYM and ABJM. It may find applications in the separation of variables and bootstrability approaches to the correlation functions. Some applications to correlators at strong coupling are also presented.

Recently, new methods became available which utilise a combination of the conformal bootstrap techniques with the insights from integrability to obtain observables beyond the spectrum [12][13][14][15]. The information about the non-perturbative spectrum combined with the conformal bootstrap method is already shown to give very accurate bounds on the structure constants in some examples [12,13].
The non-perturbative spectrum can only be studied numerically, except for some limiting cases, such as in the near-BPS regime [2,3,8,16] or in the limit of infinitely long operators, using the Asymptotic Bethe Ansatz (ABA) [17,18]. There is very little chance of an exact analytic expression for all states in the theory, which is already quite clear from the cumbersome structure of the weak coupling expansion for the states which are tractable analytically [10,11].
The numerical methods to tackle the non-perturbative spectrum of short operators was pioneered in [1] on the example of the simplest non-trivial Konishi operator, and was based on the infinite set of integral equations called the Thermodynamic Bethe Ansatz (TBA).
The TBA based approach has a number of problems: the linear convergence rate, high computational cost, and finally going beyond the simplest set of operators was proven to be complicated [19]. With the discovery of the Quantum Spectral Curve (QSC) [20,21], which is a finite system of equations, both analytic techniques [10,11,[22][23][24][25] and numerical methods [9,26] were developed to go well beyond the state of the art of what was possible at the level of TBA.
In this paper we take the previous developments of the numerical approaches to the non-perturbative spectrum to an industrial level. By adding new ideas to improve its performance at weak coupling, where the perturbative data is readily available [11], to start up the numerics for a huge number of states, and by using fast C++ implementation of our new algorithm. We computed all 219 states with bare dimension ∆ 0 ≤ 6 on one PC (with two multi-core CPUs) in a timescale of a few months, which also included the time for adjusting our algorithm.
Details on the results. In addition to making the source code available via GitHub on https://GitHub.com/julius-julius/qsc alongside the additional tools to test the C++ compilation, to generate the starting points and to control the parameters during the run.
We also generated and share the following data: for the 45 states with additional symmetries we computed the spectrum in the range g = [0, 5] 1 (λ ∈ [0, ∼ 4 000]). Their spectral plot is 1 The state 6[0 0 5 5 1 1 0 0 ]3 is computed in the range g = [0, 4.6] and will be updated. In the paper we also analyse the strong coupling behaviour by fitting of our data. When the accuracy allows, we determine analytic numbers or provide numerical expression for the strong coupling expansion otherwise.
Finally, we use the strong coupling data to find the structure constants large g coefficients, which were studied recently using conformal bootstrap [52,53], by boosting these results with the spectral data.
Structure of the paper. In section 2 we introduce the general notations for the local operator in N = 4 SYM and introduce notations for their quantum numbers; in section 3 we review the key aspect of the QSC needed for the introduction of the numerical algorithm for the QSC solver in section 4. Then in section 5 we describe the C++ implementation with the installation and usage instructions as well as the benchmarking of the performance.
In section 6 we present the data for the spectrum and its detailed analysis followed by a conclusion in section 7.

Local operators and quantum numbers
In this section we describe the notations we use through this paper for the states and for their quantum numbers. N = 4 SYM is a superconformal field theory with PSU(2, 2|4) symmetry group [54][55][56][57][58][59][60]. As such, a state in this theory is characterised by six quantum numbers. Namely, the scaling dimension: ∆, the Dynkin labels of the Lorentz group SO(3, 1): ( 1 2 ) (a.k.a. Lorentz spin labels), and the Dynkin labels of the R-symmetry group SO(6): (q 1 p q 2 ) (a.k.a. R-symmetry labels).
The Lagrangian of N = 4 SYM depends on two parameters: the Yang-Mills coupling g YM , and the rank of the gauge group SU(N ) [54,55]. The integrable structure of the theory manifests itself [61] in the planar, or large-N limit, where we take g YM → 0 and N → ∞ in such a way that the 't Hooft coupling λ ≡ g 2 YM N , is held fixed. Another scaling of the 't Hooft coupling is g ≡ √ λ 4 π , which is a standard notation in the integrability literature. We use both g and λ, preferring g for weak and finite coupling, and λ for strong coupling.
In this paper, we work exclusively with planar N = 4 SYM. As such, due to large-N factorisation [62,63] one can focus entirely on linear combinations of the single trace operators tr W 1 . . . W L , (2.1) where W I are "fields" of N = 4 SYM, defined as 3 Here, D is the covariant derivative, and Φ is a scalar, Ψ andΨ are gluinos, and F is the gluon field strength. 4 The trace is over the SU(N ) colour indices, under which all fields transform in the adjoint representation. The number of fields entering under the trace is called the length of the operator and is denoted by L, which is only well defined at zero 't Hooft coupling.
In general, the scaling dimension ∆ of a state in this theory is a non-trivial function of the 't Hooft coupling λ. As N = 4 SYM is a conformal field theory (CFT) [56,64], the 3 Here, the superscript on D is used to denote multiple applications of the covariant derivative, rather than an index associated to the symmetry group. 4 Strictly speaking, each of these have Lorentz and/or R-symmetry indices, which we have suppressed for convenience. The form of these fields with their associated indices is given in table 1. scaling dimension of an operator O can be extracted from its two-point function 5,6 O(x)O(y) = 1 |x − y| 2 ∆ , (2.3) where x and y are coordinates in the 4D space. We can separate the scaling dimension ∆ into the bare/engineering/classical dimension -denoted by ∆ 0 , which is the sum of the bare dimensions of its constituent fields -and its quantum corrections, called the anomalous dimension. The bare dimension ∆ 0 of an operator may be obtained by setting the 't Hooft coupling g = λ = 0.
Due to supersymmetry, operators organise themselves into supermultiplets which comprise all operators with the same anomalous dimension. Supermultiplets are denoted by the quantum numbers of the highest weight state (we choose the convention/grading where the state with the lowest dimension ∆ [10]) in the multiplet. All primary members of the multiplet are superdescendants of the superprimary, obtained by acting with a set of 16 supercharges on the highest weight state. Each time we act with the supercharge we either get zero or an operator with a dimension increased by 1/2.
Superconformal multiplets can be short, semi-short and long [69][70][71][72][73]. Here we focus on long multiplets as for them the scaling dimension is not protected. Semi-short multiplets require fine-tuning of the parameters and as a result are only known to exist at g = 0 [73].
To classify weak coupling solutions of the QSC, it is convenient to switch to the notation which uses the oscillator numbers following [10,[74][75][76][77]. The reason for this is that at g = 0 there is an additional "quantum number" (or rather a label): the length of a single-trace operator L. In addition to the Lorentz spin labels, and the R-symmetry labels, these quantum numbers parameterise the length L as well [10]. Therefore, we use the following set of oscillator numbers: {n b 1 , n b 2 , n f 1 , n f 2 , n f 3 , n f 4 , n a 1 , n a 2 } . (2.4) These are directly related to Lorentz spin labels and R-symmetry labels as 1 = n b 2 − n b 1 , 2 = n a 1 − n a 2 , q 1 = n f 1 − n f 2 , p = n f 2 − n f 3 , q 2 = n f 3 − n f 4 , (2.5) and parameterise the bare dimension ∆ 0 and length L of single-trace operators as 5 Assuming the correct normalisation. 6 Note that operators dressed with more indices that furnish non-trivial representations of the symmetry group PSU(2, 2|4) will have a more complicated expression for their two-point function, which are available in the standard references on CFT [65][66][67][68].
At finite coupling, the symmetry group reduces to PSU(2, 2|4) and a multiplet is described by the usual set of 6 quantum numbers {∆ , 1 , 2 , q 1 , p , q 2 } (2.5), mixing operators with different length L [78]. Even though the integer L is no longer associated with any symmetry at finite coupling, we still can use it to label the states by tracing them back to g = 0. Similarly, the bare dimension ∆ 0 is not a quantum number, but rather a useful label which can be used to distinguish different states.
Even though the oscillator numbers contain more information than the Dynkin labels, we will still find several operators sharing the same oscillator numbers. In order to distinguish them we introduce an additional label which we call the multiplicity label, it is a positive integer and denoted as sol. Therefore, in order to identify each superprimary state/operator in planar N = 4 SYM uniquely, we introduce the State ID, which comprises the set of oscillator numbers, the bare dimension and solution number. We have Another reason why oscillator numbers are particularly convenient is because they enable us to immediately associate a state/operator with its possible field content at zero coupling.
Below, we present the map between oscillator numbers and "fundamental fields". 7 We have  Konishi multiplet as The superprimary of this multiplet has the bare dimension 2 (the left subscript) and as it is the only multiplet with these quantum numbers it has unit multiplicity, and therefore sol = 1 (the right subscript).
SL(2) sector. Consider another important example of multiplets in the SL(2) sector. Each such multiplet contains an operator that is constructed out of L scalars Z, with S covariant derivatives, that can schematically 9 be represented as tr D S Z L . 10 These operators are not the highest weight states in the multiplet by our choice of convention/grading (see [10] for the details). In our notation SL(2) sector multiplets can be identified by the following form of their State ID. We have These states have quantum numbers The indices used in the r.h.s. of the second equality, are SO(6) vector indices, and take the values I = 1, . . . , 6. 9 Concretely, such operators are tr (D n 1 Z) . . . (D n 1 Z), with L k n k = S. 10 Indeed, the Konishi multiplet contains the operator tr D 2 Z 2 of this form and thus the Konishi multiplet is also an sl(2) multiplet.
Finally, we assign to each state an integer number State Number or St. No. which is just a useful label to enumerate all states in the database.

QSC: A practical user's manual
Here we describe the Quantum Spectral Curve (QSC) [20,21] a method to solve the spectral problem for planar N = 4 SYM for both analytical and numerical applications. In this section we cover the basics of the QSC and introduce key notations.

Generalities
Each state in this theory is associated with 2 8 Q-functions of a complex spectral parameter u, associated with the PSU(2, 2|4) symmetry of the theory. The Q-functions are not all independent and are related by the QQ-relations described below. In the well studied cases of the spin-chains the Q-functions (or Baxter polynomials for the compact symmetry group) are the building blocks for the wave-functions of the states in a separation of variables basis (for recent developments see [79][80][81][82][83][84][85][86][87][88]). In particular, they contain all the quantum numbers of the state in the large-u asymptotic, which includes, in particular, the exact scaling dimension ∆.
In addition to the large-u asymptotic, injecting the information about the quantum numbers, one has to impose additional analyticity constraints on the Q-functions. As functions of the complex variable u they may have branch cuts with the branch points at ±2 g + i n/2 for an integer n and g = √ λ 4π related to the 't Hooft coupling λ. The monodromy of the Q-functions through branch cuts is linked to a symmetry of the QQ-relations, which can be traced back to the crossing symmetry and which allows to constrain them to a discrete set of solutions, corresponding to the local states of the theory [21]. We describe this property below in more detail.

Algebraic properties
A Q-function is labelled by up to 8 indices, with a bar separating up to four indices from the rest. The Q-functions are absolutely anti-symmetric for each of two groups of indexes.
It is easy to see that this labelling accounts for all 2 8 Q-functions.
The advantage of this labelling of the Q-function is that the QQ-relations (a.k.a. Wronskian identities) can be written in a compact form (in the conventions of [21]): Here we use the notation f ± ≡ f (u ± i/2), lower case indices denote single indices, and upper case indices denote sets of indexes.
Additionally, we have the following determinant-type relations [21]: One defines Hodge-dual Q-functions with upper indices in the following way [21,89]: with b 1 < · · · < b 4−n and j 1 < · · · < j 4−m . Then, we see that Q ∅|∅ = Q 1234|1234 , which we further require to be equal to unity as well. 11 We can also define distinguished Q-functions with single upper-indices: Using the determinant relations (3.4) and (3.5) to write Q 1234|1234 as a 4 × 4 determinant of various Q a|i , one gets Now, since we know that Q 1234|1234 = 1, we see that Q + 1234|1234 − Q − 1234|1234 = 0. At the same time, the r.h.s. of (3.10) gives a non-trivial constraint: Combining with (3.8) and (3.9), we get Expanding the determinant (3.10) on the first row, and repackaging the 3×3 sub-determinants as Q a|i , we get From equations (3.8) and (3.13), we get 14) Combining equations (3.6) and (3.15), we obtain Equations (3.6), (3.8), (3.9) (3.12)-(3.16) will play the central role in the description of the numerical algorithm, as we describe in section 4.
The Q-functions are further constrained by the analyticity conditions, which we describe below.

Symmetries of the Q-system
The subset of the QQ-relations above (3.6)-(3.16) is invariant under the following Λ- and where x is some function of the spectral parameter u. As we show below, in order for this symmetry to remain the symmetry of the QSC, which also includes the analyticity constraints one can conclude that x(u) should be the Zhukovsky functions, defined in equation (3.30) (or its power). Note that other Q-functions outside the reduced set transform rather non-trivially under this map.
In addition one has the following linear transformation H-symmetry [21] for some constant matricesĤ andȞ. The transformation of other Q-functions can be deduced from (3.19) and can be found in [21]. We also impose detĤ = detȞ = 1 in order to preserve Q ∅|∅ = 1. We also have whereĤ a b andȞ i j are inverse and transposed matrices ofĤ b a andȞ j i .

Asymptotic properties
Like in the case of the simple spin chains, in the case of the planar N = 4 SYM, the quantum numbers of the states are hidden in the large-u behaviour of the Q-functions. Since, we can build all the Q-functions starting from the distinguished single-index ones, it is sufficient to specify the asymptotic information only for these functions. The asymptotics of P a and Q i are related to the quantum numbers of the state that they characterise. Introduce the constants powP a and powQ i , so that as u → ∞, we have where A a and B i are constants, defined below. The large-u asymptotic is given by where Λ is the ambiguity related to the symmetry (3.17), n's are the oscillator numbers of the state and the anomalous dimension γ is defined as where ∆ 0 and L can be written in terms of the oscillator numbers via (2.6). Furthermore, the asymptotics of the single-upper-index Q-functions are The constants A a and B i contain certain freedom as they can be re-defined using Hsymmetry (3.19), However the combinations A a A a and B i B i (for each a and i) are invariant under the subset of the H-symmetry, which preserves the asymptotic (3.21). Those 8 combinations can be also written explicitly in terms of the quantum numbers of the states [21].
In this paper we adopt the conventions of [11] where the individual A a and B i are fixed in terms of powP and powQ , for j = 1, 2 , , for j = 3, 4 . (3.29)

Analytic constraints
We will mainly focus on the analytic properties of P a and Q i and the Hodge dual P a and Q i . Firstly, on their main sheet, the P a and P a have only two branch points, which are located at −2g and +2g, and joined by a short branch cut, i .e. a branch cut passing through zero. As such, they can be efficiently parameterised by the Zhukovsky variable x(u) defined by As a result, we have the following convergent expansion for P which parameterises it everywhere on the main sheet and even includes some area around the cut on the next sheet: Note that explicit dependence on ∆ of the state enters through A a and A a to (3.26) and (3.27). The expansion coefficients c a,n , c a,n and ∆ are some functions of g for each particular state and contain complete information on all other Q-functions [26] as we also argue below. In particular one can reconstruct Q i and Q i . Denote the analytic continuation of P a on this sheet byP a . Due to the parametrisation (3.31) in terms of the Zhukovsky variable, we can obtain an expression forP a by just Note that the fact that the branch cut [−2 g, 2 g], which becomes a unit circle in x variable, is inside the convergence domain of the Laurent series in x guarantees fast convergence with the coefficients c a,n decaying exponentially at sufficiently large n. For numerical applications this allows to truncate the series at finite n while keeping the high precision of the parametrisation in the vicinity of the branch-cut and the whole main sheet.
One should note, that due to the H-symmetry (3.19), partially fixed by the asymptotic (3.21), there is still some freedom in re-defining the coefficients c a,n , which we fix below in a particular way, depending on a type of the state. We refer to this symmetry as gauge symmetry.
Finally, in order to constrain the infinite set of constraints c a,n , c a,n and ∆ we have to impose the following analytic conditions on the Q's: • The Q i should be analytic in the upper half-plane, and have infinitely many branch points on the lower-half plane, located at ±2g − i n, with n ∈ Z and n ≥ 0, joined by short cuts. This property is guaranteed by construction and follows from the representation of P in (3.31).
• Upon passing through the branch cut on the real axis, we obtainQ i , which are analytic in the lower half-plane, and have an infinite series of cuts on the upper half plane.
Note that the branch structure ofQ i is similar to that ofQ i (i.e. complex conjugate function), the complex conjugate of Q i . In fact, we have to impose the following "gluing" conditions [21,25,27] 12 : where the most general form of them can be written as which relateQ i andQ i . The coefficients in the gluing matrices are related as we should haveḠ i j is inverse matrix of G i j . Furthermore, G i j =Ḡ j i [27]. Therefore, only two parameters are independent, i.e.
The gluing conditions (3.32) impose strong constraints on {∆, c a,n , c a,n }, giving a discrete set of solutions (up to the gauge symmetry). The way to impose and solve this constraint numerically was initially developed in [26], and it is generalised in section 4 for all types of states. Next we review our notation for the 4 types of states, classified according to some discrete symmetries of the system. 12 This condition is essentially a reflection of the unitarity of the N = 4 SYM. However, for some applications, such as complex 't Hooft coupling analytic continuation, complex spin or Fishnet limit this may not be general enough. In general one can always construct two sets of Qi with upper and lower half analyticity and glue those sets on the cut [−2 g, 2 g] on the real axis. In the present case, when the coupling is real, the complex conjugation is simply a quick way to construct the lower-half-plane analytic set of the Q-functions for numerical applications.

Types of States and discrete Symmetries
Based on quantum numbers, we can divide the states of N = 4 SYM into four broad groups according to additional symmetries of their P a functions. For that we will use two types of symmetries: left-right symmetry and parity symmetry, which we explore below.
Left-right Symmetry. There is an additional symmetry of interchanging Q functions with their Hodge dual. This symmetry interchanges two SU(2|2) subgroups of the PSU(2, 2|4) and is usually referred to as left-right Symmetry.
The state is invariant under the left-right symmetry if there exists a combination of Hand Λ-transformations which allows to relate P a and P a in the following way where The relation (3.35) implies, via the QQ-relations that one can also assume that Q i = χ ij Q j . Thus, the anti-symmetric tensors χ ab and χ ij identifies the two PSU(2|2) subgroups reducing the symmetry to OSp(4|4) (up to a real form).
In this case we do not need to use both set of P a and P a , and we have only one set of parameters {c a,n } corresponding to the lower index, which gives considerable acceleration of the numeric.
Parity Symmetry. Another Z 2 symmetry of the local operators of N = 4 SYM is the symmetry of reversing the order of operators under the trace. In the language of the spectral parameter this corresponds to changing u to −u.
There is a distinct set of states which is invariant under this symmetry. In this case one can assume that P's have certain parity under u → −u, so that it is possible to to simplify the ansatz for P a as follows, skipping every second power of x Again this amounts to efficiently reducing the number of parameters in the system by a factor of two, considerably improving the numerical performance for these states.
Types of states. As we explained above we have two types of Z 2 symmetries: left-right symmetry and parity symmetry. Each operator in SYM can be either in a singlet or doublet representation of each of the two symmetries. Accordingly, we have 4 types of states, which require slightly different numerical treatment, as the additional reflection symmetries are beneficial for the reduction of the parameters. At the same time if the state is not invariant under a reflection then there should be another state with exactly the same dimension ∆, which may also be used to save the computational time. We introduce the notation for the four type of states as shown in the

Setting up the Numerical Problem
In this section we describe the main steps of the numerical algorithm. We setup our numerical problem essentially following [26]. The main novelty in our approach is elucidated in the next section, where we describe a method which considerably improves the performance at weak coupling, which is essential for highly excited states, that get more an more degenerate with increased bare dimension. As we describe below, this allows us to initialise the numerical procedure from weak coupling, using the perturbative QSC solver of [10,11].

Basics
We describe the main steps of the algorithm for states with no additional symmetry (type IV or general states). There are some further simplifications to the procedure if we consider states with Z 2 symmetries (left-right symmetry or parity symmetry, cf. section 3.2), and some additional subtleties appear when we relax symmetries. Here for the most general states we describe the main steps of the algorithm and will comment on deviations, depending on the type of the states, in the section 4.3. 13 There could be states invariant only w.r.t. to both reflection performed simultaneously. As these states are quire rare we have not introduce a separate code for those and they are treated as general.
Step 1: Parameters of the problem. We begin with the ansatz (3.31). We truncate the infinite sum at some fixed value of n, which we denote as cutP: c a,n x n , (4.1) Thus, the set of parameters of our problem is {∆, c a,n , c a,n }, with a = 1, . . . , 4, n = 1, . . . cutP. However, there is a nontrivial condition, following from the QQ-relations (3.12), which we use to constrain a set of cutP parameters by setting Therefore, the total number of parameters to fix from the gluing conditions (3.32) is therefore 7 cutP + 1 (which also includes ∆).
The gauge symmetry (i.e. left over H-symmetry) allows us to fix six coefficients  Step 2: Finding Q a|i at large u. Having the ansatz for P's at hand, we firstly go to the large u limit. In this limit we can find the whole set of functions P a , Q i and Q a|i asymptotically.
For instance, given the asymptotics of P a and Q i , one can read the asymptotic for Q a|i from the QQ-relation (3.16). Therefore, one can approximate Q a|i at large u as the following asymptotic expansion introducing coefficients B a|i,n and the cut-off for the number of terms in the asymptotic expansion cutQai. Our goal is to find the set of coefficients B a|i,n in terms of c a,n and c a,n . For that we use the QQ-relation (3.16) substituting the asymptotic expansions of all functions. 14 At the leading order, coefficients B a|i,0 can be fixed in terms of A a , B i , powP a and powQ i : where We can solve the linear system (4.7) to obtain the B a|i,n coefficients in terms c a,k and c a,k , for k = 1, . . . , n, and the quantum numbers of the state. As we explain in section 4.3.4 this step should be taken with extra care, as the linear system may not have a solution.
In some cases, as we will see in section 4.3, the linear system (4.7) may have zero-modes, which need to be handled with care. We explain how to do this in detail in section 4.3.
As a result of solving the linear system (4.7), we obtain Q a|i at large u up to cutQai orders in 1/u. In the next steps, we need to explore Q a|i near the branch cut.
Step 3: Approaching the cut. Equipped with Q a|i at large u, we can proceed to approximating Q a|i numerically close to the branch cut. We start by evaluating (4.5) at some particular value of u = u 0 + i QaiShift − i/2, where QaiShift is a large positive integer, so that the expansion (4.5) is reliable at this point. Then we rewrite (3.16) as This allows us to recursively reduce the magnitude of the imaginary part of the argument of Q a|i , until we reach just above the real axis. We have (with a matrix multiplication under the square brackets) (4.10) Notice that it is possible to do this because the ansatz for P a in (3.31) is convergent, and therefore we can evaluate the product P a P b at any value of u in the upper-half plane, whereas the expansion (4.5) for Q a|i is asymptotic, and only gives a good approximation at large enough Im(u).
We repeat the procedure for nPoints probe points on the cut. More precisely we use Chebyshev nodes in the interval [−2g, 2g] i.e. u 0,k = 2g cos 2k−1 2 nPoints π for k = 1, . . . , nPoints. The Chebyshev nodes are particularly convenient because: firstly, they are more dense near the branch points, where the Q-functions are the most unpredictable; secondly, they are mapped points with equal separation on the unit circle under the Zhukovsky map (3.30).
Having computed Q + a|i we can easily find both Q i and its analytic continuation under the cut on the real axisQ i from (3.15): where we remember thatP a can be obtained from (4.1) P a , by sending x →x = 1/x.
Here let us notice that one can find Q a|i directly from Q a|i by using (3.13), which states that Q a|i is minus inverse of the transpose of Q a|i . Then one constructs Q i and Q i in the same way, using In the end, we obtain lists of Q-functions: Chebyshev nodes on the cut. We denote the lists as Qidownlist, Qiuplist, Qidowntlist and Qiuptlist respectively. After we obtain the Q-functions just above the branch cut, we can impose the gluing conditions (3.32).
In the next section, we will see how the gluing conditions are used to obtain a system of equations for the free parameters of our problem which is then solved using Newton's method.

Closing the system: Newton method for harmonics on the cut
Here we describe the main novelty of the current implementation.
First let us briefly summarise the previous method by [26]. We started from the Pfunctions, solved the QQ-relations to construct Q a|i which by itself enabled us to compute the Q-functions both above and below the branch cut. Basically for any set of parameters {c} ≡ {c a,n , c a,n } and ∆, we would be able to conduct this procedure. In order to close the equations we have to impose the gluing condition, which relates where we added the parameters into arguments of the Q-functions to emphasise that they can be constructed for given values of parameters.
The method of [26] consisted in using an optimisation algorithm to adjust parameters in order to minimise the mismatch in the gluing condition evaluated at the Chebyshev nodes u m : (4.14) This works well at finite g, but loses precision quickly at very small couplings. The reason for this is that the convergence radius of the expansion (3.31) of the P-functions into inverse powers of x at weak coupling has the convergence radius 1/|x (±2 g ± i)| g, meaning that the expansion coefficients c a,n ∼ g n all vanish quickly but at a different rate. The fact that they vanish quickly of course allows for very efficient truncation, which becomes exact at a given order in small coupling, however, the difference in the scale makes the optimisation problem ill-conditioned.
In order to resolve this problem in this paper we propose an alternative method. De- such that on the solution of the QSC we should have P Updated a = P a , but "off-shell" it is a function with rather different analytic properties -it is almost guaranteed to have infinitely many cuts even on the defining sheet. Nevertheless, we can use the Zukovsky map to resolve the cut on the real axis, which produces an analyticity annulus, containing the unit circle where we can define the Laurent expansion the prefactor we introduced for convenience, without loss of generality. We notice that on-shell (i.e. when the gluing condition is satisfied) we should have all positive powers coefficients to vanish c Updated a,n = 0, n < 0, zero power coefficient should become 1 and the rest should coincide with the initial parameters c a,n .
Technically, finding coefficients c Let us note that for the Newton method we have to compute the derivative of the vector valued function F , which we do by the simplest first order finite difference formula ∂ F = 1 F (· · · + ) − F . This is of course the slowest part of the algorithm as it requires approximately 8 × cutP evaluations of the faction F . Furthermore, one should take into account the following subtlety in the technical implementation. The coefficients c 1,n and c 3,n are purely real, whereas c 2,n and c 4,n are purely imaginary. When computing numerical derivatives, these parameters need to be shifted by a small constant , as described above.
To ensure the correct direction of the shifting, for real parameters c 1,n and c 3,n the real part of equations (4.17) is taken, and for imaginary parameters c 2,n and c 4,n the imaginary part. Furthermore, this restricts the gluing coefficients α and β in (3.34) to be real.
The main steps of the algorithm are summarised in table 5.

Details of the numerical procedure for different types of states
Based on the behaviour of the states w.r.t. the two discrete symmetries (left-right symmetry and parity symmetry) we split all operators into 4 types as explained in section 3.2.
Step 1: Parameters of the problem Start with initial parameters ∆, c a,n c a,n , a = 1 . . . 4, n = 1 . . . cutP, setting some c a,n to zero by the gauge choice.
Step 2: Finding Q a|i at large u Go to the large Im u limit. Solve the linear system of coeffi- Step 3: Approaching the cut Starting from large Im u limit, shift Q a|i down to the branch cut. Get Q a|i just above the branch cut and calculate Q i ,

Newton-Fourier
• Use Steps 1-3 to find c  Accordingly we have 4 slightly different implementations of the QSC solver for each of these types and in this section we describe the specific differences.

Left-right and parity symmetric case (Type I)
This case is the most symmetric one and therefore the simplest. Here we will ojutline the simplifications in the algorithm which arise because of symmetries.
First of all, we have the least amount of independent parameters among all cases. Due to the left-right symmetry, at Step 1 we introduce coefficients only for P a : as P a can be obtained by raising the index P a = χ ab P b . Here let us note that we convert P a and P a available in the perturbative solver [10,11] into the left-right symmetric form (3.35).
The details of this conversion are provided in the appendix B. The expansion (4.18) includes only even powers because of the parity symmetry. In addition, we impose the following gauge conditions: In practice, three of the gauge conditions are already satisfied trivially, due to parity symmetry requiring that particular coefficients are zero. Because of this, we end up with only 2 cutP parameters (cf . table 4).
To start the algorithm, we need to find Q a|i at the large-u limit at Step 2. The linear system (4.7) for the coefficients B a|i,2 n where a = 1 . . . 4, i = 1 . . . 4, n = 1 . . . cutQai/2 can be solved as a function of c a,n without complications, i .e. there are no zero-modes in this case.

The step
Step 3, where we compute Q a|i does not require a modification.
The gluing conditions are simplified as they connect only Q i andQ i . The lists of Q i andQ i can be obtained as The gluing matrix (3.33) is simplified and involves the only one independent gluing constant (3.32) and noticing that β =ᾱ. To close the system and set up a numerical search we follow the procedure of section 4.1. Regarding the Newton search at the last step, the main procedure is unchanged, and is the same as that described in section 4.2, albeit with only ∆ and c a,n as parameters.

Left-right symmetric states with general parity (Type II)
For this case, the parity symmetry is relaxed but the left-right symmetry is preserved. Therefore, in Step 1 the ansatz for P a should contain both even and odd powers in the expansion: The main important difference here is due to the additional obstruction in the linear system for the asymptotic expansion of Q a|i . In the absence of the parity symmetry we notice that Q 1 and Q 2 can mix, i.e. the large-u expansion of Q 2 will be invariant under for some constant η, as well as the large-u expansion of Q 4 will be invariant under for some constant µ. In the parity symmetric case one can avoid this mixing by requiring Q i to have definite parity, which then prevent the mixing as, for instance, Q 1 and Q 2 must have opposite parity due to the asymptotic. We notice that in the left-right symmetric case, due to the gluing condition (4.21), we must have µ = −η ᾱ α , where α is the gluing constant. This freedom affects the Step 2, where we build a large-u expansion for Q a|i . The ambiguity in the definition of Q i translates into existence of a zero mode in the linear system (4.7), responsible for the redefinition of Q a|2 → Q a|2 +ηQ a|4 and Q a|4 → Q a|4 +µQ a|1 in accordance with (4.23) and (4.24). Thus the zero mode occurs at n = n η in (4.7), which, as follows from the asymptotics, should be equal to Note that we can always choose the 4 component to be 1, because the combination 2p + is always nonzero we can impose an additional constraint B 4|2,nη = 0 and solve for remaining a = 1, . . . , 3. However, in order for this reduction to be consistent, one should also impose that the r.h.s. of the linear system (4.7) shares the left null-vector with T . This is not guaranteed for arbitrary values of the parameters and imposes and additional constraint on the expansion coefficients c a,n of the P a , with n = n η .
To give the simplest example for such a constraint, for the states 4 [0 0 4 3 2 1 0 0] sol , with sol = 1, 2, the value n η = 1 and the corresponding constraint is 15 (4.28) Conceptually, this constraint means that the power-like asymptotic for Q i is consistent with the ansatz for P a , the same way as for the general type of states the coefficients A a are partially fixed by the quantum numbers.
Once Q a|i at large u is defined, we can proceed with Step 3, i.e. shifting to the cut.
Just above the cut we obtain Q i andQ i the same way as it is done for the left-right and parity symmetric case in section 4.3.1.
During the Newton search, we need to treat the coefficient fixed by the constraint as a dependent parameter. Once the independent parameters are shifted by , we need to calculate the shift in the dependent variable using the constraint (cf. in the above example, where (4.28) which can be used to fix c 1,1 ).

General and parity symmetric states (Type III)
For this case we do not consider the left-right symmetry any more, though we keep parity symmetry. The ansatz for P a and P a in Step 1 is as follows: The rest of the algorithm is the same as outlined in the section 4.1. However, we will mention a technical subtlety in implementing the Newton's search: some coefficients {c} are not independent due to the identity (4.3). Thus, we do not vary them in the Newton search and need to recalculate them after shifting the independent parameters.

General states (Type IV)
As the main steps of the algorithm have been described on the basis of this case, we just highlight two points. Complications for this case, with least symmetry are essentially a 15 There is a freedom to choose which of the coefficients are considered to be the dependent ones. In our implementation we always treat c1,n ν as a dependent parameter for all type II states. 16 There is a freedom to choose which coefficients are kept made dependent ones, in our implementation we always treat c 1,n as dependent parameters.
combination of subtleties which we have seen separately for left-right symmetric and general (i.e. for type II states), and general and parity symmetric (i.e. for type III states) cases.
In this case n η = 1 and n µ = 2. The first constraint which looks the following: and the second one is similar, but a more cumbersome expression for c 2,2 . In general, the constraint for n µ > 1 becomes nonlinear in c a,n and c a,n , however, it is linear in c a,nµ terms.
Finally, let us note that there are poles in ∆ in the constraints (4.33), (4.35). However, for states we have studied, the poles satisfy ∆ ≤ ∆ 0 and as ∆ is expected to grow monotonically with the coupling, they should not cause a problem at finite real g.

C++ Implementation
This section is devoted to give the most necessary technical details for using the C++ implementation of the QSC solver algorithm described in the previous section. We give step by step instructions on the installation and usage.
The source-code is available at GitHub alongside with the Mathematica prototype and a database for the initial points and perturbative data.

Installation
The numerical solution of the QSC equations requires higher than machine precision arithmetic. This implies, that the usual double arithmetic of C or C++ is not enough for our purposes. This is why, we used the Class Library of Numbers (CLN), which is a special library for efficient computations with all kinds of numbers in arbitrary precision. All information about this free library can be found in the library's homepage: https://www.ginac.de/CLN.
Nevertheless, for practical use of our C++ codes, knowing all subtle details is not necessary.
The CLN package is available at Unix-like operating systems 18 , but can also be used under Windows operating system, provided a Windows Subsystem for Linux (WSL) is installed 19 .
Then all linux programs can be executed from the Windows PowerShell. The only difference with respect to usual linux is that each command should be anticipated by typing wsl and a space to indicate that the WSL system is used.
The prerequisites to run the C++ codes is that, a C++ compiler and the CLN package should be installed on the computer. This requires the installation of three packages either using the package manager or from the terminal window. Namely, these are the g++ compiler 20 , and the libcln and libcln-dev packages for the usage of arbitrary precision arithmetic 21 . In Debian linux systems, they can be installed easily from terminal window by the following commands: The source codes should be compiled only once.

Parameters and outputs
The parameters of an individual run are passed as arguments to the executable code. Here under the word parameters we mean all the data, which is necessary for a proper execution.
For example: cutoff values, precision of internal computations, number of sampling points, quantum numbers of the states, initial values of ∆ and c a,n etc. The simplified syntax for an execution from a terminal window is as follows: > executable.out followed by numbers separated by space characters, such that each number corresponds to an argument of the executable code. The meaning of each argument can be read off from the example notebook files, where the string corresponding to the appropriate command for the execution is created and the C++ code is run from Mathematica as an external program. 22 Just to show an example for windows users, the equivalent command from PowerShell using the WSL, is > wsl sudo apt install libcln-dev.
Parameters. Here, we give the detailed list of arguments below, such that the numbers in front of the description of a given argument, indicates its position in the list of arguments.
1. WorkingPrecision: is an integer for specifying the precision of internal computations.
2. Integer cutoff: cutP or cutP/2 in the general parity and parity symmetric cases respectively.
5. An integer denoting the number of Chebyshev sampling points on the cut: nPoints.
6. DH: is an integer for giving the small shift = 10 −DH for computing the numerical derivatives by a first order formula.
7. precSSF: is an integer for specifying one of the quit conditions, as described below.
8. precDelta: is another integer for giving the 2nd condition for quitting the iteration cycle, as described below.
9. maxiter: is an integer giving the maximal number of iterations.
In terms of the previous three integers, one can give the condition at the fulfilment of which the program quits the iteration process and returns. This is as follows: The program returns after the n th iteration, if either n ≥ maxiter or if the following two conditions are simultaneously satisfied: a,n |c a,n − c Updated a,n | 2 + |c a,n − c a,n,Updated | 2 < 10 −precSSF , where ∆ (n) denotes the value of ∆ after the nth iteration.
10. -17. Integer oscillator quantum numbers: 18. An integer for giving the multiplicity label sol of the state, 19. A real number for the 't Hooft coupling g. with WorkingPrecision −1 zeroes followed by the digit 1. The initial value 0.1 for g would be interpreted as a machine precision number, but since in the code WorkingPrecision digit arithmetic is set, the code generates a WorkingPrecision precision number for g with a value deviating from 1/10 in the order of 10 −16 . At that random generated value of g, the code will serve with precise WorkingPrecision precision computations and results. Though, this number representation seems to be a bit cumbersome it can be simply and safely treated e.g. in Mathematica. Because of such issues, the output of the C++ code also writes out the parameters of the actual calculation.
It is also very important to note, that the C++ implementation uses a bit different from (4.1) and (4.2) conventions for the P-functions. Both in the input and the output, the coefficients c a,n and c a,n should be meant as coefficients of the series representation as follows: Thus, the C++ code requires in the list of arguments, and in the output gives, the coefficients in the convention defined by (5.1). The relation between the coefficients of the two conventions can be given by simple formulas: a,n = g powP a c a,n , c a,n,(C++) = g −powP a −1 c a,n , where c and c (C++) stands for the coefficients of (4.1)-(4.2) and (5.1), respectively. When listing the arguments, and the output of the C++ codes, for short we neglected to put the C++ label on the coefficients.
Output. The output contains numbers in a form immediately interpretable by Mathematica.
If somebody would like to get the numbers in a different form, than it can be reached after appropriate string manipulations on the output string. The structure and content of the output is explained in the example notebook files, but for completeness we summarise them here, too. 2. The value of g.

3.
Oscillator numbers and multiplicity label of the state in a vector of the form:   14. (LR 10.) A number measuring the "goodness" of the cutoffs: where Q i (u n ) QaiShift denotes the numerical value of Q i (u n ) computed with cutoff value QaiShift. This indicates how precise the solution of the QQ-system is within the given approximation.
15. (LR 11.) A quantity, which is proportional to the L 2 -norm of the deviations of the Q/Q ratios from the gluing matrix elements: This number is the indicator of the overall precision of the numerical solution, which is sensitive to both precision of the QQ-system and also measures the precision of the gluing condition fulfilment.

Management of the parameters, precision control and extrapolation
In this section we discuss our strategies for adjusting various parameters controlling precision of the solver and our adaptive strategy for selecting the step in the coupling g.
For managing and adjusting these parameters in the real time we provide a simple Python script, which is also available on GitHub .

Adjusting precision controlling parameters
The are quite a few parameters we have to adjust to ensure optimal performance within certain precision. Let us remind the main parameters: • cutP -controls the number of parameters. Bigger cutP gives better approximation of P-functions but also slows the performance as ∼ cutP 3 .
• QaiShift -controls the number of jumps from the real axis to the domain where we use the asymptotic expansion. Increasing this parameter increases the precision of the Q-functions computed for given P-functions. 25 In the left-right symmetric cases the superscripted terms are not involved.
• cutQai -controls the number of terms in the asymptotic expansion of Q a|i to use.
It also increases the precision of the Q-functions for given P-functions but due to the asymptotic nature of the large u expansion increase in this parameter has to be balanced with an increase in QaiShift.
• nPoints -number of the probe points to use. As roughly each probe points gives 8 gluing equations we can simply fix it to cutP + 2 for the general Type IV states, for Type I and III: cutP + 4 and cutP + 8 for the state of the Type II.
Step A. Optimizing QaiShift and cutQai. One has to ensure that at the given starting points for the parameters {c} and ∆ the QQ-relations are solved to a sufficiently good precision. To test this we rely on the output 14 with the goal to keep it below our targeted precision. If the output 14 is not small enough we have to increase either QaiShift or cutQai. In practice, to decide which of the parameters needs to be increased we run the C++ code with the input parameter 9 maxiter = 0, which then make the solver to stop before performing any Newton method iterations, and returns the output 14 very quickly.
Then we run it multiple times, changing QaiShift and cutQai separately to check which parameter has the biggest impact on the precision and then adjusting them to reach the desired precision of the output 14.
Step B. Optimizing cutP. This parameter needs to be changed when we are confident that the Q-system is solved to a sufficiently high precision (Step A). Then the indicator that cutP needs an increment is 1) Newton method converges with a good enough precision, but at the same time 2) the output 15, measuring the "goodness" of the gluing condition is too large. In this situation we increase cutP by 1 or by 2, depending on the type of the state, and run again.

Extrapolation and step in g
The convergence of the QSC Solver relies a lot on good set of initial points for c's and ∆.
At small g we use the perturbative QSC Solver of [10,11]. The main advantage of the current algorithm is that it is very stable at very small values of g such as 10 −5 , where the precision of the perturbative solution (which we provide for up to g 10 order) is excellent and provides good starting points for our numerical QSC Solver. This way we managed to initialise all 219 states. Then having sufficiently many numerical points (say n × 10 −5 , n = 1 . . . 10) one can use polynomial extrapolation to create starting points for the next value of g. Here we explain the strategy we used to manage the step in g.
We start from the step dg = 1/10/2 n where typically n = 7. Then if in 3 consecutive runs for different g's the convergence was achieved with less than 4 iterations of the Newton method we increase dg by factor of 2 (making sure dg is below 1/10). If the convergence failed within the limit of the iterations specified, this likely means that the starting point are too far from the true values, in this case we decrease the step by to dg → dg/2 and step back, close to the previously computed point. Finally, when increasing g we make sure that the solver runs through all values of g for which 10g is an integer.

Examples of Parameters and Benchmarking
In this section we give the details on how fast our QSC solver is for various regimes in g and different types of states. We also compare the timings for the C++ code with our Mathematica implementation for a Type I state. We were using Mac Book Pro with 2.3 GHz Quad-Core Intel Core i7 and Mathematica version 12.3.1.0 for the benchmarking. We also give typical parameters and the corresponding precision estimate for the result for ∆ for various types of states.

Data Analysis at Strong Coupling
In this section we present analysis of the numerical data obtained by the Fast QSC Solver for the first 219 states/operators -all the single trace operators with the bare dimension ∆ 0 ≤ 6. Whereas at weak coupling, the perturbative expansion can be obtained analyticonumerically by using the Perturbative QSC solver [11], at the moment, the strong coupling behaviour of these states can only be obtained by fitting the high precision numerical data.
The main content of this section thus dedicated to strong coupling analysis.
In section 6.1, we explain the fitting procedure which we use to obtain strong coupling predictions, and compare the predictions to those in the literature. In section 6.2, we analyse our database of states, to see the match to counting formulas available in the literature, and furthermore propose some predictions for strong coupling behaviour of a subset of the states we studied numerically, but also giving prediction for infinitely many states with large bare dimension. Finally, in section 6.3, we see how our spectral data can be used to initiate the bootstrability program for local operators of N = 4 SYM at strong coupling.

Strong coupling analysis of the numerical data
A convenient way to arrange the states at strong coupling was proposed in [90]. It was argued there that in the regime where the quantum numbers [ 1 2 q 1 p q 2 ] ∼ 1, the states in N = 4 SYM, should map to massive string states on AdS 5 × S 5 . In particular, in this regime, the spectrum should organise itself into "string mass levels", parameterised by a positive integer δ, as Physically, in this regime the strings carry a lot of energy so that their wave-function is localised and can only explore a small part of the curved target space. Even though the 't Hooft coupling is large, this is a quantum regime, as opposed to classical strings, for which the quantum numbers scale as [ 1 2 q 1 p q 2 ] ∼ √ λ. From (6.1) we see that the natural expansion parameter at strong coupling is the combination δ 2 λ, which indeed simplifies the structure of the expansion coefficients as we show below. Even though there is no first principle derivation based on the string theory which would allow for a systematic large λ expansion, it is expected that the structure of the asymptotic strong coupling expansion should be of the following form [1,[3][4][5][91][92][93][94]: where ∆ const and d n are undetermined coefficients. The subscript reg is used to remind that the expansion of ∆ reg /λ 1/4 has an (asymptotic) expansion in integer powers of 1/ √ λ.

Fixing ∆ const
In the absence of the general first principle derivation, in order to fix the constant term in the strong coupling expansion (6.3) we present a heuristic argument. Consider a state in N = 4 SYM, with dimension ∆, and Dynkin labels [ 1 2 q 1 p q 2 ]. The quadratic Casimir of this representation is given by [77] It is reasonable to assume that the strong coupling expansion of this operator is regular in powers of 1/ √ λ . If we substitute the form of the strong coupling expansion of ∆ from equation (6.2) into (6.4), we see that the choice ∆ const = −2 renders the strong coupling expansion of J 2 to be only in terms of ∆ 2 reg , thereby fulfilling this assumption. Based on this we conclude that for all states in planar N = 4 SYM. In the next section, we give credence to this conjecture by performing some numerical fits. Of course for some states, such as Konishi [1,[3][4][5][91][92][93][94], and those with quantum numbers [ 0 0 0] with even [52,53], this is already well established, but here we argue that should be true in general.

Fitting procedure
For each of the 219 states in our database, our goal is to make numerical fits of our data, and obtain predictions for d n in (6.3). To this end, we fit our numerical data against the truncated expansion (6.3). We also assume that ∆ const = −2. Then for each state we proceed as follows: Step 1: Determine δ. This is relatively easy to do, as for the range of 't Hooft coupling that we consider, we can clearly see how the states divide into bands with 2 √ δ λ 1/4 asymptotics, where δ = 1, 2, . . . (cf. figures 2, 3a, and 3b) Thus, we can assign a positive integral value of δ, to each state in our database. Let δ m be the value of δ associated with the state whose St. No. is m.
Step 2: Find the best hyperparameters. Let λ max; m be the highest value of the 't Hooft coupling for which we have a data point for the state m. Consider the set of 't Hooft couplings for which we have data points for this state, in between some value of λ, denoted as λ min and λ max; m . For a particular choice of λ min , we choose a range of powers of 1/ √ λ up to 1/λ nmax/2 . Thus, a choice of the tuple {λ min , n max }, defines the hyperparameters of our fit. We would like to fit our numerical data against the "model" (6.3). Let d n; m be the value of the coefficient d n , associated with the state m. Thus, we need to determine the tuple {λ min , n max }, for which the we get the best fit. For each choice of this tuple, we subtract the first two terms of (6.3), evaluated at each value of λ in between λ min and λ max; m , where we have a data point. Then we fit this "refined" data with the following powers of λ: {λ 1/4 , . . . , λ 1/4−nmax/2 }, and measure the coefficient of λ 1/4 , which should be zero, according to (6.3). We choose that tuple {λ best , n best }, for which this coefficient is has the least magnitude. As an independent verification of the precision we also usually find the sub-leading coefficient d 1 to be some simple rational numbers.
Step 4: Measure the precision of d 1 . We repeat Step 3, but use the powers: λ: In table 18, we provide the fitted values of d 1 , for every state in our database, up to the best precision, that we could obtain using this procedure.

Predictions from the literature.
Recall, from section 2, that sl(2) operators are characterised by the number of derivatives S and length L. The State ID of such operators is of the form (2.10). Their quantum numbers are given by (2.11). As described in detail in appendix C, the scaling dimension of such operators, at one-loop at weak coupling, may be obtained by solving Bethe ansatz equations. As explained there, this procedure associates a set of integral mode numbers, which can be used as an alternative to the multiplicity label to distinguish between the states with the same quantum numbers. There are 27 sl(2) states in our database.
States with n = 1. This is the most studied case. In this case, there is a prediction for d 1 in [1][2][3][4][5][91][92][93][94], a prediction for d 2 in [6], and finally a prediction for d 3 in [8], giving ∆ S, L, n=1 In our database, we found 9 states with n = 1. In table 10 below, we test the predictions for the strong coupling expansion of the scaling dimension of these states and find perfect agreement within the numerical precision of our fit. The results are given in the table 10.    (2) states in our database with n = 1. The difference between the prediction for the various d k from (6.6), and our fitted value for the same d k is given in the last three columns. In all cases, the discrepancy occurs within the error of our fit.
We can also make a fit d 4 for some of these states, in table 11 below we display those, for which we have reliable precision.
Best fit for d 4   Conjectures for states with n > 1. There are a number of conjectures about the strong coupling dimension of states with mode numbers n > 1 in the literature [2,[6][7][8]. In [6], the following expression is proposed for the dimension of a state with n = 2 or 3, though it must be noted that the authors of that article pointed out that this prediction is based on an assumption about the behaviour in S, which definitely breaks down for n > 1 state. So that is curious to see if this affects the prediction already at the low orders in 1/ √ λ, which reads Here B 2 depends on n, and we have In [7], the conjuncture (6.7) was tested for S = 2, and various values of L and n. It was found that this conjecture holds water for n = 2, and deviations were observed at n = 3.
For the states with n = 3 studied in [7], apart from the leading contribution, the formula (6.7) doesn't predict either d 1 or d 2 correctly in this case. In [7], a prediction for d 1 is given for the case that n = 3 and S = 2. We have We will now compare these results with our numerics.
States with n = 2. There are 4 states in our database with n = 2. First, let us consider states with S = 2. We display the 3 of them in table 12. For all these states our fit appear to agree within the numerical error with (6.7).   (2) states in our database with S = 2 and n = 2. The difference between the prediction for the various d k from (6.7), and our fitted value for the same d k is given in the last two columns. In all cases, the discrepancy occurs within the error of our fit.
The State Number 6, with [S L n] = [2 4 2] was studied in [6], and the states number 24 and 119, with [2 5 2] and [2 6 2] respectively, were studied in [7] and we are able to corroborate the result, that (6.7) correctly predicts the scaling dimension at strong coupling for these states.
Next, consider the 1 state in our database with S > 2 as shown in table 13.  Thus, for the State Number 204, with S = 4, we see that the formula (6.7) does not seem to work, even at sub-leading order at strong coupling. It would be interesting to understand the reason for this disagreement and find a general formula which works for all S.
States with n = 3. There is 1 state in our database with n = 3, which is displayed in table 14.  Table 14: Test of strong coupling predictions for the sl(2) state in our database with S > 2 and n = 2. The difference between the prediction for d 1 from (6.9), and our fitted value is given in the last column. The discrepancy occurs within the error of our fit.
This state, number 120, is the same state as was studied in [7]. Therefore we are only able to obtain a consistency check, rather than a confirmation of the prediction (6.9).
In general we conclude that even in the simplest sl(2) sector we still do not know an analytic way of computing even the leading non-trivial strong coupling coefficient d 1 , which would work for all states.

Counting of states and Kaluza-Klein towers
The authors of [95] perform a counting of states of N = 4 SYM, that scale as ∼ λ 1/4 at strong coupling. In order to do this, they take the flat-space limit of AdS 5 ×S 5 , and perform a counting of representations of SO(9), the massive little group of R 1,9 . In particular, they count representations of SO(4) × SO(5), corresponding to the split into AdS 5 and S 5 .
Compactifying five directions to S 5 , implies that that each representation [m n] of SO (5) gets replaced by a Kaluza-Klein (KK-) tower [96] of SO (6) representations. Introduce the following notation for the KK-tower of states with Lorentz spin labels ( 1 2 ), and SO (5) labels [m n]: [96]  States within a KK-tower will share the same δ but their ∆ will in general split already at 1/λ 1/4 order. There are infinitely many states in each KK-tower, which then implies that there should be infinitely many states for a given string level δ. However, the number of KK-towers at each δ is finite and we will give some examples below.
Thus, for each δ, one can introduce a counting function count δ , following [95], which contains the KK-towers with that value of δ. For example The count 3 is also known explicitly [97], and is given in appendix D.2. Below we identify states belonging to different KK-towers and check the counting (6.11) and (6.12).
The process of assigning a specific state to a given Kaluza-Klein tower solely based on quantum numbers can generally be ambiguous. This is due to the degeneracies in the spectrum and the absence of additional input from the wave functions. As we mentioned already, the states within the same tower could have different sub-leading coefficients d 1 .
However, we noticed that when computing the expansion of the quadratic Casimir instead of ∆, the states within the same tower have the same sub-leading coefficient, which can in turn be used to assign states to different KK-towers. Furthermore, if this observation is correct, it implies a non-trivial set relations on the d 1 coefficients in the expansion of ∆ among infinite set of states within one KK-tower. More precisely one gets the following identity: + 3 q 1 (q 1 + 4) + 3 q 2 (q 2 + 4) + 2 q 1 q 2 + j 1 2 , (6.13) where j 1 is the sub-leading coefficient in the strong coupling expansion for J 2 2 δ √ λ + j 1 , which we conjecture to be the same for all states in one tower, based on the data we have.
We also checked that the next coefficient j 2 is already different for states in a KK-tower.
States with δ = 1. There are 5 such states in our database. All of them have Lorentz spin labels (0 0). We present them in table 15.

St. No.
State ID  Clearly, all these states are members of the KK-tower [0 0; 0 0], in exact agreement with (6.11), and thus the results of [95]. Furthermore we see that for all these states j 1 = 2. All the states are sl(2) states of the form [S, L, n] = [2, p + 2, 1], as can be seen by comparing with section 6.1.3. The strong coupling dimensions of these states, at first three sub-leading orders is given by equation (6.6).
All the other states in the tower should have bigger bare dimension due to the unitarity bounds (see e.g. [98]): ∆ > 2 + max 1 + 1 2 (3 q 1 + 2 p + q 2 ) , 2 + 1 2 (q 1 + 2 p + 3 q 2 ) , (6.14) which for the current set of states would imply ∆ > p + 2 . States with higher δ. At δ > 1 we can face an ambiguity when assigning the states from our database to various KK-towers (given in (6.12) for δ = 2 and in (D.1) for δ = 3). We see that there are many states, whose quantum numbers are consistent with being a member of more than one KK-tower. But this ambiguity is lifted if we look at the j 1 coefficient. Or in other words we verified the conjecture (6.13) which relates the sub-leading coefficients in ∆ between various members of the KK-tower.
We have a total of 38 states in our database, with δ = 2. We have classified them into the 9 KK-towers given in (6.12). The assignment is described in appendix D.1, specifically in tables 20 -23. We have a total of 128 states in our database, with δ = 3. We have classified them into KK-towers given in (D.1). The assignment is described in appendix D.2, specifically in tables 24 -35. It is also shown there that j 1 is indeed a good classifier of states, into KK-towers.
In other words, where we have sufficient precision, our fits confirm the formula (6.13).
It is important to notice that, by using the formula (6.13), it is possible to obtain the d 1 coefficients for every state in the KK tower, even ones which are not practically accessible by numerical methods. 26 In the next section, we will see how this new spectral information, can be injected into the conformal bootstrap constraints obtained in [52,53].

Bootstrability
Following the philosophy of [12][13][14], the results for the spectrum, joined with the constraints from the conformal crossing relations for some correlators could give narrow bounds, or even analytic results for some structure constants. In the current set-up of the local operators our predictions at strong coupling can be fed into the conformal constraints obtained in [52,53].
One of the complications of the consideration based purely on conformal constraints is typically related to the mixing problem -when there are several states with the same quantum numbers and same δ. As we know how the degeneracy is lifted from our data, we should be able to inject this information and improve the results based on conformal bootstrap (with some input from localisation).
where Φ I is a fundamental real scalar field of N = 4 SYM, and Y I is a polarisation null vector. After extracting the kinematic factors, the four-point function of four such operators can be written in terms of a function of the conformal and R-symmetry cross ratios. We where and V are the conformal cross-ratios, and σ and τ are R-symmetry cross ratios. These are defined as Using superconformal symmetry, S can be written as [99] S(U, V ; σ, τ ) = S free (U, V ; σ, τ ) + Θ(U, V ; σ, τ ) T (U, V ) , Here S free denotes the free theory correlator, T (U, V ) is called the reduced correlator, and c is the central charge, which is given by c = N 2 −1 4 . The reduced correlator satisfies a crossing equation T (U, V ) = T (1/U, V /U ) = 1/V 2 T (U/V, 1/V ), and can be expanded using the operator product expansion (OPE). We have Here, the four-dimensional conformal block G is given by [100] In the planar limit, the superprimaries exchanged in the OPE (6.22) include single-and double-trace operators, and belong to both short and long multiplets. In [52,53] it was shown that it is possible to rewrite the crossing equation in a way that focuses on the long multiplets corresponding to the "stringy" operators exchanged in the OPE (6.22) i.e.
In order to make a comparison we introduce further notations from [52,53]. For twists Then, for the OPE coefficients, we have [52,53] The prefactor of C 2 in (6.26) is highly oscillating when λ → ∞ and becomes infinite each time the twist T crosses an even integer. At the same time the function f (λ) has a regular asymptotic expansion: The result of [52,53] includes formulas for the CFT-data of the exchanged operators at strong coupling. In order to compare we introduce an integer "Regge trajectory number" [101] t ≡ δ − /2.
In table 16, we present all the 23 states in our database, that can be exchanged in the OPE (6.22) for given δ and .  In the cases we have some or none of the exchanged states in our database, we write the St. Nos. of the states that are there, followed by a " ?" with the total number of states with that δ and in parenthesis. The counting of states in this table is from [95].
Thus, for each value of t, we have a set of states, with particular values of δ and , associated with it. These states are said to be "on the Regge trajectory t." In the papers [52,53] the mixing problem appeared to be an obstacle to disentangle various exchanged operators with the same values of δ and , only "average" formulas were obtained that compute a sum over such operators, which they denote by . . . . Namely, they compute f 0 [52], f 1 , T 1 f 0 and f 2 [53], for the states on various Regge trajectories.
Regge trajectory 1. On the leading Regge trajectory, i.e. for the states with t = 1, the degeneracy of states is trivial, up to δ = 7 [95], and therefore, the average formulas give immediately results for the coefficients in the strong coupling expansion of the CFT-data.
Correspondingly, the first 3 coefficients f 0 , f 1 and f 2 are already known analytically for these states (see [52,53] for explicit results).
Predictions on Regge trajectory 2. States on this Regge trajectory have = 2 (δ − 2). The average of the OPE coefficients is given by [52] where The average of the leading OPE coefficients, weighted by the twist is [53] Note that the number of states with a given spin contributing to (6.33) depends on the spin. For this trajectory we only have two cells in , δ table 16 with complete number of states: = 0, δ = 2 (two states) and = 2, δ = 3 (four states) according to [95] and the counting in section 6.2.
Consider first = 0 case. The two states are St. No. 3 and 4 (as defined in table 18).
For those state we have where the second index corresponds to the State Number. Then from (6.28) and (6.30) we Curiously only one of the two OPE coefficients is nonzero, which may indicate some additional simplification at strong coupling.
For = 2 there are four states with St. Nos. 196,197,198 and 199. Furthermore,198, and 199, which are Type II states, form a parity dublet, and thus their dimensions are exactly degenerate, as discussed in section 3.2. Since the external BPS operators are parity invariant, their OPE coefficients should be equal as well, which reduces the system to 3 unknowns. Using that  Note added: while this paper was at the final stage of preparation for publication we received a communication confirming that in an upcoming paper [102], analytical expressions for T 1 , T 2 and f 0 for the state numbers 3 and 4, are obtained using the analytical conformal bootstrap, and they confirm our predictions for these quantities.
For this trajectory only for = 0 and δ = 3 we have complete set of 6 states in the To further constrain the OPE coefficients one may consider correlators with non-BPS external legs and/or additional input from localisation [103,104]. Having external non-BPS operators is at the moment technically challenging, since the super-conformal blocks for these representations are not yet known.

Conclusion
We present a highly efficient numerical method and its C++ implementation for studying the spectrum of non-protected observables in planar N = 4 supersymmetric Yang-Mills theory (SYM) using the Quantum Spectral Curve (QSC) method. This work addresses the need for a more comprehensive study of a larger set of states, due to its relevance in the conformal Bootstrap approach.
Our new method facilitates efficient initialisation of numerical iterations benefiting from already available perturbative solvers at small g, overcoming a key challenge in previous computations. Furthermore, our C++ implementation delivers a 8-fold performance improvement over the most efficient existing Mathematica implementation.
We used our developed algorithm to compute all states with a bare dimension ≤ 6.
The data generated and the codes used have been made publicly available on GitHub for the wider scientific community. We provide our code as open-source and kindly request that users cite this paper when using the code or any of its original components.
By combining the strong coupling analysis of our data with the recent results from conformal bootstrap we managed to resolve the degeneracy issue in some cases and give a prediction for strong coupling expansion coefficient of a structure constant of two BPS one non-BPS operator (with bare dimension ∆ 0 = 4).
Possible future applications and directions for our approach include: • Combining our results with the constraints from conformal symmetry and numerical conformal bootstrap to establish narrow bounds on the structure constants or other observables beyond the spectrum.
• Studying other non-local observables, including light-ray operators and Regge trajectories.
• Adapting the method for use in different theories, such as ABJM.
• Modifying the approach for use in AdS 3 cases.
• Using the method to calculate the Hagedorn temperature non-perturbatively.
By making our code and data publicly accessible, we aim to aid and inspire further research in these and related systems.
A Spectral data for all states in planar N = 4 SYM with ∆ 0 ≤ 6 In this appendix we give the list of all the states which we studied, i.e. all states with bare dimension ≤ 6. In order to help with the identification of the states we provide in the table below their oscillator numbers, quantum numbers, one loop anomalous dimension and the first two strong coupling orders: sting level δ and the one loop analytic guess for the correction d 1 . The states of the types II, III and IV, which brakes some of the discrete symmetries are paired with the states with exactly the same ∆(g) in the last column. The numerical values of the ∆ can be found on GitHub repository.   Table 18: Spectral data for all states in planar N = 4 SYM. This table has the following columns.
The first two columns contain the St. No. and State ID, both of which uniquely identify a state. The third coulmn gives the quantum numbers of the state. The fourth column is the bare dimension ∆ 0 and the fifth column gives the one-loop dimension, from [11]. The sixth column gives the string mass level δ. The seventh column gives our best fit for the sub-leading dimension at strong coupling, more precisely the d 1 coefficient and where possible, the eighth column gives our prediction for this coefficient. The ninth column gives the sub-leading value of the Casimir at strong coupling: j 1 . The tenth column gives the state type, and consequently the eleventh column gives the St. Nos. of other other states degenerate to the state, if applicable.

B Conversion to the left-right symmetric form
This appendix contains technical details on how to convert the perturbative results of [10,11] into our conventions in order to be used as starting points for the numerical procedure.
The perturbative solver [10,11] does not present P a and P a in the left-right symmetric form for states which have this symmetry. In order to use the perturbative data from the solver, we need to convert it to the left-right symmetric form, i.e. so that P a = χ a b P b , by a composition of Λ-and H-transformations (cf. equations (3.17)-(3.20)).
Let P MV a and P a MV be the P-functions obtained from the perturbative QSC solver of [10,11]. Then, the Ansatz looks as follows in their conventions We use the Λ-symmetry (powP MV a → powP MV a − Λ = powP a ) in order to ensure the asymptotics of the P a and P a satisfy the left-right symmetric condition (3.35). Since for LR symmetric states we have we need to solve equations of the form and similarly for other P a to get Such a Λ can only be found if or in terms of oscillator numbers It must be noted, however, that this is a necessary and not sufficient condition for left-right where h i can be determined from (B.10) in terms of the coefficients of the P MV -functions: In the repository we plan to add a small Mathematica notebook which performs the conversion as described in this appendix.

C One-loop spectrum of sl(2) sector and mode numbers
The one-loop anomalous dimension of sl(2) operators at weak coupling may be obtained by solving the sl(2) Bethe ansatz equations (BAE) [77,105]. A multiplet is characterised by a number of derivatives S and a length L, via (2.10). To obtain the one-loop dimension, we need solve the following set of equations: the sl (2) BAE Another way of classifying solutions to the BAE is by using "mode numbers". We explain how in the sequel. Consider the logarithm of the equations (C.1). We need to specify the branch of the logarithm, and this introduces an integer constant n k , associated with the logarithm of the k th Bethe equation. We have

D Counting of states at higher δ
Here we present details of the counting of the states at strong coupling as discussed in section 6.2. In particular, we assign states in our database, with δ = 2 and δ = 3, to KK-towers available from formulas in the literature [95,97]. In some cases, where is there is an ambiguity, we use the conjecture (6.13), to break this ambiguity. Therefore, we show that the subleading Casimir j 1 , is a good classifier of KK-towers.

D.1 States with δ = 2
As we mentioned in 6.2 there are 9 KK-towers expected at the level δ = 2. For the convenience of the reader, we display the counting from [95], i.e. equation (6.12) below:    Table 20: States in our database with δ = 2 and ( 1 2 ) = (0 0). We predict that the states in rows with the same colour belong to the same KK-tower.

D.2 States with δ = 3
The KK-towers at δ = 3 are given by [97]. We have     In table 29, the entries shaded with belong to one of the 4 [2 2; 0 0] , belong to [2 2; 0 2] , belong to one of the 2 [2 2; 2 0] . We didn't find any states with Lorentz spin labels (0 3) and (3 0), exactly as expected by [97] through (D.1). The 2 states with Lorentz spin labels (1 3) and 2 states with Lorentz spin labels (3 1), are given in   In   In      In  To conclude, we observe matching with the counting at strong coupling for δ = 2 and δ = 3. In addition, we break degeneracies of KK-towers using the sub-leading Casimir j 1 which as we notice is a good classifier.