Horizon instability of the extremal BTZ black hole

We study real-time propagation of a massive scalar field on the extremal BTZ black hole spacetime, focusing on the Aretakis instability of the event horizon. We obtain a simple time-domain expression for the AdS3 retarded Green function with Dirichlet boundary conditions and construct the corresponding time-domain BTZ retarded Green function using the method of images. The field decays at different rates on and off the horizon, indicating that transverse derivatives grow with time on the horizon (Aretakis instability). We solve the null geodesic equation in full generality and show that the instability is associated with a class of null geodesics that orbit near the event horizon arbitrarily many times before falling in. In an appendix we also treat the problem in the frequency domain, finding consistency between the methods.


Introduction
In the decade since Aretakis' initial study of massless scalar fields in the extremal Reissner-Nördstrom spacetime [1,2], it has become clear that extremal horizons generically exhibit weak derivative instabilities [3][4][5]. Independent of the type of extremal black hole (and even of the type of field that perturbs it!), the evidence suggests that sufficiently highorder transverse derivatives always grow at least polynomially in advanced time along the event horizon. This has the physical consequence that infalling observers experience large field gradients [4,6,7].

JHEP05(2020)094
This general picture has emerged from a variety of techniques, whose domains of validity are largely non-overlapping. The original technique of Aretakis [3,8,9] involves a conserved quantity along the horizon that provides an obstruction to decay. This technique has been used only in what we call the "discrete case" (e.g., axisymmetric massless perturbations of Kerr [5]), and mainly when initial data extends to the horizon. We ourselves have been involved with a different technique based on frequency-domain analysis, in which the instability is associated with a singular branch point in each mode of the retarded Green function [6,[10][11][12]. This technique has been used mainly in the non-discrete case (e.g., non-axisymmetric mode perturbations of Kerr [6]), and only when initial data is confined away from the horizon. Other techniques rely on discrete inversion symmetries present only for some black holes [4,[13][14][15]. Analytic arguments involving the near-horizon geometry [5,6,16,17] together with numerical work in a variety of settings [4,18,19], have provided complimentary insight into these various regimes. However, a fully general understanding of the instability remains elusive.
There is a dissonance here: a universal phenomenon should be simple at its core, yet the state of the art presents a sprawling complexity. One suspects that a key unifying idea has not yet been identified. In such a situation, it pays to study the effect in the simplest possible setting, where the nature of the phenomenon might be revealed with a minimum of distraction. Buoyed by this hope, in this paper we initiate the study of the Aretakis instability in the exceptionally simple setting of a 2+1-dimensional black holethe Bañados-Teitelboim-Zanelli (BTZ) spacetime [20,21]. As this black hole has negative cosmological constant, the results may also help understand the holographic implications of the instability [22,23].
The BTZ black hole allows exceptional analytic control because it is locally isometric to three-dimensional Anti-de Sitter space (AdS 3 ), a maximally symmetric spacetime. We focus on the retarded Green function of a massive scalar field satisfying Dirichlet boundary conditions. Using the results of [24] (and correcting some minor computational errors), we write the AdS 3 retarded Green function in a very simple form. We then use the method of images [25] to construct the BTZ Green function.
The Aretakis instability emerges from this image sum in a beautiful way. Each term in the sum decays at the same rate, fixed by the AdS 3 spacetime. For a field point off the horizon in BTZ, the high-order images are unimportant, and the full Green function (and hence field) also decays at this AdS 3 rate. However, when the field point is chosen on the BTZ horizon, the high-order images become important and the infinite sum decays at precisely half the AdS 3 rate, 1 as expected on general grounds [5]. Since the field decays at different rates on and off the horizon, sufficiently high-order transverse derivatives on the horizon must grow -the Aretakis instability.
By solving the null geodesic equation in BTZ, we are able to identify the late-timerelevant terms in the image sum with the arrival of wavefronts that have orbited the event horizon arbitrarily many times before falling in. The number of orbits increases linearly JHEP05(2020)094 as the geodesic conserved quantity approaches that of the event horizon. We conjecture that this behavior is part of the following larger pattern for black holes: the surface gravity of the event horizon functions as a Lyapunov exponent for the deviation of nearby null geodesics; and when this exponent vanishes for an extremal black hole, the exponential deviation goes over to a power law. One can then say that wavefronts linger far longer near the event horizon of extremal black holes, as compared to analogous non-extremal black holes, providing an interpretation of the instability of extremal horizons.
Our results offer a new way to think about the Aretakis instability, but they do not solve the problem of unifying the previous different approaches. To provide a point of comparison, in the appendices we apply two other methods to the extremal BTZ spacetime: the conserved charge method (appendix B) and the mode sum method (appendix C). We cannot compare with the conserved charged method because our analysis is restricted to the non-discrete case, with initial data not extending to the horizon. We can compare with the decay rate of individual modes in the frequency domain analysis (and we find agreement), but this method does not give the behavior of the full field (see further discussion in appendix C). Thus there remains much mystery surrounding the Aretakis instability of extremal horizons. This paper is organized as follows. In section 2 we review the AdS 3 and extremal BTZ spacetimes, introducing notation. In section 3 we consider a massive scalar field and construct the retarded Green function by the method of images. In section 4 we discuss late-time decay, and in section 5 we interpret the instability in terms of null geodesics. Our signature is (− + +) and we set the AdS radius to one.

Metric
We now discuss the AdS 3 metric and its periodic identification to produce the extremal BTZ black hole. We introduce a few different coordinate patches that are useful in the analysis that follows.

Poincaré coordinates and patch
"Poincaré coordinates" for AdS 3 are given by t = sin τ cos τ − cos Ω sin χ , z = cos χ cos τ − cos Ω sin χ , x = sin Ω sin χ cos τ − cos Ω sin χ , (2.2) JHEP05(2020)094 Figure 1. AdS 3 and its patches, plotted using global coordinates (τ, χ, Ω) as cylindrical coordinates with height τ , radius χ, and angle Ω. The maximally extended spacetime has a cylindrical conformal boundary at χ = π/2, shown here in translucent orange. On the left, we show the Poincaré horizons τ ± that bound the Poincaré patch. In the middle, we show τ − along with the future horizon τ H (which is also a Poincaré horizon); these bound the "extremal patch" that becomes the BTZ exterior after identification. Finally, on the right we show the same plot from a different perspective, along with a selection of integral curves of the Killing field ξ = ∂ X . The extremal BTZ black hole is obtained by identifying discretely along these curves.

Extremal coordinates and patch
We also introduce "extremal coordinates" for AdS 3 , JHEP05(2020)094 bringing the metric into the form These coordinates cover the region R > R H , which is bounded by the past Poincaré horizon τ − and a globally translated and rotated Poincaré horizon τ H defined by 2 τ H = arcsin (sin Ω sin χ). (2.7) We refer to this patch τ ∈ (τ − , τ H ) as the extremal patch (figure 1 middle and right). We call the bounding null surfaces τ − and τ H the past and future horizons (respectively), since these will become the event horizons of the extremal BTZ black hole after the identification X ∼ X + 2π. These horizons are described by R = R H in extremal coordinates, but the patch itself is independent of R H . 3 The inverse transformation is given by

Horizon coordinates
Finally, it will be useful to consider "horizon coordinates" (V, R, Φ) defined by where the metric becomes (2.10) These coordinates cover the future horizon R = R H , with the null generators given by Φ = const. When Φ ∼ Φ + 2π, these are ingoing, corotating coordinates for the extremal BTZ black hole. The horizon-generating Killing field is given by In global coordinates, R > RH corresponds to αβ > 0 with α = sin Ω sin χ − sin τ and β = cos τ − cos Ω sin χ. The zeros of β bound a set of Poincaré patches including our choice (2.4), while the zeros of α bound an interleaving set related by the symmetries τ → τ + π/2 and φ → φ + π/2. As our Poincaré patch (2.4) has β > 0, we require α > 0 as well, giving rise to τ ∈ (τ−, τH ). 3 Note that the original transformation given by BTZ does entail a different patch for each value of RH ; these patches approach Poincaré patch as RH → 0. Our extremal patch remains distinct in the RH → 0 limit. JHEP05(2020)094

Extremal Killing field and BTZ
The extremal BTZ black hole is obtained by identifying points separated by parameter distance 2π along the integral curves of the Killing field This means we identify X ∼ X + 2π in extremal coordinates (2.6) and (equivalently) Φ ∼ Φ + 2π in horizon coordinates (2.9). We will call ξ the extremal Killing field. Its integral curves connect the two ends of the "conformal bifurcation line" where the future and past horizons meet at the boundary (figure 1).

Scalar field
We consider a massive scalar field ψ, The AdS 3 /BTZ spacetimes are not globally hyperbolic, and well-posed evolution requires specification of the behavior of the field on the boundary. We will consider so-called "Dirichlet" conditions, where ψ is required to vanish on the boundary. We define an associated retarded Green function by together with the requirements that G vanish when either point approaches the boundary or if x is not in the causal past of x. (Here δ 3 is the covariant delta distribution, equal to 1/ √ −g times the coordinate delta function). We assume µ 2 ≥ −1 so that the Dirichlet dynamics are well-posed [26,27]. This Green function can be used to construct the field from its initial value ψ 0 (R , Φ ) on the null surface V = 0 via the Kirchhoff representation 4 where we have assumed that the field point is on or outside the event horizon. In AdS 3 , the range of Φ is unbounded, while in BTZ it must be restricted to a fiducial range such as 0 < Φ ≤ 2π. In all other respects, eqs. (3.1), (3.2), and (3.3) hold equally well in AdS 3 and extermal BTZ.

AdS 3 Green function
The Dirichlet retarded Green function for AdS 3 was obtained in closed form by [24]. In appendix A we review the derivation, correct some trivial errors, and obtain a new form of the result: show the wavefront of a point source (surface of singularity of the retarded Green function), which "bounces" off the boundary. The wavefront is described by Σ = 0 (orange) and Σ = 2 (blue). On the right, we show the behavior of the retarded green function as a function of Σ. Two qualitative behaviors exist: the "generic" case (middle; we have chosen ν = 1.2) and the "exceptional" case (we have chosen ν = 0.5). In the exceptional case, the Green function vanishes for Σ > 2.
Here ν is a scaling dimension for the field defined by while Σ is a biscalar on AdS 3 given by The Green function is singular on the surfaces Σ = 0 and Σ = 2, whose union is the null wave front of the propagating field. This wave front bounces off the boundary as shown in figure 2. In the "exceptional" case where ν is a half-integer, the Green function vanishes for Σ > 2, so that the returning wavefront appears to "cancel out" the propagating field. Although the wavefronts appear conical in the plots, they are not precisely equal to cones. The distinction is most graphically evident when the source point is near the boundary, in which case the wavefronts collapse to Poincaré horizons.
We now discuss the geometric interpretation of Σ and its role in the Green function. In a Poincaré patch there is always a unique geodesic connecting two points, and Σ is related to the world function σ (one-half the squared geodesic distance) by [29] In particular, Σ > 0 is timelike separation, Σ < 0 is spacelike separation, and Σ = 0 is null separation. However, we have seen that Σ = 2 also represents a null wavefront of the propagating field. This is not a contradiction, since the bouncing behavior of the wave fronts is not reflected in the behavior of geodesics, which simply continue indefinitely toward JHEP05(2020)094 the boundary as the affine parameter increases. Thus, at least within the Poincaré patch, points for which Σ = 2 are indeed connected by only a single geodesic, which is timelike, despite the fact that they are also connected by a bouncing wave front of the scalar field. The null character of the surface Σ = 2 (fixing the prime point) can be seen directly by observing that it corresponds to Σ = 0 for "conjugate points" τ → τ + nπ, Ω → Ω + nπ for odd n (see figure 2). Finally, note that the wavefront preserves the inverse square root character of the singularity in the Green function inherited from the local Hadamard form. See e.g. [30] and references therein for rigorous results on the global propagation of singularities.

Extremal BTZ Green function
As the BTZ black hole is a periodic identification of AdS 3 , propagation in BTZ is equivalent to propagation on a portion of AdS 3 subject to periodic boundary conditions. The retarded Green function can then be simply constructed from the non-periodic one by

JHEP05(2020)094
where the last equation holds only in horizon coordinates. The notation e 2πnξ x indicates the spacetime point obtained by flowing along the integral curves of ξ for a parameter distance 2πn; when horizon coordinates are used, this simply increases Φ by 2πn. It is straightforward to check that this expression satisfies all the conditions for the Dirichlet retarded Green function, together with the additional required periodicity. From the perspective of AdS 3 , we place "image sources" at distances of 2πn along the integral curves of ξ from the original source position (figure 3), all of which emit wave fronts. From the perspective of BTZ, the arriving fronts are interpreted as having circled the black hole |n| times. This method of constructing the Green function is generally called the method of images; it has been used extensively in Euclidean signature, but apparently not for retarded propagation.
It will be convenient to define a separate Σ n for each image n. Letting then (expressing in horizon coordinates) we have The other component of the AdS 3 retarded Green function (3.4) is a theta function, Θ(t − t ), that cuts out the past of t . We will introduce the notation indicating that one is to express t and t in terms of horizon coordinates and then send Φ → Φ + 2πn. Explicitly, we have Putting everything together, the BTZ Dirichlet retarded Green function is written where Σ n and s n are given in eqs. (3.11) and (3.13), respectively. The non-uniform behavior of Σ n in the limits V → ∞, n → ±∞, R → R H and R → R H is key to the phenomena discussed below.

Late-time decay
We now use the explicit retarded Green function to investigate late-time decay of field perturbations. We will discuss decay in null directions, varying V while fixing R and Φ. When R = R H this corresponds to decay along a horizon generator Φ, measured by affine parameter V . When R > R H the null ray remains outside the horizon. JHEP05(2020)094

AdS 3
We begin with the case of AdS 3 . The AdS 3 Green function is simply the n = 0 term in the sum (3.14). At large δV we have where we exclude the measure-zero case C 0 = 0. If C 0 < 0 then Σ 0 is negative at late times and the Green function vanishes: the points of fixed R and Φ with large V are out of causal contact with the source position R , Φ , V . When C 0 > 0, the Green function at late times is If ν is a half-integer then this expression vanishes. In fact the Green function (3.14) vanishes identically in this case, where V 0 is some constant. By the Kirchhoff integral (3.3), the field sourced by generic initial data will behave the same way, where C is some function determined by the initial data, and we introduced in order to be consistent with previous work [5][6][7]. Notice that h ≥ 1/2 given our assumption on the mass µ 2 ≥ −1. We will refer to the case 2h / ∈ Z + as "generic" and the case 2h ∈ Z + as "exceptional". This exceptional case includes the case h ∈ Z + previously called "discrete".

BTZ
For a BTZ black hole, we must consider the infinite sum (3.14). Each term in the sum is an AdS 3 Green function with a different source point, so each term behaves like sin(2πν)V −2h at late times as above. However, this does not guarantee that the full sum shares this behavior, and indeed we will see that it does not. To examine the full sum, it is helpful to rewrite (3.11) as  . The slower decay is due to spikes from wavefronts that orbit many times near the black hole and cross the horizon at arbitrarily late times. Spikes also arrive at arbitrarily late times from wavefronts that spend their time mainly far away from the black hole, but they are exponentially narrow and are not visible in the plots at late times (and do not affect the decay envelope). In these plots we have chosen R H = 0.01 and ν = 1.2 such that h = 1.7, with δΦ 0 = 0 and R = 10; in the left plot we have R = R H = 0.01, while in the right plot To explore the deviation of the full sum from the behavior of an individual term we consider the regime of late times and large image numbers, i.e. large δV and large |n|. We restrict to the generic case 2h / ∈ Z + .

Both points outside the horizon
Suppose now that both the source and field point are outside the horizon, R > R H and R > R H . Then B ± are both non-zero, so in the regime of large δV and |n| we have The contribution of the n th term to the Green function is important when Σ n is positive and of order unity (see figure 2), with the initial spike occurring when Σ n = 0. Since B ± are positive, we see that the negative n terms do not contribute (Σ n is negative), while the positive n terms are important for times of order the image number (δV ∼ n). The importance of the contribution can be assessed from the derivative ∂Σ n ∂V ≈ B + e 2πnR H , (n > 0). (4.12)

JHEP05(2020)094
Since this derivative is exponentially large, the width (in V ) of the region of importance is exponentially suppressed with image number. This suggests that the large-n terms are not very important, and (aside from brief exponentially narrowing spikes) the full behavior of the sum will share the falloff V −2h of each individual term. Plotting the Green function numerically supports this conclusion ( figure 4 left).
The behavior of the field ψ follows from that of the Green function via the Kirchhoff integral (3.3). It is clear from the discussion above that the Green function decays as V −2h apart from exponentially narrowing spikes that may be safely excluded from the integral. This late-time approximation is uniformly valid over any region of R outside (but not including!) the event horizon, so the field will similarly decay as V −2h provided the initial data is confined outside the horizon. The spikes in the Green function will be smoothed out by the integral into finite oscillations about this decay. That is, we claim that the generic outside-horizon, late-time behavior of fields sourced by initial data outside the horizon is where ψ 0 is an arbitrary function of R, Φ and an O(1) function of V .

Field point on the horizon and source point outside
If the field point is now on the horizon (R = R H ) while the source is still outside (R > R H ), then B − = 0 and at large δV and |n| we have which may be compared with eq. (4.11). The formula for the positive n terms is identical, and these terms again contribute at δV ∼ n. However, now the negative n terms contribute as well, starting at even later times δV ∼ e 4π|n|R H (where Σ n becomes positive). The derivatives in these regimes are given by n < 0 : which may be compared to eq. (4.12). The positive n terms again have a large derivative, indicating an exponentially narrow contribution to the Green function. However, the newly contributing negative n terms have a small derivative, indicating an exponentially wide spike in the Green function as a function of V . Thus the negative n spikes should be very important at late times, and there is no reason to expect that the sum has the same falloff as an individual spike. Plotting the Green function (figure 4 right) confirms the exponentially wide spikes and shows that the actual decay is V −h . We can understand this rate through the following heuristic argument. Let us label the exponentially wide spikes by m = −n, so that m is a positive integer. As discussed above, JHEP05(2020)094  .8)]. Focusing on the period of order δV ∼ e 4πmR H before the next spike has arrived, we may equivalently express this in terms of m or δV as recalling that m-independent constants are dropped. This shows the (δV ) −h decay visible in figure 4. In essence, each image decays as C m (δV ) −2h , but the effective amplitudes C m grow as (δV ) h , resulting in an overall (δV ) −h decay. By similar arguments made for eq. (4.13) above, we expect that the generic behavior of fields sourced by initial data outside the horizon is thus where ψ H is an arbitrary function of Φ and an O(1) function of V . For later comparison it will be helpful to have a more precise expression for the time of arrival of the n th spike. This is obtained by letting δV n and Σ n = 0 for δV , giving If the field point is near the horizon but not exactly on it, one expects a transient period of V −h decay followed by final decay of V −2h . Plotting the Green function confirms this expectation (figure 5). The field ψ sourced by initial data outside the horizon will similarly show this transition from V −h to V −2h when evaluated near the horizon.

