Instability of BTZ black holes in parity-even massive gravity

We investigate the linearized equations of motion in three dimensional Born-Infeld gravity theory. Motivated by this model, we calculate the quasinormal modes of BTZ black hole solutions of parity-even gravity theories in three dimensions by using numerical methods. The results are classified in three families and they are so accurate such that it allows us to propose analytical form for the quasinormal frequencies. We find new quasinormal modes which have been missing in the literature of the analytical studies of three dimensional massive gravitons. These new modes do not have the known tower structure and they are purely imaginary for any value of the angular momenta. Considering the complete set of the quasinormal modes, we show that the BTZ black hole solutions are unstable for any value of the parameters of the theory. We confirm our numerical results by computing the new eigenmodes analytically at zero angular momentum.


Introduction
It is well known that three dimensional Einstein gravity has no propagating degrees of freedom in the bulk. In spite of that, Topological Massive Gravity (TMG) [1], New Massive Gravity (NMG) [2] and its generalization [3,4] and Born-Infeld (BI) gravity [5] have provided a framework to study three dimensional gravity with higher curvature terms which allow massive graviton excitations. This setup may eventually helps us to understand the quantum gravity in three dimensions. Investigating the quasinormal modes (QNM) of the classical solutions of these theories may teach us about the characteristic features of the corresponding geometries. For reviews on QNMs of black holes see e.g. [6,7].
From a holographic point of view, the AdS 3 /CFT 2 correspondence has been well studied in last three decades initiated by Brown and Henneaux [8]. On the other hand, it is known that adding higher derivative curvature terms in the classical bulk gravity theory is dual to going away from strongly coupled regime of the boundary gauge theory. Therefore, the BI-action, which can be considered as an infinite set of higher derivative terms, is a convenient framework to study the effects of the higher derivatives in the bulk theory, and maybe, to investigate the intermediate coupling behaviour, all the way to weakly coupled regime in the dual boundary theory.
In this note we focus on three dimensional BI-gravity and its QNMs structure. We show that the fourth order linearized equations of motion can be factorized into two second order differential equations, corresponding to massless and massive graviton excitations, and has similar form as the QNM equations of the NMG studied in [9]. This motivates us to study the QNMs in a more general framework where the results can be extended to other parity-even theories.
The massive mode excitation in TMG satisfies a first order differential equation which allows an analytical approach to find the quasinormal frequencies [10]. In [11] it was shown that the second order differential equation of the perturbations in parity-even theories could be decomposed into two first order differential equations. Then following the same approach as [10] the analytical solutions are computed, which are obviously solutions to the original equation. Nevertheless, there is no guarantee that calculated modes are the full solutions of the original second order differential equation. To cover this missing point in the literature of the QNMs in parity-even 3D gravities, we employ a numerical approach to compute the quasinormal modes of the original equation.
Our results show new features of the excitations in three dimensional parity-even massive gravities. While our numerical QNMs include the analytical ones, we found new modes which do not appear in the analytical studies. Due to existence of the new modes, we show that static BTZ black hole solutions of these theories are unstable for any value of the parameters of the theory and the background. On top of that, the new modes we found are unique and there is no tower of them in the QNM structure which means they do not belong to any Christmas tree structure of the modes. In zero angular momentum we could find analytical solution corresponding to the new mode which, confirms our numerical analysis.
The organization of the paper is as follows. In the next section, 2 we will shortly review the BI-gravity and its relation with NMG and its extensions to O(R 3 ). We show a regime of the parameters which admits the BTZ black hole solutions and we study its boundary theory. In section 3, we find the linearized equation of the metric perturbations and investigate the mass of graviton excitations. Section 4 is devoted to analytical attempts to find the massive QNMs. The heart of our work is section 5 in which we use numerical tools to compute the QNMs. In this section we start with vanishing angular momentum and then we turn it on for the more general situation. For the former case we support our result by using analytical investigation. We close the paper by a summary and outlook in section 6. For completeness appendixes A and B contain some technical details of the computation of boundary stress tensor and the coupled QNM equations at finite angular momentum, respectively.

