Late-time tails, entropy aspects, and stability of black holes with anisotropic fluids

In this work we consider black holes surrounded by anisotropic fluids in four dimensions. We first study the causal structure of these solutions showing some similarities and differences with Reissner-Nordstr\"om-de Sitter black holes. In addition, we consider scalar perturbations on this background geometry and compute the corresponding quasinormal modes. Moreover, we discuss the late-time behavior of the perturbations finding an interesting new feature, i.e., the presence of a subdominant power-law tail term. Likewise, we compute the Bekenstein entropy bound and the first semiclassical correction to the black hole entropy using the brick wall method, showing their universality. Finally, we also discuss the thermodynamical stability of the model.


I. INTRODUCTION
Recently, the LIGO collaboration [1] [2] started the age of gravitational wave astronomy through the detection of a gravitational signal coming from the merger of two astrophysical black holes. Such signal was strong enough to permit the observation of the ringdown phase characterized by the so-called quasinormal modes (QNMs), which carry information of the structure of the spacetime itself. In addition, the study of QNMs spectra can bring a better understanding of the stability of a given black hole solution [3][4][5][6][7]. Moreover, this question can be addressed through the scattering of a scalar field in the fixed black hole background [8][9][10][11][12][13], which can be understood as a probe field to test the (in)stability of the black hole metric.
The QNMs and its spectrum are characterized, under appropriate boundary conditions, by a set of complex frequencies and encode the linear response of the black hole geometry to an external probe field with different spin weights. The time evolution of such probe fields is divided in three main stages: the initial burst in a short interval depending on the initial conditions, followed by the damping oscillation given by the QNMs and, at late-times, a power-law or exponential tails.
Another interesting subject that black holes bring is their thermodynamics. The similarity between classical thermodynamics and the laws governing the mechanics of black holes was well established by Bekenstein and Hawking [19,20] through the identification of black hole surface gravity and event horizon area with the temperature and entropy of a thermodynamical system, respectively. This fact led to the well known Bekenstein-Hawking formula, expressed in geometrical units. Based on this novel theory Bekenstein proposed the existence of an upper bound on the entropy of any system of energy E and dimension R given by S ≤ 2πER [21]. This equation is a consequence of the validity of the generalized second law (GSL) of black hole thermodynamics. Furthermore, in an effort to include quantum aspects in the gravitational theory describing a black hole, 't Hooft [22] proposed a semi-classical method to compute the corrections to the classical entropy formula (1). This technique known as the brickwall method consists in considering a thermal bath of scalar fields living outside the event horizon. The quantization of these fields via statistical mechanics partition function leads to quantum corrections to the black hole entropy. By carrying out this calculation on a Schwarzschild black hole 't Hooft showed that the first correction is proportional to the area, as expected, having a coefficient dependent on the proper distance from the horizon to the wall. Later calculations in other solutions showed that this first correction is the same in 4-dimensional geometries.
In this work we are interested in a solution of Einstein equations discovered by Kiselev [23], which describes a spherically symmetric black hole surrounded by an anisotropic fluid [24,25]. This constitutes a line-element derived from the solutions studied in [26], the so-called dirty black holes. Studies on its stability [14][15][16][17][18] and some aspects of its thermodynamical behavior have been implemented in the last years [27][28][29][30][31]. However, a detailed description of the causal structure of the spacetime, the late-time behavior of the scalar QNMs, and other aspects related to corrections to the entropy and thermodynamical stability are absent in the literature.
The paper is organized as follows, Section II presents the metric describing the family of black holes surrounded by anisotropic fluid and its main features. In Section III we present the causal structure of this spacetime. Also, the perturbative dynamics due to probe scalar field evolution is formulated and the QNMs spectrum and late-time tails are computed.
Section IV brings a study of some aspects of black hole thermodynamics including Bekenstein entropy bound, semiclassical corrections to entropy through t'Hooft brick wall method, and thermodynamical stability tested using specific heat and Hessian matrix criteria. Finally, in Section V some final comments are given.

