# Magnetically charged regular black holes in $f(R,T)$ gravity coupled to
nonlinear electrodynamics

Takol Tangphati<EMAIL_ADDRESS>School of Science, Walailak University,
Thasala, Nakhon Si Thammarat, 80160, Thailand.  Menglong Youk
<EMAIL_ADDRESS>Strong Gravity Group, Department of Physics, Faculty of
Science, Silpakorn University, Nakhon Pathom 73000, Thailand  Supakchai
Ponglertsakul<EMAIL_ADDRESS>Strong Gravity Group, Department of
Physics, Faculty of Science, Silpakorn University, Nakhon Pathom 73000,
Thailand

###### Abstract

We construct asymptotically flat, static spherically symmetric black holes
with regular centre in $f(R,T)$ gravity coupled to nonlinear electrodynamics
Lagrangian. We obtain generalized metric function of the Bardeen and Hayward
black holes. The null, weak and strong energy conditions of these solutions
are discussed. All the energy conditions hold outside the black hole’s outer
event horizon by appropriated choices of parameters. Quasinormal mode of
massive scalar perturbation is also investigated. Quasinormal frequencies are
computed via the sixth order Wentzel-Kramers-Brillouin (WKB) with Padé
approximation. All the imaginary parts of the frequencies are found to be
negative. Finally, we provide an analysis in the eikonal limit.

##  I Introduction

The most well-known gravitational theory describing the relation between
spacetime and matter is Einstein’s general relativity (GR). For over a decade,
this theory has been well-tested by the observations and experiments in the
weak field limit like our solar system, and the highly dense binary systems
[1, 2]. However, there are numerous open questions that GR fails to provide
answers, for instance, an accelerated expansion of the universe [3, 4], and
galaxy rotation curve [5]. Rather than applying auxiliary fields to the
theory, one could construct the modification of the GR as the extension based
on the original Einstein’s general relativity.

One of the modifications of GR is the $f(R)$ gravity where the Ricci scalar
$R$ in the Einstein-Hilbert action is replaced with an arbitrary function of
$R$ [6, 7]. This modification can describe the accelerated expansion of the
Universe without relying on the exotic matters [8]. Moreover, the
generalization of the $f(R)$ gravity theory leads to extra degrees of freedom
related to curvature invariants and scalar fields, which are called Extended
Theories of Gravity (ETG) [9]. These additional degrees of freedom play a
major role as effective fluids unlike the fluids of ordinary matter which is
adopted as sources of the field equations. One class of the ETG is the $f(T)$
gravity theory where the extension of the torsional gravity with an arbitrary
function of the torsion scalar $f(T)$ plays a major role in explaining the
cosmological and astrophysics problems [10]. Additionally, another class of
the extension of the Einstein’s gravity theory is the $f(Q)$ gravity theory
constructed from the symmetric teleparallel gravity which is based on the non-
metricity scalar $Q$. The modification of this theory represents the stable
dark energy causing the accelerated universe in which the matter perturbation
remains constant [11, 12].

In addition, the $f(R,T)$ gravity theory is designed to add the matter
components into the gravitational action by applying the arbitrary function of
the Ricci scalar $R$ along with the trace of the energy momentum tensor $T$.
This is proposed in [13] where the modified field equation is derived and
cosmological solution is analysed by introducing a self-interacting scalar
field. Numerous works of the $f(R,T)$ gravity theory have been investigated.
The cosmological solutions based on a homogeneous and isotropic spacetime
through a phase-space analysis are done in [14]. In addition, several
cosmological solutions from the $f(R,T)$ theory have been exclusively explored
in refs [15, 16, 17, 18, 19]. The violation of the energy conditions is
investigated in [20, 21]. Moreover, thermodynamics properties of the the
$f(R,T)$ gravity theory are explored in [22, 23, 24, 25]. On the other hand,
within the $f(R,T)$ framework, various compact objects are constructed and
studied e.g., wormhole [26] and compact stars [27, 28].

Black holes are ones of the most fundamental objects in the universe. They
play a crucial role in almost all relativistic gravitational field theories.
The detection of gravitational waves [29] and the first capture of black
hole’s image [30] marked the beginning of black hole’s astronomy era. This
makes black holes extremely important in astrophysics research nowadays. Black
holes are solutions of relativistic gravitational field equations. According
to GR, there is an essential singularity hidden behind each black hole. The
regular black hole proposed by Bardeen [31] offers a new possibility to obtain
black hole without a singularity. Later, it is shown that regular black holes
are the solutions of Einstein’s gravity coupled to nonlinear electrodynamics
[32] and the Bardeen black hole can be regarded as a nonlinear magnetic
monopole [33]. The Bardeen black hole is later extended to include a
cosmological constant [34]. Charged regular black holes with various mass
functions are studied in [35]. In addition, a modification of the Reissner-
Nordström black hole yields regular charged black hole, and its entropy obeys
Bekenstein’s area law [36]. We refer interested readers to ref [37] for a
recent review on regular black holes with nonlinear electrodynamics sources.
Beyond GR, the regular black holes with nonlinear electrodynamics are
extensively explored e.g. in Einstein-Gauss-Bonnet theory [38, 39, 40] and
$f(R)$ gravity [41].

In $f(R,T)$ gravity, an exact black hole solution surrounded by anisotropic
fluid is explored [42]. The energy conditions for each particular equation of
state parameter $w$ are discussed in [42]. This prompts a research question
whether there are other black hole solutions in the $f(R,T)$ gravity. Thus, in
this work, we construct asymptotically flat, static spherically symmetric
regular black holes within the framework of $f(R,T)$ gravity. There are two
approaches to obtain the black hole solutions. Firstly, we specifically choose
a mass function that yields a regular black hole, and find the corresponding
the nonlinear electrodynamics Lagrangian ($L_{NED}$). Secondly, we specify the
$L_{NED}$, and find the corresponding mass function. From both approaches, we
obtain novel magnetically charged regular black holes. Remarkably, from the
second approach, we obtain a metric function that can be considered as a
generalization of Bardeen and Hayward black hole [33, 43, 44]. Then, we
analyse the null, weak and strong energy conditions of these solutions. The
quasinormal modes and the eikonal limit of these black holes are also
investigated.

This paper is organized as follows. In Sec II, we discuss the $f(R,T)$ gravity
coupled to nonlinear electrodynamics. The modified field equation is derived
and corresponding energy-momentum tensor is given. Then, modified field
equations are solved and the regular black holes are explored in Sec III.
Then, we discuss the energy conditions in Sec IV . We study quasinormal modes
and the eikonal limit in Sec V. Lastly, we summarize our results and discuss
possible extensions of this work in Sec VI.

##  II Basic equations

We consider $f(R,T)$ gravity coupled to nonlinear electrodynamics (NED). This
theory is described by

$\displaystyle S$
$\displaystyle=\frac{1}{2}\int\sqrt{-g}d^{4}x~{}f(R,T)+\int\sqrt{-g}d^{4}x~{}L_{NED},$
(1)

where $f(R,T)$ is an arbitrary function of the Ricci scalar $R$ and the trace
$T$ of the energy-momentum tensor of the matter $T_{\mu\nu}$. The nonlinear
electrodynamics Lagrangian is given by $L_{NED}(F)$ where
$F=-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}$. The Faraday-Maxwell tensor is defined in
term of the gauge potential
$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$.

Varying this action with respect to $\delta g^{\mu\nu}$, yields the modified
Einstein field equation

$\displaystyle\mathcal{G}_{\mu\nu}\equiv
f_{R}R_{\mu\nu}+\left(g_{\mu\nu}\Box-\nabla_{\mu}\nabla_{\nu}\right)f_{R}-\frac{1}{2}fg_{\mu\nu}$
$\displaystyle=T_{\mu\nu}-f_{T}(T_{\mu\nu}+\Theta_{\mu\nu}),$ (2)

where $f_{R}=\frac{\partial f}{\partial R},f_{T}=\frac{\partial f}{\partial
T}$ and $\Box=\nabla_{\alpha}\nabla^{\alpha}$. The energy-momentum tensor
$T_{\mu\nu}$ and $\theta_{\mu\nu}$ are computed from

$\displaystyle T_{\mu\nu}$
$\displaystyle\equiv-\frac{2}{\sqrt{-g}}\frac{\delta\left(\sqrt{-g}L_{NED}\right)}{\delta
g^{\mu\nu}},$ (3) $\displaystyle\theta_{\mu\nu}$ $\displaystyle\equiv
g^{\alpha\beta}\frac{\delta T_{\alpha\beta}}{\delta g^{\mu\nu}}.$ (4)

With nonlinear electrodynamics sources, the explicit forms of $T_{\mu\nu}$ and
$\Theta_{\mu\nu}$ are expressed as

$\displaystyle T_{\mu\nu}$
$\displaystyle=g_{\mu\nu}L_{NED}+L_{F}F_{\mu\gamma}{F_{\nu}}^{\gamma},$ (5)
$\displaystyle\Theta_{\mu\nu}$
$\displaystyle=-g_{\mu\nu}L_{NED}-F_{\mu\gamma}{F_{\nu}}^{\gamma}\left[\frac{L_{FF}}{2}F_{\rho\sigma}F^{\rho\sigma}+L_{F}\right],$
(6)

where $L_{F}=\partial L_{NED}/\partial F$ and
$L_{FF}=\partial^{2}L_{NED}/\partial F^{2}$. Moreover, taking the trace of
$\eqref{EFE}$ gives the following

$\displaystyle\Box f_{R}$
$\displaystyle=\frac{1}{3}\left(T-f_{T}(T+\Theta)+2f-f_{R}R\right),$ (7)

where $T\equiv g_{\mu\nu}T^{\mu\nu}$ and $\Theta\equiv
g_{\mu\nu}\Theta^{\mu\nu}$. The equation of motion of the gauge field is

