Asymptotic growth of the 4d $\mathcal N=4$ index and partially deconfined phases

We study the Cardy-like asymptotics of the 4d $\mathcal N=4$ index and demonstrate the existence of partially deconfined phases where the asymptotic growth of the index is not as rapid as in the fully deconfined case. We then take the large-$N$ limit after the Cardy-like limit and make a conjecture for the leading asymptotics of the index. While the Cardy-like behavior is derived using the integral representation of the index, we demonstrate how the same results can be obtained using the Bethe ansatz type approach as well. In doing so, we discover new non-standard solutions to the elliptic Bethe ansatz equations including continuous families of solutions for $SU(N)$ theory with $N\ge3$. We argue that the existence of both standard and continuous non-standard solutions has a natural interpretation in terms of vacua of $\mathcal N=1^*$ theory on $\mathbb R^3\times S^1$.


Introduction
For the first time it has become possible to analyze a black hole-counting index [1,2] both in a Cardy-like limit [3][4][5] and in a large-N limit [6]. In the present work we further study the asymptotics of the 4d N = 4 index of [2], finding that the two limits shed light not only on each other, but also on new black objects in the dual AdS 5 theory. Our work is motivated by the important recent discovery [6,7] that by varying the fugacity parameters of the index, its large-N asymptotics exhibits a Hawking-Page-type deconfinement transition [8,9] from a multi-particle phase-already observed in [2]-to a long-anticipated black hole phase [10][11][12][13][14][15].
Here we argue that by varying the chemical potentials in the index, its Cardy-like asymptotics displays "infinite-temperature" Roberge-Weiss-type first-order phase transitions [16] between the fully-deconfined phase associated to black holes [10][11][12][13][14][15], and confined or partiallydeconfined phases, with the latter possibly associated to new multi-center black objects. Guided by this Cardy-limit analysis, we revisit the large-N asymptotics of the index and argue that by including some previously neglected contributions in the asymptotic analysis of [6] it is possible to see the partially-deconfined phases in the large-N limit as well.
In the rest of this introduction, we present a more precise description of our framework, as well as a brief outline of our new technical results. Section 2 spells out our terminology regarding various "phases" of the index. There we outline a correspondence with N = 1 * theory, which turns out to yield surprisingly powerful insight into the Bethe Ansatz approach discussed later in the paper. In Section 3 we study the Cardy-like limit of the index using its expression as an integral over holonomy variables [2,17], extending previous partial results by one of us in [5]. For the SU(2) case we explain that varying the chemical potentials triggers an "infinite-temperature" Roberge-Weiss-type transition [16] between a confined phase where the center-symmetric (i.e. Z 2 -symmetric) holonomy configuration dominates the index, and a deconfined phase where two center-breaking holonomy configurations take over. For N = 3 we establish a similar behavior with a Z 3 center-breaking pattern, while for N = 4 we encounter a partially deconfined phase with a Z 4 → Z 2 center-breaking pattern. We also consider the SU(N > 4) cases, and in particular argue that taking the large-N limit after the Cardy-like limit should yield various partially-deconfined infinite-temperature phases. Our investigation of this double-scaling limit leads up to a conjecture for the leading asymptotics of the index as displayed in (3.19).
In Section 4 we study the index using its expression as a sum over solutions to a system of elliptic Bethe Ansatz Equations (eBAEs) [18,19]. First we review the Bethe Ansatz formula for the 4d N = 4 index. Then we study the Cardy-like asymptotics of the index in this approach. It turns out that compatibility with the same partially-deconfined infinitetemperature phases observed in the integral approach requires the existence of new solutions to the eBAEs, which are not covered in [20]. We find a vast number of such new solutions in this work; in most cases numerically, in some low-rank cases asymptotically, and in one case exactly! The exact solution comes from the remarkable correspondence with N = 1 * theory, discussed in section 2. The correspondence also gives powerful insight into continua of eBAE solutions which exist for N ≥ 3. The existence of such continua of solutions implies in fact that the Bethe Ansatz formula in its current form [18,19] as a finite sum is incomplete for N > 2, and calls for an integration with a so-far unknown measure, which we leave unresolved. Finally we move on to the large-N limit of the index, extending previous results by Benini and Milan in [6]. Section 5 summarizes our main findings by placing them in the context of recent literature, and also outlines a few important related directions for future research. The appendices elaborate on some technical details used in the main text.
Since ∆ 1,2 are defined mod Z, we can focus on a fundamental domain. It turns out to be useful [5] to take the fundamental domain to consist of the two wings 0 < ∆ 1 , ∆ 2 , 1−∆ 1 −∆ 2 < 1 (upper-right) and −1 < ∆ 1 , ∆ 2 , −1 − ∆ 1 − ∆ 2 < 0 (lower-left) of a butterfly in the ∆ 1 -∆ 2 plane. When arg β > 0, the effective pairwise potential for the holonomies in the Cardy-like limit, given explicitly in (3.6) below, is M-shaped on the upper-right wing of the butterfly, while on the other wing it is W-shaped-and conversely for arg β < 0. See Figure 1 for a representation of the fundamental domain along with the M-and W-wings. Throughout this paper, the M-wing (resp. W-wing) denotes the part of fundamental domain where the effective pairwise potential for the holonomies is M-shaped (resp. W-shaped), namely the upper-right (resp. lower-left) wing of the butterfly when arg β > 0 and the lower-left (resp. upper-right) wing when arg β < 0.
In [5] the asymptotics of the index on the M wings was obtained. On the M wings, because of the shape of the pairwise potential, the holonomies condense in the Cardy limit, and the "saddle-points" with x ij (:= x i − x j ) = 0 dominate the matrix-integral expression (1.2) for the index. As reviewed in subsection 3.1 the resulting asymptotics allows making contact with the bulk BPS black holes. The asymptotics on the W wings has been an open problem. In this work we discover a host of interesting phenomena, most importantly partial deconfinement, by investigating the W wings.
limit ourselves for simplicity to p = q, in which case the formula takes the form I(p, p, y 1,2,3 ) ? = û∈eBAEs I(û; ∆ 1,2,3 , τ ), (1.6) for some rather involved special function I spelled out in subsection 4.1, involving the elliptic gamma function. Hereû ∈ eBAEs means that we have sum over the solutions to the elliptic Bethe Ansatz equations of the SU(N ) N = 4 theory spelled out in section 4.1, involving Jacobi theta functions. Prior to the present work, only a set of isolated solutions to the eBAEs were known. These were derived in [20], and we will refer to them as "the standard solutions". They are labeled by three non-negative integers {m, n, r} subject to mn = N , 0 ≤ r < n, and correspond to perfect tilings of the torus with modular parameter τ . In this work we discuss various new solutions, which we refer to as "non-standard". In particular, we will argue that for N > 2 there are continua of such solutions. The Bethe Ansatz formula for the index then breaks down, and needs to be reformulated to incorporate such continua; hence the question mark above the equal-sign in (1.6).
Despite the said shortcoming of the Bethe Ansatz approach for N > 2, we will still utilize it in section 4 by temporarily neglecting the continua of eBAE solutions. This way we study the CKKN limit (with b = 1) of the index, and will compare the result with that obtained in section 3 from the integral expression.
Up to the same caveat, we will also utilize the Bethe Ansatz formula in section 4 to study the large-N limit of the index. Taking the Cardy limit after the large-N limit, we obtain an answer that we will compare with the asymptotics of the index in the Cardy before large-N limit, analyzed in section 3.

Outline of the new technical results
For readers interested in specific technical results, here we provide a list of the new findings of the present paper, with reference to the appropriate section where they are discussed.
• Relation between the N = 1 * theory and the N = 4 theory.
The correspondence between vacua of the N = 1 * theory and solutions to the N = 4 eBAEs is spelled out in subsection 2.1.1. It leads to Conjecture 2 in section 4 stating that for N ≥ (l + 1)(l + 2)/2, there are l-complex-dimensional continua of solutions to the SU(N ) N = 4 eBAEs.
• The Cardy-like asymptotics of the index.
We studied the Cardy-like asymptotics of the index for generic ∆ 1,2 ∈ R and arg β = 0 using the elliptic hypergeometric integral form in section 3. For N = 2, 3, 4, it is presented for the first time in subsections 3.2.1 and 3.2.2. Some earlier studies had considered only the x ij = 0 saddle-point in the matrix-integral, which is incorrect on the W wings. For large-N , a lower bound for the index is obtained in (3.17), wherein C max can be taken to infinity. This lower bound along with Lemma 1 establishes that the index is partially deconfined all over the W wings in the double-scaling limit. This finding encourages Conjecture 1 that the said lower bound is actually optimal and therefore gives the large-N after the Cardy-like asymptotics of the index. Using this asymptotic expression for the index and assuming Q a ∈ CZ, we find critical points of the Legendre transform of (logarithm of) the index yielding micro-canonical entropies S C (J 1,2 , Q a ) = S BH (J 1,2 , Q a )/C for C = 2, 3, 4, 5; these presumably correspond to entropies of new (possibly multi-center) black objects in the bulk. (What happens to the micro-canonical entropy for C > 5 is not clear to us, as commented at the end of section 3.) We studied the same index using the Bethe Ansatz form in subsection 4.2. Compatibility with the results from the elliptic hypergeometric integral form implies the existence of new eBAE solutions that were not covered in [20]. We indeed found such solutions numerically (analytically in the Cardy-like limit) for some simple cases and reproduced the lower bound (3.17). Conjecture 1 would imply that in the Bethe Ansatz approach the other eBAE solutions, which we have not fully figured out, will not contribute to the leading Cardy-like asymptotics of the index at large N .
We discuss various new SU(N ) eBAE solutions that were not covered in [20], referred to as "non-standard" solutions. In subsection 4.3.1, we employ elementary elliptic function theory to establish the existence of one such solution (two if we count the different signs) for N = 2, and present its asymptotics. In subsection 4.3.2 we discuss numerical evidence that for N = 3 there is a one-complex-dimensional continuum of eBAE solutions. We further discuss this continuum in the low-and high-temperature limits. It turns out that a member of this continuum can be captured exactly (i.e. at finite temperature). This is thanks to the correspondence of subsection 2.1.1 with the N = 1 * theory, which allows us to borrow a result of Dorey [25]. We present analytic evidence that this exact non-standard solution is indeed a member of a one-complex-dimensional continuum. This is achieved via a beautiful three-term theta function identity presented as Lemma 2, which establishes that the associated Jacobian factor of the eBAE solution vanishes. Finally we discuss numerical evidence for Conjecture 2 for N ≤ 10: in particular, we found numerical evidence for two and three complex dimensional continua of solutions to the SU(N ) eBAEs for N = 6 and N = 10 respectively. These findings imply that the Bethe Ansatz formula for the 4d N = 4 index is valid in its currently available form only for N = 2; for higher N it needs to be reformulated to take the continua of Bethe roots into account.
• The large-N asymptotics of the index.
In subsection 4.4 we estimated the large-N limit of the index, extending previous results in [6]. In particular, the leading Cardy-like asymptotics of the improved large-N limit of the index turns out to match the large-N after the Cardy-like asymptotics of the index (3.19), with a couple of subtle issues discussed in details in the main text. This suggests that the asymptotic behavior of the index in the double-scaling limit is captured by Conjecture 1 regardless of the order of the Cardy-like limit and the large-N limit.
2 High-temperature phases of the 4d N = 4 index The 4d N = 4 index (just as any Romelsberger index [1] for that matter) can be computed as a partition function on a primary Hopf surface [26] with complex-structure moduli τ, σ. To gain intuition on these moduli, we work with b, β instead, defined through τ = iβb −1 /2π, σ = iβb/2π. When b, β are positive real numbers, the Hopf surface is S 3 b × S 1 , with a direct product metric, and with β = 2πr S 1 /r S 3 , while b becomes the squashing parameter of the three-sphere. Then, in analogy with thermal quantum physics, one can interpret the S 1 as the Euclidean time circle, and hence think of β as inverse-temperature in units of r S 3 .
More generally, the complex-structure moduli of the Hopf surface could be such that β becomes complex. Then the Hopf surface is still topologically S 3 × S 1 , but metrically it is not a direct product anymore. This situation would correspond to having a "complex temperature".
From a field theory perspective, allowing β to become complex simply amounts to extending the territory of exploration, with potentially new behaviors of the index to be discovered in the extended domain. (For example, as we will recollect in subsection 5.1, general Romelsberger indices seem to exhibit a much faster and much more universal Cardy-like growth in subsets of the complex-β domain [27,28].) From a holographic perspective, on the other hand, complexifying β finds a distinctly significant meaning through its relation with rotation in the bulk-c.f. [26,29]. This relation arises because the non-direct product geometry of the boundary can be filled in only with rotating spacetimes. This observation in turn explains why studies of the Cardy-like limit of the 4d N = 4 index prior to CKKN [3] found a much slower growth than that required by the bulk black holes: earlier studies had focused on real β, while the bulk BPS black holes have rotation and require complex β. (As we will recollect in subsection 5.1, a second important novelty of the limit studied by CKKN was considering complex y k .) In this work we study various "phases of the 4d N = 4 index at high (complex) temperatures". What we mean by this is as follows. First of all, since the Hopf surface the index corresponds to is compact, in order to have a notion of "phase" (associated to various "saddle-points" dominating the partition function) we need to take some limit of the index. When a large-N limit is taken, one can speak of high-(complex-) temperature phases of the index when |β| is small enough-smaller than some finite critical value for instance. For finite N on the other hand, to have a notion of phase we go to "infinite-(complex-) temperature" |β| → 0 (with | arg β| ∈ (0, π/2) fixed). We can then classify various behaviors of the index in those limits as various "phases". The control-parameters ∆ 1,2 in turn would often allow Roberge-Weiss type [16] transitions between such phases.
There appear to be two particularly natural classification schemes in the present context, and we now explain both in some detail. In particular, two different notions of partial deconfinement arise from the following classifications. When discussing partial deconfinement in the following sections, it should be clear from the context which of the two notions we are referring to. (A "sub-matrix deconfinement", different from partial deconfinement in the senses elaborated on below, has been recently discussed in [30][31][32].)

