Solving the Hagedorn temperature of AdS5/CFT4 via the Quantum Spectral Curve: Chemical potentials and deformations

We describe how to calculate the Hagedorn temperature of $\mathcal{N}=4$ SYM theory and type IIB superstring theory on $AdS_5\times S^5$ via the Quantum Spectral Curve (QSC) -- providing further details on our previous letters 1706.03074 and 1803.04416. We solve the QSC equations perturbatively at weak 't Hooft coupling $\lambda$ up to seven-loop order and numerically at finite coupling, finding that the perturbative results can be expressed in terms of single-valued harmonic polylogarithms. Moreover, we generalize the QSC to describe the Hagedorn temperature in the presence of chemical potentials. Finally, we show that the Hagedorn temperature in certain deformations of $\mathcal{N}=4$ SYM theory (real-$\beta$ and $\gamma_i$ deformation) agrees with the one in $\mathcal{N}=4$ SYM theory at any value of $\lambda$.


Contents
Integrability in the context of the AdS/CFT correspondence has opened a unique window of non-perturbative understanding of the properties of strongly coupled gauge theories and string theories in strongly curved backgrounds, see e.g. Refs. [3,4] for reviews. It has been developed furthest for the so-called spectral problem, the problem of finding the scaling dimensions ∆ of composite operators in planar N = 4 SYM theory and thus the energies of the corresponding strings in type IIB superstring theory on AdS 5 × S 5 . The finite-coupling solution to the spectral problem is given by the thermodynamic Bethe ansatz (TBA) [5][6][7][8][9][10], an infinite set of integral equations. These equations have been subsequently brought into the form of finite-difference equations, the so-called Quantum Spectral Curve (QSC) equations [11,12], which can be efficiently solved perturbatively at weak coupling [13][14][15] and numerically at finite coupling [16,17]. Further applications of the QSC include the pomeron and BFKL regime [18][19][20], cusped Wilson lines [21][22][23][24][25], the quark-antiquark potential [26], color-twist operators [27] and integrable deformations of N = 4 SYM theory [28][29][30][31][32]. An important aspect of understanding gauge theories and string theories, which has received less attention in the context of integrability, concerns their thermodynamic properties. One such property that occurs on both sides of the AdS/CFT correspondence is Hagedorn behavior, an exponential growth of the density of states with the energy that leads to a pole in the planar partition function at the so-called Hagedorn temperature T H . For string theory in flat space, the Hagedorn temperature could be calculated long time ago [33]. In AdS 5 × S 5 , an analogous calculation has not been possible, due to the problems with quantizing string theory in curved spacetime. In the planar gauge theory, the Hagedorn temperature has been calculated in the free theory [34] and to the first order at weak coupling [35]. The physical interpretation of the Hagedorn temperature is that of a limiting temperature; it signals the breakdown of the low-energy description and a confinement-deconfinement-like transition on the gauge-theory side, which corresponds to the Hawking-Page transition between a gas of closed strings and a black hole on the string-theory side [34,[36][37][38].
In our letters [1,2], we have recently shown how to calculate the Hagedorn temperature of planar N = 4 SYM theory and type IIB superstring theory on AdS 5 × S 5 via integrability. Concretely, we have derived TBA equations for the Hagedorn temperature [1] and recast them into the form of the QSC [2]. Solving these equations perturbatively at weak coupling, we calculated the Hagedorn temperature up to three loops. Moreover, we solved the equations numerically at finite coupling, finding in particular that Hagedorn temperature of type IIB superstring theory on AdS 5 × S 5 asymptotes to the one on flat 10D Minkowski space (calculated in Ref. [33]) in the limit of strong coupling, i.e. vanishing curvature.
The first aim of the present paper is to provide several details on the Hagedorn QSC and its derivation from the TBA, which we deferred in our letters [1,2]. Moreover, we extend the weak-coupling solution from three loops to seven loops, finding interesting number-theoretic properties. In particular, we find that the perturbative results can be expressed in terms Calculating the spin-chain free energy at a finite temperature T for a chain of finite length L requires to consider the integrable model on a torus with circumferences 1/T and L (top). In the spectral problem, the temperature is taken to be zero, resulting in a cylinder with circumference L (left). In the case of the Hagedorn temperature, we have to take the limit L → ∞, resulting in a cylinder with circumference 1/T (right). The former and the latter situation are related by a double Wick rotation.
of so-called single-valued harmonic polylogarithms [39], which have previously occurred in scattering amplitudes and generalize the single-valued multiple zeta values that occur in the spectral problem [13,15,40]. Finally, we also provide numerical values for the first three corrections to the leading strong-coupling behavior, finding perfect agreement with the analytic calculation [41,42] of the first correction. 1 The second aim of this paper is to extend said results for the Hagedorn temperature to a class of integrable deformations of the maximally supersymmetric Yang-Mills theory and to include chemical potentials. The Hagedorn temperature for integrable deformations has been previously studied in Ref. [43]. The Hagedorn temperature of N = 4 SYM theory with chemical potentials has been studied previously in Refs. [44][45][46][47][48][49]. We show in this paper how to use the QSC to determine the Hagedorn temperature for a class of integrable deformations and in the presence of arbitrary chemical potentials.
Let us briefly stress some salient features of the TBA and QSC for the Hagedorn temperature [1,2], contrasting this case to the one in the spectral problem; see also Fig. 1. In the spectral problem, the integrable model is solved at zero temperature on a cylinder with a finite circumference L that takes into account the finite length of the operator, i.e. the finite number of fields the operator contains. This is related to solving the integral model on a cylinder with circumference 1/T by a double Wick rotation, such that the TBA can be used. While the case with finite temperature is known as the physical theory, its double-Wick-rotated version with finite length that is relevant in the spectral problem is called mirror theory. In order to calculate the Hagedorn temperature, we are interested in the physical theory at finite 1/T . The other dimension accounting for the finite length of an operator is sent to infinity as the Hagedorn singularity is governed by the high-energy limit of the density of states, where only operators with large length -or rather large classical scaling dimension -contribute; see Ref. [1] for details.
One effect of the double Wick rotation concerns the analytic structure on the QSC. The planar coupling g 2 = g 2 YM N 16π 2 enters the QSC via branch cuts in its fundamental functions. The effect of the double Wick rotation is to exchange so-called short branch cuts on the interval (−2g, +2g) with so-called long branch cuts on (−∞, −2g) ∪ (+2g, +∞), and vice versa. The QSC for the Hagedorn temperature thus exhibits the opposite branch-cut structure compared to the one for the spectral problem.
A second effect of the double Wick rotation concerns the boundary conditions. In a certain class of integrable deformations of N = 4 SYM theory, characterized by diagonal twists depending on the Cartan charges of psu(2, 2|4), the deformation parameters enter the TBA by imposing twisted boundary conditions along the direction with finite circumference L [50]. The twists result in exponential asymptotics of the fundamental functions of the QSC [21,28]. Similarly, the Hagedorn temperature and chemical potentials for the different Cartan charges of psu(2, 2|4) enter by imposing twisted boundary conditions along the direction with finite circumference 1/T . Again, we see that the situation for the Hagedorn temperature is related to the one for the spectral problem by a double Wick rotation, which exchanges these two circumferences. In particular, the Hagedorn temperature formally enters the QSC as a twist, with the consequence that the fundamental functions of the QSC exhibit exponential asymptotics as well.
We can already see at this heuristic level that introducing twisted boundary conditions along the direction with finite circumference L in the calculation of the Hagedorn temperature has no effect, as we take the limit where L goes to infinity. Thus, as we will discuss in more detail in Sec. 5, the Hagedorn temperature of the corresponding integrable deformations of N = 4 SYM theory coincides with the Hagedorn temperature in undeformed N = 4 SYM theory.
In the twisted spectral problem, the twists (as well as the coupling constants and the Cartan charges S 1 , S 2 , J 1 , J 2 , J 3 ) are fixed and considered as the input while the scaling dimension ∆ is kept free and becomes the output. The Hagedorn temperature is defined as the temperature for which the free energy per unit classical scaling dimension equals −1 [1]. In the spectral problem, the free energy is connected to the scaling dimension ∆. Thus, for the Hagedorn temperature, we are keeping the scaling dimension ∆ (and thus the free energy) fixed while the single twist encoding the Hagedorn temperature is kept free and constitutes the output.
Finally, it is worth noting that the TBA and QSC in the spectral problem strictly speaking compute the Witten index, while we are interested in the spin-chain partition function or free energy. This difference can be accounted for by a (fermionic) sign in the twists. The remainder of this paper is structured as follows: In Sec. 2, we describe in detail the QSC equations determining the Hagedorn temperature. We proceed in Sec. 3 with a description of how to solve these equations, first perturbatively at weak coupling and then numerically at finite coupling. Sec. 4 is devoted to the inclusion of chemical potentials, followed by a discussion of the Hagedorn temperature for integrable deformations of N = 4 SYM theory in Sec. 5. We include several appendices on technical details of the perturbative solution at weak coupling (App. A) as well as on the TBA equations for the Hagedorn temperature and their relation to the Y-system (App. B), T-system (App. C-D) and QSC (App. E-F).

Quantum Spectral Curve for the Hagedorn temperature
In this section, we will give a brief introduction to the Quantum Spectral Curve (QSC) and describe how it can be applied to the problem of determining the Hagedorn temperature T H . In particular, we provide several details which we had to omit in our letter [2].

Generalities
Some of the features of the QSC are universal, in such as they are only reflecting the psu(2, 2|4) symmetry of N = 4 SYM theory and occur for all observables in it to which the QSC has up to now been applied. Other features are specific to a particular observable. In this subsection, we will give an overview of the former, while discussing the latter in the subsequent subsections. For a more general review of the QSC, see e.g. Refs. [51][52][53].
The QSC is also known as analytic Q-system. This Q-system is an equivalent formulation of the T-system, which in turn is an equivalent formulation of the Y-system and the TBA. For more details on the transitions between these formulations, see App. B, C and E. The Q-system consists of 2 8 = 256 functions Q A|I (u) labeled by subsets A, B ⊂ {1, 2, 3, 4}. They are (multivariate) functions of the spectral parameter u ∈ C. The Q-functions satisfy the finite-difference equations where g 1 and g 2 are arbitrary functions. Moreover, the Q-system has a GL(4) × GL (4) symmetry, called H symmetry: where H 1 and H 2 are i-periodic matrices, and Using the relations (2.1)-(2.3), all functions Q A|I can be written in terms of the functions P a ≡ Q a|∅ , Q i ≡ Q ∅|i , Q a|i and Q ∅|∅ ; see e.g. Ref. [28] for the explicit construction. Moreover, we can use the first gauge transformation (2.4) to set Q ∅|∅ = 1, which we will do throughout. In terms of these functions, the finite-difference equation (2.2) reads It will be convenient to define Hodge-dual Q functions Q A|I via whereĀ (Ī) denotes the complement of A (I), and ǫ is the completely anti-symmetric tensor with four indices. In particular, we use the notation P a ≡ Q a|∅ and Q i ≡ Q ∅|i . The Q-functions and their Hodge duals satisfy the relations and For many cases of interest, the QSC possesses an additional symmetry -called leftright symmetry due to its manifestation in the T-system and Y-system. In the left-rightsymmetric case, the Q-functions with upper and lower indices satisfy another relation on top of Hodge duality. In particular, The left-right-symmetric case occurs for example in the spectral problem in well-studied closed subsectors and for the Hagedorn temperature in the absence of chemical potentials, or if these chemical potentials satisfy certain constraints to be discussed in Sec. 4.

Asymptotic solution
The asymptotic behavior of the different Q-functions, i.e. their behavior for large values of the spectral parameter u, depends on the particular observable that is to be determined by the QSC; and by specifying the asymptotics, one can specify the observable. Our starting point for determining the asymptotic behavior of the QSC is the T-system we derived in Ref. [1]. For large spectral parameter u, we have shown that it asymptotes to a constant T-system, which is quoted in Eqs. (C.29) and (C.30). Using the relation between the Q-system and the T-system discussed in App. E, this constant T-system can be reproduced via the following Q-functions: and where 14) The corresponding Q a|i (u) are given in App. F. Note that we have already used some of the H symmetry (2.5) to bring the asymptotic solution into this form. It is further convenient to use the remaining H symmetry to set In this gauge the asymptotic P a functions for a = 1, 2 and a = 3, 4 transform into each other under u → −u, while the asymptotic Q i (u) functions are either even or odd. The exponential asymptotics of the P a contain a sign, which is a consequence of the fact that the Hagedorn QSC is based of the partition function instead of the Witten index. This sign is understood to contain a small imaginary part to resolve the branch-cut ambiguity: (2. 16) We remark that the asymptotic behavior of the Hagedorn QSC is a special case of the one occurring for the spectral problem in twisted N = 4 SYM theory discussed in Ref. [28]; Figure 2. Branch cut structure of the QSC for the Hagedorn temperature. There exists a Riemann sheet for which Q i has a single 'short' cut and P a has a single 'long' cut. (Note that the branch-cut structure is the opposite of the one occurring for the spectral problem.) concretely, the asymptotics in Ref. [28] reduce to the ones above upon interchanging, in their notation, P a ↔ Q i , x a ↔ y i and A a ↔ B i , as well as setting, 2T H , y 1 = y 2 = y 3 = y 4 = 1 and λ a = ν i = 0. We will discuss this observation in more detail when discussing chemical potentials in Sec. 4.
At tree level, i.e. in the free theory, the T-system is constant and obtained by setting 3) in the asymptotic solution. Thus, the tree-level QSC is obtained by 3) in Eqs. (2.12)-(2.14) and (F.1)-(F.2). It is worth noting that the tree-level QSC is given by polynomials times exponential factors, while the tree-level QSC in the spectral problem typically contains also negative powers of u.