$\displaystyle\partial_{\mu}\left[\sqrt{-g}\left(4f_{T}L_{FF}F-L_{F}\right)F^{\mu\nu}\right]$
$\displaystyle=0.$ (8)

Now, we consider a static spherically symmetric solution. The line element
written in Schwarzschild-like coordinate is given by

$\displaystyle
ds^{2}=-A(r)dt^{2}+B(r)dr^{2}+r^{2}\left(d\theta^{2}+\sin^{2}\theta
d\phi^{2}\right).$ (9)

We also consider a purely magnetic ansatz of the Faraday-Maxwell tensor [33]

$\displaystyle F^{\theta\phi}$ $\displaystyle=\frac{q_{m}\csc\theta}{r^{4}},$
(10)

where $q_{m}$ is an integration constant that can be interpreted as the
magnetic charge of the source. With this choice, the invariant $F$ is
$-\frac{q_{m}^{2}}{2r^{4}}$. One can show that the ansatz (10) satisfies the
equation of motion (8).

##  III Solving the modified field equations

Here we consider the $f(R,T)=R+2\beta T$ where $\beta$ is an arbitrary
constant. Together with purely magnetic field strength (10), the modified
field equations i.e.,
$\mathcal{G}{{}_{\mu}}^{\nu}=T{{}_{\mu}}^{\nu}-f_{T}\left(T{{}_{\mu}}^{\nu}+\Theta{{}_{\mu}}^{\nu}\right)$,
are

$\displaystyle-\frac{1}{r^{2}}+\frac{1}{Br^{2}}-\frac{B^{\prime}}{B^{2}r}$
$\displaystyle=\left(1+4\beta\right)L_{NED}+\frac{2\beta
q_{m}^{2}}{r^{4}}L_{F},$ (11)
$\displaystyle-\frac{1}{r^{2}}+\frac{1}{Br^{2}}-\frac{A^{\prime}}{AB^{2}r}$
$\displaystyle=\left(1+4\beta\right)L_{NED}+\frac{2\beta
q_{m}^{2}}{r^{4}}L_{F},$ (12)
$\displaystyle\frac{A^{\prime\prime}}{2AB}-\frac{A^{\prime}B^{\prime}}{4AB^{2}}+\frac{A^{\prime}}{2ABr}-\frac{A^{\prime
2}}{4A^{2}B}-\frac{B^{\prime}}{2B^{2}r}$
$\displaystyle=\left(1+4\beta\right)L_{NED}+\frac{q_{m}^{2}}{r^{4}}\left(1+2\beta\right)L_{F}+\frac{2\beta
q_{m}^{4}}{r^{8}}L_{FF},$ (13)

where prime denotes derivative with respect to $r$. The first two equations
imply that $A=B^{-1}$. Therefore, the remaining field equations are

$\displaystyle\frac{A^{\prime}}{r}+\frac{A}{r^{2}}-\frac{1}{r^{2}}$
$\displaystyle=\left(1+4\beta\right)L_{NED}+\frac{2\beta
q_{m}^{2}}{r^{4}}L_{F},$ (14) $\displaystyle
A^{\prime\prime}+\frac{2A^{\prime}}{r}$
$\displaystyle=2\left(1+4\beta\right)L_{NED}+\frac{2q_{m}^{2}}{r^{4}}\left(1+2\beta\right)L_{F}+\frac{4\beta
q_{m}^{4}}{r^{8}}L_{FF},$ (15)

In addition, the Ricci scalar is

$\displaystyle R$
$\displaystyle=-\left(A^{\prime\prime}+\frac{4A^{\prime}}{r}+\frac{2A}{r^{2}}-\frac{2}{r^{2}}\right).$
(16)

Substituting this into the trace of the modified field equations (7) allows us
to eliminate $A^{\prime\prime}$ in (15). After eliminating $A^{\prime\prime}$,
we find that (14) and (15) are identical. Thus, we are left with a single
first order ordinary differential equation (recalls that $F=F(r)$)

$\displaystyle m^{\prime}(r)$
$\displaystyle=-\frac{r^{2}}{2}\left(1+4\beta\right)L_{NED}-\frac{\beta
q_{m}^{2}}{r^{2}}L_{F},$ (17)

where the mass function $m(r)$ is defined via $A(r)\equiv 1-2m(r)/r$. There
are two ways to solve this equation. Firstly, we may choose a particular form
of $m(r)$, then solving (17) for $L_{NED}$. Secondly, we fix the form of NED
Lagrangian and solve for the mass function $m(r)$.

Before solving for a new solution, let’s examine the consistency of (17). We
consider the case where $L_{NED}=-F,L_{F}=-1$. Therefore (17) can be solved as

$\displaystyle m(r)$ $\displaystyle=-2M+\frac{q_{m}^{2}}{4r},$ (18)

where $M$ is an integration constant. By letting $q_{m}=2Q_{m}$, we obtain a
special case of dyonic Reissner-Nordström black hole [45]. We remark that when
the Lagrangian matter reduces to the $U(1)$ electromagnetic, $\beta$
automatically disappears from $f(R,T)=R+2\beta T$ since the energy-momentum
tensor is traceless.

###  III.1 Fixed mass function

In this subsection, we shall solve the modified field equation (17) for the
spherically symmetric regular black hole solution. We choose the mass function
to be in the form

$\displaystyle m(r)=Me^{-\frac{q_{m}^{2}}{2Mr}},.$ (19)

Here $M$ is a constant parameter and $q_{m}$ is the charge of the regular
black hole. We remark that this form of mass function is considered to obtain
regular black holes within the context of GR [46, 35, 36] and $f(R)$ [41]
gravity coupled to NED.

Since $\displaystyle{\lim_{r\to\infty}\frac{m(r)}{r}=M}$, one can interpret
the constant $M$ as the black hole’s mass. The black hole’s event horizons are
determined by $A(r_{h})=0$, and the location of the outer event horizon is
given by

$\displaystyle r_{h}$
$\displaystyle=-\frac{q_{m}^{2}}{2M\Omega(-\frac{q_{m}^{2}}{4M^{2}})},$ (20)

where $\Omega(z)$ is the omega function or the Lambert W function. This mass
function (19) allows for three possible outcomes regarding number of the
horizons, i.e., two positive real roots (inner and outer horizon), one
degenerated root (extremal case) and no real root (horizonless case). Remark
that, throughout this work, we shall particularly focus on the first two
cases. The behaviour of $A(r)$ are shown in Fig 1 for three possible
solutions. It can be observed from the figure that at small $r$, these
solutions are finite i.e., $A(r)\sim 1$. The solutions are clearly
asymptotically flat since $A(r)\to 1$ as $r\to\infty$. As $q_{m}$ increases,
the minimum value of $A$ increases until $A(r)>0$ for all $r$.

Figure 1: The metric function $A$ is plotted against radial coordinate $r$.
The location of zeros indicates the location of black hole’s event horizons
for $M=1$ and $q_{m}=0.8$ (red), $q_{m}=1$ (green), $q_{m}=1.213$ (blue) and
$q_{m}=1.5$ (black).

The regularity of the solution can be observed by considering two curvature
scalar quantities, i.e., the Ricci scalar $R$ and the Kretchmann scalar $K$.
For the mass function (19), we obtain

$\displaystyle R$
$\displaystyle=\frac{e^{-\frac{q_{m}^{2}}{2Mr}}q_{m}^{4}}{2Mr^{5}},$ (21)
$\displaystyle K$
$\displaystyle=R_{\mu\nu\sigma\rho}R^{\mu\nu\sigma\rho}=\frac{e^{-\frac{q_{m}^{2}}{Mr}}}{4M^{2}r^{10}}\left(q_{m}^{8}-16Mq_{m}^{6}r+96M^{2}q_{m}^{4}r^{2}-192M^{3}q_{m}^{2}r^{3}+192M^{4}r^{4}\right).$
(22)

In Fig 2, we display example plots of the Ricci and the Kretchmann scalars.
The curvature scalars are finite everywhere for various values of $q_{m}$.
Moreover, $R$ and $K$ behave as $O(r^{-5}),O(r^{-6})$, respectively as
$r\to\infty.$ In addition, the maximum value of $R$ locates at
$r=\frac{q_{m}^{2}}{10M}$. On the other hand, the radius renders $K_{max}$ is
not trivial. For instance, when $M=1$ and $q_{m}=0.8$, $R_{max}$ is $1,285$ at
$r=0.064$ while $K_{max}$ is $327,982$ at $r=0.057.$ We emphasize that the
mass function (19), together with these scalar curvatures, is already
considered in GR and $f(R)$ gravity coupled to NED [35, 41].

Figure 2: Left: The Ricci scalar, Right: The Kretechmann scalar for $M=1$,
$q_{m}=0.8$ (red), $q_{m}=1$ (green), $q_{m}=1.213$ (blue) and $q_{m}=1.5$
(Black). The sub-figures show the behaviour of $R$ and $K$ at small $r$.

To obtain $L_{NED}$, we substitute (19) into (17) and solve for $L_{NED}$.
Thus, we obtain

$\displaystyle L_{NED}(F)$
$\displaystyle=\mathcal{C}F^{1+1/4\beta}+\frac{2F}{\beta}Ei_{x}(y),$ (23)