Classification via center symmetry
A first classification scheme arises if following in the footsteps of Polyakov [33] one considers patterns of center-symmetry breaking by the dominant holonomy configurations in the index.
In section 3 we will discuss various center-breaking patterns in the Cardy-like limit of the index. The Cardy-like limit is analogous to the infinite-temperature limit of thermal partition functions. We will speak of partial deconfinement in a sense similar to that of Polyakov, when the dominant holonomy configurations break the Z N center to a subgroup (possibly an approximate one for large N ). For C > 1 a divisor of N , a useful orderparameter for a single critical holonomy configuration is the C-th power of the Polyakov loop Tr P C x * = N j=1 e 2πiCx * j , which condenses (i.e. becomes nonzero) in the "C-center phase" where Z N → Z C . The dominant holonomy configurations in such a phase would then be C-centered, as they would correspond to C packs of N/C condensed (i.e. collided) holonomies distributed uniformly on the circle; summing over all critical holonomies then recovers the center symmetry as one would expect: x * Tr P C x * = 0. Partial deconfinement in this sense has been discussed earlier (see e.g. [34]) in the more conventional thermal, non-supersymmetric context with R 3 as the spatial manifold; there a sharper notion of an order-parameter exists, since on non-compact spatial manifolds cluster decomposition forbids the analog of summing over all x * .
The one-center phase would of course be the fully deconfined one where all the holonomies condense at a given value (either 0, or 1 N , ..., or N −1 N , due to the SU(N ) constraint). However, in principle this is not the only pattern for a full breaking of the Z N center. For example, a random distribution of the holonomies on the circle would also completely break the center. When the dominant holonomy configurations break the Z N center completely, but are not one-centered, we say the high-temperature phase of the index is "non-standard"; the Polyakov loop then may or may not condense. We will encounter such non-standard phases in the Cardy limit of the N = 4 index for N = 5, 6 in section 3; they correspond to the holes in Figure 4, and in the Bethe Ansatz approach they would arise when non-standard eBAE solutions take over the index in the Cardy limit.
In section 4 we will demonstrate how partial deconfinement in a similar sense can occur in the large-N limit of the index as well. The large-N limit will be analyzed for τ = σ via the Bethe Ansatz approach, where the behavior of the index depends on which solution of the eBAEs dominates the large-N limit. Such solutions can be thought of as complexified holonomy configurations, whose τ -independent parts correspond to the Polyakov loops. At finite β, they also have τ -dependent parts though, that are analogous to 't Hooft loops. There is hence also an analog of "magnetic" Z N center at finite β, which has in fact appeared in the Bethe Ansatz context already in [20]. Therefore the high-temperature phases at finite β and large N , can be classified via subgroups of Z N × Z N , in a picture that is in a sense dual to 't Hooft's classification of phases of SU(N ) gauge theories [35][36][37]. We now proceed to expand on this duality-or correspondence-below.

Correspondence with vacua of compactified N = 1 * theory
There is a regime of parameters where we expect close connection between high-temperature phases of the N = 4 index and low-energy phases of the N = 1 * theory on R 3 × S 1 . This is the regime where i) β → 0 and β ∈ R >0 , such that in the r S 1 → 0 "direct channel" one is probing high-temperature phases on S 3 × S 1 , while in the r S 3 → ∞ "crossed channel" one is probing low-energy phases on R 3 × S 1 ; ii) the chemical potentials ∆ 1,2,3 are small enough that their periodicity and balancing condition are not significant. Even then, the ∆ k are real masses for the adjoint chiral multiplets of compactified N = 4 theory, while N = 1 * theory has complex masses for its adjoint chirals. Nevertheless, one might expect at least some resemblance between high-temperature phases of the N = 4 index and low-energy phases of compactified N = 1 * theory based on the channel-crossing argument, and interestingly enough closer inspection reveals not just a resemblance, but a precise quantitative correspondence, aspects (though not all) of which extend even to finite complex β and arbitrary ∆ 1,2 ∈ R.
First, the proper identification between the complex-structure modulus τ of the N = 4 index and the complexified gauge couplingτ of the N = 1 * theory seems to be as follows: Alternatively, the electric and magnetic loops are swapped in the two pictures-i.e. the Polyakov loops in the direct channel correspond to the 't Hooft loops in the crossed channel and the 't Hooft loops in the direct channel correspond to the Wilson loops in the crossed channel. The identification (2.1) can be motivated through the crossed-channel relation between the two pictures, but a more satisfactory derivation of it would be desirable. Once the identification is accepted though, one can compare the vacua of compactified N = 1 * theory as determined via Dorey's elliptic superpotential [25], with the possible high-temperature phases of the N = 4 index as determined via solutions to the elliptic BAEs.
A particularly interesting aspect of the correspondence which survives at finite complex β and finite ∆ 1,2,3 , is the connection between the massive phases of the compactified N = 1 * , and the standard solutions to the N = 4 eBAEs: massive vacua ←→ standard eBAE solutions, (2.2) valid for arbitrary N . Specifically, the vacuum associated to the subgroup F r of Z N × Z N generated by (0, n), ( N n , r) in the Donagi-Witten terminology [38], corresponds to the ( N n , n, r) standard eBAE solution [20] spelled out in subsection 4.1 below. 2 Moreover, just as the massive phases are permuted via S duality in N = 1 * theory, the standard solutions to the N = 4 eBAEs are permuted via an SL(2, Z) acting on τ [20].
Most strikingly for our purposes, the Coulomb phase of the SU(3) N = 1 * theory corresponds to a continuous set of non-standard solutions to the SU(3) N = 4 eBAEs. This yields an exact non-standard SU(3) eBAE solution via the corresponding N = 1 * vacuum given by Dorey [25]. We will discuss this exact non-standard eBAE solution in section 4. We suspect that more generally for any N , for general complex τ , and at least for suitably chosen ∆ 1,2 , there is a correspondence Coulomb vacua ←→ continua of non-standard eBAE solutions, (2.3) though the precise map might quite non-trivially depend on the chosen ∆ 1,2 . A further bridge, due to Dorey [25], is expected to connect the vacua on R 3 × S 1 to those on R 4 . 3 Based on this correspondence and available knowledge (see e.g. [42]) on semi-classical Coulomb vacua on R 4 , we expect that for N ≥ (l + 1)(l + 2)/2 there are l-complex-dimensional continua of eBAE solutions for the SU(N ) N = 4 theory. We have numerically checked that this expectation pans out for N = 4 through 10 as well: the N = 4, 5 cases, just like for N = 3, contain one-complex-dimensional continua of non-standard eBAE solutions, while in the N = 6 case for the first time a two-complex-dimensional continuum of solutions arises. This persists for N = 7, 8, 9, and then a new three-complex-dimensional family of solutions appears at N = 10. These continua of solutions present a serious difficulty for the Bethe Ansatz formula for the 4d N = 4 index [18,19], which is derived assuming only isolated eBAE solutions. We will comment more on this point in section 4.
An example of phenomena arising for finite complex τ or large ∆ 1,2 that are outside the regime of validity of the correspondence is presented by the isolated non-standard SU(2) eBAE solution u ∆ discussed in subsection 4.3. For ∆ 1,2 so large that the lower branch of (4.32) becomes relevant, even in the β → 0 limit the non-standard solution u ∆ does not correspond to an SU(2) N = 1 * vacuum.

