Black holes in 4d $\mathcal{N}=4$ Super-Yang-Mills

We resolve a long-standing question: does the four-dimensional $\mathcal{N}=4$ SU(N) Super-Yang-Mills theory on $S^3$ at large N contain enough states to account for the entropy of rotating electrically-charged BPS black holes in AdS$_5$? Our answer is positive. We reconsider the large N limit of the superconformal index, using the Bethe Ansatz formulation, and find an exponentially large contribution which exactly reproduces the Bekenstein-Hawking entropy of the black holes of Gutowski-Reall. Besides, the large N limit exhibits a complicated structure, with many competing exponential contributions and Stokes lines, hinting at new physics.


Introduction and results
One of the fascinating aspects of black hole physics is its connection with the laws of thermodynamics. Of particular importance is the fact that black holes carry a macroscopic entropy [1][2][3][4][5], classically determined in terms of the horizon area. In the search for a theory of quantum gravity, explaining the microscopic origin of black hole thermodynamics is a fundamental but challenging test.
String theory is proposed to embed gravity in a consistent quantum system, hence it should in particular explain the black hole entropy in terms of a degeneracy of string states. This has been beautifully shown to be the case in the seminal paper [6] by Strominger and Vafa, where the Bekenstein-Hawking entropy of a class of supersymmetric asymptotically-flat black holes was microscopically reproduced.
In the case of asymptotically-AdS black holes, the AdS/CFT duality [7][8][9] constitutes a natural and wonderful framework to study their properties at the quantum level. The duality provides a non-perturbative definition of quantum gravity, in terms of a conformal field theory (CFT) living at the boundary of AdS space. The problem of offering a microscopic account of the black hole entropy is rephrased into that of counting particular states in the dual CFT. However, despite the very favorable setup, this problem in four or more dimensions has remained unsolved for many years, and only recently a concrete example was successfully studied in [10,11]. There, the Bekenstein-Hawking entropy of a class of static dyonic BPS black holes in AdS 4 was holographically reproduced in the dual CFT 3 , via supersymmetric localization [12][13][14]. Since then, the matching has been extended to many other classes of magnetically-charged BPS black holes in various dimensions [15][16][17][18][19][20][21][22][23][24][25][26][27], including the first quantum corrections [28][29][30][31][32][33].
When moving to rotating, purely electric black holes, the situation becomes more complicated. Famously, the microstate counting for BPS black holes in AdS 5 has remained a long-standing open problem, which dates back to the work of [34,35]. In this context, BPS black holes arise as rotating electrically-charged solutions of type IIB string theory on AdS 5 × S 5 [36][37][38][39][40]. Their holographic description is in terms of 1/16 BPS states of the boundary 4d N = 4 Super-Yang-Mills (SYM) theory on S 3 , which can be counted (with sign) by the superconformal index [41,35]. One would expect the contribution of the black hole microstates to the index to dominate the large N (i.e. weak curvature) expansion. However, the large N computation of the index performed in [35] showed no rapid enough growth of the number of states, and thus it could not reproduce the entropy of the dual black holes. Additionally, that result was followed by several studies of BPS operators at weak coupling [42][43][44][45][46] in which no sign of high degeneracy of states was found.
Very recently, the issue received renewed attention leading towards a different conclusion. First, the authors of [47] related the black hole entropy to the (complexified) regularized onshell action of the gravitational black hole solutions, and then compared the latter with the S 3 × S 1 supersymmetric partition function of the field theory, finding perfect agreement at leading order in large N . Second, the authors of [48] analyzed the index in a doublescaling Cardy-like limit, finding quantitative evidence that the index does account for the entropy of large BPS black holes (whose size is much larger than the AdS radius). Third, in [49] it was observed that, even at finite values of the fugacities, the index exhibits a deconfinement transition before the Hawking-Page transition related to the known AdS 5 black holes, pointing towards the existence of hairy black holes.
In this paper we offer a resolution of the issue by revisiting the counting of 1/16 BPS states in the boundary N = 4 SYM theory at large N . We approach the problem by using a new expression for the superconformal index of the theory, derived in [50,51] and dubbed Bethe Ansatz (BA) formula, which allows for an easier analysis of the large N limit. We find that the superconformal index, i.e. the grand canonical partition function of 1/16 BPS states, does in fact grow very rapidly with N -as e O(N 2 ) -for generic complex values of the fugacities. Although the BA formula of [51] can handle the general case, this is technically difficult and in this paper we restrict to states and black holes with two equal angular momenta, as in [37].
The BA formulation reveals that the large N limit has a complicated structure. There are many exponentially large contributions, that somehow play the role of saddle points. As we vary the complex fugacities, those contributions compete and in different regions of the fugacity space, different contributions dominate. This gives rise to Stokes lines, separating different domains of analyticity of the limit. The presence of Stokes lines also resolves the apparent tension with the computation of [35], that was performed with real fugacities. We show that when the fugacities are taken to be real, all exponentially large contributions organize into competing pairs that can conceivably cancel against each other. The fact that for real fugacities the index suffers from strong and non-generic cancelations was already stressed in [48,49].
Our main result is to identify a particular exponential contribution, such that extracting from it the microcanonical degeneracy of states exactly reproduces the Bekenstein-Hawking entropy of BPS black holes in AdS 5 (whose Legendre transform was obtained in [52]). This is in line with the double-scaling Cardy-like limit of [48]. Along the way, we show that the very same I-extremization principle [10,11] found in AdS 4 , is also at work in AdS 5 guaranteeing that the index captures the total number of single-center BPS black hole states.
At the same time, we step into many other exponentially large contributions: we expect them to describe very interesting new physics, that we urge to uncover. To that purpose, we study in greater detail the case of BPS black holes with equal charges and angular momenta [36]. We find that while for large black holes their entropy dominates the superconformal index, this is not so for smaller black holes. This seems to suggest 1 that an instability, possibly towards hairy or multi-center black holes, might develop as the charges are decreased. Similar observations were made in [48,49]. It would be extremely interesting if there were some connections with the recent works [53][54][55][56], and we leave this issue for future investigations.
The paper is organized as follows. In Section 2 we review the charges and entropy of BPS black holes in AdS 5 . In Section 3 we present the BA formula for the superconformal index of N = 4 SYM, and in Section 4 we compute its large N limit. Sections 5 and 6 are devoted to extracting the black hole entropy from the index.