Branch cuts
Beyond the asymptotic limit, the analytic structure of the QSC is characterized by the structure of its branch cuts. Recall that the Hagedorn QSC is based on the physical TBA, while the QSC for the spectral problem is based on the mirror TBA. Since the physical and mirror theory are related by a double Wick rotation, the branch-cut structure of the Hagedorn QSC is exactly opposite to that for the spectral problem. In particular, there exists a Riemann sheet on which the Q i have a single 'short' branch cut on the interval (−2g, +2g), while the P a have a single 'long' cut on (−∞, −2g) ∪ (+2g, ∞); see Fig. 2.
Upon analytic continuation of P a across the single long cut, one arrives on a Riemann sheet on which P a is analytic in the upper half plane (UHPA) but possesses an infinite series of shorts cuts in the lower half plane at (−2g − in, +2g − in) for n ∈ N 0 . Similarly, the analytic continuationP a of P a across the first short cut is analytic in the lower half plane (LHPA) but possesses an infinite series of shorts cuts in the upper half plane at (−2g + in, +2g + in); see Fig. 3.
The sheet with the short branch cut corresponds to the following relation between the spectral plane and the rapidity plane:  Figure 3. For the purpose of solving the QSC for the Hagedorn temperature, a different Riemann sheet for P a is advantageous, obtained from the one depicted in Fig. 2 via analytic continuation. On this sheet, P a is an UHPA function with short branch cuts (−2g − in, +2g − in) for n ∈ N 0 , whileP a is a LHPA function with short branch cuts (−2g + in, +2g + in).
thus make the ansatz where c i,n = O(g 0 ) and starts to contribute at order 1/u n . 2 In particular, we have c 4,−1 = −1+O(g 2 ). The coefficients in this ansatz can be fixed by the gluing conditions (and asymptotic conditions) discussed in the next subsection, up to some remaining gauge symmetry discussed in the subsection thereafter. Note that this procedure is exactly opposite to the procedure for the QSC in the spectral problem, where an ansatz for P a is made.

