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. 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]). 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,[14][15][16][17][18][19]), 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 Greens 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 [20]).
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 [21][22][23][24]. 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 [25][26][27][28][29]). 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 [30] for related phenomena in the case of Fermi surfaces).
As it was shown in [19,[31][32][33] 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. 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 [19] 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 [19] 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 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 poleskipping points. We conclude with a discussion in section 7. In appendix A we present explicit representations of gamma matrices that can be useful in practical applications. In appendix B we examine the form of the Green's function near a generic pole-skipping point and discuss the appearance of anomalous points. Some of the more detailed calculations omitted in the main text are collected in appendix C. In appendix D 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.

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 [19].
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 ds 2 = −r 2 f (r)dv 2 + 2dv dr + h(r)d x 2 . (2.5) The vacuum solutions (S matter = 0) with such properties are characterized by which are the BTZ black hole [14,15] if d = 1 and the planar AdS-Schwarzschild black hole solution if d ≥ 2.

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 [19], 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 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. 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 and consequently, G R OO (ω, k) is infinitely multivalued. As we show in appendix B, the Green's function has a line of poles and a line of zeros that pass through these special points. This is why these locations have been dubbed 'pole-skipping' points because the poles do not appear as they collide with the zeros [21][22][23][24]. There also exists an interesting phenomenon where we naively observe pole-skipping, but the points are anomalous, meaning that in the boundary correlator there are no intersecting lines of zeros and poles. We discuss these in more detail in appendix B.

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 [34,35] 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 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 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 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. 3 Note that the number d in ref. [10] is equal to d + 1 in our notation.
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 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 three-dimensional 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 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 2 h(r) ψ + = 0 , (4.7a) 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 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, . . . . 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 [19]. 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 ψ ± , for n > 0. Using these equations one can iteratively express all of the higher order coefficients as a linear combination of ψ (0) ± only. In this way one can explicitly construct two independent solutions that are regular at the horizon with the leading behavior (r −r 0 ) 0 . One of the solutions is parametrized by ψ − . 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 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 ψ (0) − . 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 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 χ (0) ± 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. See appendix B for more details.
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 6 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 (j) The leading order equation reads where M + . 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 + 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 [19].
However, if, in addition to ω = ω 1 , the momentum k is such that M (00) 7 To be precise M (00) + is proportional to the two-dimensional identity matrix, as is the term multiplying ψ (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  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 (4.28) 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) poleskipping 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) 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 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 (4.37) At generic values of k, the matrix M 1 (ω 2 , k) is invertible and thus (4.36) sets all the series coefficients appearing in the equation to zero. With that we see that at (4.23) the leading behavior of the solution at the horizon is (r − r 0 ). Furthermore ψ (1) c = 0 implies that for such a solution ψ (1) 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 ψ (1) ± to zero and treat the other variable as ψ (1) c . One finds that the roots of the equation (4.39) are independent of the choice of which variable we set to 0.
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 ψ (n) ± , which translates to demanding that det M (nn) (ω, k) = 0 . (4.41) 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 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 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 (ψ The equations for (ψ (−) − ) are equivalent 10 , but with k → −k and therefore we focus only on the pair (ψ (+) The equations (5.7) are essentially the same as (4.7). The only differences are the additional factor of d in one of the terms of the equations, and that in higher dimensions the spinors ψ (±) ± 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. 10 The detailed derivation of these equations together with the explicit 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 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 three-dimensional 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 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 pole-skipping 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 ψ counterpart) provides N/2 constraint equations, thus reducing the number of free parameters to N/2, which is enough to uniquely determine the boundary correlation function. At a pole-skipping point, say, (5.14), the equation  does not hold automatically at this pole-skipping point. The pole-skipping 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 space. Thus, the number of poleskipping points is halved to (2n + 1), yet, at each pole-skipping 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.

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 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 [14,15]. 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 half-integer 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 , (6.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 which agrees with the general result (4.23) and (4.24). Since Γ r is diagonal, we note that at (6.5), the two independent solutions that are regular at the horizon can be written as 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 pole-skipping 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 (ψ (1) − ) T vanishes. One finds that det M (11) = 8r 0 (3r 0 − 2iω) , (6.8) which vanishes at The easiest way to obtain the corresponding momenta is to set one of ψ (1) ± to 0 and combine (6.7) with the zeroth order equation 11 (both evaluated at ω = ω 1 ) to obtain a system of three equations for three variables. For example, setting ψ (1) The zeroth order equation reads (−ik + r0 − mr0 − 2iω)ψ (0) 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 ψ (1) + to 0 and include ψ (1) − in the matrix (6.10), the determinant switches sign. This is because the coefficients multiplying ψ (1) ± in (6.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 . . . . . .
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 [19]), 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 (6.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-halfinteger mass fermions, it is given by It has a pole whenever the argument of any of the gamma functions in the numerator hits a nonpositive integer. Similarly, it has a zero whenever an argument of any of the gamma functions 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 (6.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 so-called anomalous points (see appendix B for details) and signal that a more thorough analysis of the boundary Green's function is needed.
in the denominator is equal to a non-positive integer.
Assuming that the mass m is fixed and is not half-integer valued, we get two infinite families of lines of poles and two infinite families lines of zeros in the (ω, k) plane. The poles are located at 14) and the zeros can be found at where in all cases n = 0, 1, 2, . . .. Pole-skipping is observed whenever a line of poles and a line of zeros intersect and thus the Green's function skips a pole. This can be shown to happen precisely at for any n ∈ {0, 1, . . .} and with q 1 ∈ {0, . . . , n}, q 2 ∈ {1, . . . , n}, 12 which precisely matches our near-horizon analysis from (6.12).

Green's function at half-integer conformal dimensions
When the mass m (or equivalently the scaling dimension ∆ of the dual operators) is half-integer 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 D.
Because ∆ is half-integer valued, the arguments of the gamma functions in the denominator and numerator of (6.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 This means that the retarded Green's function will have a family of 2∆ − 2 lines of zeros, given by the equations where n ∈ 0, 1, . . . , ∆ − 3 2 and there is no solution for n = 0 in ω Z 1 . As all gamma functions cancel out, poles arise only when the argument of any of the two digamma functions is a non-positive integer. Thus there are two infinite families of lines of poles located at for n = 0, 1, 2, . . .. One can look for intersections between the lines of zeros and the lines of poles. These occur at the following values for the frequency and momentum ω n = −iπT (2n + 1), k n,q 1 = 2πiT (n + ∆ − 2q 1 − 1), k n,q 2 = −2πiT (n + ∆ − 2q 2 ), (6.22) where n ∈ {0, 1, . . .}, q 1 ∈ {0, . . . , min n, ∆ − 3 2 }, q 2 ∈ {1, . . . , min n, ∆ − 3 2 } and again there is no pole-skipping point at n = 0 for the momenta given by k 0,q 2 .
To see that even these special cases are predicted by the near-horizon behavior, we must mention the occurrence of so-called anomalous pole-skipping points. Namely, at points in the momentum space that are infinitesimally close to a pole-skipping point, the boundary Green's function takes on a certain form, which was dubbed the pole-skipping form and is given as where δω and δk are the directions in momentum space in which we move away from a poleskipping point and (δω/δk) p,z correspond to the slope of the lines of poles and lines of zeros going through the pole-skipping point.
The locations where the near-horizon analysis predicts pole-skipping, but the correlator does not take on the pole-skipping point form are called anomalous. For the fermionic field, these occur when two pole-skipping points overlap (see figure 2) and can only occur for n ≥ 1. The detailed analysis of the pole-skipping form for the fermionic field is given in appendix B.
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 are only nonanomalous pole-skipping points. For n ≥ m + 1/2, the non-anomalous pole-skipping 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} . (6.25) Therefore, all in all, the near-horizon analysis predicts that the non-anomalous pole-skipping 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 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 [19]. 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 ingoing boundary condition at the horizon fixes half of the components of a spinor, whereas at poleskipping 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 pole-skipping points for each n. These scenarios are analyzed in section 4 and 5, respectively. The fermionic case is conceptually similar to the bosonic case [19]. In both cases, the nearhorizon 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ñ = 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 multi-component 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 discuss in appendix B, 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. Poleskipping 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 2dimensional CFT, it has been shown [36] 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 [37,38]. This model contains an additional parameter which regulates the strength of the energy dissipation in the boundary theory. Ref. [23] 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 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.

A Gamma matrices in various dimensions
The gamma matrices used in the calculations satisfy the Clifford algebra relations where η ab = diag(−1, +1, +1, . . . , +1). In particular this means that the gamma matrix associated with the time direction 13 squares to −1 while the gamma matrices associated to the spatial directions square to 1.
Most of the calculations in the main text were done without referring to a particular representation of the gamma matrices. However, in practice it might be useful to choose a nice representation to easily extract the locations of the pole-skipping points. Here we present a choice we found particularly useful.
Recall that in AdS/CFT, a bulk spinor and a boundary spinor can have a different number of components, depending on the number of dimensions of spacetime. If the boundary theory is even-dimensional (d + 1 is even) and thus the bulk theory is odd dimensional, then the boundary and bulk spinors have an equal number of components. If 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 14 . 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 6.
We start with a 3-dimensional bulk spacetime. The spinors are two dimensional and so we can use the gamma matrices 13 In our case this is the Γ v matrix. 14 One can take for example γ d+2 = i − d−1 2 γ 0 γ 1 . . . γ d 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 15 .
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 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.4) 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.6). 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.
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 d-dimensional 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.11). 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.9) and (A.11) respectively) and Γ v and Γ x (given in (A.10)) 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 16 , each containing 2 degrees of freedom, similar to the BTZ case discussed in 6. 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 15) 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 pole-skipping 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 4-dimensional spinor. According to the above procedure, one can use the following representation 16) 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ŝ This means that the four-component spinor can be written as ψ(r) = (ψ (+) a denoting an independent degree of freedom with well-defined eigenvalues under Γ r and Γ vx . The Dirac equations are then split into two subsystems, one involving ψ − . Each is governed by two coupled, first order differential equations for scalars.

B Green's function near pole-skipping points and anomalous points
In the main text we have described how to obtain the location of the pole-skipping points and claimed that at such points a line of poles and a line of zeros of the boundary Green's function intersect. Here we explicitly show that this is the case by moving infinitesimally away from the pole-skipping point in momentum space. Furthermore, we show that the solution depends on the direction of the move and thus near the special locations in Fourier space, the correlator takes the pole-skipping form where ω n and k n are the frequency and momentum of a pole-skipping point respectively and (δω/δk) p,z correspond to the directions in which we need to move away from the pole-skipping point in order to obtain a normalizable or a non-normalizable solution at the boundary. As normalizable solutions correspond to poles in the Green's function, the associated direction is the slope of the line of poles, passing through the pole-skipping point. Non-normalizable solutions are related to zeros in the correlator and thus the associated direction is the slope of the line of zeros passing through the pole-skipping point.
Originally, all these calculations were performed for the energy-density component of the stress-energy tensor [23] and the minimally coupled scalar field [19]. Here we will show that analogous calculations can be done for the minimally coupled fermionic field as well.

B.1 Near the lowest Matsubara frequency
We have seen that the pole-skipping point at ω = ω 0 = −πiT is different from other points and we thus consider it separately. We also saw that this pole-skipping point comes from the interaction of the zeroth order coefficients in the spinor expansion (4.12) 17 . At the pole-skipping point ω = ω 0 = −πiT and k = k 0 = im h(r 0 ), the system has two independent solutions that are regular at the horizon and thus the boundary retarded Green's function is ill-defined. However, at any point infinitesimally close to the pole-skipping location, there exists only one independent ingoing solution. To see this, let us look at the leading order in the series expansion of the equations of motion (4.13) at where is a small dimensionless parameter. If = 0, the equations (4.14) are automatically satisfied. But at linear order in , we get a constraint relating ψ This allows us to express one of ψ (0) ± in terms of the other. The relation can be written as The relation between ψ (0) ± explicitly depends on the direction (δω/δk) in which we move away from the pole-skipping point. One can interpret this relation in a different way, as one can think of the aforementioned slope as the additional undetermined parameter of the regular solution at the pole-skipping point. The solution then has two free parameters -the overall normalization, and the direction in which we move away from the pole-skipping point in Fourier space.
One can use the relations (B.4) in the reverse way. Let us assume we found a particular solution to the bulk equations of motion, specified by certain boundary conditions, and let's expand it around the horizon. The equations (B.4) allow us to determine the slope at which the solution will approach the pole-skipping point in Fourier space. This is important because from the near-boundary analysis we know that the spinors separate into 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 Green's function (see (3.6)). By a choice of appropriate boundary conditions, we can therefore find bulk solutions which are either fully normalizable or non-normalizable at the boundary. Both can be expanded near the horizon as where ψ (n) denotes the normalizable and ψ (nn) the non-normalizable solution. Near the poleskipping point, the components of the zeroth order coefficient are then related by and we take ψ as the two free parameters associated with the normalizable and the non-normalizable solution. As mentioned above, (δω/δk) p,z correspond to the directions in which we need to move away from the pole-skipping point in order to obtain a normalizable or a non-normalizable solution at the boundary. The meaning of the subscripts will become apparent momentarily.
Let us assume that we are near the location of the first pole-skipping point at (B.2). At linear order in , the ingoing solution can be written as a linear combination of the normalizable and non-normalizable component where neither of the components is normalized. Following the above argument, both ψ (n) and ψ (nn) contain one free parameter, ψ , respectively. If both were left undetermined, the solution would have too many free parameters. However, the direction in which we move away from the pole-skipping point determines one free parameter in terms of the other. One can write the relation between the two as where R(δω/δk) is a matrix that depends on the slope. Using the prescription of [10], the boundary Green's function is then proportional to this matrix meaning that the Green's function depends on the slope as well.
To obtain the explicit form of the Green's function, we insert (B.7) into the Dirac equations, and expand them around the horizon. At linear order in , the two undetermined sets of parameters ψ .

(B.10)
Ignoring all the unimportant factors, one can see that the Green's function is proportional to which is precisely the pole-skipping form. In addition to this, the details of the location of the pole-skipping point do not enter the calculation at any point. Thus, this pole-skipping point will never be anomalous. Here, we follow the definition from [19], where a pole-skipping point was called anomalous if it appeared as a possible location from the near-horizon analysis, but the Green's function near such a point did not take the pole-skipping form (B.1). For a scalar field such anomalous points usually appeared when two pole-skipping points collided. We will shortly see that this is the case for the fermionic field as well. With that, one can understand that the pole-skipping point at ω = ω 0 can never be anomalous as it is the only pole skipping point with such frequency, hence there is no other pole-skipping point that it can collide with. Thus there will always be a line of zeros and a line of poles that will intersect at this pole-skipping point.
Finally, one might wonder if the prefactors in (B.10) cause some trouble. They are related to the prefactors in (B.6) and one notices that if they vanish, ψ  vanish as well. In that case, these particular spinor components cannot be taken as undetermined free parameters. One must rather use the other half of the spinor as the free parameter. In fact, the slopes (δω/δk) = ±r 0 /(2 h(r 0 )) denote the two special cases where the leading orders of the normalizable and non-normalizable solutions have a well defined eigenvalue under Γ r at leading order expansion around the horizon.

B.2 Near higher Matsubara frequencies
To analyze the form of the Green's function at higher fermionic Matsubara frequencies, we use the method involving second order differential equations, as it allows us to draw close comparisons with the scalar field case. For higher frequencies we find that pole-skipping points can be anomalous. For simplicity, we will again work only in the 3-dimensional bulk spacetime with two dimensional spinors. Without loss of generality, let us look at the variable ψ + . The variable ψ − is fully determined by ψ + through the first order differential equations.
Let us analyze the form of the Green's function at where ω q , k q is the location of a pole-skipping point with q > 0 and is a small parameter. If ω and k were generic points in momentum space, then we could use equations like (4.22) and (4.25) to iteratively express all ψ where M (q) is the matrix defined in (4.30) and where we assumed that ω = ω s with s < q. Now we want to evaluate the equation (B.13) in the vicinity of a pole-skipping point with ω = ω q . At the pole-skipping point ( = 0), the equation is automatically satisfied, however, at linear order in , the equation becomes 1 N (ω q ) (∂ k det M(ω q , k q ) δk + ∂ ω det M(ω q , k q ) δω) ψ Again, the direction in which we move away from the pole-skipping point determines the relation between the two coefficients.
In particular, there exist slopes associated to normalizable (ψ .

(B.17b)
At linear order in , the normalizable and non-normalizable solution have one free parameter which we take to be ψ (0) + for both solutions. If we move away now from the location of the pole-skipping in a general direction, then at linear order in , the solution will be a linear combination of the normalizable and nonnormalizable solution ψ = ψ (n) + ψ (nn) , (B. 18) where again neither of the components are normalized and hence naively the above solution has two free parameters. However, at linear order in , we get a relation between the free parameters of the normalizable and non-normalizable solution, which depends on the direction in which we move away from the pole skipping point and is given by .

(B.19)
Using the prescription (B.9), the retarded Green's function is proportional to the multiplicative factor relating the non-normalizable and normalizable component and thus has precisely the pole-skipping form (B.1).
However, in the case of higher Matsubara frequencies, we may have anomalous points. These occur whenever the determinant of the matrix (4.30) satisfies both det M(ω n , k n ) = 0 , and ∂ k det M(ω n , k n ) = 0 . (B.20) This occurs whenever we have a repeated root or in other words, when two pole-skipping points overlap. Notice that the above condition does not automatically include k n = 0 roots. This is the consequence of the determinant being a function of k and not k 2 , which is the case for the scalar field. An example of anomalous roots is given in section 6, where anomalous poleskipping points occur in the case of a fermion with half-integer mass (in units of the AdS radius) propagating in the BTZ background.

C 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.
The spin connections for this frame are given by with all other components, which are not related by symmetry to the ones above, being 0. Using these spin connections, one can calculate the Dirac equation to be (1 − f (r))h (r) 8h(r) Γ v + d r (1 + f (r))h (r) 8h(r) Γ r − m ψ(r, v, x j ) = 0 .

(C.5)
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 (C. 6) 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 ψ = ψ + + ψ − , Γ r ψ ± = ±ψ ± , P ± ≡ 1 2 (1 ± Γ r ) , (C. 7) and insert this decomposition into (C.6). This allows us to split the Dirac equations into two independent equations according to the subspaces for ψ ± , which we obtain by acting on (C.6) with the two projection operators (C.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 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 equation only contains the derivative of a single component We also observe that for d ≥ 2, the two matrices Γ r and Γ vi are independent and commuting 18 . The matrixk which is in agreement with the results obtained in [10]. Note that since the equations (C.12c) and (C.12d) are essentially the same except for k → −k, the same asymptotic behavior is observed for ψ (∓) ± as well.

D Exact fermionic Green's function for BTZ black hole
The metric of a spinning BTZ black hole [14,15] can be defined as follows where we define φ as an angular coordinate having a period of 2π. There are several parameters of the system that are given by where M is the mass, J is the angular momentum, T L and T R are the left and right moving temperature of the system respectively, and G is the Newton constant in 3-dimensions.
It is convenient to change to a new coordinate system (r, t, φ) → (ρ, T, X), in which the variables are defined as In these new coordinates the metric is written as ds 2 = − sinh 2 ρ dT 2 + cosh 2 ρ dX 2 + dρ 2 .

(D.4)
We then choose the diagonal frame, such that E T = − sinh ρ dT , E X = cosh ρ dX , E ρ = dρ . (D.5) 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 µ ) . (D.7) The momenta (ω, k) and (k T , k X ) are related via