BPS black holes in AdS 5
In this paper we study the entropy of rotating charged BPS black holes in AdS 5 [36][37][38][39][40] that can be embedded in type IIB string theory on AdS 5 ×S 5 [57]. In order to set the stage, let us briefly review such gravitational solutions. The black holes are solutions to the equations of motion of type IIB supergravity that preserve one complex supercharge [58], thus being 1/16 BPS. The metric interpolates between the AdS 5 boundary and a fibration of AdS 2 on S 3 at the horizon. Moreover, the black holes carry three charges Q 1,2,3 for U (1) 3 ⊂ SO(6) acting on S 5 , that appear as electric charges in AdS 5 , and two angular momenta J 1,2 associated to the Cartan U (1) 2 ⊂ SO(4) (each Cartan generator acts on an R 2 plane inside R 4 ). The black hole mass is fixed by the linear BPS constraint where g = −1 5 is the gauge coupling, determined in terms of the curvature radius 5 of AdS 5 (whereas charges are dimensionless). It turns out that regular BPS black holes with no closed time-like curves only exist when the five charges satisfy certain non-linear constraints. The first constraint relies on the fact that one parameterizes the solutions by four real parameters µ 1,2,3 , Ξ [40]. 2 The second constraint is Alternatively, one can have the same constraint with Ξ substituted by Ξ −1 which corresponds to exchanging J 1 ↔ J 2 . The third constraint is where the entropy S BH is defined in (2.10) below.
Charges and angular momenta of the black holes are completely determined by these four parameters µ I , Ξ with I = 1, 2, 3. Defining In [40] the authors use five real parameters µ 1,2,3 , a, b with 0 ≤ a, b < g −1 , however the black hole charges only depend on the combination Ξ = (1 − b 2 g 2 )/(1 − a 2 g 2 ). The parameters a, b are useful to write the full supergravity solutions. They are determined, in terms of µ 1,2,3 and Ξ, by the extra relation the electric charges and angular momenta are where G N is the five-dimensional Newton constant and It is easy to see that one of the charges Q I can be zero or negative. 3 There are some combinations, though, that we can bound above zero, for instance: (2.8) In particular, at most one charge can be zero or negative. Setting g = 1 for the sake of clarity, we also have for I = K = L = I. The two inequalities follow from (2.3).
The Bekenstein-Hawking entropy is proportional to the horizon area, and can be written as a function of the black hole charges [59]: (2.10) The constraint (2.4) requires the quantity inside the radical to be positive. The BPS solutions have a regular well-defined event horizon only if the angular momenta are non-zero: in other words there is no static limit in gauged supergravity.
The charges are The entropy is Once again, the constraint (2.4) requires the quantity inside the radical to be positive. 4 3 The dual field theory and its index A non-perturbative definition of type IIB string theory on AdS 5 × S 5 is in terms of its boundary dual: 4d N = 4 SYM theory with SU (N ) gauge group [7], where (3.1) The weak curvature limit in gravity corresponds to the large N and large 't Hooft coupling limit in field theory. Up to the choice of gauge group, SYM is the unique four-dimensional Lagrangian CFT with maximal supersymmetry. The field content, in N = 1 notation, consists of a vector multiplet and three chiral multiplets X, Y, Z, all in the adjoint representation of the gauge group. Besides, there is a cubic superpotential W = Tr X[Y, Z]. The R-symmetry is SO(6) R : going to the Cartan U (1) 3 , we choose a basis of generators R 1,2,3 each giving R-charge 2 to a single chiral multiplet and zero to the other two, in a symmetric way.
Considering the theory in radial quantization on R × S 3 , we are interested in the states that can be dual to the BPS black holes described in Section 2. These are 1/16 BPS states preserving one complex supercharge Q, and characterized by two angular momenta J 1,2 on S 3 and three R-charges for U (1) 3 ⊂ SO(6) R . The angular momenta J 1,2 are semi-integer and each rotates an R 2 ⊂ R 4 . Indicating with J ± the spins under SU (2) + × SU (2) − ∼ = SO(4), we set J 1,2 = J + ± J − . With respect to the N = 1 superconformal subalgebra (SCA) that contains Q, we describe the R-charges in terms of two flavor generators q 1,2 = 1 2 (R 1,2 − R 3 ) commuting with Q, and the R-charge r = 1 3 (R 1 + R 2 + R 3 ). All fields in the theory have integer charges under q 1,2 . The counting of BPS states is performed by the superconformal index [41,35] defined by the trace Here p, q, v a with a = 1, 2 are complex fugacities associated with the various quantum numbers, while the corresponding chemical potentials τ, σ, ξ a are defined by The fermion number is defined as F = 2(J + + J − ) = 2J 1 . The index is well-defined for By standard arguments [60], I only counts states annihilated by Q and Q † and is thus independent of β.
It will be convenient to redefine the flavor chemical potentials as and use y a = e 2πi∆a . (3.6) The index becomes 5 Notice that J 1 , J 2 , 1 2 F , 1 2 R 3 are all semi-integer and correlated according to It is then manifest from (3.7) that the index is a single-valued function of the fugacities.
The index (3.2) admits an exact integral representation [41,35,61]. In order to evaluate its large N limit, though, we find more convenient to recast it in a different form, called Bethe Ansatz formula [50,51] (see also [62] for a 3d analog, and [63][64][65][66][67] for similar Higgs branch localization formulas). Computing the large N limit with this formula is still challenging, and in this paper we will restrict ourselves to the case of equal fugacities for the angular momenta: Hence, let us describe the Bethe Ansatz formula with this restriction [50], 6 in the case of N = 4 SU (N ) SYM. The superconformal index reads: This is a finite sum over the solution set {û} to a system of transcendental equations, dubbed Bethe Ansatz Equations (BAEs), given by (3.11) for i = 1, . . . , N and where u ij = u i − u j . We call Q i the BA operators. The unknowns are the "complexified SU (N ) holonomies" u i subject to the identifications meaning that each one lives on a torus of modular parameter τ , and constrained by as well as a "Lagrange multiplier" λ. The function θ 0 is defined as (3.14) in terms of the q-Pochhammer symbol. Some of its properties are collected in Appendix A. The prefactor in (3.10) is defined in terms of the elliptic gamma function [68] Γ(u; τ, σ) = Γ z = e 2πiu ; p = e 2πiτ , q = e 2πiσ = ∞ m,n=0 The function Z is Finally, the Jacobian H is when evaluated on the solutions to the BAEs. Notice that both Q i , κ N , Z and H are invariant under integer shifts of τ , ∆ 1 and ∆ 2 , implying that the superconformal index (3.10) is a single-valued function of the fugacities.
Let us add some comments on how (3.11) and (3.18) are obtained from the general formalism in [51]. The maximal torus of SU (N ) is given by the matrices diag(z 1 , . . . , z N −1 , z N ) with N j=1 z j = 1 and, setting z j = e 2πiu j , is parameterized by u 1 , . . . , u N −1 . For general gauge group G, the BA operators Q i have an index i that runs over the Cartan subalgebra of G. Let us denote the BA operators of SU (N ) as Q 1 , . . . , Q N −1 , then the BAEs are Q j = 1. The BA operators of SU (N ) can be written as Q j = Q j /Q N in terms of the BA operators Q 1 , . . . , Q N of U (N ). Introducing a "Lagrange multiplier" λ, we can set Q N = e −2πiλ and write the BAEs as e 2πiλ Q j = 1 for j = 1, . . . , N (this includes the definition of λ). Absorbing e 2πiλ into Q i , we end up with (3.11).
The Jacobian H for SU (N ) is given by When evaluated on the solutions to the BAEs, we have To see the last equality, one should notice that ∂Q i /∂λ BAEs = 2πi.
The chemical potentials u j are defined modulo 1, and the SU (N ) condition implies that they should satisfy j u j ∈ Z. However, it is easy to check that the BAEs (3.11) are invariant under shifts of one of the u j 's by the periods of a complex torus of modular parameter τ , namely u k → u k + n + mτ for a fixed k. Hence the BAEs are well-defined on N − 1 copies of the torus. Consistently, both H and Z-when evaluated on the solutions to the BAEs-are invariant under shifts of u j by the periods of the torus (see [51] for the general proof).
As one could suspect at this point, the BAEs (3.11) are also invariant under modular transformations of the torus. To that purpose, it might be convenient to rewrite them in terms of the function θ(u; τ ) = e −πiu+πiτ /6 θ 0 (u; τ ) that has simpler modular properties (see Appendix A). When doing that, the term j u ij in the exponential in (3.11) disappears. One easily shows that Q i are invariant under

