Pole skipping and chaos in anisotropic plasma: a holographic study

Recently, a direct signature of chaos in many body system has been realized from the energy density retarded Green’s function using the phenomenon of ‘pole skipping’. Moreover, special locations in the complex frequency and momentum plane are found, known as the pole skipping points such that the retarded Green’s function can not be defined uniquely there. In this paper, we compute the correction/shift to the pole skipping points due to a spatial anisotropy in a holographic system by performing near horizon analysis of EOMs involving different bulk field perturbations, namely the scalar, the axion and the metric field. For vector and scalar modes of metric perturbations we construct the gauge invariant variable in order to obtain the master equation. Two separate cases for every bulk field EOMs is considered with the fluctuation propagating parallel and perpendicular to the direction of anisotropy. We compute the dispersion relation for momentum diffusion along the transverse direction in the shear channel and show that it passes through the first three successive pole skipping points. The pole skipping phenomenon in the sound channel is found to occur in the upper half plane such that the parameters Lyapunov exponent λL and the butterfly velocity vB are explicitly obtained thus establishing the connection with many body chaos.


Introduction
A large number of studies has been conducted in recent years to quantify the chaotic behavior of quantum systems with large number of degrees of freedom. Classically, the chaotic behavior of a dynamical system is characterized by a parameter λ L , known as the Lyapunov exponent. A positive value of λ L indicates an exponentially fast growth (w.r.t.

JHEP03(2021)232
time t) of separation between two phase space trajectories which were infinitesimally close at some initial time t 0 . In quantum system, the characteristics of chaotic behavior is somewhat analogous to that of the classical systems. However, the mathematical tools required to diagnose or quantitatively compute the chaotic behavior in quantum systems are different from the classical ones. For instance, in quantum systems the random matrix theory is one such commonly used tools to describe chaos [1]. However, due to the progress in the last few years it turns out that a more satisfactory diagnosis of quantum chaos can be achieved from the study of black hole, in particular using the holographic principle. 1 The gauge/gravity duality and for that matter, the AdS/CFT correspondence [3][4][5] provides a significant improvement of our understanding on quantum systems with large number of degrees of freedom at strong coupling. Using the tools of gauge/gravity duality one can perform the gravitational shock wave analysis [6][7][8][9][10] to calculate the out-of-time ordered correlation function (OTOC) which is regarded as the measure of chaos in quantum systems. The OTOC characterises the chaotic behavior in many body quantum systems in terms of two parameters, namely the lyapunov exponent λ L and the butterfly velocity v B . A series of comprehensive work has already been done over the past few years in order to establish this connection between OTOC and quantum chaos, for instance see [6,7,9,[11][12][13][14][15][16] and the references therein. More precisely, in chaotic systems the OTOC or essentially the four point correlation function shows the following exponential growth w.r.t. time and space,

1)
V and W being some generic operator. The butterfly velocity v B represents the speed at which the perturbation propagates in space. However, most recently it is observed that a particular component of a much more simple correlation function namely, the energy density retarded Green's function can provide a direct signature of quantum chaos in most of the holographic theories with Einstein gravity [17][18][19][20][21].
In strongly coupled quantum field theories at finite temperature, thermal retarded Green's function encode information about the near equilibrium physics of the system. Using the concept of gauge/gravity duality the properties of these thermal Green's functions has been investigated for strongly coupled systems in [22][23][24][25]. It is recently observed that the dispersion relation of collective excitations in the energy density Green's function is actually related to the particular form of the OTOC as given in (1.1). In particular, the parameters λ L and v B can be obtained by analysing the behavior of quasinormal modes in holographic systems. The phenomenon that sets up a direct relation between these collective excitations in the retarded Green's function and the parameters of chaos in quantum many body systems is known as the 'pole skipping' phenomenon. It is defined as the special locations in the complex (w − q) plane such that the lines of zeroes and the lines of poles of the retarded Green's function in momentum space coincide or crosse each other at that special point and hence it is not defined uniquely.
The non-uniqueness of Green's function at special location in the complex plane was explicitly shown for example in [18,19] by performing a near horizon expansion of the JHEP03(2021)232 equation of motion. To explain the basic idea with a particular example, let us consider the equation of motion for a massive scalar field φ in planner AdS black hole background with horizon radius r H [18]. Writing the metric in terms of the ingoing Eddington-Finkelstein coordinates (v, r, x) with φ = ϕ(r)e −iwv+iqx , one obtains the second order equation of motion for ϕ(r). Now for the near horizon behavior of the solution, ϕ(r) can be expanded as ϕ(r) = (r − r H ) λ with the two possible results for λ as λ 1 = 0, λ 2 = (iw/2πT ), T being the black hole temperature. At this point one imposes the ingoing wave boundary condition at the horizon and picks the exponent λ 1 which makes the solution regular at the horizon. However at special values of w n = −i2πT n, (n = 1, 2, 3, . . .), the other exponent λ 2 becomes a positive integer and hence there exist two independent ingoing regular solutions at the horizon. As a result the retarded Green's function can not be defined uniquely. In fact it turns out that the retarded Green's function actually depends on the slope with which one approaches the pole skipping points.
In order to have a clear idea about the phenomenon, it would be helpful to describe in short the procedure to obtain the pole skipping points at different orders of expansion near the horizon. As already discussed, the pole skipping points are defined as the special locations in the complex (w − q) plane at which the pole of the retarded Green's function is skipped because the numerator as well as the denominator (of the Green's function) vanishes simultaneously. These special points are represented by the set of values (w n , q n ) obtained by the near horizon analysis of the equation of motion involving the bulk fields in the dual gravity theory. In the above the index n can take positive integer values (n = 0, 1, 2, 3, . . .) that indicates the order for the near horizon expansion of the bulk equation. Let us consider an equation of motion for any bulk field Z(r) having the following general form, where r denotes the radial coordinate of the dual gravitational background and the horizon is defined at r = r H . In this analysis we will be using the ingoing Eddington-Finkelstein coordinates in which the metric for a generic gravity background takes the following form, where v is defined in terms of the tortoise coordinate r * as v = t + r * and the ellipses indicates the spatial part of the metric. To proceed further, consider the near horizon expansion of the bulk field Z(r) as, and then put it back into the equation of motion (1.2). The resultant equation can be expanded in a power series near the horizon r H and are given at different orders as,