II. BLACK HOLE SOLUTIONS
We are interested in a kind of dirty black hole whose line-element can be written as where dΩ 2 represents the metric of the 2-sphere and f (r) is given by [23] f where p r and p t represent the total energy-momentum tensor components T 11 and T 22 = T 33 , respectively. This non-zero anisotropy labels a non-quintessential fluid, different from what was stated in the first work which presented such a metric [23].
Furthermore, we can reinterpret the energy-momentum tensor of the solution as a sum of anisotropic fluids with different state parameters instead of considering a black hole surrounded by just one fluid component. By writing with c n and f n being constants, the energy-momentum tensor is linear in each 'charge' n, i.e., T = T f 1 +T f 2 +T f 3 +· · · . In such case by the proper choice of c n 's and f n 's we can easily have the charged black hole surrounded by a fluid as represented previously, meaning that the traditional components of charge and mass can be seen as fluid charges in the Kiselev picture [23].
Now the null energy condition imposes severe restrictions on the state parameter of the fluid ω f . By taking the condition of validity of the null energy statement [24,25] we have that the density gradient of the fluid is where m ′ represents the derivative of the position-dependent mass function m(r) defined as [25] 2 m(r) with K i and ω i being general coefficients and exponents of a Puiseux series. In our case we Thus, the energy condition is preserved whenever −1 ≤ ω f ≤ 0, and violated otherwise.
For this reason in this work we will study dynamical and thermodynamical aspects of the geometry within the range of validity of such condition.
In the next section we are going to characterize the causal structure of the family of solutions represented by the line-element (2) establishing the nature of the singularity and the horizons. In addition, we will check the late-time behavior of scalar QNMs in that geometry.

III. CAUSAL STRUCTURE AND PROBE SCALAR FIELD EVOLUTION
We are going to describe the causal structure for two different representative black hole solutions of the metric (2). We start by considering the behavior of the Kretschmann invariant given by where we have defined σ = 3w f + 1, p 1 = σ 4 + 2σ 3 + 5σ 2 + 4, p 2 = σ 2 + 3σ + 2 and p 3 = 3σ 2 + 7σ + 2. In the cases when w f ≤ 0 we have σ ≤ 1, so the Kretschmann invariant always diverges at r = 0 and is well behaved at the horizons and, thus, the line-element (2) has a physical singularity at the origin r = 0. In what follows, we are going to show that for two specific cases ω f = −1/2 and ω f = −2/3 with M > Q there is a range of parameters that represents a black hole with cosmological-like horizon r c , an event horizon r + , and Cauchy inner horizon r = r − covering the time-like singularity at r = 0. Such causal structure is very similar to the Reissner-Nordström-de Sitter black hole, except in the region beyond the cosmological-like horizon r > r c , where the spatial infinity (r → ∞) is light-like.
A. Black hole solution with w f = −1/2 Considering the line-element (2) with w f = −1/2 and the redefinition of the radial coordinate r = z 2 we have where the function H(z) is given in terms of three real roots z c > z + > z − denoting, respectively, the cosmological-like, event, and Cauchy horizons, and two real negative roots (z 1 , z 2 ). Thus, yields a tortoise coordinate given by which defines the usual double null system, U = t − z * and V = t + z * . Here the constants (α c , α + , α − , α 1 , α 2 ) are all positive definite and are given in terms of the horizons where the indices i and j denote the horizons (z c , z + , z − , z 1 , z 2 ).
We perform a detailed examination of the behavior of the black hole solution in the vicinity of each horizon in order to obtain the Kruskal-Szekeres extension to end up with the Penrose-Carter diagram of the entire manifold.
Near the cosmological-like horizon z = z c , the Kruskal-Szekeres coordinates U c and V c obey the following relation where the plus sign denotes the region z > z c and the negative sign corresponds to the region z < z c . Similarly, near the event horizon z + we have where the upper sign refers to z > z + and the lower sign refers to z < z + . Finally, for the region near the Cauchy horizon z ≈ z − , we have Introducing the Penrose coordinates T = 1 2 (Ũ +Ṽ ) and R = 1 2 (Ũ −Ṽ ) in each region covered by the relations (14 -16) withŨ = arctan(U) andṼ = arctan(V ), we compactified the coordinates. Furthermore, combining different overlaping coordinate patches it is possible to extend the metric through each horizon, thus, constructing the conformal diagram for the entire spacetime (10) in Fig. (1). Such diagram shows a causal structure very similar to that of a Reissner-Nordström-de Sitter black hole [32,33]. We observe an infinite sequence of structures featuring two outer horizons (event and cosmological-like), an inner Cauchy horizon, and a time-like singularity at the origin z = 0. However, the spatial infinity (z → ∞) in the black hole solution with w f = −1/2 displays a light-like structure, which is different from the Reissner-Nordström-de Sitter case, where the spatial infinity is space-like (see Fig.2 in [33]).
For an observer in region I crossing the event horizon and entering region III, we observe that the coordinate z is now time-like and the subsequent motion occurs with z decreasing.
However, after the observer crosses the Cauchy horizon, the coordinate z becomes space-like again, so it is possible for this observer to avoid the time-like singularity at z = 0 and emerge in another copy of region III.
The maximally extended black hole with w f = −2/3, and the conformal diagram is the same as in the case w f = −1/2, and can be obtained by performing the same steps as discussed here. The detailed calculation of the extension is given in the Appendix A.