Closing the equations
In the original formulation of the QSC for the spectral problem, the discontinuities of P a (u) and Q i (u) are determined via certain i-periodic functions µ ab (u) and ω ij (u) [11,12]. Our generalization to the Hagedorn temperature will instead be based on an alternative formulation [19,21,26,51], in which the analytic continuation across the cut is realized via complex conjugation. The corresponding equations are known as gluing conditions.
The gluing conditions connectP a (u) and P a (−u) on the Riemann sheet with short cuts depicted in Fig. 3. It turns out that for u ∈ (−2g, +2g), where complex conjugation is understood to change the sign of the i0 in Eq. (2.16). However, in contrast to the case of the spectral problem, the gluing conditions are not sufficient to fix all coefficients. In addition we should impose further restrictions on the large-u asymptotics of the P a functions that are not guaranteed by the ansatz (2.18) and the asymptotics of the Q a|i functions. In particular, for the perturbative analysis in Subsec. 3.1 we should impose the following behavior at u → ∞: Note that this follows from Eqs. (2.12) and (2.14) when imposing the gauge choice (2.15). One can easily modify this for other gauge choices. Similarly, one should also impose the correct asymptotics of the P a functions as an additional requirement in the numerical solution of Subsec. 3.2.

Gauge fixing
In the particular gauge (2.15) we chose, we have where ∼ denotes that the functions have the same coefficients in a large-u expansion while differing in there analytic structure. In particular, P 1 (u) and P 2 (u) are UPHA while P 3 (−u) and P 4 (−u) are LHPA. The above relations imply similar relations for Q a|i : for a = 1, 2. While this property is manifest for the leading coefficients via Eq. (2.21), it is a gauge choice that we can also impose on the higher-order coefficients to eliminate the gauge freedom at higher loop orders in the weak-coupling expansion as well as in the numeric approach at finite coupling.

Solving the Quantum Spectral Curve
Having derived the QSC equations for the Hagedorn temperature in the previous section, we now set out to solve it -first perturbatively at weak coupling in Subsec. 3.1 and then numerically at finite coupling in Subsec. 3.2. While some of the discussion provides further details of the procedure used in our letter [2], we also provide new perturbative results up to seven loops and discuss their rich number-theoretic structure, finding in particular that they can be expressed in terms of single-valued harmonic polylogarithms.

Solution algorithm
The starting point for the perturbative solution at weak coupling is the tree-level QSC presented in Subsec. 2.2 and App. F. We can determine the perturbative corrections to the tree-level QSC order by order in the coupling constant g using a slightly modified version of the general solution strategy of Ref. [19]. At each order ℓ, we start with the ansatz (2.18) for the Q i functions, which truncates at a finite order in 1/u and contains the undetermined coefficients T (ℓ) H as well as c H and similarly for c i,n (g 2 ). Combining Eqs. (2.6) and (2.8), we know that the exact solution Q a|i (u) to the QSC satisfies We define the mismatch in this equation when instead using the tree-level solution Q We can write the exact solution in terms of the tree-level solution Q Consequently, the matrix b a c satisfies the first-order difference equation At any loop order of perturbation theory at weak coupling, . with and where h is given while f is the unknown. Finite-difference equations of this type have previously been studied in the context of the twisted spectral problem [28] and of cusped Wilson lines [21]. In particular, the work [21] contains the Mathematica package TwistTools.m that implements some of the identities for solving these finite-difference equations. We derive further identities that are required in our case in App. A. Following the notation of Ref. [21], denote by Σ the operator generating f from h and z, i.e. f (u) = Σ(h(u), z). 3 To have a unique solution, we define Σ(0, z) = 0. Note that the action of Σ crucially depends on whether z = 1 or z = 1, e.g.
For h(u) given by a polynomial in u, Σ(h(u), z) is again a polynomial in u. Moreover, At subsequent loop orders, the single η functions η z a (u) can occur on the right-hand-side of the finite-difference equation (3.7). However, at any loop order at weak coupling, a solution to the twisted finite-difference equation (3.7) can be found in the space of rational functions in u multiplying the generalized multiple η functions [21]: These functions satisfy several important properties that we summaries in App. A; some of these have already been described in Ref. [21], while others are new. The finite-difference equation (3.7), also admits homogeneous solutions. For z = 1, these are given by constants, which we have to add to b with coefficients that have to be determined in the following steps. For general z, there is also the homogeneous solution It is removed by imposing P a to be UHPA. Note that in the twisted spectral problem [28], the variables z are pure phases, |z|= 1, such that z * = 1/z. In our case, however, both z = (2 + √ 3) 2 > 1 and z = 1/(2 + √ 3) 2 < 1 occur, and they are real. The functions above are understood to be defined by the given sum representations for values inside of their respective regions of convergence, and by analytic continuation outside of them; see also the discussions in Refs. [13,28].
We are now in the position to fix the various coefficients we introduced. 4 Using the orthogonality conditions (2.9), we find 3 Strictly speaking, the authors of Ref. [21] denote by f (u) = Σ(h(u)) the solution to the finite-difference where they allow for an exponential dependence in f and h. However, we have found it advantageous to factor out the exponential factor z −iu , cf. Eqs. (3.6)-(3.7), indicating this by Σ(·, z). 4 As can be seen from the discussion in Subsec. 2.2, already the tree-level coefficients contain √ 3. As discussed in Ref. [13] for the case of the spectral problem, it is thus more efficient to solve for the algebraic part of the higher-loop coefficients numerically with high accuracy, reconstructing the exact values using the PSLQ algorithm implemented e.g. in Mathematica as FindIntegerNullVector.
Imposing this fixes half of the coefficients in b stemming from the homogeneous solution.
Via Eq. (2.8), we then obtain P a from Q i and Q a|i . This allows us to impose the gluing conditions (2.19). As in Ref. [19], they enter by imposing that P a (u) +P a (u) and (P a (u)−P a (u))/ u 2 − 4g 2 are regular at u = 0. 5 In the process, the generalized η functions are evaluated at u = i, where they are proportional to multiple polylogarithms (MPLs). More precisely, For all z i = 1, the multiple polylogarithms reduce to multiple zeta values (MZV): The weight of these functions and numbers is defined as s 1 + · · · + s k . Note that MPLs have branch cuts; in particular, already classical polylogarithms Li s (z) have branch cuts in z on use the +i0 prescription (2.16). The occurrence of MPLs at this step already indicates that the Hagedorn temperature can be written in terms of these functions. Knowledge about the identities between the MPLs is crucial when solving for the undetermined coefficients, as solutions to the gluing conditions only exist when taking these identities into account. MPLs satisfy the so-called shuffle and stuffle relations, see e.g. Ref. [54] for a review. They can be used to reduce MPLs of weight less than four to classical polylogarithms. Moreover, the following inversion identity for classical polylogarithms exists for z / ∈ (0, 1]: where B n is the Bernoulli polynomial. We can numerically evaluate MPLs to arbitrary precision using GiNaC [55]. This can be used to find identities among MPLs using the PSLQ algorithms implemented e.g. in Mathematica as FindIntegerNullVector. As already indicated in Sec. 2, the gluing conditions do not suffice to fix all undetermined coefficients. In addition, we have to impose the correct asymptotics at large u via Eq. (2.20). This requires to expand η functions at large u; see App. A for details of this expansion. Similarly, an expansion is also required when imposing the final condition (2.22).