JHEP03(2021)232
Considering equations upto nth order one can construct a n × n matrix with elements as the coefficients of the above set of equations is given as, where the above matrix elements are functions of w and q. With the above matrix, the pole skipping points are determined by the solutions of the following equations, In this work, we manage to solve the above equation analytically only for the first few pole skipping points and one has to rely on numerical solutions for the results at higher orders. Several interesting aspects of quantum chaos can be realised from the phenomenon of pole skipping. The connection between hydrodynamics and chaos has been explicitly shown in several holographic theories using the pole skipping. In [26], the authors showed through numerical calculation that for the hydrodynamic sound mode, the dispersion relation provides the results of lyapunov exponent, a parameter for quantum chaos. This connection with hydrodynamics was shown to be valid even for gravitational theories with curvature squared correction in [27] using the pole skipping phenomenon. Also in [28] the author obtained the higher curvature correction to the special value of the momentum at which the pole is supposed to be skipped, with no correction to the results for the frequency. Another interesting connection between the transport coefficient and the parameter of quantum chaos was developed in [29,30]. Moreover, the diffusion constant for both charge and momentum is shown to be related to the square of the butterfly velocity in strongly coupled theory. In holographic theories, both transport coefficients and the parameters of chaos are related to the near horizon physics. So the connection can be realized from the AdS/CFT correspondence. A list of several other works recently done on pole skipping phenomenon can be found in [31][32][33][34][35][36][37][38]. 2 In this paper we have considered a gravitational background with a spatial anisotropy as obtained in [40,41]. 3 The primary goal of this paper is to find explicitly the corrections that the pole skipping points receive in the complex frequency-momentum plane due to the spatial anisotropy in the background theory parameterized by a (or a dimensionless one b = a/T ). We consider the bulk scalar, axion and metric field perturbations and in each case we make two different choices for the direction of propagation for the field fluctuations, namely (i) along the direction of anisotropy (ii) perpendicular to the direction of anisotropy. We find that the frequency at the pole skipping point receives no correction due to the anisotropy but only the momentum gets corrected. The rest of the paper is organized as follows, In section 2, we present a short discussion on the gravitational background dual JHEP03(2021)232 to a SYM plasma with spatial anisotropy. In section 3, we describe quantitatively the pole skipping phenomenon in energy density Green's function in the upper half complex plane by doing a near horizon analysis of the vv component of Einstein's equation and obtain the Lyapunov exponent and butterfly velocity associated to the phenomenon of chaos. In section 4, we first study the pole skipping in the lower half plane for scalar and axion field and then carry on the analysis for the metric perturbations. For the metric perturbations we take into account all the non-zero components for the shear and the sound channel and construct the corresponding master equations involving the gauge invariant variables (we present a detailed discussions on the construction of gauge invariant variables for different metric field perturbations in appendix A). In the same section we also solve the dispersion relation of the transverse momentum diffusion in the shear channel using numerical method and showed that it passes through the corresponding pole skipping points. Finally, we conclude in section 5.

