Quasinormal modes, greybody factors and thermodynamics of four dimensional AdS black holes in Critical Gravity

Analog to the ringing of a bell, when a black hole (BH) is perturbed, due for example to an infalling object or some disturbance in its surroundings, it responds via certain characteristic modes, denoted as quasinormal modes (QNMs) and the associated quasinormal frequencies (QNFs), that govern the time evolution of the initial perturbation [1]. These QNMs have attracted a high interest due to the detection of gravitational waves [2, 3] (for reviews, see Refs. [4, 5, 6]). As was shown firstly by Regge and Wheeler [7], as well as by Vishveshwara [8], at the moment of exploring Schwarzschild BHs, these modes depend only on parameters of these configurations, remaining independent of the nature of the perturbation and being a distinctive ’fingerprint’ of the BH. The above has allowed the exploration of QNMs considering massless topological BHs [9], topological AdS BHs [10, 11], regular rotating BHs [12] and with non-standard asymptotic behavior [13].

The consistent emergence and prominence of these distinctive perturbations have undergone rigorous testing, allowing us to study QNMs and QNFs deeply. For example, at linearized analyses [8, 14, 15], where the fields are approached as perturbations in the singular BH spacetime, numerical computations encompassing scenarios such as BH collisions [16, 17] or stellar collapse [18, 19], as well as the estimation of BH parameters, such as the mass, angular momentum, and charge, through gravitational waves [1, 2, 3].

The evolution of small perturbations within BHs, generally unfolds in three key stages: the initial outburst, the ringing of QNMs, and the eventual power law tail. During the QNM ringing stage, the QNM from the BH shows damped oscillations, characterized by a discrete set of complex frequencies. Here, the real part represents the frequency of the oscillation, while the imaginary part indicates the rate at which the oscillation fades.

Concerning the calculation methodologies for QNMs, we will consider the finite element method as well as the pseudospectral methods [21]. Together with the above, at the moment to consider quantum effect within a BH, the emission of Hawking radiation occurs [22]. Nevertheless, the presence of strong gravitational effects in its proximity acts as a barrier, implying that we observe is not a radiation spectrum, but rather a modified form known as a greybody spectrum. This spectrum is distinctly influenced by the curvature of spacetime, encapsulated by the greybody factor [23, 24].

Focused on the Anti-de Sitter/Conformal Field Theory (AdS/CFT) correspondence [20, 25, 26, 27], a large static BH asymptotically AdS (and its perturbation) corresponds to a thermal state in the CFT (perturbing this thermal state). Here, the QNFs in AdS BHs allow us to obtain a prediction for the thermalization timescale for a strongly coupled CFT [1]. In addition, from a thermodynamic point of view, the study of phase transitions (PTs) in strongly coupled field theories is one of the most relevant aspects of the gauge/gravity correspondence, where Hawking and Page, in their pioneering work [28], showed a PT between spherically symmetric AdS BHs and thermal AdS space-time, this represents confining and deconfining PTs in the dual quark-gluon plasma [26, 29]. Along with the above, PTs have been an object of high study when considering AdS-charged BHs, which show a resemblance to the Van der Waals fluid [30, 31].

On the other hand, by exploring gravity theories beyond General Relativity, Lü and Pope in [32] showed that via the action

	$\displaystyle S_{\tiny{\mbox{CG}}}[g_{\mu\nu},R_{\mu\nu\sigma\rho}]$	$\displaystyle=$	$\displaystyle\int{d}^{4}x\sqrt{-g}\mathcal{L}_{\tiny{\mbox{CG}}}$
		$\displaystyle=$	$\displaystyle\int{d}^{4}x\sqrt{-g}\left[\frac{R-2\Lambda}{2\kappa}\right.$
		$\displaystyle+$	$\displaystyle\left.\frac{1}{2\kappa}\left(-\frac{1}{2\Lambda}{R}^{2}+\frac{3}{2\Lambda}{R}_{\alpha\beta}{R}^{\alpha\beta}\right)\right],$

where $\kappa$ is a coupling constant and $\Lambda$ the cosmological constant, a four-dimensional renormalizable theory devoid of ghosts appears, being known as Critical Gravity (CG). For the theory (I), the massive scalar mode is eliminated, while the massive spin-2 particle becomes massless, and the on-shell energy of the remaining massless gravitons becomes zero. Looking at it from a thermodynamic perspective, and according to the authors, although the theory (I) supports the AdS metric as well as the well-known Schwarzschild AdS (Sch-AdS) BH, the adjustments integrated into the parameters of the action (I) yield null thermodynamic extensive quantities, being the price to pay to obtain a well-behaved theory [32]. The above result has generated the exploration of BHs, supported by CG (I), with non-vanishing thermodynamic properties. One of them is via with non-standard asymptotic behavior [33], through the inclusion of dilatonic fields and a non-minimally coupled scalar field, or asymptotically AdS configurations with planar base manifold [34, 35] with nonlinear electrodynamics (NLE). For instance, NLE can help circumvent the singularity of the field of a point particle in standard electrodynamics, attempt to describe quantum electrodynamics [36], and obtain stationary solutions [39, 38, 37]. This has given rise to well-known theories such as the Born and Infeld theory [40, 41, 42, 43, 44, 45], the Euler-Heisenberg models [46], ModMax electrodynamics [47, 48] and the construction of a NLE through a sum of infinite series of Maxwell invariant [49].