Perturbative results and their number-theoretic properties
Using the procedure described above, we have solved the Quantum Spectral Curve and determined the Hagedorn temperature up to and including seven-loop order. We now present these perturbative results and discuss their number-theoretic properties. We also attach the perturbative results in the ancillary file PerturbativeResults.m.
We write In general, we find that T . . is of mixed transcendentality, with a highest transcendental piece of transcendental degree 2ℓ−3. This is similar to the spectrum of anomalous dimensions in N = 4 SYM theory, which is also of mixed transcendentality, while for example ℓ-loop scattering amplitudes in N = 4 SYM theory are of uniform maximal transcendentality 2ℓ (see e.g. Ref. [54]).
As previously mentioned, the tree-level Hagedorn temperature is (3.19) in full agreement with Ref. [34]. At one-loop order, we find (3.20) in full agreement with Ref. [35]. At two-loop order, we have (3.21) in full agreement with Ref. [1]. At three-loop order, we find ≈ 37.22529358 . . . , (3.22) in full agreement with Ref. [2]. At four-loop order, we find While the Hagedorn temperature up to and including three-loop order is written in terms of classical polylogarithms, the four-loop Hagedorn temperature (3.23) contains the multiple polylogarithm Li 1,2 3) 2 is of weight three and can thus be expressed in terms of classical polylogarithms, this would introduce 1− 1 (2+ √ 3) 2 as argument, thus obscuring the origin of this number in the expansion of a generalized η function. Beyond four-loop order, however, it is not possible to write the Hagedorn temperature in terms of classical polylogarithms, as can be seen via the Lie cobracket test in Ref. [56]. 7 The higher-loop results for T H become increasingly lengthy when written in terms of MPLs, so we refrain from giving the full expressions here. Instead, we will now study the number-theoretic properties of these results, yielding much more compact expressions.
It follows from the discussion above that the arguments of the MPLs occurring in the Hagedorn temperature are necessarily z = (2+ √ 3) −2 to some power. Promoting the MPLs to functions of z in this way, we can use that MPLs depending only on a single variable can be expressed in terms of a smaller class of functions called harmonic polylogarithms [57], as discussed in the following. MPLs are equivalent to Goncharov polylogarithms [58][59][60] via 8 Harmonic polylogarithms (HPLs) [57] are defined as  9 We observe that the weak-coupling results can be expressed in terms of an even smaller class of functions, namely single-valued harmonic polylogarithms (SVHPLs). SVHPLs L were introduced by Francis Brown [39] as particular combinations of HPLs of complex conjugated arguments z, z * by requiring the branch cuts in the harmonic polylogarithms to cancel such that one obtains a single-valued function in the (z, z * ) plane. For example, Up to weight 6, they are explicitly given in auxiliary files of Ref. [62]. For higher weight, they can be generated using HyperLogProcedures [63], and we attach the relevant expressions up to weight 11 in the auxiliary file SVHPLreplacementsUpTo11.m. In order to express our result in terms of SVHPLs, we make an ansatz in terms of single-valued harmonic polylogarithms, reexpress them in terms of Goncharov polylogarithms and fix the coefficients by going to a fibration basis as implemented e.g. in HyperInt [61]. 10 We moreover define the shorthand notation In terms of SVHPLs, the perturbative results take a much more compact form. The first few orders read For space reasons, we refrain from showing the full six-loop and seven-loop results here; they can be found in the attached ancillary file PerturbativeResults.m. Their numeric values are T (6) H ≈ −49510.01767 , T H ≈ +625284.5652 . (3.33) While we have stopped at seven-loop order, there is no technical or conceptual obstacle to going to higher loop orders. It was previously observed for the Konishi anomalous dimension that the weak-coupling expansion can be expressed in terms of single-valued multiple zeta values (SVMZV) [13,15,40], which are SVHPLs evaluated at argument 1. Note that for cusped Wilson loops [21], weak-coupling results have only been obtained in the limit of small twist such that only classical polylogarithms occur; for these, the promotion to single-valued functions is trivial.
We observe that the perturbative results at weak coupling are palindromic in the arguments of the SVHPLs, i.e. they stay the same when sending L a 1 ,a 2 ,...,an → L an,...,a 2 ,a 1 . It would be interesting to understand the reason for this.
Moreover, the piece of T H with the highest transcendental degree seems to follow a simple structure. In the Konishi anomalous dimension, a corresponding piece could be determined in closed form for any loop order [40], and it would be interesting to develop a similar understanding here.

Numeric solution at finite coupling
In this section, we describe how to find the Hagedorn temperature T H numerically for finite values of the planar coupling g. We employ a modified version of the method of Ref. [16]; see also Refs. [51,64]. Moreover, we discuss the strong-coupling behaviour of the Hagedorn temperature.

Numeric algorithm
The key to solving the QSC numerically is to use the ansatz (2.18) for the four Q i (u) functions, but with the sum truncated so that n has a maximal value that we call K. This enables us to determine the four Q i (u) functions in terms of a finite number of parameters T H and c i,n with i = 1, 2, 3, 4 and n ≤ K. To find the values of these parameters for a given coupling g, we numerically and thus approximately solve the various conditions discussed in Sec. 2, similar to what we previously did at weak coupling.
For the 16 functions Q a|i (u), we make now the approximate ansatz which is truncated at a finite order depending on the value of N and where Notice that Q a|i (u) is exponentially convergent for N → ∞ when the imaginary part of s a u is large. Thus, this determines the region in which we can use Eq. (3.34) as a good approximation for sufficiently large N . The next step is to find the coefficients B a|i,n given the parameters T H and c i,n . A first step is to consider the leading large-u behavior of Eq. (2.6). On the LHS we use the asymptotic behavior of Eq. (3.34). On the RHS we use Eqs. (2.12) and (2.13) for the asymptotic behaviors of P a and Q i . Comparing the LHS and RHS, we find However, one cannot continue in this fashion since one needs to eliminate the P a (u) functions to find an approximate solution at large u. To this end, we notice that one can combine Eqs. (2.6), (2.8) and (2.10) to find in powers of 1/u for large u, giving the coefficients V a|i,n . In this way one finds equations that determine the B a|i,n coefficients. However, the procedure turns out to be slightly more complicated than one might naively expect. Expanding the LHS of Eq. (3.39), one finds that it goes like u 2 for large u, which is why we started the sum on the RHS with n = 1. The origin of this behavior is the last term on the LHS. This could seem surprising since we used above that in Eq. (2.6) all three terms starts at the same order. However, this is because in that case we imposed the asymptotics (2.12) on the P a (u) functions. Now, instead, we have eliminated the P a (u) functions and this means one cannot assume anymore that their asymptotics are of the form (2.12). Indeed, it is not difficult to show that with generic behavior of the B a|i,n and c i,n coefficients, the leading behavior of the P a (u) functions is for large u. Therefore, requiring the asymptotic behavior (2.12) imposes additional conditions on the B a|i,n and c i,n coefficients. 11 With this in mind, we now describe the procedure to solve the equations (3.40) at large u. One starts by choosing a numerical value for the coupling g, T H and the parameters c i,n , except for c 4,−1 and c 4,1 that are left as free variables.
At this point, all the equations V a|i,n = 0 are solved for n = 1, ..., 5. For n = 6, one gets a single independent equation which one can solve linearly for c 4,−1 . This in turn solves the equations for n = 7. Instead for n = 8 one finds a single independent equation which one can solve linearly for c 4,1 . Finally, for n = 9 one finds two equations V 1|1,9 = 0 and V 3|1,9 = 0 that one can solve linearly for B 1|1,1 and B 3|1,1 . With this, all equations V a|i,n = 0 are solved for n = 1, ..., 9. At n = 10 one has 16 equations V a|i,10 = 0 that one can solve linearly for the 16 variables B a|1,2 , B a|2,7 , B a|3,7 and B a|4,10 with a = 1, 2, 3, 4. At this point, one can look at the orthogonality relations (2.9) to obtain the two parameters B 2|1,1 and B 4|1,1 . Expanding for large u, and inserting the solutions of the variables found so far, one finds that its expansion starts at order u 3 . One can now solve for B 2|1,1 and B 4|1,1 by demanding that the u 3 terms are zero.
One proceeds now for n ≥ 11 as follows. For each successive n one has 16 equations V a|i,n = 0 that one can solve linearly for the 16 variables B a|1,n−8 , B a|2,n−3 , B a|3,n−3 and B a|4,n with a = 1, 2, 3, 4. In this way one can proceed order by order in n. To determine the approximate expression (3.34) for a given N , one needs to solve V a|i,n = 0 up to and including the order n = N + 8. Given this, one can check that the large-u expansion of Eq. (3.42) is zero to order u 5−N .
To recap, given numerical values for the coupling g, T H and the parameters c i,n (except for c 4,−1 and c 4,1 ) we have now found the approximate large-u Q a|i (3.34) to the desired accuracy given by N , and we have in addition obtained c 4,−1 and c 4,1 . One can now use this to determine the functions P a and its analytic continuationP a for a = 1, 2 on the real axis. Indeed, starting with a large and positive imaginary u, we can use Eq. (3.38) iteratively to find Q a|i closer and closer to the real axis. We choose the starting imaginary part of u, written here as U = Im(u), to be a sufficiently large and positive odd integer. We choose the real part of u such that after the iterative procedure one can find Q + a|i (u) for u ∈ I P , where I P is a suitably chosen set of P points in the real-valued interval (−2g, 2g). 12 This is possible for a = 1, 2 thanks to the above-mentioned exponential convergence. Building on this, one obtains P a (u) for u ∈ I P from the first line of Eq. (2.8) where Q i (u) is obtained from the truncated version of Eq. (2.18) as described above. One can furthermore find the analytic continuationP a (u) for u ∈ I P from P a = −Q i Q + a|i , (3.43) whereQ i is obtained from the truncated version of Eq. (2.18) by replacing x →x = 1/x. Having established a numerical procedure to compute P a andP a for a = 1, 2 and u ∈ (−2g, 2g), given values for the coupling g, T H and the parameters c i,n (except for c 4,−1 and c 4,1 ), we can construct a function F that parametrize how far we are from obeying the gluing conditions (2.19): (3.44) Given a value of the coupling g, one can now use a numerical minimization procedure based on F (T H , {c i,n }) that, given a suitable starting guess for T H and the relevant c i,n coefficients, can approach values of these parameters that minimize F (T H , {c i,n }), preferably so that it gets very close to zero. We have implemented this using the Levenberg-Marquardt algorithm, an improved version of Newton's method.