3D-BI gravity
In this section we provide a review on BI gravity in three dimension and, its relation with NMG theory [2] and its extensions [3,4]. Although the equations of motion of this theory have a diversity of solutions [12][13][14][15], we will focus on the BTZ geometry and investigate its physical properties. The BI extension of NMG introduced in [5] as, where G µν = R µν − 1 2 Rg µν is the Einstein tensor, m 2 is a positive definite dimension-full parameter, σ = ±1 is a sign coefficient, λ is a parameter related to the cosmological constant and κ 3 is related to the three dimensional Newton constant as κ 2 3 = 16πG 3 . Note that λ = 0 is excluded from the parameter space since at this point the field redefinition g µν → g µν − σ m 2 G µν makes the theory trivial [16]. Expanding this action in terms of 1/m 2 , which is practically a derivative expansion, has the following interesting feature The zeroth order expansion gives the Einstein-Hilbert action with a cosmological constant Λ = 2m 2 (1−λ) σ , and to first order it reproduces the NMG action [2] which is a minimal parity preserving massive gravity in three dimension. On top of that, it is straightforward to show that expansion of the BI action (1) up to second order reproduces the extension of NMG theory proposed in [3] by demanding the existence of a holographic c-theorem. This is a fascinating agreement between two independent approaches of investigating higher derivative corrections to Einstein-Hilbert action in three dimensional gravity. While direct extensions of NMG to higher order seems to be an extremely difficult task, one may expand the BI action eq. (1) to any order of derivative curvature invariants. One main difference between BI theory and any truncated version of that appears in their vacuum solutions. Although BI has a single vacuum solution, truncated models may have several vacua. 1 We would like to study the quasinormal modes of the BTZ black holes of the BI theory by solving the linearized equations of motion. The following formulas are useful for analytical calculation, where A is a general invertible matrix and A.A −1 = 1. Looking at action (1) one may identify a matrix as A µ ν = δ µ ν + σ m 2 G µ ν , and its variation is given by δA µ ν = σ m 2 δG µ ν . In order to write the equations of motion it is convenient to define a tensor, and an operator, Using the following relation for variation of Einstein tensor, one may show that the equations of motion have the following compact form, We are interested in BTZ geometries which are examples of locally maximally symmetric spacetimes defined by 1 We thank B. Tekin for bringing this point to our attention.
where ± corresponds to the dS/AdS respectively, ℓ is the radius of dS/AdS and for flat spacetime the right hand side vanishes. From now on, we focus on Anti-de Sitter solutions and we work with minus sign in eq. (8). For AdS geometries one may show that the tensor defined in eq. (4) simplifies to Plugging the above results in the equation of motion (7) one may obtain a relation between the AdS radius ℓ and the parameters of the theory as This shows that for given parameters of the theory, there is a unique vacuum solution. Expanding this equation to second order in 1/m 2 , one may find, which is the corresponding relation for vacua geometries in NMG theory [2]. Note that in order to have locally AdS solution, satisfying eq. (10) , there are two choices for parameters of the theory with m 2 > 0, namely So far we have reviewed some features of the 3D-BI theory under conditions of having AdS solution. According to the holography conjecture, one may expect that there is a boundary theory dual to an asymptotically AdS background. The gravity side of the theory is a resummation of infinite number of higher derivative corrections to the Einstein-Hilbert action. Therefore the field theory side might be related to weakly coupled regime of a dual gauge theory. In the following subsection we investigate this boundary theory by computing the corresponding central charge and the stress tensor.