Classification via asymptotic growth
Deconfinement in a second sense can be associated with the asymptotic growth of Re log I (or more generally Re log of a partition function, with the Casimir-energy piece removed), either as |β| → 0 for finite N , or as N → ∞. This is a more general sense as it does not depend on a center symmetry. In the case of the 4d N = 4 index, an O(1/|β| 2 ) growth as |β| → 0, or an O(N 2 ) growth as N → ∞ [9] could count as "deconfinement" in this second sense.
So far these criteria do not distinguish between the fully-deconfined and partially-deconfined phases, as classified in the first sense. A more refined classification is possible in the present context though, due to the presence of the first-order Roberge-Weiss type transitions [16]. Let us consider the specific case of the SU(4) N = 4 theory in the Cardy-like limit as an illustrative example. In Figure 3 we see three high-temperature "phases" of the index separated by first-order transitions. Each phase can be labeled according to its fastest growth, say by its s := sup( lim |β|→0 |β| 2 Re log I), (2.4) which is proportional to the maximum height of its corresponding curve in Figure 3. Then we get an ordering of the phases with s > 0: although no longer a fully-deconfined or a confined phase in an absolute sense, we still get a "maximally deconfined" and (possibly) a "minimally deconfined" phase in a relative sense, as well as other "partially deconfined" phases in the middle. Phases with s ≤ 0 might more appropriately be called "non-deconfined". Figure 3 displays a clear correlation between this and the previous sense of deconfinement: the curves corresponding to larger center-breaking have higher maxima and therefore faster maximal growth. Also, note that the blue curve which would correspond to a confined (i.e. center-preserving) phase in the Polyakov sense, is associated to a "minimally deconfining" phase in the sense of asymptotic growth. On the other hand, the non-standard phases arising for N ≥ 5, exhibit intermediate asymptotic growth, so would be partially deconfined in the second sense, even though their dominant holonomy configurations break the center completely.
In the large-N limit, besides a similar ordering of the deconfined phases viã there is also a useful notion of a "confined" phase [9] where Re log I = O(N 0 ). In similar problems, there could of course be various other phases with intermediate scaling as well.

Cardy-like asymptotics of the index
Following [5], the elliptic gamma functions can be expanded in the CKKN Cardy-like limit (1.5), so that the index (1.2) simplifies as with the integral over the N − 1 independent holonomies corresponding to the maximal torus of SU (N ). Here Note in particular that κ(x) is compatible with the unit periodicity of the holonomies. Note also that on the RHS of (3.2) every ∆ 3 can be replaced with −∆ 1 − ∆ 2 ; this is because the balancing condition y 1 y 2 y 3 = pq fixes ∆ 3 = τ + σ − ∆ 1 − ∆ 2 mod Z, and since we are interested in the leading Cardy-like asymptotics, we can neglect τ and σ in ∆ 3 . (To capture the subleading effects a generalization of κ(x) to complex domain is needed [5]-c.f. (4.15) below.) Since we are interested in the |τ σ| → 0 limit, the integral in (3.1) is dominated by the global maximum of the real part of the exponent, or alternatively the global minimum of Re(iQ h /τ σ). Moreover, this limit is well defined since κ(x) is continuous and bounded and the integration domain is compact. As a result, the leading asymptotic behavior of the index is given by where x * is the holonomy configuration corresponding to the global minimum 4 . Taking the parametrization τ = iβb −1 2π and σ = iβb 2π , and noting that Q h is a real function, we see that x * corresponds to the global minimum of (3.5) From the x-dependent part of Q h we see that minimizing V eff is equivalent to minimizing a potential of the form 1≤i<j≤N V Q (x ij ; arg β, ∆ 1,2 ), with the pairwise part explicitly reading (3.6) Figure 1 above shows the qualitative behavior of this pairwise potential.

Behavior on the M wings
As Figure 1 shows, on the M wings the pairwise potential (3.6) is minimized at x ij = 0. 5 Consequently the overall potential V eff also takes on its global minimum when all holonomies are identical. Taking SU (N ) into account, there are N possible configurations, namely all These configurations are dominant in the Cardy-like limit and completely break the Z N center. Here we have complete deconfinement [5] (see also [3,4]), as the index exhibits "maximal" asymptotic growth where the sign should be taken to be the same as that of arg β. (More precisely, it is after appropriate tuning of ∆ 1,2 that the maximal asymptotic growth is achieved in this fully deconfined phase; see Figure 2.) This "grand-canonical" asymptotics in (3.7), when translated to the micro-canonical ensemble, yields the expected entropy S BH of the bulk AdS 5 black holes [15,43]. (The original work [15] showed this last statement for the minus sign, and [43] later established it for the plus sign as well.) Our focus in this work is on the W wings though, to which we now turn. 4 If there are degenerate minima, x * can be taken to correspond to any one of them, as the added degeneracy factor is subleading in the CKKN limit. 5 This was found in [5] by numerically scanning the space of the control-parameters. In [4] an analytic proof was suggested for an M-type behavior all over the parameter-space when arg β > 0; however, as pointed out in the Added Note of [5], the proof actually applies only to the upper-right wing of Figure 1, and the oddity of the potential under ∆1,2 → −∆1,2 establishes in fact the W-type behavior on the lower-left wing when arg β > 0. The two wings of course switch places for arg β < 0.

Behavior on the W wings
The W wings are characterized by the feature that the minimum of the pairwise potential (3.6) is displaced away from x ij = 0. In particular, it is located either at x ij = 1/2 mod 1 or in a flat region around this point. The issue now is that, except for special case of SU (2), it is impossible to center the differences x ij around 1/2 mod 1 for all i and j. As a result, the global minimum of V eff cannot correspond to the individual minima of all the individual pairwise potentials, and the extremization problem then becomes quite challenging. For this reason, it has been an open problem to find the Cardy-like asymptotics of the index on the W-wings (c.f. Problem 1 in Section 5 of [5]). In this section we completely address the problem for N ≤ 4, and take steps towards addressing it for N > 4.

SU(2) and infinite-temperature confinement/deconfinement transition
The W-wing behavior is easy to determine for the SU (2) case, as the minimum at x 12 = 1/2 along with the SU (2) condition x 1 + x 2 = 0 is trivially solved by the "confining" holonomy configuration x 1 = −x 2 = 1/4. This leads to the W wing asymptotics Note that i) the chosen dominant holonomy configuration on the W wings respects the Z 2 center symmetry generated by x i → x i + 1/2, and ii) as already discussed in [5], the fastest asymptotic growth of the index on the W wings is slower than the fastest asymptotic growth on the M wings, as expected.

SU(N ) for finite N > 2
For N > 2 it is no longer possible to have all x ij equal to 1/2, and finding the global minimum of the effective potential becomes a difficult problem. Nevertheless, we can obtain lower bounds on the asymptotic growth by examining special sets of holonomy configurations.
In particular, we note that for any set of holonomies x 0 on the maximal torus. This bound is saturated when x 0 = x * , but is suboptimal otherwise. Our goal is then to pick a family of configurations {x 0,i } and optimize over this family. Because the potential is exponentiated, we can in fact write which is a convenient way to package the lower bound on the asymptotic growth. The choice of holonomy configurations to optimize over will of course determine how optimal the bound will be. As a compromise between simplicity and robustness of the estimate, we consider the family based on grouping the N holonomies into packs of d collided There are a total of N/d distinct packs, and they are then distributed uniformly on the periodic interval [−1/2, 1/2] in such a way that they satisfy the SU (N ) condition j x j ∈ Z. This latter condition gives rise to d discrete configurations (which we collectively denote by x d ), signalling a partial breaking Z N → Z N/d of the center. These special configurations also arise as the hyperbolic (or "high-temperature") reduction of the set of eigenvalue configurations found in [20]; we will comment more on this point below.
For a given divisor d, since there are C := N/d packs distributed evenly on the periodic unit interval, the spacing between packs is 1/C. As a result, the configuration x d yields .
Here the second term is the contribution of the d(d − 1)/2 collided holonomy pairs inside a single pack, with x ij = 0, and the third term is the contribution of the holonomy pairs between the different packs, with x ij = J/C. To simplify (3.11) further, we use the remarkable identity which can be derived with Mathematica's aid. Applying this to (3.11) and substituting in d = N/C then gives The symbol emphasizes that the RHS is only a lower bound on the asymptotic growth of I in the CKKN limit. While we are mainly interested in the W wings, derivation of this bound is independent of the wings. In particular, this bound is optimal on the M wings, where the optimal term is given by C = 1, corresponding to the condensation of all holonomies into a single pack. On the W wings, however, except for N = 2 where it reduces to (3.8), the bound (3.13) is not necessarily optimal, as we will argue below. Although the bound (3.13) is written as a sum over 'trial' configurations, generically, depending on where exactly we are on the W wings, only one term would dominate the sum. The question of which divisor C provides the strongest bound then boils down to the comparison of 1 for various divisors C of N . For N a prime number, there are only two divisors, namely C = 1 and C = N . In this case, the answer is simple: the "confined" C = N term dominates the Figure 2: The functions C −3 3 a=1 κ(C∆ a ) for C = 1 (brown), C = 2 (green), C = 3 (yellow), and C = 6 (blue). sum in (3.13) on the W wings, while the "fully deconfined" C = 1 term dominates on the M wings. For composite N , however, other divisors (besides 1 and N ) may give the dominant contribution to the sum in (3.13) on the W wings. For example, for N = 6, Figure 2 shows that while on the M wing the fully condensed C = 1 term is dominant, on the W wing there are regions where the other divisors take over the sum in (3.13).
As already emphasized, the sum in (3.13) gives us a lower bound, but not necessarily the true asymptotic growth of the index. Nevertheless, we conjecture that at least on subsets of the regions where various divisors become dominant in (3.13), the corresponding term in the sum actually gives the true asymptotics. In other words, that for any finite N there are confining or partially deconfining phases on the W wings.
For small values of N we can attack the extremization problem numerically of course, and make more precise statements. This is what we will do next. For example, for N = 6 we establish that on a subset of the region in Figure 2 where the green curve takes over, the dominant holonomy configuration is indeed the partially-deconfining (Z 6 → Z 2 ) configuration x d=3 , that on a subset of the "yellow region" the dominant holonomy configuration is the partially-deconfining (Z 6 → Z 3 ) configuration x d=2 , and that on a subset of the "blue region" the dominant holonomy configuration is the confining (Z 6 -symmetric) configuration x d=1 .
A remarkable surprise of the numerical investigation discussed below is that for N = 3, 4 the bound (3.13) is in fact optimal! Therefore (3.13) gives the exact leading-order asymptotics of the index in these cases. For larger values of N on the other hand, the numerical analysis shows that there are indeed regions on the W wings where the bound (3.13) is not optimal.

SU(3): infinite-temperature confinement/deconfinement transition
Our numerical investigation shows that for N = 3 the bound (3.13) is optimal to within two parts in 10 15 -which is essentially the machine precision. We therefore conclude that the exact leading asymptotics of the index in this case reads Hence, just as in the SU(2) case, we have infinite-temperature confinement/deconfinement transitions moving from the W wings to the M wings. On the M wings the first term on the RHS of (3.15) dominates, and on the W wings the second term. Before moving on to richer cases, once again we emphasize that in the present paper we are studying the asymptotics on generic points of the parameter-space. On non-generic points where ∆ a ∈ Z or τ, σ ∈ iR + , the asymptotic growth would be slower, and a more involved analysis is required; c.f. section 3 of [5].