Details of the anisotropic background
In this section we will briefly describe the supergravity solution as obtained by the authors in [40,41] which is dual to a spatially anisotropic strongly coupled SYM theory at finite temperature. In relativistic heavy ion collision the plasma that is created has been found to be locally anisotropic for a very short time period due to the pressure difference along the longitudinal and transverse direction. This motivates the authors towards a dual gravitational background with anisotropy along a spatial direction. The five dimensional action involving the metric (g), the dilaton (φ) and the axion (χ) field excitation is given as, with 2k 2 = 16πG as the five dimensional gravitational constant. The supergravity solution is given by the following five dimensional metric as [40,41], where the spatial anisotropy is considered along the x 1 direction. The above solution is static, completely regular on the horizon and also asymptotically AdS. Notice that the axion field χ is linearly proportional to x 1 , the anisotropic direction and the dilaton field φ depends on the radial coordinate r. The black hole horizon is located at r = r H with the boundary at r = ∞. The explicit form of the different metric components in (2.2) is given in [40,41], where the authors have introduced an anisotropy parameter a. Also, assuming the anisotropy to be weak (a/T 1, T being the hawking temperature), the authors have kept terms only upto quadratic order in a in the series expanded form of the metric components of (2.2). It is important to note that the anisotropy parameter a is a dimensionfull quantity, [a] = dim [Length]. The hawking temperature is given upto JHEP03(2021)232 quadratic order in a as, Varying above action (2.1) with respect to g µν , φ and χ one gets the Einstein equations as well as the equations of motion for the scalars as, However for the near horizon analysis we need to write the above five dimensional metric in terms of the ingoing Eddington-Finkelstein coordinates defined as, The metric (2.2) can be rewritten in Eddington-Finkelstein coordinate as, where the metric components including the correction due anisotropy are given as, (2.10) Using the above metric, in the following sections we will do the near horizon analysis of the EOM for different field perturbation to obtain the special points in the complex (w − q) plane where the pole skipping phenomenon will be explicit.

Pole skipping in energy density Green's function
Recently it is explicitly shown that an universal description of the chaotic behavior in many body system can be achieved by a hydrodynamical effective field theory [43]. 4 In other words, this effective field theory predicts that the exponential growth of the OTOC can be realized from the energy density retarded Green's function in a sense that it exhibit pole

JHEP03(2021)232
skipping at a particular value of frequency and momentum that is directly related to the parameters λ L , v B appearing in the exponential form of the OTOC as, So one can obtain λ L , v B from the energy density Green's function using the pole skipping phenomenon. In this section we will try to obtain the explicit form of λ L and v B from the near horizon expansion of Einstein's equation. In particular we would be interested in the corrections that these parameters receives due to the spatial anisotropy in the background theory. In turns out that only the vv component of Einstein's equation has to be computed near the horizon in order to determine the pole skipping point. We will start by considering the small perturbation of the above unperturbed metric (2.9). There are two different choices that one can make regarding the direction of propagation for the metric perturbation, (i) Perturbation along the direction of anisotropy, (ii) perturbation perpendicular to the direction of anisotropy. In the following we will consider these two cases separately to study the phenomenon of pole skipping for the metric fluctuation in the sound channel. Let us now consider the following linear perturbations of the above fields as, µν , φ 0 , χ 0 represents the background values of the fields and h µν , ϕ, ψ are their linear perturbations. For computational simplification here we will work in radial gauge such that all the components of metric perturbation which are of the form h rµ are zero for all µ. Substituting the above linear fluctuation of the bulk fields to the corresponding equations of motion as given in (2.5), (2.6), (2.7), we obtained the following linearized equations as, µν are the linearized fluctuations to the affine connection and the ricci tensor respectively defined as, JHEP03(2021)232