where $\mathcal{C}$ is an integration constant and $Ei_{x}(y)$ is the
exponential integral function where $x\equiv 1+1/\beta$ and
$y\equiv-\frac{(-1)^{3/4}q_{m}^{3/2}F^{1/4}}{2^{3/4}M}$. Since $F$ is
negative, this restricts the value of $\beta$ i.e., $1+\frac{1}{4\beta}=n$
where $n$ is an integer. In Fig 3, we illustrate the behaviour of $L_{NED}$ as
a function of invariant $F$. These plots demonstrate clearly a modification of
standard Maxwell Lagrangian. As can be seen from the plots, $L_{NED}$
approaches zero as $F\to 0$. With a given $\mathcal{C}$ and requiring that
$L_{NED}$ should be a real value, one can show that $L_{NED}\sim F+O(F^{5/4})$
at small $F$. Interestingly, the no-go theorem states that the Einstein field
equation couples to Lagrangian with the Maxwell behaviour at small $F$ (i.e.,
$L\to 0,L_{F}\to 1$ as $F\to 0$) does not admit static spherically symmetric
purely electric solution with a regular centre [47]. Let us remark that, our
attempts to find regular black holes with electric charge are not successful.
This is because with the purely electric gauge potential, the field equations
reduces to a much more complicated second order differential equation
comparing to (14).

Figure 3: Plot of $L_{NED}$ as a function of $F$ with $\mathcal{C}=1,~{}M=1.$
Left: for fix $q_{m}=0.1$ and various values of $\beta$. Right: for fixed
$\beta=\frac{1}{4}$ and various values of $q_{m}$

###  III.2 Fixed Lagrangian

Here, we solve (17) for the mass function when the matter Lagrangian is fixed.
The Lagrangian of nonlinear electrodynamics is chosen to be

$\displaystyle L_{NED}(F)$
$\displaystyle=-\frac{2a}{\alpha}\frac{\left(-4\alpha
F\right)^{(b+3)/4}}{\left(1+\left(-4\alpha F\right)^{b/4}\right)^{1+a/b}},$
(24)

where $a,b$ and $\alpha$ are positive-valued arbitrary constants. This
Lagrangian is adopted from [44, 48] where the authors construct regular black
holes in Einstein-NED and Einstein Cubic gravity, respectively. By inserting
this Lagrangian into (17), the following mass function is obtained

$\displaystyle m(r)$
$\displaystyle=M-\frac{q^{3}}{\alpha}(1+\beta)+\frac{q^{3}}{\alpha}\mathcal{Q}^{-a/b}(1+\beta)+\frac{aq^{3}}{\alpha}\mathcal{Q}^{-1-a/b}\left(\mathcal{Q}-1\right)\beta,$
(25)

where $\mathcal{Q}(r)\equiv 1+\left(\frac{q}{r}\right)^{b}$. $M$ is the
gravitational mass and $q$ is an integration constant related to the magnetic
charge $q_{m}=\frac{q^{2}}{\sqrt{2\alpha}}$. As it was pointed out in [49], we
may define the effective mass $M_{eff}$ as the difference between
gravitational mass $M$ and the magnetically induced mass
$M_{em}=\frac{q^{3}}{\alpha}\left(1+\beta\right)$, i.e., $M_{eff}=M-M_{em}$.
Then, regular black hole is obtained by letting $M=M_{em}$. Therefore, the
metric function of regular black hole in $f(R,T)$ gravity coupled to nonlinear
electrodyanmics source is

$\displaystyle A(r)$ $\displaystyle=1-\frac{2q^{3}}{\alpha
r}\mathcal{Q}^{-a/b}\left[(1+\beta)+a\beta\frac{\left(\mathcal{Q}-1\right)}{\mathcal{Q}}\right]$
$\displaystyle=1-\frac{2q^{3}}{\alpha}\left(r^{b}+q^{b}\right)^{-a/b}\left[(1+\beta)+\frac{a\beta
q^{b}}{r^{b}+q^{b}}\right]r^{a-1},$ (26)

From (26), it appears that to avoid the singularity, one must take $a\geq 1$.
But a closer investigation on the Ricci and the Kretechmann scalar reveals
that to ensure the regularity of the solution as $r\to 0$, $a$ must be equal
or greater than three ($a\geq 3$). This is demonstrated in Fig 4. We observe
that both scalar curvatures diverge as $r\to 0$ for $a<3$. The leading order
term of both scalar curvatures are

$\displaystyle R$ $\displaystyle\sim r^{a-1},~{}~{}~{}~{}K\sim r^{2a-6}.$ (27)

This agrees with the results found in [44]. For $a=4$, the maximum value of
$R$ and $K$ are $12.11$ and $24.85$ for $r=0.3$ and $r=0.29$ respectively. For
the remaining part of this work, we consider only the case where $a\geq 3$.

Figure 4: The Ricci scalar (Left) and the Kretchmann scalar (Right) as a
function of $r$ with varying $a$ and $\beta=0.1,q=0.9,\alpha=1$ and $b=2$.

Now, we consider asymptotic structures of $A,R$ and $K$. As $r\to\infty$, we
find that for $a\geq 3$

$\displaystyle A$ $\displaystyle\sim 1-\frac{2q^{3}}{\alpha
r}(1+\beta)+O\left(\frac{1}{r^{b+1}}\right),$ (28) $\displaystyle R$
$\displaystyle\sim
O\left(\frac{1}{r^{5}}\right),~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{for $b\leq
2$},$ (29) $\displaystyle\sim
O\left(\frac{1}{r^{b+3}}\right),~{}~{}~{}~{}~{}~{}\text{for $b>2$},$ (30)
$\displaystyle K$ $\displaystyle\sim O\left(\frac{1}{r^{6}}\right).$ (31)

The leading order of $A$ suggests that the solution (26) is asymptotically
flat while the others display the regularity of the scalar curvatures at large
$r$. The location of the black hole’s event horizon is subtle without
specifying $a$ and $b$. For the sake of demonstration, we consider three
particular cases, i.e., (i) $a=3,b=2$, (ii) $a=3,b=3$ and (iii) $a=4,b=2$. The
first two cases are chosen such that the Lagrangian (24) gives rise to the
Bardeen-like and Hayward-like solutions [33, 43, 44]. The regular black holes
for (i-iii) in $f(R,T)$ gravity are
(i) $a=3,b=2$

$\displaystyle A_{B}(r)$
$\displaystyle=1-\frac{2q^{3}r^{4}}{\alpha\left(r^{2}+q^{2}\right)^{5/2}}\left[(1+\beta)+\frac{q^{2}}{r^{2}}(1+4\beta)\right],$
(32)

(ii) $a=3,b=3$

$\displaystyle A_{H}(r)$
$\displaystyle=1-\frac{2q^{3}r^{5}}{\alpha\left(r^{3}+q^{3}\right)^{2}}\left[(1+\beta)+\frac{q^{3}}{r^{3}}(1+4\beta)\right],$
(33)

(iii) $a=4,b=2$

$\displaystyle A(r)$
$\displaystyle=1-\frac{2q^{3}r^{5}}{\alpha\left(r^{2}+q^{2}\right)^{3}}\left[(1+\beta)+\frac{q^{2}}{r^{2}}(1+5\beta)\right].$
(34)

Figure 5: Example plots of $A(r)$ as a function of $r$ for $\alpha=1$. Left:
the Bardeen-like solutions with $q=1.2$. Center: the Hayward-like solutions
with $q=1.2$. Right: the regular black holes with $a=4,b=2,\beta=0.5$ for
various values of $q$.

As $\beta\to 0$, the solutions $A_{B}$ and $A_{H}$ become regular black holes
in general relativity i.e., the Bardeen and the Hayward solutions,
respectively. The behaviours of these solutions are illustrated in Fig 5. The
Bardeen-like and the Hayward-like solutions are shown for varying $\beta$. For
the right figure, we fix $\beta=0.5$ and vary $q$ instead. It can be seen that
extremal black holes and horizonless solution are also possible by varying
$\beta$ or $q$. For example, in the third case, the two horizons coincide at
$r=1.157$ with $q=0.85617015$. Moreover, it is clear that they are
asymptotically flat as $A\to 1$ at large $r$.

##  IV Energy conditions

In this section, we consider null, weak and strong energy conditions (NEC,
WEC, SEC) of the solutions discussed in the previous section. To consider the
energy conditions in $f(R,T)$ gravity, let us re-write (2) as

$\displaystyle R_{\mu\nu}-\frac{1}{2}Rg_{\mu\nu}$
$\displaystyle=f^{-1}_{R}\left[T_{\mu\nu}-f_{T}\left(T_{\mu\nu}+\Theta_{\mu\nu}\right)-\left(g_{\mu\nu}\Box-\nabla_{\mu}\nabla_{\nu}\right)f_{R}+\frac{1}{2}g_{\mu\nu}\left(f-Rf_{R}\right)\right],$
$\displaystyle\equiv T^{(eff)}_{\mu\nu}.$ (35)

where we have defined the effective energy-momentum tensor
$T^{(eff)}_{\mu\nu}$. We identify ${T^{(eff)0}}_{0}=-\rho^{(eff)}$,
${T^{(eff)1}}_{1}=p^{(eff)}_{1},{T^{(eff)2}}_{2}=p^{(eff)}_{2}$ and
${T^{(eff)3}}_{3}=p^{(eff)}_{3}$. The energy conditions in $f(R)$ gravity
coupled to NED are discussed in [41]. In addition, the energy conditions of
$f(R,T)$ gravity have been addressed properly in [20, 21, 50]. These are

$\displaystyle NEC$ $\displaystyle\mathrel{\mathop{\mathchar
58\relax}}\rho^{(eff)}+p^{(eff)}_{1,2,3}\geq 0,$ (36) $\displaystyle WEC$
$\displaystyle\mathrel{\mathop{\mathchar 58\relax}}\rho^{(eff)}\geq
0,\rho^{(eff)}+p^{(eff)}_{1,2,3}\geq 0,$ (37) $\displaystyle SEC$
$\displaystyle\mathrel{\mathop{\mathchar
58\relax}}\rho^{(eff)}+p^{(eff)}_{1}+p^{(eff)}_{2}+p^{(eff)}_{3}\geq
0,\rho^{(eff)}+p^{(eff)}_{1,2,3}\geq 0.$ (38)

