Fermionic pole-skipping in holography

We examine thermal Green’s functions of fermionic operators in quantum field theories with gravity duals. The calculations are performed on the gravity side using ingoing Eddington-Finkelstein coordinates. We find that at negative imaginary Matsubara frequencies and special values of the wavenumber, there are multiple solutions to the bulk equations of motion that are ingoing at the horizon and thus the boundary Green’s function is not uniquely defined. At these points in Fourier space a line of poles and a line of zeros of the correlator intersect. We analyze these ‘pole-skipping’ points in three-dimensional asymptotically anti-de Sitter spacetimes where exact Green’s functions are known. We then generalize the procedure to higher-dimensional spacetimes and derive the generic form the boundary correlator takes near the pole-skipping points in momentum space. We also discuss the special case of a fermion with half-integer mass in the BTZ background. We discuss the implications and possible generalizations of the results.


Introduction
Despite immense progress in our understanding of quantum field theories, a complete description of strongly interacting theories is still lacking. The gauge/gravity correspondence [1][2][3] opened up a new path towards studying certain strongly coupled large-N quantum field theories by investigating their dual, weakly coupled, gravitational theories on curved backgrounds.
A basic quantity of interest in finite temperature quantum field theories is the retarded two-point function of an operator. It measures how the system in equilibrium responds to perturbations. The prescription of how to compute the correlators in holographic theories in real-time was formulated in [4] (see also [5][6][7][8][9][10][11][12][13][14]). A thermal state in the field theory corresponds to a black hole in an asymptotically anti-de Sitter (AdS) spacetime on the gravity side. Boundary operators are dual to fields in the bulk (i.e. the curved background). The AdS/CFT dictionary relates the Green's function of a boundary operator O, to solving the equations of motion for the corresponding bulk field φ. Near the black hole event horizon the second-order equation of motion has an ingoing and an outgoing solution. In order to calculate the retarded (advanced) Green's function, one should pick the ingoing (outgoing) solution [4]. This wavefunction is then evolved in the radial direction outwards to the spatial boundary of AdS where the Green's function can be read off. Using this prescription, the retarded Green's function is uniquely defined in terms of the bulk solution satisfying the prescribed boundary conditions. One of the important conditions for this is the uniqueness of the ingoing solution in the interior.
In principle the prescription for calculating the retarded Green's function G R (ω, k) is straightforward. However, evolving the ingoing solution to the boundary turns out to be computationally challenging. While it can be done explicitly in the simplest cases (e.g. the BTZ black hole [4,10,[15][16][17][18][19][20]), typically one has to use numerical methods to obtain the solutions. Generically, the retarded Green's function depends in a complicated way on the details of the state in the quantum field theory. Simplifications occur in the low-frequency and low-wavenumber limit of the correlator. In this case, the form of the retarded Green's function is dictated by near-horizon physics in the bulk and its qualitative features are independent of the rest of the geometry (see e.g. the results on shear viscosity [21]).
Recently it has been observed that certain properties of the correlators away from the ω = 0, k = 0 point in Fourier space can already be seen in the near horizon behavior of the solutions [22][23][24][25]. Initially it was observed that at special complex values of the frequency and the momentum, the retarded Green's function contained information about the chaotic behavior of the theories (see [26][27][28][29][30]). Such behavior was dubbed "pole-skipping" as it occurs where a line of poles intersects a line of zeros in the Green's function of the dual boundary operator (see also [31] for related phenomena in the case of Fermi surfaces).
As it was shown in [20,[32][33][34] pole-skipping is not limited only to the components of the energy-momentum tensor but can also be observed in other fields in the theory. Its gravitational origin stems from the fact that at these points there is no unique ingoing solution at the interior of the bulk spacetime. With this, the holographic retarded Green's function ceases to be uniquely defined and becomes multivalued, taking on a spe-JHEP07(2020)203 cial pole-skipping form [23,24]. The interesting aspect of this phenomenon is that the bulk computation is limited to the horizon and has no knowledge of the boundary. In this sense a local calculation in the bulk constrains the structure of boundary Green's functions.
In this work we build on the findings of [20] and describe the pole-skipping for minimally coupled spinor fields on asymptotically AdS backgrounds. Looking at the exact Green's function for fermions in the BTZ black hole background found in [10], one can observe that there are special points at which the lines of poles intersect the lines of zeros. They occur precisely at the fermionic Matsubara frequencies 1 ω = ω F n := −2πiT n + 1 2 , n = 0, 1, 2, . . . . (1.1) This nicely complements the fact that scalar and energy-momentum pole-skipping points occur at bosonic Matsubara frequencies ω = ω B n = −2πiT n, with n = 1, 2, 3, . . . . In light of this, it has been conjectured in [20] that this behavior has a bulk interpretation in terms of non-unique ingoing solutions. Here we will explicitly show that this is indeed the case.
We emphasize that only the energy-momentum near-horizon behavior is clearly related to chaos in holographic theories. There, one can observe pole-skipping at the first positive bosonic Matsubara frequency ω = +2πiT where the right-hand side is precisely the Lyapunov exponent that characterizes out-of-time order higher-point functions.
In other examples, such as the case of the scalar field, pole-skipping occurs on the lower-half frequency plane. Since we will encounter pole-skipping at fermionic Matsubara frequencies, it is even less likely that this phenomenon can be related to quantum chaos in a straightforward manner. However, these features might be important in holographic theories in general.
The paper is organized in the following way. In section 2 we review the pole-skipping phenomenon in the case of a minimally coupled scalar field. In section 3 we define a minimally coupled fermion field on an anti-de Sitter background and discuss spinors in holography. Then in section 4 we look at pole-skipping in 3-dimensional bulk spacetimes. The generalization to higher dimensions is given in section 5, while in section 6 we examine the form of the Green's function near a generic pole-skipping point and discuss the appearance of anomalous points. In section 7 we discuss some examples. Most notably we use the results to calculate the fermionic pole-skipping points for the BTZ black hole and compare our results with the known retarded Green's function. We examine the special cases of boundary operators with half-integer conformal dimensions and relate them to anomalous pole-skipping points. We conclude with a discussion in section 8. In appendix A we present explicit representations of gamma matrices that can be useful in practical applications. Some of the more detailed calculations omitted in the main text are collected in appendix B. In appendix C we review the calculation of the exact Green's function for the BTZ black hole. We also consider the equality of the retarded and advanced Green's function at the pole-skipping points. Finally we calculate the form of the retarded Green's function in the special cases where the mass of the fermion takes a half-integer value.

JHEP07(2020)203 2 Review of pole-skipping
In this section we present the general form of the background metric in ingoing Eddington-Finkelstein coordinates and review the systematic procedure to extract the locations of the pole-skipping points in the case of a minimally coupled scalar field, which was developed in [20].
We start by assuming that the action for the background fields is given by where Λ = −d(d + 1)/2L 2 is the cosmological constant and L is the AdS radius, which we henceforth set to L = 1. The term S matter allows for additional matter content which can also contribute to the curvature of the background. We further assume that the equations of motion for this action admit a planar black hole solution given by the metric where r is the radial direction. The boundary of spacetime is located at r → ∞. Furthermore, t denotes time and x i with i = 1 . . . d are the (flat) coordinates of the d spatial dimensions. The combination (t, x) ∈ R 1,d also denotes the Minkowski coordinates of the corresponding boundary theory. The exact form of the two functions f (r) and h(r) in general depends on the matter content of the theory. Since we want our spacetime to be asymptotically anti-de Sitter, they have to approach f (r) → 1 and h(r) → r 2 as r → ∞. We assume that the background has a horizon at r = r 0 , i.e. the emblackening factor vanishes at this radius: f (r 0 ) = 0. We also assume that the Taylor series of the functions f and h have finite radii of convergence near the horizon. The Hawking temperature of the black hole is given by In order to extract the pole-skipping points, it is convenient to introduce the ingoing Eddington-Finkelstein coordinates, defined by v = t + r * , dr * dr = 1 r 2 f (r) , (2.4) in which the background metric takes the form The vacuum solutions (S matter = 0) with such properties are characterized by which are the BTZ black hole [15,16]  JHEP07(2020)203

Minimally coupled scalar field in the bulk
The simplest instance for which one can observe pole-skipping is a minimally coupled scalar field in an asymptotically anti-de Sitter spacetime. To that end, we add to the action of the background (2.1) the action of a massive scalar in a curved background, which is given by The Green's function can be extracted by finding solutions to the equation of motion Note that this is a second order differential equation for a single scalar field. As such it has two free parameters that we need to fix with boundary conditions. The scaling dimension ∆ and the mass m of the scalar field are related via where we take the larger of the two roots to be the scaling dimension in the standard quantization.
If we wish to calculate the retarded Green's function, we need to choose the ingoing solution at the horizon [4]. To do so, we consider the ansatz ϕ = φ(r)e −iωv+i k· x and perform a series expansion of φ(r) around the horizon. We find that this boundary condition gives a unique ingoing solution to (2.8) for generic values of ω and k up to an overall normalization. The next step is to expand this solution near the boundary as to obtain the boundary retarded Green's function up to the possible existence of contact terms by (2.11)

Pole-skipping points
Here we briefly explain why imposing boundary conditions at special values of frequency ω and momentum k is not sufficient to uniquely (up to an overall factor) specify a solution ϕ to the equation (2.8) and give the locations of the pole-skipping points for a minimally coupled scalar field. We closely follow [20], where these calculations were initially performed. See that work and the references therein for more details. After performing the Fourier transform and switching to the Eddington-Finkelstein coordinate system, (2.8) becomes

JHEP07(2020)203
We look for solutions that are regular at the horizon. Such solutions can be written as a Taylor series expansion φ(r) = φ 0 + φ 1 (r − r 0 ) + . . . around r = r 0 . Near the horizon, there exist two power law solutions φ = (r − r 0 ) α with which do not depend on k and m. For generic values of ω, only the solution with exponent α 1 is regular and is therefore taken to be the ingoing solution. (Note that α 1 had to be zero, because the horizon is not a distinguished location in infalling coordinates.) However, at the special values of frequency ω n = −2iπT n with n ∈ {1, 2, . . .}, the second exponent becomes α 2 = n and naively both solutions seem to be regular at the horizon. A more detailed calculation shows that one of the solutions contains logarithmic divergences which spoil the regularity, so there is still a unique regular ingoing solution. One then finds that all such logarithmic divergences vanish for some particular values of the momentum k. This means that for finely tuned values of ω and k, there is no unique ingoing solution to (2.12), which renders G R OO (ω, k) ill-defined. To see this explicitly we expand (2.12) in a series around the horizon and solve the resulting equation order by order. At the zeroth order, one obtains a relation between the lowest two field coefficients φ 0 and φ 1 by which, for generic values of k and ω, fixes φ 1 in terms of φ 0 . Higher order terms of the series expansion of the equation of motion allow us to relate all the field coefficients φ n in terms of only φ 0 . Thus we explicitly construct a unique regular solution with an undetermined overall normalization in the form of the factor φ 0 . If the frequency takes the value of the first bosonic Matsubara frequency ω = ω 1 = −2iπT , this method fails as (2.14) reduces to For generic values of k, the above equation sets φ 0 = 0. All higher order coefficients φ n are then related to φ 1 , which can be taken as the undetermined normalization of the unique solution. However, by finely tuning both the frequency ω and the momentum k to take the values the equation (2.14) becomes trivially satisfied. In this case, both φ 0 and φ 1 remain undetermined and all higher coefficients of the series expansion of the scalar field φ n are determined in terms of both φ 0 and φ 1 . The regular solution then has two independent parameters and is thus not unique. Consequently, at (2.16), the boundary Green's function is not uniquely defined. One finds that there are pole-skipping points at higher Matsubara frequencies as well. At ω = ω n = −2πiT n, there are 2n wavenumbers k n at which we observe pole-skipping.

JHEP07(2020)203
In order to locate these points, one needs look at higher orders in the expansion of (2.12) around the horizon. Setting all the coefficients of the expansion to zero results in a coupled set of algebraic equations that can be written as where the coefficients are generically of the form M ij (ω, k 2 ) = iω a ij + k 2 b ij + c ij , with a ij , b ij , and c ij determined by the background metric.
At generic values of frequency, (2.17) is easily solved in an iterative manner. In fact, these are the equations that allow us to express all φ n as functions of φ 0 . However, at ω = ω n , it is not possible to construct an ingoing solution in this way as the coefficient of φ n vanishes in the n th row of (2.17). We then obtain a closed set of equations for the coefficientsφ = (φ 0 , . . . , φ n−1 ), which is of the form where M (n) (ω n , k 2 ) is the submatrix of M (ω, k 2 ) consisting of the first n rows and first n columns. For generic values of k, the matrix M (n) (ω n , k 2 ) is invertible, settingφ = 0. With that, φ n takes the role of the free parameter and the remaining equations in (2.17) can be used to relate φ m , with m > n, to φ n , thus obtaining a unique ingoing solution up to an overall normalization, which is now φ n . If, on the other hand, the value of k is such that the matrix M (n) (ω n , k 2 ) is not invertible, then we get an additional non-trivial ingoing solution which is parametrized by a free parameter that we can choose to be φ 0 . The regular solution has two free parameters (φ 0 and φ n ) and the boundary Green's function is again not unique. The values of k for which M (n) is not invertible are the same as the ones at which the determinant of the matrix vanishes. Pole-skipping at higher Matsubara frequencies can therefore be observed at the special locations In summary, at special points in Fourier space (2.19), imposing the ingoing boundary condition at the horizon is not enough to select a unique solution to the wave equation. It can be shown that as a consequence, near such points in momentum space, the correlator takes on the pole-skipping form [20,24], given by where (δω/δk) p is the direction in Fourier space at which we obtain a normalizable solution corresponding to the poles in the boundary Green's function, while (δω/δk) z denotes the direction of the non-normalizable solution associated with the slope of the line of zeros of JHEP07(2020)203 the correlator. This form shows that generically pole-skipping points are intersections of lines of zeros and lines of poles of the retarded Green's function. This is the motivation behind the name "pole-skipping", because at such points the expected divergence coming from the vanishing denominator of (2.20), is cancelled out by a coinciding zero of the numerator of the correlator [22][23][24][25].
There also exist anomalous points [20] which are locations that appear as possible pole-skipping points from the near-horizon analysis, but the Green's function near them does not take the pole-skipping form (2.20).

Minimally coupled fermion in the bulk
The aim of this paper is to locate the pole-skipping points for a general fermionic field in an asymptotically anti-de Sitter background. To do so, we must add to the background the action of a minimally coupled fermion field given by [35,36] where S bdy is a boundary term that does not alter the equations of motion, the fermion conjugate is defined as ψ = ψ † Γ 0 , and the covariant derivative acting on fermions is defined by with ω M denoting the spin connection one-form and Γ ab is the anti-symmetrised product of flat space gamma matrices (see (A.2)). In what follows we will denote the curved indices by upper-case Latin letters while flat space indices are denoted by lower-case Latin letters. 2 The resulting equation of motion for the spinor ψ is then the Dirac equation Recall that for a theory in d + 2 spacetime dimensions, the number of components of a spinor is given by where q denotes the highest integer that is less than or equal to q. This makes the Dirac equation (3.3) a system of coupled first order differential equations for the N components of the spinor. To fully specify the solution we thus need to impose N boundary conditions. To calculate the retarded Green's functions for spinors we follow the prescription given by [10]. We first introduce the decomposition of the spinor in terms of the eigenvectors of the matrix Γ r defined by

5)
2 A further comment on notation. A general flat space tensor has lower-case Latin letter indices, but particular values for the indices are underlined, for example v, r, or x. This is to distinguish them from curved space indices where a generic tensor has upper-case Latin letters, but a particular value is lower-case letter that is not underlined, for example u, v, or x.