Both points on the horizon
If both the source and field point are on the horizon (R = R = R H ), then we have the simple expression Σ n = 1 − cosh (R H (δΦ 0 − 2πn)) .

JHEP05(2020)094
In particular, the Green function is actually independent of V (apart from the causal factor Θ(s n )). Of course, for a source point at at any small distance off the horizon the Green function does decay, with the rate determined by whether the field point is on or off the horizon. That is, the late-time behavior of the Green function is highly non-uniform as the points approach the horizon. Determining the decay rate of fields with initial data that extends to the horizon would require a careful estimate of the Kirchhoff integral (3.3) using the full behavior of the Green function near the horizon. We leave this to future work.

Null geodesics
We now study the null geodesics of the BTZ spacetime in order to understand the behavior of the Green function that gives rise to the instability. Let p µ denote the four-momentum of the null geodesic, where p V > 0 is our time orientation. Using the null condition and the two Killing fields ζ = ∂ V (horizon Killing field) and ξ = ∂ Φ (axial Killing field), we have three conserved quantities, The quantity L is interpreted as the angular momentum, while E is the energy according to the "static" Killing field Note that E can be negative for trajectories inside the "ergoregion" R < √ 2R H , where ∂ T is spacelike [31]. The special case E = L = 0 corresponds to the horizon generators, Henceforth we will consider the region outside and including the horizon, We will further assume E = 0 and introduce which introduces an energy-rescaled affine parameter λ that increases toward the future.
The case E = 0, L = 0 can be handled straightforwardly by the limit b → ∞, and we shall see that the case E = L = 0 is recovered (more subtly) by b → 1. Solving eqs. (5.1) assuming R ≥ R H > 0 and E = 0, we find where s = sgn(E). The requirement that k R be real imposes an allowed region of R for each value of b. Taking into account the assumption (5.3), we find

