Black hole perturbation theory and multiple polylogarithms

We study black hole linear perturbation theory in a four-dimensional Schwarzschild (anti) de Sitter background. When dealing with a positive cosmological constant, the corresponding spectral problem is solved systematically via the Nekrasov-Shatashvili functions or, equivalently, classical Virasoro conformal blocks. However, this approach can be more complicated to implement for certain perturbations if the cosmological constant is negative. For these cases, we propose an alternative method to set up perturbation theory for both small and large black holes in an analytical manner. Our analysis reveals a new underlying recursive structure that involves multiple polylogarithms. We focus on gravitational, electromagnetic, and conformally coupled scalar perturbations subject to Dirichlet and Robin boundary conditions. The low-lying modes of the scalar sector of gravitational perturbations and its hydrodynamic limit are studied in detail.


Introduction
In this paper, we study analytical approaches to solve the differential equations describing black hole (BH) linear perturbation theory. The concrete problem consists in the study of Einstein equations with a cosmological constant, approximated around a particular BH solution to first order in perturbation theory. 1 Because of the symmetries of the Schwarzschild (Anti) de Sitter solution, the linearised equations separate and reduce to a second-order linear ODE of Fuchsian type; we refer to [1] for a review and a list of references. In the particular cases we study in this work, the relevant equation has four or five regular singularities. Fuchsian equations appear in many fields of theoretical and mathematical physics, and the relevance of the parametric analysis of their solutions and the corresponding connection coefficients goes well beyond the application to BH perturbation theory. In this paper, we employ two distinct, complementary strategies to analytically study such equations. First, by following modern developments in the context of the supersymmetric gauge theories, we tackle such problems using the Nekrasov-Shatashvili (NS) functions [2] (see Appendix A for their definition). These functions have been shown to be building blocks to compute quantum periods [2][3][4][5][6], eigenfunctions [7][8][9][10][11][12][13], Fredholm determinants [14,15], and connection coefficients [12] for Fuchsian differential equations and their irregular limits. These techniques were recently applied to studying spectral problems describing black hole perturbation theory. Initially introduced in [16] for the study of quasinormal modes (QNMs) in four-dimensional asymptotically flat black holes, this approach has been generalized to various gravitational backgrounds and extends beyond the QNMs computation [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34]. 2 Other interesting related results have been elaborated in .
In this paper, we further extend this approach to the framework of four-dimensional BHs in the de Sitter (dS) background. More precisely, we compute the relevant connection formulae following the methodology developed in [12], where exact connection formulae for the Heun equation were obtained from the classical limit of Virasoro conformal blocks, which is, in turn, related to NS functions. 3 The approach based on the NS functions can be applied in its simplest form only when the boundary conditions of the spectral problem are imposed at singular points of the differential equation. However, there are cases when this condition is not satisfied. For example, when considering gravitational or conformally coupled scalar perturbations of black holes in a fourdimensional anti-de Sitter (AdS) spacetime, see Sec. 4. Moreover, the NS functions expansion is potent but only in some regions of parameter space, and to explore other regions, it is sometimes necessary to resort to other analytic methods. This happens, for instance, in studying the socalled hydrodynamic limit. 4 In these situations, we analyze the equation using an alternative "polylog" method where we reduce the relevant problem to recurrence relations which we solve in terms of multiple polylogarithms. This paper is structured as follows.
In Sec. 2, we briefly describe the NS and the polylog methods and comment on their implementation. 1 Higher orders in perturbation theory can also be studied with the methods we develop here, but this is beyond the scope of this work. 2 See also [35][36][37][38][39][40] for another approach based on Painlevé equations. 3 For a discussion of these connection formulae from a mathematical viewpoint, see [65]. 4 Indeed, the natural expansion in this limit does not map to the instanton expansion in gauge theory. This obstacle may be overcome by using TBA techniques for the computation of the NS functions and the corresponding quantum periods, see [29][30][31]. However, we will not explore this path in this work.
In Sec. 3, we study conformally coupled scalar, electromagnetic, and vector-type gravitational perturbations of four-dimensional Schwarzschild de Sitter black holes (SdS), where both methods are applicable. We use the Heun connection formulae to obtain the quantization condition in Sec. 3.1, which gives the quasinormal frequencies as series expansions in the radius of the black hole horizon R h . In Sec. 3.2, we apply the polylog method and express the corresponding wave functions in terms of multiple polylogarithms in one variable. By gluing the relevant local solutions, we determine the frequencies of the QNMs, and the resulting series expansions in R h agree with the ones obtained in Sec. 3.1. In all computed orders, we find purely imaginary QNM frequencies in agreement with the earlier observations made by numerical computations [66][67][68].
In Sec. 4, we study the same class of perturbations of Sec. 3 in the case of four-dimensional Schwarzschild Anti-de Sitter black holes (SAdS) imposing Dirichlet boundary conditions at spatial infinity. For these perturbations, one of the boundary conditions is imposed at a regular point of the equation. Hence the method based on the NS function is more complicated to implement, and we apply the polylog method instead. As in the Schwarzschild de Sitter case, the two local solutions are described in terms of multiple polylogarithms in one variable. We check our analytic results against the numerical values of [69]. 5 Our results suggest that the leading order of the imaginary part of the QNM frequencies is −c R 2 +2 h in the small R h approximation, where is the angular quantum number and c is a real positive constant depending on all the quantum numbers. This is consistent with earlier numerical results obtained via the Breit-Wigner approach [71] and some earlier analytic studies in [72].
In Sec. 5, we study the low-lying modes of the scalar sector of gravitational perturbations of Schwarzschild Anti-de Sitter black holes in the big R h limit. In this Section, we use Robin boundary conditions at spatial infinity, which preserve the metric at the boundary and, as such, are more suited for holography. The corresponding differential equation has five regular singular points, and we use the polylog method to compute the relevant local solutions as Taylor expansions in 1/R h . To make this computation more efficient, we introduce three local regions. In the region near the horizon, the local solution is given in terms of multiple polylogarithms in several variables (all but the first argument are constants). The local solutions are described in terms of Laurent polynomials in the other two regions. The QNM frequencies are obtained by gluing the three local solutions, and the first two orders in 1/R h agree with the results from [73]. Theoretically, one can compute the QNM frequency up to any given order in 1/R h . However, due to the exponential growth of the number of polylogarithm functions that appear in each order in the perturbative expansion, we could determine the QNM frequency up to order 1/R 6 h . By taking the hydrodynamic limit, we can also reproduce the results from [74] and obtain four additional corrections in the expansion.
Appendix A introduces the notations and conventions used in the main text for the NS functions. Appendix B reports relevant identities between classical polylogarithms and multiple polylogarithms and relations for multiple zeta values. In Appendix C, we prove by induction on K ∈ N that the local solutions at order R K h in Sec. 3 and Sec. 4 are given in terms of multiple polylogarithms in one variable of weight at most K. Appendix D presents the linear basis of multiple polylogarithms in several variables that describes the local solution near the horizon in Sec. 5. We also show how nontrivial identities arise between multiple polylogarithms at a fixed level. In Appendix E, we write the Heun connection formula relevant to the Schwarzschild Antide Sitter case. We obtain the first order correction in R h of the QNM frequency with n = 0, differential equation with four regular singularities, for which the connection formulae are given in [12]. The canonical form of the Heun equation is given by To study the quasinormal modes equation, the gauge theory approach is useful when the boundary conditions are imposed at the singularities of the problem, like for the SdS case, analyzed in Sec. 3. Thanks to the power-like behaviour of the local solutions near the singular points, it is easy to identify the two local solutions selected by the two boundary conditions. These two solutions are then used to quantize the frequencies, taking into account the connection formula that relates them. In these cases the quantization condition is expressed in terms of the quantum periods of the underlying SW geometry, computed via NS functions (see formula (3.17) and Appendix A for the conventions used). If, as it happens in the SAdS case analyzed in Sec. 4, at least one of the boundary conditions is imposed in a regular point of the differential equation, the gauge theory approach is less effective: it is still possible to solve the problem with the connection formulae, but the quantization condition for the frequency will not be expressed in terms of NS functions only (see formula (E.3)).
By introducing an appropriate change of variables, we can always transform the perturbation equation with four regular singularities in the canonical Heun form and send the singularities in z = 0, 1, t, ∞. In all the cases in which the connection formulae are used, we will put us in a regime in which the complex modulus of t is small, |t| 1, and such that the relevant connection formula is among local solutions in z = t and in z = 1. The independent solutions of the canonical Heun equation for z ∼ t are (2.2) and the ones for z ∼ 1 are In terms of the connection matrices of hypergeometric functions where θ, θ = ±, the connection formula for small t from z ∼ t to z ∼ 1 is given by [12] (2.5) Within the examples analyzed in this paper, the gauge theory approach is particularly effective in computing quasinormal modes at large in the SdS case (see subsection 3.1.2).