Exact solutions to the BAEs
When evaluating the BA formula (3.10), the hardest task is to solve the BAEs (3.11). The very same equations appear in the T 2 × S 2 topologically twisted index [12], and one exact solution was found in [16,69]: Hereū is a suitable constant that solves the SU (N ) constraint (3.13); since all expressions depend solely on u ij , we will not specify that constant. Notice that the solution does not depend on the chemical potentials ∆ a . To prove that it is a solution, we compute To go to the second line we used the periodicity relations (A.3). Taking the product over ∆ = {∆ 1 , ∆ 2 , −∆ 1 − ∆ 2 } we precisely reproduce the inverse of the prefactor of (3.11), for every i. Furthermore, notice that the shiftū →ū + 1 N generates a new inequivalent solution that solves the SU (N ) constraint. Repeating the shift N times, because of the torus periodicities, we go back to the original solution. Therefore, (3.22) actually represents N inequivalent solutions.
The solutions (3.25) organize into orbits of P SL(2, Z) with the following action: T : {m, n, r} → {m, n, r + m} , S : {m, n, r} → gcd(n, r) , m n gcd(n, r) , m(n − r) gcd(n, r) (3.26) where the last entry of {m , n , r } is understood mod n . One can check that S 2 = 1. If {m, n, r} have a common divisor, then one can see that also {m , n , r } have that common divisor, and since T, S are invertible, it follows that gcd(m, n, r) ≡ d is an invariant along P SL(2, Z) orbits.
We can prove that if {m, n, r} have gcd(m, n, r) = 1, then they are in the orbit of {1, mn, 0}, i.e. there exists a P SL(2, Z) transformation that maps them to {1, mn, 0}. Indeed, letr = gcd(m, r). We can perform a number of T transformations to reach {m, n,r}. Necessarily gcd(n,r) = 1, therefore an S transformation gives {1, mn, m(n −r)}. Now a number of T transformations gives {1, mn, 0}. On the other hand, we observe that if gcd(m, n, r) = d > 1, then the orbit under P SL(2, Z) is in one-to-one correspondence with the one of {m/d, n/d, r/d}, which is generated by {1, mn/d 2 , 0}. This shows that the number of orbits is equal to the number of divisors d 2 of N which are also squares. Each orbit is generated by {d, N/d, 0}, and is in one-to-one correspondence with the orbit generated by {1, N/d 2 , 0}, which we can regard as the "canonical form".
At this point we recognize that the set of inequivalent solutions in the first class (3.24) (neglecting shifts ofū) is precisely the P SL(2, Z) orbit with gcd(m, n, r) = 1 in the second class (3.25). Indeed, start with a solution of type (3.24) for some N and some coprime integers a, b. Let m = gcd(a, N ) and n = N/m. We can write the solution as We can identifyk = (a/m)j mod n. Since (a/m) and n are coprime, as j runs from 0 to n − 1,k takes all values in the same range once. Moreover there exists s = (a/m) −1 mod n, such that j = sk mod n. In other words, (a/m) is invertible mod n and its inverse s is coprime with n. We can write j = sk + n (3.28) and as j runs from 0 to N − 1, covers a range of length m. Substituting the expression for j we obtain Notice that gcd(b, m) = 1. Indeed, suppose that b and m have a common factor, then this must also be a factor of a, which is a contradiction. Therefore we have the equality of sets {b mod m} = { mod m}. Finally, we set r = bs mod n and we reproduce the expression in (3.25). The values {m, n, r} obtained this way have gcd(m, n, r) = 1. Indeed, suppose they have a common factor, then this must also be a factor of a but not of (a/m), and thus it must also be a factor of b, which is a contradiction.
On the contrary, start with a solution {m, n, r} of type (3.25) with gcd(m, n, r) = 1. It is easy to see, by repeating the procedure, that it is equivalent to a solution of type (3.24) with a = m and b = r (which imply s = 1).

The large N limit
In this Section we take the large N limit of the BA formula (3.10) for the superconformal index. The first part of the Section is technical, and the uninterested reader could directly jump to Section 4.3 where the final result is presented.
In the related context of the T 2 × S 2 topologically twisted index [70,12], it was shown in [16] that the basic solution (3.22) leads to the dominant contribution in the high temperature limit. Assuming that such a solution gives an important contribution in our setup as well, we will start evaluating its large N limit. We will find that it scales as e O(N 2 ) , therefore in the following we will systematically neglect any factor whose logarithm is subleading with respect to O(N 2 ). We will also find that the solution (3.22) is not necessarily dominant in our setup, rather other solutions can compete, and we will thus have to include the contributions of some of the solutions (3.24).
First of all, consider the prefactor κ N in (3.10) and the multiplicity of the BA solutions, whose contribution does not depend on the particular solution. Each BA solution (3.25) has multiplicity N · N !, where the first factor comes from shifts ofū while the second factor from the Weyl group action. Thus, from (3.15), we find This contribution can be neglected at leading order.

