Higher curvature corrections to pole-skipping

Recent developments have revealed a new phenomenon, i.e. the residues of the poles of the holographic retarded two point functions of generic operators vanish at certain complex values of the frequency and momentum. This so-called pole-skipping phenomenon can be determined holographically by the near horizon dynamics of the bulk equations of the corresponding fields. In particular, the pole-skipping point in the upper half plane of complex frequency has been shown to be closed related to many-body chaos, while those in the lower half plane also places universal and nontrivial constraints on the two point functions. In this paper, we study the effect of higher curvature corrections, i.e. the stringy correction and Gauss-Bonnet correction, to the (lower half plane) pole-skipping phenomenon for generic scalar, vector, and metric perturbations. We find that at the pole-skipping points, the frequencies $\omega_n=-i2\pi nT$ are not explicitly influenced by both $R^2$ and $R^4$ corrections, while the momenta $k_n$ receive corresponding corrections.


Introduction
In recent years, progress in quantum many-body chaos has attracted much interest. In particular, developments in the gauge/gravity duality [1] have exhibited close connection between black hole physics and chaos in quantum many-body systems. Usually, the chaotic behavior is characterized by the out-of-time correlation functions (OTOCs), from which two characteristic parameters can be obtained, i.e. the quantum Lyapunov exponent λ L , and the butterfly velocity v B [2,3,4,5]. Within the framework of holography, the OTOC can be obtained from the shock wave analysis in the dual gravity theory [3]. In particular, black holes are argued to be the fastest scramblers [6,7], which saturate an upper bound on the Lyapunov exponent [8] λ L = 2πT .
Later, it was argued that besides the OTOCs, which are essentially four point functions, quantum chaos can also be manifested in the retarded two point functions. Numerical analysis in [9] first indicates that information of chaos previously obtained via holography from non-linear shock wave geometry can be extracted from hydrodynamic sound modes in linearized gravitational perturbation equation. More precisely, at certain imaginary values of frequency ω * and momentum k * of the sound pole of the retarded stress-energy two point function, the residue of the pole also vanishes, i.e., the pole is "skipped". At the pole-skipping point, the Lyapunov exponent can be read off from the imaginary frequency ω * = iλ L , while the butterfly velocity can be determined from the dispersion relation right at the point ω * = v B k * . In [10], this pole-skipping phenomenon was also explained in terms of the shift symmetry of an effective hydrodynamic description. The pole-skipping was also analytically studied in [11] as an universal behavior near the horizon where the timetime component of the Einstein equation (in ingoing Eddington-Finkelstein coordinates) vanishes such that the dual retarded two point function is not uniquely defined. The poleskipping phenomenon has also been checked to hold for SYK system [12] and 2D CFT with large central charge [13]. The stringy correction and Gauss-Bonnet high curvature correction to pole-skipping was investigated in [14], where imaginary frequency ω * receives no explicit correction 1 while the butterfly velocity v B does receives correction, which was shown to agree with the results obtained using the shock wave solution as in [4]. Poleskipping for CFT in hyperbolic space dual to AdS-Rindler geometry is also discussed in [15,16].
Recently, the near horizon analysis is generalized in [17] to equations of bulk fields dual to spin-0, spin-1 and spin-2 operators, and pole-skipping is found to exist in retarded two point functions of these operators. However, these pole-skipping points appear in the lower half plane of the complex frequency, in contrast to the aforementioned pole-skipping point of chaos located in the upper half plane at ω = +i2πT . This indicates that pole-skipping may not always be directly related to quantum chaos, but could be a consequence of a more general feature of near horizon bulk equations. Relevant discussions can also be found in [18,19,20].
In this paper, we continue to study pole-skipping along the line of [17] by involving the stringy correction and Gauss-Bonnet correction. It turns out that the dependence of the frequencies at the pole-skipping points remain the same as in the uncorrected case, while the momenta receives corrections. This pattern agrees with that found in [14] for pole-skipping in the upper plane. Moreover, the upper half plane pole-skipping point can also be recovered and is shown to agree with the results obtained in [14].
This paper is organized as follows. Section 2 reviews the key ideas relevant to poleskipping in the uncorrected background. In Section 3, we discuss pole-skipping in the presence of the stringy correction and obtain the corresponding imaginary values of ω and k for typical scalar operator, current operator, and stress-energy tensor, respectively. Similar analysis involving the Gauss-Bonnet term will be presented in Section 4. We conclude with a summary and discussion in the final section.