To clarify the notation, $\rho^{(eff)}+p^{(eff)}_{1,2,3}\geq 0$ is simply
$\rho^{(eff)}+p^{(eff)}_{i}\geq 0$ for $i=1,2,3$ separately. In this model,
the non-vanishing diagonal components of the effective energy-momentum tensor
are given explicitly by

$\displaystyle\rho^{(eff)}$ $\displaystyle=-\left(\beta T+L_{NED}\right),$
(39) $\displaystyle p^{(eff)}_{1}$ $\displaystyle=\left(\beta
T+L_{NED}\right),$ (40) $\displaystyle p^{(eff)}_{2}$
$\displaystyle=p^{(eff)}_{3}=L_{NED}+\frac{q_{m}^{2}}{r^{4}}L_{F}+\beta\left(T+\frac{2q_{m}^{4}}{r^{8}}L_{FF}\right).$
(41)

Therefore, the energy conditions mentioned above reduce to

$\displaystyle NEC$ $\displaystyle\mathrel{\mathop{\mathchar
58\relax}}\rho^{(eff)}+p^{(eff)}_{1,2}\geq 0,$ (42) $\displaystyle WEC$
$\displaystyle\mathrel{\mathop{\mathchar 58\relax}}\rho^{(eff)}\geq
0,\rho^{(eff)}+p^{(eff)}_{1,2}\geq 0,$ (43) $\displaystyle SEC$
$\displaystyle\mathrel{\mathop{\mathchar 58\relax}}2p^{(eff)}_{2}\geq
0,\rho^{(eff)}+p^{(eff)}_{1,2}\geq 0,.$ (44)

Overall, we have four distinct inequalities. These will be considered in the
following subsection for each regular black holes.

###  IV.1 Energy conditions I

Now, we consider the Lagrangian (23). We also assume that $q_{m}$ and $M\geq
0$. It turns out that the energy conditions demand the following

$\displaystyle NEC_{2}~{}\&~{}WEC_{3}~{}\&~{}SEC_{3}\mathrel{\mathop{\mathchar
58\relax}}\frac{e^{-q_{m}^{2}/2Mr}q_{m}}{r}\left(8Mr-q_{m}^{2}\right)\geq 0,$
(45) $\displaystyle SEC_{1}\mathrel{\mathop{\mathchar
58\relax}}\frac{e^{-q_{m}^{2}/2Mr}q_{m}}{r}\left(4Mr-q_{m}^{2}\right)\geq 0.$
(46)

Note that, SEC1 and SEC3 refer to $2p^{(eff)}_{2}\geq 0$ and
$\rho^{(eff)}+p^{(eff)}_{2}\geq 0$ respectively. The NEC1, WEC1, WEC2 and SEC2
are automatically satisfied. The NEC2, WEC3 and SEC3 share similarities.
Moreover, the SEC1 provides another constraint on radial coordinate $r$ (46).
However, all the energy conditions are satisfied simultaneously in a region
$r\geq\frac{q_{m}^{2}}{4M}$.

For parameter set chosen in Fig 1, we find that the NEC and the WEC are
violated in the region $r<\\{0.08,0.125,0.184\\}$ for $q_{m}=0.8,1$ and
$1.213$ respectively. However, these radii are much smaller than the inner
event horizons. Thus, in these cases, the NEC and the WEC are satisfied in the
exterior regions of the black holes. In contrast, the SEC1 is found to be
violated in a region between inner and outer horizon in $q_{m}=0.8$ and $1$
cases. While the SEC3 (45) for these two cases requires $r\geq 0.08,0.125$,
the inner horizon locates at $r=0.110$ and $r=0.232$ respectively. Therefore,
the SEC3 holds very well inside of the inner horizon. However, in the near
extremal limit $q_{m}=1.213$, the SEC holds between two horizons. For these
parameter sets, we find that all the energy conditions hold outside the outer
event horizon. The NEC, WEC and SEC are not met in a region deep inside the
black holes.

###  IV.2 Energy conditions II

For the Lagrangian (24), the energy conditions become complicated and lengthy
without specifying $a$ and $b$. For this reason, we explicitly discuss the
energy conditions for $a=4$ and $b=2$ case only. Similarly, the NEC2, WEC3 and
SEC3 are identical. The energy conditions are, therefore,

$\displaystyle WEC_{1}\mathrel{\mathop{\mathchar
58\relax}}r^{2}(\beta-1)+q^{2}(1+5\beta)\geq 0,$ (47) $\displaystyle
SEC_{1}\mathrel{\mathop{\mathchar
58\relax}}r^{4}(1-\beta)+10q^{2}r^{2}\beta-q^{4}(1+5\beta)\geq 0,$ (48)
$\displaystyle NEC_{2}~{}\&~{}WEC_{3}~{}\&~{}SEC_{3}\mathrel{\mathop{\mathchar
58\relax}}5r^{4}(1-\beta)+2q^{2}r^{2}(2+19\beta)-q^{4}(1+5\beta)\geq 0,$ (49)

where we replace $q_{m}$ with $q$ to match with the notation used in
subsection III.2. The NEC1, WEC2 and SEC2 are naturally satisfied. In Fig 6
and 7, we display the energy conditions as a function of $r$ for fixed $\beta$
and $q$, respectively. The energy conditions are violated if these curves
become negative. When $q=0.86$, the WEC1 holds continuously from the origin
toward the exterior region. However, at certain radius outside of the black
hole, the WEC1 is violated (x-intercept is at $r=2.28$). In contrast, the
other energy conditions hold from the certain radius inside the black hole all
the way toward black hole’s exterior. As we move away from near-extremal
scenario, the EC1 becomes negative relatively faster than the near-extremal
case (x-intercept is at $r=3.97$). These are shown in Fig 6.

Figure 6: Plots of the energy conditions against $r$ for fixed $\beta=0.5$.
The plot legends EC1,EC2,EC3 correspond to the energy conditions (47)–(49),
respectively. Left: $q=0.86$, Right: $q=1.5$. The dashed vertical lines
represent the location of inner and outer horizons.

Figure 7: Plots of the energy conditions against $r$ for fixed $q=1$. The plot
legends $EC_{1},EC_{2},EC_{3}$ correspond to the energy conditions (47)–(49),
respectively. Left: $\beta=0.5$, Centre: $\beta=1$, Right: $\beta=2$. The
dashed vertical lines represent the location of inner and outer horizons.

In addition, we explore how $\beta$ affects the energy conditions in Fig 7.
The left figure shows similar behaviours of the energy conditions as mentioned
in the previous figure. The WEC1 is violated just prior to the outer event
horizon (x-intercept is at $r=2.65$ while the outer horizon is at $r=2.69$).
While the others are positive from the inner horizon. When $\beta=1$, the EC1
becomes positive constant in $r$, therefore it is always satisfied. The other
two remain positive right after the inner horizon. Remarkably, the energy
conditions change dramatically for $\beta=2$. The EC2 and EC3 are positive for
some particular region inside the black hole’s outer horizon before rapidly
become negative. The EC1, on the other hand, holds throughout spatial
coordinate $r$.

For this particular case i.e., $a=4$ and $b=2$, we find that the NEC and SEC
are easily met at the black hole’s exterior while the WEC will be violated at
certain radius. However, appropriate selections of parameters can possibly
make all the energy conditions satisfied.

##  V Quasinormal modes

A massive scalar field $(\Phi)$ on curved spacetime is described by the Klein-
Gordon equation

$\displaystyle\nabla_{\gamma}\nabla^{\gamma}\Phi-\mu^{2}\Phi$
$\displaystyle=0,$ (50)

where $\mu$ is the scalar field’s mass. In a spherical symmetric spacetime,
the scalar field can be expressed as

$\displaystyle\Phi(t,r,\theta,\phi)$ $\displaystyle=\frac{R(r)}{r}e^{-i\omega
t}Y(\theta,\phi),$ (51)

where $Y(\theta,\phi)$ is the spherical harmonics. Under the spacetime metric
(9) (with $B=A^{-1}$), the Klein-Gordon equation takes the form

$\displaystyle\frac{d^{2}R}{dr_{\ast}^{2}}+\left(\omega^{2}-V(r)\right)R$
$\displaystyle=0.$ (52)

The effective potential is

$\displaystyle V(r)$
$\displaystyle=A(r)\left(\mu^{2}+\frac{\ell(\ell+1)}{r^{2}}+\frac{A^{\prime}(r)}{r}\right),$
(53)

where $\ell$ is the spherical harmonic index. Moreover, we have introduced the
tortoise coordinate defined by

$\displaystyle r_{\ast}$ $\displaystyle=\int\frac{dr}{A(r)}.$ (54)

The appropriate boundary conditions that lead to the quaisnormal mode are
purely ingoing at the black hole’s event horizon $r\to r_{h}$ or
$r_{\ast}\to-\infty$ and no incoming flux at infinity $r,r_{\ast}\to\infty$.
The frequencies $\omega$ corresponding to these boundary conditions will be
discrete complex number or quasinormal frequencies. This complex frequency can
be written in the form $\omega=\omega_{R}\pm i\omega_{I}$.