The multi polylog approach
When the gauge theory approach proves less effective, we solve the QNM spectral problem order by order in some suitably chosen expansion parameter κ and for fixed values of some quantum numbers. 6 For instance, we can consider perturbation theory in R h or R −1 h , R h being the radius of the black hole horizon. This is like doing Hamiltonian perturbation theory, and although the numerical implementation of this algorithm is well known (see e.g. [1,[77][78][79]), if we want an analytical answer, the calculation quickly becomes cumbersome. In particular, it is necessary to find suitable stratagems for higher orders. For the situations we consider in this work, we find that higher orders can be determined systematically using the underlying structure that involves multiple polylogarithms (similar techniques are also used to compute Feynman integrals in QCD, see e.g. [80][81][82] and references therein). In this section, we sketch the general idea while we give more details in concrete examples; see Sec. 3, Sec. 4, and Sec. 5. Mathematica notebooks are also attached.
The spectral problems we are interested in are two-point boundary value problems associated with differential equations on the sphere with n regular singularities. 7 More precisely, we will focus on the cases n = 4, 5. The boundary conditions are fixed at generic points z = z 1 and z = z 2 , not necessarily coinciding with the position of the regular singularities.
For the problems at hand, we can use the following Ansatz for the eigenfunctions in each region of the patch decomposition 8 of the n-punctured sphere Sometimes, it is useful to introduce additional regions with respect to the minimal patch decomposition to optimize the efficiency of perturbation theory. Different scalings in κ of the regular singularities of the equation under the scale redefinitions z → κ c z determine the possible number of regions. We assume that in the differential equation (and therefore in the position of the 6 Here "suitably chosen" means that the 0 th order is solvable in terms of relatively simple functions (e.g. rational functions or logarithms). In addition, we would like good convergence properties for the expansions in κ in the spectral problem. We believe that this is the case, at least for the example of Sec. 3 where κ is related to the instanton counting parameter t of the underlying gauge theory (see [25,75,76] for the study of convergence properties in gauge theory). 7 The case of irregular singularities will appear elsewhere [83]. 8 We remark that the expansions in z of the local solutions, performed around a given singularity of the problem, hold in an open disk centered in that singularity. The radius of this disk is equal to the distance between that singularity and the closest one. singularities too), only integer powers of κ appear, therefore c ∈ Z. The position of the regular singularities can depend on the perturbation parameter κ. Therefore, in the perturbative expansion in κ ∼ 0, the singularities will tend to cluster differently as we change the critical parameter c. These different clustering schemes define the different regions for the perturbative expansion and determine a finer structure in the patch decomposition, which also considers the different geometric situations describing the potential terms in the differential equation. For a nontrivial example, see Sec. 5. At each order in κ, ψ (z) is determined by a second-order equation which we solve by using the method of variation of parameters. The functions ϕ and ν in (2.7) are known, 9 and the non-homogeneous part of the differential equation η K (z) is fully determined by the solutions to the previous orders f m with m ≤ K − 1. Let f 0 , g 0 be the two solutions of the homogeneous part of (2.7). 10 Then we write the generic solution to (2.7) as 11 where W 0 is the Wronskian of the two leading order solutions In each region, the integration constants c K 's can be absorbed into a normalization of the solution, and they can be fixed to zero without loss of generality. Imposing the two boundary conditions and gluing the local solutions fixes the integration constants b K and gives the quantization of the frequency of the perturbation. If either z 1,2 has a non-trivial dependence on the parameter κ, the boundary condition is applied by expanding ψ (z 1,2 ) in powers of κ. In the following sections, we will be more detailed in describing how this expansion works case by case. In principle, the relations (2.8) allow us to compute the wave function up to any given order in κ. However, to implement this algorithm in practice, there is still a non-trivial step: explicitly compute the integrals in (2.8). In all our examples, the leading order solutions are described in terms of rational or logarithmic functions, and their Wronskian is a rational function. Hence the wave function at order κ K is described in terms of multiple polylogarithms of weight K and lower. In the cases analyzed in sections 3 and 4, up to order R 4 h , one can avoid using multiple polylogarithms due to identities (B.6)-(B.9) presented in Appendix B. However, from order R 5 h on, the expansions in multiple polylogarithms cannot be avoided to our knowledge. Here we would like to mention that a more general statement about multiple polylogarithms in several variables is well-known. According to Theorem D of [84], not every multiple polylogarithm of weight ≥ 4 can be expressed as a finite combination of classical polylogarithms. Since, in Sections 3 and 4, we only deal with multiple polylogarithms in a single variable, we can push this bound to weight 5.

Schwarzschild de Sitter black hole
The metric describing the de Sitter Schwarzschild black hole in four dimensions (SdS 4 ) is where M is the mass of the black hole and Λ > 0 is the cosmological constant. In what follows, we will fix Λ = 3, and then we suppose M to be in the range 0 < M 2 < 1/27 to have three real roots for the equation rf (r) = 0, since otherwise we would have unphysical solutions. We will denote these roots by where R h ∈]0, 1/ √ 3[ is the smallest positive real root, and R ± are real and given in terms of R h by We will study a class of linear perturbations of the SdS 4 geometry with spin s ∈ {0, 1, 2}, encoded in the following radial equation (see [1] and reference therein) where the potential is For s = 0, this equation describes conformally coupled scalar perturbations; for s = 1, electromagnetic perturbations; and for s = 2, odd (Regge-Wheeler or vector-type) gravitational perturbations.
The boundary conditions we impose on the wave function are the presence of only ingoing modes at the event horizon R h and the presence of only outgoing modes at the cosmological horizon R + . These conditions can be made explicit by introducing the tortoise coordinate r * defined by In terms of r * , the behavior of Φ near R h , R + is described by plane waves, so we ask that Φ behaves as exp(−iωr * ) for r ∼ R h and as exp(iωr * ) for r ∼ R + . The latter radial equation apparently has five regular singular points located at r = {0, R h , R ± , ∞}. However, as pointed out in [85], under the change of variable and redefinition of the wave function the singularity at infinity is removed, and the equation becomes a canonical Heun equation (2.1) with (3.11) In the z coordinate, the horizon r = R h is mapped to z = t, the cosmological horizons r = R ± are mapped to z = 1 and z = ∞, respectively, while the origin, r = 0, is mapped to z = 0. The boundary conditions described for Φ imply the following behaviors for the function ψ: We now want to obtain the analytic formula from which the quasinormal modes can be computed in the limit where t is small, 0 < t 1, or, equivalently, R h is small, R h 1. For this purpose, we write the following dictionary for the gauge parameters in terms of Heun's parameters and gravitational quantities (see appendix A for the conventions used): , (3.13)

Connection Problem
The computation of quasinormal mode frequencies is obtained by imposing purely ingoing boundary conditions at the event horizon z = t and purely outgoing at the positive cosmological horizon z = 1. The independent solutions of the Heun equation for z ∼ t are given in (2.2), and the ones for z ∼ 1 are given in (2.3). Taking into account the boundary conditions (3.12), the connection coefficient between ψ (t) + and ψ (1) + has to be set equal to zero. The connection formula [12] 12 for small t from z ∼ t to z ∼ 1 is given by (3.14) This leads us to the quantization condition for the quasinormal modes in the form which can be rewritten as (3. 16) Note that this is nothing but (see Appendix A) where F full (t) is the full NS free energy, since the ratio of Gamma functions in (3.17) represents the 1-loop corrections.

QNMs at large
The previous quantization condition gets simplified in the large limit, where we neglect nonperturbative effects in of the form R h . This regime was studied for AdS 5 black holes in [23,87], since in this limit, the quasinormal mode frequencies become real, and, via the AdS/CFT correspondence, they compute the dimensions of certain operators in the holographic conformal field theory, see [88][89][90][91][92][93][94][95] and references therein. In the dS case, in this regime, the quasinormal mode frequencies are purely imaginary, and their interpretation from the point of view of holography is, at present, less clear (at least to us).
In the leading order in R h , a ∼ ± + 1 2 . Choosing the plus sign, the quantization condition since the other term is exponentially suppressed. This condition is satisfied if and only if

20)
12 For the computation of the traces of monodromies of the Heun equation, see also [9,86].
(3.22) Higher orders can also be computed systematically, but their expressions are cumbersome; hence we do not write them explicitly. Notice that in this limit, all the odd orders ω 2k+1 seem to vanish. Moreover, these formulas are correct for finite up to order R 2 +1 h , as will be shown in section 3.2. The pure de Sitter case can be obtained by taking the limit t → 0 or, equivalently, R h → 0. As the event horizon disappears in this limit, it is enough to consider only the region near the cosmological horizon r = R + . The leading order solutions in t of the canonical Heun equation (2.1) are given by