Contribution of the basic solution
Here we consider only the contribution of the basic solution (3.22) to the sum in (3.10).
The Jacobian. We use the expression in (3.18). The derivative of Q i with respect to u j can be computed and it gives: This relation holds because we take u 1 , . . . , u N −1 as the independent variables, and fix u N using (3.13). Substituting we get When we evaluate this expression on u ij = τ (j − i)/N , we notice that-for generic values of ∆ a -the terms in the second line are of order O(1). Indeed, the distribution of points u ij generically does not hit any zeros or poles of G. Retaining only the terms in the first line, the Jacobian reads where the diagonal entries are Let us estimate the behavior of A i with N . By the same argument as above, A i contains the sum of N elements of order O(1) and thus it scales like O(N ) (or smaller). The determinant can be computed at leading order and it gives This scales as O(N N ), therefore log H = O(N log N ) and can be neglected.
The functions Γ. The dominant contribution comes from the function Z defined in (3.17).
To evaluate it, let us analyze N i =j log Γ(u ij +∆; τ, τ ) with ∆ ∈ {∆ 1 , ∆ 2 , ∆ 1 +∆ 2 } separately. Making use of the relation (A.17) proven in [68], we write The function ψ(t) is defined in (A.10), while To make progress, we perform a series expansion of log θ 0 and log ψ, evaluate this expansion on the basic solution for u ij in (3.22), and perform the sum N i =j . We define the modular transformed variables where we introduced A which denotes the following sum over ij: The series can be resummed to N log θ 0 , however we do not need that. We collect the terms into two groups: where the second term comes from = N j. For |q| < |ỹ| < 1, namely for the series converges. The second term is suppressed at large N , whereas the first term is of order O(N ) and can be neglected.
We then perform a similar analysis of log ψ, using the series expansions of the functions log and Li 2 . We find where we used that the following sum vanishes: Once again, the expression can be resummed by breaking the sum into two groups (corresponding to generic and = N j): (4.17) The first term (that comes from setting A → N ) is of order O(N ) and can be neglected. The second term is (4. 18) In the regime of convergence (4.14) this series goes to zero as N → ∞. We conclude that the only contribution at leading order in N is from the polynomial Q in (4.9).
The limit we computed is valid as long as ∆ satisfies (4.14). That inequality has the interpretation that ∆ should lie inside an infinite strip, bounded on the left by the line through −1 and τ − 1, and on the right by the line (that we dub γ) through 0 and τ (see Figure 1). On the other hand, Γ(u ij + ∆; τ, τ ) is a periodic function invariant under shifts ∆ → ∆ + 1. Therefore, unless ∆ sits exactly on one image of the line γ under periodic integer shifts, there always exists a shift that brings ∆ inside the strip. This means that we can use our computation to extract the limit for all ∆ ∈ C \ {γ + Z}.
Let us define the periodic discontinuous function Figure 1: In yellow is highlighted the domain (4.14) in the complex ∆-plane. The right boundary is the line γ, passing through 0 and τ . The left boundary is the line γ − 1, passing through −1 and τ − 1. The dashes lines are other elements of γ + Z.
The function is not defined for Im(∆/τ ) ∈ Z × Im(1/τ ). Essentially, this function is constructed in such a way that [∆] τ = ∆ mod 1, and [∆] τ satisfies (4.14) when it is defined. It also satisfies We use such a function to express the limit as This expression is, by construction, invariant under ∆ → ∆ + 1. The lines that we have dubbed γ + Z, are Stokes lines: they represent transitions between regions in the complex ∆-plane in which different exponential contributions dominate the large N limit, and along which the limit is discontinuous. 7 We do not know what is the limit along the lines, because different contributions compete and a more precise estimate would be necessary to evaluate their sum. We will elaborate on Stokes lines in Section 4.3.
The term with ∆ = 0 requires a special treatment, because it does not satisfy (4.14). We can still use the expansion (4.8). The term log θ 0 is evaluated as 7 Stokes lines divide the complex plane into regions in which the limit gives different analytic functions.
Because of their origin, Stokes lines have the property that only the imaginary part of the function can jump, while the real part must be continuous. One can indeed check that (4.21) satisfies this property.
To calculate the first term in the second line, we notice that x−e 2πij/N , and substituting x = 1 we get N = N −1 j=1 1 − e 2πij/N . At this point we can shift j by k units and multiply over k: To compute the second term we use the series expansion as before. We see that log θ 0 contributes at order O(N log N ) and can be neglected. The product of terms log ψ gives In the first equality we changed sign to u ij because it is summed over ij; to go to the second line we used (A.11). This term is of order O(N ) and can be neglected. We conclude that The expression depends on [∆ 1 ] τ , [∆ 2 ] τ and [∆ 1 + ∆ 2 ] τ . We notice the following relation: 27) The second one can be rewritten as The large N limit of the summand is then where we have introduced the following function for compactness: The two cases were defined in (4.27).
We can rewrite the function Θ in a way that will be useful in Section 6. Define an auxiliary chemical potential ∆ 3 , modulo 1, such that It is also useful to define the primed bracket The primed bracket selects the image of ∆, under integer shifts, that sits inside the strip on the right of the line γ through zero and τ , as opposed to the strip on the left. Hence Irrespective of the integer appearing in (4.31), the bracketed potentials satisfy the following constraints:  [69]. Some of those solutions might contribute at the same leading order in N A class of inequivalent solutions-particularly simple to study-that contribute at leading order in N is obtained through T -transformations: These are the solutions {1, N, r} in the notation of Section 3.1. To evaluate their contribution, simply notice that both Z in (3.17) and H are invariant under τ → τ + r, thus the contribution of (4.35) is the same as in (4.29) but with τ → τ + r. In the large N limit, r runs over Z.
We have not evaluated the contribution of all other {m, n, r} solutions, which is a difficult task. However, in order to have an idea of what their contribution could be, let us estimate the contribution from the S-transformed solution which is {N, 1, 0} in the notation of (3.25). The large N limit of κ N does not depend on the solution, and is subleading. The large N limit of log H is computed in the same way as in Section 4.1, and it gives O(N log N ) or smaller. Let us then analyze Z. In the regime |q| 2 < |y| < 1 we can directly expand log Γ in its plethystic form: If |y| is outside the range of convergence of the plethystic expansion, either above or below, we can simply shift ∆ → ∆ ± τ . This gives a shift by which can be treated in a similar way. This confirms that the estimate above is valid for all ∆'s, even outside the original regime of convergence. The case ∆ = 0 requires a special treatment. We have Thus, there is no contribution from log Z at leading order in N .
In the following we will assume that the only solutions contributing at leading order, namely O(N 2 ), are the T -transformed solutions.