Boundary theory data
Boundary theory of AdS 3 geometries has been investigated in last 3 decades initiated by Brown and Henneaux [8]. In this section we make some comments from this perspective on the boundary theory of BTZ black holes in three dimensional BI-gravity. In three dimensions the central charge for higher derivative gravities can be obtained by using the Wald formula [17,18] where G 3 is three dimensional Newton constant and L is the Lagrangian density. Since the BI action (1) is introduced as a function of Einstein tensor we use the following chain rule, to express the derivative in eq. (12) in terms of derivative with respect to Einstein tensor. By employing eq. (3) and the tensor defined in eq. (4) one may show that where in the second equality we plug the conditions for maximally symmetric solution eqs. (9) and (10). Using these formulas the corresponding central charge simplifies to This is nothing but the scaling of Brown-Henneaux central charge [8] computed for Einstein theory. This expression has been found in [19,20] and it matches the Weyl anomaly computation, but here we could rederive it with less difficulty. From holographic point of view one may expect that the central charge (15) corresponds to the zero energy of a dual (conformal) field theory. Unitarity of the theory requires that this central charge should be non-negative which leads to, On the other hand, we found a regime of the parameters of the BI action (1) in which the BTZ geometries are solutions to the equations of motion in eq. (11). Now by comparing two conditions in eqs. (11) and (16), one can see that the unitarity of the boundary theory fixes the sign parameter σ and restricts the regime of λ, namely 2 In order to find the correct energy of a solution in a theory one may follow the Hamilton-Jacobi analysis [21]. The first step in this approach is calculating the correct boundary terms so that the variational principle is well-posed. For general metric backgrounds boundary terms are known for a few number of models such as Lovelock theories [22]. But for a symmetric background the boundary terms can be computed in a systematic way. In appendix A we show that for a maximally symmetric geometry in BI theory the boundary term is given by which leads to the Brown-York stress tensor as in whichĥ ab is the boundary metric and K ab is the extrinsic curvature. Interestingly enough, comparing this stress tensor with the results for Einstein theory one see the same scaling factor (σλ) as found for the central charge in eq. (15). Therefore, the energy density computed from the time-time component of the energy-momentum tensor (19) does not modify the regime of the parameters.
Similar to Einstein theory we can compute the central charge from the energy-momentum tensor given eq. (19), which obviously leads to eq. (15), and this could be considered as a crosscheck for our analysis. As a side remark we would like to note that, similar to our discussion after eq. (10), the central charge and Brown-York stress tensor for NMG theory and its nth order extended versions can be found by expanding our results for BI action to O(1/m 2 ) and O(1/m 2n ) respectively.

Linearized equations of motion
In this section we linearize equations of motion (7), derived in previous section for BI gravity, and compare its expansion in 1/m 2 with the linearization of NMG and extended NMG. We consider the metric as g µν =ḡ µν + h µν , in whichḡ µν is the background spacetime and h µν is the perturbation on top of that. Using the rules in eq. (3) and giving some massage, the equations of motion for perturbation is found to be, The pre-factor in this equation plays a crucial role later when we want to compare with the results for NMG theories by expanding eq. (20) in 1/m 2 and therefore, we keep it to the end. The trace of above equation has a simple expression, On the other hand, the Harmonic gauge ∇ µh µν = 0 with the aid of linearized equation of motion (21), leads to the so called transverse traceless (TT) gauge for original metric perturbations, ∇ µ h µν = h = 0. In this gauge, the linearized equation of motion simplifies to, which, equivalently, can be factorized as, Therefore the equation of motion for linear perturbation splits into two parts. Similar structure has been found for TMG [10] and NMG [9] theories. The second parenthesis in eq. (23) is related to a massless graviton in BTZ geometry while, the first one corresponds to a massive graviton with the mass Note that for σ = −1, which corresponds to non-unitary boundary theory (17), there is a critical point at mℓ = 1 where the mass of graviton, the parameter λ and the central charge of the boundary theory (15) vanish but, at this critical point the BI-model becomes trivial. However, for σ = −1 and mℓ = 1 the mass of gravioton satisfies the Breitenlohner Freedman bound m 2 > m 2 BF = −1/ℓ 2 [23]. On the other hand, σ = 1, which gurantees the unitarity of the boundary theory, leads to BF-bound violation, m 2 < m 2 BF . This indicates the "bulk-boundary clash" [24] of BI-model mentioned in [25] (for a review see [26]).
At this point, we would like to discuss about how one can get the NMG data from our results. Expanding eq. (23), including the pre-factor, to the first order in 1/m 2 leads to which is the equation for perturbations in NMG theory [9] and apparently the pre-factor has an important role in the first parenthesis. One may calculate the graviton mass for NMG theory by looking at the first term in eq. (25) which is given by Now the expansion of (23) to second order in 1/m 2 leads to Note that in the second line after factoring out the coefficient of box operator, we expand the remaining terms in 1/m 2 up to the second order. Using this linearized equation the mass of graviton is given by By following this procedure of expansion to O(1/m 2n ) one can show that expression for the mass of graviton has an interesting structure where c n 's are numerical coefficients and their absolute value are less than one, |c n | < 1. Note that, apart from the first two terms in eq. (29), there is only one other term, which is proportional to 1/(m 2 ℓ 2 ) n and at the limit n → ∞, this term would disappear if we demand m ℓ > 1. Let us emphasize that, all the potentially intermediate terms in the expansion (29) disappear due to the similar factoring procedure used in eq. (27). This is a fascinating result which shows how the truncation of the expansion in 1/m 2 at any order, which leads to three terms in right hand side of eq. (29), would be different than keeping all the terms as in BI theory, which corresponds to eq. (24) and is in agreement with the first two terms in eq. (29).