Perturbation theory around dS
where ω 0 is the leading order term in the R h expansion of the frequency: Since is a non-negative integer, the hypergeometric functions get truncated to polynomials as (3.25) The boundary conditions require that the radial part of the gravitational perturbation Φ (r) is well-defined as r → 0. Using the dictionary for the wave function (3.9), we rewrite the latter requirement in terms of ψ (z): Thus, we have to pick a regular solution at z ∼ 0 and consider an additional factor of z −s+1/2 . Looking at the first solution from (3.23), we can see that is not regular at z ∼ 0 for any allowed value of . Indeed, the other combination gives the solution, which is regular at z ∼ 0: In addition, the boundary conditions at the cosmological horizon require the eigenfunction to be regular with a well-defined Taylor expansion at z = 1. This is possible only if iω 0 ∈ Z ≥0 (due to the term (1 − z) iω 0 in (3.25)). Moreover, to avoid the poles in the Gamma functions in (3.25): we must exclude all the values of iω 0 that are smaller or equal to (these poles indicate that the second expression in (3.25) have to be rewritten in terms of log (z − 1) for iω 0 = , − 1, . . . , − + 1, − ). This gives the well-known quantization condition for the QNM frequencies of the pure dS 4 : The corresponding eigenfunction is We also note that the discarded solution is The Wronskian between f L 0 and g L 0 is

Left Region
Here we call the region near the cosmological horizon r = R + left region due to the analogy with the corresponding quantum mechanical problem on the complex plane. The local variable in this region is z, and the leading order solutions in R h (and so in t) of the canonical Heun equation (2.1) are given in (3.30), (3.31). Expanding in small R h the solution and the frequencies we get for the outgoing solution ψ where ω 1 is a coefficient in the R h expansion of the frequency (3.24). Since g 0 (z) blows up as z → 0, it should not be present in the leading order of the wave function in the left region. Hence, we require ω 1 = 0. On the other hand, the incoming wave solution at the cosmological horizon is After we fix ω 1 = 0 and proceed with the general method described in Sec. 2.2, the logarithm function log (z − 1) appears in higher orders in R h (and t). The only source of this function is the incoming wave solution (3.34), and we will be canceling any contributions of log (z − 1) by fixing the coefficients b K in the perturbative expansion of the wave function (2.6), (2.8).
After establishing the boundary condition, we compute the integrals in (2.8). As we show in Appendix C, these integrals are described in terms of the multiple polylogarithms in a single variable: The latter admits for s 1 ≥ 2: and for s 1 = 1, k ≥ 2: The weight of the multiple polylogarithm Li s 1 ,...,s k (z) is s 1 + · · · + s k , and the level is k. At each order t K+1 , both integrands in (2.8) are linear combinations of the following terms with maximum weight K: where r 1,2 , i 1,2 , j 1,2 , p 1 are some non-negative integers, and 0 ≤ p 1 ≤ K, s 1 + · · · + s k ≤ K. After taking the integrals, the only new functions that appear are multiple polylogarithms of maximum weight K + 1. Moreover, both integrals in (2.8) are linear combinations of terms similar to (3.38): and terms containing new combinations of logarithms and multiple polylogarithms that were not present in (3.38): where the maximum weight is K + 1:ŝ One of the differences between (3.38) and (3.39) is that r 1,2 , i 1,2 , and j 1,2 are shifted by 1 or −1.
These shifts are specific to the left region of the SdS 4 case (and even then may be subjected to reevaluation for some values of quantum numbers n, , and s that we did not consider). In the right region, the shifts are different but can be determined on the case by case basis (the details can be found in the attached Mathematica files, where we distinguish two regimes with ≤ n and > n). Even though the optimal choice of shifts depends on the case at hand, there is a choice of big enough shifts applicable to all quantum numbers for both regions.
To summarize, we reduced the problem of solving the initial ODE in a given order in t to a system of linear equations on the coefficients in front of the functions from (3.40) and γ m , δ m . 13 The resulting corrections f L K (z) to the wave function in the left region are linear combinations of the following functions: where k 1,2 , l 1,2 , p 1 are some non-negative integers, 0 ≤ p 1 ≤ K, s 1 + · · · + s k ≤ K, and ζ L m , ξ L m are z-independent quantities.