Final result and Stokes lines
Since we end up with competing exponentials, the one with the largest real part dominates the large N limit. Assuming that solutions other than the T -transformed of the basic one are of subleading order in N , we find the final formula lim N →∞ log I(q, y 1 , y 2 ) = max (4.40) The function Θ is defined in (4.30). The meaning of max is that we should choose the value of r ∈ Z such that the real part of the argument is maximized. One of the good features of eqn. (4.40) is that it is periodic under integer shifts of τ, ∆ 1 , ∆ 2 . We already observed that Θ is periodic in ∆ 1,2 because the functions [∆ 1,2 ] τ are. Taking the max over τ → τ + r gives periodicity in τ as well. This implies that the RHS of (4.40) is actually a single-valued function of the fugacities q, y 1 , y 2 . This is a property of the index at finite N , as manifest in (3.7) and (3.10), and it is reassuring that the large N expression we found respects the same property.
The function I ∞ has a complicated structure. The full range of allowed fugacities q, y 1 , y 2 gets divided into multiple domains of analyticity, separated by Stokes lines. In each domain of analyticity, only one exponential contribution (for some value of r) dominates the large N limit: the function log I ∞ takes the form of a simple rational function given by Θ(∆ 1 , ∆ 2 ; τ + r). The Stokes lines are real-codimension-one surfaces, in the space of fugacities, that separate the different domains. When crossing a Stokes line, a different exponential contribution dominates, and log I ∞ takes the form of a different rational function. In particular, on top of a Stokes line there are two (or more) exponential contributions that compete: their exponents have equal real part. This characterizes the locations of Stokes lines. In terms of the function Θ: for some r 1,2 ∈ Z.
In fact, also the values of ∆ 1 , ∆ 2 such that Θ(∆ 1 , ∆ 2 ; τ + r) is discontinuous (for the value of r picked up by max) should be regarded as forming a Stokes line. In this case, the two competing exponents correspond to the values of Θ on the two sides of the discontinuity. There are two possible sources of discontinuity. First, one of the bracket functions, say [∆ 1 ] τ , could be discontinuous. This happens when Im(∆ 1 /τ ) ∈ Z × Im(1/τ ), namely when α ≡ lim →0 + [∆ 1 − ] τ /τ ∈ R. Taking into account that on the left of the discontinuity we are in the 1 st case, while on the right we are in the 2 nd case-in the terminology of (4.27)-and assuming that ∆ 2 is generic, we find where the limit is taken with real positive. Second, we could pass from the 1 st to the 2 nd case of the definition (4.30). This happens when [∆ 1 ] τ + [∆ 2 ] τ + 1 = α τ for some α ∈ R. Assuming that ∆ 1,2 are otherwise generic, we find (4.43) In both cases we confirm that the codimension-one surface of discontinuity is a Stokes line, because Im Θ is equal on the two sides.
When we sit exactly on a Stokes line, two (or more) exponential contributions compete, and in order to compute the large N limit we should sum them. However we do not know the relative phases, because they are affected by all subleading terms and a more accurate analysis would be required. Therefore, we cannot determine the large N limit of the index along Stokes lines.
It turns out that a value of r that maximizes the real part of the argument of max may or may not exist. We can estimate the behaviour of the real part at large r by noticing that This implies that Thus, the real part of the argument of max approaches a constant value. If there is no maximum but rather the constant value is a supremum, then our computation is not finished: All contributions from the T -transformed solutions should be summed, however for large |r| they form an infinite number of competing exponentials, whose sum crucially depends on how they interfere. In order to determine such a sum we would need more accurate information.
We conclude by stressing that-even though only the dominant exponential determines the large N limit of the index-we expect that all exponential contributions, including the subdominant ones, have some physical meaning. Each of them plays the role of a "saddle point", although our treatment is not the standard saddle-point approximation. We will make this comment more concrete in Section 6, when comparing the large N limit of the index with BPS black hole solutions in supergravity.

Comparison with previous literature
The large N limit of the superconformal index of N = 4 SYM was already computed in [35]. There, it was found that the large N limit does not depend on N , and therefore it does not show a rapid enough growth of the number of states to reproduce the black hole entropy. In this Section we would like to explain how the results here and there can be compatible.
The authors of [35] took the large N limit of the index, for real fugacities. Their result, in our notation and restricted to the case p = q, is lim N →∞ I(q, y 1 , y 2 ) = ∞ n=1 1 1 − f q n , y n 1 , y n 2 (4.46) (4.47) In particular, log I is of order O(1). On the contrary, we computed the large N limit for generic complex fugacities, and found that log I is of order O(N 2 ). It was already discussed in [48], in a double-scaling Cardy-like limit, that the large N limit of the index is completely different for real and complex fugacities, and it was observed in [49] that there exists a deconfinement transition once complex fugacities are taken into account.
The resolution we propose relies on the fact that, for complex fugacities, the limit shows Stokes lines. As we described, along those codimension-one surfaces multiple exponentials compete. In order to know what the limit is there, we would need to sum those competing exponentials, but this requires a more accurate knowledge of the subleading terms.
What we notice, though, is that the codimension-three subspace of real fugacities is precisely within a Stokes line. Therefore, although we cannot prove it, it is conceivable that the competing terms cancel exactly, leaving the O(1) result (4.46). Indeed, in Appendix B we prove the following result, which is stronger than the statement that we sit on a Stokes line. Take the angular fugacity q to be real positive, namely 0 < q < 1 and set τ ∈ iR ≥0 for concreteness, and take the flavor fugacities y 1,2 to be real. Then Θ(∆ 1 , ∆ 2 ; τ ) is along a Stokes line and is not defined, while for r > 0. On the other hand, take the angular fugacity real negative, namely −1 < q < 0 and set τ ∈ − 1 2 + iR ≥0 , and take again the flavor fugacities to be real. Then for r ≥ 0. Therefore, among the various contributions from T -transformed solutions parameterized by r ∈ Z, there is an exact pairing of all well-defined terms where, in each pair, two terms have the same real part and can conceivably cancel. In other words, not only the term with maximal real part can cancel, but also all other terms we computed at order O(N 2 ). This scenario is a strong check of our result, that makes it compatible with [35].