JHEP05(2020)094
with turning point radius (for |b| > 1 only) Thus there are three kind of trajectories outside the horizon: transits from the (white hole) horizon to the boundary (|b| < 1), transits from the boundary to the (black hole) horizon (also with |b| < 1), and transits from the white hole to the black hole (|b| > 1). There are no trajectories that begin at the boundary, turn, and end at the boundary. The sign s = sgn(E) is fixed by the time orientation k V > 0 (equivalently p V > 0) and may be expressed directly in terms of b, as follows. Since | − bR 2 Given the allowed ranges (5.6) of R, we find that In the case b = 1 the sign is indeterminate, with the limits b → 1 ± having qualitatively different behavior. It is possible to solve eqs. (5.5a) entirely for the trajectory x µ (λ). Choosing the integration constants for simplicity, we find where ± is the sign of dR/dλ (i.e. + when outgoing and − when ingoing). The three integration constants may be restored by shifting V , Φ, and λ by (separate) constant values. With the above choices, the range of λ is |b| > 1, turning at λ = 0 (5.10c) As λ → ±R H /(1 + b) the particle approaches the horizon R → R H , whereas as λ → ±∞ the particle approaches the boundary R → ∞.
We now fix a starting radius R 0 > R H and compute the total change in time (∆V = V − V 0 ) and angle (∆Φ = Φ − Φ 0 ) before the particle enters the horizon. These quantities also depends on the initial radial direction (ingoing or outgoing) when |b| > 1, since the particle will encounter a turning point if initially directed outwards. Considering only trajectories that end at the horizon, We will label initially ingoing trajectories (present for all b) with "nt" (for "no turning point") and initially outgoing trajectories (present only for |b| > 1) with "t" for "has a turning point". Eqs. (5.9) show that the lapse in time and angle is If the quantity R a is not real, then there are no geodesics linking R 0 to the horizon. This occurs when |b| < 1 and the associated turning point is smaller than R 0 . Eqs. (5.11) display divergences as b → 1, corresponding to geodesics that circle the horizon many times before falling in. However, as b → 1 + (i.e. from above) the turning point moves to the horizon and there are no longer any trajectories linking R 0 > R H to the horizon (also seen by R a becoming imaginary). These b → 1 + trajectories emerge from the past horizon and orbit arbitrarily many times near the horizon radius before falling in, and

