Numerical calculation of the quasinormal frequencies for the Dirac field in a Lifshitz black brane

In the zero momentum limit we numerically calculate the quasinormal frequencies of the massive Dirac field propagating in a Lifshitz black brane. We focus on the non-exactly solvable cases for the fermionic perturbations, so that our results are an extension of the examples already reported for the massive Klein–Gordon and Dirac fields in the zero momentum limit. Based on our numerical results, we propose an analytical approximation of the obtained quasinormal frequencies of the Dirac field and compare their behavior with those of the Klein–Gordon field. We extend the results on the Klein–Gordon quasinormal frequencies already published. Furthermore, by imposing the Dirichlet boundary condition at the asymptotic region, we are able to find more general results for the fermionic exactly solvable case previously studied.


Introduction
In the study of perturbed spacetimes it can be found a well defined set of proper oscillations known as quasinormal modes (QNM), in order to find these QNM solutions particular boundary conditions have to be fulfilled. When the spacetime has an event horizon the perturbation has to be purely ingoing near it, due to the fact that classically nothing can escape from the black hole or black brane, whereas the other boundary condition depend on the asymptotic structure of the spacetime [1][2][3][4][5][6][7][8]. Associated to the QNM it is found a complex set of frequencies, named as quasinormal frequencies (QNF), which depend on the parameters of the spacetime and the perturbation itself. Thus by computing the QNF of a black hole or a black brane we can study its physical propa e-mail: aaresdepargar0800@alumno.ipn.mx b e-mail: alopezo@ipn.mx (corresponding author) erties [1][2][3]. Moreover in the AdS-CFT correspondence the QNF are useful to compute the relaxation times of the dual field theories [2,3,5,7,8].
The spacetimes which present Lifshitz scaling at the boundaries where the parameter z is the critical exponent and t, x represent the coordinates of the spacetime boundary, recently have been used to model non-relativistic systems near critical points in the holographic context [9][10][11][12]. Therefore, it is relevant to compute the QNF spectrum in Lifshitz spacetimes, so that we can compute the relaxation times of the associated dual field theory and in Refs. [5,6,[13][14][15][16][17][18][19][20] we can look at some examples. Even though the electromagnetic, scalar and fermionic perturbations are sometimes easier to study than the metric perturbations, the equations of motion that describe the dynamics of the perturbations does not always allow us to find analytical solutions, actually there exist several cases of perturbations for Lifshitz spacetimes where the problem has to be treated numerically and usually there is not an evident analytical expression which describes the numerical QNF [5,7,17,19,20]. Fortunately, there exist several methods which allow us to compute numerically the QNF spectrum, in particular the Asymptotic Iteration Method (AIM) has good results [21][22][23].
By imposing the Dirichlet boundary condition at the asymptotic region and in the zero momentum limit, we extend the results presented in Refs. [5,6], where the QNM exact solutions are found for Klein-Gordon (KG) and Dirac fields propagating in an asymptotically Lifshitz black brane [24,25], and in this paper we give more general results than in Ref. [6] for the exact solutions of the fermionic field. Furthermore, for the KG field we numerically calculated its QNF for the non-exactly solvable cases and we analyze again the examples considered for the first time in [5], we extend it for the Dirac field and we give analytical approximations for our numerical QNF of both fields.
The paper is organized as follows. In Sect. 2 we give a general outlook of the Lifshitz black brane which we are going to work with. In Sect. 3 we write the equations of motion for the Dirac field, study again the exactly solvable example of Ref. [6] and analyze a limit of the equations of motion where we find analytical solutions. In Sect. 4 we briefly summarize the numerical method that we use in the rest of the paper, the Improved Asymptotic Iteration Method (IAIM), so that in Sect. 5 we show the numerical results for both the fermionic and scalar fields. In Sect. 6 we discuss our main results. Finally, in Appendix A we give a proof of the classical stability for the Dirac field in the zero momentum limit.

Lifshitz black brane
It is said that a D-dimensional (D ≥ 3) background is a Lifshitz spacetime if it behaves as at the asymptotic region. Motivated by the results found in Refs. [5,6], we analyze an exact solution of the Einstein-Maxwell-Dilaton equations of motion given by [24,25] where and dx 2 D−2 is the line element of a maximally symmetric manifold representing a (D − 2)-dimensional plane [24,25], thus (3) is the line element of a symmetric spacetime. The spacetime (3), (4) represents a Lifshitz black brane, whose radius of the event horizon is r h and r ∈ (r h , ∞) is the interval of interest in this work. Thus, this spacetime behaves as (2) at infinity and it presents Lifshitz scaling at the asymptotic boundary. When z = 1 the holographic theory is relativistic, that is, it falls on the solution given in Ref. [7], whereas if z > 1 the holographic theory is non-relativistic. In what follows we have set the Lifshitz radius l equal to 1.
The tortoise coordinate r * for the black brane (3) is defined by with r * ∈ (−∞, 0) when r ∈ (r h , ∞), such that for D = z + 2 it is found [5] The Hawking temperature of the black brane (3) is given by [5,25] In what follows, in the Lifshitz black brane (3), we study the oscillations of the scalar and fermionic perturbations in the zero momentum limit, which are studied for the first time in Refs. [5,6]. One purpose of this work is to compute numerically the Dirac QNF for the cases z > D − 2 which are not studied in Ref. [6] and we compare them with the KG QNF given in Ref. [5].

Dirac perturbations in the Lifshitz spacetime
Let us consider a D-dimensional (D > 3) symmetric spacetime whose line element is [26,27] where the functions F, G and H depend only on the r coordinate, and dΣ 2 D−2 is a (D − 2)-dimensional maximally symmetric base manifold. It is worth to mention that the Lifshitz black brane (3) can be written as (8) if we identify and the submanifold dΣ 2 D−2 is identified with the (D − 2)dimensional plane dx 2 D−2 . In D-dimensional symmetric backgrounds, the Dirac equation reduces to a pair of coupled partial differential equations in two variables [27,28]. In order to do so, we note that the D-dimensional symmetric spacetime (8) can be written as and by using the transformation properties of the Dirac operator under conformal transformations, that is, ifg μν = Ω 2 g μν then / ∇ψ = Ω (D+1)/2/ ∇ψ and ψ = Ω (D−1)/2ψ , we find that in the spacetime (11), conformal to the symmetric spacetime (8), the Dirac operator simplifies to [27,28] where I 2 (D−2)/2 is the 2 (D−2)/2 × 2 (D−2)/2 -dimensional identity matrix, the symbol ⊗ stands for the direct product and σ i are the Pauli matrices. As in Ref. [27], we have considered the following representation of the gamma matrices The matricesγ i are a representation of the gamma matrices for a (D − 2)-dimensional space and the operators / ∇ 2D and / ∇ dΣ stand for the Dirac operators on the two-dimensional spacetime ds 2 2D and on the (D−2)-dimensional submanifold dΣ 2 D−2 , respectively [27].
Taking the Dirac spinorψ(r, t, x i ) for the spacetime (11) as with ψ 2D (r, t) a two-spinor on a two-dimensional spacetime and with in such a way that κ and χ(x i ) are the eigenvalues and the eigenfunctions of the Dirac operator on the submanifold with line element dΣ 2 D−2 [27]. From (12) and the previous facts we obtain that the Dirac equation (10) can be reduced to the pair of coupled partial differential equations where the functions ψ 1 and ψ 2 are components of a 2Ddimensional spinor and for the Dirac field we assume that the mass m is always positive (m > 0). We point out that the system of equations (16) is valid for D = 3 [29], and the black brane solution (3) is also valid in D = 3 [4,23,24]. Therefore we include the spacetime dimension D = 3 in our study. Based on the test field approximation, in the zero momentum limit we are going to compute the fermionic QNF for the Lifshitz spacetime (3), therefore we use (16) with κ = 0 which is an allowed eigenvalue because the base manifold dx 2 D−2 is a (D − 2)-dimensional plane [30]. In what follows we set an harmonic time dependence for the components of the field, which means Thus by inserting the previous expression into the system (16) with κ = 0 and by defining y = r h /r , such that y ∈ (0, 1), we obtain withω = ω/r z h . We focus on the solution for R 2 (y), by mentioning that similar results are found for R 1 (y). Hence, if we take the first derivative of the first equation in (18) with respect to y and by inserting the second equation of (18) into it, we obtain a second order differential equation given by To find QNM solutions we are going to define them based on Refs. [5,6], that is, the solution R 2 (y) of Eq. (19) is a QNM, if satisfies the boundary conditions: • Near the horizon it must behave as an ingoing wave.
• It goes to zero as r → ∞ (Dirichlet boundary condition).
To fulfill the requirement near the horizon, that is in the neighborhood y = 1, we have to look at the solution of (19) near it, thus we analyze with solutions where A 1 , B 1 are constants. 1 From (5) we find that such that near the horizon we obtain and the solutions at this limit are given by Recalling that in Ref. [6] only for m ≥ z we found exactly a well defined QNF spectrum, since for 0 < m < z the Dirichlet boundary condition was not sufficient to determine a discrete set of QNF. We see that decoupling the system (16) with κ = 0 in a different way from the one presented in Ref. [6], we are able to set the Dirichlet boundary condition for all the positive values of the mass m by choosing the y m term. It is straightforward to find the exact solutions presented for the first time in Ref. [6] by using (19). The exactly solvable example corresponds to z = D − 2 or α = 2z, hence by inserting the following ansatz into Eq. (19), we find out that f (v) is solution of the hypergeometric differential equation [31] If we suppose that c is not an integer, the solution of (27) is [31] where F(a, b; c; v) is the hypergeometric function and because we find that near the horizon From (23) to have an ingoing wave at the horizon we must set B 3 = 0, thus (26) is written as For c − a − b different from an integer, owing to the Kummer properties for the hypergeometric solution, as r → ∞ (v → 1), the solution R 2 at this limit behaves as [31] therefore we require that the second term of (33) goes to zero, because we need that R 2 (v) → 0 at the asymptotic region. This can be done if a = −n 1 or b = −n 2 with n i = 0, 1, 2, . . ., and if we look at the expressions (28) we can conclude that in order to have QNM solutions for all the positive values of the mass it is necessary that The last expression is the QNF set found in Ref. [6] with the difference that all the positive values of the mass are allowed. We have to mention that the condition a = −n 1 gives a condition over the mass which is not physical, and even though we have supposed that neither of the parameters c nor c−a −b can be integers, when we consider that they are integers, the same results are obtained, we just have to use the corresponding solutions of the hypergeometric differential equation (27) [31].
For the previous reasons we believe that in order to study the cases z = D − 2 we first use the approach of this paper instead of the one used in Ref. [6]. Thus, as we already know the appropriate behavior of R 2 (y) at the boundaries we insert into (19) the following ansatz and it is obtained that R 2 (y) must be a solution of Although we have not yet found the analytical solutions of (36) for z = D − 2, if we suppose z >> D − 2 or equivalently α ∼ z, Eq. (36) transforms into with x = y z . Once again, the differential equation (37) is of hypergeometric type (27) with parameters thus if we proceed in a similar way as for α = 2z, to fulfill the QNM requirements at the boundaries it is necessary that When we compare (34) with (39) we see that both sets of QNF have a similar analytical form. In Sect. 5 we show that for z > D − 2 the behavior of our numerically calculated Dirac QNF is well approximated by a similar expression. Since, as we mentioned before, we have not been able to treat analytically the problem of z = D − 2, we implement the numerical method known as the Improved Asymptotic Iteration Method [21][22][23] as in Ref. [5] for the scalar perturbations.

The improved asymptotic iteration method
The Asymptotic Iteration Method is developed in Ref. [21], improved in [22,23] and has been used to compute the QNF for several gravitational systems [5,17,20]. The method is based on taking derivatives iteratively of the second order differential equation [21,22] Y By taking the n-th derivative of (40) we can always keep the form of the original differential equation if we define the functions such that for n = 1, 2, . . ., where the notation Y (n+2) (x) refers to (n + 2)-derivative of Y (x). The main point of the method is that we can build the general solution of (40) if we suppose that for a certain n the quantization condition is fulfilled [21]. To obtain the QNF by using the method described above, we have to rewrite (43) as and search its roots once we evaluate it in an appropriate value of x. Also, in Ref. [21] it is shown that the condition (43) is sufficient to find the solution of Eq. (40). Previous work shows that the main problem of the AIM [21,22], is that the computation of the derivatives of λ i (x) and s i (x) for each iteration could be a hard work issue, that is why in Ref. [23] an optimization is proposed and in that paper the Improved Asymptotic Iteration Method is developed. In this method, the functions λ n (x) and s n (x) are expanded in a Taylor series, around an appropriate point x 0 , that is, and the quantization condition (44) is translated into As we mentioned previously, to obtain the QNF with the method described in the last paragraphs we use Eq. (47), by taking into account that the boundary conditions have to be set before in order to obtain an equation in the form of (40). This fact means that the functions, λ 0 (x) and s 0 (x) have already the information of the boundary conditions and therefore the condition (47) is a polynomial in ω, whose zeros are in principle the QNF. We point out that the method is implemented by solving the condition (47) for a previously selected value of n. Because we are solving the condition in an approximate way, unphysical QNF could appear. To choose the appropriate frequencies, we use physical arguments (see Appendix A, for example) and we have to pick the ones that repeat as n increases.
By recalling that (36) has already the information of the boundary conditions for the QNM, we can implement IAIM for the cases z = D − 2, if we identify In the zero momentum limit, notice that in Appendix A we show that for the Dirac field moving in the black brane (3) the imaginary parts of their QNF satisfy Im(ω) < 0 and that is why we can actually implement the numerical method since the system is classically stable.

Numerical results
In Ref. [5] the QNF of the scalar perturbations moving in the Lifshitz black brane (3) are studied. In that paper exact QNF are calculated for D = z + 2 and by implementing the IAIM, the KG QNF were obtained numerically for the cases D = z + 2. In this paper we study again the KG QNF Also we compare them with the Dirac QNF that are computed by implementing the same numerical method. For the Dirac field we also give an analytical approximation for the QNF that are obtained numerically. We begin by calculating the QNF of the Dirac field.

Dirac QNF in a Lifshitz black brane
In Ref. [6], for D = z + 2, exact Dirac QNF are calculated by imposing the Dirichlet boundary condition at the asymptotic region and in Sect. 3 we extend the previous results to all the positive values of the mass. Based on it, we implement the IAIM to compute the QNF for the non-exactly solvable cases which correspond to D = z + 2. In order to use the IAIM, we take as a basis Eq. (36) for R 2 which has already the information of the boundary conditions and by emphasizing that similar numerical results are obtained for the R 1 function. We notice that for z ≥ D − 2 the numerical results (and the analytical results (34) and (39)) point out that the QNF of the Dirac field are purely imaginary. Considering the numerically calculated QNF and due to the form of the analytical expressions (34) and (39), we found out that our numerical results for the QNF of the Dirac field for z ≥ D − 2 are well approximated by the following formula At present time we do not have an analytical deduction of this expression. However, we notice that in the appropriate limits, it reproduces the exactly solvable cases as well. For example, once we replace α = 2z into (49), we obtain the formula (34) and if we consider the limit α ∼ z it yields the expression (39). Thus, (49) is a good approximation to our numerical Dirac QNF, but we have to point out that it is only based on the numerical output that we obtained and we do not know if it reproduces the whole spectrum of QNF for the Dirac field.
For several values of the spacetime dimension for the black brane (3), in Figs. 1, 2, 3, 4 and 5 we illustrate the behavior of the Dirac QNF when we change the critical exponent z while we keep fixed m. For the cases z ≥ D − 2, we observe that |I m(ω)| increases as z increases. In more detail, we see that when z starts to grow the behavior is approximately linear Nevertheless the behavior is not linear for values of z that are close to D − 2, while for the cases z < D − 2, as z moves away from the line z = D − 2, the QNF are not longer described by (49) since, in general, the QNF are not purely imaginary as it can be verified in Tables 1 and 2 where explicit numerical data are shown. 2 The same situation happens for 2 In all the tables the headline AN refers to the values of the QNF computed from the analytical approximations or exact expressions. The headline IAIM refers to the numerical QNF obtained by implementing the IAIM. We notice that in Tables 1 and 2 the real parts of the analytical values are not shown because the proposed expression (49) produces purely imaginary frequencies. Also, we define the quantityω 0 as the frequency of the fundamental mode normalized by r z h . This quantity is useful in what follows. the KG QNF, as it was described in Ref. [5]. Notice that in Sect. 5.2 we study the QNF of the KG field again [5].
In Figs. 6, 7 and 8 and Table 3, we show particular examples with D and z fixed, but satisfying z ≥ D − 2, to illustrate how the QNF change when we vary the mass of the field. We note that for the first three modes their dependence on the mass is similar in all the cases shown. Furthermore, in Figs. 6, 7 and 8 we see that the mass of the field changes slightly the values of the Dirac QNF, at least for n = 0, 1, 2.
It is relevant to mention that in this Section we have presented representative examples, but we have calculated numerical data for other values of the parameters, and as far as we know, the results can be extended for more values of the parameters by noticing that problems of convergence of the method appear when they take higher values. We conclude that the formula (49) is a good approximation for the behavior of the QNF obtained numerically in the cases z ≥ D − 2, however for z < D − 2 the expression (49) does not work and we are far from proposing an analytical expression that approximates the values of these QNF.

KG QNF in a Lifshitz black brane
In this section we expound some numerical data for the QNF of the KG field to extend some of the previous results in Ref. [5]. First we present an analytical formula for its QNF in the limit α ∼ z that is not previously known, and in a similar way to the Dirac field, we propose an analytical approximation for the obtained numerical QNF in the non-exactly solvable cases (see the expression (63)).
The equation of motion for a massive scalar field is with m denoting the mass of the field, = ∇ μ ∇ μ is the d'Alembert operator and ∇ μ the covariant derivative. Taking into account the results of Ref. [5] and by considering that we are in the zero momentum limit, if we insert the following ansatz into (50) with the help of the tortoise coordinate (5), for the radial function of the scalar field moving in (3) we obtain a Schrödinger-like equation with (53)   We proceed in a different way from Ref. [5] to achieve higher values of the physical parameters, by recalling that the objective of this section is to compute the QNF for the scalar field in the Lifshitz black brane (3) to extend the previous published results. Thus we are going to consider the boundary conditions imposed in Ref. [5] which are the same that we have used while studying the fermionic perturbations. Equation (52) transforms into where we have used once again the variable y = r h /r . By analyzing the solutions of Eq. (54) at the boundaries, we propose the following ansatz The parameter Δ is given by and is the scaling dimension of the operator O Δ dual to φ, as it is explained in [5] and through the Breitenlohner-Freedman bound (BF bound), there exist a constrain on the scaling dimension given by [32,33] 0 < α 2 < Δ.
As for the fermionic perturbations, we can consider the limit case α ∼ z for the scalar field. In this limit, with the help of the variable x = y z , as previously, Eq. (56) transforms into Thus, following a similar procedure to that described in Ref. [5], to satisfy the boundary requirements for the QNM, in the limit α ∼ z the KG QNF must be equal to Therefore we expect that as α ∼ z the QNF of the KG behave as (60). This analytical result is not given in Ref. [5]. The exactly solvable case (z = D − 2) was already analyzed in Ref. [5] and it was found that the QNF are Owing to the previous results, it is interesting to inspect higher values for the parameters z and D, since we want to compare them with the Dirac QNF of the previous section. Furthermore, these values extend the results of Ref. [5] where they only reached z = 6, so that in this section we will go further. As it can be verified with (22) and (23) and with the first term of (55) where it is assured that the field goes to zero at the asymptotic region, Eq. (56) has already the information of the boundary conditions. Thus, taking as a basis Eq. (56), we are going to compute the KG QNF with the help of the IAIM as in Ref. [5].
From the expression (56) for the Klein-Gordon field we identify By implementing the IAIM, taking as a basis expressions (62), we reproduce Table 2 of Ref. [5] and in Table 7 we show the numerical results where we found a minimal difference. Furthermore it is found that the obtained numerical QNF for the scalar field moving in the black brane (3) for z > D − 2 tend to behave as Owing to the BF bound (58), we choose representative examples and in Figs. 9, 10 and 11 it is shown the behavior of the QNF to compare them with Eqs. (60), (61), or (63) as appropriate. To compare in more detail the numerical results  with the produced values by the proposed expression (63) in Tables 4, 5 and 6 we show explicitly some of the results for the QNF that are used to plot the Figs. 9, 10 and 11. In Figs. 9, 10 and 11 we can see that for the KG field the |Imω 0 | goes to a constant value as z increases, this is consistent with the proposed expression (63) since by taking z >> D − 2 (α ∼ z), it is found that |Imω 0 | ∼ Δ/2. For the Dirac field, as z increases, in Figs. 1, 2, 3, 4 and 5 we visualize that the |Imω 0 | increases instead of decreasing to a constant value which is also congruent with the formula (49) proposed for the numerical QNF of the Dirac field. For the overtones of both fields, as z increases, our numerical results show that |Imω| increases. Based on the proposed expressions (49) y (63) we think that the linear dependence on the mode number n contributes to this behavior. Surprisingly, for both fields the line z = D − 2 represents a breaking point, since for z < D − 2 the QNF are not purely imaginary [5]. Even though we are far from understanding the reason of it, we  believe that we have to study in detail the case z < D − 2 in a future work. We observe in Figs. 9, 10 and 11 that for the numerical data of the KG QNF, the formula (63) does not give an approximation as precise as the proposed formula (49) for the numerical QNF of the Dirac field.

Discussion
In the zero momentum limit we study the QNF of the massive Dirac and KG perturbations in the Lifshitz black brane (3). We extend the results presented in Ref. [6], since in Sect. 3 for z = D − 2, by imposing the Dirichlet boundary condition at the asymptotic region and by considering that we are on the test field approximation, we are able to find the exact QNF of the Dirac field for all the positive values of the mass, in contrast to Ref. [6] in which a restriction on the value of the mass was imposed. Based on it, we numerically calculate the fermionic QNF in the Lifshitz black brane (3) for the non-exactly solvable cases z = D − 2. Owing to the exactly solvable example and the limit case α ∼ z, where analytical QNF were computed, in Sect. 5 we were able to give the approximate analytical expression (49) for our numerical QNF with z > D − 2. For the cases z < D − 2 we found out that, in general, the QNF are not purely imaginary, and the proposed formula (49) does not work. It will be interesting to analyze in detail the last case in a future work so that an explanation can be found.
Additionally to the results for the Dirac field, we study the QNF of the scalar perturbations in the spacetime (3) which are calculated for the first time in Ref. [5] (but see Table 7). Once again based on the case α ∼ z and the exactly solvable example of Ref. [5], for which analytical expressions for the KG QNF exist, we conclude that, for fixed D, the numerical QNF of the KG field, starting from a particular value of z, tend to behave as predicts the proposed expression (63) when z > D − 2. Even though the analytical expressions (49) and (63) approximate very well the numerical results, notice that we do not have an analytical proof of these formulas.
Taking into account the Hawking temperature (7) and by considering examples where (49) and (63) approximate the numerical results, it is obtained that the numerical results of the QNF for the Dirac and KG fields respectively can be approximated by the expressions Curiously, for both fields and for fixed D, z, m, the previous formulas predict that the spacing between consecutive QNF is given by From the numerical results and the proposed approximations (49) and (63), we can state that the fermionic QNF depend on the physical parameters in a different way that the QNF of the scalar field, as it was already commented in Ref. [6]. In Sect. 5.1 for z ≥ D − 2 we also find that when z varies the Dirac QNF behave in a similar way for the different dimensions of the spacetime that we have studied. Moreover we see that for the analyzed dimensions of the spacetime the QNF of the Dirac field change with the mass in a similar way. As far as we know these facts for the Dirac QNF in the black brane (3) are not previously described.
We can compute that the decay time from [1][2][3] where ω 0 is the fundamental QNF, as previously. By assuming that (64) and (65) are valid (at least for the values of the physical parameters studied in this work) for the Dirac and Klein-Gordon fields we obtain that their decay times are well approximated by In the limit z → ∞, for both fields the previous expressions predict that the decay times have the following behavior which is similar to the exactly solvable examples [5,6]. Therefore, in this limit and for both fields, we expect that if r h ≥ 1 then τ → 0, whereas if r h < 1 then τ → ∞.
From the results of Ref. [5] and our numerical results, we see that the fundamental frequencies of the Dirac and KG field behave in a different way as z increases. For the scalar Table 7 We show the numerical results which are slightly different from the ones obtained in Ref. [5]. In this field |Im(ω 0 )| decreases, whereas for the Dirac field |Im(ω 0 )| increases. As shown in Appendix A, for the massive Dirac field propagating in the black brane (3) and in the zero momentum limit, the imaginary part of its QNF is negative. Owing to the time dependence exp(−iωt) we assure that its QNM are classically stable, since the field decays in time [5,6]. 3 As far as we are aware the perturbations in the black brane (3) for κ = 0 have not been studied except the case z = 1 and D = 3 where analytical solutions are found [8]. The extension is natural and could be considered in the future.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authos' comment: This paper is a theoretical study of black branes. All relevant information is presented in the paper.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendix A Classical stability
In this appendix to obtain decoupled radial equations we follow the method of Ref. [6] where fermionic perturbations in the Lifshitz black brane (3) were studied for the first time. The aim of this is to show that the QNM for the massive Dirac field in the zero momentum limit are classically stable. In order to do that, as in Ref. [6], we transform (16) with κ = 0 into two Schrödinger-like differential equations. 3 See Appendix A for a detailed analytical proof of this fact for the Dirac field in the zero momentum limit.