JHEP07(2020)203
where ψ ± each contain N/2 degrees of freedom. Assuming that the metric components only depend on the r coordinate, we make the plane wave ansatz ψ = ψ(r)e −iωt+i k· x and solve the Dirac equation in Fourier space. If we want to calculate the retarded Green's function, we need to choose the solution that is ingoing at the horizon. This boundary condition usually reduces the number of free parameters in the solution to N/2. We then evolve the solution to the AdS boundary (r → ∞), where we find that in general it takes the following form 3 with the Dirac equation imposing relations between the pairs B(k), D(k) and A(k), C(k).
For m ≥ 0 the dominant contribution comes from the term multiplied by A(k), which thus is identified with the source. The response is given by D(k) as it is related to the finite term in the conjugate momentum to the field ψ + in the appropriate limit. With this identification the mass m of the spinor in the bulk and the conformal dimension ∆ of its corresponding response in the boundary spinor are related via 4 The prefactors A(k) and D(k) are spinors, and after imposing the ingoing condition one can find that they are related by a matrix R(k) as The retarded Green's function in the boundary theory is given by It might be worth stressing how choosing the ingoing solution at the horizon renders the retarded Green's function unique for both scalar and fermion fields. A scalar field has only one component, but since its dynamics is governed by a second order differential equation, we need two boundary conditions to fully determine the solution. The ingoing condition at the horizon imposes one constraint and thus the solution is effectively determined up to an overall normalization. As the correlator is a ratio between the two leading terms in the asymptotic expansion (2.10), this overall normalization cancels out and the Green's function is thus uniquely defined.
For a spinor field the procedure is conceptually the same as one can see the matrix R(k) as a generalized ratio between two terms in the asymptotic expansion. However, the calculations are more involved. The ingoing solution at the horizon fixes half of the degrees of freedom. This is usually achieved by transforming the Dirac equations into a second JHEP07(2020)203 order equation for half of the components, say ψ + , and then taking the ingoing solution.
Putting the ingoing solution into the first order Dirac equation fixes the other half of the components, in this case ψ − , in terms of the free parameters left in ψ + . Therefore the solution is completely determined up to an overall spinor with N/2 free parameters that multiplies both ψ ± . When the solution is then evolved and expanded near the boundary, both D(k) and A(k) are proportional to this overall spinor, albeit the factor of proportionality can be a matrix in spinor space. This means that R(k) does not depend on any free parameters and therefore the retarded Green's function is uniquely defined.