JHEP05(2020)094
are not relevant to the late-time behavior of fields from initial data confined outside the horizon, for which R 0 can be considered fixed. 5 For the relevant trajectories b → 1 − (i.e. from below) that originate from some fixed R 0 > R H , we may expand to obtain (dropping the label "nt") To leading order we thus have We now relate this result to the BTZ Green function between points (V , R , Φ ) and (V, R H , Φ). Recall that the BTZ Green function consists of a sum over image charges, such that the spikes in the Green function correspond to geodesics with initial values at the image charges. These initial values are given by where n is any integer. These equations also imply ∆V = δV and ∆Φ = Φ − Φ − 2πn. Substituting in, we see that (5.14) agrees exactly with the arrival times of the n → −∞ BTZ spikes (4.20), confirming that the b → 1 − geodesics are "responsible" for the Aretakis instability. In figure 6 we plot a selection of these trajectories. Eqs. (5.11) also display weaker, logarithmic divergences as b → −1 − , corresponding to geodesics that are initially outgoing and reach a turning point at very large radius. As opposed to the b → 1 − that spend a lot of time near the black hole, these b → −1 − geodesics spend a lot of time near the boundary. Expanding eqs. (5.11a) and (5.11b) (and dropping the label "t"), we find Thus to leading order we have Making the substitutions ∆V = δV and ∆Φ = Φ − Φ − 2πn (see eq. (5.15)) now with n > 0 to make ∆V positive, eq. (5.17) agrees precisely with the arrival times (4.19) of the 5 If we scale R0 → RH at the same rate as b → 1 + then we can recover the geodesics that turn arbitrarily close to the horizon. These are likely relevant to the case where initial data extends to the horizon.