B. Klein-Gordon equation
For a black hole spacetime as represented in Fig.1 the physical universe lies in region I, where we choose to integrate a scalar field that do not change the geometry.
In this domain the integration of the Klein-Gordon equation, ✷Φ = 0, will be affected by the definition of a tortoise coordinate system, dx = f −1 dr, (now in terms of r) used to fix the field propagation as ingoing plane waves crossing through the boundaries of x. In terms of this system the field equation turns to the usual simple form where Ψ represents the radial-temporal part of the Klein-Gordon field written as and V (r) plays the role of a potential for the scattered scalar waves given by In de Sitter spacetimes the tortoise coordinate places the cosmological horizon r → r c at the point x = ∞ and the event horizon of the black hole r → r + at x = −∞. This is also the case for dirty black holes with an anisotropic fluid as discussed in this paper. As a consequence, when using the above wave equation we will restrict the integration to the When studying the evolution of fields in fixed geometries, Eq.(17) establishes a master equation and for different fields (or spherical geometries) the proper V (r) must be taken.
The numerical integration in double null-coordinates for the calculus of the quasinormal modes is a well-establish method, which in general does not depend on the initial conditions.
Except for the initial burst of evolution, the quasinormal ringing phase that follows and the late-time behavior depend only on the geometry parameters. In terms of the null coordinate the Klein-Gordon equation takes the form or, written as a discrete equation, The boundary conditions in such system can be put in the form although discussions on the preservation of polar and radial symmetry (for the gravitational field) have presented Neumann boundary condition as the appropriate one.
After obtaining the field profile in time domain we can employ the Prony method [34] to acquire the quasinormal frequencies or, in the case of non-oscillatory profiles, linear regression. We will also use the WKB6 method [35][36][37] as a matter of comparison.

C. Late-time behavior and quasinormal modes
The late-time evolution of the probe scalar field brings two distinct behaviors depending on the fluid parameters c and ω f . In Fig.2 we can see different field profiles evolved from a similar initial burst as defined above. Depending on the fluid charge parameters we have an exponential decay or a power-law tail dominating the final stage of evolution. The fact comes surprisingly as a combination of two distinct behaviors already found in black holes with/without cosmological constant, being such final stage an exponential decay/power-law tail, respectively, for the Reissner-Nordström case. A second element present in the scalar field evolution of the above figures is the quasinormal modes, damped oscillations that arrive given the presence of a black hole potential barrier such as (19). In table I we list the fundamental modes for different values of fluid density. As expected, the influence of the fluid in the scalar field QNMs is very mild when its density is small (not detectable, e.g. for c ∼ 10 −6 ), no matter what the state parameter is. As c increases, the differences coming from several state parameters of the fluid increase as well. We can see that the quality factor, Q = Re(ω) − Im(ω) , decreases as we increase |ω f |. In fact, in a spacetime with an anisotropic fluid the scalar field oscillates better compared to a spacetime with cosmological constant: e. g. when M = 2Q = l = 10c/3c max = 1, we have Q = 3.30, 3.18, 3.03 and 2.88 for ω f = −1/2, −2/3, −5/6 and −1, respectively.
The results in the Table I were double checked with the WKB6 method [35]. The convergence of both calculations is as good as 0.1% for c/c max 0.5, where c max represents the maximum value of fluid density to which 3 horizons arise. Whenever the fluid density is high, higher divergences are found. This comes as no surprise as long as the WKB6 has a poor convergence for near extremal black holes.
For a large range of parameters we investigate the transitional behavior of the scalar field at late-times. Testing for the linear correlation of two different profiles written as we perform calculations for different state parameters going from ω f = −0.5 to ω f = −1.
The results are given in Table II. Observing the high values of linear correlation we state that both behaviors (exponential decay and power-law) are present in the final stage of the field evolution being one of them dominant.
We can see a small variation in the linear coefficients of the power-law series for −0.65 ω f −0.5 and an explosion after that, softening its behavior in the field composition Ψ| late times → C 1 t −a + C 2 e −αt . This makes the presence of this term subdominant in relation to the exponential decay series, which is prevalent for ω f −2/3.
This comes as an interesting result not stated until now in the available literature, e. g.
for RNdS geometries, the presence of a power-law tail term subdominant to the imaginary quasinormal mode in late-times in such spacetimes.
In the last two columns of the Table II we can see the quasinormal modes frequencies for a variety of ω f . The frequencies were obtained via Prony method with the same field profiles used in the late-time test. Again they were checked with WKB6 method with very good agreement in the results (maximum deviation of 0.1%).

