Residues, modularity, and the Cardy limit of the 4d $\mathcal{N}=4$ superconformal index

We compute the superconformal index of the $\mathcal{N}=4$ $SU(N)$ Yang-Mills theory through a residue calculation. The method is similar in spirit to the Bethe Ansatz formalism, except that all poles are explicitly known, and we do not require specialization of any of the chemical potentials. Our expression for the index allows us to revisit the Cardy limit using modular properties of four-dimensional supersymmetric partition functions. We find that all residues contribute at leading order in the Cardy limit. In a specific region of flavour chemical potential space, close to the two unrefined points, in fact all residues contribute universally. These universal residues precisely agree with the entropy functions of the asymptotically AdS$_5$ black hole and its"twin saddle"respectively. Finally, we discuss how our formula is suited to study the implications of four-dimensional modularity for the index beyond the Cardy limit.


Introduction
Recently, there has been renewed interest in the superconformal index of the fourdimensional N = 4 SU (N ) Yang-Mills theory, which captures the degeneracies (with signs) of the 1/16 th BPS spectrum. This interest is caused by the fact that, in contrast to the conclusions of an earlier investigation [1], the index is able to reproduce the entropy of supersymmetric black holes with AdS 5 asymptotics [2][3][4][5]. More precisely, the first series of works [6][7][8] aimed to reproduce a specific "entropy function" or free energy from the index, which was previously shown in [9] to lead to the black hole entropy upon a Legendre transformation. Apart from this success, the index is a rather rich mathematical object and further study has among other things suggested the existence of unknown gravitational saddles [8,10,11]. There are also indications of an interesting structure relating these various saddles, perhaps akin to the SL(2, Z) family of BTZ black holes in three dimensions [12][13][14].
However, the comparison with AdS 3 /CFT 2 is not obvious. In particular, this structure arises from modularity of the CFT 2 partition function defined on a two-torus. The present work was originally geared toward exploring analogous structures in a four-dimensional context. Specifically, it has been known for a while that certain fourdimensional supersymmetric partition functions, including the superconformal index, obey certain SL(3, Z) properties, as emphasized in [15][16][17]. These properties can be explained from the fact that the manifolds on which such partition functions are defined can all be viewed as a gluing of two solid three tori along their boundary, identified up to an SL(3, Z) transformation. The associated holomorphic block factorization [18,19] is central to the SL(3, Z) modular property, which was explained in full generality in [20]. However, it turns out that applying the general framework of [20] to the problem of interest requires a different approach to the computation of the superconformal index as compared with existing approaches. The main reason for this is that the previous works, which we briefly summarize below, all lack a certain generality. This fact prohibits the application of the main modular property of [20] to the existing expressions for the index. In the present work, we will be mainly concerned with the computation of the index through residues and a subsequent Cardy limit via the SL(3, Z) modular property. We postpone the investigation of the aforementioned structure to a forthcoming publication [21].
Let us now turn to a summary of modern existing approaches to the computation of the index. In the process, we will also indicate their respective drawbacks for our purposes. Finally, we will briefly mention how our method is able to avoid these specific drawbacks.
1. In [6], supersymmetric localization is employed to compute the index. It is shown that in the black hole background, a supersymmetric Casimir energy-like 1 prefactor precisely reproduces the required free energy. The remaining index is argued to not contribute at leading order in N in the relevant parameter regime. The main disadvantage from our perspective is the fact that certain aspects of the gravitational solution are used to obtain the result. This makes it difficult to see how this method will lead to insight into unknown gravitational solutions.
2. Various authors have considered a Cardy-like limit of the partition function [7,10,11,[24][25][26][27][28]. This limit sends the chemical potentials τ and σ, associated with the angular momenta, to zero. It is analogous to the high temperature limit of the torus partition function in two-dimensional CFT, hence the name. In this limit, the computation of the index simplifies significantly, although one still has to rely on a saddle point approximation of the gauge integral to obtain the final result.
It was pointed out in [20] that there is an interesting connection between this Cardy limit and a certain modular property of the index. However, this modular property cannot be used for all values of the gauge parameters, and therefore its use has to be justified a posteriori by the precise value of the saddle. We believe that the implications of modularity will be more transparent by first computing the gauge integral (exactly) and only then applying the modular property. In this way, one can a priori justify its use and the effect of the Cardy limit is not obscured by a subsequent saddle point approximation.
3. The Bethe Ansatz method [29][30][31] allows one to compute the gauge integral exactly and can be used to evaluate the large-N limit of the index [8,32]. Using this method, the index can be expressed as a sum over "Bethe vacua" and is reminiscent of a sum over saddles. The Bethe Ansatz method was the first to point out a rich phase structure exhibited by the large-N limit of the index, by showing that in different regimes of parameter space different Bethe vacua provide the dominant contribution to the index. However, there are some unsatisfactory aspects such as the fact that not all Bethe vacua are known (see however [26,27]). Also, one requires a certain specialization of the τ and σ chemical potentials. More precisely, the initial work [8] required τ = σ, which was later generalized to τ /σ ∈ Q [32]. But even for the more general setting, the modular property of [20] does not hold. 2 4. Another line of work has attempted a saddle point approximation of the gauge integral in the large-N limit, as opposed to the Cardy limit [11,34]. 3 Due to an elliptic extension of the gauge parameters, it is possible to find a large class of saddles with interesting interrelations. In addition, this analysis makes it possible to study the phase structure of the 1/16th BPS sector and predicts the existence of new gravitational saddles. The main drawback of this approach is the apparent lack of an a priori justification of the elliptic extension. Indeed, such ellipticity is naturally associated with the presence of large gauge transformations along two non-trivial one-cycles in the geometry, whereas S 3 × S 1 only has a single nontrivial one-cycle. In addition, this approach so far has required the specialization τ = σ, which as we already mentioned is problematic for the use of the full modular property.
5. Most recently, there have also been numerical efforts to study the superconformal index [36,37]. These methods both require specializations of the parameters and truncations, which obscure modular properties.
The main issues we have pointed out with the existing methods can be summarized as: specialization of parameters and/or working at the level of the gauge integrand. We obtain another, exact representation of the index similar in structure to the Bethe Ansatz analysis. In particular, we first evaluate the gauge integral through residues and therefore do not rely on saddle point approximations. In addition, we do not require any specialization of the parameters to perform the computation. However, as one varies the chemical potentials some poles may enter or exit the contour of integration. Although this prohibits us from finding a fully explicit expression for the index for general parameters, this is inconsequential for the study of the Cardy limit. The residue analysis is the same as the one that relates the Higgs branch localization formulas of N = 1 gauge theories with fundamental matter to the gauge integral expression of the index [18,38,39]. The term "Higgs branch localization" in the case of the N = 4 is inappropriate, since the N = 4 theory does not have a Higgs branch. We will therefore refrain from using this terminology, even though the analysis and expressions will be completely analogous to those obtained in gauge theories using Higgs branch localization.
The computation of the index is performed in Section 2, first for SU (2) gauge group in Section 2.1 and then for SU (N ) in Section 2.2. Only after this step will we take the Cardy limit of the resulting expression in Section 3, using a modular property of the index. We conclude and discuss our results in Section 4. In Appendix A we collect definitions and some useful properties of special functions that appear in the index. In Appendix B, we compute the unrefined limit of the index, which requires special care, and its Cardy limit. Finally, we collect the expression for the anomaly polynomial in Appendix C and briefly discuss how it is related to the standard expression of the entropy function used in the main text.

Computation of the N = 4 index
In this section, we will perform the computation of the superconformal index for the N = 4 SU (N ) Yang-Mills theory. We will warm-up with the SU (2) case and then proceed to general N . Before turning to the computation, we will describe the index in some detail and in the process set up notation. See [40] for a recent review.
The superconformal index can be expressed as follows: Our parametrization is equivalent to the one used in [8] upon identifying f i = y i . Furthermore, H Q is the 1/16 th BPS Hilbert space corresponding to those states on which: vanishes. The charges J 1,2 are the rotation generators along the Euler angles of S 3 and the r i correspond to the Cartan generators of SU (4). Thinking of the N = 4 theory in an N = 1 language, the R-symmetry charges of the three chiral multiplets are given by: The SO(6) Cartan generators R i used in [8] are related to the r i via: Each R i equals 2 for a single chiral multiplet and vanishes on the other two. Finally, the N = 1 superconformal U (1) R-charge r corresponds to: Since the index is independent of continuous couplings, one can compute it at weak coupling. In this case, the trace can be explicitly performed and the resulting expression is given by: The integral over the gauge fugacities x i = e 2πiu i (x ij = x i x −1 j ) ensures the projection onto gauge invariant states. Moreover, for notational convenience we have written the elliptic Γ function Γ(u; τ, σ) as a function of fugacities: where the second equation indicates how the elliptic Γ function can be thought of as a generating function for multiletter indices from a single letter index. More precisely, writing the integrand of (2.7) in terms of plethystic exponentials, the argument of the combined exponentials is 1 − f with f the 1/16 th single letter index of the N = 4 theory. The additional 1 originates from the Vandermonde determinant obtained by replacing the matrix integral with an integral over eigenvalues [1,41].
In addition, we defined f 3 = pqf −1 1 f −1 2 and x N = (x 1 · · · x N −1 ) −1 , the latter corresponding to the SU (N ) constraint. Furthermore, κ N consists of the Cartan factors of both the chiral multiplets and the vector multiplet, and is given by: where the q-Pochhammer symbol is defined in (A.2).
Finally, for convergence of the product formula for the elliptic Γ function, one should require |p|, |q| < 1. For the summation formula, one needs in addition: |pq| < f a < 1. (2.10) Notice that if the second requirement holds for f 1,2 , it automatically holds for: The domain of convergence of the Γ functions can be extended outside the unit disk. This is discussed in more detail in Appendix A.

SU (2) index
For SU (2) gauge group, the gauge integral consists of a single contour integral: . The computation of this integral by residues has already been done in [39]. More precisely, in that paper the index is computed for an N = 1 SU (2) vector multiplet coupled to N f fundamental chiral and Nf anti-chiral multiplets. Even though we consider adjoint chiral multiplets, for SU (2) gauge group the only distinction at the level of the index is the power x 2 instead of x appearing in the argument of the Γ functions. Keeping this minor distinction in mind, the following computation may be viewed as a special case of the computation in [39] for N f = 3 chiral multiplets and Nf = 0 anti-chiral multiplets.
We will compute the integral by picking up residues from the poles of the integrand inside the unit circle. From its definition (2.8), it is not difficult to see that the Γ(x) function has simple poles at x = p −k q −l for k, l ∈ Z ≥0 . Therefore, each Γ(x 2 f a ) in the numerator has poles at: whereas each Γ(x −2 f a ) factor has poles at: (2.14) Given our restriction (2.10), only the latter poles lie inside the unit circle and therefore only these will contribute to the residue sum.
The Γ functions in the denominator do not contribute poles. One way to see this is to note that: where the θ q function is defined in (A.1). Now, the statement that the Γ functions in the denominator do not contribute poles follows directly from the fact that the right hand side of (2.15) only has zeros.
For future use, we will now derive the following residue corresponding to a basic pole (k = l = 0) of the elliptic Γ function: where the contour γ(±f To compute the residues at the more general poles (2.14), we need the following two properties of the Γ function (see Appendix A for more details): . (2.17) Now we are ready to evaluate the contour integral. We deform the contour such that it splits into a sum of three "towers" of contours, encircling the poles at x 2 = f a p k q l for a = 1, 2, 3 and k, l ≥ 0. The sum of the residues of the integrand at the two basic poles x 2 = f a reads: Here, we used the shift properties of the Γ function listed above to extract the same prefactor as on the basic pole. Noting that f 1 f 2 f 3 p −1 q −1 = 1 and that the 1 (p;p)∞(q;q)∞ cancels that same part in κ 2 , upon summing all residues one arrives at the final result: . (2.20) Let us note here that the sum over k, l factorizes into a part that only contains θ q functions and a part that only contains θ p functions. In the Higgs branch localization literature, these functions are called vortex partition functions. As already mentioned in the introduction, the N = 4 theory strictly speaking does not have a Higgs branch. Despite this, we will still use the terminology of vortex partition functions to indicate these products of θ functions because of their similarity to the vortex partition functions of gauge theories with fundamental matter. The final result can then be expressed as: where: .
(2.22) 4 We use the notation x 2 = f a p k q l to denote the collection of the two poles x = ±(f a p k q l ) 1 2 . This is also the reason why in the expression for the residue, the factor of 1 2 has disappeared with respect to (2.16).
Let us make some final comments. First of all, notice that this computation is valid for values of the chemical potentials obeying (2.10). In particular, we do not need to constrain the values of τ and σ, which is required in the Bethe Ansatz formalism [31]. Secondly, the unrefined limit f 1 = f 2 = f 3 = (pq) 1 3 of our expression (2.21) is singular. This can be traced to the fact that in this limit, the integrand develops cubic instead of simple poles. Therefore, to access the unrefined limit in our formalism one has to redo the residue computation, now taking into account the higher order poles. We defer this analysis to Appendix B.

SU (N ) index
We would now like to do a similar computation for SU (N ) gauge group. The expression for the index was given in (2.7). For convenience, we implement the SU (N ) constraint such that the index can be written as: .

(2.23)
To compute such multidimensional contour integrals, one cannot in general resort to Cauchy's theorem directly. The reason for this is that poles may not factorize in their dependence on x i , as is indeed the case for the integrand at hand. Let us therefore briefly review how to deal with such multivariate residue integrals.
Interlude on multivariate residues: Let g(x) = g 1 (x), . . . , g n (x) : C n → C n and h : C n → C be holomorphic functions. We are interested in computing the residue of the meromorphic n-form ω: A pole of ω is defined as an isolated point p ∈ C n such that g(p) = 0. The residue of ω is now computed by the integral: where γ p is an n-torus centered around p.
One can evaluate the following Jacobian determinant at a pole x = p: If J p = 0, which will turn out to be true away from the unrefined limit, 5 one can perform the coordinate transformation y i = g i (x) such that the poles factorize in y coordinates. Then, the residue can be evaluated as a product of one-dimensional residue integrals: We will apply this general formula to the computation of (2.23) by first classifying all the poles of the integrand. Subsequently, we deform the contour such that it splits into a sum of (N − 1)-tori, each of which encircles a pole of the integrand p = (x 1 , . . . , x N −1 ) for which all x i lie inside |x i | = 1. Other poles will not contribute to the resulting residue sum.
Back to the index: For the same reason as in the SU (2) case, all poles of the integrand originate from the Γ functions in the numerator. The total number of Γ functions in the numerator is equal to 3(N 2 − N ). A (simple) pole of the integrand is realized at those points where N − 1 of these Γ functions have a pole. Therefore, poles of the integrand can be found by selecting N − 1 Γ functions in the numerator Γ(y i ) and subsequently solving the system of equations y i = p −k i q −l i for some k i , l i ≥ 0.
The 3(N 2 − N ) pole equations are linear equations when written in terms of the chemical potentials: for some k ij , l ij , k i , l i ≥ 0. Selecting N − 1 of these equations and solving them for u i leads to a pole of the integrand. At first sight, it may seem that for large-N this leads to a huge number of poles to analyze. However, there turns out to be a significant reduction in the number of poles that contribute non-trivially to the full residue sum as we will now argue.
First of all, for the system of N − 1 equations to be solvable, at least one of the equations has to be of the type in the second line of (2.28). In particular, the system is solvable when choosing the N −1 equations to be all of the second type. More generally, given some number of equations of the second type, the system remains solvable as long as one adds equations of the first line such that the system can be rewritten in the form of N − 1 equations of the second type. For example, say we choose the i th equation from the second line, while the rest of the equations come from the first line. Then, the system is solvable when we choose the N − 2 u i − u j equations for 1 ≤ j ≤ N − 1 and j = i, since subtracting these equations from the i th equation of the second line brings us back to the system with N − 1 equations of the second type. Now comes the first main simplification of the analysis. Suppose that at least one of the equations in the (solvable) system is of the first type. For definiteness, let this equation be: For each such equation of the first type in the system, there exists another system of N − 1 equations where in all equations the labels i and j are exchanged and such that φ a ij = φ a ji , k ij = k ji and l ij = l ji . The solutions to these two systems are identical, up to an exchange of u i and u j . Together with the fact that the residue formula (2.27) is odd under the exchange of two integration variables x i and x j , this implies that the residues corresponding to these two systems are equal but of opposite sign. Therefore, they will cancel in the sum over residues.
This leaves us with the analysis of the class of unpaired poles: the solution to the system of N − 1 equations of the second type. Let us consider first the N − 1 equations with the + sign on the right hand side, since these poles have the best chance of lying inside the unit circles. The solution of this system is as follows: Here, we indicate by p (n) with n = 1, . . . , N , the N distinct solutions for the x i .
If one wishes to replace the + sign in the i th equation with a − sign, the solution can be obtained from (2.30) by taking: Let us now check whether the pole (2.30) indeed lies within all unit circles. First of all, it is easy to see that a pole with all k i equal and all l i equal lies inside all unit circles for comparable values of the f a . This is because in this case the net power of fugacities is positive for every i. However, we will quickly see that poles with all k i equal or all l i equal have a vanishing residue. In fact, for the residue to be non-vanishing all but a few k i have to be distinct and similarly for the l i . A minimal choice of a pole with non-vanishing residue would for example be: This specific choice leads to the following absolute value: (2.33) One may convince oneself that there is a finite domain in parameter space where all x i are inside their unit circles. For example, when |p| ≈ |q|, one finds that: For comparable values of the |f a |, this point will lie inside all unit circles, while at large-N the point lies inside all unit circles for all values of the |f a |.
Furthermore, the domain for which (2.33) lies inside all unit circles can also be made parametrically large by considering the following shift: This new pole will have an additional factor of p n 1 q n 2 in the numerator, thus enlarging the domain in parameter space for which the pole lies inside all unit circles. For large enough n 1,2 , one does not need to take |p| ≈ |q| for this to be true.
A complete analysis of all poles lying inside the unit circles, including systems of equations with some + signs are replaced with − signs, is beyond the scope of this paper. In fact, for our purposes, i.e. computing the Cardy limit of the index using a modular property, it will be sufficient to know that there is at least one pole with non-vanishing residue. As argued above this is always the case, irrespective of where we are in parameter space.
We will now continue to compute the residues for the poles where all N −1 equations are taken with a + sign. Other poles originating from a set of equations including − signs can be obtained from these residues through the transformation (2.31).
The computation of the residue makes use of the formula (2.27). So, we first have to evaluate the Jacobian of the pole (2.30). For the functions g i (x), we take: Then, for each of the N poles (2.30) the Jacobian consists of: It is not difficult to check that this implies: Notice that the Jacobian is independent of n. In addition, the function h(x) in our case can also be seen to be equal on each of the N poles in (2.30). Therefore, upon summing the contributions of the N poles, the residue formula: has a trivial contribution from the Jacobian and one just has to evaluate h(x) on any of the N poles p (m) .
To finish the computation of the residue, we first note that on any of the N poles (2.30): x Now we can straightforwardly evaluate the residue: In the second line, second product, we kept the product over b complete, even though N − 1 of these factors are used in the residue. The reason for this is to keep notation simple. The price is that we have to add to the denominator a factor of Γ(p −k i q −l i ) to cancel those superfluous factors in the numerator. In addition, the θ functions in the last line originate from using the shift properties (2.17) for the Γ functions contributing to the pole in order to evaluate the residue in terms of q-Pochhammer symbols. This also results in C i , which is given by: Summing the result over k i and l i gives us the final result for the class of poles (2.30). Before stating the (form of) the final result, let us make the following remark. As alluded to above, for specific values of k i and l i , the residues may be vanishing. This is caused by the Γ functions in the denominator of the first factor of the residue. To see this explicitly, let us first note that (see (2.15)): Since the a i only takes three different values, for N > 4 there are necessarily terms where the f a i fugacities in the argument of the θ functions cancel. If in addition for these values of i and j either k i = k j or l i = l j , one of the θ functions on the right hand side has a zero. Therefore, generically one has to require distinct k i and distinct l i for the residue to be non-vanishing. In particular, this implies that the basic pole (k i , l i = 0) necessarily has a vanishing residue for N > 4.
We are now in a position to state the full residue sum for the class of poles (2.30): Before simplifying this expression, let us make some comments. First of all, the sum over (a i ) is similar to the sum over a in the SU (2) case and should be thought of as a sum over all possible N − 1 tuples (a 1 , . . . , a N −1 ) for a i = 1, 2, 3. Similarly, by (k i ), (l i ) and (0) we indicate the N − 1 tuples (k 1 , . . . , k N −1 ), (l 1 , . . . , l N −1 ) and (0, . . . , 0). Secondly, the prime on the index indicates that we have only considered the poles (2.30). Two other comments pertaining to the expression are: 1. As discussed around (2.33), not all poles considered in the sum will fall inside all unit circles for generic values of the fugacities. Poles which do not should not be included in the residue sum. This necessary omission is indicated by the primes on both sums. The explicit prescription would depend on where one evaluates the residue in fugacity space. Finding this prescription seems a complicated problem, and we will not consider it in the present paper.
2. Similarly, we have not considered any of the poles where some of the + signs are replaced with − signs in the pole equations in the second line of (2.28). Some of these poles may still lie inside all unit circles in some regime of parameter space and therefore could be included into a full expression for the index.
For our purposes, these two issues are not consequential for the following two reasons respectively: 1. To leading order in the Cardy limit of the index, we will show that there is a universal contribution (at large-N ) for all poles (in a specific parameter regime). Hence, to compute the Cardy limit one really only relies on the fact that the residue sum contains at least one non-vanishing residue.
2. Suppose there is a pole inside all unit circles after having replaced some of the + signs with − signs. Then, one can obtain its residue from the residue of the pole with all + signs upon using the transformation (2.31). Now, it is not difficult to see that the O(N 2 ) part of I N , captured by the i<j part of (2.44), is invariant under this transformation. Therefore, at large-N the residues for this more general class of poles are indistinguishable from the residues of poles with only + signs.
With the above comments in mind, we proceed to simplify the expression. We will make use again of the shift properties of the Γ functions (see (A.15)). Without further ado, we have: (2.45) The precise form of the vortex partition functions of the numerator Z V depends on the sign of k i − k j and l i − l j . For example, if both are positive or both negative for all i < j, then the vortex partition function is given by: . (2.46) Here we used, as in the SU (2) case, f 1 f 2 f 3 p −1 q −1 = 1. Moreover, note that the products in the denominator of the first line start at m = 1. This is because the m = 0 terms cancel the terms, originating from the Γ functions in the denominators of the first line of (2.45). We find it more transparent to cancel these terms against each other directly, instead of defining Z V including the relevant m = 0 terms and keeping the Γ functions in the expression for the index.
For more general signs of k i − k j and l i − l j , one has to make use also of the second and third line of (A.15). The net effect is a shuffling of θ functions between denominator and numerator. Since this is not consequential for our purposes, we refrain from providing the precise formulas. The precise form only depends on (k i ) or (l i ), and is summarized in that label on the vortex partition function.
Summing over all residues, we find a final expression for the index: 7 (2.47) To repeat, there are two main provisos that should be kept in mind when reading this expression for the index. Firstly, the precise specification of the primed summation 6 It is also possible to compute this partition function from the point of view of the vortex worldsheet theory (see [42,43] for examples of such a computation in N = 1, 2 gauge theories). In particular, it should match the elliptic genus of a specific (4, 4) GLSM appearing for example in Section 5.1 of [44]. Recently, the computation of the elliptic genus for a special example of the relevant GLSM appeared in [45]. 7 In this expression, Γ(1) N −1 is only included to cancel the Γ(1) N −1 coming from the last product in the first line. These latter factors should not be included in the first place, since these represent precisely the factors that define the pole at which we evaluate the residue. The reason to include them is purely for notational convenience. This is related to the remarks in the paragraph below (2.41).
domain depends on where one is in parameter space. This means that the final expression is not fully explicit. Furthermore, we did not include any other classes of poles, corresponding to exchanges of + sign equations with − sign equations. The first issue is problematic if one wants to study the index exactly. However, we will see that for purposes of the Cardy limit it is not consequential at leading order. The second issue can be easily resolved by inserting an additional sum over all transformations (2.31). However, also for these transformed residues, one would need to establish a summation domain, which will differ from the untransformed residue sum. As mentioned in the main part of the section, the second issue disappears in the Cardy limit if combined with a large-N limit, since the leading part of the residue at large-N is universal for any combination of + and − sign equations. Finally, as in the case of SU (2), this expression is singular in the unrefined limit. This is due to the fact that the integrand develops higher order poles, which require a separate analysis. We perform this analysis in Appendix B.

Cardy limit of the index
In this section, we will study our final expression for the index (2.47) in the Cardy limit. As commented above, even though the expression is not fully explicit, we will see that it suffices for our purposes.

Cardy limit of θ and Γ functions
In order to study the Cardy limit of the index where τ, σ → 0 +i keeping σ/τ ∈ H \ R fixed, we will make use of modular properties of θ and Γ functions. We will briefly review these properties here. For a collection of relevant formulas, see also Appendix A.
θ function: The θ function satisfies the following modular property: where B(z, τ ) is defined in (A.8). Using the summation formula for the θ function (A.3), we have in the limit that τ → 0: if This domain is illustrated in Figure 1. In this domain, the modular property implies: which we will call the Cardy limit of the θ function.
The domain where (3.4) holds can be extended to z ∈ C \ Z + γ, where γ is the line running through 0 and τ . In the figure, such domains correspond to arbitrary horizontal integer shifts of the strip. This extension is possible due to periodicity under z → z + 1 of the left hand side of (3.1), which is reflected on the right side through: The second line follows from the quasi-ellipticity of the θ function (A.4). One easily sees that the shift properties cancel in the product, thus reproducing the periodicity of the left hand side. The extension of (3.4) to z ∈ C \ Z + γ can now be performed through repeated use of (3.5). The result is that for any z ∈ C \ Z + γ, we can write: where the bracket is defined as: [z] τ ≡ z + n, n ∈ Z such that Im In words, the bracket implements a horizontal shift on z such that its image lies in the fundamental domain indicated in Figure 1. It is easy to verify that the brackets satisfy the following relations: Γ function: The elliptic Γ function also satisfies a modular property [33], as recently discussed as well in [20]. For Im(τ ), Im(σ), Im σ τ > 0, one has: 8 Here, Q(z, τ, σ) is defined in Appendix A.
Following [20], consider the limit τ, σ → 0 with τ σ / ∈ R of the Γ factor in the numerator on the right hand side of the modular property. Using the summation formula (A.12), we find: (3.10) Similarly to the θ function, this factor becomes equal to 1 when: The story is exactly the same for the Γ function in the denominator if: The intersection of these two domains is shown in Figure 2. For later convenience, we will refer to the green shaded diamond shaped domain as D 0 . We will refer to integer shifts of this domain by D n , where n ∈ Z is defined in (3.7) and indicates a horizontal translation of D 0 by −n. Thus, when the z variable sits inside the diamond shaped green domain, the modular property reduces in the Cardy limit to: Similarly to the case of the θ function, since Γ(z; τ, σ) is periodic under z → z + 1 this limit can be extended to z ∈ D n for any n ∈ Z. The periodicity of the left hand side of (3.9) is reflected on the right hand side through the following identities: (3.14) where we used the basic shift property (A.14) of the elliptic Γ function, the shift property (A.4) of the θ function, the extension (A.6) of the θ function and finally the modular property (A.9) of the θ function. In particular, the quadratic polynomial Φ(z, τ, σ) in the chemical z appears in the latter modular property and is defined in (A.10).
Repeated use of these identities allows one to extend the limit (3.13) to z ∈ D n for any n ∈ Z: lim where the bracket is defined similarly as in (3 .7): However, in this case only the first and the third relation of (3.8) are satisfied: The second relation is not satisfied since enough translations by τ (or σ) will eventually bring a point inside the diamond into one of the red regions. This will not be an issue for the remainder, since in the Cardy limit τ, σ → 0 and we can simply ignore their appearance in brackets to leading order. Notice also that this bracket is only defined for z ∈ D n . This constrains the possible values of z for which the Cardy limit results in (3.15), since the D n do not cover the full complex z-plane unless arg(τ ) = arg(σ).
Let discuss the diamond domain D 0 in a little more detail. In the following, we take z = φ a . Depending on the sign of Im(φ a ), the domain can be written as: (3.18) Notice that our choice Im σ τ implies that Re(τ ) Im(τ ) > Re(σ) Im(σ) (see Figure 2). It is then clear from these equations that the domain is maximal at Im(φ a ) = 0 and shrinks linearly for both signs of Im(φ a ), as is manifest from Figure 2. In particular, the interval shrinks to zero size when: Therefore, if |Im(φ a )| is close to this value, φ a will generically lie inside a red domain of Figure 2. In this case, the divergence of the τ, σ → 0 limit cannot be isolated inside the Q function as in (3.15). Instead, one has to keep (one of) the elliptic Γ functions on the right hand side of (3.9). 9 To avoid restrictions on Re(φ a ), we will consider the following regime in parameter space: . (3.20) This limit zooms into the part of the domain where the difference between the τ and σ strip is very small. Effectively, in this regime we may treat the diamond as a strip and for generic values of φ a in this regime: φ a ∈ D n for some n. In this limit we can rewrite the domains in (3.18) as follows: where we defined scaled normal componentsφ + andφ − to the upper and lower right boundary of the diamond respectively: (3.22) Notice that for arg(τ ) and arg(σ) close enough, (3.20) is not a strong constraint. However, we do want the difference between arg(τ ) and arg(σ) to be finite in order to use the modular property.

Cardy limit
We are now ready to compute the Cardy limit of our expression for the index (2.47).
In the Cardy limit, the θ functions diverge as exp( 1 |τ | ) or exp( 1 |σ| ), as can be seen from (3.6). On the other hand, the elliptic Γ functions diverge as exp( 1 |τ σ| ), as shown in (3.15). Therefore, to leading order we may ignore the θ functions and the difficulties associated to the sum over (k i ) and (l i ) commented upon at the end of Section 2.2. Instead, we only have to consider the part of the residues that is made up from Γ functions: Here, we temporarily reinstated the Γ functions in the denominators, which cancel against θ functions in the vortex partition functions (see the remark below (2.46)). We warn the reader that this may cause confusion for two reasons. Firstly, it falsely suggests the residues are all vanishing. This is because the f a i can only take on three values, implying that for N > 4: necessarily has a zero. However, this is not a true zero since these factors are cancelled by the m = 0 terms that we left out in the definition (2.46). In addition, the rewriting is not necessary for taking the Cardy limit, again because these terms are cancelled by the m = 0 θ functions. In spite of this, we keep these Γ functions for the moment, because they allow a nice derivation of the anomaly polynomial of the theory as we will now show.
We use the modular property (3.9) to replace all elliptic Γ functions in (3.23). In particular, every Γ function contributes a certain Q function to the overall prefactor. Collecting all those Q functions, we find: where the summation runs over the set for 1 ≤ i < j ≤ N − 1. Interestingly, the full summand of the φ summation does not depend on φ if φ 3 = τ + σ − φ 1 − φ 2 − 1. At this stage, it is not clear why we should choose φ 3 as such. Indeed, at the level of the index any integer k could have been added: φ 3 = τ + σ − φ 1 − φ 2 + k. We will derive k = −1 later when considering the Cardy limit in a specific region of parameter space. For now, let us just take k = −1 and in addition note that: Since the second line of (3.25) does not depend on φ we may set φ = 0 to find: This object is formally identical to the supersymmetric Casimir energy [22,23], although the latter is defined for k = 0 in φ 3 . This apparently small distinction played a crucial role in the derivation of the AdS 5 black hole entropy of [6] (see also [9]). Furthermore, (3.27) is very closely related to the anomaly polynomial of the N = 4 SU (N ) theory (see Appendix C), as first observed in [46]. The fact that anomaly polynomials can be derived through the use of the modular property of elliptic Γ functions was already known from several works, including [15,17,18] and discussed recently in detail in [20].
Notice that naively (3.25) seems to depend on the choice of pole, i.e. a choice of a i . However, the above shows that it is independent of this choice. The interpretation of this fact, as for example appearing in [20], is that the residue sum can be thought of physically as a sum over Higgs branch vacua and the anomaly polynomial Q tot should not depend on a specific vacuum. 10 Finally, this object is also the so-called entropy function that upon a Legendre transformation leads to the correct black hole entropy [9] (see also [6][7][8]).
However, as our analysis of the Cardy limit of the Γ function has indicated, this is not yet the end of the story. In order to evaluate the limit, we have to evaluate Q tot on bracketed potentials (see (3.15)). We now only take the Q polynomials of the Γ functions appearing in (2.47). In other words, we ignore the −Q(φ) − Q(−φ) part of (3.25), which originates from the Γ functions in the denominator of (3.23). Consistent with the identity (2.15), this part is subleading in the Cardy limit: where we made use of the bracket relations (3.8).
Thus, the total polynomial to consider is now given by: Before turning to an analysis of this object as a function of the specific residue, we will make some comments: 1. At large-N , one may just consider the summation over φ a i − φ a j since only this part contributes O(N 2 ) terms.
2. The brackets will reintroduce a dependence on the summation variable φ, in contrast with the analysis of (3.25). In principle, this implies that different residues contribute differently in the Cardy limit. Notice that this presents an opposing point of view to the analysis in [20], where it is asserted that also in the Cardy limit, the residues/vacua contribute universally.
3. We will show below that the latter point can be partially resolved by taking N large and restricting the flavor fugacities appropriately.

Equal a i
Let us start with the analysis of a residue for which a i = a j for all i, j. Focusing on the part that scales with N 2 , 11 we find: To evaluate the bracket appearing in the last Q function, we first note that we may ignore τ and σ to leading order in the Cardy limit. Now, there are two possibilities depending on whether [φ 1 ] + [φ 2 ] ∈ D 0,1 (see Figure 2): we cannot proceed, as we explained at the end of Section 3.1. The two choices lead respectively to: We can make the expression look more symmetric by using [φ 3 ] = [−φ 1 − φ 2 ]. In this case, the expression becomes: This is the expected answer for the entropy function (3.27), 12 where the two possibilities correspond to the twin saddles discussed for the first time in [6][7][8]. In particular, a Legendre transform of both expressions gives rise the correct black hole entropy [9] (see also [6][7][8]). Even though this is encouraging, we still have to analyse more general residues before we can make definite statements about the Cardy limit of the full index.

Unequal a i
To understand if there is a universal contribution from all residues to the index in the Cardy limit, or if some residues are subleading with respect to the residue studied above, we will now analyse the most general residue. Apart from the terms in the sum over φ a i − φ a j where a i = a j , there are three other possible terms contributing at leading order in N , corresponding to the pairs (a i , a j ) = (1, 2), (1, 3), (2,3). The sums over b = 1, 2, 3 in (3.29) work out for each pair respectively as: To proceed, we first have to understand how to evaluate the more general brackets appearing in the arguments of the Q functions. For a bracket [aφ 1 + bφ 2 ] with a, b ≥ 0, the following cases should be considered: (3. 36) In addition, the brackets [2φ 1,2 − φ 2,1 ] can be evaluated as follows: If any of the brackets do not fall within any diamond, we cannot proceed. As explained in Section 3.1, we can avoid this issue by working in the limit (3.20). In the following, we will always work in this limit.
We will now compute the value of the added constant integer for all brackets as a function of −1 < [φ ± 1,2 ] < 0, whereφ ± 1,2 were defined in (3.22). As functions ofφ ± 1,2 each separate bracket [aφ 1 + bφ 2 ] divides the square −1 < [φ ± 1,2 ] < 0 up into |a| + |b| parallel Notice that this corresponds to the usual unrefined point φ = 1 3 (τ + σ ± 1) (up to a translation to D 0 ) in the Cardy limit. Also, note that these points obey the requirement (3.20). At these points a i = a j is irrelevant, and for the same reasons as in the case when a i = a j for all i, j we obtain the entropy functions (3.34) (see footnote 12), now evaluated at the two unrefined points in (3.38), respectively. We should note here that naively we are not allowed to evaluate (3.29) at the unrefined points. This is due to the fact that at the unrefined points our expression for the index is singular because the integrand develops higher poles in this limit. However, at large-N , i.e. when restricting to the φ ∈ {φ a i − φ a j } terms in (3.29), we can safely take the unrefined limit of Q tot . This is because this part of Q tot does not change when taking the higher order poles into account, as we show explicitly in Appendix B. We will return to finite N at the end of this section, where we will see how the O(N ) contributions to Q tot prohibit its evaluation on the unrefined point and we instead have to resort to the expression for Q tot in Appendix B. Concluding, at least in the unrefined limits and at large-N , our expression for the index in the Cardy limit reproduces the expected entropy functions.
Let us now try to move away from the unrefined point. Before considering other domains, let us first just consider the domains that contain the unrefined points. These domains are given by the blue and yellow domain in Figure 4.
In the blue domain, the (non-trivial) brackets evaluate as follows: Evaluating the Q ab polynomials on these brackets, we find: Again, we can write this more symmetrically in terms of ] − 1 as follows: This result shows that as we move away from the unrefined point, as long as we stay in a small enough region around the unrefined point [φ 1 ] = [φ 2 ] = − 1 3 , we still obtain a universal residue up to small corrections of order (φ a − φ b ) 2 . That is, for this region in parameter space and at large-N we have: (3.42) where the pair (a, b) takes the values (1, 2), (1, 3) and (2,3). We stress that we did not have to impose φ 3 = −φ 1 − φ 2 − 1. Instead, the integer −1 emerges from a careful limit of the modular property and the associated bracketed potentials. This is very similar to the emergence of this integer in the analysis of [8].
We can repeat the above analysis for the yellow domain in Figure 4. Now, the unrefined point [φ] = − 2 3 is enclosed in the domain. For this domain, the non-trivial brackets take the following form: (3.43) Plugging these brackets into (3.29), we find: where we remind the reader that: [φ a ] = [φ a ] + 1.
Similarly to the blue region, all remainder terms are small close to the unrefined point . Therefore, close enough to the unrefined point we are also able in this case to conclude that each residue contributes a universal Q function in the Cardy limit and at large-N up to small corrections: This function coincides with the entropy function for the twin saddle (again, see footnote 12).

Other domains
Apart from the blue and yellow domains of Figure 4, there are many more domains in Figure 3. In each of these domains, the set of brackets (3.39) takes on a different value. Similarly as for the blue and yellow regions, we can ask for each of these regions whether they contain a point where all the Q ab are equal to each other. This would be indicative of the existence of a universal residue in such a region. We find that there are two possible scenarios. In the first scenario, which occurs for most domains, the point at which the Q ab associated to a certain region are equal falls outside that region. Therefore, it seems that we cannot associate a universal residue to these regions.
Another scenario occurs for those domains whose boundary touches or overlaps with the lines [ In these cases, all Q ab are equal and vanishing along the intersection. Interestingly, the entropy function for φ 1 + φ 2 ∈ D 0 vanishes along φ 1 = 0, φ 2 = 0 and φ 1 + φ 2 + 1 = 0 whereas the entropy function for φ 1 + φ 2 ∈ D 1 vanishes φ 1 + 1 = 0, φ 2 + 1 = 0 and φ 1 + φ 2 + 1 = 0. Therefore, for these special regions we find that our result is again consistent with the known entropy functions (see footnote 12). However, we find that moving away from the vanishing locus generates remainder terms for at least one of the Q ab that are of the same order as the entropy function. Therefore, we cannot keep the remainder terms for all three Q ab small while moving away from the vanishing locus, in contrast to the case of the blue and yellow region in Figure 4.
Concluding, we note that our expression for the index singles out the blue and yellow regions: only in these regions, close enough to the unrefined points, do we find a universal residue. Moreover, this universal residue is consistent with the results in the literature. An important difference with the Bethe Ansatz analysis is that in their case in the relevant parts of parameter space only a single residue dominates [8].
For the other domains we remain inconclusive because there does not seem to exist a universal residue. This could mean that, instead of being able to extract such a universal piece, one would have to sum the various residues to find the Cardy limit. This seems a very complicated task. Or one may be able to argue that some residue provides the dominant contribution to the residue sum in the respective domain, similar to [8]. If a dominant residue would correspond to a residue with all a i = a j , we would get the expected entropy functions for either the right upper triangular region or lower left triangular region in Figure 3, consistent with the literature. However, we have so far not been able to find a convincing argument for this scenario. Finally, let us note that whether a given residue will contribute to the residue sum depends in particular on the values of |f a |, as we have seen in Section 2.2. To fully understand the question of the existence of a universal residue in different regions of parameter space, one should first know which residues contribute in the first place. We do not expect this to fully resolve the issue, but further analysis of this point is required.

A brief look at finite N
At finite N , we should take into account all terms in (3.29). Explicitly, we have: (3.48) The first line gives the expected entropy function in the blue domain of Figure 4, while the second line will gives us the entropy function up to remainder terms proportional to (φ a i − φ a j ) 2 as we explained above. However, the last line gives rise to new brackets and remainder terms. There are again three possible terms: At large-N we have been able to ignore the fact that our computation of the index really requires us to stay away from the unrefined point, since the O(N 2 ) part of the Q tot remains unchanged in the proper calculation (see Appendix B). However, the Q functions in (3.49) coming in at O(N ) do not appear in the unrefined limit and should therefore not be evaluated at the unrefined point. This can also be seen from the terms with the brackets [−2φ 1,2 − φ 2,1 ] and [φ 1 − φ 2 ], which are not defined at the unrefined points. In particular, these brackets divide the blue and yellow domain of Figure 4 into six new regions. The unrefined point lies precisely at the intersection of the boundaries of regions. See Figure 5. For this reason, we need to use the results from Appendix B to study the Cardy limit of the index at finite N at the unrefined points. The finite N Q-polynomials are given in (B.28), which we repeat here for convenience: (3.50) These functions correspond to the total Q polynomials at the unrefined point in the blue and yellow domain respectively. The fact that we do not see the expected N 2 − 1 is explained at the end of Appendix B.2.
At finite N , we cannot move away as nicely from the unrefined point as we did at large-N . The reason for this is that the unrefined point lies at the intersection of the domains in the refined computation. Because in each of the six domains, the set of brackets takes on a different value, moving away from the unrefined point in this case depends on the direction one is going in. In particular, the resulting remainder terms associated to any of the domains are not small because of the discontinuous nature of the brackets.

Discussion
In this paper, we have computed the superconformal index of the N = 4 theory through residues. This method is the same as the one used to derive the Higgs branch localization formula from the index, which as far as we are aware was only applied in the literature to gauge theories with fundamental matter (see e.g. [18,38,39] and references therein). 13 We have seen that in the refined computation, the adjoint matter results in non-degenerate non-factorized poles. Because the poles are non-degenerate, it is still fairly easy to compute the residue. 14 The real complication shows up in deciding when a pole falls inside all unit circles, which depends on the precise values of the chemical potentials. This prohibits us from finding a fully explicit expression for the index at general values of the parameters. Luckily, for purposes of taking the Cardy limit we have been able to argue that this complication is irrelevant. Studying the Cardy limit of our expression, we first of all find that all residues contribute at leading order. In addition, close to the unflavoured or unrefined points in chemical potential space, all residues become equal to the entropy function of the 1/16 th BPS asymptotically AdS 5 black hole and its "twin saddle" respectively (see the remark in footnote 12).
The Bethe Ansatz method, as developed in [29][30][31], also computes the gauge integral through residues and was applied to compute the large-N limit and the study of AdS 5 black holes in [8,32]. However, our method is technically distinct from the Bethe Ansatz method and the final formulas have important differences. Let us discuss this in some detail. Firstly, the following rewriting of the gauge integrand is employed in [31]: where Z is short-hand notation for the integrand (see e.g. (2.7)). The chemical potentials coupling to the angular momenta are specialized: τ = aω and σ = bω with a, b ∈ Z. Furthermore, C is a new contour andQ i = 1 represent the so-called Bethe Ansatz equations. Inside this new contour, one only has to take into account poles coming from 13 In the context of the Hilbert series, similar residue computations were performed for SQCD in [48] and for adjoint SQCD in [49].
14 In Appendix B, we compute the index in the unrefined limit. In this case, the poles are degenerate and one needs more sophisticated techniques to compute the residue, as we discuss in detail.
the denominator, i.e. solutions to the Bethe Ansatz equations. This is a very different set of poles than the set we consider, which originate from Z. In particular, the basic solutions to the Bethe Ansatz equations do not depend on the R-symmetry chemical potentials φ a , whereas all our poles do. 15 Moreover, in our case all poles are explicitly known even though whether the associated residue contributes depends on the precise values of the chemical potentials. In contrast, in the Bethe Ansatz approach not all poles are known due to difficulties in solving the Bethe Ansatz equation. However, in their case the residue associated to a pole will always contribute. Another difference that stands out is that in the Bethe Ansatz method at large-N , there is generically a single dominant residue which captures the correct entropy functions. For us instead, it seems that generically all residues contribute at leading order in the Cardy limit, as we have discussed at length in Section 3.2. It would be very interesting to compare the two methods. Since we only have a fully explicit general expression for SU (2) gauge group, this may be a good place to start. Otherwise, for SU (N ) one would have to consider a specific regime in parameter space where our expression can be made fully explicit. We leave this as future work.
This work was originally motivated to understand what role certain modular properties of four-dimensional supersymmetric partition functions [20] play in the evaluation of the N = 4 SU (N ) superconformal index and the associated gravitational interpretation. Indeed, a very interesting question is whether in AdS 5 /CFT 4 there is an analogous story to the Farey tail expansion of the elliptic genus and the SL(2, Z) family of BTZ black holes familiar from AdS 3 /CFT 2 [12][13][14]. Such a story may also connect to the (m, n) saddles of the matrix integral found in [11] using an elliptic extension of the gauge integrand. We believe our expression for the index is suited for such a study for the following basic reason: we have performed the gauge integral before taking the Cardy limit, which firstly allows us to justify the use of the modular property a priori and secondly makes the relation between modularity and the Cardy limit more transparent. Instead, taking the Cardy limit at the level of the gauge integrand as in e.g. [7,10,24], these two consequences are less transparent. Given our expression, we can study more general modular properties obeyed by the elliptic Γ function [20]. Even though the modular property used in this paper is suited for the study of the Cardy limit τ, σ → 0, such other modular properties can be useful to study more general "Cardy" limit. In addition, they can help to explore the phase structure suggested by [11] within our formalism. This will be the subject of a future publication [21].
In a previous work [50] we studied a CFT 2 subsector of the N = 4 theory from the perspective of the superconformal index. We found that a two-dimensional Cardy formula arises from the superconformal index, which is indicative of ordinary SL(2, Z) modularity. It would be interesting to understand how this relates to the modularity mentioned in the previous paragraph. See also [16] for earlier and possibly related work. Other future directions include considering N = 1 superconformal field theories, which have been studied in the Cardy limit in [25,51], by Bethe Ansatz methods in [32,52,53], and finally by large-N saddle point approximation in [34]. Finally, it would be interesting to understand if our expression allows one to compute subleading corrections in the Cardy limit.
θ function: The θ q function, also known as the q-theta function, can be defined as an infinite product for Im(τ ) > 0: where q = e 2πiτ , x = e 2πiz and the q-Pochhammer symbol is defined as: Alternatively, there is the summation formula defined for 0 < Im(z) < Im(τ ): The θ function is quasi-elliptic under the translation z → z + mτ + n, m, n ∈ Z: Furthermore, it satisfies a reflection property: and can be extended to Im(τ ) < 0 through: .
Finally, the θ function satisfies a modular property: where: Another version of the modular property, which we also require, is given by: where: (A.10)

B Unrefined limit of the index
In this appendix, we consider the unrefined limit of the index: This limit has to be treated separately because any set of three simple poles y = f a p k q l associated to three Γ factors of the form 3 b=1 Γ(yf b ) collide into a single cubic pole. Even though the final expression for the index is different in some respects, we will find that, at large-N , the leading order expression in the Cardy limit remains unchanged. This justifies the naive unrefined limit of the Cardy limit of the refined index computed in the main text. As in the above, we treat SU (2) and SU (N ) separately.

B.1 SU (2) index
The expression for the index of the N = 4 SU (2) theory in the unrefined limit reads: Inside the unit circle, there are now cubic poles at: for k, l ≥ 0.
Let us first compute the residue of Γ(x) 3 at its basic cubic pole, x = 1. We first write: (B.4) To obtain the residue, Cauchy's integral formula tells us that we have to compute the second derivative of the function in brackets and evaluate it at x = 1. Let us call this function Γ(x). We first rewrite Γ(x) using the plethystic exponential (cf. (2.8)): The residue can be written in terms of γ as follows: Here, we have: This leads us to our final expression for the residue at the cubic pole: This residue, together with the use of shift properties of the elliptic Γ functions as in Section 2.1, allows us to find for the unrefined SU (2) index: where: (B.10)

B.2 SU (N ) index
We will now attempt a similar computation for the SU (N ) theory. The unrefined limit of the index is given by: .

(B.11)
In this case, poles arise from the intersection of any solvable set of N −1 of the following N 2 − N equations: where we defined f = e 2πiφ .
For the same reasons as in Section 2.2, we only consider poles originating from the pole equations on the second line of (B.12) with all + signs. However, in contrast with the analysis in Section 2.2, it is not possible to use the formula (2.27) for a multivariate residue. We will now briefly discuss the reason for this and provide an alternative formula that can be used.
Firstly, given the N − 1 equations of the second type with all + signs, a natural choice for the functions g i (x) is: 16 The issue with higher order poles is that in this case: This prohibits the use of the formula (2.27) to compute the residue, which is only defined for non-degenerate residues with J p = 0.
Luckily, there exist more sophisticated techniques to compute the residue of such degenerate multivariate residues. The main formula for the degenerate case is reviewed around Theorem 1 in [54], in which also additional references may be found. The basic idea is still similar to the non-degenerate case: one wants to find a transformation that factorizes the multivariate residue integral into a product of univariate ones. The main formula, Theorem 1 in [54], for the evaluation of degenerate multivariate residues is given by: This formula is the analogue of (2.27) for a degenerate pole, i.e. an isolated zero at x = p of g(x) = (g 1 (x), . . . , g n (x)) with J p = 0. In this formula, the g i (x i ) are functions that only depend on x i , and can be obtained from the g i (x) via: where A(x) = (a ij (x)) is a matrix of holomorphic polynomials. The polynomial A(x) can be determined algorithmically through a so-called Gröbner basis computation. For details, we again refer to [54] and references therein.
We have done the Gröbner basis calculation through Mathematica at low values of N . These computations suggest that the g i for general N are given by: where we have chosen to only to the computation for k i , l i = 0. We do not have to understand the more general case, since any p, q shifted pole can be first brought back to this basic pole with the help of shift properties of the elliptic Γ functions. The formula for g i (x i ) shows that the perhaps from (B.13) naively expected cubic poles are actually poles of order 2N − 1.
We now continue to compute the residue at the pole (cf. (2.30)): Firstly, we bring the residue back to x N i = f using the shift properties of the elliptic Γ function: where: Now, we use the formula (B.15) for the last line: We will not find an explicit expression for the residue. This is because of a technical reason: we have been able to find the matrix A(x) for low values of N , but it is a somewhat complicated polynomial whose generalization to arbitrary N is not clear to us. Luckily, for purposes of the Cardy limit, as we will discuss in more detail in the next section, the form of A(x) is unimportant. We will also see in the next section that the high order of the poles and the associated derivatives will not be of relevance to the Cardy limit either. Therefore, we hide all these details in a function R(f, q, p). Then, analogously to the refined analysis of Section 2.2, the expression for the index in the unrefined limit takes the form: where the Pochhammer symbols originate from the residue similarly to the SU (2) case discussed in the previous section, and details of the A(x) matrix and derivatives of it and Γ are hidden in R. As in the refined case, the precise form of the vortex partition functions of the numerator Z V depends on the sign of k i − k j and l i − l j . For both positive or both negative for all i < j, the vortex partition function is given by:

(B.23)
This is just the the specialization φ 1 = φ 2 = φ 3 of the refined vortex partition functions (2.46). All further comments made there apply here as well, so we will not repeat them.

B.3 Cardy limit of the unrefined index
As is clear from the expression (B.22), the universal part of the residue is simplified significantly in the unrefined limit of the index. We already know that to leading order in the Cardy limit, the vortex partition functions do not contribute at leading order (see Section 3). We will now argue that also the q-Pochhammer symbols and R do not contribute at leading order.
The function R(f, q, p) consists of a sum of products of derivatives of the polynomial A(x) and derivatives of γ(x) (see (B.5)). Since A(x) is a finite degree polynomial in x, the evaluation of its derivatives on the poles lead to a finite degree polynomial in the fugacities. This contributes at order O(e 2πiτ ) or O(e 2πiσ ) and can be safely ignored in the Cardy limit. For derivatives of γ(x) it is not difficult to see, for example from the expressions (B.7), they will diverge as O(e − log(τ σ) ), which is again subleading.
We conclude that also in the unrefined case to study the Cardy we only have to consider the Γ functions: Γ(f ) 3(N −1) 2 Γ(f 2 ) 3(N −1) .

(B.27)
Since f = (pq) 1 3 we have: φ = 1 3 (τ + σ − k) for some k ∈ Z. However, evaluation of the brackets, discussed in Section 3.1, reduces this choice to the independent values k = 0, 1, 2. It is easy to see that for k = 1, 2 the unrefined points lie inside the diamond D 0 in the Cardy limit (see Figure 2). Instead, for k = (0 mod 3) it will lie outside the diamond. However, also in this case it is not difficult to see that the Cardy limit of the modular property to leading order only gets a contribution from the Q polynomial. This is because both Γ functions on the right hand side of the modular property (3.9), even though they do not simplify to 1 in this case, will be convergent functions of σ τ . Therefore, their contribution will be subleading with respect to the diverging e

prefactor.
We then proceed to compute the total Q polynomial for these separate cases: (B.28) Here, we have used that the brackets evaluate to leading order as:

(B.29)
At large-N , the expressions for Q 1 and Q 2 agree with the unrefined limits of the Q polynomials computed away from the unrefined points (3.42) and (3.47). This is what we wanted to show. As a final comment, we note that the subleading pieces in N do not take the perhaps expected form of N 2 − 1. The reason for this is that the missing 3(N − 1) terms originate from the q-Pochhammer symbols, which do not contribute at leading order in the Cardy limit.
(C.5) These coefficients encode the global anomalies of the N = 4 SU (N ) Yang-Mills theory.