Review of key ideas
To elucidate the key ideas of [17] as well as [11] that are relevant to our discussion of pole-skipping, we will first consider a minimally coupled scalar field ϕ in the uncorrected background, i.e. AdS 5 black brane 2 , where f (r) = 1 − r 4 0 /r 4 with the horizon location r 0 . Note the metric has already been written in ingoing Eddington-Finkelstein coordinates. The scalar field obeys the equation of motion 3 Assuming that the perturbation depends only on x, in addition to time and the radial direction, using the Fourier transform ϕ(v, r, x) → e −iωv+ikx φ(r), the EOM becomes where a prime indicates taking derivative with respect to r. The holographical dual of the bulk scalar field is a scalar operator of dimension ∆ determined by the mass of the bulk field via ∆(∆ − 4) = m 2 . The retarded two point function of the scalar operator (in Fourier space) is given as [21,22] where A and B are coefficients in the asymptotic expansion of the scalar field near the boundary φ → Ar ∆−4 + Br −∆ . (2.5) Moreover, the field obeys the ingoing wave condition at the horizon. Then the poles of G R , i.e. A(ω, k) = 0, just defines the quasi-normal mode spectrum [21,23] for the scalar field perturbation. As argued in [17], there exists certain values (ω n , k n ), referred to as pole-skipping points, at which both A and B vanish, such that the retarded two point function is not well defined. As first indicated in [11] and further explored in [17], the pole-skipping points manifest themselves in the dual gravity theory as some special locations in ω and k for the bulk EOM. This can be understood as follows. Since we are in ingoing Eddington-Finkelstein coordinates where the metric functions are regular near the horizon r 0 , we can insert the near horizon expansion of the scalar field into EOM (2.3), and then expand the EOM near the horizon r 0 . Then a (infinite) series of perturbed EOM in the order of (r − r 0 ) can be obtained as where the coefficients C ij are functions of ω and k. For generic ω and k, one can solve for φ 1 in terms of φ 0 from the O[(r − r 0 ) 0 ] equation in (2.7), and iteratively obtain other φ i in terms of φ 0 order by order. Then the ingoing solution is uniquely determined (up to the normalization associated with φ 0 ), and the retarded two point function (2.4) is well defined.
However, when C 01 = 0, which gives ω = ω 1 ≡ −i2πT , φ 1 cannot be determined by φ 0 . Moreover, when C 00 is also vanishing, leading to certain value k = k 1 , φ 0 is also unconstrained. This gives the first pole-skipping points (ω 1 , k 1 ). Now the two free parameters φ 0 and φ 1 imply that the ingoing solution is not uniquely defined, leading to the pole-skipping phenomenon in the two point function in the dual field theory.
Similarly, other pole-skipping points with higher frequencies can be obtained. Indeed, in the O[(r − r 0 ) n−1 ] equation of (2.7), the vanishing of the coefficient C n−1 n gives ω = ω n ≡ −i2πT n, and thus implies that φ n is unconstrained. Moreover, with C n−1 n = 0 and generic values for k, the first n equations as a set of algebraic equations for the first n variables (φ 0 , . . . , φ n−1 ) imply that all of these variables should vanish, unless the momentum k takes some special values k n arising from the vanishing of the determinant of the coefficient matrix (2.8) -5 -Note, in particular, M 1 = C 00 . In sum, the two conditions, C n−1 n = 0, det M n = 0, (2.9) together determine the locations of the pole-skipping points (ω n , k n ). Note that the algebraic equation det M n = 0 in general produces n complex values for k n . In the following, we will investigate the effect of the stringy correction and Gauss-Bonnet correction to the pole-skipping phenomenon by performing the near horizon analysis as above. In particular, we will work out similar equations as (2.7) for various types of bulk fields, from which the pole-skipping points of the corresponding retarded two point functions are determined by the two conditions in (2.9).

Setup
The finite 't Hooft coupling correction in the SU (N c ) N = 4 supersymmetric Yang-Mills theory (SYM) in the large N c limit is holographically dual to the stringy α ′ correction in supergravity. More precisely, the leading finite λ correction comes at O(λ −3/2 ), corresponding to the α ′3 correction to Einstein gravity. The usual starting point for discussing the stringy correction (e.g. in [24,25]) is the 10D type IIB low-energy effective action [26,27,28,29] where κ 2 10 is essentially the 10D gravitational constant, R (10) is the 10D Ricci scalar, γ = α ′3 ζ(3)/8 ∼ λ −3/2 is the parameter for the leading order α ′ correction, W (10) is a fourth order high curvature term, which can be expressed in terms of the Weyl tensor C µναβ as Since the dilaton Φ decouples from the gravitational perturbation equation to leading order in the α ′ correction, it can be simply neglected in the following. As argued in [30], the RR 5-form F 5 and other fields are also irrelevant for our purpose. We will follow [31] (see also [14]) and only consider the dimensionally reduced 5D action with a correction term where κ 5 gives the effective 5D gravitational constant, W is just given by (3.2) with the 10D Weyl tensors replaced by the 5D ones. To our purpose in this paper, we will focus on the 5D action (3.3) and study the effect of the leading order γ correction on the pole-skipping phenomenon.
The background solution is the γ-corrected black brane [24,25]  The Hawking temperature receives the leading order correction with T 0 = r 0 /π being the uncorrected temperature. To facilitate our near horizon analysis, we change to ingoing Eddington-Finkelstein coordinates Then, the metric takes the form where Z vv = Z t , and up to O(γ 2 ) terms which are dropped.