IV. THERMODYNAMICS
In this section we are going to discuss some thermodynamical aspects of the dirty black holes under consideration.  First of all, we can rewrite the metric coefficient (3) in terms of the event horizon as In addition, using the metric (2) we can write the surface gravity as Both expressions will be useful in our next calculations.

A. Entropy Bound
Let us consider a particle in equatorial motion near a black hole. The constants of motion are given by corresponding to the energy and angular momentum of the particle, respectively. Since the energy conservation for a particle of mass m implies −m 2 = π µ π µ , using the metric (2) together with the metric coefficient (3) we can obtain a quadratic equation for the conserved energy E of the particle, whose solution becomes As the particle is approaching the black hole gradually, this process must stop when the proper distance from the body's center of mass to the black hole horizon equals the body's radius R, where r + + δ(R) represents the point of capture of the particle by the black hole. At this point the energy of the particle (30) can be evaluated and minimized with respect to the angular momentum of the particle. This results in J min = 0, such that In order to perform the integral (31), express δ in terms of R, and evaluate Eq.(32), we considered 3 cases, ω f = −1/2, −2/3, −5/6. To first order in δ the proper distance integral yields, , for ω f = −5/6 From the first law of thermodynamics we have that being A r the rationalized event horizon area A/4π and dM = E min , the change in the black hole mass due to the assimilation of the particle. Using Eqs. (27), (32), and (33) we obtain This result is independent of the black hole parameters and perfectly agrees with the universal bound found by Bekenstein [21].

B. Semiclassical corrections to black hole entropy
Following 't Hooft's brickwall method [22] we consider a thermal bath of scalar fields propagating just outside the horizon of a black hole background given by Eqs. (2) and (3).
The minimally coupled scalar field with mass µ satisfies Klein-Gordon equation, The idea is to quantize this field using the partition function of statistical mechanics, whose leading contribution comes from the classical solutions of the Euclidean action that yield the Bekenstein-Hawking formula. This scalar field will produce quantum corrections to the black hole entropy which can be calculated using the brickwall method. The 't Hooft method consists in introducing an ultraviolet cut-off near the event horizon such that Φ = 0 for r ≤ r + + ǫ. In addition, in order to eliminate infrared divergences another cut-off is introduced at a large distance from the horizon, L ≫ r + , where Φ = 0 for r ≥ L. By decomposing the scalar field as the radial part of Eq.(37) turns into where ℓ(ℓ + 1) is the variable separation constant. Then, using a WKB approximation for R(r) ∼ e iS(r) in Eq. (39), where S(r) is a rapidly varying phase, to leading order only the contribution from the first derivative of S is important. This contribution represents the radial wave number K ≡ S ′ , which can be obtained from the real part of Eq.(39) as In terms of this quantity the number of radial modes n r is quantized semiclassically as, Furthermore, the entropy of the system will be calculated from the Helmholtz free energy F of the thermal bath of scalar particles with temperature β −1 = κ/2π, where we made an integration by parts in the last step. Using Eqs. (40) and (41) and performing the integral in ℓ we obtain According to brickwall method we should study the contribution of this integral near the horizon. Thus, using Eq. (26) to write an approximate expression of the metric near the horizon and performing the integral in E we get where we rescaled some quantities as y = r/r + ,L = L/r + , andǭ = ǫ/r + . At this point it is convenient to consider different values of ω f separately. We should notice that the divergent contribution of the integral to the Helmholtz energy comes from its lower limit. Thus, the leading divergent term F ǫ is given by The corresponding entropy S ǫ = β 2 ∂F ∂β , then, becomes for ω f = −1/2 , for ω f = −5/6 (46) We can express our results in terms of the proper thickness α defined as To first order this expression can give us a relation between ǫ and α for the values of ω f considered here, Replacing these values and the corresponding expressions for the surface gravity (27) in Eq.(46) we finally obtain in the three cases, or in terms of the black hole horizon area A = 4πr 2 + , This expression is the same correction found by 't Hooft and other authors for 4-dimensional black holes, a fact that shows its universality.