Perturbation parallel to the direction of anisotropy
We first consider the perturbations to propagate along the direction of the anisotropy, x 1 so that one can use the fourier transform to write the same as, Moreover, with this particular choice of the field fluctuation one can categorize all the metric perturbations into three different modes depending on their transformation under the SO(2) rotational symmetry in the x 2 −x 3 plane [45], namely (i) Scalar mode, (ii) Vector mode, (iii) Tensor mode. In the following we will write down the nonzero components of these three modes of metric perturbation.
For the computation of special point in the complex (w − q) plane we consider only the sound modes of metric perturbation which corresponds to the retarded Green's function for the temporal component of the energy momentum tensor, G R T 00 T 00 . To proceed further we consider the near horizon expansion of the above fluctuations to have the following form, where Y represents in general the fluctuations of the metric, scalar and the axion field. The reason for the above near horizon expansion is due to the fact that the location of the special point depends on the near horizon value of the background metric. Substituting (3.8) into the linearized Einstein equation one gets the following result for the vv component in the near horizon limit as, The above equation is identically satisfied for the particular value of w and q, where, we define the dimensionless quantity b defined as b = a/T 1, in the limit a T . The Lyapunov exponent and the butterfly velocity can be obtained as (3.1),

JHEP03(2021)232
The Lyapunov exponent takes the maximum value allowed by the chaos bound even in the presence of a spatial anisotropy and only the butterfly velocity receives a correction due to the anisotropy.

Perturbation perpendicular to the direction of anisotropy
We now consider the perturbation along the x 2 coordinate, that is perpendicular to the direction of anisotropy. So the perturbation of the fields in this case can be written as, (3.12) In this case the non zero components of the metric perturbations for the scalar, vector and the tensor mode are given as, Again analysing the vv component of the linearized Einsteins equation near the horizon one gets, The corresponding value of w and q from the above equation can be obtained as, with the butterfly velocity given as, Similar to the previous case, here the Lyapunov exponent remains the same and the butterfly velocity gets corrected. The results for the butterfly velocity as obtained here in (3.11) and (3.15) matches exactly with the results obtained in [46] using gravitational shock wave analysis.

JHEP03(2021)232 4 Gauge invariant variable and pole skipping phenomenon in anisotropic plasma
In the previous section we have considered the near horizon analysis of only the (vv)component of the Einsteins equation (3.3) in the sound channel to figure out the location of the lowest order pole skipping point and from that we also obtain the Lyapunov exponent and the butterfly velocity. However even more rigorous way of doing the same is to take into account all components of the Einsteins equation and study their behavior near the event horizon. For example, in the sound channel there are different components of the metric fluctuation and hence we required to solve multiple equations simultaneously. In particular, we will construct a gauge invariant master variable for the anisotropic background so that all the components of the Einstein's equation can be put together into a single equation which is easier to deal with. In the following, we will first discuss the computation of pole skipping points at different orders for the scalar and the axion field and then we will move on to the fluctuations of the metric corresponding to both the shear and the sound channel where the gauge invariant master variable will play an important role.

Scalar field fluctuation
The equation of motion for the scalar field follows from (2.6) with the background metric components in terms of Eddington-Finkelstein coordinates as given in (2.10). As before we have considered two separate choices for the field fluctuation to propagate along the direction of anisotropy or perpendicular to that.
Parallel case. Inserting the form of the scalar field perturbation as given in (3.7), propagating along x 1 into (2.6) one gets equation of motion for the scalar field as, where the coefficients S 1 , S 2 ,S 1 andS 2 are given in appendix B. Using the near horizon power series expansion of ϕ 0 as in (1.4) one can construct the coefficient matrix C as in (1.6). The first few elements of the same is given below as,