Right Region
The right region is near the event horizon r = R h , or z = t. We introduce the local variable z R = t/z so that the horizon is at z R = 1. In the z R variable, the equation (2.1) reads (3.43) In the remaining part of this subsection, we will mostly omit the R index on the z variable (except for the cases where it could be confusing). We take as leading order solutions in R h (and so in t) of this equation The Wronskian between f R 0 and g R 0 is .
Here we would like to comment on the choice of the logarithm function log 1 − z R in the solution g R 0 . The other possible choice of the logarithm could be, for example, log z R − 1 . This choice dictates what functions will appear in higher orders in t and affects the R h expansion of the frequency ω. Throughout the paper, we work with the principal value of the complex logarithm, and thus the change in the argument affects the position of the branch cut on the complex z plane. Our wave function ψ (z) can be viewed as an analytic continuation of the physical solution on half of the real line r ≥ 0. In the de Sitter case, the coordinate transformation z (r) is (3.8) with real parameters R ± . Since we want the solution to be continuous across the real slice R h < r < R + , the branch cut should not cross the interval t < z R < 1, where t is small and positive. This leaves us with log 1 − z R , and the branch cut runs from z R = 1 to z R = +∞. The other logarithm function that appears in higher orders in t is log z R , and the corresponding branch cut runs from z R = 0 to z R = −∞ also avoiding the interval t < z R < 1 (see Figure 1).
The boundary condition near the horizon requires us to keep the solution corresponding to the incoming wave and discard the one corresponding to the outgoing wave. According to (2.2), the two solutions behave like (3.47) , both waves in the rhs of (3.47) have Taylor expansions in R h that start with 1. One must also consider the higher orders in R h to distinguish the two expansions. The incoming wave solution has a particular dependence on the logarithm function log (1 − z) in each order in R h (or t): In the leading order in R h both ψ (t) − and ψ (t) + are given by the same function f R 0 (z). Since the other function g R 0 (z) contains the logarithm, it enters ψ (t) + in the higher orders in R h . The constants b K from (2.8) are fixed by matching with the logarithmic behavior of the incoming wave solution (3.48) in each order in R h . The integrals in (2.8) are again described in terms of the multiple polylogarithms in a single variable (see Appendix C). We construct the linear basis of functions for each integral in the way it was done in the previous section for the left region. The only difference is that we need to add powers of the second logarithm function log (1 − z) to formulas (3.38), (3.39) and (3.40). In particular, the second integrand from (2.8) at order t K of the form will have a maximum weight K because the logarithm function log (1 − z) is present in the leading order solution g R 0 (z). The resulting integral, however, will be of the same weight K due to the pole structure in (3.49). Eventually, the corrections f R K (z) to the wave function in the right region are linear combinations of the following functions of maximum weight K: where k 1,2 , l 1,2 , p 1,2,3 are some non-negative integers, and 0 ≤ p 1 + p 2 ≤ K, p 3 + s 1 + · · · + s k ≤ K.

Results for QNM frequencies
The final step in the procedure described in Section 2.2 is to glue the local solutions by requiring that the wave function and its first derivative are continuous at the intersection of the two regions. There is a certain freedom in choosing the intersection point as long as it lies in the region of convergence of both local solutions. We choose the point z = t 1/2 , which is the same as z R = t 1/2 . Note that the expansions of ψ L,R z L,R are given as series expansions around z L,R = 1 up to orders t m L,R : What happens when we take z L,R ∼ t 1/2 and expand for a small t? Some terms f L K (z)t K in ψ L (z) will contribute to orders lower than t K . This could lead to a reshuffling, where, for example, f L 1 (z)t becomes the leading order contribution at z ∼ t 1/2 . This happens when ≥ 1, as seen from (3.30): However, since we are within the radius of convergence of ψ L (z), this reshuffling involves only a finite number of terms. For all values of quantum numbers we have considered, the reshuffling is superficial and goes away after the frequency is set to one of the quasinormal modes. The continuity condition can be equivalently stated as where C(t) is a normalization factor. The advantage of (3.54) is that we can use one of the equations to understand which orders in t we can trust when expanding (3.51) at z L,R = t 1/2 , then use the other one to fix the frequencies.
Using Mathematica, we compute the local solutions up to orders m L = 10 and m R = 7. This allows us to determine the R h expansion of the frequency up to order R 9 h or less depending on the value of . In all computed orders, we find the real part of the quasinormal modes is zero, which agrees with the earlier observations made by numerical computations [66][67][68]. The results for the imaginary part of the quasinormal mode frequencies ω n, ,s , starting from n = 0, are Let us also report the results for n = 1: Some of the results presented above were shortened for the reader's convenience. The full expressions and more expansions of frequencies for other choices of and s can be found in the attached Mathematica files. The irrational numbers entering these QNM frequencies are log(2) and multiple zeta values.