C. Thermodynamical stability
In order to see the influence of the fluid on the stability from a thermodynamical point of view the first step is to analyze the specific heat, which in our case becomes The plot of the specific heat for different values of ω f , displayed in the left panel of Fig.3, shows the rich structure of the geometry already noticed in the literature [27][28][29][30][31]. There are positive (stable) and negative (unstable) regions alternating with each other. These regions are separated by several points that signal first order phase transitions where C = 0 and also second order transitions whenever C becomes infinite. However, the sign of the specific heat is not enough to ensure stability. One additional criterion to verify the existence of critical points comes from the Hessian matrix of the Helmholtz free energy F related to the the black hole [42] where C is the conjugate quantity to the "charge" c related to the presence of the anisotropic fluid given by (54) Using the entropy S = πr 2 + and the temperature of the black hole T = κ/2π with κ given in Eq. (27) we find that With all this information we can calculate the determinant of the Hessian matrix. However, this determinant vanishes, what means that one of the eigenvalues of the matrix is zero.
The other eigenvalue corresponds to the trace Tr(H) of the Hessian matrix (53). Then, a necessary criterion for the model to be stable is the positivity of this quantity, i.e., Tr(H) ≥ 0.
We plotted this trace in the right panel of Fig.3. We observe that, in fact, there are regions where Tr(H) ≥ 0 for the values of ω f considered along this work. Moreover, in the left panel of Fig.4 we see that small black holes fulfill the stability criterion independent of the value of c, whose influence is only visible for bigger r + . A curious fact is that for ω f = −2/3 the trace of the Hessian matrix does not depend on c. In addition, the effect of the charge on the stability criterion can be seen in the right panel of Fig.4, small charge black holes have shorter regions of instability. Therefore, with this analysis we see that it is possible to have phase transitions for different values of the black hole and anisotropic fluid parameters.

V. FINAL REMARKS
In this paper we investigate charged black hole spacetimes surrounded by anisotropic fluids. We firstly obtained the conformal structure of the entire manifold showing that its Penrose-Carter diagram is similar to Reissner-Nordström-dS spacetime, i.e., there is a cosmological-like horizon, an event horizon, and inner Cauchy horizon. In addition, there is a time-like singularity at the origin that could be avoided by an observer crossing the inner horizon. The novelty in the spacetimes considered in the present work is the light-like structure beyond the cosmological-like horizon differently from the RN-dS black hole where this region presents a space-like structure.
Having established the causal structure of the black hole spacetime we evolve the scalar field between the event and cosmological-like horizons obtaining two interesting features.
The first one is that the late-time behavior is dominated firstly by a power-law term for small state parameter of the fluid |ω f | and, afterwards, by an exponential decay (purely imaginary quasi normal mode) for higher |ω f |. For these geometries the presence of a power-law term in the final stage comes as an interesting new result never reported before even in de Sitter black hole spacetimes where this phenomenon is also present. The second one concerns the quasinormal modes obtained in Section III. They provide the spectrum of oscillation of the black hole when perturbed by a scalar field. We show that they are very similar for different state parameters when the fluid density is small being hugely influenced when it becomes large. When varying the state parameter, the oscillations have both imaginary and real part diminished as we increase |ω f |.
Regarding the thermodynamical calculations, we considered an arbitrary particle of proper energy E in equatorial motion and captured by these black holes surrounded by anisotropic fluids. Our result shows that these geometries yield the universal bound for the entropy of the falling system originally found by Bekenstein [21]. In addition, we also considered a thermal bath of scalar fields propagating outside the event horizon of these black holes in order to find the semiclassical corrections to their entropy. Following 't Hooft's brickwall method we found the same kind of correction corresponding to 4-dimensional black holes showing the universality of this result [22]. Finally, we also analyzed the thermodynamical stability looking at the specific heat of the black hole. As an additional criterion to ensure the presence of critical points, we also calculated the trace of the Hessian matrix of the Helmholtz free energy. In this way we showed that phase transitions of first and second order are possible for different values of the black hole and anisotropic fluid parameters.