JHEP05(2020)094
weak late-time spikes n → ∞ in the BTZ Green function. Thus these weaker spikes are associated with the b → −1 − geodesics.
As discussed in detail in section 3.1, the Green function also contains spikes that are associated with wavefronts that bounce off the boundary. Each such spike is in effect associated with two geodesics (one outgoing and one ingoing), which must be glued together by some rule determined by solving the wave equation near the boundary. Since the Aretakis instability is visible from the above analysis of the ordinary geodesics alone, we do not study this phenomenon further.
To summarize, we have shown that geodesics originating from a fixed point R 0 > R H spend arbitrarily long time outside the horizon only near the special limits b → 1 − and b → −1 − . The first limit b → 1 − corresponds to geodesics that orbit many times near the horizon and contribute wide, important spikes to the BTZ Green function at late times. The second limit b → −1 − corresponds to geodesics that spend a long time near the timelike boundary of the spacetime and contribute narrow, unimportant spikes at late times.

JHEP05(2020)094
The retarded Green function vanishes for t < t by definition, so (A.4) reduces to an initial condition We resolve the remaining delta functions using the relevant completeness relations The retarded Green function must vanish for t < t , satisfy (A.7) at t = t , and be a homogeneous solution of (A.1) for t > t . In light of the homogeneous mode solutions (A.2)-(A.3), we may satisfy the latter two conditions by multiplying the right-hand side of (A.7) by −iωe −iω(t−t ) . Adjoining a step function Θ(t − t ) fulfills the first condition, giving simply time evolve this condition to get where we also use that fact that the integrand is odd in k and even in q. Since the Green function vanishes as we approach the boundary z → 0 (or z → 0), this is the Dirichlet retarded Green function as desired. Note that with our choice of Dirichlet boundary conditions, G is conformal to the retarded propagator on the upper-half plane of Minkowski space with a mirror at infinity [29] with conformal factor zz . We perform the k-integral using eq. (3.876-1) of [33], giving (A.9) The remaining q integral is resolved with eq. (6.578-8) of [33], giving the AdS 3 Green function as This is the form given in [24], correcting a couple of typos. We can simplify further using eqs. (14.5.17), (14.3.10), and (14.5.11) in [34] to obtain This form explicitly shows the inverse square root singularity at the wavefronts Σ = 0, 2.

