Holographic thermal correlators: A tale of Fuchsian ODEs and integration contours

We analyze real-time thermal correlation functions of conserved currents in holographic field theories using the grSK geometry, which provides a contour prescription for their evaluation. We demonstrate its efficacy, arguing that there are situations involving components of conserved currents, or derivative interactions, where such a prescription is, in fact, essential. To this end, we first undertake a careful analysis of the linearized wave equations in AdS black hole backgrounds and identify the ramification points of the solutions as a function of (complexified) frequency and momentum. All the equations we study are Fuchsian with only regular singular points that for the most part are associated with the geometric features of the background. Special features, e.g., the appearance of apparent singular points at the horizon, whence outgoing solutions end up being analytic, arise at higher codimension loci in parameter space. Using the grSK geometry, we demonstrate that these apparent singularities do not correspond to any interesting physical features in higher-point functions. We also argue that the Schwinger-Keldysh collapse and KMS conditions, implemented by the grSK geometry, continue to hold even in the presence of such singularities. For charged black holes above a critical charge, the energy density operator does not possess an exponentially growing mode, associated with `pole-skipping' (from one such apparent singularity). Our analysis suggests that the connection between the scrambling physics of black holes and energy transport has, at best, a limited domain of validity.


Introduction
The study of linearized perturbations of black holes, initiated in [1,2], is relevant for a broad class of physical questions. These range from astrophysical signals of black hole horizon ringdown, to mathematical relativity questions of black hole stability, and holography [3]. The linear wave equations have therefore been analyzed by many authors over the years. 1 Our interest is in the context of holography, especially in the analytic structure of the solutions, and how this informs boundary correlation functions, in particular higher-point functions.
The context for our analysis is the computation of real-time thermal correlators in a strongly coupled holographic field theory. We will focus on Schwinger-Keldysh correlators, which have been argued to be computed on a particular complexified spacetime in [7] (we refer to the resulting geometry as the grSK spacetime). As it stands, their prescription can be viewed as a particular contour choice for the computation of bulk Witten diagrams, as elaborated in [8][9][10]. The virtue of having a contour prescription is that one can fold in the information about the nature of local behaviour into the contour integral. This turns out to be especially important when considering correlation functions of conserved currents where, as we shall see, the bulk vertices involve functions with singularities which potentially interfere with the contour of integration. In turn, these new singularities could, in principle, break the general arguments advanced in the above references that correlators computed in the grSK spacetime satisfy the general properties required of Schwinger-Keldysh thermal correlators.
To be clear, the analysis of real-time correlation functions in the holographic context has a long and rich history. The original prescription for computing two-point functions was given in [11] and justified in [12]. The logic here was to argue that the retarded correlators, by virtue of boundary causality, should involve infalling boundary conditions on the future horizon and suitable boundary sources controlled by non-normalizable modes at the AdS asymptopia. 2 The boundary two-point function was obtained as the ratio of the normalizable mode to the non-normalizable mode for such solutions. Much of the subsequent literature in applications of AdS/CFT to compute thermal real-time observables, spectral functions, fluid-gravity correspondence etc., has relied on this prescription.
While the prescription of [11] captures correctly the retarded observables, owing to the presence of a single boundary condition on the asymptopia, it does not naturally lend itself to computing fluctuations (one could, however, get these by demanding that the fluctuationdissipation theorems hold at the Gaussian order). A more natural prescription therefore was the proposal of [13,14] who argued that one should take the boundary Schwinger-Keldysh contour seriously. In particular, they posited that the bulk spacetime dual to such contours should be constructed by piecewise gluing Lorentzian and Euclidean geometries across codimension-1 bulk surfaces. Within this framework, [15] demonstrated the naturalness of the infalling prescription of [11]. Taking this into account, one can indeed compute (higher-point) correlators with general Schwinger-Keldysh time-ordering, as was done already a decade ago in [16,17]. More recent efforts in this direction include [18,19] for non-thermal excited states, and [20] who computed thermal 3-point functions.
One might thus ask, what does the prescription of [7] buy us, apart from perhaps the elegance associated with writing the result in terms of a contour integral? For instance, using this prescription it was demonstrated that a general correlation function of a class of boundary operators could be written as an integral over a single copy of the black hole exterior, of an integrand that is a (multiple) discontinuity [9,10] (see also [21]). However, certain assumptions 3 were made in the process of deriving these statements. It transpires that for a class of operators, the bulk vertices for the corresponding holographic fields involves nonanalytic vertices. Having a contour prescription allows one to disambiguate these situations, and allows one to argue in its favour.
A simple example which explains the issue is a cubic interaction of two bulk-fields (say ϕ and χ) involving derivatives, e.g., χ ∇ A ϕ ∇ A ϕ. In this case, the interaction vertex written in a basis adapted to the analytic infalling modes has a pole at the location of the horizon. If we were to say that the prescription involves gluing together two copies of the black hole exterior along the codimension-1 null hypersurface (the horizon) we have a jump discontinuity in the (null) extrinsic curvature, which requires regulating. The contour integral provides a natural regulator, giving a clean prescription for the boundary correlator.
Another situation where one encounters something interesting corresponds to the computation of correlation functions of the boundary energy density operator. Now the integrands involve a function Λ k of geometric data and spatial momenta k, (2.17), whose zeros r k potentially pinch the integration contour for certain values of k ∈ C. 4 This behaviour turns out to be related to explorations in the literature regarding the contribution of the energy density operator to the out-of-time-order correlator capturing chaotic and scrambling dynamics [22] (see also [23] for an analysis in 2d CFTs). From a gravitational perspective the phenomenon of interest dubbed 'pole-skipping' in the aforementioned work, refers to loci in the complex frequency ω and momentum k space where the naive expectation of the correlator having a pole turns out to be unfounded. As observed originally in [24] and elaborated upon in [25], in the holographic context for neutral black holes, the phenomenon occurs at a codimension-2 locus in the (ω, k) space, at (ω ⋆ , k ⋆ ) = 2πi T (1, 2 (d−1) d ). The real space profile of the solution is that of an exponentially growing mode at the rate set by the maximal Lyapunov exponent λ L = 2πi T , damping out spatially with a rate set by the butterfly velocity v B = these loci the two-point correlator is ambiguous. From a mathematical perspective, at these particular values of the parameters the differential equations governing linearized fluctuations no longer have a regular singular point at the horizon. Rather the horizon becomes an apparent singularity. Physically, this means that the outgoing mode, which is generically non-analytic, 5 stops being so. Consequently, there is no branch cut in the grSK geometry and thus the retarded two-point function becomes ill-defined. Per se, this is not a problem, happening as it does in a set of measure zero. However, given that this happens when the function Λ k has a zero at the horizon (which a-priori could have led to an irregular singular point), one might ask if there are other regions where the poles of this function can interfere with the grSK contour.
We address this question and show that apart from the horizon, which is relevant for correlators of generic primaries, global currents, and the transverse tensor and vector polarizations of the stress tensor [30], the only other point where the function changes the character of the solution corresponds to the asymptotic boundary, which owes to the presence of soft modes (as was understood in [31]). In summary, we establish that for all bosonic wave equations of interest in the Schwarzschild-AdS d+1 background, the only singular points of the linearized equations are at the asymptopia, black hole singularity, and at the horizon (and other roots of the metric function f ).
While it is useful to establish that the solution to the wave equation corresponding to physical (gauge-invariant) gravitational fluctuations is well-behaved at other radial positions, it is not quite sufficient for the analysis of correlation functions. First note that the energy density operator couples to all the operators of the CFT; in the bulk this translates to the statement that the scalar graviton polarizations have a non-trivial vertex with all primaries (and contributes to ⟨T µν O O⟩). The corresponding bulk degree of freedom, Z, a scalar graviton polarization, satisfies an autonomous differential equation with only an apparent singular point at the roots of the function Λ k , and is therefore regular there. The metric functions themselves are determined not only in terms of Z and its derivatives, but crucially have a Λ k dependence that diverges as r − r k (the transformation between the metric perturbations and Z has a simple pole). Should the residue at this pole be non-vanishing for k lying along the ray described in footnote 4, we would have a problem satisfying the Schwinger-Keldysh collapse rules and the KMS condition.
Since this issue only arises upon complexification of momenta, a skeptical reader might ask whether this has any bearing on the physical correlation functions. While we are sympathetic to this view, it seems a bit bizarre that the real-time correlation functions given by the grSK prescription could fail to respect basic consistency conditions along a real codimension-1 surface in the complex (ω, k) space. Moreover, the universality of the energy density coupling, suggests that the resolution should be sufficiently general. Indeed, we verify by direct computation that not only is Z analytic at the zeros r k of Λ k , but so are the metric functions. Therefore, there are no subtleties even in the higher-point correlation functions computed using complexified grSK geometry. In other words, the analysis of these equations, provides a non-trivial consistency check of the grSK contour prescription.
We not only analyze the situation in the Schwarzschild-AdS d+1 black hole, which has been examined in part in the literature, but also undertake an analogous analysis for the Reissner-Nordström-AdS d+1 black hole. Now one has to contend with additional singular points arising from the charge dependence through the Ohmic radius introduced in [32]. These are benign lying as they do inside the horizon, and thus do not affect the grSK contour. For this charged plasma, however, the longitudinal scalar polarization of stress tensor couples to the longitudinal mode of the charge current. They can be decoupled at the linearized order, leading to an identification of the energy density and charge diffusion operators [33] (the corresponding two bulk fields are Z and V, respectively).
The decoupling again introduces an analogous function Λ k (still given by (2.17)). Since the metric function is different, the set of zeroes of Λ k is correspondingly distinct -it has two sets of zeros on two circles of radii r k1 and r k2 , with phases determined by k. Moreover, |r k1 | > r + for small charges, |r k2 | > r + for large charges, with the swap happening well away from extremality. 6 The mapping from Z and V to the physical metric and gauge field perturbations has poles at the zeros of Λ k . The energy density field Z has apparent singularities at both sets of roots, but the charge diffusion field V, while having an apparent singularity at the r k1 set, instead has a simple pole at the roots in the r k2 set. For either set of roots we can tune the momentum so that we have a potential pole along the grSK contour. However, once again, if we evaluate the metric functions in the vicinity of the Λ k zeros (either set), the metric and gauge field perturbations are manifestly analytic, ensuring that the correlation functions are well-behaved, and respect the Schwinger-Keldysh collapse rules and KMS conditions.
When we further fine-tune the momentum so that one of the roots from the r k1 set lies on the horizon, we find a behaviour analogous to the Schwarzschild-AdS d+1 case. The horizon can be made an apparent singular point at codimension-2 locus (ω ⋆ , k ⋆ ) = 2πi T (1, (d−1) r + 2π T ). 7 Translating to position domain we again have an exponentially growing mode in time, and damped in space, mimicking the scrambling mode for localized perturbations [34] (see the general expression given in [35]).
Above this critical charge, the aforementioned behaviour shifts to the charge diffusion mode. This however occurs at a different value of frequency, specifically at (ω ⋆ , k ⋆ ) = 2πi T (−1, (d−1) r + 2πT ). Now one no longer has an exponentially growing mode, even at an apparent singular point of the differential equation. Furthermore, the energy density field does not show any special features; its outgoing mode is non-analytic at the horizon as usual. This is perhaps unexpected from the viewpoint of [22], who argued for a general connec- 6 The parameterization for Reissner-Nordström-AdS d+1 we use has a charge parameter Q ∈ (0, d d−2 ] with the upper limit corresponding to extremality. The critical charge where the two sets of roots exchange dominance is lower; it is at Q * = d 3d−4 . 7 We work in units where ℓ AdS = 1, and thus the radial coordinate has units of energy. T /r+ is dimensionless, and we prefer to emphasize the dependence on the physical temperature measured in horizon size units. tion between quantum chaotic dynamics encoded in the exponential growth of the out-of-time ordered 4-point function to the energy density 2-point function. This connection fails to hold for charged black holes above a critical charge. While [36,37] have earlier analyzed the fluctuations of Reissner-Nordström-AdS 4 for connection between conserved current correlators and chaos they appear to have missed this fact. We return to what this could imply for maximally chaotic systems in §5.
The outline of the paper is as follows: In §2 we review the grSK geometry and the contour prescription for evaluating real-time thermal correlators. We also take the opportunity to give a synopsis of the features of second order linear differential equations that we are interested in, reminding the reader of some the standard terminology. This allows us to identify issues that arise in computation of real-time correlation functions using the grSK geometry. In particular, we highlight several potential failure points of the prescription, which we systematically explore in the sequel.
In §3 we summarize the upshot of our analysis of wave equations for the Green's functions necessary for the grSK contour. Throughout this work, we assume that our background geometry has time and spatial translational symmetries along the boundary (i.e., CFT) directions, which allows us to work in Fourier domain; our parameters also comprise then the frequencies and momenta. We discuss both the Schwarzschild-AdS d+1 black hole and the Reissner-Nordström-AdS d+1 black hole. Special focus will be on the fields dual to conserved currents, which have an intricate behaviour. Having established the analytic structure of the solutions of the perturbations, we then turn in § 4 to the computation of Witten diagrams using the grSK geometry. Here we discuss several scenarios: interaction terms where vertex functions have poles at the horizon, the apparent singularities in the energy density operator, and finally comment on the observations regarding 'pole-skipping' behaviour in Schwinger-Keldysh correlation functions. Some of these features are explicitly exemplified by a computation of a three-point function in a simple setting (derivative interaction of two scalar fields in the BTZ geometry).
To keep the main text simple, we have relegated various facts, and some detailed derivations to appendices. For completeness, we collate all the differential equations we analyze in Appendix A and some formulae relevant for the energy density field in Reissner-Nordström-AdS d+1 in Appendix B. A detailed analysis of singularity structure of the differential equations, which feeds into our summary in §3 can be found in Appendix C. For clarity, here we ascribe names to various singular points that occur, delineating features that are generic from those that aren't (the latter occur at higher codimension in parameter space). Appendix D contains some details relating to the evaluation of three-point functions in §4.4. Finally, Appendix E gives a quick and simple-minded discussion of deducing when a particular singular point of a second order ODE is apparent by a local Frobenius analysis (the general discussion dates back to [38]).

Preliminaries
Consider a general static, translational symmetric black hole geometry with the line element in ingoing coordinates given as (2.1) We restrict attention to non-degenerate horizons with f (r) having a simple zero. The largest real root f (r) is taken to be r = r + (the outer horizon) and one defines the Hawking temperature as

The grSK geometry
The grSK geometry (also referred to as the grSK contour) introduced in [7] is a particular complexification of the black hole spacetime (2.1), which is conjectured to be dual to the Schwinger-Keldysh contour of a thermal holographic CFT. In thermal field theory, one computes real-time observables by complexifying the time coordinate, to include the forward and backward evolution of the thermal density matrix ρ β → U (J R ) ρ β U † (J L ). Holographically, this is achieved by complexifying the radial direction in the ingoing chart specified above. Given the boundary Schwinger-Keldysh contour, we have a spacetime with two asymptotic boundaries, labeled L and R, above. We will work in the modified advanced retarded (FP) basis introduced in [39], where the sources J F and J P are the following linear combinations of the L and R Schwinger-Keldysh sources (likewise for operators) Here n B is the Bose-Einstein factor n B = 1 e βω −1 , which we note has poles at Matsubara frequencies w = i m, m ∈ Z. The advantage of this basis is that the Schwinger-Keldysh, and KMS conditions are cleanly implemented as the statement of vanishing of all correlators with purely F or P insertions. Any prescription for the bulk dual, therefore, must respect these constraints. Our goal here is to argue that the grSK geometry indeed achieves this in a self-consistent manner.
To describe the geometry, consider the mock tortoise coordinate defined as This coordinate ζ(r) has a logarithmic branch cut emanating from the zeros of f . We will be interested in the cut running from the outer horizon, which we ruin along the ray (r + , ∞). The normalization has been chosen so that ζ has unit monodromy across the cut. We take Re(ζ) = 0 on the top leg of the grSK contour and Re(ζ) = 1 on the bottom leg, as indicated in Fig. 1.
The basic data one needs to compute boundary Schwinger-Keldysh observables, is the boundary-bulk Green's function for the field dual to the QFT operator. Of primary interest Re(r) Figure 1: The complex r plane with the locations of the two regulated boundaries (with cut-off rc) and the outer horizon at r+ marked. The grSK contour is a codimension-1 surface in this plane (drawn at fixed v). The direction of the contour is taken to be counter-clockwise, encircling the branch point at the outer horizon. The cut is chosen to run out towards the boundary. Depending on the system we study, there will be other branch points, cuts, and singularities in the complex radial plane. We will discuss these in the sequel, in a case-by-case basis.
is the one regular at the future horizon, the ingoing Green's function G in (ζ, v, x). For planar black holes we exploit the Killing symmetries to work in Fourier domain (ω, k), with wavefunctions e −i ω v+i k·x and let k = |k|. We also work with a basis of the tangent space adapted to the ingoing coordinates [30]. In particular, the derivation is useful to analyze time-reversal properties. 8 Depending on the circumstance, it will also be useful to consider the frequencies measured in Matsubara units, so we introduce The ingoing Green's function G in (ζ, ω, k) satisfies suitable non-normalizable boundary conditions at the spacetime boundary, and is regular at the horizon. This is the analytic solution. From it one can recover the outgoing boundary-bulk Green's function by using the time-reversal property of the background [9]; specifically, G out (ζ, ω, k) = e −βωζ G in (ζ, −ω, k). Furthermore, the bulk-bulk Green's function may also be directly obtained from the knowledge of G in [10].
The computation of a thermal correlation function of the boundary field theory, therefore, involves convolving these boundary-bulk and bulk-bulk Green's functions, and carrying out the radial integrals. The main difference from the usual Witten diagram computation in Euclidean-AdS is that now the radial integral is carried out over the grSK geometry, using the contour depicted in Fig. 1. While the ingoing boundary-bulk Green's function being regular is insensitive to the branch cut, the outgoing boundary-bulk and bulk-bulk Green's functions pick up monodromies from the cut. The contour integral can, in the absence of any other singularities, be evaluated as an ordinary radial integral restricted to the domain [r + , ∞), with the integrand being given as a discontinuity of the convolved Green's functions, vertex factors, and radial Boltzmann weights (the e −βωζ factors) [10].

General properties of ODEs
As reviewed above, the essential ingredient in the computation of thermal real-time correlators is the ingoing boundary-bulk Green's function. These Green's functions in Fourier domain satisfy second order ordinary differential equations. The class of equations we are interested in are of the following form: The 'dilaton' factor e χ is a function of the radial coordinate. In certain cases it also depends on the spatial geometry, which will prove to be the ones of interest for us. The potential functions V (r) can be complicated, but these will be mostly irrelevant for our local analysis. We have compiled the various equations of interest in Appendix A.
Our aim is to analyze the singularity structure of the differential equations (2.7), which are second order, linear ODEs (SOLDE), in Fourier domain. These can be canonically be put in the form: We let r ∈ P 1 and wish to classify the branch points of the solutions φ(r) of (2.8). The interplay of this branch structure with the contour prescription of [7] described in §2.1 will be important for purposes of understanding higher-point thermal correlation functions. It may be possible to reduce the number of such branch points, by projecting down from the P 1 parameterized by r, say with the variable r 2 for even boundary dimensions. However, since we wish to work with the grSK geometry, which complexifies r, we will directly analyze the ODEs in this radial variable with no redefinitions. We recall some elementary facts. A point r = r 0 (including the point at infinity) is an ordinary point if p(r) and q(r) are analytic in an open region containing r 0 . If not, it can be a regular singular point if (r − r 0 ) p(r) and (r − r 0 ) 2 q(r) are analytic in an open domain of r 0 , or an irregular singular point otherwise. A SOLDE is Fuchsian if all of its singular points (abbreviated SPs) are regular (modulo a monodromy condition, which will describe below). For reasons that are not a-priori clear, all the wave equations we study for perturbations of non-extremal planar AdS d+1 black holes are Fuchsian with regular singular points. For equations of massive scalars in such backgrounds this was already noted in [40].
In the neighbourhood of a regular SP x = x 0 , the local behaviour of the solution can be understood by the Frobenius expansion. Laurent expanding the functions p(r) and q(r) about r 0 , and letting {p n , q n } be the coefficient of (r − r 0 ) n , we learn that characteristic exponents α ± are solutions to the indicial equation, These tell us that the solution φ(r) has a branch point at r 0 and can be used to set up a Frobenius expansion for φ(r) in an open neighbourhood. This local solution picks up a monodromy upon encircling r 0 , r → r 0 e 2πi , rotating by In the former case the two local Frobenius solutions φ ± are linearly independent, while in the latter case there is a logarithmic branch. 9 There is however one exceptional case: for α ± ∈ Z and α + ̸ = α − , it may be possible for there to be no monodromy -the logarithmic branch is absent. Furthermore, if both exponents are non-negative then the solution admits a regular Taylor expansion. In such a situation the putative SP r = r 0 is said to be an apparent singular point (abbreviated ASP), for it secretly is simply an ordinary point of the equation. An analysis of such ASPs, including a general criterion for their existence, can be found in [38], where they are referred to as pseudo-singular points in translation. For a Fuchsian SOLDE, the so-called Fuchs condition, requires that sum of the characteristic exponents from all the SPs, regular and apparent, (including infinity, whose exponents counted as usual with an opposite sign) is fixed in terms of the number of SP: if we have an equation with m SPs, then m i=1 α i+ + α i− = m − 2. A classic reference on this material is [41]; for modern perspective on Fuchsian ODEs see [42].
The singularity analysis of equations compiled in Appendix A, which are of the form (2.7), is straightforward. Our interest will be in ascertaining how the singular points morph as we change the parameters in the equation. Before we do so, however, let us take stock of why this is relevant for the interaction vertices. Following this, we will summarize the basic features of the equations relevant for computing conserved current correlators in Appendix C.

Bulk interactions
We will focus on two types of bulk interactions, which require attention for checking the consistency of the grSK contour prescription. The first of these, involves interactions comprising of derivatives of fields. For simplicity, consider the following cubic interaction of two scalar fields ϕ and χ While it looks innocuous, recall that the choice made for the grSK contour requires the ingoing solution to be analytic, putting all the monodromies into the outgoing solution. Adapting to the ingoing tangent basis, the vertex can be written as 10 Assuming α+ > α−, the linearly independent solutions are φ+ and φ− + γ φ+ log(r − r0). 10 The integral with subscript k is a shorthand for the momentum space integrals with the usual momentum conserving δ-functions. The d-vector k µ has components k µ = (ω, k), with k ≡ |k|. The plane wave basis in R d−1,1 is taken to be e −iωv+ik·x . We have also made clear that the radial integral is to be evaluated on the grSK contour Fig. 1. where we introduced D ± = r 2 f ∂ v ± iω, and have restricted attention to translationally invariant backgrounds. The radial kinetic term here has now a simple pole at the horizon, which needs to be dealt with. This behaviour is generic: for derivative interactions one always encounters a pole at the horizon, the order being related to the number of derivatives involved. This phenomenon was in fact noticed before in the analysis of fluctuations of a probe Nambu-Goto string [8], where the quartic coupling of the transverse string modes have a double pole at the horizon.
The second vertex of interest is the interaction between a primary O of dimension ∆, and the energy density operator. The cubic vertex arises already from the kinetic terms for the scalar field For the scalar graviton polarizations, which capture the dual of the energy density operator, we can express the metric perturbation in a suitable gauge as [31]: (2.14) From here it is straightforward to read off the couplings (2.15) The complication lies in the translation to the fields which satisfy simple ODEs since, owing to the background gauge invariance, the fields Φ E , Φ O , Φ W are not independent degrees of freedom. As explained in [31,33], the information about the perturbation can be repackaged into gauge invariant variables (following the original analysis of [43,44]). For a planar Schwarzschild-AdS d+1 black hole, we have a single energy density operator, dual to a field Z, which determines the metric functions as [31]: The function Λ k which will play an important role below is defined as The corresponding expressions for the planar Reissner-Nordström-AdS d+1 black hole are given in Appendix B, since in this case we have to also deal with the mixing between the energy density and charge diffusion modes. Despite the additional complications, one finds the translation to the gauge invariant variables to involve the function Λ k as in (2.17), with the only change being in the metric function f (r). For the neutral black hole, we can therefore write the vertex as The point to notice here is the presence of the function Λ k in the denominator, which potentially can give rise to singular contributions. We are interested in ascertaining whether its zeros are relevant, and if so whether one also has to account for the monodromy of the field Z about these singularities. This requires us to undertake a careful analysis of the differential equations, thus motivating the study in this paper. We note here that the factors of Λ k in the vertex are due to the redefinition of the metric variables in terms of Z in (2.16). Independent of applications to real-time correlators, one would also like to know if metric perturbations are regular themselves. On physical grounds we expect them to be so for real frequencies and momenta, but the complete analysis, as it will transpire, should also allow for complexified spatial momenta.

Green's functions on black hole backgrounds
We turn to an analysis of the wave equations encountered in the study of black hole perturbations. We will focus on the set of equations that captures holographically the dynamics of spinless operators of arbitrary conformal dimension and conserved currents, for the most part. For systems with multiple conserved currents (as in the case of finite density), diagonalize the kinetic terms for the corresponding bulk fields. In all cases of interest, the dynamics of conserved currents can be repackaged into equations for designer scalars with non-minimal couplings. The details of the analysis (which is pretty standard) are compiled in Appendix C. We simply summarize below the features which play a role in the computation of Witten diagrams using the grSK contour prescription.

Schwarzschild-AdS black hole
Let us begin with a class of equations that encompasses minimally coupled massive scalar fields dual to CFT primaries, probe gauge fields dual to a global current operator, and tensor Figure 2: Singular points of Schwarzschild-AdS5 designer scalar equation in the complex r. We have drawn a circle of radius r+ that separates the interior of the black hole from the exterior. The cut depicted is the one running from the horizon to the AdS5 boundary for the outgoing Hawking mode (we make the choice that the ingoing mode is analytic). The grSK contour encircles this cut counter-clockwise as shown. The cuts from other singular points are not indicated, but can run in any directions so long as they do not cross the grSK contour.
and vector polarizations of gravitons. As demonstrated in [30] all these cases are described by the dynamics of a designer scalar field (2.7) with a simple power law behaviour of the auxiliary dilaton, This system has SPs at the asymptotic boundary, at the black hole curvature singularity (the origin), and at the zeros of f for generic values of ω, k. We have depicted this generic behaviour in Fig. 2 The behaviour at infinity dictates the familiar asymptotic fall-offs. The only other SP relevant for the grSK contour, is that originating from the outer horizon r = r + , where the characteristic exponents are 0 and i βω 2π , corresponding to the ingoing and outgoing solution, respectively. Using the freedom to orient cuts from the other SPs away from the grSK contour, we see that there is no ambiguity in the prescription. Thus, in this case we only need to deal with issues arising from bulk vertices, which as noted in § 2.3 may themselves contribute singularities.
There are additional features at special values of frequency and momenta. At Matsubara frequencies w = i n with n ∈ Z ≥0 the outgoing solution can be made regular at the horizon by fine-tuning the momentum. For this choice the horizon is an ASP. While we explore these in Appendix C.2, in the context of grSK contour, their presence simply implies that the correlators are analytic at these special points.
There is one equation (A.3), that for the bulk mode Z(r) dual to the energy density operator, which has to analyzed independently. Naively, it appears that in addition to the asymptotic, curvature, and horizon SP, we have an additional d − 2 SPs from the roots of Λ k = 0. For Schwarzschild-AdS d+1 are located on circle of size 1/k at where ϖ n denotes the n th roots of unity. For arg(q 2 ) = π we even potentially have one of the roots located between the horizon and the asymptotic boundary. However, a closer analysis reveals that for generic q 2 the points r k are apparent SPs (ASPs) of (A.3).
At special values of parameters, we find new behaviour. The roots of Λ k can be made to coalesce with other SPs of the system, by tuning the momentum. This leaves the Fuchsian nature of the equation unchanged, but modifies the characteristic exponents. There are two cases of interest: r k → ∞ and r k → r + .
When k = 0, the roots of Λ k are at infinity. However, this leaves the point at infinity a regular SP, albeit with changed exponents. As was understood in [31] the change in behaviour owes to enhanced diffeomorphism symmetry in the zero momentum limit (a fact already appreciated in [43] and explained in part in [45]). 11 By choosing we can achieve r k = r + . For generic ω, the only change is that the characteristic exponents are shifted by unity: the ingoing mode has exponent +1, and the outgoing mode has exponent 1 + i βω 2π . At a particular fine-tuned values of (ω, k), setting ω = 2πi T the energy horizon coalescence renders the horizon to be an ASP. At this point profile of Z along the boundary directions takes the form: .
The result is suggestive of an exponentially growing mode in time at a rate set by the maximal Lyapunov exponent λ L = 2πT with the spatial dependence being damped at a rate given by the butterfly velocity. These are in fact, the values computed for Schwarzschild-AdS d+1 black holes using the shockwave approximation [26] (see also [34]). As noted earlier, this phenomenon was first noticed in [24] and explained carefully in [25]. In the SYK model it is the presence of such an exponentially growing mode that leads to the growth of the out-oftime-order correlation functions [46]. Inspired by this, an effective theory of maximal chaos was proposed in [22], arguing for there to be a growing mode in the energy density correlator. We have simply verified these earlier observations (trivially extending them for all d ≥ 3). At this particular value of parameters, as with the Matsubara ASP, there is an ambiguity in the definition of the boundary retarded Green's function. In the present case this is just as well, since an exponentially growing mode would have signaled an instability (the quasinormal spectrum is required to be supported in the lower half-plane with retarded Green's functions analytic above the real frequency axis). We are also wary of referring to this phenomenon as 'pole-skipping' for similar reasons, and prefer the term, apparent anti-quasinormal mode.
In summary, the singularity structure for all the ODEs of interest in Schwarzschild-AdS d+1 is captured by the designer scalar equation (3.1), as depicted in Fig. 2 for d = 4. The regular SPs are only at the roots of f , the boundary and at r = 0. The nature of these regular singular points changes due to coalescence as we tune k to some special values, but apart from that there is no other singularity to deal with.

Reissner-Nordström-AdS black hole
For the Reissner-Nordström-AdS d+1 black hole, the analysis of the differential equations is quite similar. The main differences are that depending on the equation we study there are, in addition to the SPs at infinity, the curvature singularity, and the roots of f (r), additional features at specific radial positions. One of these corresponds to the roots of f ′ (r), denoted r Q , defined by The dimensionless parameter S Q is serves as a proxy for the charge, with S Q = 0 in the neutral case, and S Q = 1 for the extremal black hole. The other is the roots of the function Λ k , which are distributed in two sets on circles of radii r k1 and r k2 , respectively. Of these r k1 is analogous to the locus r k for the Schwarzschild-AdS d+1 black hole; it is an apparent singular point. The locus r k2 , however, is a genuine SP for the charge diffusion mode. There are of course some special points in the parameter space where the nature of the SPs changes (eg., at Matsubara frequencies for fine-tuned momenta at the horizon SP). Details of these can be found in Appendix C.3, but the general picture is summarized in Fig. 3. There is one particular situation of fine-tuning that we wish to emphasize here, involving the roots of Λ k lying at the horizon, and leading to apparent SPs. The perturbations involved are the field Z dual to the energy density operator, and V dual to the charge diffusion operator, which are analyzed in Appendix C.3.3. These two satisfy decoupled wave equations with the function Λ k controlling both their kinetic terms. The roots of Λ k satisfy r k2 < r + < r k1 for S Q ∈ [0, 1 2 ), but are ordered as r k1 < r + < r k2 for S Q ∈ [ 1 2 , 1). When the dimensionless ratio S Q < 1 2 , and k chosen so that r k1 = r + , we find that the horizon is an apparent singularity for the energy density field Z, at a specific frequency ω = 2πi T . However, there is no choice of frequencies for which the field V is analytic, and so the horizon remains a regular SP for this field.
At this value of parameters, the profile of Z is as in (3.4), viz., an exponentially growing mode. We recover from this profile the Lyapunov exponent for a neutral probe. 12 Further- 12 While it has been argued that in the presence of global conservation laws the bound on the Lyapunov exponent is modified to 2πT 1− µ µc , [47], our probes are the conserved currents which are uncharged under the symmetry. Figure 3: Singular points encountered on complex radial plane for the various differential equations of interest in the Reissner-Nordström-AdS5 background. The singular points {±r+, ±r−, ±i ri are the roots of f (r), which along with the point at infinity are SPs of all the equations of interest. The black hole singularity is not a SP in all cases, but is for some of the equations. Additionally, we have some equations (transverse vector perturbations of the gauge field) for which the roots of f ′ (r) located at ±r Q are also singular. Finally, r k2 andr k2 are the roots of Λ k which are regular SPs of the field dual to the charge diffusion operator. more, from the value of the momentum at which the outgoing mode remains analytic, we recover the butterfly velocity. We conclude that by fine-tuning parameters, the energy density field can have a putative exponential growing mode, which is analytic at the horizon, with the wavefunction in the boundary directions behaving as This is indeed the expected value of the butterfly velocity for a AdS black hole [35]. For this range the charge density field has a non-analytic outgoing mode at the horizon. For S Q > 1 2 , however, something very different occurs. In this range of charges we can tune r k2 = r + by choosing the momentum appropriately. We find that energy density mode Z is always non-analytic at the horizon, for there is no choice of frequency that makes the outgoing mode regular. However, the horizon can be made into an apparent singularity for the charge diffusion mode, albeit at a different frequency, ω = −2πi T . This frequency leads to an exponentially damped mode, that cannot be identified as having to do with any kind of Lyapunov behaviour (and thus with the chaotic behaviour of the OTOC).
This change in behaviour of the fields the horizon occurs for S Q = 1 2 , which translates to a Reissner-Nordström-AdS d+1 black hole having a charge parameter The parameterization of the metric in (C.13) with 0 ≤ Q ≤ d d−2 , the upper limit corresponding to the extremal solution. When the charge of the black hole is small Q < Q * , the behaviour we find for the energy density mode Z is analogous to that for a neutral black hole, by continuity. However, there is a transition at a critical charge Q * (because there is a secondary branch of zeros). For large charges there is no exponentially growing analytic mode, in either the charge density, or in the energy density operators. While in neither case is one actually computing the OTOC directly, the essential point is that in the absence of a temporally growing mode. Therefore, there can be no meaningful association with chaos for large enough charge. As far as we can tell this point has been missed in earlier analyses [36,37]. We will examine what these features mean in the Schwinger-Keldysh correlator, in due course.

Bulk interactions and the grSK contour
Having ascertained the behaviour of solutions to the wave equations, we can now turn to the analysis of boundary correlation functions. As described in §2.3, our primary concern is with bulk vertices that potentially have singularities that interfere or pinch the grSK contour. The locations of the SPs along with the grSK contour for the equations of interest are given in Figs. 2 and 3 (for specific equations in the charged case, see Figs. 4 to 6).

Interaction terms with horizon poles
Let us being with the cubic interaction of derivative coupled scalars (2.12) where the vertex function naively has a simple pole at the horizon. To process this vertex we will need to understand which correlation function is affected by the presence of the pole. In particular, we first need to check that the SK and KMS conditions are not spoiled by the presence of the horizon pole in the vertex. Furthermore, we need to determine which of the non-vanishing correlators get a contribution from it and estimate them.
To proceed, we then recall some useful facts about the D ± derivatives that was used to relate the outgoing boundary-bulk Green's function to the ingoing one [9]. In terms of the ingoing adapted derivatives satisfy a conjugation property Using this we can argue that, given an ingoing boundary-bulk Green's function G in (ζ, k µ ), the outgoing Green's function is given by Herek µ is the frequency reversed d-momentum; if k µ = (ω, k), thenk µ = (−ω, k). In terms of this data the solution to the free wave equation for a field ϕ on the grSK geometry with sources J R and J L on the R and L boundaries is given in the FP basis by (cf., [9]) Next, we note that suitable combinations of the D + derivative and temporal derivative lead to factors of f , viz., Hence, the only case of concern is one where there is a non-trivial term without a factor of f , which occurs when we have a contribution of the form With this information, we now go over the various 3-point functions and check what the presence of this pole does. For the vertex arising from (2.12) we will use G ϕ in (ζ, k µ ) and G χ in (ζ, k µ ) to denote the Green's functions for the two fields. We also order the fields with momentum labels k µ 1 and k µ 2 assigned to the two ϕ fields, and k µ 3 assigned to χ. The arguments below are a specialization of those given in [9] (see also [10]) so we will be brief in our exposition.
• The horizon pole is benign for the FFF correlator which vanishes as desired. Here factor gives a factor of r 2 f which cancels the pole. The vanishing of the correlator follows along the lines then described in the aforementioned references. A similar argument applies to the FFP correlator.
• Using a similar property of (D + − iω)G ϕ out (ζ, k 2 ) as noted above, we can check that PPP, PPF and PFF correlators (nb: last label is for χ) are unaffected. In particular, PPP correlator vanishes, and the other two are computed as the integral of the discontinuity across the cut.
• This leaves us with FPF and FPP, which are both non-trivial since there is no factor of f from the derivatives acting on the Green's function for ϕ. The FPF correlator turns out to be computed by the integral where we have made use of the relations deduced above. Now we have an explicit horizon pole in the integrand. A similar story pertains for FPP, which requires us to look at the following integral: In both cases the first lines of the expressions have a pole at the branch point, which needs evaluation.
To get a sense of how to tackle these integrals with a pole on the branch point, let us consider the following toy integral (here C H is the keyhole or Hankel contour as in Fig. 2) The instinct is to realize that brackets in the last line involve two divergent pieces, which should cancel amongst themselves. The result is then which is manifestly analytic as a function of α.
The integral we want with the horizon pole is something more involved. Moreover, the Green's functions do not have exponential damping but only power-law behaviour in the radial direction. However, we can use the above example to inspire a simple prescription for the computation. Let us therefore consider the following integral which captures both correlators of interest (with x ∈ {F, P}): We have isolated the dependence on the frequency that enters the monodromy and suppressed dependence on all other parameters. The function in the integrand, X, is obtained from the various boundary-bulk propagators, measures etc. The only information we need from it, is that it is analytic at r = r + . By the above logic The final result involves a localized contribution from the horizon the serves to cancel the IR divergent term in the radial integral, giving thereby a sensible result for the boundary correlation function. In [8], such a localized contribution was extracted by working in a small frequency gradient expansion.
From our perspective, the contour integral prescription of [7] has the advantage of making clear how to deal with these horizon poles. In the earlier prescriptions of [14], as for example recently employed in [20], the piecewise integrals will each have IR divergences, which will need regulating. The contour prescription provides one such, and has the added benefit of being geometric. We will employ this to compute an explicit correlator in §4.4 below.

Coupling to energy density
The other coupling of interest is the coupling to the energy density operator. We have explicitly written out the cubic coupling of a primary of dimension ∆ to the energy density operator (2.18). There are several terms in S EOO , but all have a common factor of Λ −1 k . There are additionally potential horizon poles. We expect something similar for cubic self-coupling of the scalar graviton polarization with itself. For instance, there will be a Z 3 vertex (with some complicated action of derivatives), which will contain an overall Λ −3 k . While computing the contributions to the ⟨T µν OO⟩ three-point function, naively, one will have to evaluate integrals like˛d The horizon pole contributions can be understood as in §4.1 (later we will argue that some factors present in (2.18) are spurious). Let us therefore focus on the Λ k contribution for now. If the zero of Λ k , r = r k lies outside the contour of integration, then we do not need to concern ourselves, and the calculation proceeds as usual. The key issue is when the grSK contour is pinched by the pole at r k , which, as we have seen, can occur for complex momenta. For instance, restricting attention to the Schwarzschild-AdS d+1 case, taking arg(q 2 ) = π will achieve this. 13 It is also worth emphasizing that if we focus on small k and analyze the problem in a boundary gradient expansion, as has been done in earlier works [31,33], then the zeros of Λ k do not show up.
The ingoing Green's functions for the fields dual to the energy density operator (or charge diffusion) do not have a branch cut at this point from the SOLDE analysis. However, the explicit pole in the integrand (4.12) cannot be avoided. This is a serious problem for the grSK contour, unless there is a mitigating effect from the rest of the integrand. This is because, should the function X(ω 1 , ζ) have a constant limit at r k , then we will pick up this value, as a localized residue contribution at the pole. Such a contribution will violate the SK and KMS conditions, rendering the grSK prescription suspect. The easiest way to see this is to note that one can end with a non-vanishing answer for the correlator for all F-type operators with energy density insertion. While the usual Schwinger-Keldysh rules will imply vanishing in the boundary, X(ω 1 , ζ) will be a product of ingoing Green's function, which are analytic, and may have a non-vanishing residue (setting ω 2 = 0 for simplicity). The only way to avoid this would be for the residue to vanish.
In other words, consistency with the Schwinger-Keldysh rules, requires that the all interactions with the energy density operator, have suitable zeros in the numerator to cancel the inverse powers of Λ k . The origin of the Λ −1 k factors in the integrand is from the field redefinitions (2.16) between the metric fields {Φ E , Φ O , Φ W } and Z. Explicitly, carrying out a local Frobenius expansion of Z in the neighbourhood of r k , we can that the integrand has no poles for any interaction involving the energy-momentum tensor's scalar polarization. 13 For Reissner-Nordström-AdS d+1 we can attain this by requiring arg(p 2 s ) = π S Q < 1 2 or arg(p 2 s ) = 0 when S Q > 1 2 , respectively (cf., Appendix C.3.3).
As such, one would find it strange for the metric functions to diverge in the physical region between the horizon and the boundary for complex momenta. Lest the reader imagine that this would be guaranteed, we should clarify that the intuition which holds for real momenta could have failed for complex momenta. That this does not happen involves nontrivial cancellations as can be checked algebraically. Once this is established, one can directly use (2.15) to deduce that the correlation functions remain well-defined even for arg(q 2 ) = π in the Schwarzschild-AdS d+1 geometry.
In fact, a similar reasoning can be applied to the potential horizon poles in (2.18). Reverting to the metric perturbation itself, say (2.14), which has explicit factor of f in the denominator, we recall, that the natural cotangent basis in ingoing coordinates involves the one-form dr r 2 f . One can see that all the factors of 1/f are part of such a one-form, which secretly guarantees regularity of the metric perturbation at the horizon. However, in the kinetic term of the field ϕ dual to the primary O the factors of 1/f are physical, and do contribute as in (4.11).
For Reissner-Nordström-AdS d+1 we have to deal with an additional complication. The field V actually diverges at r = r k2 . This seems more serious; (B.3) and (B.4) suggest that there might be double-poles in the metric and gauge field perturbations. However, once again, the physical metric and gauge field perturbations are regular at the zeros of Λ k . This can be directly checked using the local Frobenius solution of {Z, V}.
The upshot of our analysis is that the zero loci of the function Λ k , which was introduced to find autonomous SOLDEs for the perturbations, are irrelevant for the purposes of computing any physical observable. These zeros are for the most part ASPs of the differential equation (the exception being the one set of zeros in the charged black hole for V), and not only are the decoupled fields regular there, but the physical metric and gauge field functions are themselves regular. The skeptical reader might inquire whether this was all a red herring from the beginning. Could one have dispensed with fields like Z in the Schwarzschild-AdS (or Z and V in the Reissner-Nordström-AdS case) and worked with the metric functions directly? The function Λ k was originally presented in [43,44], and a (technical) rationale was its origin was offered in [31,33]. What we can assert with some confidence is that none of the metric functions themselves satisfy an autonomous SOLDE; one just has a coupled system of first and second order equations (see the latter works for a thorough analysis of the dynamics). While there is some satisfaction in noting that there are no special loci in the geometry determined for a subset of momenta, we do not have at present more physical picture to paint for the existence of Λ k and its irrelevance in the final answers.

Apparent singularities at the horizon
The preceding analysis focused on having additional poles in the integrand either at the horizon, which generically is a regular SP, or at r k which is always an ASP of the SOLDEs. Let us now turn to the case where the horizon itself becomes an ASP, at special codimension-2 loci in the complexified frequency domain. There are two situations we should consider here: the Matsubara ASP, which pertains for all the equations analyzed herein, and the case where the energy ASP at r = r k merges with the horizon leading to energy horizon coalescence described in Appendix C.2.2 and Appendix C.3.3. In these cases, particular values of frequency and momenta render the horizon an ASP of the corresponding equation. Even in the exceptional case of charge diffusion mode V, at the coalescence of a root from the set r k2 with r + , we encounter an ASP along a codimension-2 locus.
Consider first the Matsubara ASP, where for ω = −2πi m T , with m ∈ Z + and fine-tuned values of q ensuring that both the ingoing and outgoing modes of some field are analytic. Let us first compute a 2-point correlation function. Such a correlator is computed by a pure boundary quantity, and as noted earlier has a one-parameter ambiguity. The fine-tuning involved to make the horizon an ASP is codimension-2 in the (ω, k) space. Generally, the 2-point correlators are meromorphic with simple poles along the codimension-1 quasinormal spectral curve ω = ω QN (k, m) labeled by a discrete parameter m. To determine the behaviour of the boundary Green's function at isolated points associated with our fine-tuning, we need to specify a direction of approach, which gives us our one-parameter ambiguity. Curiously, at these apparent quasinormal modes, examining the 2-point function, one might conclude that the black hole does not Hawking radiate and appears to lose its thermal character. Now consider a higher-point function arising from a bulk contact interaction without a horizon pole. Generically, the integrands have a branch-cut, coming from the outgoing Green's function, the coefficient of the advanced sources J P . So the integrand for an n-point of the form F . . . FP . . . P with (n−l) F and l P fields, can be arranged to not have a branch cut if all the P fields are at a Matsubara ASP point in parameter space. Note that with sufficient retarded (F) fields we do not have a constraint from momentum or frequency conservation. Owing to the absence of a cut in the integrand, such correlation functions vanish. This is not surprising per se; as already noted in [9] and explored further in [10], even in the generic case correlation functions can vanish at Matsubara frequencies. What is special here is that this happens only when the momenta are also fine-tuned. The one exception is when the vertex has a horizon pole as in (2.12). Then the contour picks up the residue at the pole.
Thus, unlike the 2-point functions, we have just argued that the effect of the Matsubara ASPs on higher-point functions is more benign. At most the correlation function with some advanced operators tuned to special points in the parameter space vanish.
We can argue likewise for the energy horizon ASP. Now the energy density operator, say Z in a neutral plasma, is analytic for a particular frequency ω ⋆ = 2πiT and momentum k ⋆ . So correlation functions with a number of advanced operators dual to Z will have no branch cut at (ω ⋆ , k ⋆ ). Correlators of such operators along with some other advanced operators tuned to the Matsubara ASP, and potentially any number of retarded operators, will produce an integrand with no branch cut emanating from the horizon. Such correlators will vanish, unless there is a horizon pole, which allows one to pick up a localized residue contribution. Note that argument does not affect the retarded field. Moreover, it also doesn't guarantee an extended domain of analyticity for the corresponding advanced field, since we only control the behaviour at an isolated point. The advanced correlator of Z(ω ⋆ , k ⋆ ) continues to be analytic in the lower-half plane, and potentially free of singularities at this codimension-2 locus.
The fact that the correlation function of advanced energy density operators at (ω ⋆ , k ⋆ ) gets only a localized contribution is reminiscent of the argument of [48]. They argued for localized contributions for the out-of-time-ordered correlator at the horizon due to high relative boost (in the shockwave approximation). While we are only computing the Schwinger-Keldysh ordered correlators, the assertion is that exactly at the particular locus where the spatiotemporal dependence of the advanced field takes the form of an exponentially growing mode, we pick up a contribution from the horizon. As noted before, this motivated a proposal for an effective field theory for computing out-of-time-order correlators in maximally chaotic systems presented by [22]. They however required that the growing mode was part of the retarded field, but as we have seen, there is nothing special for the retarded field at (ω ⋆ , k ⋆ ). The special aspects of this locus pertain to the advanced operator, but all they do is force the correlator to vanish for isolated points. Our analysis, however, does not indicate whether, there is a direct feature associated with these special points in parameter space for out-of-time-order correlators.

Three-point correlator with horizon pole contribution
To illustrate the general considerations of the preceding discussion, we will examine a threepoint correlator arising from derivative interactions in the bulk. Our aim is to explicitly capture the contributions from the horizon pole of the vertex function, as described in §4.1. To keep things simple, we will evaluate I FPF for the interaction (2.11) in the BTZ geometry. This is a simple setting, where 3-point functions for non-derivative interactions can be obtained analytically [9] (four-point functions were computed in [10]).
Let us start with the simplified expression for the FPF correlator (4.6). Expanding out the derivatives D ± in (4.6), we have (4.13) The last term in the above expression has an explicit pole at the horizon, while all other terms are manifestly regular (since the ingoing propagators are analytic). Naively, one might expect this to lead to an IR divergence, but as we have argued, based on a local analysis, this is not the case. To evaluate the integral we employ a simple change of variables, r = r + /z, and rewrite the above as (4.14) The regular piece, I reg FPF , has been reduced to an integral over a single sheet, by picking up the discontinuity across the horizon branch cut. This term is evaluated exactly as in [9].
The ingoing bulk-to-boundary propagators have analytic expressions in the BTZ geometry. For a minimally coupled scalar primary of conformal dimension ∆, one has c ; z is the regularized hypergeometric function, and Furthermore, in the present case, we have To evaluate the integrals, we let ∆ 1 = ∆ 2 = ∆ ϕ and ∆ 3 = ∆ χ , for simplicity, Further introducing, as in [10], the function in terms of which, the boundary retarded Green's function of a primary with dimension ∆ is given by .
The rest of the evaluation proceeds by employing the Mellin-Barnes representation of the hypergeometric function. For the part of the integrand which is regular on the horizon, we find We have written the expression by factoring out the boundary two-point function for the fields, which makes manifest, as in [9,10] the analytic structure of the correlators, for the form factor J reg is an analytic function in the Fourier domain. Since the expression for this function is a bit cumbersome, it is relegated to Appendix D. The localized contribution, which is new, instead gives upon using the Mellin-Barnes representation, the following: (4.21) The key point to note is that the IR divergence from the radial integral, precisely cancels against the localized contribution (the last term in the final line). This is indeed in accord from the local analysis of §4.1. Evaluating the integral, we find at the end of the day, The form factor J loc , which is also given in Appendix D, is manifestly regular. Thus, all the singularities of the three-point function are isolated in the pre-factors, which are products of two-point functions. Combing the two pieces, we find The form factor J can be found in (D.6). This result is similar to that obtained in [9,10], where there were no vertices with poles at the horizon. In particular, the analytic structure of the three-point function is unchanged from the term containing the horizon pole. There are quasinormal modes for the two F operators in the lower half of the ω 1 and ω 3 plane, and anti-quasinormal modes in the upper half ω 3 plane. This piece has been isolated in the first line (which can be rewritten in terms of the product of the two-point function of the three fields). There are some zeros of the correlator at Matsubara frequencies, from the Boltzmann factor 1 − e 2πw 2 (though half of these cancel the poles of Γ(1 + i w 2 )). The reader can also verify that there are no features due to the apparent singularities at the horizon, as argued in §4.3.
While this example illustrates how the grSK contour deals with singularities at the horizon, writing down an example with the Λ k terms in the vertex is a bit more involved. We have however identified a situation, involving the Chern-Simons coupling of a U (1) gauge field in Einstein-Maxwell-Chern-Simons theory, where the computation of anomaly induced terms in the thermal correlator requires understanding both the horizon poles, discussed herein, and the potential singularities at the zeros of Λ k . This analysis will appear in a separate publication.

Discussion
Our goal in this paper was to provide further evidence that real-time thermal correlation functions in holographic field theories are sensibly computed on the grSK geometry introduced in [7]. While it has been clear from earlier works [8][9][10] that the prescription determines the correlators consistent with Schwinger-Keldysh and KMS conditions, there were some specific features associated with conserved currents that had not been addressed hitherto. Most of the features we discussed are not always directly visible in a low-energy gradient expansion, which owing to the absence of closed form solutions, we generically have to resort to for analytic expressions.
A related motivation for us was to close the gaps in our recent discussion regarding the analytic structure of holographic thermal correlators [10]. There we had assumed that the bulk vertices did not have any singularities which interfered with the grSK contour. As we have seen here, there are indeed such vertices, either when the bulk fields have derivative interactions, or when we consider the gauge invariant combination that corresponds to the energy density field (and charge diffusion field at finite density).
One could have argued that such singularities, especially the latter, might have originated from our poor choice of variables. However, for reasons described in [30] (and amplified in [31][32][33] for charged plasmas) the bulk fields we have worked with are natural. They correspond to gauge invariant data, have good time-reversal transformation properties, all of which allows for determination of bulk Green's functions without ambiguities. The price one pays is the relative complexity of the differential equation obeyed by the fields. Therefore, the first task we undertook here was to clarify the nature of the differential equations involved and elucidating their general structure.
We find it curious that all the differential equations of import in static, translationally invariant AdS black hole backgrounds analyzed here are Fuchsian. This has been known for the massive scalar wave equations [40] (and in fact known already in [3] for massless scalars), and we verify it to be so even for conserved currents. While we don't a-priori see a rationale for this structure, as argued in the preceding reference one could relate the connection problem for finding quasinormal modes to instanton partition functions in supersymmetric field theory. This has indeed been successfully used to extract thermal correlators of scalar primary operators of arbitrary dimension in [49] in 4 dimensional holographic CFTs (and generalized to include charge and rotation in [50]). These analyses work with a different variable z(r 2 ) that reduces the number of regular singular points to four, allowing for use of known data for Heun equations. Since we wished to use the grSK contour in the radial variable, we have refrained here from adopting such a change of variables in our analysis. It is clear that these techniques can be adapted to deduce current correlators as well.
One corollary of our analysis was a refined understanding of the purported connection between energy transport and maximal scrambling behaviour of black holes. The effective field theory description of maximally chaotic systems studied in [22,23] has been justified in holographic context by invoking the phenomena of 'pole-skipping'. As has been argued by others, this phenomenon is not just restricted to the energy density operator, but acquires special import in this case, since the wave function in the boundary direction resembles a scrambling mode: a temporally growing mode, which is damped out spatially at a rate fixed by the butterfly velocity. We prefer to view this phenomenon in terms of apparent quasinormal modes (for reasons already explained in [10]). In particular, the mode in question arises when the zeros of a particular function Λ k , that naturally appears in the bulk ODE, coincides with the location of the horizon.
While for a neutral black hole the apparent quasinormal mode behaves like a scrambling mode, it fails to do so for a charged black hole past a critical value of charge. The specific value of charge does not appear to be singled out from the geometry in any way apart from being well away from extremality. Past this critical charge there is no mode that has exponential growth in time in the energy-momentum tensor components. This suggests to us that the link between such apparent quasinormal modes of the energy density and chaotic dynamics is perhaps coincidental.
Finally, our analysis also revealed one further surprise for charged black holes. While the black hole singularity is a regular singular point for all the fields analyzed herein in a neutral black hole background, it generically isn't for the charged black hole. Per se, this is not surprising, since the causal structure of the two spacetimes are different. However, it is curious that the momentum diffusion and the transverse photon fields around a Reissner-Nordström-AdS black hole actually are sensitive to the black hole singularity. This deserves further investigation, to understand whether such operators can be used to extract information about the spacetime in the vicinity of the singularity along the lines investigated in [51][52][53].
speaking the equations of interest reduce to those of a massless scalar equation with nontrivial kinetic terms, and a radially varying potential, schematically of the form given in (2.7). We will present the equations as they were introduced in the analysis of thermal correlators in [30,31] for the Schwarzschild-AdS d+1 perturbations, and in [32,33] for the Reissner-Nordström-AdS d+1 case. We will write them in a time-reversal invariant form, using the ingoing adapted radial derivative

A.1 Perturbations of Schwarzschild-AdS d+1
The class of equations we need to analyze for the perturbations of planar Schwarzschild-AdS d+1 with the geometric data given in (C.1) are as follows: Designer scalar equation: Massive designer scalar with index M, which is captured by the ODE The interesting cases here are M = d − 1 and m 2 ̸ = 0 for a regular massive scalar, dual to a conformal primary of weight ∆ of the boundary CFT. A special case is the massless scalar, which is also the equation for transverse tensor gravitons (which comprises of  [30].

Scalar graviton perturbations:
The one other equation of interest is the scalar graviton polarization, which characterizes the dynamics of the energy density field, and takes the form [31] with Λ k given in (2.17), reproduced here for convenience These equations have been analyzed in Appendix C.2.1 and Appendix C.2.2, respectively.

A.2 Perturbations of Reissner-Nordström-AdS d+1
The class of equations we analyze for the perturbations of planar Reissner-Nordström-AdS d+1 with line element (C.13) are the following: Designer scalars: The perturbations in the tensor sector, or generic massive perturbations of the black hole are captured by (A.2) with the replacement of f (r) to (C.13) [32].

Vector polarizations of gravitons and gauge field:
The vector perturbations of metric and gauge field (parity preserving) are captured by equations for two fields, one dual to the momentum current, and the other dual to transverse photons, denoted X and Y, respectively. 14 These fields have a non-trivial dilatonic modulation of their kinetic terms, and a potential which depends on the spatial momenta. The equations as re-derived in [32] are The parameter C is defined as which in turns defines a deformed spatial momentum variable p v Note that we have added a v subscript to this momentum variable (relative to [32]) to distinguish it from the deformation that appears in the scalar sector.

Scalar polarizations of gravitons and gauge field:
The scalar sector equations are a bit more complicated. They are given in terms of two fields V and Z which correspond to the charge diffusion and energy density mode, respectively. These two fields obey the following decoupled equations as re-derived in [33] In this case we define the deformed momentum parameter p s (note the subscript) as Note that the function Λ k continues to be given by the expression (2.17). The potentials are complicated expressions involving various combinations of background functions and the deformed momentum. For the field V we find the following expression for the potential V V : On the other hand the potential for Z takes the form As noted earlier, these expressions were derived in [44], but our presentation and rewriting of them should make the structure more transparent. It is important that the modulation function Λ k only appears in the potential terms in the dynamics of V.
We can rewrite these equations in the standard form (2.8). The coefficient functions for the SOLDE for V are given to be (A.12) For the field Z we instead find the coefficient functions . (A.13)

B Scalar metric and gauge field perturbations of the charged black hole
We had promised in the main text to relate these fields to the physical metric and gauge field perturbations. Writing the metric and gauge field for the Einstein-Maxwell theory as with the subscript '(0)' referring to the background solution. The metric perturbations are parameterized by fields {Φ E , Φ O , Φ W } as in (2.14), while the gauge field perturbation is captured by another field V These four functions are given in terms of V and Z through the following relations (we have chosen not to solve for the fields directly). First we redefine the metric functions using a field similar to Z encountered in the Schwarzschild-AdS d+1 example, taking into account the mixing between the gauge field and graviton perturbations: The fields V and Z themselves are finally given in terms of V and Z by a functional basis change in a set of coupled ODEs From our analysis in Appendix C.3.3 we noticed that while Z is regular at either set of roots r k1 and r k2 of Λ k , the function V was only regular at the former, but had a simple pole at the latter. From here we conclude that Z is regular, but V potentially has a simple pole.
However, upon using the local Frobenius solution for Z and V, we find however that both Z and V are regular at all the roots of Λ k . Furthermore, the fields {Φ E , Φ O , Φ W } themselves are regular with the residue at the poles coming from the Λ k factors canceling out. Apart from noting that this is indeed physically sensible -metric functions and gauge fields should not have divergences -we do not have an explanation for the conspiracy in the coefficients which achieves this. The best we can offer is to note that we could have anticipated on physical grounds that the function V has to be singular like Λ −1 k . The fact that it is only singular at the r k2 set of roots has to do with the relative coefficients in the linear combination (recall that the roots are controlled by p 2 s and p 2 s + 2).

C Analysis of black hole wave equations
As we undertake our analysis we will have the need to refer to various singular points (abbreviated SP). Some of them will occur for generic values of parameters corresponding to the background geometry and the linearized field, while some appear upon some fine-tuning. We will distinguish the two cases for convenience, isolating also the features that are specific to certain probes.

C.1 Terminology for singular points
The SPs of interest that appear at generic values of the parameters occur for the most part at geometric loci. The obvious locations are the AdS boundary, the black hole singularity, and zeros of the function f (including the horizon). We find it convenient to refer to these using the following terminology: • Asymptotic SP, refers to the point at the AdS boundary r → ∞.
• Curvature SP, associated with the black hole singularity.
• Horizon SP, located at the horizon for a non-degenerate black hole, and other roots of the metric function f .
In addition, there are some singular points that are specific to equations governing conserved currents, after repackaging the information into gauge invariant variables (as described above). We encounter • Energy density ASP: This is the potential SP from the zero of the function Λ k , (2.17) encountered in the analysis of scalar polarizations of gravitons (gauge fields).
• Ohmic SPs: These are specific to the Reissner-Nordström-AdS d+1 solution, and are determined by the Ohmic radius r Q introduced in [32] (related to zeros of f ′ ) typically lying inside the black hole horizon.
The SPs listed above exist for generic parameter choices, viz., for generic values of black hole parameters and at generic frequencies and momenta of probes. However, additional features arise as we vary parameters, either in the form of confluence, where two regular SPs merge and become irregular, or coalescence, where the merger remains a regular SP. 15 Moreover, either for the geometric SPs or at confluence/coalescence, we may be able to convert the SP to be an ASP by a further fine-tuning of the parameters. For these, our terminology will be the following, defined in decreasing order of genericity: • Extremal confluence: The extremal limit where the horizon SPs become irregular. 16 • Matsubara ASP: These refer to special fine-tuned values of frequencies, specifically ω = −2πi n T with n ∈ Z + , where the horizon SP has some special features. At this locus the ingoing and outgoing solution both become analytic (at certain specific values of the momenta), resulting in the horizon becoming an ASP. This was noticed in [25,27] and explored in some generality in [28].
• Energy asymptotic coalescence: The equation of the scalar graviton polarization has additional features for translational homogeneous modes (k = 0) owing to the presence of soft gravitons [31]. Here the asymptotic characteristic exponents change, whilst however retaining the regular singular nature of the SOLDE.
• Energy horizon coalescence: Refers to the situation when the horizon SP merges with the energy density ASP. This occurs at a codimension-2 locus, at fixed ω, k, and has been analyzed hitherto in the context of chaotic dynamics of black holes [25].

C.2 Schwarzschild-AdS black hole wave equations
As remarked above, some of the general features we explore below have been uncovered in the literature. There are a few situations primarily involving the charged Reissner-Nordström-AdS d+1 black hole where there is a qualitatively new behaviour. For completeness, however, we will describe the neutral case first, phrasing the results in a manner which resonates 15 The adjective confluence is used traditionally in the theory of differential equations to characterize mergers of singular points which lead to essential singularities in the solution (hence irregular SP). When the merger does not change the nature of the SP we choose to refer to it as a coalescence (even when it changes the characteristic exponents). 16 This may also involve a merger with the Ohmic singular point.
with our understanding of the physics. The background solution is the planar Schwarzschild-AdS d+1 black hole whose metric takes the form (2.1) with

C.2.1 Designer scalar ODEs
We begin with the simplest equation (3.1) capturing the dynamics of a designer scalar field. It has the geometric singular points at r = ∞, r = 0, and at f (r) = 0. This comprises d + 2 singular points for a probe of Schwarzschild-AdS d+1 , rendering it as a Fuchsian equation with the given number of regular singular points [40]. The local behaviour of the differential equation in the neighbourhood of the singular points can be summarized as follows: • Asymptotic SP 17 • Horizon SP 18 We have written the result in dimensionless frequency and momenta, w and q, natural from a bulk perspective. 19 The singular points are depicted in the complex radial plane in Fig. 2. At the asymptotic SP, the characteristic exponents are determined by where, by parameterizing the mass appropriately in terms of a 'dimension' ∆, we have identified the behaviour with that of a scalar field in a AdS spacetime of effective dimension d eff = M + 1. This system has been analyzed in [10]. 17 We have directly expanded the coefficient functions to determine the behaviour at the asymptotic SP. A more natural way to deduce this would have employed a fractional linear map ϱ = r + r to bring the point to the origin. 18 When we write the local behaviour at the horizon SP, we are going to do so only for the largest positive real root of f . If we wish to examine the behavior at a different zero of f , say r = r+ϖ d , for some d th root of unity, then the expression (C.4) holds provided we also rotate our dimensionless frequencies and momenta appropriately, viz., replace r+ → ϖ d r+ whilst also defining (ω, k) = ϖ d r+ (w, q). 19 For the Schwarzschild-AdS d+1 solution the two dimensionless frequencies are simply related as w = 2 d w.
At the curvature SP, the indicial equation implies that the characteristic exponents are 0 and d − 1 − M. Minimally coupled scalar with M = d − 1 always has a branch cut at the singularity (since both exponents are equal). For d > 3, designer scalars with d − 1 − M ∈ Z + naively could find the origin to be an apparent singularity. Examples of such fields are scalar and vector polarizations of a probe Maxwell field in the background, which have M = 3−d and M = d − 3, respectively, and vector graviton polarizations which have M = 1 − d. However, one can check that despite the exponents being integral, there is always a logarithmic branch cut, and thus the black hole singularity is a regular SP of the SOLDE.
Finally, the horizon SP, are at the zeros of f , which depends on the geometry in question. For the planar Schwarzschild-AdS d+1 black hole one has d-roots at r = r + ϖ d , where ϖ d are the d th roots of unity. 20 At the horizon SP, the characteristic exponents are 0 and i w. The former is the analytic ingoing mode, while the latter is the outgoing mode. While this suffices for the local solution around the horizon SP, if we were to pick some closed form function representation for the solution (e.g., the hypergeometric function usually employed in the BTZ geometry), one should ensure that any branch cuts inherent there should be run in a direction away from the ray [r + , ∞). From these results, it follows that (3.1) is indeed a Fuchsian SOLDE with d + 2 regular SP. We recall here the Fuchsian sum rule, which demands that the sum of the characteristic exponents at all the SPs equal two less the number of SPs.
Matsubara ASP: If we consider the horizon SP, for iw = n with n ∈ Z ≥0 , then the second mode function, viz., the outgoing mode, can also be made regular at the horizon (this was noticed a long time ago in [54,55]). With this choice we make the two characteristic exponents non-negative integers, but a further fine-tuning is required to ensure the absence of any monodromy. This was explored in [28], and can be understood as follows. The Frobenius expansion for the outgoing mode despite starting out with an integral power, will generically include a logarithmic branch, rendering the mode non-analytic. This can be avoided, since we can tune the momentum q. Assuming one of the characteristic exponents to be zero (which can be done by factoring out a power), and the difference n ≥ 1, we will need to Laurent expand p φ and q φ to order n − 2. The coefficient γ of the logarithmic mode can be determined in terms of these Laurent coefficients, see Appendix E for details. Thus, setting a particular combination of the coefficients to vanish will suffice to ensure that we have an ASP.
For the horizon SP we have, for given n, a polynomial of degree 2n in q determining the non-trivial monodromy. Choosing q to be a root of this polynomial will ensure that we have a horizon ASP. Once this is done, there are no quasinormal modes, as there is no quantization condition to solve for (we have fixed both w and q); imposing normalizability at the boundary is ineffective since any linear combination of the modes on the horizon is acceptable. This is the reason why we have apparent quasinormal modes. Likewise, attempting to find the boundary Green's function leads to an ambiguity. This has been explained in [25,[27][28][29], so we will not elaborate further. The reader may find wave equations on the BTZ background studied recently [10] helpful to see how the ambiguity arises (since analytic expressions are available in this case).

C.2.2 The energy density operator ODE
Let us now look at (A.3) for the planar Schwarzschild-AdS d+1 black hole. This has been examined in d = 4 in Appendix B of [25] for some particular aspects, but we will highlight some additional features, which will be of relevance to our discussion. The differential equation is for a variable Z(r), with coefficient functions The local behaviour of the differential equation in the neighbourhood of the geometric singular points can be summarized as follows: • Asymptotic SP The residue of q Z at the horizon is (C.10) While this particular value is not relevant for the purposes of computing the characteristic exponents, it will be useful when we consider fine-tuned values of parameters.
From the asymptotic SP the characteristic exponents for Z for generic (ω, k) are 0 and d − 4, respectively, which was used in [31] to argue for its non-Markovian character. At the curvature SP, we have a degenerate pair of exponents (both zero), which implies that there is a logarithmic branch in the solution for Z near the black hole singularity. The behaviour at the horizon is characterized by an ingoing mode that is analytic (exponent 0) and an outgoing mode which is not -due to the exponent being iw -as is the case for any non-extremal black hole at generic values of frequency and momenta.

Energy density ASP:
The zeros of Λ k are potential SPs of the equation for Z(r). They are located at r k , given in (3.2), and turn out to be ASPs.
To see this we note that the characteristic exponents are 0 and 3, determined by the fact that p Z ∼ − 2 r−r k , while q Z starts off with a simple pole. Working out the coefficients to O(r − r k ) we find that the relation obtained in (E.3) always holds. Hence, as advertised in Appendix C.1 it is appropriate to refer to the roots of Λ k as energy density ASP. Thus, the singular points of (A.3) are d + 2 in number, coinciding for generic (ω, k) with the picture we had for the designer scalar wave equations Fig. 2. Taking the energy density ASPs into account, we find the sum of the characteristic exponents to be 4−d+3 (d−2) = 2 d−2, which is indeed what is demanded by the Fuchsian condition for an equation with d + 2 regular SPs and (d − 2) ASPs.
There are additional features at fine-tuned values of the parameters. We again have a discretuum set of Matsubara ASPs. Here we still tune w = −i n with n ∈ Z ≥0 . The tuning of k involves also the function Λ k leading to a more complicated condition than before, but one can ascertain that suitable choices of k can be made. Their physical interpretation follows along the lines described in Appendix C.2.1, so we will not elaborate further. Fine-tuning k to particular values, however, leads to some new features that are unique to the energy density ODE. We will describe two aspects of coalescence that are interesting below.
Energy asymptotic coalescence: While modes with k 2 ̸ = 0 have the asymptotic SP determined by (C.7), translationally invariant modes, with k 2 = 0, have a very different behaviour. This has to do with the function Λ k (r) which has non-trivial momentum dependence. The explicit k 2 dependence in the potential does not affect this argument just as the momentum dependence does not enter the asymptotic fall-offs of fields in AdS spacetimes. In any event, setting k 2 = 0 we find the locus r = ∞ is still a regular SP, albeit with a modified set of characteristic exponents, r 0 and r −d , respectively. This behaviour is akin to that of a minimally coupled, massless scalar field in the geometry. The Fuchsian nature of the equation at this locus is more obvious, reducing as it does to the earlier analysis of designer scalar equation.

Energy horizon coalescence:
We have a set horizon SPs arising from the roots of f , and a set of energy ASP from Λ k . We can ask when two of these coalesce. As with the energy asymptotic coalescence, it is reasonable to expect that the horizon remains a regular SP, perhaps with some changed exponents. Focusing specifically on the largest real zero of f , we see that Λ k will have a root at r + provided we set With this choice we then notice that the coefficient functions become The main change here is the double pole in q Z which arises from the merger of the two zeros (there is also a shift in the residue of p Z ). The indicial equation for the power law (r − r + ) α , near the singular point, now indicates that the characteristic exponents shift to α + = 1 and α − = 1 + iw. This is all that happens for generic frequencies.
However, one might imagine that as with the Matsubara ASP, we can tune the frequency to make the α − mode regular at the horizon, say by choosing w = −i (n − 1) for n ∈ Z ≥0 . Note however, that the momentum has been fixed already in (C.11), which implies that we don't have too much fine-tuning freedom. Not only, do we have to make sure that the characteristic exponent is integral, but we also have to check that the Laurent coefficients work out to remove the logarithmic branch. This turns out to be possible, but only at a particular frequency, w = i, i.e., ω = 2πiT . At this point, both modes are analytic, and the horizon locus is an ASP, as we can already see from (C.12). For any other choice of n, the solutions have a logarithmic branch cut, and the outgoing solution picks up a monodromy.

C.3 Reissner-Nordström-AdS black hole wave equations
We now turn to the discussion of wave equations in a charged black hole background corresponding to a boundary CFT plasma at non-zero chemical potential. The background line-element is still given by the translationally invariant (2.1), now with .
The physical parameters, temperature and chemical potential, are related to r + and Q through We also introduce the Ohmic radius r Q , Ohmic function h, and the dc conductivity parameter S Q following [32], for they will play an important role in our analysis. These are given by The wave equations of interest are collected in Appendix A.2.

C.3.1 Designer scalar ODEs
The designer scalar ODE in the Reissner-Nordström-AdS d+1 black hole background is similar to the that studied in Appendix C.2.1. It is relevant for the study of massive scalar probes, tensor graviton perturbations of the black hole (which again are massless minimally coupled scalars), both of which have M = d − 1. We will however analyze the equations for arbitrary M.
The analysis of the designer ODE (3.1) in the planar Reissner-Nordström-AdS d+1 is similar to our earlier discussion. The key difference is that the function f (r) has 2(d − 1) zeros from (C.13). Let us first consider a non-degenerate horizon where the roots of f (r), and in particular the largest real-root, are non-degenerate. The asymptotic SP is unaffected by the charge, and the analysis at the horizon SP is similar to the Schwarzschild case (accounting for the charge dependence of the temperature). Specifically, • Asymptotic SP • Curvature SP 21 In this case using the Matsubara normalized frequency results in cleaner expressions, since T (r + ) is more involved. 22 The reader can check that (C.4) and (C.19) are the same accounting for the translation between horizon size and temperature. Note that we are again focusing on the largest positive root of f to define the horizon SP behaviour. As indicated in footnote 18 we can carry out an analogous exercise for the other roots, but we now will have to work with the locally defined (complex) temperatures if we wish to write the expressions in terms of the dimensionless variables.
Since the Reissner-Nordström-AdS d+1 black hole has a different causal structure, in particular a timelike singularity, the behaviour at the curvature SP is different. In particular, the characteristic exponents are 0 and 2d − 3 − M, respectively. For a minimally coupled massive scalar, with M = d − 1, not only are both exponents are integral, but the logarithmic branch can also be checked to be absent. Thus, in this case, we have a curvature ASP, with the fields not really being sensitive to the black hole singularity -however, some components of the conserved currents have duals that are cognizant of the singularity, cf., Appendix C.3.2. We believe this should hold for a range of M, but have not checked in detail whether it does. One can again check that the designer scalar wave equation in Reissner-Nordström-AdS d+1 background is a Fuchsian SOLDE, with 2d SPs (we include the curvature SP in this count even when it is an ASP). Essential features of this equation are summarized in Fig. 4. 21 We are assuming the limit is taken with Q ̸ = 0. The order of limits is important, as higher terms in the Laurent expansions have singular behaviour in the limit Q → 0 (which makes sense since the function f diverges more rapidly for the charged solution). 22 Now w and w are related by a factor which is d−(d−2) Q 2

2
, which we prefer to avoid writing, and have hence chosen to work with both normalizations. Extremal confluence: There are other additional features of the Reissner-Nordström-AdS d+1 solution that have to do with degenerate horizons, and mergers of zeros of f (r). For an extremal black hole when the inner and outer horizon radii are equal 23 the horizon singular point becomes irregular due to a confluence of two regular singular points (say at r = r + and r = r − ). 24 We will not analyze this limit any further.
We next turn to the vector and scalar perturbations of the metric and gauge field in Einstein-Maxwell theory about a Reissner-Nordström-AdS d+1 background. The equations have been compiled in Appendix A.2.

C.3.2 Momentum diffusion ODE
In the vector sector we have two equations, one for the momentum diffusion mode, X and the other for the transverse photons Y, cf., (A.5). The former is of the form of a designer scalar with an exponent M = 1 − d but with a more involved dependence on the spatial momentum through p v defined in (A.7). We denote the coefficient functions appearing in the standard form of the SOLDE as {p X , q X } and {p Y , q Y }, respectively.
The generic SPs of these equations have the following behaviour: • Asymptotic SP: X behaves as a designer field with index M = 1−d, while Y behaves as one with index d − 3. The potentials do not affect this analysis (but enter in the determination of counterterms). 23 Note, however, that the proper distance between the inner and outer horizon remains finite in the limit. 24 If we play formal games of working with Q ∈ C then we can find scenarios where the complex roots of f (r) coalesce. A simple example seems to be d = 4 and Q = i 2 , for which the roots are at ±r+, ±i r + √ 2 , each of the latter two being two-fold degenerate. The conventions are as in Fig. 4; we have now made clear that the origin is a regular SP, and indicated the Ohmic SPs for this equation as well. Note that these lie between the inner and outer horizons [32].
• Asymptotic SP: Both V and Z behave as designer fields with index 3 − d. The potentials again are irrelevant for this analysis, and these fields behave as non-Markovian modes with characteristic exponents of 0 and d − 4.
• Curvature ASP: The black hole singularity continues to be a singular point of the Z and V SOLDE. The coefficient functions at most have simple poles. Specifically, The characteristic exponents are integral and there is no logarithmic branch cut in the solutions. Hence, black hole singularity is an ASP of the equations for both Z and V.
• Horizon SP: The behaviour of p V and p Z is as for any other field in a non-extremal background; a simple pole with residue 1 − iw. The functions q V and q Z also only have simple poles, but with a residue that is depends non-trivially on p s and S Q (or Q).
• Ohmic SP: The fields V and Z have Ohmic SP, since their kinetic terms are modulated by h. The functions p V and p Z have a simple poles with residue ∓2, respectively, while q V and q Z start off with a double pole, as can be read-off directly from (A.12) and (A.13).

Energy density (A)SP:
In the Schwarzschild-AdS d+1 analysis, we saw that the roots of Λ k , which naively appear to be SPs of the SOLDEs, are actually not. We will confirm this to be the case for the Reissner-Nordström-AdS d+1 equations as well, albeit in a slightly complicated manner. Notice that the vanishing loci of Λ k are located at the following 2(d − 2) points of the complex radial plane. (C.28) The first set here is continuously connected to the locations where the corresponding function for the neutral black hole has roots, while the second set is unique to the charged black hole.
To deduce this, we ask which of the two roots can be made to lie on the ray [r + , ∞) along the real axis by choosing an appropriate phase for p s . Since p 2 s ∼ q 2 for small momentum it follows that the first set of roots connect to the neutral case. Therefore, for small q the first set of roots lies outside the circle |r| = r + , while the second set lies within. At large charge however, the situation is reversed, with the exchange between them occurring at a critical charge When S Q < 1 2 , we need to ascertain whether the first set of zeros of Λ k lies between the horizon and the boundary, but the question switches over to the second branch for larger charges, with S Q = 1 2 being a critical point. For the present however, we can simply analyze the two sets of roots in turn, not worrying about whether they lie within or outside the circle of radius r + .
Let us start with the particular root r k1 = r Q − p 2 s 2 − 1 d−2 (fixing the phase). For the field Z we find that the function p Z has a simple pole with residue −2. The function q Z also has only a simple pole, leading to characteristic exponents 0 and 3, respectively. We then expand the functions to linear order around this putative SP to check if there is a potential logarithmic term, and again find that it is absent, rendering the singular point apparent. The analysis for V is a bit simpler. At this zero of Λ k one can check that the functions {p V , q V } are, in fact, manifestly regular, thus rendering the zero set to be ordinary points of the SOLDE.
When we pick the root r k2 = r Q 1 + p 2 s 2 − 1 d−2 (with phase fixed to unity), we find that p Z has a simple pole with residue −2, while q Z has a double pole with residue 2. The characteristic exponents are 1 and 2, respectively. Examining the constant term in p Z we find that it is negative of the residue of q Z at its simple pole term. This guarantees the absence of any monodromy around this SP, rendering it also to be an apparent SP. For the field V however, we find a double pole in q V at this root, with residue −2. This implies that the characteristic exponents are −1 and 2. We have checked that there is no monodromy around this point as well. However, the solution for the field V does have a simple pole at these roots of Λ k . Ergo, the set of roots r k2 is a genuine SP of the SOLDE. Figure 6: Singular points of the longitudinal scalar mode V, which is dual to the charge diffusion mode photon in Reissner-Nordström-AdS5 in the complex r plane. The conventions are as in Fig. 5. The new elements here are the energy density SPs from the d − 2 roots of Λ k in the set r k2 , which we have labeled as {r k2 ,r k2 }. For the latter, we have chosen non-zero value for arg(k 2 ) to aid with the depiction, but note that as we rotate k we can make the SP pinch the grSK contour.
The above completes the description of the singularity structure for generic values of parameters for non-extremal Reissner-Nordström-AdS d+1 black holes. There are again features specific to fine-tune loci in parameter space. Some of these are as before, viz., • Extremal confluence: When we take the extremal limit, demanding that f has a double zero at r = r + , we convert r = r + to an irregular SP of the equations.
• Energy asymptotic coalescence: As in the Schwarzschild-AdS d+1 case we have a special situation when the spatial momentum vanishes for the energy density mode Z. Owing to the presence of Λ k in the kinetic term, at k 2 = 0, we find coalescence of the roots of Λ k which are now at the asymptotic boundary. This leaves the asymptotic singular point regular, but changes the characteristic exponents to 1 and −d, respectively. Note that this behaviour is not present for the charge diffusion field V, since its kinetic term does not have a corresponding factor of Λ k .
• Matsubara ASP: As in previous equations we can tune to specific frequencies iw = n with n ∈ Z ≥0 . At these values the outgoing mode can be made regular, provided we remove the monodromy due to the logarithmic branch of the solution in the outgoing mode. To achieve this we have to tune k 2 , or equivalently p s , which we have used to parameterize both the Z and V equations.
Energy horizon coalescence: Now that we have understood the behaviour at the horizon, and at the roots of Λ k , let us turn to the case where we tune the momenta for the singular points arising from these two to merge. This will happen whenever we have tuned the momenta to be such that the root of Λ k lies at r = r + . It is easy enough to determine when this happens, we need either based on the two branches of zeros. One can check that this is consistent with fixing the physical momentum using At this point it becomes clear that the jump is due to the fact that there is branch structure in the momentum dependence, owing to the presence of p s . Let us start with the energy density equation for Z and examine the two roots in (C.30) in turn. When p 2 s = −2 S Q we find the local behaviour of the coefficient functions to be On the other hand for the root p 2 Now, for w = −i we see that the residue of the simple pole in q V vanishes, allowing for this coalescence to turn itself into an apparent singularity. Finally, at S Q = 1 2 , we simply encounter a regular SP, despite the double-root of Λ k (r), as noted in (C.35). Now we find (C.40) For w = −i we again find an apparent singularity -in fact, this behaviour seems to smoothly connect to the second branch of Λ k roots unlike what happened for the Z equation above.

D Form factors in the FPF correlator
We collect in this appendix the formulae for the functions J reg and J loc , which enter into our analysis of the three-point function evaluated in §4.4. Let us start with J reg entered the regular piece of the three-point function. This function involves completing the Mellin-Barnes contour integrals, which, we recall, were introduced to allow one to complete the radial integral over a product of hypergeometric functions. There are three propagators, and hence, three contour integrals. These are evaluated by standard residue calculus, and lead to the following infinite sum representation equation. We will focus on a particular form inspired by the equations we want to study, and w.l.o.g. place the singular point the origin. The equation we study is therefore q m z m y(z) = 0 (E.1) We have expanded the coefficient functions in a Laurent series, fixing the polar terms based on the structure we encounter for Fuchsian equations. The characteristic solutions are z 0 and z n , and so the general solution is y(z) = Ay n (z) + B(y 0 (z) + γ y n (z) log z) .

(E.2)
Here y n (z) and y 0 (z) are the naive Frobenius solutions y n (z) = z n ∞ j=0 a j z j and y 0 (z) = ∞ j=0 b j z j , respectively. We wish to understand the conditions under which the log term can be made absent, or equivalently how the monodromy matrix set to be the identity matrix. An early analysis can be found in [38].
One strategy is to truncate the Laurent expansions for the coefficient functions and attempt to solve the SOLDE; this leads to closed form solutions for n = 1 and n = 2. The former is reliable, but the latter misses out some admixture from higher order terms in the coefficient functions. A more straightforward approach is the following: setting a 0 = 1, fix a 1 by solving the SOLDE. Then use the data to solve for the second independent function, and extract the condition for γ = 0, as a constraint on the series coefficients. Carrying out the exercise we find the following condition for the monodromy γ to vanish: n = 1 : q −1 = 0 , n = 2 : q 2 −1 + p 0 q −1 + q 0 = 0 , n = 3 : q 3 −1 + 2 p 0 q 2 −1 + 2 q −1 (p 2 0 + p 1 + 2 q 0 ) + 4 (p 0 q 0 + q 1 ) = 0 (E. 3) The larger the difference in the characteristic exponents, the higher one needs to go in the series solution. Correspondingly, one is sensitive to more terms in the Laurent expansion of the coefficient functions about the singular point. This point has been explained in [28] for some specific cases; they express the result as a vanishing condition for a determinant of Frobenius series coefficients. The analysis presented here can be adapted to any regular singular point.