Perturbations of Anti-de Sitter black holes in four dimensions
The metric describing the AdS 4 Schwarzschild black hole is given by (3.1), with Λ < 0. We denote the roots of rf (r) = 0 by where, for Λ < 0 , R ± are complex conjugate and given by in terms of the BH horizon R h ∈ R >0 . We will fix Λ = −3 and study the same perturbations of the Schwarzschild de Sitter case, described by equation (3.5). According to AdS 4 /CFT 3 holography, the conformally coupled scalar field is dual to scalar operators of conformal dimension ∆ = 1 or ∆ = 2, from the relation µ 2 = ∆(∆ − 3). The main difference with the SdS 4 case lies in the boundary conditions we impose on the solution. Indeed, we will still require the presence of only ingoing modes near the horizon, but we will impose the vanishing Dirichlet boundary condition at the AdS boundary. 14 With the following change of variables and redefinition of the wave function the singularity at infinity is removed, and the equation (3.5) becomes canonical Heun equation In these coordinates, the horizon is at z = t while the boundary is at We will also consider the small black hole limit, R h 1. The Dirichlet boundary conditions in terms of the ψ function are given by Notice that the AdS boundary (z = z ∞ ) is not a singular point of the perturbation equation. This makes the approach based on the Seiberg-Witten theory less effective. One can write the quantization condition using the connection formulae between Heun functions, but in this case, an expansion of the Heun functions in R h is needed. We will report some results in this direction in appendix E.

QNMs in pure AdS 4
The pure AdS 4 case can be recovered in the limit t → 0 or, equivalently, R h → 0. In this limit, the z variable is given by 9) and the AdS boundary is at z = 2. The leading order solutions in t of the canonical Heun equation (2.1) are given by where ω 0 is the leading order term in the R h expansion of the frequency: As in the de Sitter case, these hypergeometric functions reduce to (3.25), where we replace −iω 0 by ω 0 . The first boundary condition from (4.8) tells us that the wave function ψ (z) is regular at z = 0. This singles out the second solution from (4.10). Then, the second boundary condition at z = 2 requires the following expression to vanish: which gives the quantization condition for the QNM frequencies of the pure AdS 4 Here we have two branches of frequencies, positive and negative, and one is related to another by the complex conjugation of the radial part of the perturbation Φ (r).
In the following subsections, we will perturb around the pure AdS case to obtain the corrections for the Schwarzschild Anti-de Sitter small black holes. Following the same logic as in the de Sitter case, we will divide the space into two regions: left (L) and right (R). The left region describes the physical space near the AdS boundary with r → ∞, and the right one is the space near the horizon r = R h . After having determined the expansion of the solution ψ(z) in each region up to certain orders in the expansion parameter t, we require that the function ψ (z) and its first derivative are continuous in a point in the intersection of two regions, which we can fix at z = t 1/2 (other values of z are possible as long as they lie inside the convergence radius of the two solutions).

Left Region
The local coordinate in the left region is z, and the AdS boundary is at z ∞ , which has the following expansion in R h : (4.14) The wave function in the left region ψ L (z) satisfies the same Heun equation (2.1). The form of the leading order solutions depends on which branch of frequencies we choose in (4.13). For the negative branch ω 0 = − − 2n − 2, we have and for the positive branch ω 0 = + 2n + 2: (4.16) For both branches, the Wronskian can be written in terms of ω 0 as We will apply the perturbative method described in Section 2.2 to both positive and negative values of ω 0 , but the final result is straightforward. The only difference between the two branches is the sign of the real part of the frequency expansion (4.11), which again corresponds to complex conjugation of Φ (r). The boundary condition in the left region is simply ψ L (z ∞ ) = 0. Since f L 0 (2) = 0 and g L 0 (2) = 0, we get the following perturbative expansion for the wave function in the left region: where f L K (z) are given by (2.8). The constants b K in (2.8) are fixed by expanding ψ L (z ∞ ) in powers of t and requiring the coefficients in this expansion to vanish.
As we explain in Appendix C, the integrals in (2.8) are described in terms of the multiple polylogarithms in a single variable (3.35). Since the weights of the multiple polylogarithms appearing at order t K are less or equal to K, we can construct a linear basis of functions in which the integrals in (2.8) can be expanded. We take the same steps (3.38)-(3.40) as we did in the SdS 4 case to do this. The only difference is that we add the second logarithm function log (z − 1) to (3.38). To be more precise, the integrands in (2.8) at order t K+1 are given by the linear combination of the following functions: where r 1,2 , i 1,2 , j 1,2 , p 1,2,3 are some non-negative integers and p 1 + p 2 ≤ K, p 3 + s 1 + · · · + s k ≤ K. The reasoning behind our choice of the branches of the logarithm functions log (z) and log (z − 1) is the same as in Section 3.2.3. We want the wave function ψ (z) to be continuous across the real slice R h < r < +∞. In the SAdS 4 case, the coordinate transformation z (r) is given by (4.3) with complex parameters R ± . Taking into account that r and R h are real, we have from (4.3): (4.20) Thus, the real slice is approximately half the circle with the center in z = 1 on the complex z plane (see Figure 2). It starts at z = t and ends at z = z ∞ . Simple analysis shows that Im (t) > 0 and Im (z ∞ ) > 0 when R h > 0. This justifies our choice of logarithm functions since both branch cuts do not cross the real slice. On the other hand, if one picks log (1 − z) instead of log (z − 1), the corresponding branch cut would touch the real slice at the point z = 2 when evaluating ψ L (z ∞ ). This, in turn, would lead to incorrect results for QNM frequencies.

Right region
In the right region, we introduce local coordinate The horizon is now situated at z R = 1. The wave function in the right region ψ R z R satisfies the following equation in terms of z R : Suppressing the R index on z R , the two leading order solutions are given by where the constants a s m , b s m can be determined for any ≥ s ≥ 0 as The expressions in (4.23) are independent of which branch of frequencies we choose in (4.13) because the leading order of (4.22) does not contain ω 0 . The Wronskian of f R 0 and g R 0 is given by . (4.24) The boundary condition in the right region tells us that ψ R is regular at z R = 1. Thus, we can write the following perturbative expansion: where f R K (z) are computed using (2.8). Unlike in the left region, the choice of the logarithm function in g R 0 (z) is unimportant. This is due to the boundary condition that requires canceling contributions of log (1 − z) in each order t K . The resulting corrections f R K (z) are linear combinations of the following functions of maximum weight K: where k 1,2 , l 1,2 , p 1 are some non-negative integers, and 0 ≤ p 1 ≤ K, s 1 + · · · + s k ≤ K.