B Massless axisymmetric perturbations: Aretakis' method
Using the technique employed originally by Aretakis for extreme black holes in four dimensions [1,2], in this appendix we show that massless linear perturbations arising from axisymmetric data with support extending to the horizon conserve a charge on the extremal BTZ horizon. As in the four dimensional case, the associated conservation law prevents derivative decay, and higher-order derivatives grow polynomially. Assuming decay of the perturbation itself, we obtain rates for the derivatives by a hierarchy argument.
To begin we express the massless wave operator in "radially inverted" ingoing cororating coordinates (V, ρ, Φ), where ρ = R H /R, Here, f is a 2nd-order differential operator having smooth coefficients in a neighborhood of ρ = 1 (the horizon) whose precise form is irrelevant for the arguments to follow. Integrating the wave equation on a cylindrical section of the horizon ρ = 1 extending from an initial time V 1 to some final time V 2 leads to the conservation law where we have used that S 1 ∂ Φ ψ = 0 for sufficiently smooth ψ. The conservation law is stated locally as The conservation law implies that if ψ decays to zero on the horizon as V → ∞, the radial derivative must asymptote to a constant. Compactly, Furthermore, if both ψ and its tangential V derivatives along the horizon decay asymptotically as V → ∞, then higher-order radial derivatives grow at rates where c n is an order one constant depending on the specific choice initial data through H.