Scalar field
Let us begin by considering pole-skipping in the case of a generic scalar operator. In the N = 4 SYM theory with leading finite 't Hooft coupling correction at O(λ −3/2 ), a scalar operator is dual to a bulk scalar field in the above background (3.9). As shown in [17], the retarded two point function exhibits pole-skipping at frequencies ω n = −i2πT n with n = 1, 2, 3, . . . , and corresponding complex momenta k n . Here we further explore this phenomenon by performing the near horizon analysis of the scalar EOM in the presence of the stringy correction. Compared to the uncorrected case, the EOM in the background (3.9) receives a γdependent source term, and equation (2.3) receives a γ-dependent source term as (3.11) where the source S 1 is given in Appendix A. Inserting the near horizon expansion (2.6), the above EOM (3.11) leads to a series of equations of the form (2.7). The first few coefficients C ij are listed in Appendix A. In particular, the coefficients in the leading O[(r − r 0 ) 0 ] equation become C 00 = −k 2 − r 0 m 2 r 0 + 3iω + γ120 k 2 + m 2 r 2 0 , (3.12) where f ′ 0 denotes f ′ (r 0 ). The two conditions (2.9), i.e. C 00 = 0 and C 01 = 0 in the present case, can be used to find the correction to the first pole-skipping point. Inserting the temperature (3.7), one can easily see that the coefficient of φ 1 vanishes at the frequency (3.14) It should be emphasized that T is the γ-corrected temperature in (3.7). At the same time, the coefficient of φ 0 vanishes at which can be written in terms of field theory quantities as In particular, a compact expression, perturbative in γ, for k 2 2 can be solved as where the first term recovers the result of [17] in the absence of the stringy correction, and the second term is the γ-correction to k 2 2 . The analytic expression for k 2 3 is too lengthy to be listed here, and higher k 2 n in general must be solved numerically. Thus, in the following, except for the case of vector perturbations with stringy corrections (3.29), only compact expressions for k 2 1 and k 2 2 will be presented. In sum, the near horizon analysis reveals the pole-skipping points at ω n = −2πnT i and the corresponding complex k n , for generic scalar operators. Compared with the uncorrected result in [17], although the temperature dependence of the frequencies remains the same, the relations between ω n and k n receive O(γ) corrections. This is similar to the result in [14] for the pole-skipping point in the upper half plane of complex frequency. There, the modification to k * leads to γ-corrected butterfly velocity v B . Note that in our case here, ω n /k n at the pole-skipping point is in general not directly related to v B .