The search for new AdS configurations supported by CG with non-vanishing thermodynamic properties and motivated by gauge/gravity duality, has led to the consideration of NLE as a matter source. In the present work, we consider a nonlinear behavior of the electromagnetic field $A_{\mu}$ with field strength $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$, and the antisymmetric tensor $\mathcal{P}_{\mu\nu}$ (known as Plebánski tensor), which reads [50, 51]:

	$\displaystyle S_{\tiny{\mbox{NLE}}}[A_{\mu},P_{\mu\nu}]$	$\displaystyle=$	$\displaystyle\int{d}^{4}x\sqrt{-g}\mathcal{L}_{\tiny{\mbox{NLE}}}$
		$\displaystyle=$	$\displaystyle\int{d}^{4}x\sqrt{-g}\left(-\frac{1}{2}\mathcal{P}^{\mu\nu}F_{\mu\nu}+\mathcal{H}(\mathcal{P})\right).$

The introduction of $\mathcal{P}_{\mu\nu}$ arises from the need to establish a relationship between a standard electromagnetic theory with Maxwell’s theory of continuous media, together with a structure-function $\mathcal{H}(\mathcal{P})$ depending on the invariant formed with $\mathcal{P}_{\mu\nu}$ ($\mathcal{P}:=\frac{1}{4}\mathcal{P}_{\mu\nu}\mathcal{P}^{\mu\nu}$) where the linear Maxwell scenario is naturally recovered when $\mathcal{H}(\mathcal{P})\propto\mathcal{P}$. The action (I) has been studied extensively in the context of gravitational theories, regular non-rotating BHs, and charged rotating BHs configurations (see Refs. [52, 53, 54, 55, 56, 38, 39]).

Given the aforementioned details, in the present paper, our first objective is (i) to show that the planar extension of the work present in [34] can be extended to an arbitrary topology base manifold, allowing us to explore asymptotically AdS BHs, with one or even more locations of the event horizon. Additionally, via the Wald formalism as procedure [57, 58], we present that these charged BHs configurations enjoy non-zero thermodynamic properties, wherein the topology plays a crucial role. In fact, via Gibbs free energy and comparing with respect to the thermal AdS space-time, we will have first-order PTs as well as situations where these BHs are the preferred configuration. With this information, the second objective is (ii) to delve into the study of QNMs as well as the calculation of greybody factors, considering the spherical situation.

The novelty of the present work is based on two key dimensions. Firstly, by combining the gravity theory (I) with the matter source (I), it opens avenues to explore charged BHs configurations beyond the planar scenario, as well as their thermodynamic properties. Secondly, it introduces an opportunity to explore QNMs for the spherical case, where the influence of quadratic corrections stemming from CG (I) and NLE as a matter source becomes a significant factor.

The structure of this paper unfolds in the following manner: In Section II we present the charged BH solution, exploring the existence of horizons independent of the topology of the event horizon, while in Section III the thermodynamic properties are explored. In Section IV the QNMs and the greybody factor are computed and analyzed. Finally, Section V is devoted to our conclusions and discussions.

To perform a complete exploration and analysis for this charged BH configurations, the equation of motions that result from the variation of the action

	$\displaystyle S[g_{\mu\nu},R_{\mu\nu\sigma\rho},A_{\mu},P_{\mu\nu}]$	$\displaystyle=$	$\displaystyle S_{\tiny{\mbox{CG}}}[g_{\mu\nu},R_{\mu\nu\sigma\rho}]$
		$\displaystyle+$	$\displaystyle S_{\tiny{\mbox{NLE}}}[A_{\mu},\mathcal{P}_{\mu\nu}],$

are required, which are given by

	$\displaystyle\mathcal{E}^{\nu}_{F}$	$\displaystyle=$	$\displaystyle\nabla_{\mu}P^{\mu\nu}=0,$		(4a)
	$\displaystyle\mathcal{E}^{\mu\nu}_{P}$	$\displaystyle=$	$\displaystyle-F^{\mu\nu}+\left(\frac{\partial\mathcal{H}}{\partial\mathcal{P}}\right)\mathcal{P}^{\mu\nu}=0,$		(4b)
	$\displaystyle\mathcal{E}_{\mu\nu}$	$\displaystyle=$	$\displaystyle{\mathcal{G}_{\mu\nu}^{CG}}-\kappa T_{\mu\nu}^{NLE}=0,$		(4c)
and the tensors ${\mathcal{G}_{\mu\nu}^{CG}}$ and $T^{NLE}_{\mu\nu}$ are defined as follows:
	$\displaystyle{\mathcal{G}_{\mu\nu}^{CG}}$	$\displaystyle=$	$\displaystyle G_{\mu\nu}+\Lambda g_{\mu\nu}+\frac{3}{\Lambda}\Bigl{(}R_{\mu\rho}R_{\nu}^{\,\rho}-\frac{1}{4}R^{\rho\sigma}R_{\rho\sigma}g_{\mu\nu}\Bigr{)}$
		$\displaystyle-$	$\displaystyle\frac{1}{\Lambda}\Bigl{(}RR_{\mu\nu}-\frac{1}{4}R^{2}g_{\mu\nu}\Bigr{)}+\frac{3}{2\Lambda}\Bigl{(}\Box R_{\mu\nu}+\frac{1}{2}\Box Rg_{\mu\nu}$
		$\displaystyle-$	$\displaystyle 2\nabla_{\rho}\nabla_{\left(\mu\right.}R_{\left.\nu\right)}^{\,\rho}\Bigr{)}-\frac{1}{\Lambda}(g_{\mu\nu}\Box R-\nabla_{\mu}\nabla_{\nu}R),$
	$\displaystyle T_{\mu\nu}^{NLE}$	$\displaystyle=$	$\displaystyle\left(\frac{\partial\mathcal{H}}{\partial\mathcal{P}}\right)\mathcal{P}_{\mu\alpha}\mathcal{P}_{\nu}^{\,\alpha}-g_{\mu\nu}\left(2\mathcal{P}\left(\frac{\partial\mathcal{H}}{\partial\mathcal{P}}\right)-\mathcal{H}\right).$