Statistical interpretation and I-extremization
We wish to extract the number of BPS states, for given electric charges and angular momenta, from the large N limit of the exact expression (3.10) of the superconformal index. Since the latter counts states weighted by the fermion number (−1) F , one may worry that strong cancelations take place and that the total number of states is not accessible. However, one can argue [10,11] that the index (3.2) or (3.7) is equal to I(p, q, y 1 , y 2 ) = Tr e iπR trial (τ,σ, where the trace is taken in the IR N = 2 super quantum mechanics (QM) obtained by reducing the 4d theory on S 3 , R trial is a trial R-symmetry, and C 1,2,3,4 are the charges appearing in (3.7): Indeed, because of the relations (3.8), we can represent the fermion number as (−1) F = e iπR 3 . Substituting in (3.7) and separating the chemical potentials into real and imaginary part, we obtain the expression (5.1) with From the point of view of the super QM, R 3 is an R-symmetry while the other four operators are flavor charges, hence R trial is an R-symmetry. We see in (5.1) that only the first exponential can produce possibly-dangerous phases, while the other two are real positive. Now, for a single-center black hole in the microcanonical ensemble, the near-horizon AdS 2 region is dual to an N = 2 superconformal QM. The black hole states are vacua of the su(1, 1|1) 1d SCA. Since we are in the microcanonical ensemble, each of those states is invariant under the global conformal algebra su(1, 1) ∼ = so(2, 1) (because AdS 2 is) as well as under the fermionic generators (because the black hole is supersymmetric). This necessarily implies that those states are invariant under the superconformal R-symmetry u(1) sc ⊂ su(1, 1|1), i.e. that they have vanishing IR superconformal R-charge R sc . Thus, when R trial is tuned to R sc , the index counts the black hole states with no extra signs or phases (this is similar to [71,72]). Of course, in a given charge sector there will be more BPS states than just the single-center black hole, but assuming that the single-center black hole dominates, the index captures its entropy. It remains to understand how to identify R sc . At large N , the entropy is extracted from the index with a Legendre transform, and this operation can be argued to effectively select R sc among the R trial 's (this large N principle was dubbed I-extremization in [10,11]).

(5.4)
Here the trace is over states with {Q, Q † } = 0, and we have identified Q I = R I /2 (for I = 1, 2, 3) with the electric charges in supergravity. We recognise that the black hole angular momenta J 1,2 are associated with the chemical potentials τ, σ and the charges Q 1,2,3 with ∆ 1,2,3 . The microcanonical degeneracies at fixed quantum numbers are extracted by computing the Fourier transform of (5.4). However, since ∆ 3 is not an independent variable, what we obtain are the degeneracies for fixed values of the four charge operators appearing in (3.7), summed over Q 3 . Using the supergravity notation, those four fixed charge operators are Thus, what we can compute is where d(J, Q) are the (weighted) degeneracies with all charges J 1,2 and Q 1,2,3 fixed.
Nevertheless, we can take advantage of the fact, reviewed in Section 2, that the charges of BPS back holes are constrained, and for fixed C 1,2,3,4 there is at most one black holefor a certain value of the fifth charge Q 3 . We can then use (5.6) to extract its degeneracy d(J, Q) = exp S BH (J, Q) at leading order because the latter will dominate the sum over Q 3 .
In the large N limit, the integral (5.6) reduces by saddle point approximation to a Legendre transform with respect to the independent variables {τ, σ, ∆ 1 , ∆ 2 }: where hatted variables denote the critical point. In this approach, Q 3 can be determined as the unique value that makes the entropy S BH (J, Q) real [10,11].
In the particular case of 4d N = 4 SYM, the large N limit of the index is a function with multiple domains of analyticity, separated by Stokes lines. This makes things more interesting. In each domain we should perform the Legendre transform, and whenever the critical point falls inside the domain itself, we obtain a self-consistent contribution to the total entropy. Even more generally, we have written the index as a sum of competing exponentials (one for each Bethe Ansatz solution) and we can compute the Legendre transform of each of those exponentials-irrespective of which one dominates. We expect each contribution to represent the entropy of some classical solution-very similarly to a standard saddle pointeven when the entropy is smaller than that of the dominant solution.
6 Black hole entropy from the index In this Section we show that the contribution of the basic solution (3.22) to the superconformal index at large N , in the domain of analyticity that we called "1 st case" in (4.27), given by precisely reproduces the Bekenstein-Hawking entropy (2.13) of single-center black holes in AdS 5 (this is in line with the result of [48] in a double-scaling Cardy-like limit). It amounts to show that the Legendre transform of (6.1) is the black hole entropy (this will be reviewed below), and that the critical point involved in the Legendre transform consistently lies within the domain of analyticity in which (6.1) holds.
Recall that the contribution of the basic solution corresponds to the r = 0 sector in (4.40). For black holes with large charges, i.e. for black holes that are large compared with the AdS 5 scale, that is indeed the dominant contribution to the index. However, intriguingly enough, as we reduce the charges the contribution of the single-center black hole may cease to dominate. We will highlight this phenomenon in Section 6.1 in the very special case of black holes with equal charges. This seems to suggest that, below a certain threshold, the BPS black holes may develop instabilities, possibly towards hairy or multi-center black holes. Indications that this is the case have also been given in [48,49]. It would be nice if there was a connection between this observation and recently constructed hairy black holes in AdS 5 [53][54][55][56].
The entropy function. The Legendre transform of the black hole entropy (2.10) in the general case, also called entropy function, was obtained in [52]. Let us review it, following the detailed discussion in Appendix B of [47]. The entropy function is and with the constraint Because of the constraint, S is really a function of four variables. The entropy S BH is the Legendre transform of S with its constraint. We can compute it as the critical point of in which the constraint is imposed with a Lagrange multiplier Λ. The equations for the critical point are 5) and the constraint (6.3). In details, It follows that with It turns out that we can find the value of S at the critical point without knowing the exact solution for the critical point. We use the fact that S is homogeneous of degree 1 (it is a monomial), and thus Substituting into (6.4) we find Since Λ is the solution to the cubic equation (6.7), it looks like there are three possible values for the entropy. However, since for real charges the cubic equation has real coefficients, we either find 3 real roots or 1 real and 2 complex conjugate roots for Λ. Imposing that the entropy be real positive, we require that there is 1 real and 2 imaginary conjugate roots, then only one of them-the one along the positive imaginary axis-leads to an acceptable value for the entropy. Since (Λ − β)(Λ − iα)(Λ + iα) = Λ 3 − βΛ 2 + α 2 Λ − βα 2 , we obtain the following constraint on the charges: One can check that the parameterization (2.6) automatically solves the first equation. Then the roots of (6.7) are Λ ∈ {−p 2 , ±i √ p 1 }. The physical solution is which is precisely eqn. (2.10). We stress that the conditions (6.11) are necessary, but not sufficient, to guarantee that the supergravity solution is well-defined. 8 It is not difficult to write the values of the chemical potentials at the critical point. To simplify the notation, let us define and use an index A = 1, . . . , 5. The equations (6.6) imply that Φ A P A are all equal for A = 1, . . . , 5 . (6.14) Implementing the constraint (6.3), the solution is Since, even for real charges, the P A 's are complex, the solutions Φ A are in general complex.
Equal angular momenta. Let us specialize the formulas to the case J 1 = J 2 ≡ J, and determine useful inequalities satisfied by the chemical potentials at the critical point. First of all, from the constraint (6.3) it immediately follows At the critical point (6.15) one finds As an example, take Q 1 = Q 2 = Q 3 ≡ Q and J 1 = J 2 ≡ J. The first equation in (6.11) is solved by J = −3Q − 1 ± (2Q + 1) 3/2 , and both branches are covered by the parameterization (2.6) as Q = µ + µ 2 /2, J = 3µ 2 /2 + µ 3 . Then p 1 = 3Q 2 + 6Q + 2 ∓ 2(2Q + 1) 3/2 and one can check that, for Q > 0, both branches have p 1 > 0. However, only the branch with upper sign satisfies also (2.3)-here µ > 0-and corresponds to well-defined supergravity solutions, while the branch with lower sign does not.
To obtain the last inequality we used that Q a + J > 0 for the BPS black holes, as we showed in (2.9). This implies that Using the explicit parameterization (2.12) presented in Section 2 (and setting g = 1 for the sake of clarity), one can also show that , where p 1 = ν 2 (1 + γ 1 )γ 3 − 1 4 γ 2 2 . In particular, the first equation shows that Entropy from the index. Finally, we compare the contribution to the index from the basic solution in the 1 st case, given in (6.1), with the entropy function S in (6.2). The latter, after eliminating X 3 with the constraint (6.3) and restricting to equal angular fugacities, reads We see that it is exactly equal to (6.1), as long as we can identify τ = ω , [∆ a ] τ = X a for a = 1, 2, 3 .