Pole-skipping in asymptotically AdS 3 spaces
We start with the simplest low-dimensional example, where the bulk theory is threedimensional and the boundary theory has two spacetime dimensions. In this case both bulk and boundary spinors have two components. We will observe pole-skipping and develop a systematic approach to extract the location of the points in Fourier space for any three-dimensional background.
Let the background metric be given by where for now we leave f (r) and h(r) unspecified, apart from the properties described in section 2. Let us choose the following frame for which We choose this frame firstly because neither the vielbein components nor any of their derivatives diverge at the horizon (assuming h(r) is regular at r = r 0 ). Secondly, we avoid any square roots of the emblackening factor f (r) in the equations. Furthermore, this vielbein reduces to a frame for AdS 3 at the leading order in the near-boundary limit r → ∞. In this frame, the spin connections are given by with all other components, which are not related by symmetry to the ones above, vanishing. In this frame the Dirac equation is given by (4.5)

JHEP07(2020)203
Since the metric is independent of the coordinates v and x i , we introduce the plane wave ansatz ψ(r, v, x) = ψ(r)e −iωv+i k. x . Furthermore, we separate the spinors according to their eigenvalues of the Γ r matrix. We define the two independent spinor components associated with these eigenvalues as The spinors ψ ± are two component objects, but contain only one independent degree of freedom each. We insert this decomposition into (4.5) and act on the equation with the projection operators defined in (4.6). After some algebra one can write the two resulting equations as Above we have used the fact that the set of matrices (1, Γ v , Γ x , Γ r ) forms a complete basis for all 2 × 2 matrices, hence Γ vx can be rewritten as a linear combination of the matrices from the set. In fact, Γ vx = ±Γ r and we choose a representation such that Γ vx = Γ r . For more details on gamma matrices and explicit representations, see appendix A. It is straightforward to transform (4.7) into two decoupled second order ordinary differential equations for the spinors ψ ± . Using these second order differential equations one can look for the leading behavior of the spinors at the horizon. In practice this is achieved by introducing an ansatz where ξ + is a constant spinor satisfying Γ r ξ + = ξ + , and expanding the second order differential equations around the horizon r = r 0 . One then finds that the equations are solved at first order for the exponents 5 One can repeat the procedure for the ψ − spinor and obtain the same exponents as in the case of ψ + . Recall that in order to obtain the retarded Green's function, we are supposed JHEP07(2020)203 to select the ingoing solution at the horizon and evolve the solution towards the boundary. In ingoing Eddington-Finkelstein coordinates, this translates to taking the solution with α 1 . However, naively both solutions are ingoing if ω is such that α 2 is a positive integer which happens at ω = ω n ≡ −2πiT n + 1 2 , n = 1, 2, 3, . . . . (4.10) These are precisely the fermionic Matsubara frequencies (1.1), with the exception of the lowest frequency ω = ω 0 ≡ −iπT , which appears to be missing. Choosing such frequencies is not enough to produce two independent ingoing solutions. Similar to the scalar field, a more thorough analysis shows that logarithmic divergences appear in the expansions, making one of the solutions irregular. If, in addition, we also tune the momentum k to values such that these logarithmic divergences vanish, then there will be two independent ingoing solutions at the horizon. In this case the corresponding Green's function will show pole-skipping, as the ingoing solution and therefore the Green's function is not unique.

Pole-skipping at the lowest Matsubara frequency
In the case of the minimally coupled scalar field, the lowest Matsubara frequency is given by ω = 0. No pole-skipping has been observed at this frequency [20]. For the fermionic field, the lowest Matsubara frequency is given by ω 0 = −iπT . The exponents (4.9) suggest that there is no pole-skipping at this frequency. However, this is not the case as we will soon see.
Pole-skipping at the lowest frequency occurs if there exist two independent ingoing solutions that behave as (r − r 0 ) 0 at the horizon. For the scalar field this actually implies that the two independent solutions are of the form with C and D being the free parameters associated with the two independent solutions. φ i (k) are coefficients fixed by the equation of motion. Unlike for any other bosonic Matsubara frequency ω n = −2πiT n , n ∈ Z + , we cannot choose any value for k that would give a vanishing prefactor multiplying the logarithmic term. The upshot of this is that for α = 0, there is only one solution that is regular at the horizon, and thus no pole-skipping can be observed at this frequency. The spinor in d ≥ 1 is a multicomponent object which allows for pole-skipping to occur at the lowest Matsubara frequency. Let us introduce a series expansion for both spinor components ± are constant spinors with definite Γ r eigenvalues. We put these expansions into (4.7) and expand the equations in a series around the horizon as In the above definitions, S + and S − are the horizon expansions of the equations (4.7a) and (4.7b) respectively and S (j) ± are series coefficients that can in principle depend on both ω and k. This dependence will be suppressed in the following.
We solve the equations (4.13) order by order. For the first instance of pole-skipping we only need to look at zeroth order coefficients. These are We can immediately notice that and thus equations (4.14) actually represent only a single constraint. This is not surprising, as the zeroth order should fix one of the components in terms of the other so that we get a unique ingoing solution, up to an overall constant. If (4.14a) and (4.14b) were two independent equations they would completely fix ψ (0) ± leaving it with no free parameters. To locate the pole-skipping points, we need the scalar coefficients multiplying ψ (0) ± to vanish. This happens precisely at which is precisely the zeroth fermionic Matsubara frequency and the associated momentum. Here we have used the definition of the Hawking temperature (2.3). At such points, equations (4.14) are automatically satisfied and thus ψ (0) ± both remain free and independent coefficients.
One can then take a look at the equations at higher orders in (4.13). These relate the expansion coefficients ψ  − . Therefore, at (4.15), the retarded Green's function is not uniquely defined. 6 Dealing with logarithmic divergences. Finally, one may ask what happens to the logarithmic terms that one observes in the scalar field expansion at ω = 0. As can be shown, such divergences also appear in the fermion field expansion and are the reason why for generic values of the momentum we do not find two independent ingoing solutions. This highlights the fact that one needs to tune both the frequency and momentum to obtain two non-divergent ingoing solutions, even in the fermionic case.
To see this explicitly, we are interested in the near-horizon solutions to Dirac equations at the frequency ω = ω 0 = −πiT . To leading order, the solutions to the equations take the following form ± being constant spinors of definite chirality. We insert the expansion into (4.7) and expand the equations in a series around the horizon. The equations now take the form and we solve them iteratively. At leading order we get the following 4 equations The last two are not independent and are related via In addition to that, inserting χ (0) − , which is the solution of (4.18c), into (4.18a) also gives (4.20) 6 Note that the solution parametrized by, for example, ψ ± are already both non-vanishing. So, the leading component in the expansion has a well defined eigenvalue under Γ r , but, as soon as we move away from the horizon the two components will start to mix. The same is true for ψ

JHEP07(2020)203
This means that for a generic value of k there are only two independent equations in (4.18) and there exist solutions with χ (0) ± = 0. We have to set these coefficients to zero if we want a regular solution at the horizon. Thus for a general value of k there is still a unique ingoing solution, as the other solution contains logarithmic divergences.
If we set k to (4.15), then (4.18c) and (4.18d) are automatically satisfied. Furthermore, the remaining two equations (4.18a) and (4.18b) are now independent and in fact the first terms in both equations vanish. The solution to (4.18) is then given by χ ± being undetermined. Thus we see explicitly that at the location of the pole-skipping point (4.15), the logarithmic terms vanish and the two independent solutions are both regular at the horizon.