Let’s first consider the effective potential (53) more explicitly. All the
solutions considered in this work are asymptotically flat, therefore
$V\to\mu^{2}$ as $r\to\infty$. Unless $A^{\prime}(r)<0$, the location where
$V$ vanished is only determined by the roots of $A(r)$. The effective
potential for several types of regular black holes are shown in Fig 8-9. For
the mass function (19), the effective potential is illustrated in Fig 8. As
the charge $q_{m}$ increases, the height of $V$ increases. We observe that the
zeroth of $V$ occurs at the location of the black hole’s outer horizon. More
precisely, these potentials have another zeroth located at smaller $r$ which
corresponds to the inner horizon. However, these are not explicitly displayed
in the plots. In the extremal case (black solid line in the left figure), the
potential possesses only one root. The central figure of Fig 8 demonstrates
the effect of harmonic index $\ell$ to the height of the effective potential.
As $\ell$ decreases, the peak of $V$ decreases. The last figure illustrates
that the peak of $V$ increases with scalar field’s mass. Asymptotic value of
$V$ approaches $\mu^{2}$ as expected. Remark that, similar plots are already
explored in the case of scalar perturbations on the Bardeen solution [51].

Figure 8: Plots of effective potential with mass function (19) against $r$. In
these plots, $M$ is fixed to $1$. Left: $\ell=2,\mu=0.1$ for various $q_{m}$,
Centre: $q_{m}=0.5,\mu=0.1$ for various $\ell$, Right: $q_{m}=0.2,\ell=3$ for
various $\mu$.

Figure 9: Plots of effective potential of three types of regular black holes
against $r$. Left: fixed $\ell=2,\mu=0.1,\alpha=1,q=1.4$ for various $\beta$,
Centre: fixed $\beta=0.45,\mu=0.1,\alpha=1,q=0.9$ for various $\ell$, Right:
$\ell=2,\beta=0.2,\alpha=1,q=1.2$ for various $\mu$. The ansatz 1, 2 and 3 are
illustrated by solid, dashed and dot-dashed lines.

Now, we consider the effective potential of the solutions reported in
subsection III.2. We name the Bardeen-like, the Hayward-like and $a=4,b=2$
solution as ansatz 1, 2 and 3, respectively. The potentials are illustrated in
Fig 9. In these plots, the ansatz 1, 2 and 3 are represented by solid, dashed
and dot-dashed lines, respectively. The effects of $\ell$ and $\mu$ on the
effective potential are qualitatively similar to the previous case as
demonstrated in the central and right figures. We notice that the differences
between the ansatz 1, 2 and 3 become more significant as $\ell$ or $\mu$
increases. In contrast, the differences are less apparent as $\beta$
increases. This is shown in the left figure. Moreover, the height of the
potential decreases as $\beta$ increases. Remark that, in the case that
$\beta=0$, the effective potential of the ansatz 1 (red solid line) is plotted
in [52].

###  V.1 The Padé averaged WKB approximation method

To calculate the quasinormal frequencies, we employ the sixth order WKB
approximation technique. With this method, it is possible to obtain the
quasinormal frequencies $\omega$ via the following expression (up to sixth
order) [53]

$\displaystyle\frac{i\left(\omega^{2}-V_{max}\right)}{\sqrt{-2V^{\prime\prime}_{max}}}-\Lambda_{2}-\Lambda_{3}-\Lambda_{4}-\Lambda_{5}-\Lambda_{6}$
$\displaystyle=n+\frac{1}{2},$ (55)

where $V_{max},V^{\prime\prime}_{max}$ are the effective potential and second
derivative of the effective potential with respect to the tortoise coordinate
evaluated at the maximum point of the potential. The overtone number is
denoted by $n$. Iyer and Will find the correction terms up to the third order
(hence $\Lambda_{2},\Lambda_{3}$) [54]. While later, Konoplya find three more
correction terms $\Lambda_{4},\Lambda_{5}$ and $\Lambda_{6}$ which are defined
in [53]. To improve the numerical accuracy, the WKB approximation is extended
to thirteen order including the Padé averaged in [55] where the quasinormal
frequencies of the Schwarzschild and Reissner-Nordström black hole are
reproduced. It turns out that with the Padé averaged technique, many known
results can be reproduced with great accuracy [56]. The Mathematica code for
calculating quasinormal frequencies up to thirteen order WKB with improved
Padé averaged is provided in [57]. Thus, we will implement the code for
computing quasinormal frequencies in this work. We remark that, throughout
this section, the parameter $\alpha$ is substituted by
$\frac{q^{3}}{M}(1+\beta)$.

Table 1: The quasinormal frequencies computed by sixth order WKB with the Padé averaged for $M=1$. Upper: the $\omega$ for solution (19) for $\ell=1$ and $\ell=2$ (in parentheses) with fixed $\mu=0.1$. Lower: the $\omega$ for ansatz 1 for $\beta=0,q=0.76$ and $\ell=2$. Here, we only display $n=0$ mode. $q_{m}$ | Padé averaged WKB | Error estimation | Results from [51]
---|---|---|---
$0.2$ | $0.299367-0.095199i$ | $1.534997\times 10^{-7}$ | $0.2993-0.0953i$
| $(0.490046-0.095900i)$ | $(2.466144\times 10^{-7})$ | $(0.4900-0.0959i)$
$0.4$ | $0.305515-0.095912i$ | $9.141781\times 10^{-8}$ | $0.3055-0.0960i$
| $(0.500232-0.096563i)$ | $(2.412731\times 10^{-7)}$ | $(0.5002-0.0966i)$
$0.6$ | $0.316811-0.097033i$ | $2.341719\times 10^{-6}$ | $0.3168-0.0971i$
| $(0.518964-0.097606i)$ | $(2.100739\times 10^{-7})$ | $(0.5190-0.0976i)$
$0.8$ | $0.335444-0.098326i$ | $7.426319\times 10^{-6}$ | $0.3354-0.0984i$
| $(0.549924-0.098795i)$ | $(1.141331\times 10^{-7})$ | $(0.5499-0.0988i)$
$\mu$ | Padé averaged WKB | Error estimation | Results from [52]
$0.1$ | $0.553916-0.079097i$ | $2.211318\times 10^{-7}$ | $0.553910-0.0791170i$
$0.2$ | $0.560827-0.077603i$ | $9.480210\times 10^{-7}$ | $0.560821-0.0776240i$
$0.3$ | $0.572506-0.074963i$ | $3.869737\times 10^{-7}$ | $0.572499-0.0749856i$
$0.4$ | $0.589209-0.070920i$ | $5.633078\times 10^{-7}$ | $0.589201-0.0709443i$
$0.5$ | $0.611315-0.065028i$ | $1.348216\times 10^{-6}$ | $0.611305-0.0650541i$

As a consistency check, we list the $n=0$ quasinormal frequencies in Table 1.
In this table, we reproduce the results already obtained in refs [51, 52]
which are shown in the rightmost column. We implement the sixth order WKB with
Padé average to compute the quasinormal frequencies of massive scalar
perturbation with spacetime background given by the mass function (19) (the
upper table) and ansatz 1 with $\beta=0$ (the lower table). The error
estimation denotes the root mean square error corresponding to the sixth order
WKB with Padé approximation. The black hole’s mass is set to unity. The upper
table displays $\omega$ as a function of $q_{m}$ for $\ell=1$ and $\ell=2$ (in
parentheses). Both real and imaginary parts increase as the black hole’s
charge increases (in magnitude). The lower table investigates the effect of
scalar field’s mass $\mu$ on the $\omega$ with fixed $q=0.76$ and $\ell=2$.
The real part of $\omega$ increases with $\mu$ whereas the imaginary part
decreases with $\mu$. Apparently, the sixth order Padé averaged WKB
approximation method agrees with those results found earlier.

Now, we shall turn our attention to the QNMs of regular black holes of ansatz
1, 2 and 3, (32)–(34). In Table 2–4, the quasinormal frequencies with
$\ell=0-2$ are displayed as a function of $q$ for Bardeen-like, Hayward-like
and $a=4,b=2$ solutions respectively. For comparison, in these tables, we fix
$M=1,\beta=0.1$ and scalar field’s mass is $0.1$. Since the WKB approximation
works very well when $\ell>n$ [51, 58], we consider $\ell=1,n=0$, $\ell=2,n=0$
and $n=1$ cases. Despite the $\ell=n=0$ case might not be well-approximated by
the WKB method, we will include them in the tables since they are the most
fundamental modes. It turns out that the quasinormal frequencies of the ansatz
1–3 share similar trends. As the black hole’s charge $q$ increases, the real
part of $\omega$ increases while the imaginary part becomes less negative.
Various studies on QNMs of regular black holes also report the similar trend
[52, 59, 51, 60]. With increasing angular index $\ell$, the real part
increases. In contrast, the effect of $\ell$ on $\omega_{I}$ is non-trivial.
At first, the imaginary part decreases (in magnitude) when $\ell$ moves from
zero to one. Later, the imaginary part increases once again as $\ell=2$.
Lastly, both $\omega_{R}$ and $\omega_{I}$ decrease as the overtone number $n$
is larger. We observe that the quasinormal frequencies of these regular black
holes (ansatz 1–3) marginally differ from each other. This should not be
surprised because the effective potentials of these ansatzes (Fig 9) are
nearly identical. Therefore, for the remaining part of this article, we will
particularly focus only on the ansatz 3 for the sake of presentation.