(6.22)
This is not obvious, but we can check that it is indeed possible. First of all, X 1 and X 2 should satisfy the strip inequalities that [ · ] τ does, at least in a neighbourhood of the critical point. This is precisely what we proved in (6.18). Second, the fugacities at the critical point should also satisfy the inequalities (4.27) that define the 1 st case. Because of the constraint, this is the same as requiring that also X 3 satisfies (6.18), which is true. Thus, this concludes our proof. Let us stress that, in our approach, the constraint (6.3) with the correct constant term simply comes out of the large N limit.
One could wonder what is the physics described by the domain of analyticity named 2 nd case in (4.30). It appears that it reproduces the very same black hole entropy as the 1 st case. Indeed, as apparent from (4.33), in the two cases Θ takes almost the same form, the only difference being that [ · ] τ and [ · ] τ satisfy opposite strip inequalities and a constraint with opposite constant term. It was already observed in [47] that the entropy function S reproduces the black hole entropy with either one of the two constraints imposed. We leave for future work to understand what is the role of such a twin contribution.

Example: equal charges and angular momenta
In order to make some of the previous statements more concrete, we now study in detail a very special case in which the index counts states with equal charges Q 1 = Q 2 = Q 3 ≡ Q and angular momenta J 1 = J 2 ≡ J. This will be instructive to elucidate the structure of Stokes lines.
Let us first quickly summarize the properties of black holes and their entropy in this case [36]. We set ν = 1 (all charges are in "units" of ν) so that 23) and the charge constraint is This is quadratic in J and potentially leads to two branches of solutions. However only one of them satisfies (2.3) when parameterized in terms of µ (we also set g = 1): The entropy is positive for Q > 0, and in this range J > 0.
The extremization problem (6.4) simplifies because we only have two chemical potentials, X ≡ X 1 = X 2 = X 3 and ω with the constraint (6.3). The critical point is Let us mention that in the alternative extremization problem in which the constraint (6.3) is modified by changing +1 into −1, the critical values of X and ω are given by the same expressions, however the critical value of the Lagrange multiplier becomes Λ = −i √ p 1 .
We now turn to the index. Given the identifications X a = [∆ a ] τ and ω = τ , we can restrict to chemical potentials such that [ where ∆ 3 is defined through the general constraint (4.31). Up to integer shifts, this amounts to The critical points (6.26) indeed satisfy this relation. We have thus reduced to a single independent chemical potential τ . Notice that the function I ∆(τ ); τ = I 2τ −1 3 ; τ is periodic under τ → τ + 3, therefore we will restrict to 0 ≤ Re τ < 3.
We study the large N formula (4.40) for the index, in particular we want to determine the structure of the leading contributions as τ is varied, and where the Stokes lines are. To do so, we need the values of the bracketed potentials [∆] τ +r for r ∈ Z. We find (6.28) In the second case the bracket is not defined because Im ∆/(τ + r) ∈ Z × Im 1/(τ + r) , i.e. because ∆ sits exactly on the boundary of a strip defined by τ + r. We can however consider [∆] τ +r for values of ∆ that are a bit off the boundary of the strip in the real direction. We consider the values ∆ (±) = ∆ ± with infinitesimal > 0 and find Using these formulas, the values of Θ(∆, τ + r) are easily computed. 9 In particular, the imaginary parts of Θ computed on ∆ (±) are the same.
The dominant contribution to the index is determined by comparing the absolute values of exp −πiN 2 Θ(∆; τ + r) -or equivalently the imaginary parts of Θ-as we vary r. When there is a particular value r for which Im Θ(∆; τ + r) is maximum, there is one dominant contribution which leads to a concrete estimate of the leading behavior of the index. When, instead, there is no maximum, we are left with an infinite number of competing contributions and more detailed information would be needed to resum them. We obtain the following values for the imaginary part of Θ: if r = 2 mod 3 .

(6.30)
Notice that the limiting value for large |r| (equal to the value for r = 1 mod 3) is as in (4.45). If there is a value of r that maximizes Im Θ, it must come from the first or third case. In particular, there exists r with r = 0 mod 3 if τ satisfies the following relation: Re τ + r > 3 |τ + r | 2 with r = 0 mod 3 . This corresponds to the interior of another semi-circle of radius 1/6, centered at τ = −1/6− r.
The two inequalities (6.31) and (6.32) define two semi-circles in the fundamental range 0 ≤ Re τ < 3, for r = 0 and r = −1 respectively, as well as all their images under the periodicity τ → τ + 3. On the other hand, outside the two regions there is no dominant contribution because, for all values of r, Im Θ is smaller than the limiting value. In Figure 2 we provide plots of Im Θ(∆; τ + r) as r is varied, both for τ insides the semi-circle (6.31), inside the semi-circle (6.32), and outside those two.
In Figure 3 we represent the fundamental range 0 ≤ Re τ < 3 of the upper half τ -plane, dividing it into regions according to the dominant contribution. In Figure 4 we represent the same information in the q-plane, using q 1/3 as the variable. The red semi-circle (6.31) corresponds to the values of τ in which r = 0, while the green semi-circle (6.32) corresponds to r = −1. These are two different domains of analyticity. The remaining "no max" region, in blue, corresponds to values of τ for which there is no dominant contribution. The three regions are separated by Stokes lines (in black).
Inside the red semi-circle (6.31) the large N limit of the superconformal index is This expression exactly matches the entropy function (6.2) of black holes with equal charges Q and angular momenta J, and its Legendre transform selects the critical points (6.26). We represent the line of critical points, as µ > 0 is varied, by a blue solid line in Figures 3 and  4. As we see from there, for µ > µ * the blue line lies inside the red semi-circle, meaning that the entropy of the single-center black hole is the dominant contribution to the index. This seems to confirm that "large" BPS black holes, with Q > Q * or equivalently J > J * , are stable. On the contrary, for 0 < µ < µ * the blue line plunges into the "no max" region. We can still identify the black hole entropy with the contribution of the basic solution (3.22) to the index, however such a contribution is no longer dominant. This suggests that "small" BPS black holes with Q < Q * might be unstable towards other supergravity configurations. : Stokes lines in the q-plane, where the variable is q 1/3 . The notation is the same as in Figure 3.
We find the following values at the transition point: where Q, J, S are in units of ν. It would be nice to derive these values from supergravity.
The green circle in Figures 3 and 4 corresponds to values of τ for which the r = −1 contribution dominates. In this domain we find It is interesting to draw the subspace where both fugacities q and y are real and the computation of [35] applies. We include this subspace both in Figure 3 and, in terms of q 1/3 , in Figure 4. We see that the real subspace does not intercept the black hole lines: it only asymptotically reaches them, at the tail that describes black holes much smaller than the AdS radius.
To derive it, one can relate θ 0 to the Dedekind η and Jacobi θ 3 functions: The function θ 3 is also called ϑ 00 in the literature.
Function θ. This is a modification of the function θ 0 , defined as θ(u; τ ) = e −iπu+iπτ /6 θ 0 (u; τ ) . (A.7) Its periodicity relations are The modular transformation is Function ψ. Following [68], we define the function Within Im t < 0, the definition is analytic and single-valued. The branch of the logarithm is determined by the series expansion log(1 − z) = − ∞ k=0 z k /k, whereas Li 2 (z) = ∞ k=1 z k /k 2 . One can show that the branch cut ambiguities of the logarithm and the dilogarithm, that appear for Im t ≥ 0, cancel in the definition of ψ(t). This means that the latter function can be analytically continued to the whole complex plane yielding a meromorphic function.
Two useful properties of ψ(t) are: valid for any t ∈ C.
Function Γ. Setting now p = e 2πiτ , q = e 2πiσ , the elliptic gamma function [68] is defined for |p|, |q| < 1. The function Γ(z; p, q) is meromorphic in z, with simple zeros at z = p m+1 q n+1 and simple poles at z = p −m q −n , for m, n ∈ Z ≥0 . A plethystic representation is which is convergent for |pq| < |z| < 1.
In the degenerate case τ = σ the identities (A.16) do not apply. However, there exists an alternative version [68]: The function ψ(t) is defined in (A.10).

B Real fugacities
In Section 4 we evaluated, in the large N limit, the contribution of some of the solutions to the BAEs to the sum in (3.10). In particular we found that all T -transformed solutions (4.35) of the basic solution, parameterized by the integer r, contribute at the same order in N , and their contributions are the arguments of max in the final formula (4.40): in terms of Θ defined in (4.30) and with r ∈ Z.
Here we show that when we take the fugacities q, y 1 , y 2 to be all real, we end up precisely on a Stokes line. More precisely, we show that all contributions for r ∈ Z organize into pairs, except for those elements that already sit on a Stokes line determined by the discontinuity of one of the functions [ · ] τ . In each pair, the two contributions have equal real part and compete. We cannot compute the sum of the two terms, as this would require more accurate information about the subleading corrections. Yet, this makes our result compatible with the result of [35]. There it was found that, for real fugacities, the index scales as O(1) at large N , implying that all O(N 2 ) contributions cancel out. This point was also stressed in [48,49].
Real fugacities corresponds to chemical potentials whose real part is either zero or −1/2 modulo 1. We distinguish the various possibilities into two major cases: the case that 0 < q < 1, corresponding to τ ∈ iR ≥0 , and the case that −1 < q < 0, corresponding to τ ∈ − 1 2 + iR ≥0 . Each case is further divided into subcases, according to the number of positive flavor fugacities y 1,2 .
B.1 The case 0 < q < 1 We start with the case of positive angular fugacity, 0 < q < 1. We take τ ∈ iR ≥0 and write We distinguish three different subcases, corresponding to y 1 , y 2 being both positive, one positive and one negative, or both negative.
If one of the flavor fugacities-that we call y-is real positive, we set the corresponding chemical potential We immediately see that [∆] τ is not defined, because the argument sits precisely along one of the lines of discontinuity. On the other hand, for r > 0 and generic δ, the functions [∆] τ ±r are well-defined and we would like to evaluate them. We can precisely determine their values by splitting the imaginary axis in the ∆-plane into a series of intervals Assuming that δ ∈ I k , we see that ∆ can be brought inside the strip corresponding to τ + r by shifting it by k, and inside the strip corresponding to τ − r by shifting it by −k − 1. In formulas: If δ is equal to an extremum of I k , i.e. if δ = k t/r for some k ∈ Z, then [∆] τ ±r are not defined.
On the other hand, if one of the flavor fugacities-that we keep calling y-is real negative, we set its chemical potential For r > 0 and generic δ, both functions [∆] τ ±r are well-defined. As before, we can determine their values by splitting the imaginary axis in intervals. This time the intervals are We then find If δ = (2k − 1)t/2r for some k ∈ Z, then [∆] τ ±r are not defined.
We now proceed to applying these formulas to the three subcases.
Let us now consider r > 0. For generic δ a , the functions [∆ a ] τ ±r are well-defined. Precisely, for δ a ∈ I ka the functions [∆ a ] τ ±r are given by (B.4). Turning to ∆ 1 + ∆ 2 , we have two possibilities: In the first case, given by (B.9), we compute In the second case, given by (B.10), we compute From these expression we see that, in both cases,
B.1.2 The subcase y 1 < 0 < y 2 with 0 < q < 1 We take one flavor fugacity to be positive and one negative, say y 1 < 0 < y 2 (recall that the index is symmetric in the two flavor fugacities). Correspondingly, we set Similarly to the previous case, [∆ 2 ] τ is not defined and the contribution r = 0 is already along a Stokes line.
B.2 The case −1 < q < 0 Now we move to the case of negative angular fugacity, −1 < q < 0, and set Once again, we distinguish three different subcases corresponding to y 1 , y 2 being both positive, one positive and one negative, or both negative. First, let us discuss the new intervals we need.