Analytical analysis of QNM's
In this section we are following [10] to investigate the quasinormal frequencies of the BTZ black hole in BI theory. The massive mode equation given in eq. (23) fits in a more general equation namely, where parameter a depends on the details of the theory and the background data. 3 In this way, we can compare the massive QNM's of BTZ black holes in different theories. Specially we want to investigate the effect of truncation of BI theory at a given order of expansion in 1/m 2 . In fact in three dimension one may further decompose eq. (30) to two first order equations. To this aim, consider the following operator with a constant M. Using this operator it is easy to show that In TT-gauge this relation is equivalent to eq. (30) by setting M = √ 3 − a ℓ 2 . Therefore every solution of the equation is a solution of the second order eq. (30). It is easy to show that the TT-gauge condition could be deduced from eq. (33). But one must note that there is a critical point where M = 0 at which, logarithmic solutions may appear. We do not consider this case in our analysis while it has been studied in various cases [27][28][29].
In order to solve eq. (33) one may consider the global coordinate system for BTZ black hole, such that the left/right moving solutions are given by [10]  The ingoing boundary condition at the horizon selects the positive signs in h r,l and therefore the least damped quasinormal frequencies are , where periodicity along φ leads to discrete angular momenta k.
Introducing the operator (L −1L−1 ), where L k are Killing vectors associated with background BTZ solution, one may find an infinite tower of quasinormal modes, and the corresponding frequencies with overtone number n, It is easy to see that M > 1 leads to stable frequencies. There is an important point regarding the solutions (35) and (36). As mentioned in [9], these two modes are not orthogonal, so they can not be interpreted as independent massive modes in the theory and they might be some missing solutions to the linearized eq. (30). We will come back to this point in next section. By comparing eq. (30) with the massive mode equation in different theories one may find the relation between the auxiliary positive parameter M and the parameters of the theory. For NMG, M 2 = −σ m 2 + 1 2ℓ 2 while for the generalized NMG it will be given by M 2 = −σ m 2 − 3 8ℓ 4 m 2 . The general structure for the truncated BI-action (generalized NMG to higher orders) can be found as On the other hand for the BI-gravity, by comparing eqs. (23) and (30), it is easy to see which, could be found by taking the limit n → ∞ from the results for extended NMG theory (41) by demanding m ℓ > 1.
At the end of this section we would like to point out the main linearized equation that one has to solve for the BI theory and other parity even massive gravity models is a second order differential equation given in (30) while the analytical investigation is based on factorized first order equations (33). Although it is clear that the solutions of to the first order equations (33) are solutions to the second order one eq. (30), it is not obvious whether it works other way around. In other words, the set of solutions to the later could be larger than the former. Of course, the same argument applied when we factorized the fourth order differential equations to one massless and one massive excitations in eq. (23). But, in this work we just focus on the massive mode and not the full linearized equation. In next section we use numerical tools to tackle linearized equation (30) with given boundary conditions corresponding to the quasinormal modes and we will compare the results of these two approaches.

Numerical analysis of QNM's
In this section, by using a numerical method, we compute the quasinomal frequencies of BTZ black holes as solutions of BI theory or any 3D theory of gravity which has a massive mode corresponding to eq. (30). To this aim, we consider the static BTZ black hole in Eddington-Finkelstein coordinate where ℓ is the radius of AdS 3 spacetime, and the asymptotic region is r → 0. Linearizing equation of motion on BTZ geometry in this coordinate has a twofold advantages. First, the ingoing boundary condition at the event horizon is automatically satisfied. Second, the vanishing source boundary condition at the asymptotic boundary (which is the definition of QNM) is easy to impose as we will discuss.
It is common to use the plane wave Ansatz for the fluctuations and work in momentum space, h µν (r, t, φ) = Z µν (r)e −i ω r 0 t+i k r 0 φ . The transverse traceless gauge, not only reduces the linearized equations to eq. (23) but also imposes four constraints on the components of h µν which are useful to simplify the QNM equations.