Results for QNM frequencies
To determine the QNM frequencies, we use the continuity condition in the form (3.54): where ψ L,R z L,R are computed up to orders m L,R in t around z L,R = 1: (4.28) Similarly to the SdS 4 case, the reshuffling of terms (3.52) occurs in ψ L (z) when we take z ∼ t 1/2 . For all values of quantum numbers we have considered, this reshuffling is superficial and goes away after the frequency is set to one of the quasinormal modes.
Using Mathematica, we compute the local solutions up to orders m L = 7 and m R = 8 (sometimes even up to m L = 9 and m R = 10). This allows us to determine the R h expansion of the frequency up to order R 7 h or less depending on the value of . In all computed cases, the imaginary part does not appear before order 2 + 2 in R h : (4.29) As mentioned, the results computed for negative and positive branches of ω 0 only differ by the sign in the real part of the frequency expansion. Below are the results for the real and imaginary parts of the quasinormal mode frequencies ω n, ,s corresponding to the positive branch, starting from n = 0: For n = 1 we have: Some of the results presented above were shortened for the reader's convenience. The full expressions and more expansions of frequencies for other choices of n, , and s can be found in the attached Mathematica files. From these, one can see that the irrational numbers entering these QNM frequencies are log(2), π, and Euler sums. Analytically computing f L 1 (z) from (4.18), we can also determine the subleading term in the QNM frequency expansion with n = 0 and ≥ 1: For small enough values of R h , our results agree with the numerical ones obtained earlier in [69].
Since the frequency expansions in higher orders in R h include multiple zeta values (B.11), we use different identities of the form (B.12)-(B.16) to compute the corresponding numerical values. Tables 1-3 present the numerical results from the frequency expansions truncated at R 7 h (in the scalar case with n = l = 0, the expansion was computed up to order R 6 h and truncated at the same order). In these tables, bold digits are the ones that are stable and agree with the numerical results obtained directly from the Heun function and the continuity condition (3.53

Scalar Sector of Gravitational Perturbations -The low-lying modes
Following [96], one can consider a subdivision of gravitational perturbations in different sectors (scalar, vector, or tensor), whose distinction comes from the expansions in scalar, vector, or tensor spherical harmonics on the S 2 component of AdS 4 . In Sec. 4 we considered the vector sector of gravitational perturbations (s = 2). We will now focus on the scalar sector and impose a new boundary condition at the AdS boundary, namely a Robin boundary condition [73,[97][98][99][100][101], see also [102] for very recent developments. This choice of boundary condition is motivated by the AdS/CFT correspondence, and it ensures that the perturbations do not deform the metric on the boundary of AdS. From the point of view of the dual CFT, these boundary conditions are related to double-trace deformations, see for instance [103][104][105] and references therein. In particular, we will analyze the so-called low-lying quasinormal frequencies, which, according to AdS/CFT duality, are related to hydrodynamic modes of the 3d thermal CFT on the boundary [74,101,[106][107][108][109][110][111][112][113]. We will therefore expand our quasinormal frequencies for large values of R h , R h 1, differently from the previous sections. Defining This equation has five regular singularities, located at r = 0, R h , R ± , R 5 , where The new singularity R 5 , coming from the potential V S (r), is in the unphysical region r < 0.
Similarly to the previous cases, we introduce the change of variables z (r) = R h r (5.5) and the new wave function ψ (z) = r −1 e iωr * Φ (r) . (5.6) The master equation (5.2) then becomes and M is related to R h via The boundary conditions in terms of the ψ function are given by (5.10) The five regular singularities of the equation (5.2) have three different scalings with R h → ∞. The singularity at r = 0 doesn't scale, the singularities R ± and R h scale linearly, and R 5 scales as R 3 h . Hence, we will divide the space into three different regions and apply the perturbative method described in Section 2.2.
The three local variables are x = R 3 h /(mr) + 1/3 for the left region (near the AdS boundary), y = R 2 h /r for the middle region, and z = R h /r for the right one (near the BH horizon). 15 Here the regions are labeled left and right as they appear on the complex z plane (see Figure 3). From the point of view of the complex z plane, the left and middle regions represent two zoomings close to the origin, with different scalings. Considering the normal form of the differential equation The two rescalings x ∼ R 2 h z m and y = R h z are such that, in both variables, the differential equation in normal form has a potential, V x (x) and V y (y), respectively, with non-vanishing leading order Out of the three, the right region is the one in which it is more challenging to expand the solution of the differential equation. In particular, the solution involves multiple polylogarithms in several variables, which we analyze in Appendix D.
Since we will work with R h 1, the small parameter is α = 1/R h , and the frequency expansion can be written as 14) The intersections of the three regions and the boundary points r = R h , ∞ determine three intervals in which the wave function should be continuous:  From the point of view of x and y, the first two intervals have infinite lengths, their left endpoints are at finite values and their right endpoints are chosen to meet the next region (and so they become infinite because of the different scalings of the local variables in powers of R h ). Finally, we will derive the low-lying QNM frequencies by requiring that the wave function and its first derivative are continuous at the intersection points y = 1 and z = α 1/4 . As we explain later, the second intersection point z = α 1/4 is chosen to avoid the reshuffling of terms in the wave function expansion (5.32).

Left Region
The left region represents the region close to the AdS boundary, where we impose the Robin boundary condition. The local variable in this region is 16) and the AdS boundary is at x = 1/3. The master equation in the left region is obtained by applying the coordinate transformation z = α 2 m (x − 1/3) to (5.7) and substituting ψ (z) with ψ L (x). In the leading order in α, we get The two leading order solutions are Since f L 0 satisfies the Robin boundary condition d dx the following perturbative expansion for the wave function in the left region can be written: We do not use (2.8) to compute f L K (x) as they are simple Laurent polynomials in x. The form of these polynomials depends on whether K is even or odd. The following general result holds for the first 30 computed orders: a 2K−1,s x s , where the coefficients a K,s depend on the parameters m and ω i . For example, we have for K = 1, 2, 3, 4: In each order in α, the contribution of g L 0 is fixed by the Robin boundary condition. The contribution of f L 0 is arbitrary and can be absorbed into a normalization of the wave function ψ L (x). We choose the normalization so that f L 0 is only present in the leading order.

Middle Region
To match the wave function expansions in the left and right regions, we introduce an intermediate region with the local variable The master equation in the middle region is obtained by applying the coordinate transformation z = α y to (5.7) and substituting ψ (z) with ψ M (y). In the leading order in α, we get The two leading order solutions are Strictly speaking, there is no boundary condition in the middle region. However, there is a way to use the expansion of the wave function in this region and apply the boundary condition near the horizon y ∼ α −1 . This requires a resummation of infinitely many terms, and the results agree with the ones obtained using three regions instead of just two. Here we focus on the procedure with three regions as it allows us to get more orders in the QNM frequency expansion. To justify our choice of functions f M 0 and g M 0 , we can either use the gluing procedure or look at the behavior near the horizon. In the first couple of orders in α, there is no resummation of terms in the wave function ψ M (y) when we take y ∼ α −1 . Since near the horizon g M 0 (y) ∼ α −3 , it can only appear in orders α 3 and higher. This leads to the following perturbative expansion of the wave function: Similarly to the left region, the corrections f M K (y) are Laurent polynomials of the form where coefficients b K,s also depend on the parameters m and ω i . Starting from order α 3 , the gluing procedure fixes the contribution of g M 0 , so we keep the corresponding integration constants c M K in the expressions for f M K , K ≥ 3. Out of the 27 computed orders, we present the first 4: (5.28)

Right Region
The local variable in the right region is z, and the event horizon is at z = 1. The leading order in α of (5.7) is The two leading order solutions are The Wronskian between these solutions is According to the boundary conditions (5.10), the wave function in the right region is regular at z = 1. The corresponding perturbative expansion of the wave function is then and for s 1 = 1, k ≥ 2, The weight and level of Li s 1 ,...,s k (z 1 , . . . , z k ) are s 1 + · · · + s k and k. When taking the integrals in (2.8) with the input from this section, we will only encounter multiple polylogarithms with s 1 = s 2 = · · · = s k = 1 (see Appendix D for more details). In this case, the weight and level are the same. Moreover, all arguments z i with i ≥ 2 are constants and can take one of the three possible values: 1, u 1 , and u 2 . These constants are the third roots of unity that arise in the following decomposition of g R 0 (z): Similarly to the previous cases with multiple polylogarithms, the corrections f R K (z) at order α K are described in terms of functions Li s 1 ,...,s k (z 1 , . . . , z k ) of weight K and lower. This allows us to construct a linear basis of functions, in which f R K (z) can be expanded: where i 1,2 , j 1,2 , k 1,2 , l 1,2 , p j are non-negative integers, and 0 ≤ p 1 + p 2 + p 3 ≤ K, 0 ≤ p 4 + p 5 + p 6 + k ≤ K. Since the first argument in Li {1} k (z 1 , z 2 , . . . , z k ) can take one of the three possible forms z 1 = z, z 1 = u 1 z, or z 1 = u 2 z, (5.39) we have 3 k functions that can enter the basis at level k ≥ 2. However, this number is reduced due to the identities that involve multiplication by ordinary logarithm functions log (1 − z), log (1 − u 1 z), and log (1 − u 2 z) (see Appendix D). These identities allow us to use only two forms of the first argument z 1 = u 1 z and z 1 = u 2 z. The reduced number of multiple polylogarithms that enter the basis is 8 × 3 k−3 for k ≥ 3, and just 3 for k = 2: Using Mathematica, we compute 7 corrections f R K (z); the first two are We estimate the following behavior of f R K (z) as z → 0 based on the obtained results: Thus, to avoid the reshuffling of terms, we choose the gluing point between the middle and the right region to be z = α 1/4 .

Results for QNM frequencies
We need two continuity conditions to determine the QNM frequencies, at z = α 1/4 and z = α: The first condition in (5.43) is used to fix the integration constants c M K , and the second one gives the coefficients ω k in the QNM frequency expansion (5.14). The first seven computed orders of the wave function expansion in the right region allow us to determine ω k up to k = 6: (3) , (3) , where we shortened the results for ω 4 and ω 5 for readers convenience. The full results, including the result for ω 6 , can be found in the attached Mathematica files. Notice that, as compared to the QNM frequencies computed in Sec. 3 and Sec. 4, here the frequencies involve different irrational numbers, for instance, log 3, √ 3, as well as colored multiple zeta values of level 3. Upon taking the scaling limit where q stays constant, we reproduce the results for the QNM frequencies of the M2-brane in the AdS 4 background (see Table IV in [74]) which are directly linked to hydrodynamics [106][107][108][109]. Also, the following rescaling of the frequency is needed: Applying this limit to (5.44), we obtain an expansion of w in q: where w 1 , w 2 , and w 3 agree with the results from [74], and the new results are where we shortened the results for w 5 and w 6 for readers convenience. The full results, including the result for w 7 , can be found in the attached Mathematica files. The numerical values of these coefficients are  These alternate between real and imaginary parts, precisely as predicted in [101,110]. 16

Conclusions
This paper focuses on analytical aspects of spectral problems associated with perturbation theory for four-dimensional (A)dS black holes. We explore these problems using two analytic strategies: one based on the NS functions and one based on a recursive structure involving multiple polylogarithms. Thanks to these tools, we can compute the quasinormal mode frequencies and their eigenfunctions analytically in various regimes. For instance, we can obtain the series expansion at large R h , small R h , or large spin (R h being the BH horizon).
We use the approach based on the NS functions in the context of four-dimensional dS Schwarzschild black holes. In this setup, the NS functions allow us to compute the large expansion of QNMs systematically. We find that, up to non-perturbative effects in , the QNMs are (negative) imaginary numbers that are even functions of R h . To include non-perturbative effects in the spin, switching to the polylog approach is convenient. Once non-perturbative effects are included, QNMs are no longer even in R h . But we still find a branch of purely imaginary modes, thereby providing analytical confirmation of the results obtained through numerical studies in [66][67][68]114]. Exploring the interplay between the NS and polylog approaches would be interesting. In particular, the appearance of multiple polylogarithms and multiple zeta values may be related to the behavior of the NS functions close to their singular points, see e.g. [115][116][117].
We extend the polylog method to study conformally coupled scalar, electromagnetic, and vector-type gravitational perturbations in asymptotically AdS 4 Schwarzschild black holes. The NS functions are less effective for these perturbations because the point at spatial infinity is not a singular point of the equation. 17 Hence, we switch to the polylog method for Dirichlet and Robin boundary conditions. As an application, we use this technique to study the low-lying modes of the scalar sector of gravitational perturbations and compute several orders in the 1/R h expansion. Even in the hydrodynamic expansion, this allowed us to go beyond the results presently available in the literature. From the point of view of holography, the polylog method presents finite spin predictions for the dual 3d CFT. It would be interesting to explore this further in higher spacetime dimensions and make contact with past and recent developments in the study of holographic CFTs [23, 87-95, 103-105, 109, 119].
The technical result we obtained about the perturbation theory of second-order linear differential equations with Fuchsian (or irregular 18 ) singularities points to the existence of a recursive structure for their solution involving multiple polylogarithms. It raises the question of whether there exists a deeper algebraic structure beyond this that could improve the algorithm. This should allow us to have a higher level of analytic control over the problem at hand. For instance, it would be interesting to quantify the precise analytic properties of QNM frequencies as functions of the BH radius and/or other relevant parameters to understand their physical meaning better. For example, this would allow to detect phase transitions and/or (in)stabilities. These considerations become especially interesting when considering rotating and/or charged black holes. Indeed, these exhibit a richer structure with intriguing (in)stability features. It would be interesting to revisit these problems within the approaches presented in this paper.
Let us remark that the polylog method we developed shares similarities with the techniques used to compute Feynman integrals in Quantum Chromo-Dynamics, see e.g. [80][81][82] and references therein. 19 Similar techniques also recently appeared in studying higher curvature corrections to the effective low-energy gravitational theory arising from string scattering diagrams [121][122][123]. [124][125][126][127]. Although these results are directly related to hyperbolic trajectories, one can use the data extracted from the amplitudes to determine the parameters of the effective one-body potential of [128] to be used to describe the gravitational bound states. From the computational viewpoint, the resulting polylog expansion in such approximation is naturally obtained by com- 17 If we consider massive scalar perturbation instead, the underlying equation has five regular singular points and spatial infinity is mapped to one of them. In this case, one can use the NS functions for an SU (2) × SU (2) linear quiver [118]. In addition, for generic mass, the leading order solution does not reduce to a rational function; hence the polylog approach presented can not be applied straightforwardly. 18 A detailed analysis of the case with irregular singularities, which is relevant for asymptotically flat black holes, will appear in [83]. 19 During the writing of this paper, we were informed that Saso Grozdanov is also exploring similar ideas in the context of AdS5 black holes (particularly the hydrodynamic limit)[120].
puting the relevant Feynman multiloop integrals in the proper kinematic regime. On the contrary, in the QNMs regime, the appearance of multiple polylogarithms does not seem to have a direct interpretation in terms of Feynman multiloop integrals. Moreover, different types of special functions arise for other gravitational backgrounds and/or other perturbations. For example, when considering asymptotically flat black holes, there is the appearance of multiple polyexponential functions as well [83]. It would be interesting to understand this better.
Finally, one of the most challenging and exciting questions would be to go beyond linear perturbation theory. The NS and polylog methods allow for the computation of the eigenfunctions and the Green functions, which are essential inputs to go beyond the linear theory.