Vector field
Consider a U (1) current operator J µ , dual to a Maxwell vector field A µ in the background (3.9), described by the EOM where Φ controls the effective coupling of the gauge field 4 . In the spirit of [32,23], the vector perturbations can be classified according to the O(2) symmetry in the plane normal to the direction of the momentum, chosen to be the x direction here. The perturbations A y and A z as O(2) vectors are in the transverse channel, whereas A v , A r and A x as O(2) scalars are in the longitudinal (or, diffusive) channel. EOMs of perturbations in different channels decouple. Since EOMs in the transverse channel are two decoupled equations similar to that of the above minimally coupled scalar field, the analysis and results in this channel are similar to the above results. So we will not discuss this channel in detail, and only focus on the longitudinal channel, where there is a hydrodynamic diffusion mode. We will also use the radial gauge A r = 0.
In the longitudinal channel, the perturbations are coupled to each other. However, A v and A x can form a gauge invariant variable, i.e. the electric field E, defined by Then the three equations for A v and A x can be combined into a single equation for E, where the coefficients A E and B E are given in Appendix B. Note that the two coefficients depend on γ, and will be expanded to O(γ) in the following calculation.
To perform the near horizon analysis, one can insert the expansion into (3.22) and expand it near the horizon. Analyzing each order in (r − r 0 ), one can obtain a set of equations of the same form as (2.7). In particular, applying the conditions (2.9), the leading order equation gives the first pole-skipping points at Again, the stringy effect only produces a γ-correction to T , but the T -dependence of ω 1 is not modified. However, k 1 receives an explicit γ-correction. To focus on the effect of the stringy correction, we take Z = 1 for simplicity. Then 26) and the momenta corresponding to ω 2 and ω 3 are determined from The expressions for k 2 2 can be solved as In this case, the three solutions for k 2 3 also take compact forms It is easy to see from (3.26) that k 2 1 is positive, or, k 1 is real, in contrast to the case of scalar operator (3.15 are positive, corresponding to real k 2 and k 3 . Since the γ-corrections are perturbations, which should not change the sign of the leading order k 2 n , these solutions for k 2 and k 3 are real in the presence of the stringy correction. In general, it is expected that k n has n values, of which at least one is real. The real values of k n are related to the diffusion mode in this channel. It is wellknown that [32,23] in the hydrodynamic limit ω ≪ T and k ≪ T , the diffusion mode has a pole in the retarded two point function, ω = −iD R k 2 with D R the R-charge diffusion constant, which receives the string correction, c.f. [33]. As argued in [9,11,17], the pole-skipping phenomenon places nontrivial constraints on the dispersion relation ω(k) at |ω| ∼ T , beyond the hydrodynamic region. In other words, the dispersion relation ω(k) of the hydrodynamic diffusion mode approaches (ω, k) = (0, 0) in the form of the diffusion pole, and passes through the pole-skipping points (ω n , k n ) for k large relative to T .
By comparing the magnitude of the numerical coefficient of the O(γ) correction relative to that of the leading O(γ 0 ) term in the expressions for k 2 1 , k 2 2 and k 2 3 in (3.26), (3.28) and (3.29), one can see that the ratio becomes larger for higher k 2 n . Indeed, for k 2 1 , the ratio is 135 in (3.26). For k 2 2 , the largest ratio in the two solutions in , the largest ratio in the three solutions in (3.29) is approximately 3132. Recall that these results are all obtained with γ treated as a perturbative parameter. So, for the O(γ) terms to be legitimate perturbations, γ should be constrained by an upper bound γ 1 ≡ 1/135 for k 1 , γ 2 ≡ 1/779 for k 2 , and γ 3 ≡ 1/3132, with tighter bounds for higher k n being expected. 5 In other words, higher k 2 n becomes more sensitive to γ-corrections. 6 Similar issue was also discussed in the study of the finite coupling corrections to quasinormal modes [34,35], where the upper bound on γ for the quasinormal modes is significantly increased by an effective resummation of a subset of higher order corrections arising solely from the first order O(γ) correction. We will not pursue a possible resummation scheme here, but leave it for future work.

Metric perturbation
In order to study pole-skipping of the retarded two point function of energy momentum tensor, we consider metric perturbations to the background (3.9), g µν + h µν . We focus on the Fourier transform h µν (v, r, x) → e −iωv+ikx h µν (r). For simplicity, we assume the radial gauge h rµ = 0. Then the perturbations can be classified by the O(2) symmetry along the yz plane into three decoupled channels: • O(2) tensor, scalar channel: h yz ; • O(2) vector, shear channel: h vα and h xα , α = y, z; In Einstein gravity, the gauge invariant variable h z y = h yz /r 2 in the scalar channel obeys the same equation as a minimally coupled massless scalar field in the same background geometry. In the presence of higher curvature corrections, the EOM of h z y is not exactly the same as that of the scalar field. However, the qualitative features of the pole-skipping results are not significantly different from that of the scalar field. Moreover, there is no hydrodynamic mode in this channel [23,36]. Therefore, in this paper, we will not present the detailed results in this channel, and only focus on the shear and sound channels where there are interesting hydrodynamic modes.

Shear channel
In the shear channel, we consider the metric perturbations with only h vy and h xy nonvanishing. To obtain the linearized equations in the presence of the stringy correction, following [37,33], it is more convenient to insert the metric ansatz into the action (3.3), which is then expanded to quadratic order in h µν to give an effective action for the perturbations, from which the linearized equations for h µν follow. 7 The two perturbations can be combined into one gauge invariant variable, also referred to as "master field", which obeys a single second order differential equation where the coefficients A, B, M 0 and M 1 are given in Appendix C.2. Its derivation is rather tedious, and a schematic strategy of derivation is given in Appendix C.1. The near horizon analysis by inserting into (3.31) and expanding in (r −r 0 ) leads to a series of equations of the same form as (2.7), with φ i replaced by Z 1i . The conditions (2.9) give again ω n = −i2πT n with k n receiving explicit γ-corrections. The first three k 2 n are determined by Compact expressions for k 2 1 and k 2 2 are As in the longitudinal channel of vector perturbations, here the real solutions for k n correspond to nontrivial constraints of pole-skipping on the momentum diffusion mode beyond the hydrodynamic range. Unlike the case of vector perturbations, however, here the γ-corrections can cause the originally positive O(γ 0 ) solutions k 2 1 = 6r 2 0 and k 2 2 = 4 √ 6r 2 0 to become negative, unless the parameter γ < 6/4422 ≈ 0.0014 for real k 1 , and γ < √ 6/(1248 + 3485 √ 6) ≈ 0.00025 for real k 2 . However, these are also the conditions for the O(γ) terms to be legitimate perturbations. Therefore, as long as γ is treated as a perturbative parameter, k n always have real solutions which recover the hydrodynamic dispersion relation at small k.

Sound channel
In the sound channel, following [37,33] again, the relevant perturbations can also be combined into one single master field The equation for Z 2 also takes the same form (3.31), with the coefficients given in Appendix C.3.
The near horizon analysis again leads to a series of equations. Again, we have ω n = −i2πT n, and k n receive γ-corrections. In particular, the first three k 2 n are determined by equations arising from det M n = 0 .

(3.38)
Compact expressions for k 2 1 and k 2 2 can be solved as The sound channel includes the metric perturbation h vv , which is dual to energy T 00 in the field theory. In contrast to the above pole-skipping points at the lower half plane of complex ω, the energy retarded two point function exhibits pole-skipping at the upper half plane ω * = +i2πT , as was originally studied in [9,10,11]. In the current setup, the upper half plane pole skipping point can also be identified by analyzing the equation for Z 2 in the sound channel, which is of the same form as (3.31), as will be discussed in section 5 and Appendix G.1.

Setup
In the above section, we studied the stringy correction which is essentially a fourth order curvature correction ∼ R 4 . In particular, the γW term arises as a top-down correction from a specific string theory (type IIB) to the supergravity action [26,27,28,29]. This form of correction is just one of a very few known corrections from specific string theories.
Without being restricted to specific known string theory corrections, one may take a pragmatic way to consider generic corrections, usually starting from quadratic curvature corrections The first two terms of couplings α 1 and α 2 can be eliminated by a field redefinition of the metric [38,39,40], leaving only the α 3 term. The higher curvature terms in general produce higher than second order derivatives in the EOM, and therefore the theory suffers from Ostrogradsky instability and other pathologies [41,42,43,44]. Thus, as the above stringy correction parameterized by γ, these corrections should only be regarded as perturbations, i.e. |α i | ≪ 1. However, for specific combinations of the coefficients, one may obtain the Gauss-Bonnet term (or, the Lovelock term [45] for general higher curvature terms), which still leads to second order EOM. Thereby, the theory is expected to circumvent the above difficulties plaguing generic higher curvature theories, and the coupling λ GB can be regarded as non-perturbative 8 . However, the range of λ GB is limited to due to causality violation and other issues in the dual boundary theory [39,47,48]. 9 Moreover, it was later argued in [51] that even for the bulk theory itself, there are bulk causality violation in generic higher curvature gravity, including the Gauss-Bonnet gravity and Lovelock gravity, unless an infinite set of higher spin fields are added. 10 Then the low energy effective theory obtained by integrating out these higher spin fields would modify the action like (4.2) with additional higher derivative terms, eventually making the EOM higher than second order, and bringing back the difficulties like Ostrogradsky instability. See [54] for more detailed discussions. Besides, there are other instability problems for the Gauss-Bonnet theory, such as the so-called eikonal instability (see [55] and the references therein).
Despite the above issues, many features of the Gauss-Bonnet theory are well-behaved for non-perturbative λ GB (at least classically). In particular, exact solutions [56,57] to the second order EOM, and the exact form of the Gibbons-Hawking boundary term [58] are known. Therefore, we still formally treat λ GB as a non-perturbative parameter in our discussion. Generically, λ GB can be regarded as a function of both λ and N c . In particular, as a perturbative parameter, it can be interpreted as λ GB ∼ 1/N c for λ ≫ N 3/2 c ≫ 1 as in the theory of [38]. More detailed discussions about the holographic dictionary relating λ GB to field theory parameters can be found in [40,50] as well as [14,54].
The background solution relevant here is the Gauss-Bonnet black brane [57] where the constant N GB is related to the Gauss-Bonnet coupling λ GB as and In general, for λ GB ≤ 1/4, N 2 GB ≥ 1/2. 11 If, as mentioned above, causality violation is taken into account, then (4.3) implies or approximately, 0.9000 ≤ N 2 GB ≤ 1.1667. The temperature of the black brane is To perform the near horizon analysis, we change to ingoing Eddington-Finkelstein coordi- where the metric that we will use takes the form (4.10)

Scalar field
Consider a scalar field with mass m determined by the EOM (2.2) in the background (4.4). The scalar field EOM becomes (4.11) Again, inserting the near horizon expansion of φ in (2.6) into (4.11) and performing the near horizon expansion lead to (2.7), where the first few coefficients C ij are listed in Appendix D for comparison. We find the similar pattern on the pole-skipping points, i.e., the frequencies ω n = −i2πT n exhibit no explicit N GB -dependence, while the momenta k n receive N GBcorrections. For example, the first pole-skipping point is The dependence of k 2 1 on r 0 is the same as that in the uncorrected case [17], whereas the dependence on N GB arises from the relation between r 0 and T in (4.8).
Momenta corresponding to ω 2 and ω 3 are given by from which k 2 2 can be solved as 14) and k 2 3 can also be easily solved, but the expressions are cumbersome and not illuminating, so will not be presented.