Results
We used our numerical procedure to capture two different regimes of the coupling g. We previously reported these results in Ref. [2]. In Fig. 4    The numerical data for T H (g) in the range 0 ≤ √ g ≤ 1.8 exhibits an approximately linear behavior where the uncertainties in c 0 and c 1 are in the last digits while the uncertainties in c 2 and c 3 are 0.0005 and 0.005, respectively. The coefficient c 0 was previously reported in Ref. [2]. The approximate leading linear behavior (3.45) is quite remarkable, as it agrees well with the Hagedorn temperature of type IIB string theory in flat space. Indeed, for large coupling g we expect that the part of the spectrum of the dual type IIB string theory on AdS 5 × S 5 that we probe corresponds to a flat space spectrum [2]. The Hagedorn temperature of type IIB string theory in flat space is 1/( √ 8πl s ) [33]. Using the AdS/CFT dictionary, this becomes Since 1/ √ 2π ≃ 0.39894 we find agreement with c 0 in Eq. (3.46). Thus, we have found a way to probe flat-space physics of ten-dimensional string theory within N = 4 SYM theory.
Moreover, since the first appearance of the present paper on the arXiv, the coefficients c 0 and c 1 have been analytically calculated by considering the inverse Hagedorn temperature as the radius for which a winding mode becomes massless [41,42]: Numerically 1/(2π) ≃ 0.1592, which thus fits with our numerical data (3.46). It would be very interesting to obtain further subleading coefficients as well.

Chemical potentials
In this section, we first discuss in Subsec. 4.1 the general relation between the Hagedorn temperature and the thermodynamic limit of Gibbs free energy per unit classical scaling dimension when turning on chemical potentials, generalizing Ref. [1]. Subsequently, in Subsec. 4.2, we use this to generalize the construction of the Quantum Spectral Curve presented in the previous sections to include non-zero chemical potentials. Finally, in Subsec. 4.3, we show how our integrability-based approach to the Hagedorn temperature is related to the Pólya-theory approach for zero 't Hooft coupling.

Hagedorn temperature and Gibbs free energy
We now generalize the relation found in Ref. [1] between the Hagedorn temperature and the Gibbs free energy per unit classical scaling dimension. This relation forms the basis for applying integrability-based methods to the calculation of the Hagedorn temperature. We write the full refined partition function of N = 4 SYM theory on R × S 3 as (4.1) Here, β = 1/T is the inverse temperature and D is the dilatation operator on flat Minkowski space R 1,3 , which is the image of the Hamiltonian of N = 4 SYM theory on R × S 3 under a conformal map. J 1 , J 2 , J 3 are the three su(4) R-charges and their chemical potentials are denoted by Ω 1 , Ω 2 and Ω 3 , respectively. Moreover, S 1 and S 2 are the two angular momenta on the S 3 , with corresponding chemical potentials Ω 4 and Ω 5 . Define the following charges associated to the psu(2, 2|4) spin chain: where with D 0 the classical scaling dimension and δD the anomalous scaling dimension. In terms of these charges, we have In the planar theory, all states can be written as products of single-trace states, which in turn can be mapped to spin chains. We can thus write the single-trace partition function as is the spin-chain free energy per unit classical scaling dimension for fixed D 0 = m 2 . A multi-trace state in the planar theory is given by a product of single-trace states, in which each bosonic single-trace factor can occur with multiplicity 0, 1, 2, 3, . . . and each fermionic single-trace factor can occur with multiplicity 0, 1. Moreover, the energy and other charges of the multi-trace state are given as sum of the contributions from the individual single-trace factors. The multi-trace partition function Z(T ) is thus given by where the alternating exploits the fact that single-trace states with even m = 2D 0 are bosons, while those with odd m = 2D 0 are fermions. The Hagedorn temperature T H is the lowest temperature above which the planar partition function diverges. Since F m is a monotonically decreasing function of the temperature, a divergence first occurs when the n = 1 contribution to the multi-trace partition diverges. The n = 1 contribution is (4.8) Define the thermodynamic limit of the Gibbs free energy per unit classical scaling dimension of the psu(2, 2|4) spin chain: In terms of this, one sees from Eq. (4.8) using the Cauchy root test that the Hagedorn temperature is determined as [1] F (T H , Ω i ) = −1 + Ω 1 , (4.10) since we have exp(− 1 2 β(1 − Ω 1 + F (T, Ω i ))) > 1 when T > T H . One can now in principle find T H for any chemical potentials Ω i and any coupling g by computing F (T, Ω i ). As described in App. B.1, this can be done by solving the TBA equations (B.2)-(B.6) with boundary conditions (B.11) to obtain the Y-functions Y a,s (u) and from this computing F (T, Ω i ) by Eq. (B.12). To simplify this, one can reformulate the TBA equations as Y-system equations, see App. B.2 as used in Ref. [1] for the case of zero chemical potentials. Finally, one can alternatively read off F (T, Ω i ) from Eq. (B.25), using the asymptotic behavior of the Y-functions Y a,s (u). However, as is clear from Sec. 2 and 3, a considerably more powerful and efficient approach to find F (T, Ω i ) is through a set of QSC equations.

Quantum Spectral Curve
We now consider the generalization of the QSC for the case of non-zero chemical potentials. A main part of the QSC does in fact remain unchanged, namely the general structure of the QSC equations and their branch cuts as described in Subsec. 2.1 and 2.3.
The main difference compared to the case without chemical potentials lies in the asymptotics of the P a and Q i functions. To find these asymptotics, we can consider the asymptotic (and thus constant) solution to the T-system. As discussed in App. C.2, it is given by the psu(2, 2|4) character solution of Ref. [65] upon identifying the parameters x a and y i of Ref. [65] with T H and the chemical potentials Ω i as Fortunately for us, the asymptotics of the P a and Q i functions that reproduce the psu(2, 2|4) character solution of Ref. [65] have already been identified in the context of the twisted spectral problem [28]. However, we have to make the interchange P a ↔ Q i together with x a ↔ y i and A a ↔ B i compared to Ref. [28]. The transformation P a ↔ Q i , x a ↔ y i and A a ↔ B i is due to the fact that we consider the direct physical theory rather than the mirror theory, as discussed in the Introduction (Sec. 1) as well as in Subsec. 2.3 and App. B.1. 13 This gives where (no sum over a) , (4.13) with (4.15) The asymptotics above change discontinuously when certain chemical potentials vanish or become equal to each other, resulting in certain x a and y i becoming equal. In the generic case of all chemical potentials being nonvanishing and unequal, the equivalence between the asymptotic T-system and the asymptotic above is illustrated in App. E. Left-right symmetry occurs in the case Ω 3 = Ω 5 = 0. Finally, it is easy to see that the asymptotics above reproduce Eqs. (2.12)-(2.14) in the case where all chemical potentials vanish.
We leave the determination of the gluing conditions, which close the equations, for future work. We now turn to the special case of the free theory, where the branch cuts vanish and the gluing conditions are thus not required. 14