Pole-skipping at higher Matsubara frequencies
There are two equivalent ways to locate the pole-skipping points associated with higher Matsubara frequencies. The first method is similar to the procedure used for the scalar field, as one uses the second order differential equations for half of the components. This method is useful to determine the positions of the so-called anomalous points, which are the locations of coinciding pole-skipping points.
The second method is inspired by the lowest frequency pole-skipping point and uses only the first order Dirac equation. This method completely bypasses the computational difficulties of obtaining a decoupled second order equation, however at the expense of working with higher dimensional systems of algebraic equations.
One can show that both methods yield the same results and we will show in section 7 that they exactly locate the points of intersection between the lines of poles and the lines of zeros for Green's function in a BTZ black hole background. There we will also illustrate the use of the procedure using the first order Dirac equation.
Here we present both methods in turn. In both cases, we initially look at the lowest frequency pole-skipping location before generalising the procedure for arbitrary frequencies.

Using the second-order differential equations
The first method mimics the procedure of the scalar field reviewed in section 2. As mentioned above one can use (4.7) to obtain decoupled second order differential equations for the components of one of the spinors. Without loss of generality, we work with ψ + . The first order Dirac equations then completely determine the components of ψ − in terms of ψ + .
We begin by expanding the second order differential equation of ψ + around the horizon. This can be schematically written as where D (j) + can in principle depend on both ω and k. These terms also depend on the expansion coefficients of ψ + defined in (4.12). We solve (4.21) perturbatively by solving D + . In this manner one explicitly constructs an ingoing solution which is unique up to an overall spinor and whose leading behavior at the horizon is (r − r 0 ) 0 .
The above procedure fails if the frequency matches the first fermionic Matsubara frequency given by 23) as in this case the coefficient of ψ + is left undetermined. One can use the latter as the free parameter and again explicitly construct a regular solution that is determined up to an overall factor, ψ + . The leading behavior at the horizon of such a solution is (r − r 0 ), as the (r − r 0 ) 0 solution includes logarithmic divergences, as discussed in [20].
However, if, in addition to ω = ω 1 , the momentum k is such that then the equation (4.22) is automatically satisfied and both ψ (1) + . We hence construct two distinct regular solutions at the horizon, one with leading behavior (r − r 0 ) 0 and one with (r − r 0 ). Consequently, the retarded Green's function is not unique. Note that in general (4.24) is a third order polynomial in k and thus one expects three complex solutions for k.
One can go further in the series (4.21). At the linear order in the expansion coefficient one gets At generic values of ω and k, (4.25) combined with (4.22) fix ψ (1) + and ψ (2) in terms of ψ (0) + and one can repeat the general procedure of obtaining a unique ingoing solution, as discussed above.
At the second fermionic Matsubara frequency vanishes. Namely, at such points, equations (4.22) and (4.25) are not independent and allow us to express ψ + , meaning that we again have two independent regular solutions at the horizon.
The procedure for finding the locations of pole-skipping points associated to higher Matsubara frequencies is easily generalized. The equation (4.21) at order (n − 1) is The pole-skipping point is obtained when the coefficient multiplying ψ where M (n) Then ψ (0) + and ψ (n) + are the two independent free parameters that can be used to explicitly construct the two regular solutions at the horizon. As det M (n) + (ω n , k) is in general an (2n + 1)-degree polynomial, we can expect the same number of complex roots and thus (2n + 1) pole-skipping locations associated to the frequency ω = ω n .
So far we have not specified the representation for the gamma matrices. Thus ψ + and all ψ (k) + are two-component objects. Therefore, all entries in (4.30) are 2 × 2 matrices. However, the second order differential equations for ψ + are diagonal and consequently the entries of (4.30) are proportional to two-dimensional identity matrices. The determinant (4.29) can then be calculated as if the coefficients were scalars. This is not surprising. If we choose a gamma matrix representation in which the Γ r matrix is diagonal, the equations (4.7) reduce to scalar equations and all the entries in (4.30) become scalar functions as well.