Vector field
In the Gauss-Bonnet background, the gauge invariant variable E defined in (3.21) in the shear channel obeys an equation of the same form as (3.20), with the coefficients given in Appendix E. The leading order near horizon analysis gives the first pole-skipping points where T is given by (4.8) with higher curvature correction. Similar to the scalar case, we find the same dependence of k 2 1 on r 0 as that in the uncorrected case [17], with the N GB -dependence entering through the relation between r 0 and T in (4.8).
The momenta corresponding to ω 2 and ω 3 are Again, with Z set to unity, in addition to k 2 1 = 2r 2 0 , a compact expression can be obtained for the two solutions for k 2 It is easy to see that k 2 1 is always positive, and the upper (+) solution for k 2 2 is positive except for N 2 GB = 1/2 where k 2 2 vanishes. Moreover, although the expressions for the k 2 3 solutions are too lengthy to be listed here, simple numerical analysis of equation (4.18) (with Z = 1) indicates that k 2 3 always has a positive solution. So there are real solutions for k 1 , k 2 and k 3 , which correspond to the nontrivial constraints imposed by pole-skipping on the hydrodynamic mode.

Metric perturbation
Unlike the theory with stringy correction (3.3), the EOM of Gauss-Bonnet gravity is still of second order. One can simply insert g µν + h µν into the EOM to obtain the linearized EOM for h µν on the black brane background (4.10).
As mentioned in Section 3.4, The linearized equations decouple according to the symmetry in the plane normal to the direction of propagation, which is taken to be the xdirection. Again, we study the perturbations in the shear and sound channels by inserting corresponding Fourier transform h µν (v, r, x) → e −iωv+ikx h µν (r) into the linearized EOM, assuming the radial gauge h rµ = 0.