SU(4): infinite-temperature partial deconfinement
In this case as well, the numerical investigation shows that the bound (3.13) is optimal to within two parts in 10 15 . We therefore conclude that the exact leading asymptotics of the index for N = 4 reads The first and the third terms on the RHS come respectively from the fully-deconfined (C = 1) and the confined (C = N = 4) holonomy configurations. But here we have also the first instance of infinite-temperature partial deconfinement in the superconformal index: the middle term on the RHS of (3.16) takes over on the middle triangle of the W wings as shown in Figure 3. This term corresponds to C = 2, and signals a Z 4 → Z 2 breaking of the center symmetry in the Cardy-like limit. This qualifies as a partially deconfined phase, not only because of its partial-breaking pattern of the center symmetry, but also because of its "partial liberation of the constituents" as signified by the fact that the height of the green curve lies between those of the blue (C = 4, confined) curve and the brown (C = 1, fully-deconfined) curve.

SU(5) and SU(6): insufficiency of the divisor configurations
In these cases the numerical analysis shows that there are regions on the W wings where none of the divisor configurations x d minimizes V eff . Fixing arg β > 0 for concreteness, we see from Figure 4 that in the SU(5) case there is a relatively large such region, but for SU (6) there are rather small subsets of the W wing where this happens. Hence the bound (3.13) Figure 3: The functions C −3 3 a=1 κ(C∆ a ) for C = 1 (brown), C = 2 (green), and C = 4 (blue). The take-over of the green curve signifies the partially deconfined phase in that region when N = 4. seems much more efficient in the SU(6) case. This is to be expected of course, as there are three contributing trial configurations (C = 2, 3, 6) on the W wings when N = 6, while there is only one such configuration (C = 5) when N = 5.

Taking the large-N limit
The bound on the asymptotic growth of the index, (3.13), was derived using a family of holonomy configurations based on divisors C of N . This bound can of course be improved by enlarging the family of trial configurations. One way to do this is to divide the N holonomies into C collided packs with the packs evenly distributed on the periodic interval for all integer C = 1, 2, . . . , N . In general, each pack cannot have the same number of holonomies unless C is a divisor of N . Nevertheless, we can make the packs nearly uniform by first distributing N/C holonomies into each of the C packs. This leaves N mod C holonomies left over, which can then be distributed in some prescribed manner in the packs. This set of trial configurations would in principle improve the bound given in (3.13). However, the resulting bound would be sensitive to the particular distribution of the left over N mod C holonomies, and can no longer be expressed in such a compact manner.
Although the refined bound that is obtained by splitting the eigenvalues into C packs Figure 4: The difference (scaled by a factor of 12) between the numerically maximized Q h , and the Q h maximized over the divisor configurations, on the arg β > 0 W wing, for N = 5 on the left, and for N = 6 on the right. When the result is zero, it means the divisor configurations are maximizing Q h (hence minimizing V eff ). Note the big hole in the middle for N = 5, and the small holes for N = 6, signalling the failure of the divisor configurations to maximize Q h (hence to minimize V eff ).
for all integers C does not admit a simple expression for finite N , it nevertheless simplifies in the large-N limit, at least for the leading order growth of the index. The idea here is that, instead of taking C = 1, 2, . . . , N , we cut off the set of trial configurations at some large but finite C max that is independent of N . For a given C, we then start with C packs of N/C holonomies and compute V eff for this subset of C N/C holonomies. This is of course incomplete, but we can add in the remaining pairwise potentials, (3.6), between these uniform holonomies and the N mod C remaining ones (as well as those among the remaining holonomies themselves). These interactions between O(C) objects and O(N ) objects (as well as those among the remaining O(C) objects) add a correction of at most O(N ) since we keep the cutoff on C fixed. Alternatively, starting with C N/C instead of N holonomies also leads to a correction of the same order. As a result, the leading O(N 2 ) behavior of the index is captured by (3.13) with the modification that the sum is taken over all integers up to the cutoff: Since we have dropped terms of O(N ) or smaller, this asymptotic bound is only valid when considering the O(N 2 /τ σ) growth of the index. In this case, in fact because the bound applies for any finite C max ∈ N, we can remove the C max cutoff and instead take the sum to infinity. Note that the large-N bound, (3.17), confirms that the finite N bound, (3.13), is not optimal in general! For N a large prime, for example, the finite N bound would consist of a sum over only the C = 1 and C = N terms, with the C = N (or "confining") term winning on the W wings. But we know (c.f. Figure 2) that for large enough N , at least on subsets of the W wings, the C = 2, 3, . . . terms in the large-N bound (3.17) dominate over the confining term. This is of course a simple result of enlarging the set of trial configurations to include more general C collided packs of holonomies for all integer C, whether C divides N or not.
Returning to the large-N analysis, we see that as long as at least one term in (3.17) has a positive real part in the exponent, the index will exhibit O(N 2 ) growth in the Cardy-like limit. This corresponds to either full deconfinement when the C = 1 term dominates, or partial deconfinement when some C > 1 term dominates. As discussed above, the C = 1 term always dominates in the M wings, even at finite N , with the resulting behavior given by (3.7). On the other hand, the situation is more elaborate in the W wings. Let us fix arg β < 0 for concreteness; then the W wing consists of all ∆ 1,2 subject to 0 The question then becomes whether for any such ∆ 1,2 we can find a C ∈ N such that a κ(C∆ a ) < 0. We now argue that this is the case! In fact since κ(C∆ a ) is periodic under ∆ 1,2 → ∆ 1,2 + 1/C, we can simply focus on the square 0 < ∆ 1,2 < 1/C. Now, it follows from the scaling ∆ a → ∆ a /C that on this square the sign of a κ(C∆ a ) is positive (resp. negative) if the representatives {C∆ 1,2 }/C of ∆ 1,2 on the square 0 < ∆ 1,2 < 1/C lie on the lower triangle with vertices (0, 0), (0, 1/C), (1/C, 0) (resp. the upper triangle with vertices (0, 1/C), (1/C, 0), (1/C, 1/C)). Hence the question boils down to whether we can find a C such that the representatives are on the upper triangle where {C∆ 1 }/C + {C∆ 2 }/C > 1/C. (The interested reader might find that a simple drawing of the said triangles would make the previous sentences obvious.) The following lemma answers this question in the positive.
Lemma 1 For every pair of real numbers x, y subject to 0 < x, y, 1 − x − y < 1, there exists a natural number C > 1 such that {Cx} + {Cy} > 1.
An elementary proof of this lemma can be found in the appendix. 6 Here we instead point out that it follows from a much stronger result, often associated 7 to the names Kronecker and Weyl, that if there is no integer relationship between x, y (i.e. no solution to ax + by + c = 0 in integers other than (a, b, c) = (0, 0, 0)) then the points ({Cx}, {Cy}) are dense in the unit torus. (In our case, even if there is such an integer relationship between ∆ 1 , ∆ 2 , we can always establish our desired result by applying the Kronecker-Weyl theorem to ∆ 1 , ∆ 2 + , with a small-enough chosen such that there is no integer relationship between ∆ 1 , ∆ 2 + .) Similar arguments apply when arg β > 0. We thus conclude that for all points strictly inside the W wings, there exists a natural number C > 1 such that the exponent of the "C-th 6 We learned the proof, as well as the following remark regarding the Kronecker-Weyl theorem, from David E Speyer, a mathematician at University of Michigan. 7 See https://mathoverflow.net/questions/162875/reference-for-kronecker-weyl-theorem-in-full-generality.
bound" in (3.17) has positive real part, and hence the index is partially deconfined. A "non-deconfined" behavior (i.e. o(N 2 )/τ σ growth for log I as N → ∞ after the Cardylike limit) might appear in the non-generic situations where arg β = 0 (c.f. section 3 of [5]), or ∆ a ∈ Z. In such cases, subdominant terms of O(N ) or smaller may be important in order to fully pin down the behavior of the index.
With some optimism, this genericity of partial deconfinement on the W wings can be taken as a sign that it would be consistent to conjecture that the large-N bound (3.17), with the cut-off removed, gives not just a lower bound but the actual leading asymptotics of the index.

Conjecture 1
The leading asymptotics of the superconformal index (1.2) of the 4d N = 4 theory with SU(N ) gauge group, in the CKKN limit (1.5), simplifies as N → ∞ to with the error such that logarithms of the two sides differ by o(N 2 /τ σ).
This conjecture is motivated in part by the following two observations: i) that in the N → ∞ limit there are infinitely many trial configurations, and hence increasing chance of their sufficiency; ii) that already for N as small as 6, as witnessed by Figure 4, the divisor configurations go a long way towards minimizing V eff on the W wings.

Entropy of the partially deconfined phases
Let us now focus on the C-th term in (3.19), and see what entropy it bears. 8 First, we rewrite the C-th term explicitly as where we have defined C∆ a as (3.21) The corresponding entropy S C (J 1,2 , Q a ) is obtained by performing a Legendre transform of (3.20), which requires adding −2πi(σJ 1 + τ J 2 + Σ a ∆ a Q a ) in the exponent, and then extremizing. The first step can be written explicitly as where we have replaced C∆ a with C∆ a in the last term; this replacement is allowed assuming Q a ∈ CZ, which we do for simplicity, because then since the difference between C∆ a and C∆ a is an integer it would not change the exponential when multiplied by −2πiQ a /C. Extremizing the functionŜ C (J 1,2 , Q a ; ∆ a , σ, τ ) in (3.22) with respect to the chemical potentials ∆ a , σ, τ under the constraint Σ a ∆ a − σ − τ ∈ Z determines the entropy as where ∆ C a , σ C , and τ C denote the critical points under the aforementioned constraint. With an appropriate relation between J 1,2 , Q a the resulting entropy will be a real number [3], and hence acceptable.
Note that the expression inside the round bracket in (3.22) with general C reduces to the expression with C = 1 under C∆ a → ∆ a , Cσ → σ, and Cτ → τ . In other words, this simple replacement maps the present problem to that of the C = 1 (single-center) black hole entropy problem. We thus find S C (J 1,2 , Q a ) = S BH (J 1,2 , Q a )/C, (3.24) as alluded to in section 1.
Using the same property, we can also figure out the critical points with general C directly from the known ones with C = 1. Since our analysis is restricted to real ∆ 1,2 , however, we should focus on the equal-charge case where all Q a 's are equal to each other (c.f. the end of section 2 in [5]). In that case, the critical points with C = 1 have been already known as with the signs the same as that of arg β. The critical points with general C are then determined as where j, k are arbitrary integers. For C = 2 as an example, fixing arg(β) > 0 for concreteness, we have the critical point at (∆ 1 , ∆ 2 ) = (− 1 3 , − 1 3 ) on the W-wing. The interested reader is encouraged to locate the C = 3 critical points in Figure 2.
An interesting question is whether at the critical point the C-th term in (3.17) is indeed dominant; otherwise the entropy derivation would not be self-consistent. Curiously, a numerical investigation shows that the answer is positive for C ≤ 5, and negative for general C > 5. The interpretation of this result is not yet clear to us.