Using the first-order equations
One can obtain pole-skipping points at higher frequencies directly from the first order equations (4.7) without having to transform them into second order equations. Not only does this method provide an alternative to the previously mentioned one, but it is also the direct generalization of the method used to find the pole-skipping point at the lowest Matsubara frequency. Using this method one can thus find all pole-skipping points for the fermionic field.
We previously looked at the series expansion (4.13) at zeroth order where we found the first pole-skipping point (4.15). To obtain the locations with higher frequencies we look at the higher order terms in the expansions of the Dirac equations around the horizon. We begin by looking at the linear terms. The two equations at this order can be written in a matrix form as where M are 2 × 2 matrices whose elements are commuting 2 × 2 matrices. 8 For example while M (10) depends on k but is independent of ω. Its explicit form is not very illuminating, so we do not present it here. As S − , there are two independent equations at linear order in the series expansion of (4.13). This is expected, as for generic values of ω and k these equations fully determine ψ (1) ± in terms of the coefficient left undetermined in (4.14). By repeating the procedure at higher orders we explicitly construct a solution that is regular at the horizon and determined up to an overall factor that contains half a spinor's worth of free parameters.
The above procedure fails if the equations (4.31) do not provide two independent constraints on ψ (1) ± . This is the case if one cannot rearrange (4.31) to express (ψ (1) (11) is not invertible. Thus we are looking for values of the frequency at which the determinant of the matrix multiplying (ψ (1) which vanishes precisely at ± . In this case the combination ψ (1) c constrained by the equations is given by Combining (4.14) and (4.31) evaluated at ω = ω 1 thus yields a system of three independent equations for three variables, which can be schematically written as (4.36) The elements of the matrix M 1 are the appropriate coefficients from the equations (4.14) and (4.31) evaluated at the first Matsubara frequency. As such they are still commuting matrices, and their k dependence has been suppressed. Elements M (11) ± are given by and thus we again obtain a unique ingoing solution that has half a spinor's worth of free parameters. With that, one can take, for example ψ + as the free parameter and then use the higher order equations to determine other coefficients and thus perturbatively construct a regular solution.
One obtains two independent regular solutions if the matrix M 1 is not invertible. In that case not all three equations in (4.36) are independent. One obtains another free parameter, for example ψ  As a consistency check, the above equation yields the same roots as (4.24). The determinant is a cubic function of k so we expect three complex roots. These are the pole-skipping points associated with the frequency (4.23). As a side note, to find the locations of the pole-skipping points, in practice it is easier to simply set one of the ψ

JHEP07(2020)203
Pole-skipping points associated to higher Matsubara frequencies are located in a similar manner. We take the equations at order n in the expansion (4.13) and write them schematically as  with all M (jk) being matrices whose elements are commuting matrices. Only the leading coefficient M (nn) depends on both the frequency and momentum, while the remaining coefficients depend only on the momentum. To get pole-skipping at ω = ω n we require that the equations (4.40) provide only one independent constraint for ψ One finds that for any n, the matrix M (nn) has the form This vanishes at the fermionic Matsubara frequencies given by The corresponding momenta at which pole-skipping occurs are then found by constructing the analogue of the equation (4.36). We start by evaluating (4.40) at (4.44). Again, only a particular linear combination of the n-th order coefficients is constrained by the equations and is given by One then combines all the equations at lower orders and evaluates them at the Matsubara frequency. These can be written in a schematic form as  We note that each M entry is a linear function in k and (4.46) is a system of 2n + 1 equations, meaning that the determinant is an order 2n + 1 polynomial in k and has that many complex roots and thus for each frequency ω n , we find 2n + 1 pole-skipping points.

Pole-skipping in higher dimensions
We now generalize the procedure presented in the previous section to higher dimensional spacetimes. We work in d + 2 bulk spacetime dimensions, which means that the boundary theory is formulated in d + 1 dimensions. Thus, the bulk spinor has N = 2 d+2 2 degrees of freedom and the boundary spinor has half as many.
We find that the equations split up into two decoupled subsystems both of which are related to the lower-dimensional case presented in the previous section. We also find that for generic values of k the number of pole-skipping points at ω = ω n is doubled to 2(2n + 1) and that the locations are in general different for the two different subsystems.
We work with the background metric in ingoing Eddington-Finkelstein coordinates (2.5). The orthonormal frame is taken to be so that 1, 1, . . . , 1) .
This frame is the direct generalization of the frame (4.2) and shares all of its special properties. The spin connections are given by with all other components not related by symmetry to the ones above being 0. The calculation of the Dirac equation is conceptually the same as in section 4, so we don't repeat it in full. We exploit again the fact that the metric does not depend on v and x i and solve the equation in Fourier space by introducing ψ(r, v, x j ) = ψ(r)e −iωv+ik i x i . The Dirac equation then reads

JHEP07(2020)203
In general, this is a system of N first order coupled ordinary differential equations for the components of the spinor. We want to decouple them in a way that makes the pole-skipping mechanism manifest. We begin by introducing the decomposition where each component ψ ± contains N/2 free parameters. Furthermore, notice that for d ≥ 2 the two matrices Γ r and k i Γ vi are independent and commuting. 9 Therefore, we can introduce an additional decomposition where a = ±, and we have used We now have divided the initial spinor ψ that with N degrees of freedom into four independent spinors ψ (±) ± which each contain N/4 independent degrees of freedom. Each ψ (±) ± has a set of definite eigenvalues under the action of Γ r and k i Γ vi . The Γ r matrix projects the spinor components along the radial direction and can be considered as the chirality projection, especially with respect to the boundary theory. We thus refer to ψ ± as positive or negative chirality spinors. Similarly k i Γ vi can be considered as a projection of the components of the spinor along the direction of the momentum. In that way, it has a similar effect as a helicity projection. We thus refer to spinors ψ (±) as positive (negative) helicity spinors. As an example ψ (−) + is a spinor with positive chirality but negative helicity. Using this decomposition into 4 independent components in the Dirac equation, one notices that the equations separate into two decoupled subsystems. One for (ψ (5.7b) 9 For d = 1 which is the asymptotically AdS3 case, we have Γ vi = ±Γ r , as (1, Γ v , Γ i , Γ r ) provide a complete basis for any 2 × 2 matrix.

JHEP07(2020)203
The equations for (ψ First, one can eliminate one of the spinors from (5.7) to obtain a diagonal second order differential equation for the other. One can expand the second order equations around the horizon region. Using the ansatz ψ (+) where ξ (+) + is a constant spinor with definite chirality and helicity, one finds that the second order equations are solved at leading order by the following exponents with the same behavior being observed for ψ (−) − . We recall that the same exponents have been found in lower dimensional case (4.9). The exponent α 1 is for generic values of the frequency associated with the ingoing solution. The exception are the cases where the frequency is such that the second exponent is equal to a positive integer. This is where we expect pole-skipping. These special values for the frequency are again the fermionic Matsubara frequencies ω n = −2πiT n + 1 2 , n = 1, 2, 3, . . . , (5.10) which again does not include the zeroth Matsubara frequency ω 0 = −iπT . For pole-skipping to occur the momentum also needs to be set to special values. In the following, we briefly discuss how to locate the pole-skipping points in general dimensions. Since the equations governing the spinors are essentially the same as in the threedimensional case, the procedure of finding the locations is also the same. In order not to repeat too much of section 4, we only outline the procedure and focus mainly on how to obtain the locations of the pole-skipping points.

Pole-skipping at the lowest Matsubara frequency
We begin by expanding the spinors in a series around the horizon as The detailed derivation of these equations together with the explicit equations for (ψ

JHEP07(2020)203
We insert these expressions into the Dirac equations (5.7) and expand them around the horizon. Schematically these expansions are written as In order to see pole-skipping at the lowest frequency, we look at the zeroth order equations in (5.12) and find that there is only one independent equation at this order. It can be written as Hence, combining these two results, one sees that for d ≥ 2, the first occurrence of poleskipping is at Comparing the above results to those in asymptotically AdS 3 spaces we see that in higher dimensions there exists an additional pole-skipping point with negative imaginary momentum. This feature repeats itself with all other pole-skipping points. Fermions in higher dimensions have twice as many pole-skipping points than fermions in 3-dimensional spacetimes. These additional locations appear due to the interaction between fermions whose chiralities are opposite to their helicities. In two dimensions such fermions are absent which explains why we observe only half as many pole-skipping points as in the general case.
Finally, each of ψ does not hold automatically at this pole-skipping point. The poleskipping point associated to this subsystem has the opposite value of k, meaning that at zeroth order we will get a constraint equation for these two coefficients. This means that although we are at a pole-skipping point, the Dirac equations still provide some constraints on the spinor, and the ingoing solution will thus have only 3N 4 free parameters. There is, however, a notable exception to this rule -the case of the massless fermion. Taking m = 0 in (5.16), we notice that the two pole-skipping points merge into one, located at ω = ω 0 and k = 0. At this point in momentum space, the ingoing condition does not impose any constraints on the spinors. This is unlike the scalar case, where at any pole-skipping point the ingoing condition does not impose any constraints on the field regardless of the mass of the field.

Higher order pole-skipping
To get the higher frequency pole-skipping points one can either use the first order or the diagonal second order differential equations, as they give the same locations. The equations in higher dimensions (5.7) have the same form as the ones in three spacetime dimensions (4.7). Therefore both methods readily generalize to higher dimensions.
Due to this similarity, we do not repeat the methods here. We only mention some of the differences. The first difference is that in higher dimensions, spinors have N components and separate into spinors containing N/4 degrees of freedom ψ (±) ± . With that in mind, all the factors multiplying ψ (±) ± in the expansions around the horizon are N × N dimensional matrices. Thus the analogues of (4.30), (4.42) and (4.46) will be matrices whose elements are (commuting) N × N dimensional matrices.
The second difference is that the equations split into two independent subsystems which in general yield two independent sets of pole-skipping points. However, it is enough to find the locations for one of the subsystems as the pole-skipping points of the other are obtained by k ↔ −k, meaning that at the same frequency the two subsystems have opposite pole-skipping points.
Doing the explicit calculations, we find that all pole-skipping points occur at the fermionic Matsubara frequencies (1.1) regardless of the dimension of spacetime. At each frequency ω = ω n we get, for generic values of the mass, 2(2n + 1) pole-skipping points.
As in the case with the lowest Matsubara frequency, at a generic pole-skipping point, only N/4 components of the spinor are constrained by the equations. Again this is related to the two subsystems experiencing pole-skipping at different locations in momentum space. However, explicitly working out the locations for the first few pole-skipping points, one again notices that the massless fermion is an exception. In that case, one finds that the pole-skipping points associated to the n-th Matsubara frequency for one of the subsystems are given schematically at ω = ω n and k = {0, ±k 1 , ±k 2 . . . , ±k n }. As we can see, this set of pole-skipping points is invariant under the reversal of the momentum and therefore the other subsystem will experience pole-skipping at the exact same locations in momentum JHEP07(2020)203 space. Thus, the number of pole-skipping points is halved to (2n + 1), yet, at each poleskipping point we are left with an entire spinor's worth of free parameters.
While we are currently lacking a proof that this pattern continues for arbitrary n, we find no reason why this feature would cease to hold after the first few pole-skipping points, for which this was checked explicitly.

Green's function near the pole-skipping points
In this section we show that the near-horizon equations predict that in the neighbourhood of the pole-skipping points in momentum space the fermionic Green's function takes the pole-skipping form This signals that at such locations a line of poles and a line of zeros of the boundary correlator intersect. We perform the calculations in an asymptotically AdS 3 spacetime, as the results of the previous section suggest that the analysis in higher dimensions is analogous.

Near the lowest Matsubara frequency
At the pole-skipping point (4.15), the equations of motion have two independent ingoing solutions, hence we cannot find a unique boundary Green's function by using the prescription [4,10]. However, one of the solutions becomes divergent as soon as we move to any infinitesimally nearby point in Fourier space. Consider the series expansion of the equations of motion (4.13), evaluated at where is a small dimensionless parameter, while δω and δk denote the directions in which we move away from the pole-skipping point (4.15) in momentum space. At zeroth order in the expansion, the equations (4.14) are automatically satisfied, as they are evaluated at the pole-skipping values. However, at linear order in we get a non-trivial constraint which allows us to write These relations depend explicitly on the slope (δω/δk), which can be interpreted as the additional free parameter appearing in the regular solution if the equations of motion are evaluated directly at the pole-skipping point.

JHEP07(2020)203
A generic solution to the bulk equations of motion contains a normalizable part, which is related to the poles of the Green's function, and a non-normalizable part, related to the zeros of the correlator. However, there exist solutions that are fully normalizable ψ (n) , and fully non-normalizable ψ (nn) , which we can expand near the horizon as where we have used the decomposition (3.5), and the dots represent higher order terms in the expansion, as in (4.12). The zeroth order coefficients are related by where (δω/δk) p,z denote to the directions in which we need to move away from the poleskipping point to obtain a normalizable (subscript p) or a non-normalizable (subscript z) solution.
As already mentioned, a generic ingoing solution is a combination of the normalizable and non-normalizable components ψ = ψ (n) + ψ (nn) , (6.7) which initially both contain a degree of freedom that we choose to be ψ respectively. If these are left unconstrained, the solution has too many free parameters to uniquely determine the Green's function, which happens at the pole-skipping point. Away from this location, the Dirac equation relates the coefficients through where R(δω/δk) depends on the slope of the displacement in momentum space. Due to the prescription of [10], the boundary Green's function depends on this direction as well, as The relation (6.8) is obtained by inserting (6.7) into the Dirac equations evaluated at the horizon. Near the pole-skipping point, at linear order in , we find the relation , (6.10)

JHEP07(2020)203
from which we can easily read off the proportionality factor R. Ignoring all unimportant factors, 11 one finds G R (ω n + δω, k n + δk) ∝ δω − δω δk z δk δω − δω δk p δk , (6.11) and thus the Green's function for fermion fields also takes on the pole-skipping form near the special locations. This form shows that the Green's function is infinitely multi-valued at the pole-skipping points. Namely, at any finite values of δω and δk, the slope of the displacement from the special locations (δω/δk) uniquely relates the components of the spinor through (6.10), and similarly fixes the value of the correlator by (6.11). However, directly at the pole-skipping point this slope is ill-defined and can be considered as an additional free-parameter of the solution, with the singularity structure of the Green's function given explicitly by the pole-skipping form.
From (6.10) it is clear that at (4.15) the Green's function always takes the pole-skipping form, or equivalently, is never anomalous. This can be heuristically justified by recalling that for the scalar field anomalous points usually appeared when two pole-skipping points coincided [20]. For the fermion field, there is only one such point at the frequency ω = ω 0 , so there exists no other pole-skipping point with which it can collide. Therefore this poleskipping point is never anomalous and subsequently there is always a line of zeros and a line of poles of the Green's function that intersect at (4.15).

Near higher Matsubara frequencies
To derive the form of the Green's function near the pole-skipping points at higher Matsubara frequencies, we employ the method that involves the second order differential equations due to the similarity with the scalar field case, which was analyzed in [20].
Without loss of generality, we work with the component ψ + . If we evaluate the equations (4.21) at a generic point in momentum space, we can express all coefficients appearing in (4.12) in terms of ψ (0) + only. Solving the equations D (p) + = 0 order by order, we find a generic relation where q = 1, 2, 3, . . . , the matrix M (q) + is defined in (4.30), and we have assumed that ω does not equal any of the fermionic Matsubara frequencies, i.e. N (q) does not vanish. Now evaluate this relation in the vicinity of a pole-skipping point ω = ω q + δω , k = k q + δk , (6.14) 11 The prefactors in (6.10) vanish or diverge if the slope takes the values (δω/δk) = ±r0/(2 h(r0)).
These correspond to the case where the (non-)normalizable solution has a definite Γ r chirality at the horizon (see (6.4)) and thus we are not able to take ψ as the free parameter, because they vanish. However, the pole-skipping form (6.11) generically persists even in this case.

JHEP07(2020)203
where ω q , k q is the location of a pole-skipping point with q = 1, 2, 3, . . . , and is a small parameter. At linear order in , the equation (6.12) reads (6.15) where N (ω q ) ≡ (q − 1)!(2πT ) q−1 , from which we can obtain a relation ψ (q) which again depends explicitly on the direction in which we move away from the poleskipping point. In particular, there exist slopes associated to normalizable (ψ . (6.17b) A generic solution contains both a normalizable and a non-normalizable part (6.7) which are related through the Dirac equation. As before, the retarded Green's function is proportional to the multiplicative factor relating the normalizable degree of freedom to the non-normalizable one ψ . This is similar to (6.10), only that in this case there exist additional prefactors because we are considering two positive chirality components (see [10]). At linear order in we obtain , (6.18) and by reading off the prefactor, we can see that the correlator has the pole-skipping form (2.20) near higher-order pole-skipping points as well.
A key condition in obtaining the pole-skipping form is that the relation (6.16) depends on the slope (δω/δk), which is not the case if ∂ k det M

Examples
The methods presented in the previous sections might be a bit abstract and an alert reader will realize that we have restrained from calculating any second order differential JHEP07(2020)203 equations or determinants of matrices like (4.39). While the methods are straightforward, the expressions quickly become rather long. To illuminate the procedure and show that our analysis matches the known results, we will consider some concrete examples.
First, we consider the case of a minimally coupled fermion on a non-rotating BTZ black hole background [15,16]. In this case the fermionic Green's function is known explicitly and we use it to verify our results. We then consider the special case of a fermion with halfinteger mass (or equivalently a dual operator with a half-integer conformal dimension). In that case, the correlation function takes a special form and we show that the near-horizon analysis still agrees with the exact result. Finally, we briefly present the case of a massless fermion propagating in the planar Schwarzschild black hole in anti-de Sitter spacetime in general dimension and show that the locations of the pole-skipping points pair up, as discussed in the previous section.

Pole-skipping points
For the non-spinning BTZ black hole, two functions that appear in the metric (4.1) are given by and h(r) = r 2 , (7.1) which implies that the Hawking temperature is given by T = r 0 /2π. We use the following set of gamma matrices where σ i are the Pauli matrices. In this case Γ r is diagonal and thus the two Weyl fermions are given by where ψ ± (r) are scalar functions. Because of the favorable choice of gamma matrices, the Dirac equation can be reduced to two coupled scalar differential equations that read To get the first pole-skipping point, we then follow the procedure given in section 4. The point in momentum space where the first pole-skipping occurs is

JHEP07(2020)203
which agrees with the general result (4.23) and (4.24). Since Γ r is diagonal, we note that at (7.5), the two independent solutions that are regular at the horizon can be written as ψ(r) = χ + χ − 1 + χ 1 (r − r 0 ) + χ 2 (r − r 0 ) 2 + . . . , (7.6) where χ ± are the two free parameters and χ n are two dimensional spinors whose components are fully determined in terms of χ ± . In general, χ n are not eigenstates of the chirality matrix.
We then use either of the two procedures presented in section 4 to find the locations of other pole-skipping points. Here, we will explicitly calculate the locations of the poleskipping points associated with the next lowest frequency using the first order differential equations, while merely stating the locations of the higher frequency pole-skiping points.
The equations (4.31) for this example read with the two matrices being All of the elements of the matrices are in this case scalar functions. The frequency of the next pole-skipping point is given by the value at which the determinant of the coefficient in front of (ψ − ) T vanishes. One finds that det M (11) = 8r 0 (3r 0 − 2iω) , (7.8) which vanishes at The easiest way to obtain the corresponding momenta is to set one of ψ ± to 0 and combine (7.7) with the zeroth order equation 12 (both evaluated at ω = ω 1 ) to obtain a system of three equations for three variables. For example, setting ψ (1) 12 The zeroth order equation reads (−ik + r0 − mr0 − 2iω)ψ

JHEP07(2020)203
The momentum values for the pole-skipping points are obtained by looking for the values at which the determinant of the matrix vanishes. One finds that which vanishes at If we set ψ − in the matrix (7.10), the determinant switches sign. This is because the coefficients multiplying ψ (1) ± in (7.7), evaluated at ω = ω 1 , only differ by a sign. This obviously does not change values of the momenta. The fact that we can simply set one of ψ (1) ± to 0 and calculate the pole-skipping points in the above way is a consequence of the fact that at any fermionic Matsubara frequency, only the combination (4.45) is constrained, which in our case is simply ψ (1) Finding the locations of other pole-skipping points follows the same pattern. We find that first few pole-skipping points are then located at (see figure 1) . . . . . .
As a final remark, one can notice that unlike in the case of a minimally coupled scalar in the BTZ black hole background (see eq (4.6) of [20]), the pole-skipping points do not occur in pairs of positive and negative imaginary momenta for a general mass. The exception is when m = 0, where one can see from (7.12) that one of the momenta vanishes and the others form pairs of the form k = ±ik n , just as in the scalar case. This is the same phenomenon observed in higher dimensions, where it is also associated with a decreased number of pole-skipping points and an increase of undetermined parameters from the near-horizon analysis.

Comparison with the exact Green's function
The exact retarded Green's function for the BTZ black hole was derived in [10]. For non-half-integer mass fermions, it is given by 4πT .

(7.13)
It has a pole whenever the argument of any of the gamma functions in the numerator hits a non-positive integer. Similarly, it has a zero whenever an argument of any of the gamma functions in the denominator is equal to a non-positive integer.

Green's function at half-integer conformal dimensions
When the mass m (or equivalently the scaling dimension ∆ of the dual operators) is halfinteger valued, the near-boundary expansion contains logarithmic terms and therefore the boundary retarded Green's function takes a different form. We focus on the case of m > 0 or equivalently ∆ > 1. The boundary retarded Green's function is then given by where ψ(z) is the digamma function and we have written the mass m in terms of the scaling dimension ∆. For a more explicit derivation of the Green's functions for half-integer mass fermions, see appendix C. Because ∆ is half-integer valued, the arguments of the gamma functions in the denominator and numerator of (7.17) differ pairwise by an integer. Thus, we can expand the ratio of the gamma functions into a product of finitely many terms as The ratio of the other two gamma functions can be found in a similar way to be 13 For n = 0, there are no solutions in the kn,q 2 branch of (7.16). 2 . The gray points correspond to the momentum written with a positive sign in and the hollow points correspond to the momenta with a negative sign as written in (7.12). We see that at half-integer values of the mass, some of the locations overlap (black circles with gray filling). These cases correspond to anomalous points and signal that a more thorough analysis of the boundary Green's function is needed.
The locations where the near-horizon analysis predicts pole-skipping, but the correlator does not take on the pole-skipping point form are anomalous. For the fermionic field, these occur when two pole-skipping points overlap (see figure 2) and can only occur for n ≥ 1.
Let us assume that the mass of the bulk fermionic field is a half-integer number and focus on m > 0. The analysis of anomalous points shows that for n < m + 1/2, there JHEP07(2020)203 are only non-anomalous pole-skipping points. For n ≥ m + 1/2, the non-anomalous poleskipping points are given by where for m = 1/2, there are no solutions in the second branch. This implies that the anomalous points are given by k n,q 1 = 2πiT (m + n − 2q 1 ) , q 1 ∈ {m + 1/2, m + 3/2, . . . n} . (7.24) Therefore, all in all, the near-horizon analysis predicts that the non-anomalous poleskipping points are located at ω n = −iπT (2n + 1), k n,q 1 = 2πiT (m + n − 2q 1 ),

Schwarzschild black hole in AdS d+2
Let the background metric be the Schwarzschild-AdS black hole in d + 2 dimensions. In this case, the functions determining the metric are given by

JHEP07(2020)203
with the Hawking temperature defined by (d + 1)r 0 = 4πT . For a convenient choice of gamma matrices in any dimension and the discussion of how to calculate the pole-skipping points in practice, see appendix A. Following the procedure outlined in section 5, the pole-skipping points at the lowest frequency are located at These include the locations for both subsystems that we discussed. Notice that if we set m = 0, these two points merge into a single point with momentum given by k = 0. Now, let us focus on the case of a massless fermion (m = 0). The next pole-skipping points are located at πT , An interesting observation is that the non-zero momenta at ω 1 for the massless fermion field coincide with the pole-skipping momenta associated with first bosonic Matsubara frequency ω = ω B 1 = −2πiT for the massless bosonic field in the Schwarzschild-AdS background [20]. This is not the case if the fields are massive. Furthermore, this ceases to hold when one compares higher frequencies.

Discussion
In this paper we have investigated the near-horizon behavior of a minimally coupled fermion in asymptotically anti-de-Sitter spacetimes. The thermal Green's function of the dual fermionic operator exhibits an ambiguity: at certain values of the frequency and the momentum, there exist multiple independent solutions to the Dirac equations that are ingoing at the horizon. As a consequence, the Green's function is not uniquely defined at these points. A pole and a zero of the Green's function collides which results in the pole not appearing. Hence, this phenomenon was termed 'pole-skipping' in the literature. The special frequencies where this happens are precisely the negative fermionic Matsubara frequencies where n is a non-negative integer. At each of these frequencies, there are in general 2(2n+1) associated values of the momentum, at which pole-skipping takes place. Generically, the JHEP07(2020)203 ingoing boundary condition at the horizon fixes half of the components of a spinor, whereas at pole-skipping points, the ingoing condition only fixes a quarter.
Interesting exceptional cases include that of a spinor in three-dimensional spacetime and the case of a massless spinor field in any dimension where there are only 2n + 1 poleskipping points for each n. These scenarios are analyzed in section 4 and 5, respectively.
The fermionic case is conceptually similar to the bosonic case [20]. In both cases, the near-horizon behavior of the fields determines the behavior of the boundary field theory correlators away from the origin in Fourier space. Furthermore, there is a similarity in that the higher the frequency of the pole-skipping point, the farther we probe into the spacetime, away from the horizon. This is manifested in the fact that the special momenta depend on higher and higher derivatives of metric functions evaluated at the horizon.
In addition, we see that the pole-skipping points in general have a similar structure in both cases. The frequency is determined purely by the temperature of the black hole and thus the surface gravity at the horizon. The momentum has two general contributions, one coming from the mass term and the other which is independent of the mass.
Despite all the similarities, there are also some differences between the bosonic and the fermionic cases. The first one is that in the scalar case, we have fewer pole-skipping points: at any strictly negative imaginary bosonic Matsubara frequency, i.e. ω = −2πiñT , with n = 1, 2, 3 . . . , there are 2ñ values for the momentum where the Green's function exhibits pole-skipping.
Another interesting difference is the existence of the pole-skipping point at the zeroth Matsubara frequency for fermions. This is due to the fact that the spinors are multicomponent objects, and thus there can exist two linearly-independent solutions which have the same behavior near the horizon. Higher (bosonic or fermionic) pole-skipping points depend in some way on the derivatives of metric functions. Since the zeroth-order fermionic pole-skipping point is independent of these derivatives, it is the most localized probe at the horizon. Furthermore, as we discussed in section 6, this pole-skipping point can never be anomalous and thus it is a robust feature of holographic Green's functions of fermionic operators for any value of the conformal dimension.
There are a few potential pathways in which one could generalize the above results. Pole-skipping has now been observed and analyzed for both bosonic and fermionic fields. It should be possible to extend the analysis to the case of a gravitino field. Furthermore, using a 2-dimensional CFT, it has been shown [37] that pole-skipping is also seen for frequencies which are non-integer multiples of πiT . These are neither bosonic nor fermionic Matsubara frequencies and could be associated with non-half integer spin particles: anyons. It would be interesting to see whether there is a corresponding bulk object, whose near-horizon behavior would explain pole-skipping at such frequencies.
Our hope is that one can get a better understanding of pole-skipping by considering more complicated, yet soluble models, such as the axion model [38,39]. This model contains an additional parameter which regulates the strength of the energy dissipation in the boundary theory. Ref. [24] discusses pole-skipping in this model for the energy density function and finds that the pole-skipping point does not change as the dissipation is increased and correctly predicts the dispersion relation of the collective excitations in the JHEP07(2020)203 boundary for both the weakly and strongly dissipating regime. It would be interesting to see whether such statements could be translated to the scalar or spinor field case.
Another point of interest might be the interpretation of the anomalous points. Anomalous points occur whenever two pole-skipping points overlap. From the example of the BTZ black hole we see that such points correspond to the locations in momentum space where two lines of poles overlap. It would be interesting to see if there is some additional physics that happens at such points.
The detailed analysis of the Green's function revealed that at (bosonic or fermionic) Matsubara frequencies, the retarded and advanced Green's function are equal. Another interesting aspect worth looking into is to see how this is manifested in the boundary theory.
Finally, we have added to the literature of properties of the boundary theories that are encoded in the near-horizon region. One may wonder if there are other universal properties of holographic theories that can be seen from simple near-horizon analysis of bulk fields.

JHEP07(2020)203
the boundary theory is odd-dimensional (d + 1 is odd) and the dual bulk theory is even dimensional, then the bulk spinor has twice as many components as the boundary theory.
Our choice of gamma matrices should reflect this counting. Following [10], we choose the representations of the bulk gamma matrices Γ a in terms of the boundary matrices γ a in the case of even (d + 1) as where γ d+2 is the analogue of the usual γ 5 matrix in flat space quantum field theories. 15 In the case of odd d + 1, the bulk theory has spinors with twice as many components. We can make the choice In the latter case, the Γ r matrix is explicitly diagonal, which is not necessarily the case in the former. Here, we go a step further. We construct a representation in which the Γ r matrix is diagonal in any dimension and in which the matrixk i Γ vi is also diagonal. In this way, we can show that any fermionic system can be effectively reduced not only to 2 subsystems, each involving N/2 degrees of freedom, but to N/2 subsystems each containing only 2 degrees of freedom. In that way every system can be reduced to solving equations similar to the BTZ example in section 7. We start with a 3-dimensional bulk spacetime. The spinors are two dimensional and so we can use the gamma matrices where σ i are the usual Pauli matrices given by We see that Γ r is a diagonal matrix, despite the bulk theory being odd-dimensional. Furthermore, we see that in this case, we have Γ v Γ x = Γ r , which is something we have demanded in the main text. One could also change the sign of any of the three matrices and still get a possible representation. 16 Let us now look at a 4-dimensional bulk theory, where the spinors have four components and the associated boundary theory has spinors with only two components. In order to make the Γ r matrix diagonal, we choose the following representation 15 One can take for example γ d+2 = i − d−1 2 γ 0 γ 1 . . . γ d . 16 However, the relation Γ v Γ x = Γ r would hold only if we simultaneously change the sign of two of the matrices in (A.5). If we change the sign of only one or all three, then Γ v Γ x = −Γ r .

JHEP07(2020)203
In this case, Γ r is diagonal, while all other matrices are obtained by tensor multiplying (from the left) the gamma matrices in the 3-dimensional theory (A.5) with σ 1 . In fact this is the representation used by [10].
In a bulk theory with 5-dimensions, both the boundary and the bulk spinors have four components. We get a representation for the 5-dimensional case by adding the ±Γ 5 matrix to the 4-dimensional set of matrices (A.7). In this case, we can choose to add Using the above definition, we can easily see that which means that we have constructed a representation where Γ r is again diagonal and again the analogue of the Γ 5 matrix. The generalization to higher dimensional cases is straightforward. When constructing the gamma matrices for a bulk theory in even dimensions, we pick where Γ a are the gamma matrices from the bulk theory in one dimension lower. In particular, notice that in this case the Γ v and Γ x matrices have the following forms and consequently When constructing gamma matrices for a bulk theory in odd dimensions, we pick We see that this way, in both even and odd dimensions, Γ r , Γ v and Γ x have the same form.

JHEP07(2020)203
As we saw in section 5, in higher dimensional cases we needed to split the spinor according to two projections, Γ r andk i Γ vi which are independent for d ≥ 2. In practice, we can simplify these conditions by using the symmetry under the rotations in the ddimensional subspace. This means that we can always rotate the system in such a way that the momentum points along the x-direction, or in other words, k x = k and k i = 0, for i = x.
In such a case,k i Γ vi ⇒ Γ vx , with Γ vx being given in (A.12). Thus, using symmetry and a clever choice of gamma matrices, both projection matrices become diagonal. Furthermore, the only four matrices that are of importance are the projection matrices Γ r and Γ vx (given in (A.10) and (A.12) respectively) and Γ v and Γ x (given in (A.11)) that mix up the components.
Using these matrices, the Dirac equations do not only separate into 2 subsystems, each containing half of the degrees of freedom, as was the generic case presented in section 5, but rather, the equations separate into N/2 subsystems, 17 each containing 2 degrees of freedom, similar to the BTZ case discussed in 7. The decoupled subsystems of two degrees of freedom are the first and the last component, the second and second-to-last, the third and third-to-last and so on. Effectively, one "peels" the matrices off layer by layer.
This introduces two-dimensional "effective" gamma matrices for the subsystems. One notices that the odd numbered subsystems have the same effective matrices, which differ from the even numbered subsystems, whose matrices in turn are also all the same. For the odd numbered layers (e.g. first and last, third and third-to-last) the effective gamma matrices are given by while for the even numbered layers, the two gamma matrices are given by Γ r e = −Γ vr e = σ 3 , Γ v e = −iσ 2 , Γ x e = σ 1 , (A. 16) and the subscript denotes either even or odd. Thus, in practice, solving the equations of motion always reduces to solving a system of two coupled first order ordinary differential equations for scalar functions. This does not mean that the number of pole-skipping points or the number of free parameters at a poleskipping point changes as, although we have N/2 independent subsystems of equations, N/4 of those produce the same pole-skipping points, while the other N/4 have pole-skipping points at the same frequency, but the opposite momenta.
As an example, we can look at the case of a 5-dimensional bulk theory with a 4dimensional spinor. According to the above procedure, one can use the following representation (A.17) 17 Where N is the total number of degrees of freedom in the spinor.

JHEP07(2020)203
and Γ y = σ 1 ⊗ σ 3 and Γ z = σ 2 ⊗ 1. We choose the momentum to be along the x-direction (k x = k and k y = k z = 0) so that the helicity matrix becomeŝ 18) This means that the four-component spinor can be written as ψ(r) = (ψ − . Each is governed by two coupled, first order differential equations for scalars.

B Details of the calculations
Here we present the detailed calculations that lead towards pole-skipping in asymptotically AdS d+2 spacetimes. In principle, the same equation also applies in the asymptotically AdS 3 case and we point out where the two cases differ. We repeat some of the steps from the main text, in order to make this calculation more or less self-contained.

JHEP07(2020)203
Since the metric is independent of the coordinates v and x i , one can insert the plane wave ansatz ψ(r, v, x j ) = ψ(r)e −iωv+ik i x i . The Dirac equation in Fourier space then reads (1−f (r))h (r) 8h(r) Γ r r(1+f (r)) 2 ∂ r − iω r + 1+f (r)+rf (r) 4 + d r (1+f (r))h (r) 8h(r) In general, this is a system of N first order coupled ordinary differential equations. In order to proceed, we want to decouple them in a way that makes the pole-skipping mechanism manifest.
We start by separating the spinors according to the eigenvalues of the Γ r matrix. Since (Γ r ) 2 = 1 and Tr(Γ r ) = 0, this implies that exactly half of the eigenvalues are +1 while the other half are −1. Therefore, we introduce and insert this decomposition into (B.6). This allows us to split the Dirac equations into two independent equations according to the subspaces for ψ ± , which we obtain by acting on (B.6) with the two projection operators (B.7). Notice, however, that Γ r Γ a ψ ± = ∓Γ a ψ ± for a = r, meaning that any action of a gamma matrix that is not Γ r changes the subspace in which the spinor lives. The two independent equations then read (1 + f (r)) 2 ∂ r − iω r + 1 + f (r) + rf (r) 4 + dr(1 + f (r))h (r) 8h(r) − m ψ + = 0 , r(1 + f (r)) 2 ∂ r − iω r + 1 + f (r) + rf (r) 4 + dr(1 + f (r))h (r) 8h(r) + m ψ − = 0 .

(B.8b)
We can see that these equations in general have a linear combination of derivatives of different spinor components. However, we can transform them into a form where each

JHEP07(2020)203
It is convenient to change to a new coordinate system (r, t, φ) → (ρ, T, X), in which the variables are defined as r 2 = r 2 + cosh 2 ρ − r 2 − sinh 2 ρ , T + X = (r + − r − )(t + φ) , T − X = (r + + r − )(t − φ). (C.3a) In these new coordinates the metric is written as ds 2 = − sinh 2 ρ dT 2 + cosh 2 ρ dX 2 + dρ 2 . (C.4) We then choose the diagonal frame, such that The spin connections in this frame are given by In the new coordinates, the metric depends only on the coordinate ρ (in the old coordinates the metric depends only on r), and hence we can expand the solutions in the basis of plane waves as ψ(T, X, ρ) = e −ik T T +ik X X ψ(ρ, k µ ) = e −iωt+ikφ ψ(ρ, k µ ) . (C.7) The momenta (ω, k) and (k T , k X ) are related via The Dirac equations in Fourier space are then given by We can choose a matrix representation such that Γ ρ = σ 3 , Γ T = iσ 2 , Γ X = σ 1 and write ψ = (ψ + , ψ − ) T . Note that the subscript of the components denotes the eigenvalue under the action of the Γ ρ matrix. We rescale the two degrees of freedom by introducing ψ ± ≡ cosh ρ ± sinh ρ cosh ρ sinh ρ (χ 1 ± χ 2 ) , z = tanh 2 ρ . (C.10) In the coordinate z, the asymptotic boundary is located at z = 1 and the horizon of the black hole is at z = 0. In these coordinates, after some algebra, the Dirac equations can be written as

JHEP07(2020)203
It is straightforward to transform these two equations into second order differential equations for a single variable. One finds that the solutions to these equations are the hypergeometric functions 2 F 1 (a, b, c; z). We first wish to calculate the retarded correlator, so we choose the solutions that are ingoing at the horizon. Such solutions are of the form where the constants are given by Inserting the solutions (C.12) into (C.10) and expanding them around the asymptotic boundary, one finds that the two spinor components behave as It is important to note that when 0 ≤ m < 1 2 , every term in ψ ± is normalizable. Therefore, either A or D can be chosen to be the source, and the other as the corresponding response.
Recall that in general bulk dimensions, the mass m of the fermionic field and the scaling dimension ∆ of the dual operator are related via ∆ = d + 1 2 + m. (C. 16) In the case of the BTZ black hole, d = 1, and we get the relation If m > 0, the source is taken to be A and the expectation value is D and the retarded Green's function in this case is given by their ratio (C. 18)