JHEP03(2021)232
Solving (1.7), the location of the first and the second order special points can be obtained analytically and they are given as The pole skipping points appearing at higher orders can be calculated numerically for a given value of the dimensionless parameter b.
Perpendicular case. Considering the perturbation propagating along x 2 direction we get the following equation of motion for the scalar field, where the coefficientsS 1 andS 2 are given in appendix A. In this case the first few elements of the coefficient matrix C are given as (2)) (2)) (2)) (2)) (4.5) Again only the first and the second ordered pole skipping points can be solved analytically and they are given below as, In figure 1a we manage to plot the first four pole skipping points for the scalar field in the complex (w−q) plane for both the cases with perturbation propagation parallel (denoted by the blue dots) and perpendicular (denoted by the orange star) to the direction of anisotropy. We see from the plot that the results for the parallel and the perpendicular case differ by very small amount. Also note that, the value of w at the special points at different order does not receive any correction due to non zero anisotropy but only the value of q gets a finite correction. Next we compute the pole of the retarded Green's function corresponding to some scalar operator which is dual to the bulk scalar field considered above. Using the numerical  the pole has to appear at W = −1i and W = −2i respectively for W = w 2πT which is clearly evident from figure 2.

Axion field fluctuation
The equation of motion for the axion field fluctuation χ is almost similar to that of the scalar field discussed above. Also after performing the near horizon analysis, the components of the matrix C in (1.6) turns out to be very similar to those obtained for the scalar field fluctuation. In fact the leading ordered terms appearing in C are exactly the same. So we will not write them down again but only mention the final results for the pole skipping points at different order. Again the pole skipping points at first and second order can be solved exactly as given below for the two different cases with the perturbation being parallel and perpendicular to the direction of anisotropy.
Parallel case.
Perpendicular case.
The higher order pole skipping points are solve them by numerical methods. We have shown the locations of those points in figure 1b with the blur colored dots and the orange stars indicating the parallel and perpendicular cases respectively.

Metric perturbation
Let us now discuss the metric field perturbations with two different modes of perturbations, the vector and scalar modes which corresponds respectively to shear and sound channel. The non zero components are discussed in section 3 for both parallel and perpendicular case. However in order to make the calculations simpler we will work in a particular gauge where all the metric perturbations h rµ for all µ is set to zero. Einstein's equation for the two modes of perturbations can be cast into a closed form in terms of a single equation involving the gauge invariant variable. The constructions of the gauge invariant variables are discussed in appendix A.

Shear channel
In this section we consider the vector modes of metric perturbation. Similar to the scalar and axion field, we will consider the following two separate cases,

JHEP03(2021)232
Parallel case. In the particularly chosen gauge we have only two non zero components for the field perturbation in this case. They are defined in the fourier space as, with H vx 2 = h vv /g vv and H x 1 x 2 = h x 1 x 2 /g 11 . The two independent Einstein's equations can be written in the following form, where all the coefficients F x 1 x 2 , G x 1 x 2 , J x 1 x 2 . . . are functions of w, q and r H . We will not write their exact form as the expressions are too long. Using the gauge invariant variable Z v as given in (A.23), the above two equations can be clubbed into a single equation involving Z v which is given as, where the coefficients N ,Ñ , P andP in the above equation are given in appendix C. We consider a near horizon expansion for Z v as given in (1.4) and substitute it into equation (4.11) to get the first few components of the matrix C in (1.6),  (2)) (4.12)

JHEP03(2021)232
Solving (1.7), the first two pole skipping points are obtained as, Notice that here the pole skipping points corresponds to real momentum along with the imaginary solutions. These real momentum puts nontrivial constrains to the transverse momentum dispersion relation which will be discussed in the next subsection.

Transverse momentum diffusion in shear channel
In this section we wish to evaluate explicitly the location of the diffusion poles in the complex (w−q) plane which according to the phenomenon of pole skipping is constrained to pass through the pole skipping points as obtained in the previous subsection. In particular, we are interested in the dispersion relation that arises from the pole of the retarded Green's function associated to the transverse momentum density. The non zero components for the perturbed fields in this case again will be of vector type which in (t, u = r 2 H r 2 , x) coordinates are h tx 2 and h x 1 x 2 (again we consider a particular gauge such that h µu = 0, for all µ) with the perturbation propagating along the anisotropic direction The Einstein's equation in this case can be put in a closed form in terms of the gauge invariant variable, Z v = wh tx 1 + qh x 1 x 2 as, with the coefficients R, S,R,S as given in (C.5). Here as before, we define dimensionless frequency W and momentum Q as W = (w/2πT ) and Q = (q/2πT ). Now, to determine the dispersion relation we follow the numerical approach as given in [47,48]. First, the behavior near the horizon is determined by inserting the ansatz (1 − u) α into equation (4.18). Two possible solutions for α is obtained as α = ±iW/2 in which the solution with the negative sign is chosen to impose the incoming wave boundary condition at the horizon. Then the final solution can be written as the following power series, Imposing the following Dirichlet boundary condition at the boundary (u = 0) one determines the quasinormal modes,