Comparison with the Bethe Ansatz type approach
It was argued in [18,19] that the index (1.1) can be rewritten as a Bethe Ansatz type formula. One advantage of this reformulation is that the integral over the Coulomb branch in (1.2) is replaced by a sum over solutions to a set of Bethe Ansatz like equations. This was the approach used in [6] to obtain the black hole microstate counting in the large-N limit. Here we briefly review the Bethe Ansatz approach (BA approach) to the index and then demonstrate how the partially deconfined phases identified in the previous section emerge in this approach.

The Bethe Ansatz type expression for the index
For simplicity, we restrict to p = q (i.e. τ = σ = iβ 2π ) in the index. In this case, the Bethe Ansatz type formula reads [6,19] (1−e 2πikτ ) 2(N −1) andû = {u 1 , · · · , u N } denotes all possible solutions to the following system of elliptic Bethe Ansatz equations (eBAEs), 2) under the SU (N ) constraint Σ N i=1 u i ∈ Z + τ Z and the abbreviation u ij = u i − u j . The third chemical potential ∆ 3 is constrained via ∆ 3 = 2τ −∆ 1 −∆ 2 (mod Z) as in the previous section. Note that λ is a free parameter independent of i. Here we have introduced Z(û; ∆ a , τ ) and H(û; ∆ a , τ ) as , (4.3b) and the elliptic functions θ 0 (u; τ ) andΓ(u; τ, σ) are defined as (z = e 2πiu , q = e 2πiτ , p = e 2πiσ ) Γ(u; σ, τ ) = Γ(z; p, q). Note that H(û; ∆ a , τ ) has to be evaluated at the solutions to the eBAEs after taking the partial derivatives of Q i 's with respect to u i 's. Since the right-hand side of (4.1) is summed over all possible solutions to the eBAEs (4.2), the first step towards the computation of the index (4.1) is to find the most general solutions to the eBAEs (4.2). Here we find it convenient to use the relation (1 − e 2πikτ )θ 0 (u; τ ), (4.6) to rewrite the eBAEs (4.2) in terms of the Jacobi theta function θ 1 (u; τ ) as These eBAEs are a set of highly non-linear equations, and it seems rather challenging to find the most general solutions. However, in the form (4.7), these equations coincide with those for the topologically twisted index of N = 4 SYM on T 2 × S 2 [44]. In [20] oddity and quasi-periodicity of θ 1 (u; τ ) with respect to the first argument were used to find a large set of solutions to (4.7). These solutions are denoted in terms of three non-negative integers {m, n, r} with N = mn and r = 0, 1, . . . , n − 1, and have the u i 's distributed aŝ u {m,n,r} = uĵk =ū + nĵ + rk N +k n τ N = mn, 0 ≤ r < n, (ĵ,k) ∈ Z m × Z n , (4.9) whereū is determined by the SU (N ) constraint ĵ k uĵk ∈ Z + τ Z. In essence, these {m, n, r} solutions correspond to regular distributions of the N holonomies over the fundamental domain of the torus specified by (1, τ ). We will refer to these as the standard solutions to the eBAEs (4.2) and refer to the u i 's as holonomies in analogy to the holonomies x i in the integral representation of the index. It turns out, however, that the standard solutions, (4.9), are in fact not the most general solutions to the eBAEs. In a way, this is not particularly surprising because of the non-linear nature of the equations. We will refer to the additional solutions that do not fall into the class of (4.9) as non-standard solutions. Such solutions will not correspond to a periodic tiling of the fundamental domain and moreover may depend on the chemical potentials ∆ k . The Bethe Ansatz form for the index is then a sum over standard and non-standard solutions, which we denote schematically as I(q, q, y 1,2,3 ) = (4.10) Note that the contribution from the standard solutions to the index, dubbed as I {m,n,r} (∆ a , τ ) in (4.10), can be written explicitly as by substituting (4.9) into (4.1, 4.3a).

The Cardy-like limit of the index
In Section 3, we were able to bound the asymptotic growth of the index in the Cardy-like limit based on a set of trial holonomy configurations with the N holonomies distributed among C packs of collided holonomies. For finite N , we took C to be a divisor of N so that all packs contain an equal number N/C of holonomies and arrived at the bound (3.13). This bound can be improved in the large-N limit by removing the requirement that N/C is an integer; the O(N 2 ) behavior of the index is then governed by (3.17). In this subsection we reproduce the same results in the BA approach.

Standard solutions and the asymptotic bound (3.13)
At first it may not be obvious how (3.13) can be obtained in the BA approach. However, the connection can be made by taking the τ → 0 limit of the standard solutions for the holonomies, given in (4.9). They reduce in this hyperbolic (or "high-temperature") limit tô This distributes the N holonomies in the unit interval spaced in multiples of 1/N . In fact, it is not difficult to see that periodicity of the holonomies ensures that (4.12) is equivalent to the N holonomies distributed evenly into C distinct packs of N/C collided holonomies with C = N gcd(n, r) .
(4.13) (In the special case r = 0 we define gcd(n, 0) := n.) This indeed corresponds directly to the holonomy configurations considered in Section 3 which led to the asymptotic bound, (3.13), for the finite-N index in the Cardy-like limit. To derive (3.13) in the BA approach more explicitly, we compute the contributions from the standard solutions to the index, namely I {m,n,r} (∆ a , τ ) given in (4.11), in the Cardy-like limit. The result, whose derivation we now outline, confirms that they saturate the asymptotic bound in (3.13). First, since we are interested in the 1/|τ | 2 leading behavior of log I, the contribution log α N (τ ) from the prefactor is negligible. Next, the Jacobian determinant H(û {m,n,r} ; ∆ a , τ ) involves derivatives of theta functions which are of order O(e −1/|τ | ) [20] in general 9 and therefore − log H(û {m,n,r} ; ∆ a , τ ) does not contribute to the 1/|τ | 2 leading order of log I. Hence the leading contribution to log I comes entirely from log Z {m,n,r} (∆ a , τ ), which can be estimated using the following asymptotic formula for the elliptic gamma function: where κ τ (u) is defined as This κ τ (x) is a generalization of κ(x) introduced in (3.3) to complex domain [5]. Note that {x} τ and κ τ (x) reduce repectively to {x} and κ(x) for x ∈ R. The 1/|τ | 2 -leading order behavior of the index is then obtained from the asymptotic expansion of (4.3a) where we made use of the relation κ τ (u) + κ τ (−u) = 0. The connection to the integral approach to the index is now manifest, as this reproduces Q h (x; ∆) in (3.2) obtained in the CKKN limit of (1.2). Of course, here, the holonomies are not integrated over, but rather are taken to satisfy the eBAEs, (4.2) or equivalently (4.7). In the hyperbolic limit, the standard solutions all reduce to (4.12), which is equivalent to dividing them into C distinct packs of collided holonomies. Substituting any one of these solutions into (4.16) then necessarily gives the same result that was previously obtained in the CKKN limit. In particular, taking ∆ 1 , ∆ 2 ∈ R and noting that κ τ (x) reduces to κ(x) for real arguments, we obtain where C = N/ gcd(n, r). In the Cardy-like limit, the Bethe Ansatz form of the index, (4.10), then takes the form where we have made explicit the degeneracy factor d(C) counting the number of distinct standard solutions that gives rise to a given C, though it of course does not contribute to the 1/|τ | 2 leading order. If we discard the contributions of any possible non-standard solutions, then this saturates the asymptotic bound obtained earlier as (3.13). Recall that, in the integral approach, the bound was obtained by taking a set of trial configurations corresponding to C equal packs of collided holonomies. Here, the Bethe Ansatz approach uses the same set of holonomy configurations, so it is no surprise that the final asymptotic expression for the index is the same. However, in this approach, the holonomy configurations are not just trial configurations, but are exact solutions to the eBAEs, and remain exact solutions even away from the Cardy-like limit.

Non-standard solutions and the improved asymptotic bound (3.17)
As demonstrated in the CKKN limit of the elliptic hypergeometric integral, the basic asymptotic bound, (3.13), can be improved by enlarging the set of trial configurations to encompass C nearly uniform packs of collided holonomies for all integers C, not necessarily dividing N .
In the BA approach, however, we cannot arbitrarily choose trial configurations; rather we are limited to solutions of the eBAEs. Since the standard solutions only allow for values of C that divide N , they cannot generate the improved bound, (3.17). This strongly suggests that the set of standard solutions is in fact incomplete and that we need non-standard solutions resembling arbitrary C packs where C does not divide N .
As an example, consider the case C = 2. When N is even, this is realized by a standard solution, which corresponds in the Cardy-like limit to N/2 holonomies taking the value u i = 0 and another N/2 holonomies taking the value u i = 1/2 (up to a constantū shift enforcing the SU (N ) condition). Away from the Cardy-like limit, this splits into two eBAE solutions, the first corresponding to {m, n, r} = {2, N/2, 0}, and the second to {1, N, N/2}. In both cases, the N/2 holonomies in each pack are evenly distributed along the periodic τ direction, although in the second solution the two packs are offset by τ /N while in the first solution they are not.
The more interesting case is how C = 2 is realized when N is odd. The expectation here is that there must be a non-standard solution where the holonomies are split into two packs. Although we cannot equally divide an odd number of holonomies, we can imagine grouping (N + 1)/2 of them into one pack and (N − 1)/2 of them in the other. While the first pack has one additional holonomy, this should become unimportant in the large-N limit. Nevertheless, we demonstrate that such a non-standard solution exists, even for finite N .
In appendix C we establish a non-standard solution in the Cardy-like limit with the chemical potentials satisfying ∆ 1 + ∆ 2 ≤ 1 2 . It is given explicitly aŝ This asymptotic non-standard solution satisfies the eBAEs in the Cardy-like limit (as displayed in (C.22)) up to exponentially suppressed terms. Note that this solution corresponds to both packs having holonomies equally spaced along the periodic τ direction and hence can be viewed as a natural odd-N version of the standard C = 2 solution. Although this solution only satisfies the eBAEs asymptotically, we have demonstrated that similar solutions continue to exist at finite τ by numerically solving the exact eBAEs (4.7). As an example, we present a numerical solution with N = 11 in Figure 5a. This strongly implies that there is an exact non-standard solution to the eBAEs (4.7), whose asymptotic form is given as (4.19). Unlike the standard solutions, however, this non-standard C = 2 solution does not have the holonomies uniformly distributed on the torus, and moreover generally depends on the potentials ∆ a , except in the Cardy-like limit.
It is now straightforward to insert the asymptotic solution, (4.19), into (4.16) to obtain the contribution to the index which can be compared with the standard C = 2 solution obtained by taking C = 2 in (4.17). Although these expressions still demonstrate an even/odd effect at finite N , the leading O(N 2 /τ 2 ) behavior is identical in the large-N limit.
We now turn to the second possibility, where the chemical potentials satisfy ∆ 1 + ∆ 2 > 1 2 . Here the non-standard solution takes the asymptotic form