JHEP05(2020)094
To demonstrate the blow-up of the first derivative, take the radial derivative of the wave equation ∂ ρ 2ψ = 0, pull back to the horizon, and average over the circular cross section: Using our decay assumptions and results above for the horizon rates ψ → 0, ∂ V ψ → 0, and ∂ ρ ψ → −H/2 as V → ∞, integration of (B.7) gives linear growth The derivation of the rates for higher derivatives proceeds analogously, and straightforward induction yields (B.6).
C Mode approach

C.1 Preliminaries
The starting point for our mode calculations is the introduction of corotating, constant proper volume coordinates given by In these coordinates the separated mode solutions resemble those found for perturbations of AdS 2 with an electic field at small frequencies, whose properties have been exhaustively analysed in previous studies of the Aretakis instability [5]. Decomposing into modes, the field takes the form of a Fourier-Laplace integral For Dirichlet data, the field modes ψ mω are determined by a radial convolution integral containing mode decomposed initial data and an integral kernel called the "transfer function". The range of the integral spans the support of the data, which we choose to be bounded away from the event horizon and infinity. Since we are considering generic modeevolution of this data, we choose to work directly with the transfer function, denoted here as g m (y, y ; ω). It satisfies the inhomogeneous equation where prime denotes ordinary y-differentiation. Here again, h is as defined in (4.8).
The proceeding arguments may be understood directly from the BTZ wave equation and Kirchhoff representation of the mode decomposed retarded propagator