Sound channel
The non zero components of the metric perturbations for the scalar modes are already given in equation (3.1) and (3.2) respectively for the perturbation propagating along or normal to the direction of anisotropy. Here, to make the discussion simpler we will consider a particular gauge so that all the metric fluctuations of the form h rµ for all µ will be set to zero.
Parallel case. Considering the scalar modes of metric perturbation the full set of fluctuations are given as, where in the above we define 22 . Now, using the above form of different fluctuations into the linearized equations (3.3), we obtain the following four linearly independent coupled differential equations for the metric perturbations as, where, H mn = {H vx 1 , H x 1 x 1 , H x 2 x 2 } and all the coefficients in the above equation namely, A mn , B mn , C mn . . . etc. are too lengthy to write in the paper but are functions of (r, w, q).
We can see that all the equations written above in (4.23) are coupled which makes it difficult to solve them. However constructing gauge invariant variables [48][49][50][51][52][53]  Here we are interested to find the pole skipping points in the upper half complex plane and for this we closely follow the analysis done in [18,28].
Following the results as obtained in section 3, the location of the special point in the absence of the anisotropy (a = 0), is given from the near horizon analysis as, In order to proceed with the near horizon analysis of equation (4.24), we must check its singularity structure near r = r H which changes at the special location as given in the above equation. In particular at q = 3/2w, in the near horizon limit M and L in the above equation is dominated by terms proportional to 1/(r − r H ) and 1/(r − r H ) 2 respectively. In presence of the anisotropy which is considered in a perturbative approximation (a T or b 1), equation (4.24) must abide by the above mentioned regularity condition at the special point. In other words, any term that appears in the near horizon expansion of M + a 2M and L + a 2L which is proportional to (1/(r − r H ) p ) with p ≥ 2 and p ≥ 3 respectively must be equated to zero [28].
Turning on the anisotropy, we expect the special point to get shifted from the value mention in the above equation (4.25). Let us assume that the coordinate of the shifted point in the complex (w − q) plane is given as, where we required to determine w 1 and q 1 . Inserting the above choice for w, q into (4.24) we obtain the following near horizon expansion for the coefficient of Z s in (4.24) as, Equating the above to zero, we obtain the following result for q 1 as, On the other hand, the near horizon expansion for the coefficient of Z s in (4.24) again with the shifted w and q is already dominated by term O (r − r H ) −2 . Using the obtained

JHEP03(2021)232
result for q 1 and keeping only the most dominating terms for the near horizon expansion of M + a 2M and L + a 2L one gets, Now, we consider the following power series ansatz for Z s , Inserting the above in (4.24) with the coefficients of Z s , Z s as given in (4.29) and solving the indicial equation for λ, one gets the following two solutions, Solving for w 1 such that λ 1 = 0, we get the final result for w and q from (4.26) as, Perpendicular case. Let us now consider the perturbation along the x 2 direction. The non zero components of the perturbations in this case are given as, Again with the above perturbations we will get four linearly independent equations similar to (4.23) which can be put in a closed form using the gauge invariant variable obtained in (A.18). The equation in terms of the gauge invariant variable Z s⊥ is given as, The expressions ofM andL are given in (D.5) and (D.7) respectively. The analysis towards the final results for w and q are exactly similar to what we have done in the previously corresponding to the perturbation that is parallel to the direction of anisotropy.
In particular in this case also we found the coefficient of Z s⊥ to behave near the horizon as ∼ 1/(r − r H ) 2 which must vanish in order for the perturbative analysis to be consistent.

JHEP03(2021)232
Considering the same ansatz for w and q as given in (4.26), the near horizon behavior of M + a 2M is given as, The above vanishes exactly for the value of q 1 given as, Again considering the power series ansatz for Z s⊥ similar to the one as given in (4.30) and solving the indicial equation the exponent λ can be solved as, such that, the value of w 1 remains the same as in (4.31), with the final results given as, The results for the pole skipping points in the sound channel as obtained in this section using the gauge invariant approach is exactly matches with the results in section 3 for both parallel and perpendicular case.