(4.22)
This corresponds to two isolated holonomies along with packs of (N − 3)/2 and (N − 1)/2 holonomies, respectively. Numerically, we can see that, for sufficiently large τ , the first two holonomies are actually part of a single pack of (N + 1)/2 holonomies, just as in case 1. However, as Im τ approaches zero, this pair first moves towards the real axis and then finally towards ±∆ 1 /2 in the Cardy-like limit. The exponentially small deviation away from this endpoint is given by provided the quantity in the parentheses is positive. This is true for a generic value of ∆ 1 + ∆ 2 > 1 2 (not asymptotically close to 1 2 ) and a sufficiently large N . Although (4.22) only holds asymptotically, numerical solutions of this type can be obtained for finite τ , and we provide an example in Figure 5b.
Inserting the non-standard solution (4.22) into the asymptotic expression (4.16) then gives the contribution to the index

(4.24)
Although this expression is more complicated than the corresponding one for Case 1, it reduces in the large-N limit to the same leading order behavior as all the other C = 2 cases, namely with C = 2. While we have focused on non-standard solutions for C = 2, numerical investigations confirm that similar non-standard solutions exist for arbitrary values of C and N , at least for N sufficiently large so that it can be divided into nearly equal packs of holonomies. These non-standard solutions are similar to that with C = 2 in that they are sensitive to the choice of chemical potentials ∆ 1 and ∆ 2 , with the simpler configuration, corresponding to Case 1 above, occurring only when ∆ 1 + ∆ 2 ≤ 1/C. In this case, the non-standard solution in the Cardy-like limit is given bŷ (4.26) where N = C N/C + D (D = 1, · · · , C − 1). This satisfies the eBAEs in the Cardy-like limit (as displayed in (C.22)) up to exponentially suppressed terms. In principle, this solution can be inserted into (4.16) to obtain the contribution of a given C non-standard solution to the index. The resulting expressions simplify in the large-N limit, and reduce to (4.25) as expected.
Although this extension of the Case 1 solution to arbitrary values of C only holds for sufficiently small ∆ 1 + ∆ 2 , we expect that generalizations of the Case 2 solution (where pairs of holonomies may be pulled away from the main packs) exist for other values of the chemical potentials. We thus conjecture that, at least in the Cardy-like limit, solutions to the eBAEs exist for all values of C and N with N/C > d. Here, d is a number of order unity corresponding to the minimum number of holonomies in a single pack that allows the solution to be categorized as a set of packs instead of individually distributed holonomies. When C divides N , the solution is standard, but otherwise it is non-standard. For each value of C, the contribution to the index then takes the form (4.25), regardless of whether the solution is standard or non-standard. As a result, the Bethe-ansatz type approach to the index reproduces the improved asymptotic bound (3.17) found above.

Additional non-standard solutions
As we have seen explicitly, the non-standard solutions (4.19) and (4.22) to the eBAEs (4.7) for C = 2 and odd N contribute to the index in much the same way as the standard {2, N/2, 0} solution for even N , at least in the large-N limit. Moreover, such contributions were crucial to reproduce the improved asymptotic bounds (3.17) with arbitrary integer C that may or may not divide N . This makes us conclude that we cannot ignore the non-standard solutions in the Bethe ansatz approach to the index.
While we have focused on non-standard solutions that are similar to the {C, N/C, 0} solutions since all standard solutions reduce to this form in the Cardy-like limit, we have also found numerical evidence for the existence of non-standard solutions at finite τ . This suggests that many if not all standard {m, n, r} solutions have generalized counterparts as non-standard solutions. In addition, there may be additional non-standard solutions that do not fall into any particular classification. Thus it would be nice to understand the full set of solutions to the eBAEs.
In principle, we would just solve the eBAEs (4.7) to obtain a complete set of standard and non-standard solutions. However, in practice, this is a rather challenging problem, even in the asymptotic limits. Therefore, to make the problem tractable, we focus on the N = 2 and N = 3 cases. Even so, much of the analysis is rather technical, so we relegate the details to Appendix C and only highlight the main results here.

Non-standard solutions for N = 2
Since the eBAEs (4.7) depend only on the difference u ij = u i − u j of the holonomies, and since the free parameter λ is unconstrained, the N = 2 eBAEs reduce to a single equation for a single variable u 21 Here for convenience we have used the real chemical potentials∆ 1,2,3 , which are essentially the same as ∆ 1 , ∆ 2 , −∆ 1 − ∆ 2 , up to simple shifts and reflections such that See Appendix C for more details.
Since each theta function has a first-order zero, this fraction has six zeros. Furthermore, since it is elliptic, the fraction takes all values, including unity, six times, and thus the eBAE has six solutions. Four of them are the familiar ones It seems challenging to find u ∆ in closed form, except in the asymptotic limits, so part of the investigation will be numerical. We denote the 'low-temperature' limit to be |τ | 1 and the 'high-temperature' limit to be |τ | 1 where arg τ ( = 0, π) is held fixed for both cases. In either limit, the Jacobi theta function θ 1 (u; τ ) has a straightforward asymptotic expansion, so the equation (4.27) can be directly solved. The asymptotic analysis is presented in Appendix C, and we summarize the N = 2 results here. Figure 6: Numerical plots of the non-standard N = 2 solution u ∆ with arg(τ ) = π/4. The figure on the left corresponds to∆ 3 < 1/2 while that on the right corresponds to∆ 3 > 1/2. Note that the vertical axis is given in units of (complex) τ .
As indicated above, the standard solutions, (4.29), are exact solutions for any∆ a and for all τ , so we only focus on the non-standard solution specified by u ∆ . In the low-temperature limit the solution-given in (C.11)-takes the form The high-temperature asymptotic solution splits into two cases, depending on whether∆ 3 is less than or greater than 1/2. The solution is given by (C. 26) and (C.30), which can be summarized as (4.32) There does not appear to be a simple analytic solution away from these asymptotic limits. However, numerical investigations demonstrate that the non-standard solutions (4.31) and (4.32) are indeed continuously connected. Figure 6 shows explicit examples of the numerical solution for u ∆ with both possibilities of∆ 3 .
To be complete, it should be noted that for∆ 3 ≥ 1/2, the high-temperature limit also admits a continuous set of solutions which interestingly correspond to the plateaus of Q h in the integral approach of the previous section. However, these solutions do not appear to be continuously connected to any solutions away from the high-temperature limit, so we believe they are only an artifact of the hightemperature asymptotics, and are not true solutions to the eBAEs once subleading terms are included.

Non-standard solutions for finite N > 2
The structure of the eBAEs is considerably harder to analyze for N > 2 as there are more holonomy pairs u ij to consider. For any N , note that the product of the N individual eBAEs (4.7) gives the condition e 2πiN λ = 1 so that e 2πiλ is a root of unity (not necessarily primitive). Focusing next on a single equation, say i = 1, then gives the condition (4.34) For N = 3, this means we can solve for u 31 in terms of u 21 where ω = e 2πi/3 is a primitive cube root of unity and k = 0, 1, 2. We thus look for nonstandard solutions by picking a given u 21  The two roots are related by the map u ij → −u ij , corresponding to taking z ij → 1/z ij . The important feature of this solution is that the eBAEs reduce to a single condition on two complex holonomies. As a result, we end up with a continuous family of non-standard solutions.
To further explore the nature of the non-standard solutions, we consider, for simplicity, the case where all the chemical potentials are identical to each other, namely∆ a = 1/3. In this case, the right-hand side of (4.36) vanishes, and we find simply Since the second case can be obtained from the first by taking u i → −u i , we focus on the first case. The expression 1 + z 21 + z 31 = 0 has a simple interpretation in the complex plane as a triangle with sides 1, z 21 and z 31 . Since the solution is restricted to |q| ≤ |z 21 | ≤ 1 and |q| ≤ |z 31 | ≤ 1, the space of solutions is given by the intersection of two disks of radius one centered at 0 and 1, respectively. The continuous family of non-standard N = 3 solutions is shown schematically in Figure 7. The analysis of the high-temperature asymptotic solutions is rather involved, so we again focus on the case∆ a = 1/3. In addition to the standard solutions, the continuous family of non-standard solutions survives in the high temperature limit, and is given by (C.34) in Appendix C, which we rewrite in terms of u 21 = x 21 + y 21 τ and u 31 = x 31 + y 31 τ u 21 , u 32 ∈ u 21 , u 32 0 ≤ x 21 mod 1 < 1 3 , u 21 − 2u 31 = Z + 1 3 τ (Z + 1 2 ) ∪ u 21 , u 32 This indicates that, up to discrete shifts, there is a one-complex dimensional family of hightemperature solutions. Although the family appears to be split into four branches, the first and last subset in (4.38) is part of the same branch, as shown in Figure 8.
We have verified numerically in several examples that the asymptotic low and high temperature solutions are continuously connected with each other. This suggests that the continuous family of non-standard solutions is in fact generic for any value of τ . In general, the which is∆ a -independent. This is an exact non-standard solution, whose validity can be checked using (4.8). We were actually led to it via the correspondence of subsection 2.1.1 with N = 1 * theory: it corresponds to a Coulomb vacuum of the compactified N = 1 * theory with SU(3) gauge group, listed in Dorey's paper [25] as the "fifth" vacuum.
Despite its simple form somewhat reminiscent of the standard solutions, the solution (4.39) is part of the continuous family, and not isolated. To see this, we look for vanishing eigenvalues of the Jacobian matrix (∂Q i /∂u j ) where the Q i represent the eBAE expressions as given in (4.7). In general, there are N eBAEs along with N holonomies u j . However, the constraint i Q i = 1 along with the SU (N ) constraint j u j ∈ Z + τ Z effectively reduces this to an (N − 1) × (N − 1) matrix.
Lemma 2 For any τ in the upper-half plane, and any complex∆ 1,2,3 subject to a∆ a ∈ Z, we have This lemma, combined with (4.42), establishes that the Jacobian matrix indeed has a vanishing eigenvalue. This establishes the existence of a "zero mode" taking (4.39) to nearby eBAE solutions. Proof 10 of this lemma can be found in appendix B.
Although we have yet to perform an analytic investigation of the non-standard solutions for N > 3, we have investigated some cases numerically. For N = 4 and 5, in addition to the isolated standard solutions, we find solutions on which ∂Q i /∂u j has a single vanishing eigenvalue, suggesting that they are part of one-complex dimensional families of solutions. More interestingly, for N = 6, 7, 8, and 9 we find evidence for both one and two complex dimensional families of solutions, the latter signalled by two vanishing eigenvalues for ∂Q i /∂u j . Similarly for N = 10 we find evidence for one, two, and three complex dimensional families of solutions. These are in agreement with the expectation that the continuous families of solutions correspond to Coulomb vacua of N = 1 * theory. Based on these numerical results and the putative correspondence with N = 1 * theory, we conjecture that at least the dimensionality of the space of the N = 4 eBAE solutions is captured correctly by the semi-classical formula for the rank of the Coulomb vacua in SU(N ) N = 1 * theory.
The existence of such flat directions is somewhat unusual, nevertheless, in that it necessarily demands the existence of non-trivial identities among products of Jacobi theta functions.
Of particular significance is that the existence of continuous solutions to the eBAEs poses a serious issue to the full validity of the Bethe ansatz type approach as established in Ref. [19]. The reason is that the rewriting of the index, (1.2), into the Bethe ansatz form relies on the fact that the general solutions to the eBAEs are isolated so that Cauchy's formula may be applied [19]. Hence, if there are continuous families of non-standard solutions, the Bethe ansatz type approach must be reformulated to incorporate such continuous Bethe roots into the sum over solutions to the eBAEs, (4.10).