k = 0
We start with zero angular momentum k = 0 which leads to simpler equations. It is easy to show that in Eddington-Finkelstein coordinate there are two sets of metric perturbations which are coupled among each other and decoupled from the other ones. These two sets are (h tφ , h rφ ) and (h rr , h tr , h tt , h φφ ) which are orthogonal to each other and contain all of the metric perturbations. But using the constraints of the TT-gauge, one can show that the linearized equation (30) will reduce to two decoupled equations for metric perturbations h rφ and h φ . At the practical level, the following redefinitions are natural to have proper boundary behaviour, It is easy to show that using the following redefinitions reduce the equations in terms of new functions as Note that the same parameter M defined in eq. (45) has been used in the analytical studies introduced in eq. (31). Interestingly enough, not only M is the only parameter in both equations but also, there is no temperature dependency since we define the frequency and angular momentum in r 0 unit. Also, it is easy to see that in both equations the near boundary solution to the eqs. (46) and (47) is From this expansion it is obvious that we are interested in the cases in which parameter M is positive which means aℓ 2 < 3. In principle, for 3D massive gravity theories the parameter M could be related to the mass of the graviton in BTZ background and the parameters of the theory. Note that M = 0 is not in the domain of our interest since it corresponds to the logarithmic cases and M = 1 is associated with the massless graviton mode. According to the definition of quasinormal modes in asymptotically AdS backgrounds, the mode should be source free (A i = 0) and should satisfy the ingoing boundary condition at the horizon. On the other hand, according to the holography dictionary, the QNMs correspond to the poles of the retarded Green's function of the dual operators [30]. Therefore, our studies may shed some light on a field theory dual to the BI-gravity which might be considered as an infinite summation of the higher derivative corrections to the 3D Einstein-Hilbert action.
For our numerical calculation, we use Chebyshev spectral method to discretize the radial coordinate between the horizon (u = 1) and the asymptotic boundary (u = 0) and solve the eigenvalue eqs. (46) and (47) to find the QNM's. This method has been widely used in QNM literature and implementing the method practically has been reviewed in [31]. We should emphasize that by using this numerical method we are able to find the lowest QNMs such that, the higher the number of grid points the more number of lowest QNMs one could compute with higher accuracy. Nonetheless, accurate results may help us to propose analytical structure for the frequencies.
Here we summarize our numerical results with high precision (in our numerical calculation we use the number of significant decimal digits between 100 to 1000) and large number of grid points along the radial coordinate (between 100 to 250). The results are so accurate which shows following analytic forms: 1. The frequencies for the Z 1 are classified in two sets, Obviously, the first two sets in Z 2 channel coincide with the two sets of Z 1 channel and, the first set is in perfect agreement with the analytical expressions given in eq. (40) for zero angular momenta. But here we find the third set in Z 2 channel with a single mode which has been missing in the analytical investigations. On top of that, the degeneracy between the first and second sets for some modes was not apparent in our analytical results. Note that the dependency of the modes on the parameter M is completely different in the third set which motivate us to put this unique mode in a separated set. With these results some comments are in order: • Since the largest power of ω in equation eq. (46)/(47) is two/three, one may expect that there are two/three sets of QNM's in Z 1/2 channel a priory.
• For 0 < M < 1, the zeroth mode in the first family, ω (i,1) 0 , is unstable and all other modes are stable. Figure 1: The accuracy of the first three numerical QNM frequencies from the first set in comparison with analytical formula for M = 0.5 for Z 1 (left panel) and Z 2 (right panel). The horizontal axis shows the number of grid points along the radial coordinate and green, red and blue dots correspond to the n = 0, 1, 2 QNMs respectively. It is interesing that in this case the first mode in the Z 1 family has higher accuracy than the zeroth mode.
• For M > 1, the single mode in the third family ω (2,3) is unstable and all other modes are stable.
Therefore, while the Z 1 channel is stable for M > 1, the Z 2 channel is always unstable for any value of the parameter M, in the regime of our interest. Note that this is true for any 3D gravity theory which leads to the linearized equation in TT-gauge, (30) with a massive mode. This is a fascinating result which we found by implementing the numerical analysis for the first time. In fig. 1 we show the accuracy of our numerical results comparing with the analytical QNMs given in eqs. (49) and (50) for the first three modes in the first set.
By looking at the Eigen-vectors of the QNMs, we can justify our results and show that there is no pathology in the new modes that we found. In fig. 2 we show the real and imaginary parts of the Z 2 Eigen-vectors for the zeroth mode of the first set of QNMs (left panels) and the single mode in the third set of QNMs (right panels) at zero angular momentum for M = 0.5 (upper panels) and M = 1.1 (lower panels). According to our results given in eq. (58), for M = 0.5 the unstable mode is ω (1) 0 but for M = 1.1 the instability is due to the ω (3) . The expected behaviour of the Eigen-vectors at the horizon (regularity condition) and at the asymptotic region (Z ∼ u 2M ) could strongly support our claim about instability of BTZ solutions in parityeven 3D-gravity models.
It is very interesting that one can check analytically if the third set frequency corresponds to true QNM eigenmode. It is straightforward to see that plugging ω = −i (1 − M 2 ) in eq. (47) reduces to 4 As mentioned earlier, QNM Eigen-mode is defined with source free boundary condition at the asymptotic boundary which means A = 0. Also since the eq. (51) is a linear equation, one may rescale the Eigen-mode and the QNM mode is given by, We would like to point that one can find the Eigen-modes of other QNM frequencies in the first and second channels by plugging eqs.
which are fascinating results supporting our numerical analysis. While we found these expressions one by one, it might be possible to find a general analytical solution for nth frequency.