Shear channel
In the shear channel, the relevant perturbations are h xy and h vy and the rest are decoupled from them. Following [23], the gauge invariant master field can be introduced which obeys a single second order differential equation where A 3 and B 3 are given in Appendix F.1. The near horizon analysis leads to pole-skipping points ω n = −i2πT n and corresponding k n . The momenta corresponding to the first three ω n are given from from which compact expressions can be obtained for k 2 1 and k 2 2 as (4.23)

.(4.24)
To ensure the existence of a real solution for k 1 , one must have 3 + 8N 2 GB − 8N 4 GB > 0, which implies N 2 GB < (2 + √ 10)/4 ≈ 1.2906. Compared with the range arising from the argument of causality violation, one can see that (4.7) ensures that k 1 has a real solution, corresponding to the hydrodynamic diffusion mode. An easy analytic analysis of (4.24) indicates that the upper (+) branch of the k 2 2 solutions can give real k 2 for N 2 GB < (2 + √ 6)/4 ≈ 1.1124. So this can be regarded as an upper bound which is tighter than the causality upper bound in (4.7), if the existence of real solutions for k 2 is required a priori. Furthermore, numerical analysis of the equation for k 3 in (4.22) suggests that k 3 can be real for N 2 GB < 1.0632, which is an even tighter upper bound. Based on these observations, one can expect that, in general, requiring the existence of real solutions for k n imposes a n-dependent upper bound on N 2 GB , and that this bound becomes tighter for larger n. Moreover, k n approaches zero for N 2 GB approaching its upper bound for n, such that at this particular pole-skipping point, |ω n | ∼ T , but |k n | ≪ T . This is in contrast to the generic pole-skipping phenomenon without higher curvature corrections, where both |ω n | ∼ T and |k n | ∼ T [9,11,17]. In addition, it would be interesting to further explore the physical implication when k n has no real solutions but N GB is still within the range in (4.7).

Sound channel
In the sound channel, the gauge invariant variable constructed using the relevant perturbations is given by The equation for Z 4 is again of the form where A 4 and B 4 are given in Appendix F.2.
The near horizon analysis in this case again gives the pole-skipping frequencies ω n = −i2πnT . The momenta corresponding to the first three ω n are 27) from which compact expressions can be obtained for k 2 1 and k 2 2 as 12 implying that no pole-skipping occurs at this point. Note that this value of N 2 GB is outside the range given in (4.7). In contrast, it is not hard to see, by inserting the above N 2 GB value into (4.27), that k 2 2 and k 2 3 always have finite solutions.
Again, the upper half plane pole-skipping location can be extracted from equation (4.26), as will be discussed in the next section and Appendix G.2.