The large-N limit of the index revisited
While we have been motivated by the Cardy like limit of the index, the appearance of partially deconfined phases as well as non-standard solutions suggests that we re-examine the large-N limit at finite τ as investigated in Ref. [6]. In this large-N limit, one or at most a few solutions (if there are any degeneracies) to the BAE would be expected to dominate in the sum over solutions, (4.10). The focus is then on identifying the dominant solution, much as one would search for a dominant saddle point contribution.
The main emphasis in [6] was on the basic solution, corresponding to I {1,N,0} (∆ a , τ ). In the large-N limit, the contribution was found to be where Θ(x a , τ ) is defined as where {x} τ and κ τ (x) are defined in (4.15). Using the constraint ∆ 3 = 2τ − ∆ 1 − ∆ 2 , we can rewrite Θ(∆ a , τ ) explicitly as and therefore (4.44) is in fact a function of ∆ 1 , ∆ 2 and τ only. While it is clear that the basic solution yields a contribution of (N 2 ), there are potentially many other sources of contributions at the same order. These include other standard {m, n, r} solutions as well as isolated and continuous families of non-standard solutions. The conjecture of Ref. [6] is that only the set of T -transformed solutions,û {1,N,r} , would yield additional contributions at O(N 2 ). The resulting large-N index is then argued to take the form log I(p, p, y 1,2,3 ) = max log I {1,N,r} (∆ a , τ ) : (N log N ), (4.47) unless the ∆ a 's are located along Stokes lines where the asymptotic expansions of the elliptic functions fail to converge or where different contributions may compete with each other. We note that there is some tension between the conjectured large-N index, (4.47), and our conjectured leading asymptotics in the CKKN limit, namely (3.19). In particular, the latter, as well as the more rigorous bound, (3.17), include configurations where the holonomies are split into C nearly equal packs distributed evenly on the circle, and this is not visible in (4.47). Of course, it is possible that the Cardy like limit and the large-N limit do not commute. Nevertheless, this motivates us to ask whether additional solutions to the eBAEs, either standard or non-standard, may contribute at large N .

A new parametrization of the standard solutions
In the Cardy like limit, we have seen that standard solutions obtained from C identical packs of N/C collided holonomies will contribute at O τ −2 (N 2 /C 3 ) . Here we are assuming that C divides N , so this is a standard solution. If the limit is smooth, then we expect this behavior to persist even as we move away from the Cardy limit. As in (4.13), the standard {m, n, r} solutions that reduce to C packs of holonomies have gcd(n, r) = N/C. This can be parametrized as the set of standard solutions Hereâ labels the pack andb labels a particular holonomy within that pack. In the large-N limit, provided C ∼ O(1), the holonomies become dense along the τ circle, and hence the shift between packs by a fraction s/N does not qualitatively change the appearance of the solution. This suggests that we can compute I {C,s} in the large-N limit in a universal manner, provided C ∼ O(1) so that N/C remains large. To do so, we start with the expression for the standard solution given in (4.11). Ignoring the subleading contributions from the prefactor α N (τ ) and the Jacobian H(û {C,s} ; ∆ a , τ ) −1 as in [6] then gives log I {C,s} (∆ a , τ ) where the primes in the summation symbols indicate that the (â,b) = (â ,b ) case is excluded. This can be rewritten as where on the first line theb =b case is excluded. In fact, the first line is similar to the expression (4.44) for the basic {1, N, 0} solution, but with N replaced by N/C and with an extra prefactor of C. Note that the replacement N → N/C in the large-N formula is valid since we have assumed C ∼ O(1). The second line can be simplified using the formulas (for ∆ = 0) [6], along with (4.44). The result is provided that ∆ a is not an integer multiple of 1/C. In particular, this demonstrates O(N 2 ) scaling of the standard {C, s} solution for C ∼ O(1), even away from the the Cardy-like limit. Moreover, this leading-order behavior is independent of the offset s along the τ circle between adjacent packs of holonomies, thus confirming the intuition that once the packs are sufficiently dense, the distribution of holonomies within each pack becomes unimportant, at least at leading order. The analysis of the large-N limit of the standard solutions is still incomplete, as we have not considered the case where C scales with N . Nevertheless, the contributions (4.54) suggests that the conjecture (4.47) must be refined by enlarging the holonomy configurations that need to be considered. In terms of the {C, s} labels, the basic {1, N, 0} solution is included as {1, 0}, while the T-transformed {1, N, r} solutions fall under {C, s} with r = qN/C where q and C are relatively prime. Note that this presents a bit of a puzzle as (4.54) and (4.47) appear to be in conflict except for r = 0.
The resolution of this puzzle is based on two observations. The first is that, for r ∼ O(1) in the large-N limit, we must have C ∼ N , in which case the expression (4.54) breaks down. The second is that if instead C ∼ O(1) then r must be large. However, in this case, the map from the first to the second line of (4.47) breaks down as it is only valid for r scaling as O(N 0 ). This is because the derivation of Θ(∆ a , τ + r) from the large-N asymptotics of the standard solution (4.11) with {m, n, r} = {1, N, r} was based on making the shift τ → τ + r. However, shifts with r = O(N ) may lead to a modification to log I {1,N,r} (∆ a , τ ) that is not captured by simply shifting τ in (4.44).
Provided we only take standard solutions into account, the large-N limit of the index will be given by maximizing over all {C, s} solutions. Although we are unable to provide an analytic expression for C ∼ O(N ), we nevertheless expect the large-N index to take a form along the rough lines of log I(q, q, y 1,2,3 ) = max where C divides N for standard solutions. The shift of τ → τ + r in the above allows us to include T-transformed versions of the {C, s} solutions, and the restriction of this expression to C = 1 reproduces the conjecture (4.47) of [6].
Of course, this is still expected to be incomplete, as we have not yet accounted for the non-standard solutions. We know from the Cardy limit that isolated non-standard solutions exist for all integer values of C. Numerically, these solutions are continuously connected to low temperature solutions, and hence should exist for arbitrary τ . Again, in the large-N limit, the isolated non-standard solutions correspond to C packs of either N/C or N/C + 1 holonomies distributed along the τ circle. Since we take C ∼ O (1), the difference between packs shows up as a O(1/N ) correction, and hence we expect the leading O(N 2 ) behavior to be correctly captured by (4.55) where we now drop the restriction that C divides N .
Finally, we are able to connect the large-N limit of the index, (4.55), to the Cardy like limit by taking τ → 0. Looking only at the leading O(1/τ 2 ) order, it is easy to see that all terms with r = 0 will not contribute. Making use of along with the identity (3.12) then gives log I(q, q, y 1,2,3 ) = max (4.57) in perfect agreement with the conjecture (3.19) for the leading large-N asymptotics in the Cardy like limit. This certainly suggests that the asymptotic behavior of the index is independent of the order of limits between the Cardy-like limit and the large-N limit. Of course, in order to remove the conjectures and be more rigorous in the Bethe Ansatz approach, we have to consider solutions with C ∼ O(N ), as well as the issue of continuous non-standard solutions.