Zero-coupling limit and single-particle partition functions
Previously in the literature, the Hagedorn temperature for N = 4 SYM theory with nonzero chemical potentials has only been computed via Pólya theory [44][45][46][47][48][49] based on the methods introduced for the case of vanishing chemical potentials in Refs. [34,38] at tree level and in Ref. [35] at one-loop order. Below we check our tree-level results that emerge from the TBA analysis with the previously derived results of Refs. [44,45,47] using Pólya theory. This is an important consistency check on our methods.
At zero coupling g = 0, the full refined partition function (4.1) can be computed using Pólya theory [34,38]: where η B (T, Ω i ) and η F (T, Ω i ) are the single-particle partition functions for the bosonic and fermionic modes on the three-sphere, computed for N = 4 SYM theory with chemical 14 Note that the zero-coupling solution is given by exponential factors times polynomials; no (generalized) η functions occur, in contrast to the spectral problem.
potentials in Refs. [44,45,47]. Defining the total single-particle partition function as Thus, at zero coupling g = 0 we have two rather explicit sets of methods to compute T H . Either we determine it from the single-particle partition function via Eq. (4.18). Or, we use the integrability methods laid out in this paper. Since one cannot solve Eq. (4.18) explicitly for general chemical potentials Ω i , a pertinent question is: can we reproduce Eq. (4.18) from the methods of this paper? The answer is affirmative, and a simple route to this is to use the T-system instead of the QSC. At zero coupling g = 0, the T-system is known to be constant [1], and hence the general asymptotic T-system reviewed in App. C.1 and C.2 should hold for all u, i.e. T a,s (u) = T ∞ a,s . As shown in App. D, one can derive from the requirement of a constant Y-system that for g = 0. In App. D, we compute T 1,0 and find that we can identify Hence, the two methods are equivalent for g = 0.

Deformations
The maximally supersymmetric Yang-Mills theory admits several deformations that preserve integrability but break some or all of supersymmetry and or Poincaré symmetry. The best-understood class of these deformations contains so-called diagonal twists and was studied at the level of the Bethe equations in Ref. [66]. This class encompasses the N = 1 supersymmetric one-parameter real-β deformation, which is a special case of the N = 1 supersymmetric Leigh-Strassler deformations [67], as well as the non-supersymmetric threeparameter γ i deformation [68]. These deformations change the interaction vertices of the theory -leaving the tree-level partition function and tree-level Hagedorn temperature trivially the same. In Ref. [43], the one-loop corrections to these quantities were calculated in the real-β as well as γ i deformation. It was found that while the one-loop partition function depends on the deformation parameters, the one-loop Hagedorn temperature is the same as in N = 4 SYM theory. 15 We will now argue that the Hagedorn temperature of these deformations agrees with that of N = 4 SYM theory for any value of the coupling. At the heuristic level, we mentioned already in the Introduction that the deformations are encoded in twisted boundary conditions along the direction that becomes infinite in the thermodynamic limit, and thus they do not affect the final result.
This can equally be seen at the technical level. Deriving TBA equations for the Hagedorn temperature of the deformed theories with diagonal twists following Ref. [1], the starting point are the asymptotic Bethe equations of Ref. [66]. These contain the deformation parameters as a separate factor. Upon taking the logarithm and then the derivative, this factor drops out completely. This shows that the Hagedorn temperature in the deformations with diagonal twist, and in particular in the real-β and γ i deformation, is the same as in N = 4 SYM theory at all orders in the 't Hooft coupling.
A related deformation of N = 4 SYM theory is the integrable, conformal planar fishnet theory [72][73][74]. It arises by taking a double-scaling limit of the γ i deformation, taking γ 1 → i∞, λ → 0 withĝ = e −iγ 1 λ fixed. In this limit, all fields except two complex scalars decouple. Naively taking the same double-scaling limit of our result for the Hagedorn temperature would result in T H (λ = 0), since the Hagedorn temperature of the γ i deformation is independent of the deformation parameters, as discussed above. However, the Hagedorn temperature depends on all fields in the theory, whether they decouple or not, and the conformal fishnet theory is conventionally defined without the decoupled fields. Using the single-particle partition function of the conformal fishnet theory, 4 , the tree-level Hagedorn temperature is given as the solution to the equation η(e −1/T tree,fishnet H ) = 1, T tree,fishnet H = 0.508028 . . . ; it clearly differs from the result in N = 4 SYM theory. It would be very interesting to apply our method to calculate the Hagedorn temperature of the conformal fishnet theory at any value of the respective couplingĝ, in analogy to what we have done in the present paper for N = 4 SYM theory. 16

Conclusion and outlook
In this paper, we have derived a Quantum Spectral Curve for the Hagedorn temperature, providing several details deferred in our letters [1,2].
The Hagedorn QSC can be efficiently solved perturbatively at weak coupling as well as numerically at finite coupling. We have extended our previous perturbative results [1,2] up to and including seven-loop order; we have attached these results in the ancillary file PerturbativeResults.m. Our perturbative results show interesting number-theoretic properties, namely being expressible in terms of single-valued harmonic polylogarithms (SVHPLs). This is similar to the situation for the spectrum of local composite operators, which is expressible in terms of single-valued multiple zeta values (SVMZVs) [13,15,40], which are SVHPLs evaluated at 1. In our case, however, the SVHPLs are evaluated at At the technical level, the Hagedorn temperature enters the QSC as a twists. We expect the tools used here, as well as the number-theoretic observations, to be useful also for the systematic solution of the QSC in other cases with twists, such as deformations of N = 4 SYM theory [28][29][30][31], cusped Wilson loops [21][22][23][24] and color-twisted operators in N = 4 SYM theory [27], which were recently investigated in the context of structure constants.
Finally, we generalized the Hagedorn QSC to include also chemical potentials, as well as to a class of integrable deformations of N = 4 SYM theory. In the letter case, we found that the Hagedorn temperature is identical to the one in N = 4 SYM theory for any value of the coupling.
Let us end discussing a number of interesting future directions. From our numeric solution, we have extracted the leading and subleading behavior of the Hagedorn temperature at large λ. It would be interesting to extract also further subleading orders, or to develop a systematic approach to solving the Hagedorn QSC perturbatively at strong coupling. It would also be interesting to solve the Hagedorn QSC in the presence of chemical potentials.
The integrability-based approach to the Hagedorn temperature should also be applicable to further integrable theories. These include deformations of N = 4 SYM theory that are not given by diagonal twists, such as the η deformation, for which a TBA and QSC has been developed in Refs. [29,77]. Further examples of integrable theories to which our approach should be applicable are the ABJM and ABJ theory in the context of AdS 4 /CFT 3 [78,79], for which a QSC has equally been studied [80][81][82][83][84][85]. Additional examples occur in the context of AdS 3 /CFT 2 , for which a QSC has recently been proposed [86,87]. In particular, the theory considered in Ref. [88] is free of wrapping corrections, making it an interesting starting point to investigate an integrability-based approach to the full partition function.
Moreover, it has been pointed out in Ref. [89] that four-point functions of determinant operators exhibit a critical behavior that bears resemblance to Hagedorn behavior and it would be interesting to calculate the critical configuration via a similar integrability based approach as applied here for the Hagedorn temperature.
Similar techniques to the ones employed for the Hagedorn temperature could also be used to calculate critical exponents, which describe how exactly the partition function diverges when approaching the Hagedorn temperature In this paper, we have calculated the Hagedorn temperature, which plays the role of a limiting temperature, signaling a phase transition. The low-energy phase ceases to exist at the Hagedorn temperature. The confinement-deconfinement phase transition, dual to the Hawking-Page transition, is expected to occur at a lower temperature [38]. It has recently been shown that the Hagedorn behavior at infinite N is replaced by Lee-Yang behavior at large but finite N [90]. It would be extremely interesting to develop an integrability-based approach also for the confinement-deconfinement temperature.
grants number DFF-4002-00037, by the ERC advance grant 291092, by the ERC starting grant 757978 and by the research grants 00015369 and 00025445 from Villum Fonden.