Discussion
In this paper, we have studied the effect of the stringy correction (∼ R 4 ) and the Gauss-Bonnet correction (∼ R 2 ) to the pole-skipping phenomenon of typical scalar, vector and tensor operators dual to corresponding bulk fields. Of course, as one can easily check, all of our results recover the known results in the uncorrected case previously studied in [17]. Something new here is that, the general feature of these pole-skipping points, i.e., the locations of the frequencies are all given by ω n = −i2πT n, with the corrections only modify the expression of the temperature. On the other hand, the momenta k n receive explicit stringy or Gauss-Bonnet corrections. The similarity in this qualitative feature for these higher curvature corrections is in keeping with the results discussed in [36], where the quasinormal spectra of metric perturbations are shown to exhibit similar behavior regardless of the R 2 and R 4 corrections.
Moreover, the way these corrections affect the frequencies and momenta is similar to the pole-skipping point of chaos in the upper half complex ω plane as studied in [14]. For example, there, pole-skipping occurs at ω * = +i2πT , and k * = i √ 6πT (1 − γ23/2) for the stringy correction, and k * = i √ 6πT /N GB for the Gauss-Bonnet correction. It is in this way that the butterfly velocity v B = ω * /k * receives correction. This suggests that at the poleskipping points, the dependence of frequency on temperature exhibits certain universality, which is robust against the finite N c and finite 't Hooft coupling corrections which are holographically dual to the typical R 2 and R 4 corrections studied here and in [14]. Of course, it would be interesting to further investigate the robustness of this universality under more general higher order curvature corrections.
In fact, the upper half plane pole-skipping point can also be obtained by studying a special point of the sound channel equations (3.31) and (4.26), following the argument given in [17] for the uncorrected case. The special point here refers to a particular relation between k and ω at which the pole structures of the coefficients A and B in (3.31) or (4.26) change. The technical details of calculation are given in Appendix G.1 and G.2. Our results indicates that, instead of analyzing the vv component of near horizon Einstein equation, the pole-skipping point of chaos in the upper half plane of complex ω can also be obtained by analyzing the equation for the gauge invariant variables in the sound channel, even in the presence of typical R 2 and R 4 higher curvature corrections. This further lends support to the expectation that pole-skipping is a universal phenomenon holographically encoded by near horizon physics.
In the absence of any higher curvature correction, it was argued in [20] that two important parameters of chaos, i.e. λ L and v B , can be recovered, irrespectively of the channel of metric perturbations, as where ω * and k * in different channels are, Note that the results in the last two channels are just the first pole-skipping points (ω 1 , k 1 ). 13 In the presence of the higher curvature corrections studied here, ω * remains the same, but k * changes differently in the three channels. From the detailed results listed in Appendix H, one finds that only the results (c.f. Appendix G) of pole-skipping in the upper half plane in the sound channel agree with the shockwave analysis of the OTOC [14]. So this suggests that, when the higher curvature corrections are taken into account, λ L in the united expression in (5.1) still holds for all channels, whereas v B can only be obtained from the sound channel, not from the other two channels. We end this paper by noting an interesting open question worthy of further investigation. In the study of the gravitational quasinormal modes in the presence of higher curvature corrections, it has been found that there is a series of modes with pure imaginary frequencies which are non-perturbative in γ or λ GB , and absent in the Einstein gravity limit at γ = 0 or λ GB = 0 [36] (see also [54,55,62,63] for related discussions). In particular, in the shear channel, in addition to the gapless hydrodynamic diffusion mode, there exist other modes having pure imaginary frequencies in the lower half plane with real momenta. Moreover, the first of these non-perturbative modes can interact with the hydrodynamic mode when they are colliding at certain critical value γ c or λ c GB for a fixed momentum (or, equivalently, at certain momentum k c for fixed γ or λ GB ), after which the latter acquires real parts and the hydrodynamic description breaks down. Since the pole-skipping in the shear channel studied in this paper also occurs in the lower half plane, and it has been found that pole-skipping imposes nontrivial constraints on the hydrodynamic mode, it would be interesting to study the relation between the non-perturbative modes and the pole-skipping points. This requires high precision numerical methods. Hopefully, the results will be reported in a future publication.
After the completion of this paper, [61] appears in arXiv, which has some overlapping with the discussion here of the Gauss-Bonnet correction. It can be checked that the results agree where they overlap.

Acknowledgments
The calculations performed in this work were facilitated by two Mathematica packages: xAct, in particular, xCoba by David Yllanes and José M. Martín-García, and RGTC by Sotirios Bonanos. This work was in part supported by the NSFC grant 11647088, and a 333 grant from the North University of China.

A. Stringy correction to scalar field
The γ-dependent source term in the scalar EOM (3.11) is