Table 2: The Bardeen-like solution with $M=1,\beta=0.1$ and $\mu=0.1$. $\ell$ | $n$ | $q=0.1$ | $q=0.3$ | $q=0.5$ | $q=0.7$
---|---|---|---|---|---
$0$ | $0$ | $0.113977-0.096346i$ | $0.115682-0.095390i$ | $0.119425-0.093134i$ | $0.123845-0.086839i$
$1$ | $0$ | $0.297824-0.094875i$ | $0.301256-0.094145i$ | $0.308793-0.092224i$ | $0.322230-0.087376i$
$2$ | $0$ | $0.487471-0.095589i$ | $0.493002-0.094832i$ | $0.505246-0.092863i$ | $0.527654-0.087942i$
| $1$ | $0.466262-0.292988i$ | $0.472838-0.290381i$ | $0.487197-0.283678i$ | $0.512242-0.267288i$
Table 3: The Hayward-like solution with $M=1,\beta=0.1$ and $\mu=0.1$. $\ell$ | $n$ | $q=0.1$ | $q=0.3$ | $q=0.5$ | $q=0.7$
---|---|---|---|---|---
$0$ | $0$ | $0.113739-0.096406i$ | $0.113850-0.096230i$ | $0.114325-0.095620i$ | $0.115884-0.094445i$
$1$ | $0$ | $0.297418-0.094953i$ | $0.297639-0.094813i$ | $0.298484-0.094270i$ | $0.300409-0.092945i$
$2$ | $0$ | $0.486816-0.095670i$ | $0.487162-0.095532i$ | $0.488481-0.094994i$ | $0.491519-0.093689i$
| $1$ | $0.465483-0.293268i$ | $0.465938-0.292810i$ | $0.467658-0.291022i$ | $0.471516-0.286672i$
Table 4: The ansatz 3 ($a=4,b=2$) with $M=1,\beta=0.1$ and $\mu=0.1$. $\ell$ | $n$ | $q=0.1$ | $q=0.3$ | $q=0.5$ | $q=0.7$
---|---|---|---|---|---
$0$ | $0$ | $0.114032-0.096285i$ | $0.116337-0.095010i$ | $0.121850-0.091256i$ | $0.122825-0.080233i$
$1$ | $0$ | $0.297963-0.094847i$ | $0.302605-0.093826i$ | $0.313203-0.090768i$ | $0.333047-0.080207i$
$2$ | $0$ | $0.487694-0.095560i$ | $0.495182-0.094504i$ | $0.512499-0.091388i$ | $0.547741-0.080482i$
| $1$ | $0.466529-0.292887i$ | $0.475418-0.289264i$ | $0.495503-0.278719i$ | $0.529632-0.243894i$

Figure 10: The real part and imaginary part of quasinormal frequencies of
regular black holes (ansatz 3) as function of $\beta$ for various values of
$q$. With $M=1,\ell=2,n=0$ and $\mu=0.05$.

Figure 11: The real part and imaginary part of quasinormal frequencies of
regular black holes (ansatz 3) as function of $\ell$ for various values of
$\beta$. With $M=1,q=0.7,n=0$ and $\mu=0.1$. For the imaginary part, the
$\ln(|\omega_{I}|)$ is plotted instead of $\omega_{I}$.

Figure 12: The real part and imaginary part of quasinormal frequencies of
regular black holes (ansatz 3) as function of $\mu$ for various values of $n$.
With $M=1,\beta=0.1,q=0.4$ and $\ell=5$.

We explore how the coupling constant $\beta$ affects the quasinormal
frequencies $\omega$ in Fig 10. In this plot, we choose four particular values
of black hole’s charge $q=0.1,0.5,1.0$ and $1.4$. As $\beta$ increases, the
real part and imaginary part of $\omega$ decrease. The decrease of frequencies
becomes less obvious at larger $\beta$. In addition, the change in $\omega$
can be clearly seen as $q$ increases. As $q$ gets larger, both $\omega_{R}$
and $\omega_{I}$ are smaller (in magnitude). It is worth mentioning that at
$q=1.0$ and $q=1.4$, there are no regular black holes for $\beta<0.76$ and
$\beta<6.2$ respectively.

In Fig 11, the dependence of $\omega$ on spherical harmonic index $\ell$ is
illustrated. The real part of $\omega$ increases monotonically with $\ell$. In
addition, the differences in $\omega_{R}$ between each fixed $\beta$ become
more evident at large $\ell$. To demonstrate the change in $\omega_{I}$, we
express them with $\ln(|\omega_{I}|)$. At small $\ell$, $\omega_{I}$ varies
drastically with $\ell$, but as $\ell$ gets bigger, the change in $\omega_{I}$
becomes less significant. Remark that, the imaginary parts of quasinormal
frequencies are more negative as $\ell$ increases. These trends are also
observed for QNMs of regular black holes [51], Bardeen black holes [52] and
Bardeen de-Sitter black holes [59]. It can be seen from the plots that
$\omega_{R}$ and $\omega_{I}$ decrease as $\beta$ increases which are in
agreement with what was discussed in the previous figure.

In Fig 12, we demonstrate the effect of scalar field’s mass $\mu$ on the
quasinormal frequencies. As the field’s mass increases, the real part of the
frequencies increases, while the imaginary part becomes smaller as $\mu$
increases. The lowest overtone mode ($n=0$) has larger $\omega_{R}$ and
$\omega_{I}$ comparing to the higher overtone modes. We notice that
$\omega_{R}$ increases monotonically along with the $\mu$. In contrast, for
the first overtone ($n=1$) of QNMs of Bardeen black hole [52], its real part
approaches a certain maximum value and then decreases with the scalar field
mass. In addition, a study on massive scalar perturbation on Reissner-
Nordström black hole reveals that it is possible to have arbitrarily long live
modes or quasi-resonance modes as the scalar field’s mass increases [61, 62].
Similarly, our results are expected to respect this behaviour. However, the
WKB approximation method is not sufficient to accurately provide the quasi-
resonance modes [57].

###  V.2 The eikonal limit

When solving for quasinormal frequencies of black holes, various numerical
schemes are applicable. Nevertheless, there is an approximation that provides
a useful formula for quasinormal frequencies with great accuracy. To calculate
QNMs of black holes, one can consider the so-called geometric-optics or
eikonal limit as suggested by Mashhoon and Ferrari [63, 64, 65]. In the
eikonal limit ($\ell\to\infty$), the effective potential (53) is simply

$\displaystyle V_{eik}(r)$ $\displaystyle\approx\frac{A\ell^{2}}{r^{2}}$ (56)

This greatly simplifies the radial wave equation (52). The reduced radial
equation can be solved given that the effective potential $V_{eik}$ satisfies
the quantization condition

$\displaystyle\frac{i\left(\omega^{2}-V_{max}\right)}{\sqrt{-2V^{\prime\prime}_{max}}}$
$\displaystyle=n+\frac{1}{2}.$ (57)

where $V_{max}$ now denotes the maximum point of $V_{eik}$ i.e.,
$\frac{dV_{eik}}{dr_{\ast}}|_{r=r_{0}}=0$. It turns out that the eikonal QNMs
can be expressed as the first order of the WKB formula (55) (e.g. see [66]).
The higher order terms $\Lambda_{2}-\Lambda_{6}$ can be considered as
correction terms to the eikonal limit. Remarkably, it is pointed out in [67]
that QNMs in the eikonal limit can be related to unstable circular null orbit
around black holes in any dimensions. The real part of quasinormal frequency
is determined by the angular velocity at the unstable null geodesics
($\Omega$). The imaginary part is related to the Lyapunov exponent
($\lambda_{L}$) which corresponds to an inverse of instabiltity timescale of
the null orbit. In [68], an upper bound of the Lyapunov exponent of a particle
near the horizon is considered. The upper bound is determined by the surface
gravity at the horizon [68]. Very recently, the violation of the Lyapunov
exponent bound is found for Kerr-Newmann de Sitter black hole [69].

The approximation formula of the quasiormal frquencies in the eikonal limit
can be expressed as [67]

$\displaystyle\omega_{eik}$
$\displaystyle=\Omega\ell-i\left(n+\frac{1}{2}\right)|\lambda_{L}|,$ (58)

where

$\displaystyle\Omega$
$\displaystyle=\left.\sqrt{\frac{A}{r^{2}}}\right|_{r=r_{0}},$ (59)
$\displaystyle\lambda_{L}$
$\displaystyle=\frac{1}{\sqrt{2}}\left.\sqrt{A\left(A^{\prime\prime}-\frac{6A^{\prime}}{r}+\frac{6A}{r^{2}}\right)+A^{\prime
2}}\right|_{r=r_{0}}.$ (60)

Both the angular velocity and the Lyapunov exponent are evaluated at $r_{0}$
and the prime refers to derivative with respect to $r$. We have checked and
confirmed that (58) agrees with the Padé averaged WKB in the limit $\ell\gg
1$.

In Fig 13, we illustrate the behaviour of the angular velocity $(\Omega)$ and
the Lyapunov exponent $(\lambda_{L})$ as a function of $\beta$ and $q$. As can
be seen from the plots, the angular velocity decreases with $\beta$. At the
lower $\beta$, the angular velocity drops significantly comparing to the
larger $\beta$. Notice that, $\Omega$ changes rapidly with $\beta$ at higher
black hole’s charge $q$. In addition, as $\beta$ increases, the Lyapunov
exponent becomes larger before approaching a certain asymptotic value. From
the bottom figure, we observe that the angular velocity (the Lyapunov
exponent) increases (decreases) monotonically with $q$. The behaviour of
$\Omega$ against $q$ is in agreement with those found earlier in [52, 51] for
the Reissner-Nordstrom, the Bardeen black holes and the regular black holes
with exponential mass function. In contrast, our results on the Lyapunov
exponent plotted against $q$ resemble the Bardeen black holes but
substantially different from the Reissner-Nordström and the regular black
holes [52, 51]. Remark that, the Lyapunov exponent of the Reissner-Nordström
black hole increases with $q$ until it reaches its maximum value at some
certain $q$ and then it decreases.

Figure 13: The angular velocity and the Lyapunov exponent as a function of
$\beta$ (top) and $q$ (bottom). The black hole’s mass is fixed to one.

##  VI Conclusions