A NS functions
This appendix reports the notations and conventions used in Sec. 3.1, where the gauge theory approach is applied to the Heun connection problem. The relevant theory is N = 2 SU (2) gauge theory with N f = 4 fundamental hypermultiplets.
If Y is a Young diagram, we denote with (Y 1 ≥ Y 2 ≥ . . . ) the heights of its columns and with (Y 1 ≥ Y 2 , . . . ) the lengths of its rows. For every Young diagram Y and for every box s = (i, j), we denote the arm length and the leg length of s with respect to the diagram Y as Note that we do not require s to be in Y : if this is the case, the arm length and the leg length are non-negative quantities, but this is not true in general. We now introduce the main contributions coming into play for the definition of the instanton partition function of N = 2 SU (2) gauge theory with fundamental matter. Let us denote with Y = (Y 1 , Y 2 ) a pair of Young diagrams and with | Y | = |Y 1 |+|Y 2 | the total number of boxes. We denote with a = (a 1 , a 2 ) the v.e.v. of the scalar in the vector multiplet and with 1 , 2 the parameters characterizing the Ω-background. We define the hypermultiplet and vector contribution as [129,130] z hyp a, Y , m = k=1,2 (i,j)∈Y k .

(A.2)
We will always take 1 = 1 and a = (a, −a). Let us denote with m 1 , m 2 , m 3 , m 4 the masses of the four hypermultiplets and let us introduce the gauge parameters a 0 , a t , a 1 , a ∞ satisfying Moreover, we denote with t the instanton counting parameter t = e 2πiτ , where τ is related to the gauge coupling by In the Mathematica programs available at https://github.com/GlebAminov/BH_PolyLog, we also use the redefined masses M i , which are related to m i via The instanton part of the NS free energy is then given as a power series in t by In the text, we will also refer to the full NS free energy, which contains not only the instanton part but also the classical and one-loop contributions. This is explicitly given by The gauge parameter a is expressed in a series expansion in the instanton counting parameter t, obtained by inverting the Matone relation [131,132] where the parameter u (0) is the complex moduli parametrizing the corresponding SW curve. Explicitly, the expansion reads as follows (A.10)