Conclusion
In the current manuscript we have done a detailed analysis of the very recently observed phenomenon called "Pole skipping" in a strongly coupled plasma with anisotropy along a spatial direction from the near horizon analysis of the equation of motions for different bulk field perturbations. We have also shown that this phenomenon helps us determine the parameters of chaos for the same anisotropic quantum theory. To this end we wish to list the following new aspects/results that we have obtained after doing the above analysis.
• In this paper we have explicitly computed the occurrence of the pole skipping points in the complex plane for an anisotropic plasma using the corresponding dual holographic set up. We have obtained the pole skipping points from the near horizon analysis of the equation of motion for three different bulk field fluctuations: scalar, axion and the metric field perturbations. For the metric field we considered both the shear modes and the sound modes of the field fluctuations. We find that only the momentum JHEP03(2021)232 value receives a correction due to the spatial anisotropy as parameterized by a or the dimensionless ratio b = a/T . For scalar, axion and vector modes of the metric perturbation the pole skipping happen to appear in the lower half of the complex plane. However in the sound channel it occurs in the upper half plane and is related to the parameters of chaos. So in this regard the current paper provides a complete description of the above phenomenon for the anisotropic plasma which is one of our primary motivation.
• As discussed earlier, the pole skipping phenomenon constraints the dispersion relation to pass through the special pole skipping points in the complex plane. In this work we have explicitly shown that the numerically obtained poles of the retarded green's function for the diffusion of transverse momentum exactly passes through the first three successive pole skipping points. Also the same kind of exact overlapping is shown to happen for the scalar field green's function.
• The connection between the quantum chaos and the modes of collective excitations is remarkably established by the phenomenon of pole skipping such that instead of four point functions of generic single trace operators in QFT one needs to find the points (at different orders) at which the associated energy density two point correlation function has zeroes in both the numerator and the denominator. The point at lowest order gives the butterfly velocity of quantum chaotic spread and also the Lyapunov exponent. One of the most important result of this paper is that even in the presence of a spatial anisotropy, the pole skipping phenomenon correctly produce the Lyapunov exponent and the butterfly velocity where the Lyapunov exponent saturates the chaos bound as expected and the butterfly velocity receives the anisotropic correction.
• In this manuscript, for the first time (to the best of our knowledge), we have constructed the gauge invariant variables regarding the metric perturbations in both shear and scalar channel for a gravitational background dual to the anisotropic plasma. In appendix A, we have discussed this construction in details. Using this gauge invariant variable one can write the Einstein's equation in a simple closed form which turns out to be very useful in determining the pole skipping points.
• Finally, as discussed in [17], a satisfactory understanding of the pole skipping phenomenon has been achieved in the gravity side in terms of a particular component of Einstein's equation which becomes trivial at the lowest order pole skipping point near the horizon. However from the field theory point of view, the reason for this phenomenon is yet to be understood. In particular, there are infinite number of special points in the complex plane where the two point correlation function has zero over zero form among which only the lowest order point has a connection to the parameters of quantum chaos. However the physical meaning of the other points is still unclear from the perspective of quantum field theory. So we hope this work would be a valuable contribution to this field of research in future.

JHEP03(2021)232
Before closing, we must mention that in [29] the author, using the gravitational shock wave analysis obtains a direct connection between the coefficient of momentum diffusion and the butterfly velocity for anisotropic background. For the special anisotropy there exists two different diffusion coefficient and hence two butterfly velocity, one along the direction of anisotropy and the other which is perpendicular to direction of anisotropy. In this paper we also compute two different butterfly velocities. Physically these two velocities indicated the speed with which the momentum diffuse in the x 1 (direction of anisotropy in our case) and x 2 direction. In [29] it was shown that the ratio of these two butterfly velocities are given by the following relation, v 2 With the metric as given in (2.2), the r.h.s. of the above equation can be easily calculated as, The above analysis can be repeated for a gravitational background which is deformed by the presence of uniformly distributed heavy quark in the dual field theory [57]. 6 Moreover it would be interesting to see how the Lyapunov exponent and the butterfly velocity modify due to the non zero quark density.