In this work, we study the $f(R,T)$ gravity coupled to nonlinear
electrodynamics Lagrangian. With purely magnetic component of the gauge field,
asymptotically flat static spherically symmetric black holes with regular
centre are constructed. The black hole solutions are obtained via two
approaches: I.) Fix the mass function and solve for $L_{NED}$ II.) Fix the
$L_{NED}$ and solve for the mass functions. The first approach yields a
functional form of a novel $L_{NED}$. Figure 3 clearly shows the difference
from the standard $U(1)$ electromagnetic Lagrangian. From the second approach,
we find a generalized metric function that can be reduced to the Bardeen and
the Hayward black holes in an appropriate limit. From both approaches, we find
that these charged black holes possess two event horizons without essential
singularities. These are shown in Fig 2 and Fig 4 where the Ricci scalar and
the Kretchmann scalar are plotted.

The energy conditions (null, weak and strong) of these solutions are also
explored. For the regular black holes obtained via the first approach, all the
energy conditions considered here hold in exterior region of the black holes.
The black hole solutions from the second approach respect the null and strong
energy conditions outside the black hole’s outer horizon for small value of
$\beta$. As $\beta$ increases, the null and strong energy conditions are not
guaranteed to hold. In contrast, the weak energy condition is violated at some
certain radius inside the outer horizon for small $\beta$. As $\beta$
increases, the weak energy condition will be satisfied outside the black
hole’s outer horizon.

We investigate a massive scalar perturbation on these regular black holes. The
corresponding quasinormal frequencies are computed via the Padé average WKB
method. For all cases considered in this work, the imaginary parts of the
frequencies are all negative. We find that the real parts of the frequencies
increase with $q,\ell$ and $\mu$ while they decrease with $\beta$. In
addition, the imaginary parts of the frequencies become less negative as $q$
and $\mu$ increase, and become more negative as $\beta$ and $\ell$ increase.
In the eikonal limit, the angular velocity (the Lyapunov exponent) decreases
(increases) with $\beta$. Furthermore, the dependence of $\lambda_{L}$ on $q$
is different from those of the Reissner-Nordström black holes and the regular
black holes with exponential mass function.

There are several ways to extend this work further. Despite what the no-go
theorem states in [47], it is crucial to show whether the $f(R,T)$ admits
electrically charged regular black holes. It is interesting to consider the
thermodynamics properties of regular black holes in $f(R,T)$ gravity. This
problem could be challenging since the first law of thermodynamics in $f(R,T)$
is violated [24]. Moreover, the photon motions around the black holes
discussed in this article are also important since this could lead to the
study of optical appearances of these black holes. Additionally, there are
various forms of the function $f(R,T)$, and it is of great interest to explore
whether they admit an exact black hole solution.

###### Acknowledgements.

This work (Grant No. RGNS 64-217) was supported by Office of the Permanent
Secretary, Ministry of Higher Education, Science, Research and Innovation (OPS
MHESI), Thailand Science Research and Innovation (TSRI) and Silpakorn
university. M. Youk was supported by the Faculty of Science, Silpakorn
University, Thailand, through grant SCSU-STA-2566-11.