A Generalized η functions
In this appendix, we discuss two important properties of the generalized η functions (3.11): Many properties of the generalized η functions were already discussed in Ref. [21], such as their (shuffle) product and their behavior under shifts; we refer the reader to App. F of Ref. [21] for details. Moreover, these relations are conveniently implemented in the Mathematica package TwistTools.m accompanying Ref. [21]. We will mainly need two additional properties.

Relations for η with vanishing arguments
Generalized η functions can be further simplified in the case that one of their indices vanishes, s i = 0. This simplifications requires us to distinguish whether the corresponding z i = 1 or z i = 1; the case z i = 1 was already worked out in the Mathematica package TwistTools.m accompanying Ref. [21].
Let us first consider z i = 1. We have from which it follows that Using partial fractioning, we find and a similar expression for the index i + 1. Thus, An important special case is s 1 = · · · = s k = 0. In particular, Note that the sum in η z 0 is divergent for |z|> 1; the right hand side is the appropriate analytic continuation.

Expansion at infinity
In the perturbative algorithm, it is required to expand the η functions around u = ∞. Consider an individual term in Eq. (3.11). We have Thus, (A.9) In this expansion, we in particular encounter the formally divergent quantities Li −n (1) for n ∈ N 0 , which we regularize as Li −n (1) = ζ −n .

B TBA equations and Y-system
In this appendix, we provide further details on the TBA and Y-system equations for the Hagedorn temperature, which we deferred in our letter [1].

B.1 TBA equations
In the following, we review the TBA equations for the psu(2, 2|4) spin chain of N = 4 SYM theory at finite temperature and in the presence of chemical potentials. These are obtained in complete analogy with the TBA equations of the spectral problem for N = 4 SYM theory [5][6][7][8][9][10]. The only subtle difference is that in our case we use the direct physical theory, based on the direct association between temperature of N = 4 SYM and the temperature of the spin chain, whereas one uses instead the so-called mirror theory in the spectral problem, since in that case one should perform a double Wick rotation as explained in the Introduction.
To obtain TBA equations for the psu(2, 2|4) spin chain of N = 4 SYM theory, one starts with the asymptotic Bethe ansatz of Refs. [91,92]. One then assumes the string hypothesis, which is a specific assumption to what constitutes the complete set of solutions to the asymptotic Bethe ansatz equations. In particular, these solutions are organized in terms of so-called strings. One considers an ensemble of configurations of strings of various types and lengths in the continuum limit by sending the length of the spin chains to infinity. This is described by densities of values of the centers of the strings realized by the particular configuration of strings in the ensemble as well as the density of center values not realized. From this one can define the entropy per unit classical scaling dimension s, the energy per unit classical scaling dimension e, and the charges per unit classical scaling dimensioñ q (i) with i = 1, ..., 5. Imposing the first law of thermodynamics δe = T δs + i Ω i δq (i) at temperature T and for chemical potentials Ω i , one obtains the TBA equations in terms of Y-functions and an associated free energy per unit classical scaling dimension. The Yfunctions are defined as the ratio of the density of values not realized over the density of values that are realized by the configurations in the ensemble. Following the notation of Ref. [1], one can define Y-functions Y a,s (u) with (a, s) in the set M , They form a so-called hook, illustrated in Fig. 6.
In terms of the Y-functions Y a,s (u), the TBA equations for the direct theory take the following form: where n ≥ 1 and u ∈ R \ (−2g, 2g) using the kernels ǫ, K, Σ, Θ and a defined in App. B.4. We define here the convolutions for the functions f (u) and g(v, u). The above TBA equations for the direct theory were previously considered in Refs. [10,77] but in different thermodynamic limits. Note that we have assumed Concretely, this is important to get the right asymptotics for large n. 17 One finds the large n asymptotics (B.11) 17 Note that one can always reparametrize the given charges to obtain this.
The Gibbs free energy per unit classical scaling dimension is given by

B.2 Y-system
While the TBA equations in principle solve the problem of determining the Gibbs free energy per unit classical scaling dimension F (T, Ω i ) -and thus the Hagedorn temperature -in practice they are quite difficult to work with. A first step towards a simplification is to recast them in terms of so-called Y-system equations, as done for the case of the spectral problem in Refs. [6][7][8][9][10]. In case of zero chemical potentials, we used this in Ref. [1] to determine the Hagedorn temperature up to and including order g 4 (two loops). Analytically extending the Y-functions Y a,s (u), one finds from the analytic properties of the TBA equations (B.2)-(B.6) that they are analytic in the strip with Im(u) < 1 2 |a−|s||. One can now derive the Y-system equations with a n (u) defined in Eq. (B.26). For the exception (a, s) = (1, 0), we find with the source term ρ(u) given by This is important primarily with respect to reformulating the TBA equations in terms of a QSC [2] but it is also useful for solving the Y-system at zero coupling g = 0 [1]. As we describe in App. C, the asymptotic Y-system Y ∞ a,s can be deduced from a constant Tsystem, and this in turn provides the seed for all-important boundary conditions for the asymptotic behavior of the P a (u) and Q i (u) in the QSC, as described in Subsec. 2.2 and App. E. A crucial piece in this is how we can infer the Gibbs free energy per unit classical scaling dimension F (T, Ω i ) from the asymptotic behavior of Y-functions in the Y-system, since this in turns makes it possible to find it from the asymptotics of the T-functions, and therefore also from the asymptotics of the P a (u) and Q i (u) functions in the QSC. Below we provide this piece, by showing that one can infer the Gibbs free energy per unit classical scaling dimension F (T, Ω i ) directly from the asymptotic Y-system Y ∞ a,s . Indeed, this provides the argument behind Eq. (C.31), which we used in Ref. [2].
From Eq. (B.38), we have . (B.20) We consider u ∈ R\(−2g, 2g). For large u, we have x(u) = u+O (1/u). Considering the first term, one observes that since g 2 /x(u) goes to zero for u → ∞; thus, one obtains −θ n (v) as one can infer from the definition (B.27). The second term instead goes to zero if v is kept fixed for u → ∞. However, if v − u is finite, it is non-zero for u → ∞. Since this requires v → ∞, one has x(v + i 2 n) ≃ v; hence, the second term gives a n (v − u), as one can infer from the definition (B.26). Thus, we have derived For the second term, one observes that one picks up only contributions for large v in the integral since u → ∞. Hence, A slight rewriting of the Gibbs free energy per unit classical scaling dimension (B.12) reveals Therefore, combining this with Eq. (B.23) one finds that

B.4 Definitions of functions and kernels
In this appendix, we have collected the various definitions for families of functions and kernels that are needed for formulating the TBA equations (B.2)-(B.6). We define for a positive integer n and u ∈ C the functions a n (u) = i 2π We extend these three families of functions to n = 0 by taking the limit n → 0 + : We turn now to the kernels. For n = m, we define the kernel K mn (v, u) via For n = m, we define Moreover, we define the kernels where we defined We also define the kernels , (B.42) which is defined in terms of the dressing factor [92]

C Asymptotic T-system
In this appendix, we determine the asymptotic values of the T-system functions, from which one can infer the the asymptotic values of the Y-functions Y ∞ a,s defined by Eq. (B.19). The purpose is to generalize the asymptotic T-system found previously in Ref. [1] for zero chemical potentials to the case of non-zero chemical potentials. This is crucial for identifying the correct asymptotic behavior of the P a (u) and Q i (u) functions in the QSC with generic chemical potentials in App. E and furthermore to the case of general chemical potentials presented in Subsec. 4.2.

C.1 T-system
To set up the framework for our analysis of the asymptotic values of the T-system functions, we review first very briefly what a T-system is.
One can translate the The T-functions should obey the Hirota equations In addition to this, the T-functions should obey certain analyticity properties listed in Ref. [93]. Note that there are certain gauge freedoms of the T-system functions T a,s (u) relating different T-systems that correspond to the same Y-system; see for instance Ref. [93]. We impose the following gauge conditions on the T-functions: T 2,n = T n,2 and T 2,−n = T n,−2 for n ≥ 2 . (C.4)