JHEP05(2020)094
To ensure causal evolution, we choose the transfer function corresponding to the retarded solution. This choice dictates our selection of homogeneous radial functions subject to boundary conditions corresponding to ingoing waves at the horizon R in and Dirichlet falloff infinity R D ∼ y h as y → 0. These solutions exist provided µ 2 > −1 [27]. Imposing a C 0 match of the two homogeneous solutions at the support of the delta function allows the transfer function to be written as where y < = min(y, y ) and y > = max(y, y ). Here we have introduced the conserved Wronskian of the two solutions W [R in , We express the homogeneous solutions in terms of Whittaker functions as These solutions are linearly independent provided ν is non-integer. This special case will not be considered here. Using eq. (13.14.27) of ref. [34] for the Wronskian of these two functions, we find that (C.5) can be written as The coordinates (T, y, η) do not extend to the horizon of the black hole, and thus fail to characterize Aretakis' instability. For a suitable extension to R H we now readopt the "horizon coordinates" V and Φ (see eq. (2.9)) and invert y, Explicitly, R * and Φ are given in terms of Y by In these coordinates the horizon sits at Y = 0. On the horizon, ∂ V is tangent and Φ labels generators. While at infinity, V = T and Φ = η. From the T − T and η − η dependence in the phases e −iω(T −T ) and e im(η−η ) of the mode decomposition (C.4), we see that, under the coordinate transormation (C.9), the radial Green function g mω is transformed by Both the primed and unprimed phases may be determined from the expression by complex conjugating at the primed point.

C.2 Modes at late times near the horizon
To obtain the late-time behavior of a fixed m mode, we apply the asymptotic theory of Fourier-Laplace transforms [35] wherein the functional behavior as T → ∞ is determined by the leading-order non-analytic term in an asymptotic expansion of the Laplace transform about its uppermost singular point in the complex ω plane.
Inspection of (C.8) reveals that the transfer function has its uppermost singular point at ω = 0. To determine the late-behavior near the horizon, we expand g mω in an asymptotic series about ω = 0, keeping ω/Y held fixed. This limit motivates the definition of the nearfar late-time transfer function fixing ω/Y, (C. 13) which excludes all terms analytic in ω with the exception of those involving ω/Y . In the late-time limit the M function simplifies via eq. (13.14.14) of ref. [34], and we find (C.14) Here we have defined To perform the inverse transform, we use the identity eq. (9) in section 5.20, of ref. [36]. With this, the late-time result, stated terms of the shifted time coordinate introduced previously, is (C.17) We conclude this section with a few remarks regarding our result for g nf,lt m . On the horizon, the decay of a fixed m-mode is given by δV −h . Higher transverse derivatives ∂ n R | H exhibit the Aretakis instability, growing at rate δV −h+n . These decay and instability rates are consistent with those determined from the sum over images using the method of images above in the main text. However, at large m, the magnitude of g nf,lt m grows as m 2h−1 , indicating that the mode sum does not converge pointwise (recall that h ≥ 1/2). The high-m regime of the transfer function probes the roughest part of the initial data's Φ-dependence. If the data is sufficiently regular, with differentiability exceeding 2h − 1, then its Fourier coefficients will decay in m at a faster rate than the m-growth exhibited by the coefficients in the transfer function, such that formal application of the Kirchhoff integral (3.3) provides a finite result for the field. However, the field will then be less regular than the initial data, in contradiction with rigorous results of Warnick [37]. Thus our results for each m-mode cannot be straightforwardly promoted to results about the whole field, indicating some nonuniformity in the late-time and large-m limits. Nevertheless, the rates predicted by the mode expansion do agree with those of the image sum. We hope to explore these issues further in future work.

JHEP05(2020)094
C.3 Non-periodic limit -AdS 3 The late-time behavior (C.17) of each angular mode displays growth of transverse derivatives, consistent with the full extremal BTZ spacetime possessing the Aretakis instability. However, the angular modes also apply to the AdS 3 spacetime, which has no such instability. From the perspective of the mode approach, the only difference between AdS 3 and BTZ is that the latter admits only a subset of modes compatible with the periodicity of Φ (i.e., m is continuous for AdS 3 and quantized for BTZ). That is, formally speaking we have We have already noted that the mode sum (C.18) does not converge due to g nf,lt m behaving as m 2h−1 at large m, and the integral (C.19) is likewise divergent. However, the following formal manipulation assigns it a value that agrees with AdS 3 expectations. Using eq. (C.17) and expressing the integral (C. 19) in terms of q 0 defined in eq. (C.15), we have Assuming the generic case 2h / ∈ Z + , the integrand contains poles at q 0 = −i(h + n) with n ∈ Z + . By Cauchy's theorem we may express the integral as a sum over these poles, together with a contribution from a semicircular arc at large radius in the lower-half plane. This arc contribution does not vanish, and in fact it is infinite on account of the behavior q 2h−1 0 of the integrand at large |q 0 |. Ignoring this infinite contribution provides a natural regulator that leaves a finite result for the integral, namely the residue sum. This sum can be done in closed form using the formula The result for G nf,lt AdS 3 , regulated by dropping the arc contribution, is then The dependence on Y δV in the mode result (C.17) has now disappeared, and correspondingly eq. (C.22) shows no Aretakis instability (all transverse derivatives decay). In fact, the second term in square brackets is subleading in the limit we consider (δV → ∞ fixing Y δV ), 6 so that to leading order we have simply

JHEP05(2020)094
Although there are many unresolved issues with convergence and regulators, the discussion in this appendix supports the following general picture: each angular mode in BTZ or AdS 3 displays the Aretakis instability, but these are only promoted to an instability arising from compactly supported initial data in the BTZ spacetime, where a subset of modes is selected out by the periodic identification. This is consistent with the simple idea that a single angular mode in BTZ has compact spatial support, while a corresponding mode in the AdS 3 spacetime does not.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.