A Construction of the gauge invariant variable
In this section we present a detailed discussion towards the construction of the gauge invariant variable for the anisotropic background. 7 The five dimensional metric with anisotropy (2.9) can be rewritten following [20] as, where x a = (v, r) and x i = (x 1 , x 2 , x 3 ). Note that all the metric components in the above equation depends only on the redial coordinate r. Also unlike g ij , the matrixform of g ab has nonzero off-diagonal components. Now given the linear perturbation of the above background metric, the perturbation can be decomposed in the following way, 6 Also see [58] for the study of different entanglement measures on the same back reacted background. 7 We find the lecture as given in [59] very useful in computing the gauge invariant variable in our case.

JHEP03(2021)232
In the above decomposition A 1 , A 2 , A 3 , C are scalars, b ai is a vector and e ij is a symmetric traceless tensor. The vector b ai can be decomposed as, where B a is a scalar for a = (v, r) and B ai is a divergence free vector, that is D i B ai = 0 with D i denoting the covariant derivative with respect to the metric g ij . Furthermore the tensor e ij in (A.2) can also be decomposed as, In the above decomposition E is a scalar, E j is a vector and E ij is a symmetric traceless tensor quantity. Now consider the infinitesimal transformation of the coordinates as, x m → x m + ξ m , where m = (a, i). To study the gauge invariance of the metric perturbation under the above transformation of the coordinate x m we first note the following definition, where δ ξ denotes the gauge invariant transformation and £ ξ is the lie derivative. The infinitesimal transformation ξ m can again be decomposed as,

A.1 Scalar modes of metric perturbation
The gauge invariant transformation for the scalar modes of the metric perturbation is obtained using (A.5) as, δ ξ E = 2L. (A.7) Substituting the final relation into the second one of the above equation we get, Also using the above expression for T a in the first and the third relation of (A.7) gives, From the above equation we see that H ab and H L are gauge invariant.

A.1.1 Perturbation along the direction of the anisotropy
For the perturbation along the direction of the anisotropy, the perturbation is written as the plane wave form: h µν = e −iwv+iqx 1 h µν . So in this case different components of the perturbation can be evaluated from (A.2) as, (A.10) where we define h ss = h x 2 x 2 + h x 3 x 3 , that is the trace part of the perturbation in the (x 2 − x 3 ) plane. The above equation can be solved for C and E to get, Substituting the above expression for C and q in (A.9) we obtain the following two gauge invariant metric perturbation for the scalar mode as, Finally from the above equation the gauge invariant variable Z s for the scalar modes of metric perturbation is obtained as, Z s = 2q 2 (g vv H vv ) + 4wq (g 11 H vx 1 ) + 2w 2 (g 11 H x 1 x 1 ) − 2 g 11 w 2 + 3q 2 g vv g ij ∂ r g ij H x 2 x 2 .
(A.13) The corresponding gauge invariant variable for the dilaton (Z d ) and axion (Z a ) in a similar way can be obtained as, (A.14)

A.1.2 Perturbation perpendicular to the direction of the anisotropy
Let us now take the plane wave like perturbation to propagate along the other direction say x 2 such that it can be written as, h µν = e −iwv+iqx 2 h µν . The non zero components for JHEP03(2021)232 the scalar modes of metric perturbation are computed using (A.2) as, (A. 15) where in this case h ss = h x 1 x 1 + h x 3 x 3 . Again the above can be solved for C and E as, With this the gauge invariant metric perturbations are obtained as, The gauge invariant variable Z s⊥ is given as, Z s⊥ = 2q 2 (g vv H vv ) + 4wq (g 22 H vx 2 ) + 2w 2 (g 22 H x 2 x 2 ) − 2 g 22 w 2 + 3q 2 g vv g ij ∂ r g ij H x 1 x 1 .
(A.18) Again, in this case the gauge invariant variable for the dilaton and the axion field is obtained as,

A.2 Vector modes of metric perturbation
For the vector modes the gauge invariant transformations are given as, Substituting the second relation into the first one in the above equation one gets, Hence in this case H ai is the gauge invariant variable.

A.2.1 Perturbation along the direction of the anisotropy
The nonzero components of the vector modes of metric perturbation can be expressed as, So the gauge invariant variable for the vector modes with perturbation along the anisotropic direction is given as, (A.23)

A.2.2 Perturbation perpendicular to the direction of the anisotropy
Again considering the perturbation along the x 2 direction with the same plane wave form the nonzero components of perturbation can be written as, yielding the gauge invariant variable as,