C.2 General asymptotic T-system
We turn now to the asymptotic values of the T-functions T a,s (u). We define these as The for (a, s) ∈M . We inherit the gauge choice (C.4) and impose the additional gauge condition T ∞ 0,s = 1 for s ∈ Z. This is possible to impose since it involves only a constant transformation of the T-system, which is consistent with Eq. (C.4). Considering now the formula (B.25) for the Gibbs free energy per unit classical scaling dimension in terms of Y ∞ a,s , we can translate this into a relation between the free energy and the asymptotic T-system: Using Eq. (4.10), this reveals a direct connection between the Hagedorn temperature T H and the asymptotic T-system. The general solution to the constant Hirota equations (C.7) is the psu(2, 2|4) character solution of Ref. [65]. The solution is presented in terms of the eight variables x 1 , x 2 , x 3 , x 4 , y 1 , y 2 , y 3 and y 4 , where one can think of x 1 , x 2 , x 3 , x 4 as associated with the su(2, 2) subalgebra and y 1 , y 2 , y 3 , y 4 , as associated with the su(4) subalgebra. For a ≥ |s|, the solution is From the above, one finds T ∞ a,s for s ≥ 0. For negative s, one can find T ∞ a,s via the general relation (C.13) Note that one has manifestly T ∞ 0,s = 1 in the above solution. The gauge choice (C.4) corresponds to x 1 x 2 x 3 x 4 = y 1 y 2 y 3 y 4 . (C.14) One can check that making the rescaling x a → αx a and y i → αy i results in T ∞ a,s → α as T ∞ a,s which corresponds to a gauge transformation of the asymptotic T-system. Thus, one can freely make a rescaling of this kind to impose y 1 y 2 y 3 y 4 = 1, thus x 1 x 2 x 3 x 4 = y 1 y 2 y 3 y 4 = 1 .
We impose the following ordering in the case of non-zero chemical potentials: The x i 's are negative due to the boundary conditions of the fermions in the partition function. Indeed, from (E.6) we see that negativity of the x i 's is in accordance with the signs in the (−e −1/(2T H ) ) ±iu factors in Eq. (2.12). The inequalities (C. 16) assume This means the above solution is not fully general as it does not work for the full range of chemical potentials constrained by the bounds (B.10) that one can always assume without loss of generality. However, one can describe the cases in which some or all of the more general bounds (B.10) are saturated by taking limits of the above solution. For instance, as we review in App. C.4, one finds the case of zero chemical potentials by a limit of the above [1]. Thus, in this sense, we can think of the above solution for the asymptotic T-system as fully general.
Combining the above solution for T ∞ a,s with the Y-system asymptotics (B.11), we find From this, we find the y i 's in terms of the chemical potentials and temperature: Inserting now the relation to the Hagedorn temperature (4.10), one finds Eq. (4.11). Once inserted into the above solution for T ∞ a,s , this provides the asymptotic T-system for a Hagedorn temperature T H with chemical potentials obeying the bounds (C. 17). From this, one can approach the cases in which some or all of the more general bounds (B.10) are saturated by taking limits of this solution.

C.3 Explicit form of the asymptotic T-system
For connecting to the Q-system, and hence the QSC, it is useful to rewrite the general solution (C.9)-(C.13) for the asymptotic T-system T ∞ a,s in a more explicit form. We explicitly present the generic case (C.17) here. The cases in which some of the inequalities (B.10) are saturated can be obtained by taking limits of the expressions below, as illustrated for the case of vanishing chemical potentials in the subsequent subsection. We make the following definitions: 18 , For the right band s ≥ a ≥ 0, we record For the left band s ≤ −a ≤ 0, we record Finally, for the upper band a ≥ |s|, we record In addition, one finds T ∞ a,0 from the general relation . (C.28) One can check that the above formulas give the same result as Eqs. (C.9)-(C.13). The above explicit form (C.23)-(C.27) for T ∞ a,s is useful in that one can directly infer the dependence on a and s. This dependence is instead hidden in the character solution (C.9)-(C.13). We use this explicit form below in App. E to connect to the large-u behavior of the QSC.

C.4 Asymptotic T-system for zero chemical potentials
One can find the asymptotic T-system for zero chemical potentials by first setting Ω 3 = Ω 5 = 0 and subsequently taking the limits Ω 1 → 0, Ω 2 → 0 and Ω 4 → 0 of the solution (C.23)-(C.27) of App. C.3. One finds for a ≥ |s|, and for |s|≥ a. This is the T-system given in our letter [1], where it was written in terms of the parameter Imposing the Hagedorn temperature condition (4.10) with zero chemical potentials further fixes F (T H ) = −1.

D Zeroth order T-system and the Hagedorn temperature
In this appendix, we provide details on the relation between the integrability-based method and the Pólya-theory method, which we mentioned in Subsec. 4.3. For vanishing 't Hooft coupling λ = 0, the Y-system is constant. Thus, it is equal to the asymptotic Y-system Y a,s | λ=0 = Y ∞ a,s | T H =T H (λ=0) . Similarly, for the T-system For vanishing 't Hooft coupling, one can see from Eq. (B.15) along with the definition of the kernels in App. B.4 that a constant Y-system implies Y 1,1 Y 2,2 = Y 1,−1 Y 2,−2 = 1. Using Eq. (C.2) together with the gauge (C.4) as well as T ∞ 0,1 = 1, we derive the condition Imposing Eq. (C.15), we compute from the general asymptotic T-system (C.9)-(C.10) where we defined η(x, y) = η s (x, y) + η v (x, y) + η f 1 (x, y) + η f 2 (x, y) , (D.4) with η s (x, y) = (x 3 x 4 − x 1 x 2 )(y 1 y 2 + y 3 y 4 + y 1 y 3 + y 2 y 4 + y 2 y 3 + y 1 y 4 ) ( . (D.8) Noting that Eq. one can easily check that η(x, y) = η(T H , Ω i ) is the single-particle partition function for N = 4 SYM theory at zero 't Hooft coupling, with η s originating from the scalar fields, η v from the gauge boson and η f 1 as well as η f 2 from the fermions. As derived in Refs. [34,38], the condition for the Hagedorn temperature T H at zero 't Hooft coupling is η(T H , Ω i ) = 1. This corresponds precisely to the requirement (D.2). Therefore, we have shown that the Hagedorn temperature at zero 't Hooft coupling for any value of the chemical potentials correspond to the one obtained using the single-particle partition function for N = 4 SYM theory in Refs. [34,38].

E Asymptotic Q-system with generic chemical potentials
In this appendix, we use the results for the asymptotic T-system T ∞ a,s recorded in App. C.2 and C.3 to obtain the large-u asymptotic behavior of the Q-system for the case of generic chemical potentials obeying the bounds (C. 17). From this, one can infer the asymptotic behavior of Q i (u) and P a (u). In Subsec. 4.2 this is extended to the fully general case. Note that this treatment is a special case of the one considered in Ref. [28], as discussed in Subsec. 4.2.
Consider the asymptotic T-system T ∞ a,s appropriate for obtaining the Hagedorn temperature. This is given by Eqs. (C.23)-(C.28) with x i and y i given by Eq. (4.11). This T-system is associated with a Q-system through the relation (4.11) of Ref. [94]. For the left and right bands, we have  for a, j = 1, 2, 3, 4. For the constants A a , A a , B j and B j we impose the constraints where A a and B i are defined by Eq. (C.22) and where there are no sums over a and i on the left-hand sides of the two expressions. As mentioned in Subsec. 4.2, this ansatz is related to a special case of the large-u asymptotics of the twisted QSC in Ref. [28]. From Eqs. (E.6) and (E.7), the Q-system can now be inferred from the QQ-relations (2.1)-(2.3) as well as the definitions of the Hodge dual Q-functions (2.7). In particular, we record for a ≥ 1. Combining these relations with Eq. (E.10), we reproduce the general asymptotic T-system (C.25)-(C.27) for the upper band. It is not necessary to check that T ∞ a,0 matches, since that is determined by the relation (C.28).

F Asymptotic Q a|i for vanishing chemical potentials
In this appendix, we give the functions Q a|i that correspond to the asymptotic T-system in the case of vanishing chemical potentials; we deferred these in Subsec. 2.2. In particular, these Q a|i functions give Q We have