k = 0
Now let's turn to the finite angular momentum, k = 0. In this case, by using the TT-gauge constraints, one can show that the linearized equation (30) will reduce to two coupled equations for the metric perturbations h rφ and h φφ . Hiring the same redefinitions introduced in eqs. (44) and (45) one can find the linearized equations in new functions which are too lengthy and given in appendix B. Similar to the k = 0 case, the only free parameter in the equations is M and there is no explicit temperature dependency. But at finite angular momentum the equations are coupled. The near boundary analysis shows the same expansion of the Z i , given in eq. (48). Again, our numerical results are as accurate such that the following analytical forms of the frequencies, which are classified in three sets, can be proposed The third set has completely different dependency not only in the parameter M but also in angular momentum k such that they are purely imaginary for any value of the parameters in the domain of our interest. Note that there is a degeneracy for any value of angular momentum k and to our knowledge, this is the first example in asymptotically AdS black holes that there are unique modes which do not belong to any tower of the QNMs. On the other hand, for any value of parameter M and finite angular momentum k = 0, there are two unstable modes either in the first set or in the third set of the frequencies given in eq. (58). In this case all the metric fluctuations are coupled which means the static BTZ black hole for any theory, which has massive mode satisfying the general equation (30), is dynamically unstable with respect to any non zero angular momentum perturbation.
While the structure of the first two sets is similar to what we found from analytical studies in eq. (40), there is an obvious difference in the real part of the modes. Our numerical modes can move in both left and right directions since they have both signs in their real part. This is the second new results we found using the numerical approach to solve the second order linearized equation (30). Unfortunately, in contrast to the previous subsection, we could not find the analytical Eigen-modes of linearized coupled equations at non-zero angular momentum. But the numerical results are as accurate as the zero angular momentum case.