Summary and relation to previous work
Several papers investigating asymptotic growth of supersymmetric indices of the 4d N = 4 theory have appeared in the last year. Below we outline how our findings complement those of the most closely related recent work.
• Cardy-like asymptotics: Choi-Kim-Kim-Nahmgoong (CKKN) [3], Honda [4], and Ardehali [5]. All three of these papers investigated the Cardy-like limit of the 4d N = 4 index I(p, q; y k ) using its integral representation, and in the limit where y k approach the unit circle (for y k not approaching the unit circle the problem is still open-see Problem 2 in [5] and subsection 5.2 below). CKKN [3] took actually also a large-N limit-after the Cardy-like limit-to simplify the analysis, but in [4,5] N was left arbitrary. CKKN [3] assumed that the dominant holonomy configurations in the Cardylike limit correspond to equal holonomies x ij = 0; it was later realized in [4,5] that this assumption fails in essentially half of the parameter-space (see the Added Note of [5] for the relation between the findings of [4] and [5]). As in [5] we divide the parameter-space into complementary M wings and W wings. While on the M wings the dominant holonomy configurations indeed correspond to x ij = 0 and the asymptotics has been understood, on the W wings finding the Cardy-like asymptotics has been an open problem-see Problem 1 of [5]. In section 3 we solved this problem for N = 2, 3, 4, and also conjectured the formula (3.19) in the large-N limit. In particular, for N = 4, as well as in the large-N limit, we have discovered regions on the W wings corresponding to partially deconfined phases in the Cardy-like limit of the 4d N = 4 index.
• Large-N asymptotics: Benini-Milan [6]. This reference studied the large-N limit of the 4d N = 4 index using its Bethe-Ansatz representation. In the present paper we pointed out some difficulties with the Bethe-Ansatz representation for N ≥ 3: there seem to be continuous families of Bethe roots, calling for an integration measure which is so far not understood. Setting this difficulty aside, we pointed out that some of the discrete Bethe roots (found in [20]) are not negligible in the large-N limit of the index as assumed in [6]. We moreover discovered new discrete Bethe roots that play an important role in the large-N asymptotics of the index. The new non-negligible contributions that we have found demonstrate partially deconfined phases in the large-N limit of the index as well.
• The density-distribution approach to the large-N limit: CKKN [7]. In this reference, following the original work of Kinney-Maldacena-Minwalla-Raju [2], the large-N limit of the index was analyzed using its integral representation, by rewriting the integral in terms of ρ n := N j=1 e 2πinx j , rather than x j . These new variables are Fourier coefficients of the density distribution ρ(x) for the holonomies in the large-N limit. The drawback of this approach is that the range of ρ n does not appear to be simple to derive. This difficulty hinders an accurate analysis of the saddle-points of the integral over the ρ n variables. We have thus avoided this approach in the present work.
• Cardy-like asymptotics on higher Riemann sheets: Cabo Bizet-Cassani-Martelli-Murthy [28] and Kim-Kim-Song [27]. Since the fugacities of the index satisfy y 1 y 2 y 3 = pq, one can "turn off the flavor fugacities" by setting y k = (pq) 1/3 ; the resulting function I(p, q; (pq) 1/3 ) is the usual N = 1 superconformal index of the N = 4 theory, and is known not to have fast asymptotic growth in the usual Cardy-like limit p, q → 1, due to bose-fermi cancelations. [More precisely, for SU(N ) N = 4 theory it can be shown from the results in [23] that I(p, q; (pq) 1/3 ) grows like β N −1 , not like e 1/β 2 as expected from the bulk black holes.] These statements are not the end of the story however, as they really apply only to the fundamental Riemann sheet of the function. For p, q inside the unit open punctured disc, the function I(p, q; (pq) 1/3 ), unlike the more well-behaved I(p, q; y k ), is not meromorphic; therefore besides turning on flavor fugacities there is another way to get fast growth from it-or "obstruct its bose-fermi cancelations" if you will-and that is to go to its higher Riemann sheets. Because of the power 1/3 for pq, the non-meromorphic N = 1 index I(p, q; (pq) 1/3 ) has in fact three inequivalent sheets, which can be identified along branch cuts at arg p = π, and be labeled by n = −1, 0, +1. The papers [27,28] show that on the n = ±1 sheets it is possible to get the fast exponential growth associated to the bulk black holes-simply note that p → pe 2πin introduces nontrivial phases in (pq) 1/3 , equivalent to "turning on flavor fugacities" and setting y k = (pq) 1/3 e 2πin/3 on the n = 0 sheet, corresponding to ∆ a = n/3 in the CKKN limit, which can obstruct the bose-fermi cancelations as familiar from CKKN's work [3]. Now, while for sign(arg β) = ± the Cardy-like asymptotics of the index has been obtained from the fully-deconfining holonomy configurations on the n = ±1 sheets, on the other (n = ∓1) sheets we expect that the partially deconfined configurations become significant! Note that here the control-parameters are n, sign(arg β). In particular, for the SU (4) N = 4 theory it follows from our results in section 3 that it is the partiallydeconfining Z 4 → Z 2 holonomy configurations which take over on the n = ∓1 sheets.
Similarly, if our conjecture (3.19) is correct, in the large-N limit the partially-deconfining Z N → Z N/2 holonomy configurations take over on the n = ∓1 sheets.
• Higher Riemann sheets with path-integration: Cabo Bizet-Cassani-Martelli-Murthy [43]. This reference is also related to analytic continuation of the N = 1 index to the n = 0 sheets, and is in fact the pioneering work on investigation of such higher sheets. However, it adopts a Lagrangian (path-integral) approach, whose analytic continuation is not fully understood. In particular, the path-integral computation in [43] suffers from subtleties regarding regularization of the analytically continued supersymmetric Casimir energy: as discussed in footnote 13 of [43], already before analytic continuation the regularization scheme used there does not give the correct result when applied to general theories, and it is somewhat of a coincidence that it works for the N = 4 theory. There seems to be no reason to believe that the scheme keeps being coincidentally correct in the analytically continued case, which is relevant to black hole counting.
[Such subtleties can be overcome when working with the path-integral version of I(p, q; y k ) in the CKKN limit though, which would be relevant to the analysis here; see section 4 of [5] where this issue was addressed.] • The topologically-twisted 4d N = 4 index: Hosseini-Nedelin-Zaffaroni [44] and Hong-Liu [20]. Besides the superconformal index, there is the rich and interesting topologically-twisted index that one can compute exactly for the N = 4 theory. This index also admits a Bethe-Ansatz expression with precisely the same elliptic Bethe-Ansatz equations that feature in the 4d N = 4 superconformal index. Our results in this paper thus imply that the expressions analyzed in [44] and [20] too, for rank≥ 2, suffer from the difficulty of continuously connected Bethe roots. Furthermore, in some regions of the parameter-space this index exhibits a well-understood, fast asymptotic growth of type e 1/β in the Cardy-like limit, which is associated to black strings in AdS 5 [44].
In the rest of the parameter-space finding the asymptotics is more challenging however [20]. In those regions, we expect that the new discrete solutions that we introduced in this paper are of significance in the Cardy-like asymptotics of the topologically twisted index.

Future directions
We conclude by outlining some of the interesting open problems and exciting prospects related to this work.
• The Bethe Ansatz approach. There are several important open questions regarding the Bethe Ansatz approach. Already at rank 1 (i.e. for SU(2)) where the Bethe Ansatz formula seems valid, not all the solutions of the eBAE are known; as discussed in section 4 besides the four standard solution [20] 0, 1/2, τ /2, (1 + τ )/2, there are two more solutions of the form ±u ∆ , and it would be nice to find u ∆ in closed form.
For rank ≥ 2 we have given numerical evidence that there is a continuously connected set of solutions, undermining the Bethe Ansatz formula in its current form as a finite sum over eBAE vacua; it would be nice to have a proof for existence of uncountably many Bethe roots for rank ≥ 2. It would also be important to find ways of reformulating the BA approach to incorporate the continuous set of Bethe roots into account.
It would moreover be necessary for black hole-counting applications to characterize all the Bethe roots contributing to the leading asymptotics of the index in the large-N limit.
• General Cardy-like asymptotics and black holes with unequal charges. As remarked above, the Cardy-like asymptotics of the 4d N = 4 index for general complex ∆ 1,2 is still an open problem. This problem is relevant to the microstate counting of the bulk black holes with unequal charges, because it is only for the equal-charge black holes that ∆ 1,2 ∈ R as we have assumed [5]. While it might be quite challenging to address this problem for arbitrary N , in the SU(2) case this does not seem out of reach: finding the asymptotic form of u ∆ would lead to the desired asymptotics using the Bethe Ansatz formula.
• The gravity dual of the partially deconfined phases. The bulk duals of the partially-deconfined phases have not been constructed as far as we are aware. They should be new black objects and their accurate interpretation is not clear to us at this point. It is tempting to speculate that they might be interpreted as multi-center black holes, which in cases with lens-space horizon topology are sometimes referred to as black lenses [45,46].
• The black hole operators. One of the most exciting aspects of the recent advances in AdS 5 /CFT 4 microstate counting is the prospect they open for explicit construction of the operators dual to the bulk BPS black holes. The operators dual to the bulk KK particles have long been known of course: they are the multi-trace operators-nicely reviewed in Table 7 of [47] for instance. The black hole operators on the other hand, have been elusive even in AdS 3 /CFT 2 . We hope that the emerging refined understanding of the asymptotic growth of the 4d N = 4 index can guide future searches for these long-sought operators in the N = 4 theory.
of Energy under grant DE-SC0007859. JH is supported in part by a Leinweber Graduate Fellowship from the University of Michigan.

A Proof of Lemma 1
Cover the torus R 2 /Z 2 with balls of radius ε/2.
As discussed in the main text, the lemma proves that the SU(3) eBAEs have a "zero mode" at the exact non-standard solution (4.39).

C The elliptic Bethe Ansatz equations in the asymptotic regions
In this appendix, we investigate the eBAEs (4.7) in the asymptotic regions, including both the low-temperature (|τ | 1) and the high-temperature (|τ | 1) limits. Then we look for asymptotic, non-standard solutions for N = 2 and N = 3.

(C.2)
This is to avoid using a complex ∆ 3 = 2τ − ∆ 1 − ∆ 2 and to keep all the chemical potentials ∆ a real. We now use quasi-periodicity and the oddness of θ 1 (u; τ ), namely (4.8), to rearrange the chemical potentials∆ a and the holonomies u k for convenience when writing down the asymptotic expansions. Note that we always assume ∆ 1,2 ∈ R, or equivalently∆ a ∈ R.
We start with shifting∆ a by integers using (4.8) so that 0 ≤ Re∆ a < 1 is satisfied. This yields a∆ a ∈ {0, 1, 2}. Then since the eBAEs (C.1) are invariant under∆ a → 1 −∆ a and λ → −λ due to (4.8), redefining∆ a as∆ a → 1 −∆ a does not change the eBAE solutions {u i }. Therefore, we can redefine∆ a as∆ a → 1 −∆ a whenever a∆ a = 2 and consequently we have a∆ a ∈ {0, 1}. Finally, we assume∆ a ∈ R, then 0 ≤ Re∆ a < 1 together with a∆ a = 0 leads to∆ a = 0, which is forbidden for the index to converge. So we have 0 <∆ a < 1, The final expression of∆ a shifted from (C.2) is given in terms of ∆ a as where as usual {·} = · − · , and we have assumed ∆ 1,2 ∈ R \ Z. We can also specify the range of holonomies u k using the same properties of θ 1 (u; τ ) given in (4.8). The key is that the eBAEs (C.1) are invariant under u k → u k + p k + q k τ for arbitrary integers p k , q k ∈ Z. Using this invariance, we can set u k as u k = x k + y k τ (0 ≤ x k , y k < 1), (C.5) with y i ≤ y j for i ≤ j without loss of generality. Following the setup (C.3) and (C.5), now we investigate the eBAEs (C.1) in asymptotic regions and look for non-standard solutions there.

C.1 Low-temperature asymptotic solutions
We start with the low-temperature (|τ | 1) asymptotic region. First we rewrite the infinite product form of θ 1 (u; τ ) (4.6) as (C.9) These are the eBAEs reduced under the low-temperature limit, which yield N −1 independent algebraic equations. They are still involved to solve in general, however, so we consider simple cases: N = 2 and N = 3. where ω = e 2πi/3 , as well as a continuous family of solutions satisfying a single condition (1 + z 21 + z 31 ) 1 + 1 z 21 + 1 z 31 = 3 + 2 a cos 2π∆ a , (C. 16) which reduces to the finitely many non-standard solutions for N = 2 given in (C.11) under z 31 = z 21 . Starting with the discrete solutions, the first three are based on cube roots of unity and correspond to the high temperature limit of the standard solutions. In particular, the first entry in (C.15) corresponds to the C = 1 case, and the next two correspond to the C = 3 cases. The final solution in (C.15) is novel, as it is a non-standard solution that is nevertheless independent of the chemical potentials∆ a . Note that, strictly speaking, this solution violates the assumption made above that |z 31 | > |q|. However, it is in fact easy to see that this extends to an exact non-standard solution (u 21 , u 31 ) = (1/2, τ /2), (C. 17) for arbitrary τ . Perhaps more interestingly, we have found a continuous family of solutions given by (C. 16). This one complex dimensional family of solutions is obviously non-standard since it depends on the chemical potentials∆ a . In addition, the limiting form of (C.16) for z 21 → −1 and z 31 → 0 encompasses the non-standard solution (C.17). This suggests that (C.17) is not a discrete solution but rather a part of the continuous family. We confirm this explicitly in Appendix B by demonstrating that the transformation matrix ∂Q i /∂u j (i, j = 1, 2) where Q i represents the i-th BAE contains a zero mode corresponding to a flat direction.
(C. 21) where the curly bracket denotes mod Z as in (3.3). Under the high-temperature limit, the contributions from k ≥ 2 are exponentially suppressed and therefore (C.21) reduces to ) .

(C.22)
These are the eBAEs reduced under the high-temperature limit, which yield N −1 independent algebraic equations. They are still involved to solve in general, however, so we consider simple cases: N = 2 and N = 3.  It is still difficult to solve (C.31) for general chemical potentials and therefore we assume that all the chemical potentials∆ a are identical as∆ a = 1/3. Then the complete set of solutions to (C.31) is given as