Here, eq. (4a) stands as the nonlinear rendition of Maxwell’s equations, encapsulating the constitutive relations within (4b). Meanwhile, (4c) represents the Einstein equations. For the model (I)-(II), we consider the following four-dimensional metric Ansatz:

	$\displaystyle ds^{2}$	$\displaystyle=$	$\displaystyle-\left(\epsilon+\frac{r^{2}}{l^{2}}f(r)\right)dt^{2}+\left(\epsilon+\frac{r^{2}}{l^{2}}f(r)\right)^{-1}{dr^{2}}$		(5)
		$\displaystyle+$	$\displaystyle r^{2}d\Omega^{2}_{2,\epsilon},$

where $t\in\,(-\infty,+\infty),r>0$, $l$ is the AdS radius, and the metric function must satisfy the asymptotic condition

$${\lim_{r\rightarrow+\infty}f(r)=1,}$$

while that $d\Omega^{2}_{2,\epsilon}$ represents the line element given by

$\displaystyle d\Omega^{2}_{2,\epsilon}=\left\{\begin{array}[]{ll}d\theta^{2}+\sin^{2}(\theta)d\rho^{2},&\mbox{for $\epsilon=1,$}\\ d\theta^{2}+\sinh^{2}(\theta)d\rho^{2},&\mbox{for $\epsilon=-1,$}\\ d\theta^{2}+d\rho^{2},&\mbox{for $\epsilon=0,$}\end{array}\right.$

(9)

indexed through the constant $\epsilon$. For our analysis, we consider purely electrical configurations; this is $\mathcal{P}_{\mu\nu}=2\delta^{t}_{[\mu}\delta_{\nu]}^{r}G(r)$ where, if we replace it in the nonlinear Maxwell equation (4a), we obtain

$$\mathcal{P}_{\mu\nu}=2\delta_{[\mu}^{t}\delta_{\nu]}^{r}\frac{M}{r^{2}},$$

where $M$ is an integration constant. Therefore, the electric invariant $\mathcal{P}$ is negative definite since we only consider purely electrical configurations, which reads

$$\mathcal{P}=\frac{1}{2}{\mathcal{P}_{rt}\mathcal{P}^{rt}}=-\frac{M^{2}}{2r^{4}}<0.$$

(10)

With the above information, we note that the difference between the temporal and radial diagonal components of Einstein’s equations (4c) (this is $\mathcal{E}_{t}^{t}-\mathcal{E}_{r}^{r}=0$) is proportional to a fourth-order Cauchy-Euler ordinary differential equation, where the metric function $f$ enjoys the structure $f(r)=M_{0}+\sum_{i=1}^{3}M_{i}\left(\frac{l}{r}\right)^{i}$, with $M_{i}$’s integration constants, while that via the combination of $\mathcal{E}^{t}_{t}-\mathcal{E}_{\theta}^{\theta}=0$ (or $\mathcal{E}^{t}_{t}-\mathcal{E}_{\rho}^{\rho}=0$), and expressing the radial coordinate $r$ in terms of the electric invariant $\mathcal{P}$ from eq. (10), one can determine $\partial\mathcal{H}/\partial\mathcal{P}$, which can be later integrated to obtain $\mathcal{H}(\mathcal{P})$, which reads [34]

	$\displaystyle\mathcal{H}(\mathcal{P})$	$\displaystyle=$	$\displaystyle\frac{(\alpha_{2}^{2}-3\alpha_{1}\alpha_{3})l^{2}\mathcal{P}}{3\kappa}{-\frac{2\alpha_{1}(-2\mathcal{P})^{1/4}}{l\kappa}}$		(11)
		$\displaystyle+$	$\displaystyle{\frac{\alpha_{2}\sqrt{-2\mathcal{P}}}{\kappa}},$

where $\alpha_{1}$, $\alpha_{2}$ and $\alpha_{3}$ are coupling constants, together with a metric function

	$\displaystyle f(r)$	$\displaystyle=$	$\displaystyle 1-\alpha_{1}\sqrt{M}\left(\frac{l}{r}\right)+\alpha_{2}M\left(\frac{l}{r}\right)^{2}$		(12)
		$\displaystyle-$	$\displaystyle\sqrt{M}\left(\alpha_{3}M+\frac{2\epsilon\alpha_{2}}{3\alpha_{1}}\right)\left(\frac{l}{r}\right)^{3},$

after a redefinition of the integration constants $M_{i}$’s. Finally, to satisfy the equations of motion (4a)-(4c), the cosmological constant must be fixed as

$\displaystyle\Lambda=-\frac{3}{l^{2}},$

(13)

where we note that for the planar situation ($\epsilon=0$), the configuration obtained in [34] is naturally recovered. It is important to note that the uncharged Sch-AdS case can be obtained from the structure-function (11) as well as the metric function (12) by first fixing $\alpha_{2}=0$ and then $\alpha_{1}=0$. This implies that even the linear Maxwell case is not allowed, which emphasizes the importance of including this matter source.

On the other hand, regarding the line element (5) and the metric function (12), in order to obtain a BH configuration, some conditions need to be established. First that all, the scalar curvature $R$ reads

	$\displaystyle R$	$\displaystyle=$	$\displaystyle-\frac{12f}{l^{2}}-\frac{8rf^{\prime}}{l^{2}}-\frac{r^{2}f^{\prime\prime}}{l^{2}},$		(14)
		$\displaystyle=$	$\displaystyle-\frac{12}{l^{2}}+\frac{6\alpha_{1}\sqrt{M}}{rl}-\frac{2\alpha_{2}M}{r^{2}},$

where for our notations (${}^{\prime})$ denotes the derivative with respect to coordinate $r$, and showing us a curvature singularity located at $r_{s}=0$. With respect to the existence of event horizon, this is $r_{h}>0$ such that

$$\Big{(}\epsilon+\frac{r^{2}}{l^{2}}f(r)\Big{)}\Big{|}_{r=r_{h}}=0,\mbox{ and }\Big{(}\epsilon+\frac{r^{2}}{l^{2}}f(r)\Big{)}^{\prime}\Big{|}_{r=r_{h}}>0,$$

we note that in the limit, as $r$ approaches positive infinity, we have that $\epsilon+\frac{r^{2}}{l^{2}}f(r)\simeq\frac{r^{2}}{l^{2}}$. While that when $r\rightarrow 0^{+}$, $\epsilon+\frac{r^{2}}{l^{2}}f(r)\simeq-\sqrt{M}\left(\alpha_{3}M+\frac{2\epsilon\alpha_{2}}{3\alpha_{1}}\right)\frac{l}{r}$, allowing us to split the analysis in two cases:

	$\displaystyle\mbox{{\bf{Case 1:}}}{\mbox{ If }}\alpha_{3}M+\frac{2\epsilon\alpha_{2}}{3\alpha_{1}}>0$	$\displaystyle:$	$\displaystyle\epsilon+\frac{r^{2}}{l^{2}}f(r)\mbox{ starts}$		(15)
			increasing from
			$-\infty$ when
			$r$ increases.

	$\displaystyle\mbox{{\bf{Case 2:}}}{\mbox{ If }}\alpha_{3}M+\frac{2\epsilon\alpha_{2}}{3\alpha_{1}}<0$	$\displaystyle:$	$\displaystyle\epsilon+\frac{r^{2}}{l^{2}}f(r)\mbox{ starts}$
			decreasing from
			$+\infty$ when
			$r$ increases,

as is shown in Figure 1.

For these situations, we can conduct a study of the extreme values, denoted as $r_{\text{ext}}>0$, for $\epsilon+\frac{r^{2}}{l^{2}}f(r)$, via the first derivative with respect to the radial coordinate $r$. For this analysis, we observe that these values satisfy

	$\displaystyle\Big{(}\epsilon+\frac{r^{2}}{l^{2}}f(r)\Big{)}^{\prime}{\Big{|}}_{r=r_{\text{ext}}}$		(16)
	$\displaystyle=\Big{[}2\Big{(}\frac{r}{l}\Big{)}-\alpha_{1}\sqrt{M}+\frac{l^{2}}{r^{2}}\sqrt{M}\Big{(}\alpha_{3}M+\frac{2\epsilon\alpha_{2}}{3\alpha_{1}}\Big{)}\Big{]}{\Big{|}}_{r=r_{\text{ext}}}$
	$\displaystyle=0.$

Furthermore, to determine whether these extreme values correspond to a maximum or minimum, we examine the sign of the expression

$\displaystyle\Big{(}\epsilon+\frac{r^{2}}{l^{2}}f(r)\Big{)}^{\prime\prime}=2\left[1-\frac{l^{3}}{r^{3}}\sqrt{M}\left(\alpha_{3}M+\frac{2\epsilon\alpha_{2}}{3\alpha_{1}}\right)\right],$

(17)

through the second derivative, where when evaluated at $r=r_{\text{ext}}$, the above equation can be written as

$\displaystyle 2\left(3-\alpha_{1}\sqrt{M}\frac{l}{r_{\text{ext}}}\right).$

(18)

For Case 1, we note that when eq. (17) is zero, an inflection point arises (denoted as $r_{\text{inf}}>0$). For the special conditions

	$\displaystyle\Big{(}\epsilon+\frac{r^{2}}{l^{2}}f(r)\Big{)}^{\prime}{\Big{|}}_{r=r_{\text{inf}}}=\left(\frac{9\sqrt{M}(3\alpha_{3}M\alpha_{1}+2\alpha_{2}\epsilon)}{\alpha_{1}}\right)^{\frac{1}{3}}$
	$\displaystyle-\alpha_{1}\sqrt{M}<0,$		(19)
	$\displaystyle\left(3-\alpha_{1}\sqrt{M}\frac{l}{r_{\text{ext}^{(1)}}}\right)<0,\qquad\left(3-\alpha_{1}\sqrt{M}\frac{l}{r_{\text{ext}^{(2)}}}\right)>0,$

the existence of two extreme values is ensured, as depicted in Figure 2.

Under this scenario, $r_{\text{ext}^{(1)}}$ corresponds to a local maximum, while $r_{\text{ext}^{(2)}}$ represents a local minimum. It is important to note that $r_{\text{ext}^{(1)}}<r_{\text{inf}}<r_{\text{ext}^{(2)}}$, and the determination of whether there are one or three horizons is contingent upon the sign of

$\displaystyle\displaystyle{\Big{(}\epsilon+\frac{r^{2}}{l^{2}}f(r)\Big{)}{\Big{|}}_{r=r_{\text{ext}^{(i)}}},\mbox{ with }i=\{1,2\}.}$

(20)

All these situations can be represented graphically in Figure 3.

In contrast, for Case 2, it can be observed from eq. (17) that there are no inflection points for $r>0$. Additionally, examining its first derivative from eq. (16), we note the existence of only one extreme value, as depicted in Figure 4. This implies that $\epsilon+{r^{2}f(r)}/{l^{2}}$ behaves without any change in concavity and the existence of horizons depend on the sign of (20) at $r=r_{\text{ext}}$, as it shown in Figure 5.

Building upon the analysis and framework established for this new four-dimensional charged BH configuration, and considering the interplay between the integration constant $M$ and the coupling constants $\alpha_{i}$’s, the subsequent section will delve the thermodynamic quantities.

In this section, we derive the thermodynamic quantities of the charged solution (5)-(13), and to simplify our calculations, we consider

$$r_{h}=\zeta\sqrt{M}l,$$

(21)

with $\zeta$ is a positive constant. Among these quantities, the first one of interest is the Hawking Temperature (denoted as $T$), which can be calculated as follows:

	$\displaystyle T$	$\displaystyle=$	$\displaystyle\frac{1}{4\pi}\Big{(}\epsilon+\frac{r^{2}}{l^{2}}f(r)\Big{)}^{\prime}\Big{|}_{r=r_{h}}$		(22)
		$\displaystyle=$	$\displaystyle\frac{3r_{h}}{4\pi l^{2}}-\frac{\alpha_{1}\sqrt{M}}{2\pi l}+\frac{\alpha_{2}M}{4\pi r_{h}}+\frac{\epsilon}{4\pi r_{h}}$
		$\displaystyle=$	$\displaystyle\frac{r_{h}\Psi_{1}}{4\pi\zeta^{2}l^{2}}+\frac{\epsilon}{4\pi r_{h}},$

where

$\displaystyle\Psi_{1}=3\zeta^{2}-2\alpha_{1}\zeta+\alpha_{2}.$

(23)

Here, $r_{h}$ represents the position of the event horizon and for some special condition, we obtain $T>0$. Indeed, when $\epsilon=0$ (the planar case), the Hawking temperature is positive when $\Psi_{1}>0$, where occurs when $\alpha_{1}>0$ and $3\alpha_{2}-\alpha_{1}^{2}>0$. On the other hand, for $\epsilon\neq 0$, via the derivative of $T$ with respect to $r_{h}$, we have that:

$$\frac{dT}{dr_{h}}=\frac{\Psi_{1}}{4\pi\zeta^{2}l^{2}}-\frac{\epsilon}{4\pi r_{h}^{2}},\quad\frac{d^{2}T}{dr_{h}^{2}}=\frac{\epsilon}{2\pi r_{h}^{3}},$$

where $T$ reaches a local extremum at

$\displaystyle r_{h}^{*}=\zeta l\sqrt{\frac{\epsilon}{\Psi_{1}}},$

(24)

and the concavity depends on the sign of $\epsilon$. To obtain a real $r_{h}^{*}$, we can to separate the cases for both spherical ($\epsilon=1$) and hyperbolic situation ($\epsilon=-1$). For $\epsilon=1$ we have that

$$\Psi_{1}>0,\quad\frac{d^{2}T}{dr_{h}^{2}}>0,\quad T^{*}=T\big{|}_{r=r_{h}^{*}}=\frac{\sqrt{\epsilon\Psi_{1}}}{2\pi l\zeta}>0,$$

ensuring positivity for $T$. For $\epsilon=-1$, $r_{h}^{*}$ is real when $\Psi_{1}<0$ and a local minimum is obtained. Nevertheless, it is important to note that here $T<0$. For the remaining situation, $T$ is an increasing (for $\Psi_{1}>0$ and $\epsilon=-1$) or decreasing function (for $\Psi_{1}<0$ and $\epsilon=1$) on $r_{h}$, as shown in Figure 6.

Additionally, the remaining thermodynamic parameters can be determined through the Wald formalism [57, 58]. The first step in this formalism involves considering the variation of the action (II):

$\displaystyle\delta S=\mathcal{E}_{\mu\nu}\delta g^{\mu\nu}+\mathcal{E}^{\nu}_{F}\delta A_{\nu}+\mathcal{E}^{\mu\nu}_{P}\delta\mathcal{P}_{\mu\nu}+\partial_{\mu}{\cal{J}}^{\mu}.$

(25)

For our notations, $\mathcal{E}^{\nu}_{F}$, $\mathcal{E}^{\mu\nu}_{P}$ and $\mathcal{E}_{\mu\nu}$ represent the equations of motions (4a), (4b) and (4c) respectively, while that ${\cal{J}}^{\mu}$ is a surface term, given by

	$\displaystyle{\cal{J}}^{\mu}$	$\displaystyle=$	$\displaystyle\sqrt{-g}\Bigg{[}2\left(P^{\mu(\alpha\beta)\gamma}\nabla_{\gamma}\delta g_{\alpha\beta}-\delta g_{\alpha\beta}\nabla_{\gamma}P^{\mu(\alpha\beta)\gamma}\right)$		(26)
		$\displaystyle+$	$\displaystyle\frac{\delta\mathcal{L}_{\tiny{\mbox{NLE}}}}{\delta(\partial_{\mu}A_{\nu})}\delta A_{\nu}\Bigg{]},$

where

	$\displaystyle P^{\alpha\beta\gamma\delta}$	$\displaystyle=$	$\displaystyle\frac{\delta\mathcal{L}_{\tiny{\mbox{CG}}}}{\delta R_{\alpha\beta\gamma\delta}}=\frac{1}{4\kappa}\,\Big{(}g^{\alpha\gamma}g^{\beta\delta}-g^{\alpha\delta}g^{\beta\gamma}\Big{)}$
		$\displaystyle-$	$\displaystyle\frac{1}{4\kappa\Lambda}\,R\left(g^{\alpha\gamma}g^{\beta\delta}-g^{\alpha\delta}g^{\beta\gamma}\right)$
		$\displaystyle+$	$\displaystyle\frac{3}{8\kappa\Lambda}\,\big{(}g^{\beta\delta}R^{\alpha\gamma}-g^{\beta\gamma}R^{\alpha\delta}-g^{\alpha\delta}R^{\beta\gamma}$
		$\displaystyle+$	$\displaystyle g^{\alpha\gamma}R^{\beta\delta}\big{)},$
	$\displaystyle\frac{\delta\cal{L}_{\tiny{\mbox{NLE}}}}{\delta(\partial_{\mu}A_{\nu})}$	$\displaystyle=$	$\displaystyle-\mathcal{P}^{\mu\nu}.$		(27)

Through the surface term (26), we define a $1$-form ${\cal{J}}_{(1)}={\cal{J}}_{\mu}dx^{\mu}$ and its Hodge dual ${\Theta}_{(3)}=-(*{\cal{J}}_{(1)})$. After employing the equations of motion ($\mathcal{E}_{\mu\nu}$, $\mathcal{E}^{\nu}_{F}$, and $\mathcal{E}^{\mu\nu}_{P}$), we obtain

$${\cal{J}}_{(3)}={\Theta}_{(3)}-i_{\xi}*(\mathcal{L}_{\tiny{\mbox{CG}}}+\mathcal{L}_{\tiny{\mbox{NLE}}})=d(*{\cal{J}}_{(2)}).$$

(28)

Here, $i_{\xi}$ denotes a contraction of the vector field $\xi^{\mu}$ on the first index of $*(\mathcal{L}_{\tiny{\mbox{CG}}}+\mathcal{L}_{\tiny{\mbox{NLE}}})$. The equation (28) allows to define a $2-$form ${Q}_{(2)}=*{\cal{J}}_{(2)}$ such that ${\cal{J}}_{(3)}=dQ_{(4)}$, given by

	$\displaystyle Q_{(2)}$	$\displaystyle:=$	$\displaystyle Q_{\alpha_{1}\alpha_{2}}=\epsilon_{\alpha_{1}\alpha_{2}\mu\nu}Q^{\mu\nu}$
		$\displaystyle=$	$\displaystyle\epsilon_{\alpha_{1}\alpha_{2}\mu\nu}\Big{[}2P^{\mu\nu\rho\sigma}\nabla_{\rho}\xi_{\sigma}-4\xi_{\sigma}\nabla_{\rho}P^{\mu\nu\rho\sigma}$
		$\displaystyle-$	$\displaystyle\frac{\delta\cal{L}}{\delta(\partial_{\mu}A_{\nu})}\xi^{\sigma}A_{\sigma}\Big{]},$

with $\epsilon_{\alpha_{1}\alpha_{2}\alpha_{3}\alpha_{4}}$ representing the Levi-Civita tensor. Subsequently, we consider the vector field $\xi^{\mu}$ as a time-translation vector, representing a Killing vector that becomes null at the event horizon location (denoted as $r_{h}$ as before).

Taking all these elements into account, we can now express the variation of the Hamiltonian $\delta\mathcal{H}$ in the following form:

	$\displaystyle\delta\mathcal{H}$	$\displaystyle=$	$\displaystyle\delta\int_{\mathcal{C}}{\cal{J}}_{(3)}-\int_{\mathcal{C}}d\left(i_{\xi}\Theta_{(3)}\right)$		(30)
		$\displaystyle=$	$\displaystyle\int_{\Sigma^{(3)}}\left(\delta{Q}_{(2)}-i_{\xi}{\Theta}_{(3)}\right),$

with

	$\displaystyle\delta{Q}_{(2)}-i_{\xi}{\Theta}_{(3)}=\left[-\frac{r^{2}l^{2}f}{6\kappa}\delta(f^{\prime\prime\prime})-\frac{rl^{2}(2f-f^{\prime}r)}{12\kappa}\right.$
	$\displaystyle\times\delta(f^{\prime\prime})-\frac{l^{2}(f^{\prime}r-2f)}{6\kappa}\delta(f^{\prime})$
	$\displaystyle\left.+\frac{(r^{3}l^{2}f^{\prime\prime\prime}+4l^{2}f-2l^{2}rf^{\prime}+12r^{2})}{12r\kappa}\delta f\right]\Omega_{2,\epsilon}$
	$\displaystyle-\Omega_{2,\epsilon}A\,\delta\left(r^{2}A^{\prime}\left(\frac{\partial\mathcal{H}}{\partial\mathcal{P}}\right)^{-1}\right)^{\prime},$		(31)

where $\Omega_{2,\epsilon}$ is the finite volume of the compact base manifold, $\mathcal{C}$ represents a Cauchy Surface with boundary $\Sigma^{(2)}$. According to the Wald formalism, the variation (30) can be decomposed into two main components: one located at infinity $\mathcal{H}_{\infty}$ and the other at the horizon $\mathcal{H}_{+}$, and the first law of black hole thermodynamics is a consequence of

$\displaystyle\delta\mathcal{H}_{\infty}=\delta\mathcal{M}=T\delta\mathcal{S}+\Phi_{e}\delta\mathcal{Q}_{e}=\delta\mathcal{H}_{+}.$

(32)

Here, $\mathcal{M}$, $\mathcal{S}$, and $\mathcal{Q}{e}$ represent the mass, entropy, and electric charge. Meanwhile, $T$ and $\Phi_{e}$ take the role of temperature and electric potential, respectively. For our situation, given the solution (5)-(13), we have that at the infinity

$\displaystyle\delta{\cal{H}}_{\infty}$

$\displaystyle=$

$\displaystyle\left(\frac{\alpha_{1}\alpha_{2}r_{h}^{2}\Omega_{2,\epsilon}}{3\kappa l^{2}\zeta^{2}}\right)\delta r_{h},$

(33)

while that at the horizon

	$\displaystyle\delta{\cal{H}}_{+}$	$\displaystyle=$	$\displaystyle T\Omega_{2,\epsilon}\left(\frac{4\pi r_{h}\alpha_{1}}{\zeta\kappa}-\frac{8\pi r_{h}\alpha_{2}}{3\zeta^{2}\kappa}+\epsilon\Gamma\right)\delta r_{h}$		(34)
		$\displaystyle+$	$\displaystyle\Phi_{e}\Omega_{2,\epsilon}\left(\frac{2r_{h}}{l^{2}\zeta^{2}}\right)\delta r_{h},$

where, as before, $T$ is the Hawking temperature (22), the electric potential $\Phi_{e}$ is defined as

	$\displaystyle\Phi_{e}$	$\displaystyle=$	$\displaystyle-A_{t}(r)\big{|}_{r=r_{h}}$		(35)
		$\displaystyle=$	$\displaystyle\left(-\frac{\alpha_{1}r^{2}}{2\sqrt{M}l\kappa}+\frac{\alpha_{2}r}{\kappa}-\frac{Ml^{2}(3\alpha_{3}\alpha_{1}-\alpha_{2}^{2})}{3\kappa r}\right)\Big{|}_{r=r_{h}},$
		$\displaystyle=$	$\displaystyle\frac{r_{h}}{\kappa}\left(\frac{\alpha_{1}\alpha_{2}}{6\zeta}-\frac{\Psi_{1}\Psi_{2}}{6\zeta^{2}}\right)-\frac{\epsilon l^{2}\Psi_{2}}{3r_{h}\kappa},$

where

$\displaystyle\Psi_{2}$

$\displaystyle=$

$\displaystyle 3\alpha_{1}\zeta-2\alpha_{2},$

(36)

while that $\Gamma$ represents the contribution arising from the event horizon’s topology, its expression can be written as follows:

$\displaystyle\Gamma=-\frac{8\pi l^{2}}{3\kappa r_{h}}+\frac{4\zeta l^{2}(2\epsilon\zeta l^{2}+6r_{h}^{2}\zeta-r_{h}^{2}\alpha_{1})\pi}{3r_{h}\kappa(3\zeta^{2}r_{h}^{2}+\alpha_{2}r_{h}^{2}-2\alpha_{1}\zeta r_{h}^{2}+\epsilon\zeta^{2}l^{2})}.$

Finally, identifying (32), (33) and (34), the mass $\mathcal{M}$, entropy $\mathcal{S}$ and electric charge $\mathcal{Q}_{e}$ read

	$\displaystyle\mathcal{M}$	$\displaystyle=$	$\displaystyle\frac{\alpha_{1}\alpha_{2}r_{h}^{3}\Omega_{2,\epsilon}}{9\kappa l^{2}\zeta^{2}},$		(37)
	$\displaystyle\mathcal{S}$	$\displaystyle=$	$\displaystyle\frac{2\pi\Psi_{2}\Omega_{2,\epsilon}}{3\kappa\zeta^{2}}\left[{r_{h}^{2}}+\frac{\epsilon\zeta^{2}l^{2}}{\Psi_{1}}\ln\left(\Psi_{1}r_{h}^{2}+\epsilon\zeta^{2}l^{2}\right)\right],$		(38)
	$\displaystyle\mathcal{Q}_{e}$	$\displaystyle=$	$\displaystyle\frac{r^{2}_{h}\Omega_{2,\epsilon}}{l^{2}\zeta^{2}},$		(39)

where the $\Psi_{i}$’s are given in (23) and (36) respectively.

Let us notice that, analogous to the Hawking temperature $T$, for some suitable election of the constants $\alpha_{i}$’s and $\zeta$, the thermodynamic parameters (35) and (37)-(39) are positive. Firstly, the charge $\mathcal{Q}_{e}$ (39) is always a positive quantity, due to that the principal components are the event horizon $r_{h}>0$, the AdS radius $l>0$, the finite volume element of the compact base manifold $\Omega_{2,\epsilon}>0$, and from eq. (40) $\zeta>0$, while that $\mathcal{M}>0$ only if $\alpha_{1}\alpha_{2}>0$. For the case of the entropy $\mathcal{S}$, the presence of the topology of the event horizon (represented by $\epsilon$) yields a logarithmic behavior, where for $\Psi_{1}>0$, $\Psi_{2}>0$ as well as $\Psi_{1}r_{h}^{2}+\epsilon\zeta^{2}l^{2}>1$ we ensured that $\mathcal{S}>0$. Additionally, for $\Phi_{e}>0$, we consider:

$\displaystyle\frac{d\Phi_{e}}{dr_{h}}=\frac{1}{\kappa}\left(\frac{\alpha_{1}\alpha_{2}}{6\zeta}-\frac{\Psi_{1}\Psi_{2}}{6\zeta^{2}}\right)+\frac{\epsilon l^{2}\Psi_{2}}{3r^{2}_{h}\kappa},\,\,\frac{d^{2}\Phi_{e}}{dr_{h}^{2}}=-\frac{\epsilon l^{2}\Psi_{2}}{3\kappa r_{h}^{3}},$

and the concavity depends on the sign of $-\epsilon\Psi_{2}$. Here, the electric potential $\Phi_{e}$ reaches an extremum at $r_{h}^{*}=\sqrt{-(2\epsilon\zeta^{2}l^{2}\Psi_{2})/(\alpha_{2}\alpha_{1}\zeta-\Psi_{1}\Psi_{2})},$ where when $\epsilon\Psi_{2}<0$ and $\alpha_{2}\alpha_{1}\zeta-\Psi_{1}\Psi_{2}>0$ we have that

$$\frac{d^{2}\Phi_{e}}{dr_{h}^{2}}>0,\quad\Phi_{e}\big{|}_{r=r_{h}^{*}}=\sqrt{-\frac{2\epsilon l^{2}\Psi_{2}(\alpha_{1}\alpha_{2}\zeta-\Psi_{1}\Psi_{2})}{9\kappa^{2}\zeta^{2}}}>0,$$

ensuring that $\Phi_{e}>0$. On the other hand, when $\epsilon\Psi_{2}>0$ and $\alpha_{2}\alpha_{1}\zeta-\Psi_{1}\Psi_{2}>0$, the electrical potential behaves as an increasing function, as shown in Figure 7.

On the other hand, we can examine the system’s response to small perturbations around equilibrium by considering these thermodynamic quantities. To achieve this, the extensive thermodynamic quantities $\mathcal{M}$ and $\mathcal{Q}_{e}$ can be expressed as functions of the intensive ones $T$ and $\Phi_{e}$. From eqs. (22)-(23) and (35)-(36), we note that the location of the event horizon $r_{h}$ can be expressed in the following way

$\displaystyle r_{h}=F(T,\Phi_{e})=\frac{2\zeta^{2}(4\pi l^{2}\Psi_{2}T+3\kappa\Phi_{e})}{(\alpha_{2}\alpha_{1}\zeta+\Psi_{2}\Psi_{1})},$

(40)

which is positive for $\Psi_{2}>0$,$\Psi_{1}>0$ and $\alpha_{1}\alpha_{2}>0$. Now, the mass $\mathcal{M}$ as well as the electric charge $\mathcal{Q}_{e}$ take the form

	$\displaystyle\mathcal{M}(T,\Phi_{e})$	$\displaystyle=$	$\displaystyle\frac{\alpha_{1}\alpha_{2}\Omega_{2,\epsilon}F^{3}(T,\Phi_{e})}{9\kappa\zeta^{3}l^{2}},$
	$\displaystyle\mathcal{Q}_{e}(T,\Phi_{e})$	$\displaystyle=$	$\displaystyle\frac{\Omega_{2,\epsilon}F^{2}(T,\Phi_{e})}{l^{2}\zeta^{2}},$

and with this, we are in a position to analyze the local thermodynamic (in)stability under thermal and electrical fluctuations, represented via the specific heat $C_{\Phi_{e}}$ and the electric permittivity $\epsilon_{T}$, given by

	$\displaystyle C_{\Phi_{e}}$	$\displaystyle=$	$\displaystyle\left(\frac{\partial\mathcal{M}}{\partial T}\right)_{\Phi_{e}}=\frac{\alpha_{1}\alpha_{2}\Omega_{2,\epsilon}F^{2}}{3\kappa\zeta^{3}l^{2}}\left(\frac{\partial F}{\partial{T}}\right)$		(41)
		$\displaystyle=$	$\displaystyle\frac{8\pi\alpha_{1}\alpha_{2}\Psi_{2}\Omega_{2,\epsilon}F^{2}}{3\zeta\kappa(\alpha_{2}\alpha_{1}\zeta+\Psi_{2}\Psi_{1})},$
	$\displaystyle\epsilon_{T}$	$\displaystyle=$	$\displaystyle\left(\frac{\partial\mathcal{Q}_{e}}{\partial\Phi_{e}}\right)_{T}=\frac{2\Omega_{2,\epsilon}F}{l^{2}\zeta^{2}}\left(\frac{\partial F}{\partial{\Phi_{e}}}\right)$
		$\displaystyle=$	$\displaystyle\frac{12\kappa\Omega_{2,\epsilon}F}{l^{2}(\alpha_{2}\alpha_{1}\zeta+\Psi_{2}\Psi_{1})},$

where the sub-index stands for constant electric potential $\Phi_{e}$ and constant temperature $T$ respectively. It is important to note that the non-negativity of equations given in (41) ensure local stability under thermal fluctuations ($C_{\Phi_{e}}\geq 0$) and electrical fluctuations ($\epsilon_{T}\geq 0$). For example, for $\alpha_{1}\alpha_{2}>0$, $\alpha_{2}\alpha_{1}\zeta+\Psi_{2}\Psi_{1}>0$, $\Psi_{1}>0$, and $\Psi_{2}\geq 0$ we have that $C_{\Phi_{e}}\geq 0$ and $\epsilon_{T}\geq 0$, as shown in Figure 8. Together with the above, with eq. (40) and in the grand canonical ensemble where the temperature and electric potential are fixed, we can compute the Gibbs free energy ${G}(T,\Phi_{e})=\mathcal{M}-T\mathcal{S}-\Phi_{e}\mathcal{Q}_{e}$:

	$\displaystyle G$	$\displaystyle=$	$\displaystyle\frac{\alpha_{1}\alpha_{2}\Omega_{2,\epsilon}F^{3}}{9\kappa\zeta^{3}l^{2}}-\frac{2\pi\Psi_{2}\Omega_{2,\epsilon}TF^{2}}{3\kappa\zeta^{2}}-\frac{\Phi_{e}F^{2}\Omega_{2,\epsilon}}{l^{2}\zeta^{2}}$		(42)
		$\displaystyle-$	$\displaystyle\frac{2\pi\Omega_{2,\epsilon}\epsilon l^{2}\Psi_{2}T}{3\kappa\Psi_{1}}\ln\big{(}\Psi_{1}F^{2}+\epsilon\zeta^{2}l^{2}\big{)}.$

Although the previous equation is a priori unwieldy, via some simulations, we can obtain interesting cases when the coupling constants $\alpha_{i}$’s are fixed. For the spherical situation, we note a first-order PT for some special relation between $T$ and $\mathcal{Q}_{e}$ (see Fig. 9, up panel). For the hyperbolic case, we note that ${G}<0$, and the charged BH has lower free energy than the thermal AdS space-time, being the preferred state (see Fig. 9, down panel).

The evolution of small perturbations in BHs generally goes through three stages: the initial outburst, the QNM ringing, and the final power law tail. Here, the QNM ringing stage corresponds to the BH’s QNM, which typically exhibits damped oscillations in the form of a discrete set of complex frequencies, where the real part represents the frequency of the oscillation and the imaginary part represents the rate at which the oscillation decays. As was shown previously in the introduction, the QNMs are characterized by being independent of the initial perturbation and only dependent on the BH parameters, making them an important tool for studying BHs and gravity theories.

In the following lines, we will divide the content into four parts. The initial three segments elucidate various calculation methodologies, the first one delineates the numerical calculation of the dynamic evolution of a scalar field through the finite element method, then in the second part we elaborate on computing the QNMs across the complex plane using pseudospectral methods. For the third part, we introduce the concept of the greybody factors, and finally, we dedicate it to the presentation and analysis of the calculation results.