B Useful facts about multiple polylogarithms in a single variable
There are many identities between polylogarithms and multiple polylogarithms. Below is the list of identities that are relevant in our case. First, for multiple polylogarithms of the form Li 1,s 2 ,...,s k (z), we have: Taking derivatives and using (3.36) and (3.37), it is easy to show by induction that n ≥ 1 : Generalizing the last two identities to an arbitrary level, one gets the following identity, which we use to express Li 1,s 1 ,...,sn (z) in terms of multiple polylogarithms Li r 1 ,...,r n+1 (z) with r 1 ≥ 2: Up to weight 4, all multiple polylogarithms in a single variable can be expressed as ordinary polylogarithms by combining the above identities and the following ones [133,134]: There are identities for weight higher than 4, but not enough to express all multiple polylogarithms as ordinary polylogarithms. For example, we have for weight 5: (B.10) The latter can be checked by taking a derivative and using identity B.1. Throughout the paper, we choose not to use the powers of polylogarithms in any basis, which reduces the number of relevant identities.
Again, the integration constants c m,n can be computed with the help of the known identities: In the same way, one can derive the expansion for Li m,1 by consecutively integrating Li 1,1 :

C Solving integral recurrence relations
In sections 3.2.2 and 4.2, we claimed that the wave functions ψ L (z) at order t K (or, equivalently, R K h ) are described in terms of multiple polylogarithms of weight K and lower. In this section, we will prove this claim, but first, let us clarify the terminology. The notion of weight is related to the power of a logarithm function, as seen in the following identity: Thus, we will ascribe weight to the ordinary logarithm functions as follows. For any product of two logarithms m, n ≥ 0 : log (z) m log (z − 1) n , (C.2) the weight equals n + m ≥ 0. For the product of a logarithm and a multiple polylogarithm k ≥ 1, n ≥ 0 : log (z − 1) n Li s 1 ,...,s k (1 − z) , the weight is n + s 1 + · · · + s k > 0. Here we do not consider the other possible product log (z) m Li s 1 ,...,s k (1 − z) because, due to the identities of the form (B.4), this product can always be rewritten as a linear combination of multiple polylogarithms. Some simple examples are: In general, multiple polylogarithm functions can not be rewritten as powers of ordinary logarithm functions. We will use both logarithms and multiple polylogarithms of a certain weight to build a linear basis in which the wave function can be expanded at a certain order in t. In what follows, all powers of logarithms are non-negative integers.
We are going to prove our claim by induction. In the first order in t, the integrands in the recurrence relations are just rational functions of the form Step 1. Four types of integrals increase the maximum weight. In each case the integrand has a factor of z −1 or (z − 1) −1 . For the product of two logarithms of weight n + m, n, m ≥ 0 we have: (C.11) where in the last integral m ≥ 1 and Li n,{1} 0 ≡ Li n . The resulting weight after the integration is 1 + n + m. In the more general case of integrals involving multiple polylogarithms, we have (n ≥ 0): (C.13) Again, after the integration, the weight was increased by 1 from n+s 1 +· · ·+s k to 1+n+s 1 +· · ·+s k . The above identities were obtained by repeated integrations by parts.
(D. 20) Let us remark that with the previous procedure, one can find several identities at a fixed level involving the same multiple polylogarithm. To choose which elements to add to the basis, we followed the criterium that we omit the multiple polylogarithms with the first argument z. This criterium comes from the regularity condition on the wave function at z = 1. Moreover, when possible, we tried to include the same number of multiple polylogarithms with the first argument u 1 z and with the first argument u 2 z (for example, it is not possible at level k = 2).