Discussion and Outlook
The BI theory of gravity is interesting for different reasons. It could be considered as a resummation of infinite number of higher derivative corrections to the Einstein-Hilbert action while by expanding the BI action in parameter 1/m 2 , one may get NMG theory or its extended versions to any order. Therefore, it could provide a proper setup to investigate the effects of higher derivative corrections to the physical properties that one would like to study.
In this paper we found regimes of the parameters of the BI theory in which static BTZ black holes are solutions to the equations of motion. Then we investigate the quasinormal modes of these solutions to see if there is a range of parameters which leads to stable black holes.
We use the known analytical approach [11] to decompose the second order differential equation (30), which specifies the quasinormal modes, to the first order differential equations (33), and the results for left/right moving modes are given in eq. (40). While the solutions to the later equations clearly satisfy the former one, there is no guarantee that all the solutions to the original equation could be found in this way. Therefore, we use the numerical tools to tackle the equation (30) with ingoing boundary condition at the horizon and source free boundary condition at the asymptotic region.
Our numerical analysis is based on a general form of the linearized equation (30) in TTgauge which let us to make some comments on any 3D parity even massive gravity theories including NMG [2], its extensions [3,4] and ZDG [32,33]. On the other hand, the numerical results are accurate which let us to propose analytic expression for the QNM frequencies given in eq. (50), and they are promising to have analytical solutions for the second order linearized equation. By plugging the QNM frequencies eqs. (49) and (50) in to the eqs. (46) and (47) we found analytical solutions for some modes in vanishing angular momentum which confirms our numerical results. Finding a general form of the analytical eingenmodes could be an interesting open question which needs further investigations. We didn't study the quasinormal frequencies of the solutions at the critical point at which Logarithmic modes are new solutions to the equations. This has been studied from different point of views in various cases in [28,29,34,35].
As we know, BTZ black hole is a solution to the three dimensional Einstein-Hilbert and its massive extensions, e.g. TMG, NMG and BI-model. There is a common feature in different dimensions that the solutions to the Einstein-Hilbert action are solutions to its extended massive gravity theories as well. The stability of these solutions in massive gravity models have been studied in various cases. For example, in NMG theory it was shown that the warped AdS geometry is unstable in all range of parameters [36] while in four dimensional gravity theories instability of Schwarzschild black holes and Kerr geometries have been studied in [37,38]. Therefore this kind of instability may have the same origin and investigating this point could help us to understand the physics of massive spin-2 excitations.
What we found shows two sets of quasinormal frequencies, with a structure similar to the analytical results, and new unique modes which do not belong to any tower of the QNMs in asymptotically AdS backgrounds. On top of that, the dependency of the modes on the only parameter of the linearized equations is such that it leads to unstable quasinormal modes in any regime of the parameters in any 3D massive theory which has a linearized equation given in eq. (30).
Although we show instability of BTZ black holes in parity-even massive theories, we did not study what it decays to. The following steps has to be taken to answer this question. It is known that solutions to three dimensional gravity theories with higher derivative terms are richer than Einstein gravity and have not been completely classified [12][13][14][15]. Therefore as the first step, one has to find all the solutions (including time dependent geometries) of the equations of motion of the corresponding model. Then, one may study the full time evolution of the perturbed BTZ within the theory. This is a very interesting open question which is beyond the scope of this paper and needs further studies.
Last but not least, we would like to compare our results with the QNMs of 5D gravity theories with negative cosmological constant and at the presence of the higher derivative curvature terms in the action studied in [39,40]. It was shown that at large but finite coupling, the branches of the QNMs become more dense and lift up towards the real axis, which may lead to forming a branch cut in the limit of zero coupling. In contrast, our analysis reveals a different structure either for the full BI-action or for its truncated expansion at any order in 1/m 2 , in which the QNMs have finite distances in the frequency complex plane and there is no evidence of forming any branch cut. Comparing these two structures in three and five dimensional theory of gravity in asymptotically AdS geometries needs more investigations.

A Stress tensor for BI gravity
In order to find the equation of motion (7) one has to apply integration by part. The total derivative of this integration leads to a boundary term which has the following well-known form where the expression for J α is given by, Although the expression seems too complicated and finding the proper boundary terms for general metric background is impossible, for a maximally symmetric spacetime background eq. (60) reduces to J a = λσ(∇ α δg β β − ∇ β δg α β ) The terms in the parenthesis of this surface integral is nothing but what we get from variation of the Einstein-Hilbert action and the pre-factor σλ is the fingerprint of the BI theory. It is well-known that using these terms one can obtain the Gibbons-Hawking-York boundary term and the Brown-York stress tensor, (see e.g. [43] ). Therefore, the boundary terms and stress tensor for the three dimensional BI theory of gravity on maximally symmetric backgrounds is a multiplication of the Einstein-Hilbert ones with σλ, namely To compute the free energy of a solution one should evaluate the full on-shell action which is L T = L BI + L boundary + L ct , where L ct is the boundary local counterterm [44]. By comparing the total action eq. (63) with the three dimensional Einstein gravity theory with negative cosmological constant, it is easy to show that, for a locally AdS geometry, the free energy has the same pre-factor structure as we found for the central charge and the boundary stress tensor given in eqs. (15) and (19) respectively.

B Non-vanishing angular momenta
The coupled linearized equations with nonzero angular momentum are more involved. In this appendix we show the explicit form of these equations for completeness. As we expected, at finite angular momentum, these are coupled equations given by