B. Stringy corrections to vector field
The coefficients A E1 and B E1 in equation (3.22) are given from

C.1 Derivation of the equation for the master field of metric perturbations
In coordinates other than ingoing Eddington-Finkelstein coordinates, the equations for the master fields in the shear and sound channels have been discussed in [37,33,64,65,36,66,35]. The basic form of the equations is Or, equivalently, inserting Z = Z (0) + γZ (1) leads to Here we present the basic strategy to derive the equations in ingoing Eddington-Finkelstein coordinates for our discussion of pole-skipping. The EOM of the perturbations h µν are of the form where Ψ denotes h µν for notational simplicity, and the γ-dependent source G involves higher derivatives arising from the stringy correction γW in (3.3). Inserting Ψ = Ψ (0) + γΨ (1) , the above equation can also be written as The equation for the master field can be obtained as follows 1. As discussed in [66,35], we can use (C.5) to substitute the higher derivatives in G in terms of Ψ ′ 0 and Ψ (0) . Then (C.6) becomes 2. Insert the expression for the master field 14 where A, B, M 0 and M 1 are to be determined.

D. Gauss-Bonnet corrections to scalar field
The coefficients of the equations (2.7) with Gauss-Bonnet corrections are

E. Gauss-Bonnet corrections to vector field
The coefficients of the equation (4.15) for the gauge invariant variable are

F.1 Gauss-Bonnet corrections in the shear channel
The coefficients in the equation for Z 3 (4.21) are GB ω 2 f 5 + 2r 9 ω 2 (ω + 2iN GB r), (F.1) When the stringy correction is taken into account, the above special point is expected to be shifted by a γ-correction of the general form where k 1 is to be determined. The results of the pole-skipping points in the lower half plane in the main text are obtained under the implicit assumption that (G.1) does not hold, even though k 1 is unknown. Indeed, one cannot see k 1 directly from the expressions for A, B, M 0 and M 1 listed in Appendix C.3, as it is essentially a perturbative effect. However, further investigation reveals that k 1 can be extracted by considering the following requirement. Namely, as a perturbative effect, the stringy correction should not change the fact that r = r 0 is a regular singularity.
To extract k 1 , one can first redefine the coefficients A and B by absorbing the γdependent terms, and rewrite the equation (of the same form as (3.31)) for Z 2 in the sound channel as Z ′′ 2 +ÃZ ′ 2 +BZ 2 = 0, (G.2) whereÃ = A − γM 1 andB = B − γM 0 . Then, the above requirement implies thatÃ and B should not have poles higher than (r − r 0 ) −1 and (r − r 0 ) −2 , respectively. In particular, • Inserting (G.1) intoÃ, there are two contributions from A and γM 1 , To ensure that the pole at r = r 0 ofÃ is still a regular singularity, the (r − r 0 ) −2 terms must cancel, from which k 1 is determined as • Inserting (G.1) intoB leads to two contributions (G.5) The requirement of being a regular singularity at r = r 0 now means the (r − r 0 ) −3 terms must cancel. This is trivially satisfied by inserting the expression (G.4) for k 1 . So, solutions regular at the horizon with ρ = 0, 1, 2 are given by ω = i2πT, 0, −i2πT , respectively. Note that in this case, ω = −i2πT does not give pole-skipping, in that G R T 00 T 00 is independent of δω/δk, the parameter characterizing the way the point is approached .
The ω = 0 case corresponds to the hydrodynamic mode, which is already characterized by the pole in the two point function, and therefore should not concern us. One can show that ω = +i2πT ≡ ω * is indeed the desired pole-skipping point (dependent on δω/δk) in the upper plane. Inserting this value into (G.4) leads to k 1 = −i √ 6 23 2 r 0 , (G. 13) and the butterfly velocity v B = ω * 3/2ω * + γk 1 = 2 3 1 + 23 2 γ , (G.14) which agrees with the result obtained by analyzing the vv component of Einstein equation in [14].

G.2 Gauss-Bonnet correction
For metric perturbations in the sound channel in the GB corrected background 15 , the near horizon structure of the differential equation (4.26) depends on whether k 2 = 3ω 2 /(2N 2 GB ), as can be easily seen from the results in Appendix F.2, and in particular, the expression for α 4 in (F.6). The results obtained in the main text are implicitly under the assumption that k 2 = 3ω 2 /(2N 2 GB ). On the special point k 2 = 3ω 2 /(2N 2 GB ), however, the coefficients in (4.26) have different singular structures Then the indicial equation gives from which the argument in the previous subsection immediately leads to the conclusion that ω = +i2πT is the desired pole-skipping point in the upper half plane. Inserting this value into k 2 = 3ω 2 /(2N 2 GB ) leads to 20) and the butterfly velocity v B = 2/3N GB , which agrees with the result of [14].