## References

  * Nojiri and Odintsov [2011] S. Nojiri and S. D. Odintsov, Unified cosmic history in modified gravity: from F(R) theory to Lorentz non-invariant models, Phys. Rept. 505, 59 (2011), arXiv:1011.0544 [gr-qc] .
  * Coley and Leon [2019] A. Coley and G. Leon, Static Spherically Symmetric Einstein-aether models I: Perfect fluids with a linear equation of state and scalar fields with an exponential self-interacting potential, Gen. Rel. Grav. 51, 115 (2019), arXiv:1905.02003 [gr-qc] .
  * Riess  _et al._ [1998] A. G. Riess _et al._ (Supernova Search Team), Observational evidence from supernovae for an accelerating universe and a cosmological constant, Astron. J. 116, 1009 (1998), arXiv:astro-ph/9805201 .
  * Perlmutter  _et al._ [1999] S. Perlmutter _et al._ (Supernova Cosmology Project), Measurements of $\Omega$ and $\Lambda$ from 42 high redshift supernovae, Astrophys. J. 517, 565 (1999), arXiv:astro-ph/9812133 .
  * Rubin and Ford [1970] V. C. Rubin and W. K. Ford, Jr., Rotation of the Andromeda Nebula from a Spectroscopic Survey of Emission Regions, Astrophys. J. 159, 379 (1970).
  * Sotiriou and Faraoni [2010] T. P. Sotiriou and V. Faraoni, $f(r)$ theories of gravity, Rev. Mod. Phys. 82, 451 (2010).
  * De Felice and Tsujikawa [2010] A. De Felice and S. Tsujikawa, f(R) theories, Living Rev. Rel. 13, 3 (2010), arXiv:1002.4928 [gr-qc] .
  * Nojiri  _et al._ [2017] S. Nojiri, S. D. Odintsov, and V. K. Oikonomou, Modified Gravity Theories on a Nutshell: Inflation, Bounce and Late-time Evolution, Phys. Rept. 692, 1 (2017), arXiv:1705.11098 [gr-qc] .
  * Capozziello  _et al._ [2015] S. Capozziello, F. S. N. Lobo, and J. P. Mimoso, Generalized energy conditions in Extended Theories of Gravity, Phys. Rev. D 91, 124019 (2015), arXiv:1407.7293 [gr-qc] .
  * Cai  _et al._ [2016] Y.-F. Cai, S. Capozziello, M. De Laurentis, and E. N. Saridakis, f(T) teleparallel gravity and cosmology, Rept. Prog. Phys. 79, 106901 (2016), arXiv:1511.07586 [gr-qc] .
  * Beltrán Jiménez  _et al._ [2018] J. Beltrán Jiménez, L. Heisenberg, and T. Koivisto, Coincident General Relativity, Phys. Rev. D 98, 044048 (2018), arXiv:1710.03116 [gr-qc] .
  * Khyllep  _et al._ [2023] W. Khyllep, J. Dutta, E. N. Saridakis, and K. Yesmakhanova, Cosmology in f(Q) gravity: A unified dynamical systems analysis of the background and perturbations, Phys. Rev. D 107, 044022 (2023), arXiv:2207.02610 [gr-qc] .
  * Harko  _et al._ [2011] T. Harko, F. S. N. Lobo, S. Nojiri, and S. D. Odintsov, $f(R,T)$ gravity, Phys. Rev. D 84, 024020 (2011), arXiv:1104.2669 [gr-qc] .
  * Shabani and Farhoudi [2013] H. Shabani and M. Farhoudi, f(R,T) Cosmological Models in Phase Space, Phys. Rev. D 88, 044048 (2013), arXiv:1306.3164 [gr-qc] .
  * Houndjo and Piattella [2012] M. J. S. Houndjo and O. F. Piattella, Reconstructing $f(R,T)$ gravity from holographic dark energy, Int. J. Mod. Phys. D 21, 1250024 (2012), arXiv:1111.4275 [gr-qc] .
  * Jamil  _et al._ [2012a] M. Jamil, D. Momeni, M. Raza, and R. Myrzakulov, Reconstruction of some cosmological models in f(r,t) cosmology, The European Physical Journal C 72, 10.1140/epjc/s10052-012-1999-9 (2012a).
  * Baffou  _et al._ [2015] E. H. Baffou, A. V. Kpadonou, M. E. Rodrigues, M. J. S. Houndjo, and J. Tossa, Cosmological viable $f(R,T)$ dark energy model: dynamics and stability, Astrophys. Space Sci. 356, 173 (2015), arXiv:1312.7311 [gr-qc] .
  * Sharif and Zubair [2014] M. Sharif and M. Zubair, Cosmological reconstruction and stability in $f(R,T)$ gravity, Gen. Rel. Grav. 46, 1723 (2014).
  * Mishra  _et al._ [2018] B. Mishra, S. Tarai, and S. K. Tripathy, Anisotropic cosmological reconstruction in f(r, t) gravity, Modern Physics Letters A 33, 1850170 (2018).
  * Alvarenga  _et al._ [2013] F. G. Alvarenga, M. J. S. Houndjo, A. V. Monwanou, and J. B. C. Orou, Testing some f(R,T) gravity models from energy conditions, J. Mod. Phys. 4, 130 (2013), arXiv:1205.4678 [gr-qc] .
  * Sharif and Zubair [2013a] M. Sharif and M. Zubair, Energy Conditions Constraints and Stability of Power Law Solutions in f(R,T) Gravity, J. Phys. Soc. Jap. 82, 014002 (2013a), arXiv:1210.3878 [gr-qc] .
  * Sharif and Zubair [2012] M. Sharif and M. Zubair, Thermodynamics in f(R,T) Theory of Gravity, JCAP 03, 028, [Erratum: JCAP 05, E01 (2012)], arXiv:1204.0848 [gr-qc] .
  * Jamil  _et al._ [2013] M. Jamil, D. Momeni, and R. Myrzakulov, Wormholes in a viable f(T) gravity, Eur. Phys. J. C 73, 2267 (2013), arXiv:1212.6017 [gr-qc] .
  * Jamil  _et al._ [2012b] M. Jamil, D. Momeni, and R. Myrzakulov, Violation of First Law of Thermodynamics in f(R,T) Gravity, Chin. Phys. Lett. 29, 109801 (2012b), arXiv:1209.2916 [physics.gen-ph] .
  * Sharif and Zubair [2013b] M. Sharif and M. Zubair, Thermodynamic behavior of particular $f(R,T)$-gravity models, J. Exp. Theor. Phys. 117, 248 (2013b), arXiv:1409.7378 [gr-qc] .
  * Banerjee  _et al._ [2023] A. Banerjee, T. Tangphati, and A. Pradhan, Wormholes and energy conditions in $f(R,T)$ gravity, (2023), arXiv:2305.07541 [gr-qc] .
  * Tangphati  _et al._ [2022] T. Tangphati, S. Hansraj, A. Banerjee, and A. Pradhan, Quark stars in f(R,T) gravity with an interacting quark equation of state, Phys. Dark Univ. 35, 100990 (2022).
  * Tangphati  _et al._ [2023] T. Tangphati, G. Panotopoulos, A. Banerjee, and A. Pradhan, Charged compact stars with color–flavor-locked strange quark matter in f(R,T) gravity, Chin. J. Phys. 82, 62 (2023).
  * Abbott  _et al._ [2016] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Observation of Gravitational Waves from a Binary Black Hole Merger, Phys. Rev. Lett. 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Aki [2019] 875, L1 (2019).
  * Bardeen [1968] J. Bardeen, Non-singular general-relativistic gravitational collapse, in conference proceedings of gr5 (1968).
  * Ayon-Beato and Garcia [1998] E. Ayon-Beato and A. Garcia, Regular black hole in general relativity coupled to nonlinear electrodynamics, Phys. Rev. Lett. 80, 5056 (1998), arXiv:gr-qc/9911046 .
  * Ayon-Beato and Garcia [2000] E. Ayon-Beato and A. Garcia, The Bardeen model as a nonlinear magnetic monopole, Phys. Lett. B 493, 149 (2000), arXiv:gr-qc/0009077 .
  * Fernando [2017] S. Fernando, Bardeen–de Sitter black holes, Int. J. Mod. Phys. D 26, 1750071 (2017), arXiv:1611.05337 [gr-qc] .
  * Balart and Vagenas [2014] L. Balart and E. C. Vagenas, Regular black holes with a nonlinear electrodynamics source, Phys. Rev. D 90, 124045 (2014), arXiv:1408.0306 [gr-qc] .
  * Culetu [2015] H. Culetu, On a regular charged black hole with a nonlinear electric source, Int. J. Theor. Phys. 54, 2855 (2015), arXiv:1408.3334 [gr-qc] .
  * Bronnikov [2022] K. A. Bronnikov, Regular black holes sourced by nonlinear electrodynamics, (2022), arXiv:2211.00743 [gr-qc] .
  * Ghosh  _et al._ [2018] S. G. Ghosh, D. V. Singh, and S. D. Maharaj, Regular black holes in einstein-gauss-bonnet gravity, Phys. Rev. D 97, 104050 (2018).
  * Chatterjee [2021] B. Chatterjee, Black holes in Einstein-Gauss-Bonnet gravity: dynamical and 4-dimensional novel stationary black hole, (2021), arXiv:2110.13850 [gr-qc] .
  * Kumar  _et al._ [2022] A. Kumar, R. K. Walia, and S. G. Ghosh, Bardeen Black Holes in the Regularized 4D Einstein–Gauss–Bonnet Gravity, Universe 8, 232 (2022), arXiv:2003.13104 [gr-qc] .
  * Rodrigues  _et al._ [2016] M. E. Rodrigues, E. L. B. Junior, G. T. Marques, and V. T. Zanchin, Regular black holes in $f(R)$ gravity coupled to nonlinear electrodynamics, Phys. Rev. D 94, 024062 (2016), [Addendum: Phys.Rev.D 94, 049904 (2016)], arXiv:1511.00569 [gr-qc] .
  * Santos  _et al._ [2023] L. C. N. Santos, F. M. da Silva, C. E. Mota, I. P. Lobo, and V. B. Bezerra, Kiselev black holes in f(R, T) gravity, Gen. Rel. Grav. 55, 94 (2023), arXiv:2301.02534 [gr-qc] .
  * Hayward [2006] S. A. Hayward, Formation and evaporation of regular black holes, Phys. Rev. Lett. 96, 031103 (2006), arXiv:gr-qc/0506126 .
  * Fan and Wang [2016] Z.-Y. Fan and X. Wang, Construction of Regular Black Holes in General Relativity, Phys. Rev. D 94, 124027 (2016), arXiv:1610.02636 [gr-qc] .
  * Chen and Zhang [2013] B. Chen and J.-j. Zhang, Electromagnetic duality in dyonic reissner-nordström/conformal-field-theory correspondence, Phys. Rev. D 87, 081505 (2013).
  * Culetu [2013] H. Culetu, On a regular modified Schwarzschild spacetime, (2013), arXiv:1305.5964 [gr-qc] .
  * Bronnikov [2001] K. A. Bronnikov, Regular magnetic black holes and monopoles from nonlinear electrodynamics, Phys. Rev. D 63, 044005 (2001), arXiv:gr-qc/0006014 .
  * Lessa and Silva [2023] L. A. Lessa and J. E. G. Silva, Regular black holes in Einstein cubic gravity, (2023), arXiv:2305.18254 [gr-qc] .
  * Toshmatov  _et al._ [2018] B. Toshmatov, Z. Stuchlík, and B. Ahmedov, Comment on “Construction of regular black holes in general relativity”, Phys. Rev. D 98, 028501 (2018), arXiv:1807.09502 [gr-qc] .
  * Rajabi and Nozari [2021] F. Rajabi and K. Nozari, Energy condition in unimodular f(R, T) gravity, Eur. Phys. J. C 81, 247 (2021).
  * Panotopoulos and Rincón [2019] G. Panotopoulos and A. Rincón, Quasinormal modes of regular black holes with non linear-Electrodynamical sources, Eur. Phys. J. Plus 134, 300 (2019), arXiv:1904.10847 [gr-qc] .
  * Fernando and Correa [2012] S. Fernando and J. Correa, Quasinormal Modes of Bardeen Black Hole: Scalar Perturbations, Phys. Rev. D 86, 064039 (2012), arXiv:1208.5442 [gr-qc] .
  * Konoplya [2003] R. A. Konoplya, Quasinormal behavior of the d-dimensional Schwarzschild black hole and higher order WKB approach, Phys. Rev. D 68, 024018 (2003), arXiv:gr-qc/0303052 .
  * Iyer and Will [1987] S. Iyer and C. M. Will, Black-hole normal modes: A wkb approach. i. foundations and application of a higher-order wkb analysis of potential-barrier scattering, Phys. Rev. D 35, 3621 (1987).
  * Matyjasek and Opala [2017] J. Matyjasek and M. Opala, Quasinormal modes of black holes. The improved semianalytic approach, Phys. Rev. D 96, 024011 (2017), arXiv:1704.00361 [gr-qc] .
  * Matyjasek and Telecka [2019] J. Matyjasek and M. Telecka, Quasinormal modes of black holes. II. Padé summation of the higher-order WKB terms, Phys. Rev. D 100, 124006 (2019), arXiv:1908.09389 [gr-qc] .
  * Konoplya  _et al._ [2019] R. A. Konoplya, A. Zhidenko, and A. F. Zinhailo, Higher order WKB formula for quasinormal modes and grey-body factors: recipes for quick and accurate calculations, Class. Quant. Grav. 36, 155002 (2019), arXiv:1904.10333 [gr-qc] .
  * Cardoso  _et al._ [2004] V. Cardoso, J. P. S. Lemos, and S. Yoshida, Quasinormal modes of Schwarzschild black holes in four-dimensions and higher dimensions, Phys. Rev. D 69, 044004 (2004), arXiv:gr-qc/0309112 .
  * Wahlang  _et al._ [2017] W. Wahlang, P. A. Jeena, and S. Chakrabarti, Quasinormal modes of scalar and Dirac perturbations of Bardeen de-Sitter black holes, Int. J. Mod. Phys. D 26, 1750160 (2017), arXiv:1703.04286 [gr-qc] .
  * Liu  _et al._ [2021] P. Liu, C. Niu, and C.-Y. Zhang, Linear instability of charged massless scalar perturbation in regularized 4D charged Einstein-Gauss-Bonnet anti de-Sitter black holes, Chin. Phys. C 45, 025111 (2021), arXiv:2005.01507 [gr-qc] .
  * Ohashi and Sakagami [2004] A. Ohashi and M.-a. Sakagami, Massive quasi-normal mode, Class. Quant. Grav. 21, 3973 (2004), arXiv:gr-qc/0407009 .
  * Konoplya and Zhidenko [2005] R. A. Konoplya and A. V. Zhidenko, Decay of massive scalar field in a Schwarzschild background, Phys. Lett. B 609, 377 (2005), arXiv:gr-qc/0411059 .
  * Mashhoon [1985a] B. Mashhoon, Stability of charged rotating black holes in the eikonal approximation, Phys. Rev. D 31, 290 (1985a).
  * Mashhoon [1985b] B. Mashhoon, Erratum: Stability of charged rotating black holes in the eikonal approximation, Phys. Rev. D 31, 2697 (1985b).
  * Ferrari and Mashhoon [1984] V. Ferrari and B. Mashhoon, New approach to the quasinormal modes of a black hole, Phys. Rev. D 30, 295 (1984).
  * Ponglertsakul  _et al._ [2018] S. Ponglertsakul, P. Burikham, and L. Tannukij, Quasinormal modes of black strings in de Rham–Gabadadze–Tolley massive gravity, Eur. Phys. J. C 78, 584 (2018), arXiv:1803.09078 [gr-qc] .
  * Cardoso  _et al._ [2009] V. Cardoso, A. S. Miranda, E. Berti, H. Witek, and V. T. Zanchin, Geodesic stability, Lyapunov exponents and quasinormal modes, Phys. Rev. D 79, 064016 (2009), arXiv:0812.1806 [hep-th] .
  * Hashimoto and Tanahashi [2017] K. Hashimoto and N. Tanahashi, Universality in Chaos of Particle Motion near Black Hole Horizon, Phys. Rev. D 95, 024007 (2017), arXiv:1610.06070 [hep-th] .
  * Park and Gwak [2023] J. Park and B. Gwak, Bound on Lyapunov exponent in Kerr-Newman-de Sitter black holes by a charged particle, (2023), arXiv:2312.13075 [gr-qc] .