Dark matter signatures of black holes with Yukawa potential

A. A. Araújo Filho dilto@fisica.ufc.br Kimet Jusufi kimet.jusufi@unite.edu.mk B. Cuadros–Melgar bertha@usp.br Genly Leon genly.leon@ucn.cl

Abstract

This study uses a nonsingular Yukawa–modified potential to obtain a static and spherically symmetric black hole solution with a cosmological constant. Such Yukawa–like corrections are encoded in two parameters, $\alpha$ and $\lambda$ , that modify Newton’s law of gravity in large distances, and a deformation parameter $\ell_{0}$ , which plays an essential role in short distances. The most significant effect is encoded in $\alpha$ , which modifies the total black hole mass with an extra mass proportional to $\alpha M$ , mimicking the dark matter effects at large distances from the black hole. On the other hand, the effect due to $\lambda$ is small for astrophysical values. We scrutinize the quasinormal frequencies and shadows associated with a spherically symmetric black hole and the thermodynamical behavior influenced by the Yukawa potential. In particular, the thermodynamics of this black hole displays a rich behavior, including possible phase transitions. We use the WKB method to probe the quasinormal modes of massless scalar, electromagnetic, and gravitational field perturbations. In order to check the influence of the parameters on the shadow radius, we consider astrophysical data to determine their values, incorporating information on an optically thin radiating and infalling gas surrounding a black hole to model the black hole shadow image. In particular, we consider Sgr A* black hole as an example and we find that its shadow radius changes by order of $10^{-9}$ , meaning that the shadow radius of a black hole with Yukawa potential practically gives rise to the same result encountered in the Schwarzschild black hole. Also, in the eikonal regime, using astrophysical data for Yukawa parameters, we show that the value of the real part of the QNMs frequencies changes by $10^{-18}$ . Such Yukawa–like corrections are, therefore, difficult to measure by observations of gravitational waves using the current technology.

keywords:

Quantum–corrected Yukawa–like gravitational potential , Dark matter , quasinormal frequencies , Black Holes shadows

^†^†journal: Physics of the Dark Universe

\affiliation

[1]organization=Departamento de Física Teórica and IFIC, Centro Mixto Universidad de Valencia–CSIC. Universidad de Valencia, addressline=Burjassot–46100, Valencia, Spain \affiliation[2]organization=Departamento de Física, Universidade Federal da Paraíba, addressline=Caixa Postal 5008, 58051-970, João Pessoa, Paraíba, Brazil \affiliation[3]organization=Physics Department, University of Tetova, addressline= Ilinden Street nn, 1200, Tetova, North Macedonia \affiliation[4]organization=Engineering School of Lorena - University of Sao Paulo (EEL-USP), addressline=Estrada Municipal do Campinho Nº 100, Campinho, CEP: 12602-810, Lorena, SP - Brazil \affiliation[5]organization=Departamento de Matemáticas, Universidad Católica del Norte, addressline= Angamos 0610, Casilla 1280 Antofagasta, Chile \affiliation[6]organization=Institute of Systems Science, Durban University of Technology, addressline=PO Box 1334, Durban 4000, South Africa

1 Introduction

In modern cosmology our understanding posits that the universe displays a fundamental property of homogeneity and isotropy on its largest scales. Additionally, cosmologists propose the existence of cold dark matter as a mysterious form of matter that remains invisible and interacts solely via the gravitational force González et al. (2023); Bond and Efstathiou (1984); Trimble (1987). Despite persistent efforts, direct detection of dark matter particles has proven elusive with its existence inferred fundamentally from its gravitational effects on galaxies and other expansive cosmic structures. On the other hand, another intriguing component known as dark energy is introduced to explain the observed accelerated expansion of the universe Carroll et al. (1992). This proposition garners strong support from a wealth of observational evidence Perlmutter et al. (1998); Riess et al. (1998); Perlmutter et al. (1999).

The $\Lambda$ CDM paradigm has emerged as the preeminent model in modern cosmology characterized by its remarkable ability to account for a wide array of cosmological observations with excellent economy in terms of its parameter set Aghanim et al. (2020). Yet, mysteries persist within fundamental physics notably concerning the elusive nature of dark matter and dark energy. In addition, the role of scalar fields in the physical description of the universe, particularly in the context of inflationary scenarios, has garnered considerable attention Guth (1981). Expanding upon the $\Lambda$ CDM framework, a quintessence scalar field has been incorporated in a generalized model to delve into these cosmic enigmas Ratra and Peebles (1988); Parsons and Barrow (1995); Rubano and Barrow (2001); Saridakis (2009); Cai et al. (2010); Wali Hossain et al. (2015); Barrow and Paliathanasis (2016). Additionally, multi–scalar field models have emerged as versatile tools capable of elucidating various cosmic epochs Elizalde et al. (2004, 2008); Skugoreva et al. (2015); Saridakis and Tsoukalas (2016); Paliathanasis (2019); Banerjee et al. (2023); Santos et al. (2023). Moreover, efforts have been made to forge a unified description encompassing matter–dominated and dark energy–driven epochs leading to the development of scalar–torsion theories with a Hamiltonian foundation Leon et al. (2022).

Nonetheless, it is essential to acknowledge that an alternative current of thought exists in the literature advocating for a departure from this conventional path. This perspective suggests that by modifying Einstein’s equations novel theories of gravity may provide compelling explanations for observed phenomena challenging the traditional framework Akrami et al. (2021); Leon and Saridakis (2009); De Felice and Tsujikawa (2010); Clifton et al. (2012); Capozziello and De Laurentis (2011); De Felice and Tsujikawa (2012); Xu et al. (2012); Bamba et al. (2012); Leon and Saridakis (2013); Kofinas et al. (2014); Bahamonde et al. (2015); Momeni and Myrzakulov (2015); Cai et al. (2016); Krssak et al. (2019); Dehghani et al. (2023). When addressing the enigma of dark matter, a vital component for explaining the curious flatness observed in galaxy rotation curves Salucci (2018), one of the initial theories posited to account for this phenomenon was the Modified Newtonian Dynamics (MOND) introduced by Milgrom Milgrom (1983). This theory brings about alterations to Newton’s classical law of gravitation Ferreira and Starkmann (2009); Milgrom and Sanders (2003); Tiret et al. (2007); Kroupa et al. (2010); Cardone et al. (2011); Richtler et al. (2011). In addition, there are other intriguing proposals in the realm of dark matter. These include concepts such as superfluid dark matter Berezhiani and Khoury (2015) and the intriguing notion of a Bose–Einstein condensate Boehmer and Harko (2007), among others.

On the other hand of the cosmic spectrum, black holes are captivating astronomical phenomena capable of testing gravitational theories under extreme conditions. Among the many fascinating aspects of black holes, one that particularly intrigues scientists is their shadowy silhouette. This outline emerges due to the immense gravitational pull exerted by the black hole distorting the paths of light rays that venture too close.

In essence, photons emitted from a brilliant source in the vicinity of a black hole face two distinct destinies: they may either succumb to the irresistible gravitational force spiralling inexorably toward the event horizon, or they may be deflected away embarking on an eternal journey into the vast cosmic expanse. Within this intricate dance critical geodesic trajectories delineate the boundary between these two outcomes known as unstable spherical orbits. By meticulously observing these pivotal paths of photons against the backdrop of the cosmos, we gain the remarkable ability to capture the elusive image of a black hole’s shadow Takahashi (2004); Hioki and Maeda (2009); Brito et al. (2015); Cunha et al. (2015); Ohgami and Sakai (2015); Moffat (2015); Abdujabbarov et al. (2016); Cunha and Herdeiro (2018); Mizuno et al. (2018); Tsukamoto (2018); Psaltis (2019); Amir et al. (2019); Gralla et al. (2019); Bambi et al. (2019); Cunha et al. (2019); Khodadi et al. (2020); Perlick and Tsupko (2022); Vagnozzi et al. (2023); Saurabh and Jusufi (2021); Jusufi et al. (2020); Tsupko et al. (2020). In recent breakthroughs the Event Horizon Telescope (EHT) collaboration achieved a monumental milestone by confirming the presence of black hole shadows for two supermassive black holes, M87 and Sgr A* Akiyama et al. (2019); Psaltis et al. (2020); Akiyama et al. (2021, 2022).

In this paper we investigate the phenomenology of Yukawa–like corrections on the black hole spacetime. The Yukawa effects are associated with the presence of a graviton mass leading to modifications in Newton’s law of gravity, particularly at significant distances. These phenomena have revealed compelling implications in cosmology establishing connections between dark matter, dark energy, and baryonic matter Jusufi et al. (2023); González et al. (2023). Theoretical and observational aspects of the Yukawa potential, dGRT gravity, and dark matter in cosmology have been studied extensively. Some examples of these studies can be found in these references, Loeb and Weiner (2011); Anderson et al. (2017); Kanzi et al. (2020); Cedeño et al. (2017); Berezhiani et al. (2009); Sakallı et al. (2018).

Our goal is to examine the effects of the Yukawa potential term in the newly discovered solution presented in this paper. The metric (10) is an original black hole solution due to the contribution of the deformation parameter $\ell_{0}$ . This solution extends the black hole solution presented in González et al. (2023) and includes short- and long-range modifications. However, we will focus more on discussing the solution’s phenomenological aspects for astrophysical black holes where the impact of $\ell_{0}$ is minimal. The accuracy of the approximate solution will be examined using procedures belonging to the realm of perturbation methods in differential equations. We will first check the metric’s thermodynamic stability by analyzing its properties such as temperature, entropy, and specific heat. Next, we will test the geometry’s dynamic stability by considering various perturbations and determining the corresponding quasinormal modes (QNM) of the system. Next, we will examine how the black hole geometry affects test particles by calculating the geodesic trajectories focusing on the lightlike case that produces the shadow of a black hole. Furthermore, we will investigate the connection between shadows and QNMs in the eikonal limit. As we proceed, we will use observational data to constrain the parameter $\alpha$ , which characterizes the Yukawa potential contribution, and compare our findings to the Schwarzschild solution.

The paper is organized as follows. We present the black hole solution with Yukawa potential in $\S$ 2, pointing out its correspondence with de Rham, Gabadadze, and Tolley (dRGT) gravity and discussing the energy conditions obeyed by the apparent fluid emerging from the model. In $\S$ 3 we explore the accuracy of the approximate solution using procedures belonging to the realm of perturbation methods in differential equations. In $\S$ 4 we study the thermodynamic properties of the system. In $\S$ 5 we study QNMs for scalar, electromagnetic, and gravitational field perturbations. Geodesic trajectories are treated in $\S$ 6. Furthermore, $\S$ 7 investigates the shadow image using an infalling gas model. In addition, $\S$ 8 elaborates on checking the eikonal QNMs, the correspondence with the shadow radius, and the numerical values for the eikonal QNMs. Finally, we discuss our results in $\S$ 9.

2 Black hole solution with Yukawa potential

We shall begin this section by pointing out the main physical motivation to include a Yukawa-modified potential in the context of black holes. Our first argument comes from cosmology. Namely, as was shown in Jusufi et al. (2023) and González et al. (2023), starting from the Yukawa potential,

\Phi(r)=-\frac{Mm}{r}\left(1+\alpha\,e^{-\frac{r}{\lambda}}\right).

(1)

One can effectively obtain the $\Lambda$ CDM model in cosmology as follows. There exists an expression for the dark matter [units $c=\hbar=G=1$ ] given by,

\Omega_{DM}=\frac{\sqrt{2\alpha\Omega_{B,0}}}{\lambda H_{0}\,(1+\alpha)}\,{(1+z)^{3}}\,,

(2)

here the subscript ‘0’ indicates quantities evaluated at the present time, namely $z=0$ . This implies that dark matter can be understood as a consequence of the modified Newton law quantified by $\alpha$ and $\Omega_{B}$ . If the contribution of baryonic matter disappears, i.e., $\Omega_{B,0}=0$ , automatically the effect of dark matter also vanishes, i.e., $\Omega_{DM}=0$ . In other words, dark matter can be viewed as an apparent effect. Furthermore, one can find an expression that relates baryonic matter, effective dark matter, and dark energy as,

\Omega_{DM}(z)=\sqrt{2\,\Omega_{B,0}\Omega_{\Lambda,0}}{(1+z)^{3}},

(3)

where we introduced the definition,

\Omega_{\Lambda,0}=\frac{1}{\lambda^{2}H^{2}_{0}}\frac{\alpha}{(1+\alpha)^{2}}.

(4)

Using the observational data that were reported in González et al. (2023), it was shown that one can get $\Omega_{\Lambda,0}\simeq 0.69$ and $\Omega_{DM}\simeq 0.26$ from cosmological observations. This of course is a very interesting result as one can effectively obtain the $\Lambda$ CDM parameters from Yukawa cosmology. A second motivation is to naturally extend the discussion of dark matter to the black hole geometry. It might be possible, as in the case of cosmology, to analyze the effect of a dark matter component for black holes using the same physics, namely through the Yukawa potential.

First of all, let us review the black hole solution with the effect of dark matter recently studied in González et al. (2023). The gravitational potential we considered is modified by the regular Yukawa–type potential

\Phi(r)=-\frac{Mm}{\sqrt{r^{2}+\ell_{0}^{2}}}\left(1+\alpha\,e^{-\frac{r}{\lambda}}\right).

(5)

Notice that the wavelength of massive graviton is represented by $\lambda=\frac{\hbar}{m_{g}c}$ and $\ell_{0}$ is a deformed parameter of Planck length size. Let us see how the spacetime geometry around the black hole is modified in this theory. The general solution in the case of a static, spherically symmetric source reads,

ds^{2}=-f(r)dt^{2}+\frac{dr^{2}}{f(r)}+r^{2}(d\theta^{2}+\sin^{2}\theta d\phi^{2}).

(6)

The energy density of the modified matter can be computed from $\rho(r)=\frac{1}{4\pi}\Delta\Phi(r)$ . After considering it, we obtain the following relation for the energy density,

\displaystyle\rho(r)=\frac{e^{-\frac{r}{\lambda}}M\alpha\mathcal{X}}{4\pi r\lambda^{2}(r^{2}+\ell_{0}^{2})^{5/2}}+\frac{3M\ell_{0}^{2}}{4\pi(r^{2}+\ell_{0}^{2})^{5/2}},

(7)

where we defined $\mathcal{X}=(2\lambda-r)\ell_{0}^{4}+(3\lambda^{2}r+2\lambda r^{2}-2r^{3})\ell_{0}^{2}-r^{5}$ . The energy density, therefore, consists of two terms: the first one is proportional to $\alpha$ and it is important in large distances; the second one is proportional to $\ell_{0}$ and plays an important role in short distances instead. In particular, if we neglect the long–range modification, i.e., as a special case of our results, only the second term remains consistent with Nicolini et al. (2019).

To find a black hole solution, let us expand only the first factor in Eq. (7) in a series around $\ell_{0}$ . After taking into account the leading order terms in $\alpha$ and $\ell_{0}$ , we get

\displaystyle\rho(r)=-\frac{M\alpha}{4\pi r\lambda^{2}}e^{-\frac{r}{\lambda}}+\frac{3M\ell_{0}^{2}}{4\pi(r^{2}+\ell_{0}^{2})^{5/2}}+\mathcal{O}(\ell_{0}^{2}\alpha).

(8)

The first term coincides with the result obtained in González et al. (2023). Moreover, it is worth mentioning that the negative sign reflects that the energy conditions are violated inside the black hole. On the other hand, we assume that the Einstein field equations with a cosmological constant still hold, i.e., $G_{\mu\nu}+\Lambda g_{\mu\nu}=8\pi T_{\mu\nu}$ , so that the effects of the effective dark matter are encoded in the stress–energy tensor. Then, from the gravitational field equations we obtain,

\frac{rf^{\prime}(r)+f(r)-1}{r^{2}}+\Lambda-\frac{2M\alpha}{r\lambda^{2}}e^{-\frac{r}{\lambda}}+\frac{6M\ell_{0}^{2}}{(r^{2}+\ell_{0}^{2})^{5/2}}=0,

(9)

yielding the solution

f(r)=1-\frac{2Mr^{2}}{(r^{2}+\ell_{0}^{2})^{3/2}}-\frac{2M\alpha(r+\lambda)e^{-\frac{r}{\lambda}}}{r\lambda}-\frac{\Lambda r^{2}}{3}+\frac{c_{1}}{r},

(10)

where $c_{1}$ is an integration constant. To our knowledge, this is a new solution and generalizes that one encountered in Ref. González et al. (2023). The second term modifies the geometry due to $\ell_{0}$ and plays an important role in short–range distances to cure the black hole singularity and if we neglect the third, fourth, and last terms, our result matches with Nicolini et al. (2019). The third term is due to the apparent dark matter effect, while the fourth term is the contribution due to the cosmological constant. Since $\ell_{0}$ is of Planck length order, i.e., $\ell_{0}\sim 10^{-35}\,m$ Nicolini et al. (2019), in large distances and astrophysical black holes with large mass $M$ , we must have $r\gg\ell_{0}$ , then, we can neglect the effect of $\ell_{0}$ . Setting $c_{1}=0$ the solution turns out to be

\displaystyle f(r)=1-\frac{2M}{r}-\frac{2M\alpha(r+\lambda)e^{-\frac{r}{\lambda}}}{r\lambda}-\frac{\Lambda r^{2}}{3},

(11)

which is in agreement with González et al. (2023). However, by setting $\ell_{0}=0$ directly in (9) and integrating it, we obtain

\displaystyle f(r)=1+\frac{c_{2}}{r}-\frac{2M\alpha(r+\lambda)e^{-\frac{r}{\lambda}}}{r\lambda}-\frac{\Lambda r^{2}}{3}.

(12)

For consistency, we must recover (11) by using the freedom to choose $c_{2}$ to be $c_{2}=-2M$ . Therefore, we have to take $\ell_{0}=0$ and integrate the results in the same solution. If we desire to study the properties of small black holes, i.e., black holes with small mass $M$ , we can neglect the effect of $\alpha$ . In the present work we will focus on astrophysical black holes instead. On the other hand, $\lambda$ is of Mpc order and it is important on large distances, thus, we can perform a series expansion around $\varepsilon=1/\lambda$ , yielding the result González et al. (2023),

\displaystyle f(r)=1-\frac{2M(1+\alpha)}{r}+\frac{M\alpha r}{\lambda^{2}}-\frac{\Lambda r^{2}}{3}.

(13)

It is, therefore, natural to define the true or the physical mass of the black hole to be $\mathcal{M}=M(1+\alpha)$ and write the solutions in terms of the physical mass; then, the exact solution is written as

f(r)=1-\frac{2\mathcal{M}r^{2}}{(1+\alpha)(r^{2}+\ell_{0}^{2})^{3/2}}-\frac{2\mathcal{M}\alpha(r+\lambda)e^{-\frac{r}{\lambda}}}{r(1+\alpha)\lambda}-\frac{\Lambda r^{2}}{3}.

(14)

and the approximate solution becomes González et al. (2023),

\displaystyle f(r)=1-\frac{2\mathcal{M}}{r}+\frac{\mathcal{M}\alpha r}{(1+\alpha)\lambda^{2}}-\frac{\Lambda r^{2}}{3}.

(15)

Refer to caption — Figure 1: The metric functions in terms of $r$ for $\mathcal{M}=1$ , $\lambda=10^{5}$ , $\Lambda=10^{-5}$ , $\ell_{0}=0$ , and $\alpha=1$ .

The natural interpretation is that the extra mass term $M_{\rm dark\,matter}=M\alpha$ can be linked to the dark matter effect, which is only an apparent effect in our model. We point out that the interesting reference Laurentis and Salucci (2022) argued that in order to obtain the accurate mass distribution of M87, the mass of the M87 black hole should grow due to dark matter. They of course used observational data and empirical models for dark matter, but this result is consistent with our paper, namely, we quite naturally obtain a shift in the mass of the system as we pointed out in the last equation and the interpretation that $\alpha$ encodes dark matter effect makes sense. The fact that the mass of the system is increased in our model yields an interesting consequence that fits a well known argument for the existence of dark matter according to which the deflection of light is increased by a point mass (which can be a black hole or a galaxy). The deflection angle yields,

{\rm Deflection\,angle}=4\mathcal{M}/b=4M(1+\alpha)/b.

This relation, in principle, should also hold for galaxies and again we get an increase of the deflection angle in terms of an apparent dark matter.

In Fig. 1 we have plotted the metric functions $f(r)$ given by Eq. (15) (approximated solution) and Eq. (14) (complete solution). Further, if we neglect the term $\mathcal{O}(\alpha/\lambda^{2})$ in (15), we obtain the Kottler spacetime, i.e., a Schwarzschild black hole with a cosmological constant,

\displaystyle f(r)\simeq 1-\frac{2\mathcal{M}}{r}-\frac{\Lambda r^{2}}{3}.

(16)

The expression above displayed in Eq. (15) yields three distinct solutions representing an equal number of horizons. However, the most physically relevant to us is the event horizon given by,

r_{+}=\frac{1}{\lambda^{2}\Lambda}\left[\alpha\mathcal{M}-\frac{\sqrt[3]{2}\left(\lambda^{4}\Lambda+\alpha^{2}\mathcal{M}^{2}\right)}{\sqrt[3]{\mathcal{A}}}-\frac{\sqrt[3]{\mathcal{A}}}{\sqrt[3]{2}}\right],

(17)

where

\displaystyle\mathcal{A}

\displaystyle=

\displaystyle-2\alpha^{3}\mathcal{M}^{3}+\sqrt{\mathcal{B}}+6(\alpha+1)\lambda^{6}\Lambda^{2}\mathcal{M}-3\alpha\lambda^{4}\Lambda\mathcal{M},

and

	$\displaystyle\mathcal{B}$	$\displaystyle=$	$\displaystyle\left(2\alpha^{3}\mathcal{M}^{3}+3\lambda^{4}\Lambda\mathcal{M}\left(\alpha-2(\alpha+1)\lambda^{2}\Lambda\right)\right)^{2}$		(18)
		$\displaystyle-$	$\displaystyle 4\left(\lambda^{4}\Lambda+\alpha^{2}\mathcal{M}^{2}\right)^{3}.$		(18)

For the metric (6) the Kretschmann scalar $K=R_{abcd}R^{abcd}$ is given by,

\displaystyle K=\frac{G^{2}}{c^{4}}\left(f^{\prime\prime}(r)^{2}+\frac{4f^{\prime}(r)^{2}}{r^{2}}+\frac{4(f(r)-1)^{2}}{r^{4}}\right).

(19)

Replacing the exact solution (14) and using units where $\frac{G}{c^{2}}=1$ , we obtain,

	$\displaystyle K\|_{f(r)=\eqref{exact-sol}}$
	$\displaystyle=\frac{4}{r^{2}}\Bigg{[}\frac{4\ell_{0}^{2}\mathcal{M}r}{(1+\alpha)\left(\ell_{0}^{2}+r^{2}\right)^{5/2}}$
	$\displaystyle-\frac{2\mathcal{M}}{(1+\alpha)r^{2}}\left(\frac{r^{5}}{\left(\ell_{0}^{2}+r^{2}\right)^{5/2}}+\frac{\alpha e^{-\frac{r}{\lambda}}\left(\lambda^{2}+r^{2}+\lambda r\right)}{\lambda^{2}}\right)+\frac{2\Lambda r}{3}\Bigg{]}^{2}$
	$\displaystyle+\frac{4}{r^{4}}\left[\frac{2\mathcal{M}}{r(1+\alpha)}\left(\frac{r^{3}}{\left(\ell_{0}^{2}+r^{2}\right)^{3/2}}+\frac{\alpha e^{-\frac{r}{\lambda}}(\lambda+r)}{\lambda}\right)+\frac{\Lambda r^{2}}{3}\right]^{2}$
	$\displaystyle+\Bigg{[}\frac{2\Lambda}{3}+\frac{2\mathcal{M}}{(1+\alpha)}\Bigg{(}-\frac{11\ell_{0}^{2}r^{2}}{\left(\ell_{0}^{2}+r^{2}\right)^{7/2}}+\frac{2r^{4}}{\left(\ell_{0}^{2}+r^{2}\right)^{7/2}}$
	$\displaystyle+\frac{2\ell_{0}^{4}}{\left(\ell_{0}^{2}+r^{2}\right)^{7/2}}+\frac{\alpha e^{-\frac{r}{\lambda}}(\lambda+r)\left(2\lambda^{2}+r^{2}\right)}{\lambda^{3}r^{3}}\Bigg{)}\Bigg{]}^{2}.$		(20)

However, we use the approximated solution (15) such that,

\displaystyle K|_{f(r)=\eqref{approximated-solution}}=\frac{8\Lambda^{2}}{3}+\frac{8\mathcal{M}^{2}\left(\frac{\alpha^{2}r^{4}}{(\alpha+1)^{2}\lambda^{4}}+6\right)}{r^{6}}+\frac{8\alpha\Lambda\mathcal{M}}{(\alpha+1)\lambda^{2}r}.

(21)

Evaluating in $r_{+}$ we acquire a nonzero value for $K$ as follows,

	$\displaystyle K_{+}=\frac{8}{3}\Lambda^{2}\Bigg{[}\frac{3\alpha\mathcal{M}}{(\alpha+1)\left(-\sqrt[3]{\frac{2}{\mathcal{A}}}\left(\lambda^{4}\Lambda+\alpha^{2}\mathcal{M}^{2}\right)-\sqrt[3]{\frac{\mathcal{A}}{2}}+\alpha\mathcal{M}\right)}$
	$\displaystyle+\frac{3\lambda^{12}\Lambda^{4}\mathcal{M}^{2}\left(\frac{\alpha^{2}\left(\sqrt[3]{\frac{2}{\mathcal{A}}}\left(\lambda^{4}\Lambda+\alpha^{2}\mathcal{M}^{2}\right)+\sqrt[3]{\frac{A}{2}}-\alpha\mathcal{M}\right)^{4}}{(\alpha+1)^{2}\lambda^{12}\Lambda^{4}}+6\right)}{\left(\sqrt[3]{\frac{2}{\mathcal{A}}}\left(\lambda^{4}\Lambda+\alpha^{2}\mathcal{M}^{2}\right)+\sqrt[3]{\frac{\mathcal{A}}{2}}-\alpha\mathcal{M}\right)^{6}}+1\Bigg{]}.$		(22)

As $r$ approaches to zero, we have

K\simeq\frac{48\mathcal{M}^{2}}{r^{6}}+\mathcal{O}\left(\frac{1}{r}\right),

(23)

reproducing the Schwarzschild’s value at $r=0$ .

2.1 Correspondence with the black hole solution in dRGT gravity

Up to now we have seen how the black hole solution could be obtained by assuming the form of Einstein field equations with a cosmological constant, where the stress–energy was further modified due to Yukawa potential. From a field equation point of view – or better to say from the level of action – one can obtain a similar result. As we shall see, we shall point out an interesting link between our black hole solution and the spherical black hole solutions within the framework of dRGT massive gravity, elucidated in the work of de Rham et al. de Rham et al. (2011). In this theoretical framework it was shown that the Einstein field equations must satisfy de Rham et al. (2011),

\displaystyle G^{\mu}_{\nu}+m^{2}_{g}X^{\mu}_{\nu}=8\pi GT^{\mu(m)}_{\nu}.

(24)

Here the stress–energy tensor $T^{\mu(m)}_{\nu}$ is computed by using the matter Lagrangian. This equation is derived from the action of dRGT massive gravity, which has the form,

S=\frac{1}{16\pi G}\int d^{4}x\sqrt{-g}\left[R+m_{g}^{2}\,\mathcal{U}(g,f)\right]+S_{m}\,.

(25)

In this context $R$ represents the Ricci scalar, $m_{g}$ denotes the graviton mass, and $\mathcal{U}$ stands for the self–interacting potential of graviton. It is essential to structure $U(g,f)$ in a specific form to prevent the emergence of the Boulware–Deser ghost de Rham et al. (2011)

	$\displaystyle\mathcal{U}\equiv\mathcal{U}_{2}+\alpha_{3}\mathcal{U}_{3}+\alpha_{4}\mathcal{U}_{4}\,,$
	$\displaystyle\mathcal{U}_{2}\equiv[\mathcal{K}]^{2}-[\mathcal{K}^{2}]\,,$
	$\displaystyle\mathcal{U}_{3}\equiv[\mathcal{K}]^{3}-3[\mathcal{K}][\mathcal{K}^{2}]+2[\mathcal{K}^{3}]\,,$
	$\displaystyle\mathcal{U}_{4}\equiv[\mathcal{K}]^{4}-6[\mathcal{K}]^{2}[\mathcal{K}^{2}]+3[\mathcal{K}^{2}]^{2}+8[\mathcal{K}][\mathcal{K}^{3}]-6[\mathcal{K}^{4}]\,,$

where the tensor $\mathcal{K}^{\mu}_{\nu}$ is defined as de Rham et al. (2011); Ghosh et al. (2016),

\mathcal{K}^{\mu}_{\nu}\equiv\delta^{\mu}_{\nu}-\sqrt{g^{\mu\lambda}\partial_{\lambda}\varphi^{a}\partial_{\nu}\varphi^{b}f_{ab}},

(26)

in which

\displaystyle[\mathcal{K}]=\mathcal{K}^{\mu}_{\mu}\,\,\,\,\text{and}\,\,\,(\mathcal{K}^{i})^{\mu}_{\nu}=\mathcal{K}^{\mu}_{\rho_{1}}\mathcal{K}^{\rho_{1}}_{\rho_{2}}\ldots\mathcal{K}^{\rho_{i}}_{\nu},

(27)

the physical metric is denoted by $g_{\mu\nu}$ while $f_{\mu\nu}$ serves as a reference (fiducial) metric and $\varphi^{a}$ represents the Stückelberg fields. In this study we adopt the unitary gauge, where $\varphi^{a}=x^{\mu}\delta^{a}_{\mu}$ , thus de Rham et al. (2011),

\sqrt{g^{\mu\lambda}\partial_{\lambda}\varphi^{a}\partial_{\nu}\varphi^{b}f_{ab}}=\sqrt{g^{\mu\lambda}f_{\lambda\nu}}.

(28)

The massive graviton tensor $X^{\mu}_{\nu}$ is given by de Rham et al. (2011),

$\displaystyle X^{\mu}_{\nu}$	$\displaystyle=\mathcal{K}^{\mu}_{\nu}-[\mathcal{K}]\delta^{\mu}_{\nu}-\bar{\alpha}\left[(\mathcal{K}^{2})^{\mu}_{\nu}-[\mathcal{K}]\mathcal{K}^{\mu}_{\nu}+\frac{1}{2}\delta^{\mu}_{\nu}([\mathcal{K}]^{2}-[\mathcal{K}^{2}])\right]$
	$\displaystyle~{}~{}+3\beta\left[(\mathcal{K}^{3})^{\mu}_{\nu}-[\mathcal{K}](\mathcal{K}^{2})^{\mu}_{\nu}+\frac{1}{2}\mathcal{K}^{\mu}_{\nu}([\mathcal{K}]^{2}-[\mathcal{K}^{2}])\right.$
	$\displaystyle~{}~{}\left.-\frac{1}{6}\delta^{\mu}_{\nu}([\mathcal{K}]^{3}-3[\mathcal{K}][\mathcal{K}^{2}]+2[\mathcal{K}^{3}])\right]\,,$	(29)

where $\alpha_{3}=\frac{\bar{\alpha}-1}{3}$ and $\alpha_{4}=\frac{\beta}{4}+\frac{1-\bar{\alpha}}{12}$ . The terms of order $\mathcal{O}(\mathcal{K}^{4})$ disappear under the fiducial metric ansatz $f_{\mu\nu}=\operatorname{diag}(0,0,C^{2},C^{2}\sin^{2}{\theta})$ , where $C$ is a positive constant. In this sense the static and spherically symmetric solution is given by de Rham et al. (2011); Ghosh et al. (2016),

\displaystyle f(r)=1-\frac{2\mathcal{M}}{r}-\frac{\Lambda r^{2}}{3}+\gamma\,r+\zeta,

(30)

for $\zeta=0$ and $\gamma=\mathcal{M}\alpha/\lambda^{2}$ along with $\Lambda\sim m_{g}^{2}$ .

2.2 Energy conditions

Let us consider now the modifications introduced by the Yukawa potential as an apparent type of matter described by a fluid. If we calculate the stress–energy tensor corresponding to this fluid using Eq. (15), we obtain,

	$\displaystyle T^{t}_{t}=T^{r}_{r}$	$\displaystyle=$	$\displaystyle\frac{{\cal M}\alpha}{4\pi(1+\alpha)\lambda^{2}r}\,,$		(31)
	$\displaystyle T^{\theta}_{\theta}=T^{\phi}_{\phi}$	$\displaystyle=$	$\displaystyle\frac{{\cal M}\alpha}{8\pi(1+\alpha)\lambda^{2}r}\,.$		(32)

As stated in Simpson and Visser (2019), a necessary aspect to accomplish the null energy condition (NEC) is that both $\rho+p_{||}\geq 0$ and $\rho+p_{\perp}\geq 0$ are valid. For the apparent fluid we are modelling we identify $\rho=-T^{t}_{t}$ , $p_{||}=T^{r}_{r}$ , and $p_{\perp}=T^{\theta}_{\theta}=T^{\phi}_{\phi}$ and using the expressions (31) and (32), we obtain

	$\displaystyle\rho+p_{\|\|}$	$\displaystyle=$	$\displaystyle 0\,,$		(33)
	$\displaystyle\rho+p_{\perp}$	$\displaystyle=$	$\displaystyle-\frac{{\cal M}\alpha}{8\pi(1+\alpha)\lambda^{2}r}\,,$		(34)

in which the latter relation clearly violates NEC as long as $\alpha>0$ . As mentioned in Visser (1995), this result is enough to assume that weak, strong, and dominant energy conditions are equally violated.

We can also check if the exact solution given in Eq. (14) can alleviate this violation. Calculating both conditions we have,

	$\displaystyle\rho+p_{\|\|}$	$\displaystyle=0\,,$		(35)
	$\displaystyle\rho+p_{\perp}$	$\displaystyle=-\frac{\mathcal{M}\alpha(r+\lambda)e^{-r/\lambda}}{8\pi\lambda^{3}(1+\alpha)r}+\frac{15\mathcal{M}\ell_{0}^{2}r^{2}}{8\pi(1+\alpha)(\ell_{0}^{2}+r^{2})^{7/2}}.$		(36)

We will mathematically show that it is possible to achieve a fine–tuning such that Eq. (36) becomes zero or positive. However, physically it is practically impossible due to the order of magnitude of $\ell_{0}$ , as we shall see.

Indeed, for small values of the radial coordinate, we can approximate Eq. (36) as

\rho+p_{\perp}\simeq-\frac{\alpha\mathcal{M}}{8\pi(1+\alpha)\lambda^{2}r}+\mathcal{O}\left(r\right).

(37)

Such that $\rho+p_{\perp}\geq 0$ only for $\alpha\leq 0$ at the interior of the black hole. Nevertheless, classical energy conditions can be fulfilled by this solution for large $r$ , since we have

\rho+p_{\perp}\simeq\frac{15\ell_{0}^{2}\mathcal{M}}{8\pi(1+\alpha)r^{5}}-\frac{\alpha\mathcal{M}e^{-\frac{r}{\lambda}}(\lambda+r)}{8\pi\lambda^{3}(1+\alpha)r}.

(38)

Such that the negative term can be exponentially suppressed and $\rho+p_{\perp}\geq 0$ . The required condition is

15\ell_{0}^{2}\geq\frac{\alpha r^{4}e^{-\frac{r}{\lambda}}(\lambda+r)}{\lambda^{3}}.

(39)

Therefore, the strategy is to show that the maximum of

h(r)=\frac{\alpha r^{4}e^{-\frac{r}{\lambda}}(\lambda+r)}{\lambda^{3}}

(40)

is less or equal than $15\ell_{0}^{2}$ . The extreme values of $h(r)$ satisfy the following relation,

-\frac{\alpha r^{3}e^{-\frac{r}{\lambda}}\left(-4\lambda^{2}+r^{2}-4\lambda r\right)}{\lambda^{4}}=0.

(41)

The possible physical solutions are $r=0$ , $r=2\left(1+\sqrt{2}\right)\lambda$ , or $\alpha=0$ . Because $h^{\prime\prime}(2\left(1+\sqrt{2}\right)\lambda)=-32\left(10+7\sqrt{2}\right)e^{-2\left(1+\sqrt{2}\right)}\alpha<0$ , the maximum is obtained at $r=2\left(1+\sqrt{2}\right)\lambda$ and the required condition turns to be

\lambda\neq 0,\alpha\leq\frac{15e^{2+2\sqrt{2}}\ell_{0}^{2}}{16\left(99+70\sqrt{2}\right)\lambda^{2}}\lesssim\frac{0.591938\ell_{0}^{2}}{\lambda^{2}}=a_{\star}.

(42)

As expected, $a_{\star}$ is practically zero. In this manner, classical energy conditions are not fulfilled by the apparent fluid emerging from the Yukawa black hole solution – unless that $\alpha<0$ . At this point we should emphasize that the stress–energy tensor that we obtained corresponds to an apparent effect (it mimics the dark matter effect due to the modification of the gravitational potential, i.e., it is a screening effect in large distances) and not to real matter. Alternatively, we can consider this system as made of a composition of several fluids as Refs. Visser (2020); Boonserm et al. (2020) do and remark that the resultant does not necessarily preserves NEC. Nevertheless, the violation of NEC by this composite – but not physical fluid representing the Yukawa modified black hole – does not rule out the study of its interesting properties, especially in cosmological contexts where $\alpha$ has been estimated around $0.4$ González et al. (2023). Moreover, as we will see in $\S$ 7, we can constrain $\alpha$ using EHT observations and it can also assume negative values preventing NEC violation.

In the next section we will discuss the accuracy of the approximate solution with which we will work throughout this paper.

3 Accuracy of the approximate solution

The approach described here involves identifying an order parameter $\varepsilon={\hbar}^{-1}{m_{g}c}$ to tackle a problem that depends on two different scales: (i) a “slow” variable that has a solution similar to the outer solution in a boundary layer problem, and (ii) a “fast” system that describes rapid evolution occurring over shorter radii, which can be solved by changing the scale of the system. This approach is similar to the inner solution of a boundary layer problem and involves seeking the solution of each subsystem by using regular perturbation expansion.

For singularly perturbed problems the subsystems will have simpler structures than the whole problem. Combining the results of these two limiting problems, we can obtain information about the dynamics for small values of $\varepsilon$ . These procedures belong to the realm of perturbation methods in differential equations Fenichel (1979); Fusco and Hale (1989); Berglund and Gentz (2006); Verhulst (2006); Holmes (2012); Kevorkian and Cole (2013); Awrejcewicz (2014) and there can be used to construct uniformly valid approximations of the solutions of perturbation problems using those ones that satisfy the original expressions in the limit of $\varepsilon\rightarrow 0$ .

For this analysis we use the black hole physical mass $\mathcal{M}$ and we properly write the corresponding solution. For simplicity we set $\ell_{0}=0$ . Hence, Eq. (9) is recast as,

\displaystyle P(\varepsilon)[f(r)]

\displaystyle:=\frac{rf^{\prime}(r)+f(r)-1}{r^{2}}+\Lambda-\frac{2\mathcal{M}\alpha}{(1+\alpha)r}\varepsilon^{2}e^{-\varepsilon r}=0,

(43)

where we introduce $\varepsilon=1/\lambda$ . Associated with this problem we have the zeroth–order problem,

\displaystyle P(0)[f(r)]

\displaystyle:=\frac{rf^{\prime}(r)+f(r)-1}{r^{2}}+\Lambda=0.

(44)

For convenience we redefine $f(r)=g(r)/r$ and pose the problems $\tilde{P}(\varepsilon)[g(t)]=r^{2}P(\varepsilon)[g(r)/r]$ and $\tilde{P}(0)[g(t)]=r^{2}P(0)[g(r)/r]$ , namely,

\displaystyle\tilde{P}(\varepsilon)[g(r)]

\displaystyle:=g^{\prime}(r)-1+\Lambda{r^{2}}-\frac{2\mathcal{M}\alpha}{(1+\alpha)}\varepsilon\left(\varepsilon r\right)e^{-\varepsilon r}=0,

(45)

and

\displaystyle\tilde{P}(0)[g(r)]

\displaystyle:=g^{\prime}(r)-1+\Lambda{r^{2}}=0.

(46)

For both problems $P(\varepsilon)[g(r)]=0$ and $P(0)[g(r)]=0$ we choose the initial layer condition as,

\lim_{r\rightarrow 0^{+}}g(r)=-2\mathcal{M}.

(47)

3.1 Inner solution

In our case Eq. (45) is written as,

\displaystyle g^{\prime}(r)=1-\Lambda{r^{2}}+\frac{2\mathcal{M}\alpha}{(1+\alpha)}r\varepsilon^{2}e^{-\varepsilon r},\;g(0)=-2\mathcal{M}.

(48)

By replacing

\displaystyle g_{\text{inner}}(r)=g_{0}(r)+\varepsilon g_{1}(r)+\varepsilon^{2}g_{2}(r)+\varepsilon^{3}g_{3}(r)+\varepsilon^{4}g_{4}(r)+\ldots

(49)

we have

	$\displaystyle g_{0}^{\prime}(r)+\varepsilon g_{1}^{\prime}(r)+\varepsilon^{2}g_{2}^{\prime}(r)+\varepsilon^{3}g_{3}^{\prime}(r)+\varepsilon^{4}g_{4}^{\prime}(r)+\ldots$
	$\displaystyle=1-\Lambda{r^{2}}+\frac{2\mathcal{M}\alpha}{(1+\alpha)}r\left[\varepsilon^{2}-r\varepsilon^{3}+\frac{r^{2}\varepsilon^{4}}{2}+\ldots\right],$		(50)
	$\displaystyle g_{0}(0)=-2\mathcal{M},\quad g_{j}(0)=0,\,\,\,\,j\geq 1.$		(51)

Now comparing equal powers of $\varepsilon$ , we obtain the following systems,

	$\displaystyle g_{0}^{\prime}(r)=1-\Lambda{r^{2}},\;g_{0}(0)=-2\mathcal{M}$
	$\displaystyle\implies g_{0}(r)=-2\mathcal{M}+r-\frac{\Lambda}{3}r^{3},$		(52)
	$\displaystyle g_{1}^{\prime}(r)=0,\;g_{1}(0)=0\implies g_{1}(r)=0,$		(53)
	$\displaystyle g_{2}^{\prime}(r)=\frac{2\mathcal{M}\alpha}{(1+\alpha)}r,\;g_{2}(0)=0\implies g_{2}(r)=\frac{\mathcal{M}\alpha}{(1+\alpha)}r^{2},$		(54)
	$\displaystyle g_{3}^{\prime}(r)=-\frac{2\mathcal{M}\alpha}{(1+\alpha)}r^{2},\;g_{3}(0)=0\implies g_{3}^{\prime}(r)=-\frac{2\mathcal{M}\alpha}{3(1+\alpha)}r^{2},$		(55)
	$\displaystyle g_{4}^{\prime}(r)=\frac{\mathcal{M}\alpha}{(1+\alpha)}r^{3},\;g_{4}(0)=0\implies g_{4}(r)=\frac{\mathcal{M}\alpha}{4(1+\alpha)}r^{4}.$		(56)

Plugging all together we recover the inner solution

	$\displaystyle g_{\text{inner}}(r)$	$\displaystyle=r-2\mathcal{M}-\frac{1}{3}\Lambda r^{3}-\left(\frac{\alpha\mathcal{M}r}{\alpha+1}\right)\varepsilon^{2}$
		$\displaystyle+\frac{2\alpha\mathcal{M}r^{3}}{3(\alpha+1)}\varepsilon^{3}-\frac{\left(\alpha\mathcal{M}r^{4}\right)}{4(\alpha+1)}\varepsilon^{4}.$		(57)

Then, using the relation $\varepsilon=1/\lambda$ and $f=g/r$ we get

	$\displaystyle f_{\text{inner}}(r)$	$\displaystyle=1-2\frac{\mathcal{M}}{r}-\frac{1}{3}\Lambda r^{2}-\left(\frac{\alpha\mathcal{M}r}{\alpha+1}\right)\left(\frac{1}{\lambda}\right)^{2}$
		$\displaystyle+\frac{2\alpha\mathcal{M}r^{2}}{3(\alpha+1)}\left(\frac{1}{\lambda}\right)^{3}-\frac{\left(\alpha\mathcal{M}r^{3}\right)}{4(\alpha+1)}\left(\frac{1}{\lambda}\right)^{4}.$		(58)

3.2 Outer solution

Defining

G(\tau)=g(\tau/\varepsilon),

(59)

from (45) we obtain the “slow system”,

\frac{dG(\tau)}{d\tau}=-\Lambda\epsilon^{-3}\tau^{2}+\epsilon^{-1}+\frac{2\mathcal{M}\alpha}{(1+\alpha)}\tau e^{-\tau},

(60)

attained after the rescaling $\tau=\varepsilon r$ , where we introduced the small parameter $\varepsilon=1/\lambda$ . Using the substitution

\displaystyle G(\tau)=\varepsilon^{-3}G_{-3}(\tau)+\varepsilon^{-2}G_{-2}(\tau)+\varepsilon^{-1}G_{-1}(\tau)+G_{0}(\tau),

(61)

we yield

	$\displaystyle\varepsilon^{-3}G_{-3}^{\prime}(\tau)+\varepsilon^{-2}G_{-2}^{\prime}(\tau)+\varepsilon^{-1}G_{-1}^{\prime}(\tau)+G_{0}^{\prime}(\tau)$
	$\displaystyle=-\Lambda\epsilon^{-3}\tau^{2}+\epsilon^{-1}+\frac{2\mathcal{M}\alpha}{(1+\alpha)}\tau e^{-\tau}.$		(62)

Now, comparing equal powers of $\varepsilon$ we have these systems,

	$\displaystyle G_{-3}^{\prime}(\tau)=-\Lambda{\tau^{2}}\implies G_{-3}(\tau)=-\frac{\Lambda\tau^{3}}{3}+c_{-3},$		(63)
	$\displaystyle G_{-2}^{\prime}(\tau)=0\implies G_{-1}=c_{-2},$		(64)
	$\displaystyle G_{-1}^{\prime}(\tau)=1\implies G_{-1}(\tau)=\tau+c_{-1},$		(65)
	$\displaystyle G_{0}^{\prime}(\tau)=\frac{2\mathcal{M}\alpha}{(1+\alpha)}\tau e^{-\tau}$
	$\displaystyle\implies G_{0}(\tau)=-\frac{2\alpha\mathcal{M}e^{-\tau}(\tau+1)}{\alpha+1}+c_{0}.$		(66)

Replacing all the terms we arrive at,

	$\displaystyle G(\tau)=\varepsilon^{-3}\left(-\frac{\Lambda\tau^{3}}{3}+c_{-3}\right)+\varepsilon^{-2}c_{-2}+\varepsilon^{-1}\left(\tau+c_{-1}\right)$
	$\displaystyle+\left(-\frac{2\alpha\mathcal{M}e^{-\tau}(\tau+1)}{\alpha+1}+c_{0}\right).$		(67)

We have the matching condition given by,

	$\displaystyle\lim_{r\rightarrow\infty}g_{\text{inner}}(r)=\lim_{\tau\rightarrow 0}G(\tau)$
	$\displaystyle=\varepsilon^{-3}c_{-3}+\varepsilon^{-2}c_{-2}+\varepsilon^{-1}c_{-1}+\left(c_{0}-\frac{2\alpha\mathcal{M}}{\alpha+1}\right).$		(68)

Comparing terms of the same power in $\varepsilon$ with (57), we require that

c_{-3}=c_{-2}=c_{-1}=0,\quad c_{0}=-2\mathcal{M}+\frac{2\alpha\mathcal{M}}{\alpha+1}.

(69)

Thus,

\displaystyle G(\tau)=-\varepsilon^{-3}\frac{\Lambda\tau^{3}}{3}+\varepsilon^{-1}\tau-2\mathcal{M}-\frac{2\alpha\mathcal{M}e^{-\tau}\tau}{\alpha+1}.

(70)

Replacing

\displaystyle g(r):=G(\varepsilon r)=r-2\mathcal{M}-\frac{\Lambda r^{3}}{3}-\frac{2\alpha\mathcal{M}e^{-\varepsilon r}\varepsilon r}{\alpha+1}

(71)

in $f(r)=g(r)/r$ and using $\varepsilon=1/\lambda$ we obtain,

\displaystyle f_{\text{outer}}(r)=1-\frac{2\mathcal{M}}{r}-\frac{\Lambda r^{2}}{3}-\frac{2\alpha\mathcal{M}e^{-\frac{r}{\lambda}}}{\lambda(\alpha+1)}.

(72)

3.3 Multiple scales - Derivative expansion method

The method of multiple scales determines solutions to perturbed oscillators by suppressing resonant forcing terms that would yield spurious secular terms in the asymptotic expansions. It assumes that the solution can be expressed as a function of multiple (just two for our purposes) radial variables, which are introduced to keep a well–ordered expansion, $x(r)=X(r,\tau)$ , where $r$ is the regular (or “fast”) radial variable and $\tau$ is a new variable describing the “slow–radial” dependence of the solution. The idea is to use any freedom that is in the $\tau$ -dependence to minimize the approximation’s error and, whenever possible, to remove secular terms Fenichel (1979); Fusco and Hale (1989); Berglund and Gentz (2006); Verhulst (2006); Holmes (2012); Kevorkian and Cole (2013); Awrejcewicz (2014).

Eq. (45) depends explicitly on two scales, $r$ and $\varepsilon r$ . In this way, for $\varepsilon\ll 1$ , we set the fast $\tau_{A}=r$ and slow $\tau_{B}=\varepsilon r$ scales. Thereby,

\frac{\mathrm{d}}{\mathrm{d}r}=\frac{\partial}{\partial\tau_{A}}+\varepsilon\frac{\partial}{\partial\tau_{B}}.

(73)

We denote the partial derivatives with respect to $\tau_{A}$ and $\tau_{B}$ with the indexes $A$ and $B$ , namely, $\frac{\partial f}{\partial\tau_{A}}=f_{,A}$ . Substituting the expression

g=g_{0}(r,\varepsilon r)+\varepsilon g_{1}(r,\varepsilon r)+\varepsilon^{2}g_{2}(r,\varepsilon r)+\ldots

(74)

in the Eq. (45) and identifying the dependencies on $\tau_{A}$ and $\tau_{B}$ , the differential operator is written as,

	$\displaystyle\tilde{P}(\varepsilon)[g_{0}+\varepsilon g_{1}+\varepsilon^{2}g_{2}+\ldots]:=$
	$\displaystyle\left[\frac{\partial}{\partial\tau_{A}}+\varepsilon\frac{\partial}{\partial\tau_{B}}\right]\left(g_{0}+\varepsilon g_{1}+\varepsilon^{2}g_{2}+\ldots\right)$
	$\displaystyle-1+\Lambda{\tau_{A}^{2}}-\frac{2\mathcal{M}\alpha}{(1+\alpha)}\varepsilon\tau_{B}e^{-\tau_{B}}$
	$\displaystyle=g_{0,A}-1+\Lambda{\tau_{A}^{2}}$
	$\displaystyle+\varepsilon\left(g_{1,A}+g_{0,B}-\frac{2\mathcal{M}\alpha}{(1+\alpha)}\tau_{B}e^{-\tau_{B}}\right)$
	$\displaystyle+\varepsilon^{2}\left(g_{2,A}+g_{1,B}\right)+\varepsilon^{3}g_{2,B}+\ldots\,\,\,.$		(75)

The initial condition is given by,

	$\displaystyle g(0)=g_{0}(0,0)+\varepsilon g_{1}(0,0)+\varepsilon^{2}g_{2}(0,0)+\ldots$	(76)
$\displaystyle\implies$	$\displaystyle-2\mathcal{M}+\varepsilon\cdot 0+\varepsilon^{2}\cdot 0+\varepsilon^{3}\cdot 0+\ldots$
	$\displaystyle=g_{0}(0,0)+\varepsilon g_{1}(0,0)+\varepsilon^{2}g_{2}(0,0)+\ldots\,\,\,.$	(77)

With it we have the following problems.

Order $\varepsilon^{0}$ :

\displaystyle g_{0,A}-1+\Lambda{\tau_{A}^{2}}=0,\;g_{0}(0,0)=-2\mathcal{M}.

(78)

with solution

g_{0}(\tau_{A},\tau_{B})=\tau_{A}+\frac{\Lambda\tau_{A}^{3}}{3}+c_{1}(\tau_{B}),\quad c_{1}(0)=-2\mathcal{M}.

(79)

Order $\varepsilon^{1}$ :

	$\displaystyle g_{1,A}+g_{0,B}-\frac{2\mathcal{M}\alpha}{(1+\alpha)}\tau_{B}e^{-\tau_{B}}=0,$
	$\displaystyle g_{0}(0,0)=-2\mathcal{M},\quad g_{1}(0,0)=0.$		(80)

In such a way that,

	$\displaystyle g_{1,A}+c_{1}^{\prime}(\tau_{B})-\frac{2\mathcal{M}\alpha}{(1+\alpha)}\tau_{B}e^{-\tau_{B}}=0,$
	$\displaystyle c_{1}(0)=-2\mathcal{M},\quad g_{1}(0,0)=0.$		(81)

Order $\varepsilon^{2}$ :

\displaystyle g_{2,A}+g_{1,B}=0,\quad g_{1}(0,0)=0,\quad g_{2}(0,0)=0.

(82)

We can eliminate higher–order terms (and possibly resonant terms) by choosing $g_{1}=g_{2}=\ldots=0$ , and solving

\displaystyle c_{1}^{\prime}(\tau_{B})-\frac{2\mathcal{M}\alpha}{(1+\alpha)}\tau_{B}e^{-\tau_{B}}=0,\quad c_{1}(0)=-2\mathcal{M}.

(83)

With

\displaystyle c_{1}(\tau_{B})=\frac{2\alpha\mathcal{M}e^{-\tau_{B}}\tau_{B}}{\alpha+1}+\frac{2\alpha\mathcal{M}e^{-\tau_{B}}}{\alpha+1}-\frac{2(2\alpha+1)\mathcal{M}}{\alpha+1}.

(84)

we have,

	$\displaystyle g(\tau_{A},\tau_{B})=-\frac{2(2\alpha+1)\mathcal{M}}{\alpha+1}+\tau_{A}+\frac{\Lambda\tau_{A}^{3}}{3}$
	$\displaystyle+\frac{2\alpha\mathcal{M}e^{-\tau_{B}}\tau_{B}}{\alpha+1}+\frac{2\alpha\mathcal{M}e^{-\tau_{B}}}{\alpha+1}.$		(85)

Finally,

	$\displaystyle rf(r)$	$\displaystyle=-\frac{2(2\alpha+1)\mathcal{M}}{\alpha+1}+r+\frac{\Lambda r^{3}}{3}$
		$\displaystyle+\frac{2\alpha\mathcal{M}r\varepsilon e^{-r\varepsilon}}{\alpha+1}+\frac{2\alpha\mathcal{M}e^{-r\varepsilon}}{\alpha+1}.$		(86)

In this case with multiple scales we provide the exact solution of the model. The regular expansion method, on the other hand, provides the solution,

\displaystyle rf(r)

\displaystyle=\frac{1}{3}\left(-6\mathcal{M}+\Lambda r^{3}+3r\right)-\frac{\alpha\mathcal{M}r^{2}\varepsilon^{2}}{\alpha+1}+\frac{2\alpha\mathcal{M}r^{3}\varepsilon^{3}}{3(\alpha+1)},

(87)

with the error term,

	$\displaystyle r(f(r)\|_{\eqref{eq.(47)}}-f(r)\|_{\eqref{eq.(48)}})=-\frac{2\alpha\mathcal{M}-2\alpha\mathcal{M}e^{-r\varepsilon}}{\alpha+1}$
	$\displaystyle+\frac{2\alpha\mathcal{M}r\varepsilon e^{-r\varepsilon}}{\alpha+1}+\frac{\alpha\mathcal{M}r^{2}\varepsilon^{2}}{\alpha+1}-\frac{2\alpha\mathcal{M}r^{3}\varepsilon^{3}}{3(\alpha+1)}.$		(88)

From Eq.(87) the expansion (57) is recovered. Finally, the main topic of discussion was centred around how the approximated solution for $f(r)$ , given by Eq. (15), accurately represents the exact solution (14) as the parameters $\varepsilon=1/\lambda$ and $\ell_{0}$ approach zero. This is particularly significant when both horizon solutions are close. The error is bounded by

\displaystyle\Big{|}f(r)|_{\eqref{exact-sol}}-f(r)|_{\eqref{approximated-solution}}\Big{|}\leq C_{3}\left|\frac{2\alpha\mathcal{M}r^{2}\varepsilon^{3}}{3(\alpha+1)}\right|,

(89)

where $C_{3}$ is a constant.

Once we have established that the approximation is reliable, we can further investigate the solution characteristics with $f(r)$ as defined by Eq. (15). We will begin by checking the thermodynamic stability of the solution through the study of its main properties in the next section.

4 Thermodynamic properties

4.1 Hawking temperature

The metric defined in Eq. (15) possesses a timelike Killing vector denoted as $\xi=\partial/\partial t$ . This metric gives rise to a conserved quantity inherently associated with $\xi$ . We can construct this conserved quantity using the Killing vector as follows,

\nabla^{\nu}(\xi^{\mu}\xi_{\mu})=-2\kappa\xi^{\nu}.

(90)

Here $\nabla_{\nu}$ represents the covariant derivative and $\kappa$ remains constant along the orbits defined by $\xi$ . This implies that the Lie derivative of $\kappa$ with respect to $\xi$ is identically zero,

\mathcal{L}_{\xi}\kappa=0.

(91)

Remarkably, $\kappa$ retains its constancy across the horizon and is called surface gravity. In the coordinate basis the timelike Killing vector components are expressed as $\xi^{\mu}=(1,0,0,0)$ . The surface gravity for the metric ansatz in Eq. (15) is given by,

\kappa={\left.\frac{f^{\prime}(r)}{2}\right|_{r={r_{+}}}}.

(92)

Hawking’s groundbreaking work, as presented in Refs. Hawking (1975a, b), revealed that black holes emit radiation. This one is associated with a temperature known as the Hawking temperature, expressed as $T=\kappa/2\pi$ Filho et al. (2023a); Sedaghatnia et al. (2023); Filho et al. (2023b, c). In our case, it can be expressed as follows

T=\frac{1}{4\pi}\left(\frac{\alpha\mathcal{M}}{(\alpha+1)\lambda^{2}}+\frac{2\mathcal{M}}{r_{+}^{2}}-\frac{2\Lambda}{3}r_{+}\right),

(93)

or in terms of the event horizon,

T=\frac{1}{4\pi{\cal C}}\left(\frac{\cal J}{r_{+}}-\frac{\Lambda r_{+}{\cal K}}{3}\right)\,,

(94)

with

$\displaystyle{\cal J}$	$\displaystyle=$	$\displaystyle 2(\alpha+1)\lambda^{2}+\alpha r_{+}^{2},$	(95)
$\displaystyle{\cal K}$	$\displaystyle=$	$\displaystyle 6(\alpha+1)\lambda^{2}-\alpha r_{+},$	(96)
$\displaystyle{\cal C}$	$\displaystyle=$	$\displaystyle 2(\alpha+1)\lambda^{2}-\alpha r_{+}^{2}\,.$	(97)

In Figs. 2 and 3, we plot the Hawking temperature as a function of mass $\mathcal{M}$ and event horizon $r_{+}$ , respectively, considering $\lambda=10^{5}$ , $\Lambda=10^{-5}$ , and $\alpha=1$ for both Schwarzschild and Yukawa–corrected black holes. From Fig. 3 we can see that for small event horizons Yukawa black hole temperature mimics Schwarzschild case. However, whereas Schwarzschild temperature goes asymptotically to zero, Yukawa black hole temperature reaches zero and turns to be negative at certain value of $r_{+}$ . According to the third law of thermodynamics, we can interpret this point as the value where the event horizon meets the cosmological horizon, i.e., the extreme black hole solution, after which the temperature loses its physical meaning.

4.2 Entropy

Certainly, when we consider $f(r_{+})=0$ , a precise expression for $\mathcal{M}$ can be obtained and it is given by,

\mathcal{M}=\frac{(\alpha+1)\lambda^{2}r_{+}\left(\Lambda r_{+}^{2}-3\right)}{3\alpha r_{+}^{2}-6(\alpha+1)\lambda^{2}}.

(98)

Furthermore, following the first law of thermodynamics, the Hawking temperature is defined as,

	$\displaystyle\tilde{T}_{f}=\frac{\mathrm{d}\mathcal{M}}{\mathrm{d}S}=\frac{1}{2\pi r_{+}}\frac{\mathrm{d}\mathcal{M}}{\mathrm{d}r_{+}}$
	$\displaystyle=\frac{1}{2\pi r_{+}}\left(\frac{(\alpha+1)\lambda^{2}\left(6(\alpha+1)\lambda^{2}+\alpha\Lambda r_{+}^{4}+\mathcal{Y}\right)}{3\left(\alpha r_{+}^{2}-2(\alpha+1)\lambda^{2}\right)^{2}}\right).$		(99)

where $\mathcal{Y}=3r_{+}^{2}\left(\alpha-2(\alpha+1)\lambda^{2}\Lambda\right)$ . This temperature contradicts the value presented in Eq. (93). In line with the approach detailed in Ref. Ma and Zhao (2014), under the assumption of the validity of the area law, we employ the corrected first law of thermodynamics to determine the temperature for regular black holes. The following first law determines the corrected temperature,

\Upsilon(r_{+},\mathcal{M},\alpha,\lambda,\Lambda)\mathrm{d}\mathcal{M}=\overset{\nsim}{T}_{f}\,\mathrm{d}S,

(100)

where $\overset{\nsim}{T}_{f}$ represents the corrected version of the Hawking temperature obtained through the first law of thermodynamics and $S$ denotes the entropy.

It is important to emphasize that the function $\Upsilon(r_{+},\mathcal{M},\alpha,\lambda,\Lambda)$ , dependent on the terms of the mass function, not only defines the first law for regular black holes in the specific scenario discussed here, but also plays a crucial role in establishing the first law for various other classes of regular black holes (see Maluf and Neves (2018)).

As indicated in Ref. Ma and Zhao (2014), the general formula for $\Upsilon(r_{+},\mathcal{M},\alpha,\lambda,\Lambda)$ is as follows,

\Upsilon(r_{+},\mathcal{M},\alpha,\lambda,\Lambda)=1+4\pi\int^{\infty}_{r_{+}}r^{2}\frac{\partial T^{0}_{0}}{\partial\mathcal{M}}\mathrm{d}r.

(101)

Here the notation $T^{0}_{0}$ corresponds to the stress–energy component associated with the energy density. Therefore, Eq. (101) can explicitly be evaluated as given below,

	$\displaystyle\Upsilon(r_{+},\mathcal{M},\alpha,\lambda,\Lambda)$
	$\displaystyle=\frac{3\left(\alpha r_{+}^{2}-2(\alpha+1)\lambda^{2}\right)^{2}\left(\frac{\alpha\mathcal{M}}{(\alpha+1)\lambda^{2}}+\frac{2\mathcal{M}}{r_{+}^{2}}-\frac{1}{3}2\Lambda r_{+}\right)}{4\pi(\alpha+1)\lambda^{2}\mathcal{Z}},$		(102)

where $\mathcal{Z}=6(\alpha+1)\lambda^{2}+\alpha\Lambda r_{+}^{4}+3r_{+}^{2}\left(\alpha-2(\alpha+1)\lambda^{2}\Lambda\right)$ . After these remarks all Hawking temperatures are in agreement with each other

\overset{\nsim}{T}_{f}=T=\Upsilon(r_{+},\mathcal{M},\alpha,\lambda,\Lambda)\tilde{T}_{f}.

(103)

In this manner, we may write the entropy $S$ as given below,

S=\int\frac{\Upsilon(r_{+},\mathcal{M},\alpha,\lambda,\Lambda)}{\overset{\nsim}{T}_{f}}\mathrm{d}\mathcal{M}=\pi r_{+}^{2}=\frac{A}{4}\,,

(104)

and we recover the usual area law.

4.3 Heat capacity

Indeed, when it comes to exploring additional thermodynamic quantities, the heat capacity is of particular interest. It can be expressed as,

C_{V}=T\frac{dS}{dT}=2\pi r_{+}^{2}\frac{\cal D}{\cal E}\,,

(105)

where

$\displaystyle{\cal D}$	$\displaystyle=$	$\displaystyle 12(1+\alpha)^{2}(\Lambda r_{+}^{2}-1)\lambda^{4}-8\Lambda\alpha r_{+}^{4}(1+\alpha)\lambda^{2}$	(106)
		$\displaystyle+\alpha^{2}r_{+}^{4}(\Lambda r_{+}^{2}+3)\,,$	(106)
$\displaystyle{\cal E}$	$\displaystyle=$	$\displaystyle 12(1+\alpha)^{2}(\Lambda r_{+}^{2}+1)\lambda^{4}-24\alpha r_{+}^{2}(1+\alpha)\lambda^{2}$	(107)
		$\displaystyle+\alpha^{2}r_{+}^{4}(\Lambda r_{+}^{2}-3)\,.$	(107)

In Figs. 4 and 5 we show the heat capacity as a function of the event horizon radius for small and large values of the event horizon, respectively. From Fig. 4 we notice that, for small black holes, the Yukawa corrected heat capacity follows Schwarzschild solution, however, as $r_{+}$ grows, it departs and displays a very rich behavior (see Fig. 5) which includes first order phase transitions at the points where $C_{V}=0$ . Moreover, in the regions where $C_{V}<0$ the solution is thermodynamically unstable.

Although our results show some instabilities and phase transitions in the Yukawa corrected black hole, they say nothing about the dynamic stability of the object. Thus, in the next section, in order to test this type of stability, we should perturb the background solution and obtain the corresponding quasinormal frequencies.

5 Quasinormal frequencies

During the ringdown phase the remarkable phenomenon of quasinormal modes arises showcasing unique oscillation patterns unaffected by initial perturbations. These modes reflect the system’s inherent characteristics stemming from the natural oscillations of spacetime, irrespective of specific initial conditions. Contrary to normal modes, which are related to closed systems, quasinormal modes correspond to open systems. Consequently, they gradually lose energy by emitting gravitational waves. From a mathematical perspective they are described as poles of the complex Green function.

Solutions to the wave equation within a system governed by the background metric $g_{\mu\nu}$ must be found to ascertain their frequencies. Nonetheless, obtaining analytical solutions for these modes often poses significant challenges. Various methodologies for getting solutions for these modes have been explored in scientific literature. The WKB (Wentzel–Kramers–Brillouin) method is among the most prevalent. This approach traces its roots to the pioneering contributions of Will and Iyer Iyer and Will (1987); Iyer (1987). Later enhancements extending up to sixth order were achieved by Konoplya Konoplya (2003) and up to thirteenth order by Matyjasek and Opala Matyjasek and Opala (2017).

5.1 Scalar perturbations

For our calculations we specifically examine perturbations via the scalar field incorporating the Klein–Gordon equation in the setting of curved spacetime,

\frac{1}{\sqrt{-g}}\partial_{\mu}(g^{\mu\nu}\sqrt{-g}\partial_{\nu}\Phi)=0.

(108)

Though the investigation of backreaction effects in this context is compelling, this manuscript steers its attention elsewhere. We concentrate primarily on examining the scalar field as a minor perturbation. The existing spherical symmetry in the scenario enables us to decompose the scalar field in a particular way detailed further below,

\Phi(t,r,\theta,\varphi)=\sum^{\infty}_{l=0}\sum^{l}_{m=-l}r^{-1}\Psi_{lm}(t,r)Y_{lm}(\theta,\varphi),

(109)

where the spherical harmonics are denoted by $Y_{lm}(\theta,\varphi)$ . By inserting the decomposition of the scalar field, as illustrated in Eq. (109) into Eq. (108), the equation takes a Schrödinger–like form. This transformation imbues the equation with wave-like attributes rendering it especially apt for our analysis,

-\frac{\partial^{2}\Psi}{\partial t^{2}}+\frac{\partial^{2}\Psi}{\partial r^{*2}}+V_{eff}(r^{*})\Psi=0.

(110)

The potential $V_{eff}$ is commonly referred to as the Regge–Wheeler potential or the effective potential. It encapsulates essential details about the black hole’s geometry. We also introduce the tortoise coordinate $r^{*}$ , which covers the full expanse of spacetime stretching as $r^{*}\rightarrow\pm\infty$ . This is represented by $\mathrm{d}r^{*}=\sqrt{[1/f(r)^{2}]}\mathrm{d}r$ . With some algebraic rearrangements we can express the effective potential as

V_{eff}=f(r)\left(\frac{l(l+1)}{r^{2}}+\frac{\frac{\alpha\mathcal{M}}{\lambda^{2}(1+\alpha)}+\frac{2\mathcal{M}}{r^{2}}-\frac{2\Lambda r}{3}}{r}\right).

(111)

To better comprehend Eq. (111), we provide Fig. 6 for $\mathcal{M}=\alpha=1$ , $\lambda=10^{5}$ , and $\Lambda=10^{-5}$ . We observe that a barrier-like shape forms when considering positive values for the cosmological constant $\Lambda$ , $\lambda$ , and $\alpha$ . Notice that, as $l$ increases, the height of $V_{eff}$ also rises.

In order to calculate the quasinormal modes, we focus on the WKB method. Next, our main goal is to obtain stationary solutions for the system. To achieve this we propose that the function $\Psi(t,r)$ can be represented as $\Psi(t,r)=e^{-i\omega t}\psi(r)$ , where $\omega$ signifies the frequency. By adopting this representation, we can seamlessly isolate the time-independent aspect of Eq. (110) as outlined below,

\frac{\partial^{2}\psi}{\partial r^{*2}}-\left[\omega^{2}-V_{eff}(r^{*})\right]\psi=0.

(112)

To effectively solve Eq. (112) it is essential to consider the relevant boundary conditions meticulously. For our particular scenario solutions meeting these conditions are distinguished by their purely ingoing behaviour near the horizon,

\psi^{\text{in}}(r^{*})\sim\begin{cases}C_{l}(\omega)e^{-i\omega r^{*}}&(r^{*}\rightarrow-\infty)\\ A^{(-)}_{l}(\omega)e^{-i\omega r^{*}}+A^{(+)}_{l}(\omega)e^{+i\omega r^{*}}&(r^{*}\rightarrow+\infty).\end{cases}

The complex constants $A^{(+)}_{l}(\omega)$ , $C_{l}(\omega)$ , and $A^{(-)}_{l}(\omega)$ are crucial for our following analysis. They are foundational in examining the quasinormal modes of a black hole, characterized by frequencies $\omega_{nl}$ that meet the condition $A^{(-)}_{l}(\omega_{nl})=0$ . These modes showcase a distinctive behaviour, manifesting as purely outgoing waves at spatial infinity and exclusively ingoing waves near the event horizon. The integers $n$ and $l$ delineate the overtone and multipole numbers, respectively. Furthermore, it is worth mentioning that the quasinormal modes spectrum is anchored on the eigenvalues of Eq. (112). To delve into these frequencies we utilize the WKB method, a semi-analytical approach reminiscent of quantum mechanics.

Furthermore, the WKB approximation, initially introduced by Schutz and Will Schutz and Will (1985), has emerged as an indispensable method for determining quasinormal modes, especially when studying particle scattering around black holes. This technique has been refined over the years with significant contributions by Konoplya Konoplya (2003, 2004). It is crucial, however, to recognize that the applicability of this method is contingent on the potential assuming a barrier–like form and levelling off to constant values as $r^{*}\to\pm\infty$ . By aligning the solution power series with the peak potential turning points, the quasinormal modes can be accurately derived Santos et al. (2016). Given these premises the sixth order WKB formula is expressed as,

\frac{i(\omega^{2}_{n}-V_{0})}{\sqrt{-2V^{\prime\prime}_{0}}}-\sum^{6}_{j=2}\Lambda_{j}=n+\frac{1}{2}.

(113)

In other words Konoplya’s formulation for the quasinormal modes encompasses various components. Specifically, the term $(V^{\prime\prime}_{0}$ denotes the second derivative of the potential, calculated at its zenith $(r_{0}$ ). Additionally, the constants $\Lambda_{j}$ are influenced by the effective potential and its derivatives at this peak. Remarkably, recent progress in this domain has unveiled a 13th-order WKB approximation pioneered by Matyjasek and Opala Matyjasek and Opala (2017), which markedly elevates the precision of quasinormal frequency computations.

Notice that the quasinormal frequencies associated with the scalar field have a negative imaginary component. This defining characteristic suggests that these modes experience exponential decay over time, representing the energy dissipation through scalar waves. This feature agrees with previous studies investigating scalar, electromagnetic, and gravitational perturbations in spherically symmetric configurations Konoplya and Zhidenko (2011); Berti et al. (2009); Heidari et al. (2023a); Chen et al. (2023).

Table 1 shows quasinormal frequencies calculated via $6$ -th order WKB approximation for spin-zero particles considering different values of $\lambda$ and fixing the other parameters as $\mathcal{M}=\alpha=1$ and $\Lambda=10^{-10}$ . We can notice that as $\lambda$ grows, the real and imaginary parts of the frequency converge to a stable value, which corresponds to the Schwarzschild quasinormal frequency.

spin 0	$l=1,n=0$	$l=2,n=0$	$l=2,n=1$
$\lambda$	$\omega_{n}$	$\omega_{n}$	$\omega_{n}$
$10^{1}$	0.300356 - 0.1005030 $i$	0.495007 - 0.0995914 $i$	0.474429 - 0.304355 $i$
$10^{2}$	0.292984 - 0.0977892 $i$	0.483756 - 0.0967944 $i$	0.463953 - 0.295714 $i$
$10^{3}$	0.292910 - 0.0977619 $i$	0.483643 - 0.0967664 $i$	0.463848 - 0.295628 $i$
$10^{4}$	0.292910 - 0.0977616 $i$	0.483642 - 0.0967661 $i$	0.463847 - 0.295627 $i$
$10^{5}$	0.292910 - 0.0977616 $i$	0.483642 - 0.0967661 $i$	0.463847 - 0.295627 $i$

Table 1: Quasinormal frequencies calculated via

6

-th order WKB approximation for spin-zero particles varying

\lambda

with fixed parameters,

\mathcal{M}=\alpha=1

and

\Lambda=10^{-10}

5.2 Electromagnetic perturbations

In this section we move forward with the examination of the electromagnetic field’s propagation. To accomplish this we recall the wave equations obeyed by a test electromagnetic field,

\frac{1}{\sqrt{-g}}\partial_{\nu}\left[\sqrt{-g}g^{\alpha\mu}g^{\sigma\nu}\left(A_{\sigma,\alpha}-A_{\alpha,\sigma}\right)\right]=0.

(114)

The four-potential, denoted as $A_{\mu}$ , can be expanded in 4-dimensional vector spherical harmonics in the following manner Toshmatov et al. (2017),

	$\displaystyle A_{\mu}\left(t,r,\theta,\phi\right)$
	$\displaystyle=\sum_{l,m}\begin{bmatrix}f(t,r)Y_{lm}\left(\theta,\phi\right)\\ h(t,r)Y_{lm}\left(\theta,\phi\right)\\ \frac{a(t,r)}{\sin\left(\theta\right)}\partial_{\phi}Y_{lm}\left(\theta,\phi\right)+k(t,r)\partial_{\theta}Y_{lm}\left(\theta,\phi\right)\\ -a\left(t,r\right)\sin\left(\theta\right)\partial_{\theta}Y_{lm}\left(\theta,\phi\right)+k(t,r)\partial_{\varphi}Y_{lm}\left(\theta,\phi\right)\end{bmatrix}.$		(115)

In this expansion $Y_{lm}(\theta,\phi)$ represents the spherical harmonics. It is worth noting that the first term on the right-hand side exhibits a parity of $\left(-1\right)^{l+1}$ (referred to as the axial sector), while the second term has a parity of $\left(-1\right)^{l}$ (known as the polar sector). When we straightforwardly substitute this expansion into the Maxwell equations, we can derive a second-order differential equation for the radial component (see in particular Toshmatov et al. (2017)),

\frac{\mathrm{d}^{2}\Psi\left(r_{\ast}\right)}{\mathrm{d}r_{\ast}^{2}}+\left[\omega^{2}-V_{E}\left(r_{\ast}\right)\right]\Psi\left(r_{\ast}\right)=0.

(116)

For both the axial and polar sectors we obtain a second-order differential equation for the radial part with the relation $r_{\ast}=\int f^{-1}(r)dr$ denoting the tortoise coordinate. The mode $\Psi(r_{\ast})$ is a linear combination of the functions $a(t,r)$ , $f(t,r)$ , $h(t,r)$ , and $k(t,r),$ but the specific functional dependence varies according to the parity. For the axial sector the mode is expressed as follows,

a(t,r)=\Psi\left(r_{\ast}\right),

(117)

whereas for the polar sector it is given by,

\Psi\left(r_{\ast}\right)=\frac{r^{2}}{l(l+1)}\left[\partial_{t}h(t,r)-\partial_{r}f(t,r)\right].

(118)

The corresponding effective potential in our case is found to be,

\displaystyle V_{eff}(r)=f(r)\left(\frac{l(l+1)}{r^{2}}\right),

(119)

with

\displaystyle f(r)=1-\frac{2\mathcal{M}}{r}+\frac{\mathcal{M}\alpha r}{(1+\alpha)\lambda^{2}}-\frac{\Lambda r^{2}}{3}.

In Fig. 7 we display the effective potential, considering the vectorial perturbations for different values of $l$ . As we did before, we have considered $\mathcal{M}=\alpha=1$ , $\lambda=10^{5}$ , and $\Lambda=10^{-5}$ . We can see from this figure that the barrier grows as $l$ increases.

Table 2 shows quasinormal frequencies performed via $6$ -th order WKB approximation for spin-one particles for several values of $\lambda$ keeping the fixed parameters $\mathcal{M}=\alpha=1$ and $\Lambda=10^{-10}$ . As mentioned in the scalar case, we again notice that as $\lambda$ increases, both parts of the frequencies tend to the Schwarzschild asymptotical value.

spin 1	$l=1,n=0$	$l=2,n=0$	$l=2,n=1$
$\lambda$	$\omega_{n}$	$\omega_{n}$	$\omega_{n}$
$10^{1}$	0.253368 - 0.0953556 $i$	0.467604 - 0.0978352 $i$	0.445768 - 0.299454 $i$
$10^{2}$	0.248244 - 0.0926643 $i$	0.457694 - 0.0950395 $i$	0.436627 - 0.290816 $i$
$10^{3}$	0.248192 - 0.0926373 $i$	0.457594 - 0.0950114 $i$	0.436535 - 0.290729 $i$
$10^{4}$	0.248191 - 0.0926370 $i$	0.457593 - 0.0950112 $i$	0.436534 - 0.290728 $i$
$10^{5}$	0.248191 - 0.0926371 $i$	0.457593 - 0.0950112 $i$	0.436534 - 0.290728 $i$

Table 2: Quasinormal frequencies performed via

6

-th order WKB approximation for spin-one particles varying

\lambda

with fixed parameters,

\mathcal{M}=\alpha=1

and

\Lambda=10^{-10}

5.3 Tensorial perturbations

Let us now investigate the QNMs for tensorial (gravitational) perturbations. In particular, we shall follow closely Ref.Kim (2004), where the gravitational perturbations were investigated. In general, one can write the axially symmetric spacetime as,

	$\displaystyle\mathrm{d}s^{2}$	$\displaystyle=$	$\displaystyle-e^{2\nu}\mathrm{d}t^{2}+e^{2\psi}(\mathrm{d}\phi-q_{1}\mathrm{d}t-q_{2}\mathrm{d}r-q_{3}\mathrm{d}\theta)^{2}$		(120)
		$\displaystyle+$	$\displaystyle e^{2\mu_{2}}\mathrm{d}r^{2}+e^{2\mu_{3}}\mathrm{d}\theta^{2}.$		(120)

In the case of unperturbed black hole spacetime one can write

e^{2\nu}=f(r),\quad e^{-2\mu_{2}}=\left(1-\frac{2m(r)}{r}\right)=\frac{\Delta}{r^{2}},

(121)

along with

\Delta=r^{2}-2m(r)r,\quad e^{\mu_{3}}=r,\quad e^{\psi}=r\sin\theta,

(122)

and

q_{1}=q_{2}=q_{3}=0,

(123)

where the metric (15) has been rewritten as,

\displaystyle f(r)=1-\frac{2m(r)}{r},

(124)

such that the mass function is given by,

\displaystyle m(r)=\mathcal{M}-\frac{\mathcal{M}\alpha r^{2}}{2\lambda^{2}(1+\alpha)}+\frac{\Lambda r^{3}}{6}.

(125)

It is well known that axial perturbations are characterized by $q_{1}$ , $q_{2}$ , and $q_{3}$ . Note that for the linear perturbations $\delta\nu,\delta\psi,\delta\mu_{2},\delta\mu_{3}$ one has polar perturbations with even parity. This type of perturbation will not be considered in the present paper. Now from Einstein’s equations we have,

(e^{3\psi+\nu-\mu_{2}-\mu_{3}}Q_{23})_{,3}=-e^{3\psi-\nu-\mu_{2}+\mu_{3}}Q_{02,0},

(126)

where $x^{2}=r,x^{3}=\theta$ and $Q_{AB}=q_{A,B}-q_{B,A},Q_{A0}=q_{A,0}-q_{1,A}$ Kim (2004). In particular, this can further be written as,

\frac{\sqrt{f(r)}}{\sqrt{\Delta}}\frac{1}{r^{3}\sin^{3}\theta}\frac{\partial Q}{\partial\theta}=-(q_{1,2}-q_{2,0})_{,0},

(127)

with $Q$ given by

Q(t,r,\theta)=\Delta Q_{23}\sin^{3}\theta=\Delta(q_{2,3}-q_{3,2})\sin^{3}\theta.

(128)

Having in mind another important equation,

(e^{3\psi+\nu-\mu_{2}-\mu_{3}}Q_{23})_{,2}=e^{3\psi-\nu+\mu_{2}-\mu_{3}}Q_{03,0},

(129)

it can be shown that

\frac{\sqrt{f(r)}\sqrt{\Delta}}{r^{3}\sin^{3}\theta}\frac{\partial Q}{\partial\theta}=(q_{1,3}-q_{3,0})_{,0}.

(130)

We can further show that by using $Q(r,\theta)=Q(r)C^{-3/2}_{l+2}(\theta)$ , in which we have the Gegenbauer function $C^{\nu}_{n}(\theta)$ that satisfies the following differential equation Kim (2004),

\left[\frac{\mathrm{d}}{\mathrm{d}\theta}\sin^{2\nu}\theta\frac{\mathrm{d}}{\mathrm{d}\theta}+n(n+2\nu)\sin^{2\nu}\theta\right]C^{\nu}_{n}(\theta)=0,

(131)

then,

r\sqrt{f(r)\Delta}\frac{\mathrm{d}}{\mathrm{d}r}\left(\frac{\sqrt{f(r)\Delta}}{r^{3}}\frac{\mathrm{d}Q}{\mathrm{d}r}\right)-\mu^{2}\frac{f(r)}{r^{2}}Q+\omega^{2}Q=0,

(132)

with $\mu^{2}=(l-1)(l+2)$ . Finally, if we introduce $Q=rZ$ along with $\frac{\mathrm{d}}{\mathrm{d}r_{*}}=\sqrt{f(r)\Delta}\frac{1}{r}\frac{\mathrm{d}}{\mathrm{d}r}$ , we obtain

\left(\frac{\mathrm{d}^{2}}{\mathrm{d}r^{*2}}+\omega^{2}-V(r)\right)Z=0,

(133)

with the following potential,

V_{eff}(r)=f(r)\left(\frac{l(l+1)}{r^{2}}-\frac{6m(r)}{r^{3}}+\frac{2m^{\prime}(r)}{r^{2}}\right),

(134)

V_{eff}(r)=f(r)\left(\frac{l(l+1)}{r^{2}}-\frac{6\mathcal{M}}{r^{3}}+\frac{\mathcal{M}\alpha}{\lambda^{2}(1+\alpha)r}\right).

(135)

Notice the exciting result that the cosmological constant does not affect the potential for gravitational perturbations. In addition, in the case $\alpha\to 0$ we get the well-known potential for the Schwarzschild black hole.

In Fig. 8 we display the effective potential concerning tensorial perturbations for diverse values of $l$ , considering $\mathcal{M}=\alpha=1$ , $\lambda=10^{5}$ , and $\Lambda=10^{-5}$ . We notice the general behavior of an increasing peak as the multipole number $l$ grows.

spin 2	$l=2,n=0$	$l=2,n=1$	$l=3,n=0$
$\lambda$	$\omega_{n}$	$\omega_{n}$	$\omega_{n}$
$10^{1}$	0.381934 - 0.0913639 $i$	0.353395 - 0.281224 $i$	0.612734 - 0.0953797 $i$
$10^{2}$	0.373703 - 0.0889158 $i$	0.346368 - 0.273557 $i$	0.599577 - 0.0927294 $i$
$10^{3}$	0.373620 - 0.0888912 $i$	0.346297 - 0.273481 $i$	0.599445 - 0.0927028 $i$
$10^{4}$	0.373619 - 0.0888910 $i$	0.346297 - 0.273480 $i$	0.599443 - 0.0927025 $i$
$10^{5}$	0.373619 - 0.0888910 $i$	0.346297 - 0.273480 $i$	0.599443 - 0.0927025 $i$

Table 3: Quasinormal frequencies computed via

6

-th order WKB approximation for spin–two particles with varying

\lambda

and fixed parameters

\mathcal{M}=\alpha=1

and

\Lambda=10^{-10}

Table 3 presents the numerical values for the real part and the imaginary part for some gravitational modes obtained through the 6-th order WKB method. For smaller values of $\lambda$ the value of the real part is shown to be more significant compared to the Schwarzschild case; however, the more precise astrophysical values we use for $\lambda$ , the closer the Schwarzschild case is obtained. This means that the Yukawa black hole corrections and the Schwarzschild case are indistinguishable.

Our results in this section show that Yukawa black hole metric is stable against scalar, electromagnetic, and gravitational perturbations since the imaginary part of all the frequencies we found is always negative. Following our aim to study the outcome of the gravitational field of this black hole, in the next section we will discuss the effect on the test particle motion in the neighborhood of this geometry focusing in the lightlike trajectories.

6 Geodesic trajectories

6.1 The complete case

In this section our primary focus lies in understanding the behavior resulting from the geodesic equation. This involves expressing it as follows,

\frac{\mathrm{d}^{2}x^{\mu}}{\mathrm{d}s^{2}}+\Gamma\indices{{}^{\mu}_{\alpha}{}_{\beta}}\frac{\mathrm{d}x^{\alpha}}{\mathrm{d}s}\frac{\mathrm{d}x^{\beta}}{\mathrm{d}s}=0.

(136)

Taking into account Eq. (14) the above equation leads to four partial differential equations

	$\displaystyle\frac{\mathrm{d}\theta^{\prime}}{\mathrm{d}s}=\sin(\theta)\cos(\theta)\left(\varphi^{\prime}\right)^{2}-\frac{2\theta^{\prime}r^{\prime}}{r},$	(137)
	$\displaystyle\frac{\mathrm{d}\varphi^{\prime}}{\mathrm{d}s}=-\frac{2\varphi^{\prime}\left(r^{\prime}+r\theta^{\prime}\cot(\theta)\right)}{r},$	(138)
$\displaystyle\begin{split}&\frac{\mathrm{d}t^{\prime}}{\mathrm{d}s}=-\frac{1}{\frac{2Mr^{2}}{\left(l_{0}^{2}+r^{2}\right)^{3/2}}+\frac{2\alpha Me^{-\frac{r}{\lambda}}(\lambda+r)}{\lambda r}+\frac{\Lambda r^{2}}{3}-1}\\ &\times r^{\prime}t^{\prime}\left[M\left(\frac{4r}{\left(l_{0}^{2}+r^{2}\right)^{3/2}}-\frac{6r^{3}}{\left(l_{0}^{2}+r^{2}\right)^{5/2}}-\frac{2\alpha e^{-\frac{r}{\lambda}}}{r^{2}}\right.\right.\\ &\left.\left.-\frac{2\alpha e^{-\frac{r}{\lambda}}}{\lambda^{2}}-\frac{2\alpha e^{-\frac{r}{\lambda}}}{\lambda r}\right)+\frac{2\Lambda r}{3}\right],\end{split}$		(139)

\displaystyle\begin{split}&\frac{\mathrm{d}r^{\prime}}{\mathrm{d}s}=\frac{e^{-\frac{r}{\lambda}}}{6\lambda\left(l_{0}^{2}+r^{2}\right)^{3/2}}\\ &\times\left\{\frac{1}{r\left(\frac{2Mr^{2}}{\left(l_{0}^{2}+r^{2}\right)^{3/2}}+\frac{2\alpha Me^{-\frac{r}{\lambda}}(\lambda+r)}{\lambda r}+\frac{\Lambda r^{2}}{3}-1\right)^{2}}\right.\\ &\left.\times\left[\left(2M\left(-\frac{2r}{\left(l_{0}^{2}+r^{2}\right)^{3/2}}+\frac{3r^{3}}{\left(l_{0}^{2}+r^{2}\right)^{5/2}}+\frac{\alpha e^{-\frac{r}{\lambda}}}{r^{2}}\right.\right.\right.\right.\\ &\left.\left.\left.\left.+\frac{\alpha e^{-\frac{r}{\lambda}}}{\lambda^{2}}+\frac{\alpha e^{-\frac{r}{\lambda}}}{\lambda r}\right)-\frac{2\Lambda r}{3}\right)\right.\right.\\ &\left.\left.\times\left[\left(-l_{0}^{2}\sqrt{a^{2}+r^{2}}\left(6\alpha M(\lambda+r)+\lambda re^{r/\lambda}\left(\Lambda r^{2}-3\right)\right)\right.\right.\right.\right.\\ &\left.\left.\left.\left.-r^{2}\left(6M\left(\alpha\lambda\sqrt{l_{0}^{2}+r^{2}}+\alpha r\sqrt{l_{0}^{2}+r^{2}}+\lambda re^{r/\lambda}\right)\right.\right.\right.\right.\right.\\ &\left.\left.\left.\left.\left.+\lambda r\sqrt{a^{2}+r^{2}}e^{r/\lambda}\left(\Lambda r^{2}-3\right)\right)\right)r^{\prime 2}\right]\right.\right.\\ &\left.\left.-\frac{1}{r}\times\left[2M\left(-\frac{2r}{\left(l_{0}^{2}+r^{2}\right)^{3/2}}+\frac{3r^{3}}{\left(l_{0}^{2}+r^{2}\right)^{5/2}}+\frac{\alpha e^{-\frac{r}{\lambda}}}{r^{2}}\right.\right.\right.\right.\\ &\left.\left.\left.\left.+\frac{\alpha e^{-\frac{r}{\lambda}}}{\lambda^{2}}+\frac{\alpha e^{-\frac{r}{\lambda}}}{\lambda r}\right)-\frac{2\Lambda r}{3}\right]\right.\right.\\ &\left.\left.\times\left[-l_{0}^{2}\sqrt{l_{0}^{2}+r^{2}}\left(6\alpha M(\lambda+r)+\lambda re^{r/\lambda}\left(\Lambda r^{2}-3\right)\right)\right.\right.\right.\\ &\left.\left.\left.+\lambda r\sqrt{l_{0}^{2}+r^{2}}e^{r/\lambda}(\Lambda r^{2}-3)\right]\right]t^{\prime 2}\right.\\ &\left.-2\times\left[l_{0}^{2}\sqrt{l_{0}^{2}+r^{2}}\left(6\alpha M(\lambda+r)+\lambda re^{r/\lambda}\left(\Lambda r^{2}-3\right)\right)\right.\right.\\ &\left.\left.+r^{2}\left(6M\left(\alpha\lambda\sqrt{l_{0}^{2}+r^{2}}+\alpha r\sqrt{l_{0}^{2}+r^{2}}\right.\right.\right.\right.\\ &\left.\left.\left.\left.+\lambda re^{r/\lambda}\right)+\lambda r\sqrt{l_{0}^{2}+r^{2}}e^{r/\lambda}\left(\Lambda r^{2}-3\right)\right)\right]\theta^{\prime 2}\right.\\ &\left.-2\left[l_{0}^{2}\sqrt{l_{0}^{2}+r^{2}}\left(6\alpha M(\lambda+r)+\lambda re^{r/\lambda}\left(\Lambda r^{2}-3\right)\right)\right.\right.\\ &\left.\left.+r^{2}\left(6M\left(\alpha\lambda\sqrt{l_{0}^{2}+r^{2}}+\alpha r\sqrt{l_{0}^{2}+r^{2}}+\lambda re^{r/\lambda}\right)\right.\right.\right.\\ &\left.\left.\left.\lambda r\sqrt{l_{0}^{2}+r^{2}}+e^{r/\lambda}\left(\Lambda r^{2}-3\right)\right)\right]\sin^{2}(\theta)\varphi^{\prime 2}\right\}.\end{split}

(140)

In Fig. 9, we display the light path for $l_{0}=0$ . Here the more $\lambda$ increases, the larger the deviation becomes, represented by the orange lines (light path). The black circle represents the event horizon and the dashed points are the photon sphere of the black hole.

6.2 The approximated case

Let us discuss the case when the geodesic trajectories are obtained from the approximated solution of $f(r)$ in Eq. (15) when the parameters $\varepsilon=1/\lambda$ and $\ell_{0}$ approach zero. With these approximations the equations are simplified significantly, as we can see in what follows,

	$\displaystyle\frac{\mathrm{d}\theta^{\prime}}{\mathrm{d}s}=\sin(\theta)\cos(\theta)\left(\varphi^{\prime}\right)^{2}-\frac{2\theta^{\prime}r^{\prime}}{r},$	(141)
	$\displaystyle\frac{\mathrm{d}\varphi^{\prime}}{\mathrm{d}s}=-\frac{2\varphi^{\prime}\left(r^{\prime}+r\theta^{\prime}\cot(\theta)\right)}{r},$	(142)
	$\displaystyle\frac{\mathrm{dt^{\prime}}}{\mathrm{d}s}=\frac{r^{\prime}t^{\prime}\left(\frac{\alpha\mathcal{M}}{(\alpha+1)\lambda^{2}}+\frac{2\mathcal{M}}{r^{2}}-\frac{2\Lambda}{3}r\right)}{\frac{\Lambda r^{2}}{3}-\frac{\alpha\mathcal{M}r}{(\alpha+1)\lambda^{2}}+\frac{2\mathcal{M}}{r}-1},$	(143)
$\displaystyle\begin{split}&\frac{\mathrm{d}r^{\prime}}{\mathrm{d}s}=\\ &\frac{\left(r^{\prime}\right)^{2}\left(-6(\alpha+1)\lambda^{2}\mathcal{M}+2(\alpha+1)\lambda^{2}\Lambda r^{3}-3\alpha\mathcal{M}r^{2}\right)}{2r\left(6(\alpha+1)\lambda^{2}\mathcal{M}+(\alpha+1)\lambda^{2}\Lambda r^{3}-3\alpha\mathcal{M}r^{2}-3(\alpha+1)\lambda^{2}r\right)}\\ &-\frac{1}{18(\alpha+1)^{2}\lambda^{4}r^{3}}\\ &\times\left\{\left[6(\alpha+1)\lambda^{2}\mathcal{M}+(\alpha+1)\lambda^{2}\Lambda r^{3}-3\alpha\mathcal{M}r^{2}-3(\alpha+1)\lambda^{2}r\right]\right.\\ &\left.\times\left[\left.-6(\alpha+1)\lambda^{2}\mathcal{M}+2(\alpha+1)\lambda^{2}\Lambda r^{3}-3\alpha\mathcal{M}r^{2}\right)t^{\prime 2}\right]\right.\\ &\left.\times\left[+6r^{3}\alpha+1\lambda^{2}\left(\theta^{\prime 2}+\sin^{2}(\theta)\varphi^{\prime 2}\right.\right]\right\}.\end{split}$		(144)

In Fig. 10 we display the behavior of the light path close to the black hole. Intriguingly, after around $\lambda>30$ the deflection of light is practically not altered for such a parameter.

Light paths corresponding to the set of last unstable orbits around the compact object, known as photon sphere, give rise to an interesting phenomenon that will be studied in the next section, namely the shadow of a black hole.

7 Photon spheres and Shadows: the case of optically thin radiating and infalling gas surrounding a black hole

To fully grasp particle dynamics and light ray patterns near black holes, it is essential to understand photon spheres. These spheres are crucial in discerning the shadows, shedding light on the effects of the Yukawa potential in the spacetime under consideration.

Understanding the significance of the photon sphere, commonly labelled as the critical orbit, in the context of black hole dynamics requires, for instance, the application of the Lagrangian method for calculating null geodesics. We aim to ascertain the effects of parameters, specifically $M,\Lambda,\alpha,\lambda$ , on the photon sphere. To elucidate further, we state

\mathcal{L}=\frac{1}{2}g_{\mu\nu}\dot{x}^{\mu}\dot{x}^{\nu}.

(145)

When fixing the angle at $\theta=\pi/2$ , the previously mentioned expression simplifies to

g_{00}^{-1}E^{2}+g_{11}^{-1}\dot{r}^{2}+g_{33}^{-1}L^{2}=0,

(146)

where $L$ represents the angular momentum and $E$ denotes the energy. Subsequently, Eq. (146) is given by,

\dot{r}^{2}=E^{2}-\left(\frac{\alpha\mathcal{M}r}{(\alpha+1)\lambda^{2}}-\frac{2\mathcal{M}}{r}-\frac{\Lambda r^{2}}{3}+1\right)\left(\frac{L^{2}}{r^{2}}\right),

(147)

with

\mathcal{V}\equiv\left(\frac{\alpha\mathcal{M}r}{(\alpha+1)\lambda^{2}}-\frac{2\mathcal{M}}{r}-\frac{\Lambda r^{2}}{3}+1\right)\left(\frac{L^{2}}{r^{2}}\right)\,,

(148)

being the effective potential. To determine the critical radius, the equation $\partial\mathcal{V}/\partial r=0$ must be solved, which can be expressed as,

\displaystyle r_{\rm c}=\frac{\sqrt{(1+\alpha)\left(6\mathcal{M}^{2}\lambda^{2}\alpha+\lambda^{4}(1+\alpha)\right)}-\lambda^{2}(1+\alpha)}{\mathcal{M}\alpha}.

(149)

Here a remarkable feature can be singled out, the cosmological constant does not affect the critical orbits.

Moreover, a unique physical solution exists for this equation, denoted as $r_{c}$ – the photon sphere radius. Recent literature, such as Refs. Filho (2023), has tackled similar studies within dark matter contexts. Also, it is noteworthy to highlight that the emergence of two-photon spheres has been recently documented in the Simpson–Visser solution context Tsukamoto (2021, 2022) among others Guerrero et al. (2022). The left side of Table 4 shows various results for the shadow radius corresponding to different values of $\lambda$ while keeping $\alpha=1$ . In this scenario as $\lambda$ increases, the shadow radius decreases. Given that this parameter is directly tied to the Yukawa potential, our findings suggest that for astrophysical values of the parameters the overall influence on the shadow radius is very small. In terms of the graviton mass $\lambda=\hbar/(m_{g}c)$ , this means the smaller the graviton mass is, the stronger the effect becomes. As the graviton mass decreases, the shadow radius also decreases.

Examining shadows in the context of Yukawa potential and black hole structures is of paramount importance. Defined by the distinct outline of a black hole set against a bright background, these shadows offer perspectives on spacetime geometry and gravitational phenomena near the black hole. By analyzing these shadows, we can glean valuable information to test and enhance theoretical models further authenticating the nature of gravity. To streamline our analysis we introduce two new parameters as follows,

\xi=\frac{L}{E}\text{ and }\eta=\frac{\mathcal{K}}{{E^{2}}},

(150)

where $\mathcal{K}$ is the so–called Carter constant. After performing algebraic manipulations we obtain,

{\xi^{2}}+\eta=\frac{r_{c}^{2}}{f(r_{c})}.

(151)

In the endeavor to deduce the radius of the shadow we shall make use of the celestial coordinates $X$ and $Y$ Singh and Ghosh (2018); Heidari et al. (2023a); Filho et al. (2023d); Heidari et al. (2023b); Filho et al. (2023e) as follows, $X=-\xi$ and $Y=\pm\sqrt{\eta}$ . It is natural to define the true or the physical mass of the black hole to be $\mathcal{M}=M(1+\alpha)$ . Furthermore, in the case of large distances with a cosmological constant it has been shown that the size of the BH shadow can explicitly depend on the radial coordinate of the distant observer noted as $r_{O}$ . Following González et al. (2023) we have

\displaystyle\sin^{2}\alpha_{\rm sh}=\frac{r^{2}_{c}}{f(r_{c})}\frac{f(r_{O})}{r_{O}}\,.

(152)

In the physically relevant small-angle approximation it is easy to see that the shadow size is given by,

\displaystyle R_{\rm SH}=r_{c}\sqrt{\frac{f(r_{O})}{f(r_{c})}}\,.

(153)

It is easy to show that the shadow radius can be written as González et al. (2023),

R_{\rm SH}=\frac{r_{c}}{\sqrt{f(r_{c})}}\sqrt{1-\frac{2GM_{\rm BH}}{c^{2}r_{O}}+\frac{GM_{\rm BH}\alpha r_{O}}{c^{2}(1+\alpha)\lambda^{2}}-\frac{1}{3}\Lambda r_{O}^{2}},

(154)

where $r_{O}$ is the distance to the black hole.

Now let us consider a straightforward model, specifically an optically sparse accretion flow that emits radiation around the object. We will employ a numerical method called Backward Raytracing to uncover the shadow created by this radiating flow. To determine the intensity distribution within the emitting region, we must make certain assumptions about the radiative processes and emission mechanisms in play. In particular, the observed specific intensity $I_{\nu 0}$ at the observed photon frequency $\nu_{\text{obs}}$ at the point $(X,Y)$ of the observer’s image (usually measured in $\text{erg}\text{s}^{-1}\text{cm}^{-2}\text{str}^{-1}\text{Hz}^{-1}$ ) is given by, Bambi (2013); Saurabh and Jusufi (2021)

\displaystyle I_{obs}(\nu_{obs},X,Y)=\int_{\gamma}\mathrm{g}^{3}j(\nu_{e})dl_{\text{prop}},

(155)

In this analysis we are examining a simplified scenario involving accreting gas. Our assumption is that the gas undergoes radial free fall and its four-velocity, under static and spherically symmetric conditions, simplifies to the following Bambi (2013); Saurabh and Jusufi (2021),

	$\displaystyle u^{t}_{e}$	$\displaystyle=$	$\displaystyle\frac{1}{f(r)},\;u^{r}_{e}=-\sqrt{\frac{g(r)}{f(r)}\left(1-f(r)\right)},$
	$\displaystyle u^{\theta}_{e}$	$\displaystyle=$	$\displaystyle 0,\;u^{\phi}_{e}=0,$		(156)

where

\displaystyle f(r)=g(r)=1-\frac{2\mathcal{M}}{r}+\frac{\mathcal{M}\alpha r}{(1+\alpha)\lambda^{2}}-\frac{\Lambda r^{2}}{3}.

(157)

Using the four-velocity for the photons we can establish a connection between the radial and temporal components of the four-velocity,

\displaystyle\frac{k^{r}}{k^{t}}=\pm f(r)\sqrt{g(r)\bigg{(}\frac{1}{f(r)}-\frac{b^{2}}{r^{2}}\bigg{)}}.

(158)

The sign +(-) represents the direction in which the photon is either approaching or moving away from the massive object, respectively. As a result, the redshift function $\mathrm{g}$ is expressed as follows,

\displaystyle\mathrm{g}=\frac{1}{\frac{1}{f(r)}\pm\frac{k_{r}}{k_{t}}\sqrt{\frac{g(r)}{f(r)}\bigg{(}1-f(r)\bigg{)}}}.

(159)

Regarding the specific emissivity, we adopt a straightforward model in which the emission is monochromatic, having an emitter’s rest-frame frequency denoted as $\nu_{\star}$ . Additionally, the emission exhibits a radial profile that follows an inverse square law, i.e., $1/r^{2}$ ,

\displaystyle j(\nu_{e})\propto\frac{\delta(\nu_{e}-\nu_{\star})}{r^{2}},

(160)

where $\delta$ is the Dirac delta function along with the proper length which can be written as,

\displaystyle dl_{\text{prop}}=k_{\alpha}u^{\alpha}_{e}d\lambda=-\frac{k_{t}}{\mathrm{g}|k^{r}|}dr.

(161)

By integrating the intensity across all observable frequencies, we arrive at the observed flux,

\displaystyle F_{obs}(X,Y)\propto-\int_{\gamma}\frac{\mathrm{g}^{3}k_{t}}{r^{2}k^{r}}dr.

(162)

$\mathcal{M}$	$\alpha$	$\Lambda$ [m^-2]	$\lambda/\mathcal{M}$	$R_{SH}/\mathcal{M}$
1.0	1.0	$10^{-52}$	$10^{5}$	6.363961027
1.0	1.0	$10^{-52}$	$10^{6}$	5.209126608
1.0	1.0	$10^{-52}$	$10^{8}$	5.196153722
1.0	1.0	$10^{-52}$	$10^{13}$	5.196152422
1.0	1.0	$10^{-52}$	$10^{15}$	5.196152422

Table 4: The shadow radius is presented for varying values of mass and the parameter

\lambda

. We have set

r_{0}/\mathcal{M}=10^{10}

In order to estimate the shadow radius for the Sgr A* black hole, we used real astrophysical values, $\Lambda\simeq 10^{-52}\rm m^{-2}$ . Also, for this black hole we can take (here we have written the black hole mass in geometric units using the fact that $M_{\rm Sun}=1.5$ km),

\displaystyle M_{\rm BH}^{\rm SgrA}=4\times 10^{6}M_{\rm Sun}=6\times 10^{9}\rm m.

(163)

For the distance we have $r_{O}=8.3$ kpc and $\lambda\sim 2693$ Mpc or

\displaystyle\frac{\lambda}{M_{\rm BH}^{\rm SgrA}}=1.382406667\times 10^{16},

(164)

along with $\alpha=0.416$ González et al. (2023). We work in geometric units; hence, we use $M_{\rm Sun}=1.5\,km$ . For the shadow radius we obtain

\displaystyle R_{\rm SH}/M_{\rm BH}^{\rm SgrA}=5.196152423.

(165)

For the vacuum solution or the Schwarzschild black hole we know that

\displaystyle R^{\rm Schwarzschild}_{\rm SH}/M_{\rm BH}^{\rm SgrA}=3\sqrt{3}\simeq 5.196152424.

(166)

This means that the change in the shadow radius compared to the Schwarzschild black hole is of the order

\displaystyle\delta R_{\rm SH}=\frac{R^{\rm Schwarzschild}_{\rm SH}-R_{\rm SH}}{M_{\rm BH}^{\rm SgrA}}\sim 10^{-9}.

(167)

In Fig. 11 and Fig. 12 we have shown the intensity of the infalling gas and the shadow images of the Yukawa black hole. Finally, let us briefly mention that recent studies Jusufi et al. (2024); Jusufi and Sheykhi (2024) pointed out the possibility of varying graviton mass. This means that $\lambda=\hbar/(m_{g}c)$ can fluctuate from Mpc order in the case of cosmological scales to kpc order in the case of galactic scales. Using kpc order for $\lambda$ , we observe a similar result in the shadow radius, i.e., a change of $|\delta R_{\rm SH}|\sim 10^{-8}$ – $10^{-9}$ compared to the Schwarzschild black hole shadow.

7.1 Can we constrain $\alpha$ ?

Note that in our discussion above we have identified the black hole mass with $\mathcal{M}=M_{\rm BH}^{\rm SgrA}$ , since as we have argued $\mathcal{M}=M+\alpha M$ . This means that the total mass, as measured by an observer located far away from the black hole, is modified due to the extra term which mimics the effect of dark matter, i.e., $M_{\rm dark\,matter}=\alpha M$ . We can interpret, on the other hand, $M$ as the “bare” mass of the black hole, while the extra term $M_{\rm dark\,matter}$ is only an apparent effect which arises from the modification of Newton’s law of gravity in large distances but changes $M\to\mathcal{M}$ as measured by an observer located far away. A natural question arises: can we probe the “bare” mass $M$ ? First, we can easily see that the shadow radius can be approximated as $R_{SH}=3\sqrt{3}M(1+\alpha)+\mathcal{O}(\Lambda,\alpha,\lambda^{2})$ , and by keeping only the first term, we can apply the same approach as in Vagnozzi et al. (2023). Using observations near the black hole, we can effectively probe $M$ and not $\mathcal{M}$ , for example, by studying the motion of S-stars orbiting Sgr A*. In addition, there is uncertainty about the mass, which does not disappear when we later fix $M=1$ , it is part of what allows us to constrain the dark matter effect $\alpha$ : so we can fix the bare mass and ask how much of a change in the total mass/energy of the system we can tolerate when varying the parameter $\alpha$ . Using the EHT observations for Sgr A*, it was shown that within $2\sigma$ constraints Vagnozzi et al. (2023),

\displaystyle 4.21\lesssim R_{sh}/M\lesssim 5.56\,,

(168)

which gives the interval $-0.19\lesssim\alpha\lesssim 0.07$ (within 2 $\sigma$ ). On the other hand, for the M87 black hole it was shown within $2\sigma$ constraints Allahyari et al. (2020),

\displaystyle 4\lesssim R_{sh}/M\lesssim 7,

(169)

which gives the interval $-0.23\lesssim\alpha\lesssim 0.35$ (within 2 $\sigma$ ). This fact is also shown in Fig. 13. We see that $\alpha$ can be negative too, but in the present paper we considered only $\alpha>0$ . This shows that nearly $0.35M$ (for M87) and $0.07M$ (Sgr A*) of the total mass can be attributed to the apparent dark matter effect.

In the next section we will make the connection between our results demonstrating in an analytical way that the Yukawa corrected black hole also obeys the correspondence between shadow radius and QNMs in the eikonal limit.

7.2 The galactic flat rotating curves

Let us elaborate briefly on how the modification of the gravitational potential can lead to interesting results related to the flat rotating galactic curves. In such a case, the black hole mass does not play an important role and should be replaced by the mass distribution of the surrounding matter. Assuming a constant mass distribution for simplicity, in the limit of large radial distance we can write,

\displaystyle f(r)\simeq 1+\frac{\mathcal{M}\alpha r}{(1+\alpha)\lambda^{2}}

(170)

where we have neglected the second term $2\mathcal{M}/r$ , which is small for a very large radial distance, and the last term $\Lambda r^{2}/3$ since the cosmological constant is small and does not play an important role in galactic scales. The second term in the last equation plays the role of dark matter. For the case of a spherically symmetric spacetime ansatz with pure dark matter in Schwarzschild coordinates the tangential velocity using the radial function $f(r)$ can be calculated from the following equation,

\displaystyle v_{0}^{2}=\frac{r}{\sqrt{f(r)}}\frac{d}{dr}\left(\sqrt{f(r)}\right)

(171)

which describes a region with a nearly constant tangential speed of the star. From the last equation we get Saurabh and Jusufi (2021),

\displaystyle f(r)=\left(\frac{r}{r_{0}}\right)^{2v_{0}^{2}}\simeq 1+2v_{0}^{2}\ln(r/r_{0}),

(172)

where $r_{0}$ is some integration constant. For the outer part of the galaxy we can take $\ln(r/r_{0})\sim 1$ and $r\sim\lambda$ , that leads to,

\displaystyle v_{0}=\sqrt{\frac{G\mathcal{M}\alpha}{2(1+\alpha)\lambda}}=\rm constant,

(173)

where we have inserted the Newton’s constant $G$ . We already pointed out that $\lambda$ in cosmological scales was found to be of Mpc order. However, as was pointed out recently in Jusufi et al. (2024); Jusufi and Sheykhi (2024), the graviton mass related to $\lambda$ by $\lambda=\hbar/(m_{g}c)$ can fluctuate –for galactic scales $\lambda$ indeed can be of kpc order– and, hence, the assumption for the outer part of the galaxy $r\sim\lambda$ is justified. As an example, one can take the total mass enclosed inside the galaxy to be $\mathcal{M}\sim 10^{12}$ solar masses along with $\lambda\sim 1$ kpc order and $\alpha\sim 0.5$ which yields $v_{0}\sim 200$ km/s. The idea of a mass-varying mechanism for the graviton is yet to be explored in the near future.

8 Eikonal quasinormal modes and shadow radius correspondence

As it is well known, QNMs are characterized by complex frequencies, $\omega=\omega_{\Re}-i\omega_{\Im}$ , with a real part representing the oscillation of the wave and an imaginary part proportional to the damping of a given mode. Moreover, both numbers depend only on the features of the metric being perturbed and contain essential information related to the kind of perturbation being used. In the last decades several discussions around the nature of these vibrations emerged. But Goebel pointed out a possible relation with the photon sphere trajectories for the first time Goebel (1972). Later on, Cardoso et al. Cardoso et al. (2009) showed that the real part of the QNMs in the eikonal regime is related to the angular velocity of the outermost photon orbit $\Omega_{c}$ , while the imaginary part is proportional to the Lyapunov exponent $\tilde{\lambda}$ , which determines the instability timescale of the orbit,

\omega_{QNM}=\Omega_{c}l-i\left(n+\frac{1}{2}\right)|\tilde{\lambda}|\,.

(174)

This result was proven to be valid not only for static spherically symmetric spacetimes but also for equatorial orbits in the geometry of rotating black holes solutions. On the other hand, Stefanov et al. Stefanov et al. (2010) found a connection between black-hole eikonal quasinormal modes and gravitational lensing observables (flux ratio $\tilde{r}$ and minimal impact angle $\theta_{\infty}$ ) in the strong-deflection regime,

\Omega_{c}=\frac{c}{\theta_{\infty}D_{OL}}\,,\qquad\tilde{\lambda}=\frac{c\,\ln\tilde{r}}{2\pi\theta_{\infty}D_{OL}}\,,

(175)

where $D_{OL}$ is the distance between the observer and the lens. Afterwards, a connection with the shadow radius of a black hole was put forward in Jusufi (2020) and analytically proven in Cuadros-Melgar et al. (2020). This correspondence shows that the real part of the quasinormal frequency is proportional to the inverse of the shadow radius of the same black hole.

Nevertheless, we should stress that there are counterexamples to the correspondence. Even for high multipole numbers $l$ , the correspondence is not guaranteed for any field. As stated in Konoplya and Stuchlík (2017), the correspondence works only for test fields whose effective potential is positive definite, has a barrier, and goes to a constant value both at the event horizon and at infinity (or dS horizon), such that a WKB method can be applied. However, even under these conditions, some frequencies could be missed when using the standard WKB method Konoplya (2023).

To see if the correspondence holds for the black hole corrected by the Yukawa potential, we are going to calculate the quasinormal frequencies in the eikonal limit. Following Cuadros-Melgar et al. (2020), we consider the effective potential for the scalar perturbation given in Eq. (111). Using the 6^th order WKB method we can write the complex frequency $\omega$ as,

$\displaystyle\omega^{2}$	$\displaystyle=$	$\displaystyle V_{0}+\frac{V_{4}}{8V_{2}}\left(\nu^{2}+\frac{1}{4}\right)-\frac{(7+60\nu^{2})V_{3}^{2}}{288V_{2}^{2}}$
		$\displaystyle+i\nu\sqrt{-2V_{2}}\left[\frac{1}{2V_{2}^{2}}\left(\frac{5V_{3}^{4}(77+188\nu^{2})}{6912V_{2}^{3}}\right.\right.$
		$\displaystyle-\frac{V_{3}^{2}V_{4}(51+100\nu^{2})}{384V_{2}^{2}}+\frac{V_{4}^{2}(67+68\nu^{2})}{2304V_{2}}$
		$\displaystyle\left.\left.+\frac{V_{5}V_{3}(19+28\nu^{2})}{288V_{2}}+\frac{V_{6}(5+4\nu^{2})}{288}\right)-1\right]\,,$

where $V_{i}$ represents the $i^{th}$ derivative of the potential with respect to $r$ , evaluated at the position of the maximum value of the potential $r=r_{p}$ and $\nu=n+1/2$ , being $n$ the overtone number. If we expand this expression for large $l$ , i.e., we are taking the eikonal limit, we obtain a simple expression for the real part of the quasinormal frequency,

\omega_{\Re}=\left(l+\frac{1}{2}\right)\frac{\sqrt{f(r_{p})}}{r_{p}}\,.

(177)

Now, to find the position of the peak of the potential, we notice that in the eikonal limit Eq. (111) has the same form as the lightlike geodesic potential shown in Eq. (148), provided that $l(l+1)$ is identified with $L^{2}$ . Thus, the position of the maximum of the scalar perturbation potential corresponds to the photon sphere radius found in Eq. (149), $r_{p}=r_{c}$ . Then, since the shadow radius $R_{\rm SH}$ is given by Eq. (153), we can rewrite Eq. (177) as

\omega_{\Re}=\left(l+\frac{1}{2}\right)\frac{\sqrt{f(r_{O})}}{R_{\rm SH}}\,.

(178)

Therefore, the correspondence still holds, but we have an additional factor that depends on the radial coordinate of the distant observer $r_{O}$ . Furthermore, if we use real astrophysical values for the Sgr A^∗, for example, as was done in the previous section, $f(r_{O})$ approaches unity, and we obtain the standard result found in Cuadros-Melgar et al. (2020); Jusufi (2020).

$\mathcal{M}$	$\alpha$	$\Lambda$ [m^-2]	$\lambda/\mathcal{M}$	$l$	$\omega_{\Re}$	$\delta\omega_{\Re}$
1.0	1.0	$10^{-52}$	$10^{10}$	1	0.288675134	$6.4952\times 10^{-21}$
1.0	1.0	$10^{-52}$	$10^{10}$	2	0.481125224	$1.0825\times 10^{-20}$
1.0	1.0	$10^{-52}$	$10^{10}$	3	0.673575314	$1.5155\times 10^{-20}$
1.0	1.0	$10^{-52}$	$10^{10}$	10	2.020725942	$4.5466\times 10^{-20}$
1.0	1.0	$10^{-52}$	$10^{10}$	20	3.945226840	$1.5155\times 10^{-20}$
1.0	1.0	$10^{-52}$	$10^{10}$	50	9.718729530	$2.1867\times 10^{-19}$
1.0	1.0	$10^{-52}$	$10^{10}$	100	19.34123402	$4.3518\times 10^{-19}$
1.0	1.0	$10^{-52}$	$10^{10}$	200	38.58624300	$8.6819\times 10^{-19}$
1.0	1.0	$10^{-52}$	$10^{10}$	500	96.32126990	$2.1672\times 10^{-18}$
1.0	1.0	$10^{-52}$	$10^{10}$	1000	192.5463148	$4.3323\times 10^{-18}$

Table 5: Numerical values for the eikonal QNMs for different

l

using Eq. (177) as well as

\delta\omega_{\Re}

, which gives the difference compared to the Schwarzschild black hole case,

\delta\omega_{\Re}=\omega_{\Re}-\omega_{\Re}^{\rm Schwarzschild.}

In Table 5, we present our numerical result when computing the real part of eikonal QNMs. The results show that the effect of Yukawa corrections on astrophysical data is extremely small. Specifically, from this table we can see that comparing the QNMs to the Schwarzschild black hole case, the signature of Yukawa corrections that might mimic the dark matter effect is very small, e.g, an increase for the QNMs of order $10^{-18}-10^{-21}$ . In other words, from a practical point of view, it is almost impossible to detect such corrections using gravitational waves.

9 Conclusion

In this paper we obtained a black hole solution with non–singular Yukawa–like potential given by (10), which generalizes the metric presented in González et al. (2023). The second term of this solution modifies the geometry due to $\ell_{0}$ and plays an important role in short–range distances to cure the black hole singularity. The third term is due to the apparent dark matter effect and the fourth term is the contribution due to the cosmological constant. Since $\ell_{0}$ is of Planck length order, i.e., $\ell_{0}\sim 10^{-35}\,m$ Nicolini et al. (2019), in large distances and for astrophysical black holes with large mass $M$ we must have $r\gg\ell_{0}$ , then, we can neglect the effect of $\ell_{0}$ . In short, the spacetime geometry was modified in short distances due to the deformed parameter $\ell_{0}$ and, most importantly, in large distances due to both parameters $\alpha$ and $\lambda$ . The large distance corrections were of particular interest since $\alpha$ modifies Newton’s law of gravity and could mimic the dark matter effect. We also investigated the accuracy of the approximate solution (15) using perturbation methods in differential equations.

The Yukawa corrections could be classified in two categories: i) corrections due to $\alpha$ could be significant since the $\alpha$ parameter modifies the total mass of the black hole via $\mathcal{M}=M+\alpha M$ , or ii) corrections due to $\lambda$ were very small and could practically be neglected. Keeping this in mind, we subsequently explored the phenomenological aspects of the solution in large distances for astrophysical black holes.

We performed a thermodynamic analysis by calculating Hawking temperature, entropy, and heat capacity and we showed that this black hole solution underwent first–order thermodynamical phase transitions. In addition, we studied the quasinormal frequencies for scalar, electromagnetic, and gravitational field perturbations. For doing so, we utilized the WKB method to probe the quasinormal modes numerically. Since it was more natural to measure $\mathcal{M}$ , we analyzed the QNMs using specific values of the parameters that dictate the structure of the black hole, namely, $\mathcal{M}$ , $\alpha$ , $\lambda$ , and $\Lambda$ .

Also we used astrophysical data for the parameters and an optically thin radiating and infalling gas surrounding a black hole to model the black hole shadow image. We considered as an example the Sgr A^∗ black hole and we found that the shadow radius changes by order of $10^{-9}$ , meaning that the shadow radius of a black hole with Yukawa potential practically gave the same results as those ones encountered in Schwarzschild black hole.

On the other hand, in the eikonal regime, using astrophysical data for Yukawa parameters, we showed that the value of the real part of the QNMs frequencies changes by $10^{-18}$ . As we pointed out, the major effect via $\alpha$ was encoded in $\mathcal{M}$ . However, through astrophysical observations, it is simpler to probe $\mathcal{M}$ than $M$ . In principle, one can probe $M$ with specific observations near the black hole such as the motion of S–stars. In that case we can constrain $\alpha$ and we can calculate how much the total mass changes by working in units of $M$ . In general, we have a screening effect due to $\alpha$ and if we work in units of $\mathcal{M}$ , the Yukawa–like corrections are difficult to measure by observations with the current technology, for example, using gravitational waves.

CRediT authorship contribution statement

All the authors equally contributed to Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Writing – original draft, Writing – review & editing, Visualization, Investigation, Formal analysis, Funding acquisition, Supervision, and Project administration.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

A. A. Araújo Filho would like to thank Fundação de Apoio à Pesquisa do Estado da Paraíba (FAPESQ) and Conselho Nacional de Desenvolvimento Cientíıfico e Tecnológico (CNPq) – [200486/2022-5] and [150891/2023-7] for the financial support. Genly Leon thanks Vicerrectoría de Investigación y Desarrollo Tecnológico (VRIDT) at Universidad Católica del Norte for financial support through Resolución VRIDT N°026/2023, Resolución VRIDT N°027/2023, Núcleo de Investigación Geometría Diferencial y Aplicaciones (Resolución VRIDT N°096/2022), Proyecto de Investigación Pro Fondecyt Regular 2023 (Resolución VRIDT N°076/2023), Resolución VRIDT N°09/2024 and Agencia Nacional de Investigación y Desarrollo (ANID) through Proyecto Fondecyt Regular 2024, Folio 1240514, Etapa 2024. Genly Leon would like to express his gratitude towards faculty member Alan Coley and staff members Anna Maria Davis, Nora Amaro, Jeanne Clyburne, and Mark Monk for their warm hospitality during the implementation of the final details of the research in the Department of Mathematics and Statistics at Dalhousie University.

Data availability

No data was used or created for the research described in the article.

References

González et al. (2023) E. González, K. Jusufi, G. Leon, and E. N. Saridakis, Phys. Dark Univ. 42, 101304 (2023), arXiv:2305.14305 [astro-ph.CO] .
Bond and Efstathiou (1984) J. R. Bond and G. Efstathiou, Astrophys. J. Lett. 285, L45 (1984).
Trimble (1987) V. Trimble, Ann. Rev. Astron. Astrophys. 25, 425 (1987).
Carroll et al. (1992) S. M. Carroll, W. H. Press, and E. L. Turner, Ann. Rev. Astron. Astrophys. 30, 499 (1992).
Perlmutter et al. (1998) S. Perlmutter, G. Aldering, M. D. Valle, S. Deustua, R. Ellis, S. Fabbro, A. Fruchter, G. Goldhaber, D. Groom, I. Hook, et al. (Supernova Cosmology Project), Nature 391, 51 (1998), arXiv:astro-ph/9712212 .
Riess et al. (1998) A. G. Riess, A. V. Filippenko, and P. Challis (Supernova Search Team), Astron. J. 116, 1009 (1998), arXiv:astro-ph/9805201 .
Perlmutter et al. (1999) S. Perlmutter, G. Aldering, G. Goldhaber, R. Knop, P. Nugent, P. G. Castro, S. Deustua, S. Fabbro, A. Goobar, D. E. Groom, et al. (Supernova Cosmology Project), Astrophys. J. 517, 565 (1999), arXiv:astro-ph/9812133 .
Aghanim et al. (2020) N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. Banday, R. Barreiro, N. Bartolo, S. Basak, et al. (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
Guth (1981) A. H. Guth, Phys. Rev. D 23, 347 (1981).
Ratra and Peebles (1988) B. Ratra and P. J. E. Peebles, Phys. Rev. D 37, 3406 (1988).
Parsons and Barrow (1995) P. Parsons and J. D. Barrow, Class. Quant. Grav. 12, 1715 (1995).
Rubano and Barrow (2001) C. Rubano and J. D. Barrow, Phys. Rev. D 64, 127301 (2001), arXiv:gr-qc/0105037 .
Saridakis (2009) E. N. Saridakis, Phys. Lett. B 676, 7 (2009), arXiv:0811.1333 [hep-th] .
Cai et al. (2010) Y.-F. Cai, E. N. Saridakis, M. R. Setare, and J.-Q. Xia, Phys. Rept. 493, 1 (2010), arXiv:0909.2776 [hep-th] .
Wali Hossain et al. (2015) M. Wali Hossain, R. Myrzakulov, M. Sami, and E. N. Saridakis, Int. J. Mod. Phys. D 24, 1530014 (2015), arXiv:1410.6100 [gr-qc] .
Barrow and Paliathanasis (2016) J. D. Barrow and A. Paliathanasis, Phys. Rev. D 94, 083518 (2016), arXiv:1609.01126 [gr-qc] .
Elizalde et al. (2004) E. Elizalde, S. Nojiri, and S. D. Odintsov, Phys. Rev. D 70, 043539 (2004), arXiv:hep-th/0405034 .
Elizalde et al. (2008) E. Elizalde, S. Nojiri, S. D. Odintsov, D. Saez-Gomez, and V. Faraoni, Phys. Rev. D 77, 106005 (2008), arXiv:0803.1311 [hep-th] .
Skugoreva et al. (2015) M. A. Skugoreva, E. N. Saridakis, and A. V. Toporensky, Phys. Rev. D 91, 044023 (2015), arXiv:1412.1502 [gr-qc] .
Saridakis and Tsoukalas (2016) E. N. Saridakis and M. Tsoukalas, JCAP 04, 017 (2016), arXiv:1602.06890 [gr-qc] .
Paliathanasis (2019) A. Paliathanasis, Gen. Rel. Grav. 51, 101 (2019), arXiv:1907.12261 [gr-qc] .
Banerjee et al. (2023) S. Banerjee, M. Petronikolou, and E. N. Saridakis, Phys. Rev. D 108, 024012 (2023), arXiv:2209.02426 [gr-qc] .
Santos et al. (2023) F. F. Santos, B. Pourhassan, and E. Saridakis, “de Sitter versus anti-de Sitter in Horndeski-like gravity,” (2023), arXiv:2305.05794 [hep-th] .
Leon et al. (2022) G. Leon, A. Paliathanasis, E. N. Saridakis, and S. Basilakos, Phys. Rev. D 106, 024055 (2022), arXiv:2203.14866 [gr-qc] .
Akrami et al. (2021) Y. Akrami et al. (CANTATA), Modified Gravity and Cosmology: An Update by the CANTATA Network, edited by E. N. Saridakis, R. Lazkoz, V. Salzano, P. Vargas Moniz, S. Capozziello, J. Beltrán Jiménez, M. De Laurentis, and G. J. Olmo (Springer, 2021) arXiv:2105.12582 [gr-qc] .
Leon and Saridakis (2009) G. Leon and E. N. Saridakis, JCAP 11, 006 (2009), arXiv:0909.3571 [hep-th] .
De Felice and Tsujikawa (2010) A. De Felice and S. Tsujikawa, Living Rev. Rel. 13, 3 (2010), arXiv:1002.4928 [gr-qc] .
Clifton et al. (2012) T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Phys. Rept. 513, 1 (2012), arXiv:1106.2476 [astro-ph.CO] .
Capozziello and De Laurentis (2011) S. Capozziello and M. De Laurentis, Phys. Rept. 509, 167 (2011), arXiv:1108.6266 [gr-qc] .
De Felice and Tsujikawa (2012) A. De Felice and S. Tsujikawa, JCAP 02, 007 (2012), arXiv:1110.3878 [gr-qc] .
Xu et al. (2012) C. Xu, E. N. Saridakis, and G. Leon, JCAP 07, 005 (2012), arXiv:1202.3781 [gr-qc] .
Bamba et al. (2012) K. Bamba, S. Capozziello, S. Nojiri, and S. D. Odintsov, Astrophys. Space Sci. 342, 155 (2012), arXiv:1205.3421 [gr-qc] .
Leon and Saridakis (2013) G. Leon and E. N. Saridakis, JCAP 03, 025 (2013), arXiv:1211.3088 [astro-ph.CO] .
Kofinas et al. (2014) G. Kofinas, G. Leon, and E. N. Saridakis, Class. Quant. Grav. 31, 175011 (2014), arXiv:1404.7100 [gr-qc] .
Bahamonde et al. (2015) S. Bahamonde, C. G. Böhmer, and M. Wright, Phys. Rev. D 92, 104042 (2015), arXiv:1508.05120 [gr-qc] .
Momeni and Myrzakulov (2015) D. Momeni and R. Myrzakulov, Astrophys. Space Sci. 360, 28 (2015), arXiv:1511.01205 [physics.gen-ph] .
Cai et al. (2016) Y.-F. Cai, S. Capozziello, M. De Laurentis, and E. N. Saridakis, Rept. Prog. Phys. 79, 106901 (2016), arXiv:1511.07586 [gr-qc] .
Krssak et al. (2019) M. Krssak, R. J. van den Hoogen, J. G. Pereira, C. G. Böhmer, and A. A. Coley, Class. Quant. Grav. 36, 183001 (2019), arXiv:1810.12932 [gr-qc] .
Dehghani et al. (2023) A. Dehghani, B. Pourhassan, S. Zarepour, and E. N. Saridakis, “Thermodynamic schemes of charged BTZ-like black holes in arbitrary dimensions,” (2023), arXiv:2305.08219 [hep-th] .
Salucci (2018) P. Salucci, Found. Phys. 48, 1517 (2018), arXiv:1807.08541 [astro-ph.GA] .
Milgrom (1983) M. Milgrom, Astrophys. J. 270, 365 (1983).
Ferreira and Starkmann (2009) P. G. Ferreira and G. Starkmann, Science 326, 812 (2009), arXiv:0911.1212 [astro-ph.CO] .
Milgrom and Sanders (2003) M. Milgrom and R. H. Sanders, The Astrophysical Journal 599, L25 (2003).
Tiret et al. (2007) O. Tiret, F. Combes, G. W. Angus, B. Famaey, and H. Zhao, Astron. Astrophys. 476, L1 (2007), arXiv:0710.4070 [astro-ph] .
Kroupa et al. (2010) P. Kroupa, B. Famaey, K. S. de Boer, J. Dabringhausen, M. S. Pawlowski, C. M. Boily, H. Jerjen, D. Forbes, G. Hensler, and M. Metz, Astron. Astrophys. 523, A32 (2010), arXiv:1006.1647 [astro-ph.CO] .
Cardone et al. (2011) V. F. Cardone, G. Angus, A. Diaferio, C. Tortora, and R. Molinaro, Mon. Not. Roy. Astron. Soc. 412, 2617 (2011), arXiv:1011.5741 [astro-ph.CO] .
Richtler et al. (2011) T. Richtler, R. Salinas, I. Misgeld, M. Hilker, G. K. T. Hau, A. J. Romanowsky, Y. Schuberth, and M. Spolaor, Astron. Astrophys. 531, A119 (2011), arXiv:1103.2053 [astro-ph.CO] .
Berezhiani and Khoury (2015) L. Berezhiani and J. Khoury, Phys. Rev. D 92, 103510 (2015), arXiv:1507.01019 [astro-ph.CO] .
Boehmer and Harko (2007) C. G. Boehmer and T. Harko, JCAP 06, 025 (2007), arXiv:0705.4158 [astro-ph] .
Takahashi (2004) R. Takahashi, J. Korean Phys. Soc. 45, S1808 (2004), arXiv:astro-ph/0405099 .
Hioki and Maeda (2009) K. Hioki and K.-i. Maeda, Phys. Rev. D 80, 024042 (2009), arXiv:0904.3575 [astro-ph.HE] .
Brito et al. (2015) R. Brito, V. Cardoso, and P. Pani, Lect. Notes Phys. 906, pp.1 (2015), arXiv:1501.06570 [gr-qc] .
Cunha et al. (2015) P. V. P. Cunha, C. A. R. Herdeiro, E. Radu, and H. F. Runarsson, Phys. Rev. Lett. 115, 211102 (2015), arXiv:1509.00021 [gr-qc] .
Ohgami and Sakai (2015) T. Ohgami and N. Sakai, Phys. Rev. D 91, 124020 (2015), arXiv:1704.07065 [gr-qc] .
Moffat (2015) J. W. Moffat, Eur. Phys. J. C 75, 130 (2015), arXiv:1502.01677 [gr-qc] .
Abdujabbarov et al. (2016) A. Abdujabbarov, M. Amir, B. Ahmedov, and S. G. Ghosh, Phys. Rev. D 93, 104004 (2016), arXiv:1604.03809 [gr-qc] .
Cunha and Herdeiro (2018) P. V. P. Cunha and C. A. R. Herdeiro, Gen. Rel. Grav. 50, 42 (2018), arXiv:1801.00860 [gr-qc] .
Mizuno et al. (2018) Y. Mizuno, Z. Younsi, C. M. Fromm, O. Porth, M. De Laurentis, H. Olivares, H. Falcke, M. Kramer, and L. Rezzolla, Nature Astron. 2, 585 (2018), arXiv:1804.05812 [astro-ph.GA] .
Tsukamoto (2018) N. Tsukamoto, Phys. Rev. D 97, 064021 (2018), arXiv:1708.07427 [gr-qc] .
Psaltis (2019) D. Psaltis, Gen. Rel. Grav. 51, 137 (2019), arXiv:1806.09740 [astro-ph.HE] .
Amir et al. (2019) M. Amir, K. Jusufi, A. Banerjee, and S. Hansraj, Class. Quant. Grav. 36, 215007 (2019), arXiv:1806.07782 [gr-qc] .
Gralla et al. (2019) S. E. Gralla, D. E. Holz, and R. M. Wald, Phys. Rev. D 100, 024018 (2019), arXiv:1906.00873 [astro-ph.HE] .
Bambi et al. (2019) C. Bambi, K. Freese, S. Vagnozzi, and L. Visinelli, Phys. Rev. D 100, 044057 (2019), arXiv:1904.12983 [gr-qc] .
Cunha et al. (2019) P. V. P. Cunha, C. A. R. Herdeiro, and E. Radu, Universe 5, 220 (2019), arXiv:1909.08039 [gr-qc] .
Khodadi et al. (2020) M. Khodadi, A. Allahyari, S. Vagnozzi, and D. F. Mota, JCAP 09, 026 (2020), arXiv:2005.05992 [gr-qc] .
Perlick and Tsupko (2022) V. Perlick and O. Y. Tsupko, Phys. Rept. 947, 1 (2022), arXiv:2105.07101 [gr-qc] .
Vagnozzi et al. (2023) S. Vagnozzi et al., Class. Quant. Grav. 40, 165007 (2023), arXiv:2205.07787 [gr-qc] .
Saurabh and Jusufi (2021) K. Saurabh and K. Jusufi, Eur. Phys. J. C 81, 490 (2021), arXiv:2009.10599 [gr-qc] .
Jusufi et al. (2020) K. Jusufi, M. Jamil, and T. Zhu, Eur. Phys. J. C 80, 354 (2020), arXiv:2005.05299 [gr-qc] .
Tsupko et al. (2020) O. Y. Tsupko, Z. Fan, and G. S. Bisnovatyi-Kogan, Class. Quant. Grav. 37, 065016 (2020), arXiv:1905.10509 [gr-qc] .
Akiyama et al. (2019) K. Akiyama et al. (Event Horizon Telescope), Astrophys. J. Lett. 875, L1 (2019), arXiv:1906.11238 [astro-ph.GA] .
Psaltis et al. (2020) D. Psaltis et al. (Event Horizon Telescope), Phys. Rev. Lett. 125, 141104 (2020), arXiv:2010.01055 [gr-qc] .
Akiyama et al. (2021) K. Akiyama et al. (Event Horizon Telescope), Astrophys. J. Lett. 910, L13 (2021), arXiv:2105.01173 [astro-ph.HE] .
Akiyama et al. (2022) K. Akiyama et al. (Event Horizon Telescope), Astrophys. J. Lett. 930, L12 (2022).
Jusufi et al. (2023) K. Jusufi, G. Leon, and A. D. Millano, Phys. Dark Univ. 42, 101318 (2023), arXiv:2304.11492 [gr-qc] .
Loeb and Weiner (2011) A. Loeb and N. Weiner, Phys. Rev. Lett. 106, 171302 (2011), arXiv:1011.6374 [astro-ph.CO] .
Anderson et al. (2017) D. E. Anderson, E. A. Bergin, G. A. Blake, F. J. Ciesla, R. Visser, and J.-E. Lee, The Astrophysical Journal 845, 13 (2017).
Kanzi et al. (2020) S. Kanzi, S. H. Mazharimousavi, and I. Sakallı, Annals Phys. 422, 168301 (2020), arXiv:2007.05814 [hep-th] .
Cedeño et al. (2017) F. X. L. Cedeño, A. X. González-Morales, and L. A. Ureña López, Phys. Rev. D 96, 061301 (2017), arXiv:1703.10180 [gr-qc] .
Berezhiani et al. (2009) Z. Berezhiani, F. Nesti, L. Pilo, and N. Rossi, JHEP 07, 083 (2009), arXiv:0902.0144 [hep-th] .
Sakallı et al. (2018) I. Sakallı, K. Jusufi, and A. Övgün, Gen. Rel. Grav. 50, 125 (2018), arXiv:1803.10583 [gr-qc] .
Nicolini et al. (2019) P. Nicolini, E. Spallucci, and M. F. Wondrak, Phys. Lett. B 797, 134888 (2019), arXiv:1902.11242 [gr-qc] .
Laurentis and Salucci (2022) M. D. Laurentis and P. Salucci, The Astrophysical Journal 929, 17 (2022).
de Rham et al. (2011) C. de Rham, G. Gabadadze, and A. J. Tolley, Phys. Rev. Lett. 106, 231101 (2011), arXiv:1011.1232 [hep-th] .
Ghosh et al. (2016) S. G. Ghosh, L. Tannukij, and P. Wongjun, Eur. Phys. J. C 76, 119 (2016), arXiv:1506.07119 [gr-qc] .
Simpson and Visser (2019) A. Simpson and M. Visser, JCAP 02, 042 (2019), arXiv:1812.07114 [gr-qc] .
Visser (1995) M. Visser, Lorentzian Wormholes: From Einstein to Hawking, Computational and Mathematical Physics (American Inst. of Physics, 1995).
Visser (2020) M. Visser, Class. Quant. Grav. 37, 045001 (2020), arXiv:1908.11058 [gr-qc] .
Boonserm et al. (2020) P. Boonserm, T. Ngampitipan, A. Simpson, and M. Visser, Phys. Rev. D 101, 024022 (2020), arXiv:1910.08008 [gr-qc] .
Fenichel (1979) N. Fenichel, Journal of Differential Equations 31, 53 (1979).
Fusco and Hale (1989) G. Fusco and J. K. Hale, Journal of Dynamics and Differential Equations 1, 75 (1989).
Berglund and Gentz (2006) N. Berglund and B. Gentz, Noise-Induced Phenomena in Slow-Fast Dynamical Systems (Springer-Verlag, 2006).
Verhulst (2006) F. Verhulst, Methods and Applications of Singular Perturbations: Boundary Layers and Multiple Timescale Dynamics, Vol. 50 (Springer Science & Business Media, 2006).
Holmes (2012) M. Holmes, Introduction to Perturbation Methods, Texts in Applied Mathematics (Springer New York, 2012).
Kevorkian and Cole (2013) J. Kevorkian and J. Cole, Perturbation Methods in Applied Mathematics, Applied Mathematical Sciences (Springer New York, 2013).
Awrejcewicz (2014) J. Awrejcewicz, Applied Non-Linear Dynamical Systems, Springer Proceedings in Mathematics & Statistics, Vol. 93 (Springer International Publishing AG, Cham, 2014).
Hawking (1975a) S. W. Hawking, in Euclidean quantum gravity (World Scientific, 1975) pp. 167–188.
Hawking (1975b) S. W. Hawking, Commun. Math. Phys. 43, 199 (1975b), [Erratum: Commun.Math.Phys. 46, 206 (1976)].
Filho et al. (2023a) A. A. A. Filho, S. Zare, P. J. Porfírio, J. Kříž, and H. Hassanabadi, Phys. Lett. B 838, 137744 (2023a), arXiv:2212.03134 [hep-th] .
Sedaghatnia et al. (2023) P. Sedaghatnia, H. Hassanabadi, A. A. A. Filho, P. J. P. Porfírio, and W. S. Chung, “Thermodynamical properties of a deformed Schwarzschild black hole via Dunkl generalization,” (2023), arXiv:2302.11460 [gr-qc] .
Filho et al. (2023b) A. A. A. Filho, J. Furtado, H. Hassanabadi, and J. A. A. S. Reis, “Thermal analysis of photon-like particles in rainbow gravity,” (2023b), arXiv:2305.08587 [gr-qc] .
Filho et al. (2023c) A. A. A. Filho, J. Furtado, J. A. A. S. Reis, and J. E. G. Silva, “Thermodynamical properties of an ideal gas in a traversable wormhole,” (2023c), arXiv:2302.05492 [gr-qc] .
Ma and Zhao (2014) M.-S. Ma and R. Zhao, Class. Quant. Grav. 31, 245014 (2014), arXiv:1411.0833 [gr-qc] .
Maluf and Neves (2018) R. V. Maluf and J. C. S. Neves, Phys. Rev. D 97, 104015 (2018), arXiv:1801.02661 [gr-qc] .
Iyer and Will (1987) S. Iyer and C. M. Will, Phys. Rev. D 35, 3621 (1987).
Iyer (1987) S. Iyer, Phys. Rev. D 35, 3632 (1987).
Konoplya (2003) R. A. Konoplya, Phys. Rev. D 68, 024018 (2003), arXiv:gr-qc/0303052 .
Matyjasek and Opala (2017) J. Matyjasek and M. Opala, Phys. Rev. D 96, 024011 (2017), arXiv:1704.00361 [gr-qc] .
Schutz and Will (1985) B. F. Schutz and C. M. Will, APJL 291, L33 (1985).
Konoplya (2004) R. A. Konoplya, J. Phys. Stud. 8, 93 (2004).
Santos et al. (2016) V. Santos, R. V. Maluf, and C. A. S. Almeida, Phys. Rev. D 93, 084047 (2016), arXiv:1509.04306 [gr-qc] .
Konoplya and Zhidenko (2011) R. A. Konoplya and A. Zhidenko, Rev. Mod. Phys. 83, 793 (2011), arXiv:1102.4014 [gr-qc] .
Berti et al. (2009) E. Berti, V. Cardoso, and A. O. Starinets, Class. Quant. Grav. 26, 163001 (2009), arXiv:0905.2975 [gr-qc] .
Heidari et al. (2023a) N. Heidari, H. Hassanabadi, A. A. A. Filho, J. Kuríuz, S. Zare, and P. J. Porfírio, “Gravitational signatures of a non–commutative stable black hole,” (2023a), arXiv:2305.06838 [gr-qc] .
Chen et al. (2023) H. Chen, T. Sathiyaraj, H. Hassanabadi, Y. Yang, Z. W. Long, and F. Q. Tu, Indian Journal of Physics (2023), 10.1007/s12648-023-02734-8.
Toshmatov et al. (2017) B. Toshmatov, C. Bambi, B. Ahmedov, Z. Stuchlík, and J. Schee, Phys. Rev. D 96, 064028 (2017), arXiv:1705.03654 [gr-qc] .
Kim (2004) S.-W. Kim, “Gravitational perturbation of traversable wormhole,” (2004), arXiv:gr-qc/0401007 .
Filho (2023) A. A. A. Filho, “Implications of a Simpson-Visser solution in Verlinde’s framework,” (2023), arXiv:2308.04939 [gr-qc] .
Tsukamoto (2021) N. Tsukamoto, Phys. Rev. D 104, 064022 (2021), arXiv:2105.14336 [gr-qc] .
Tsukamoto (2022) N. Tsukamoto, Phys. Rev. D 105, 084036 (2022), arXiv:2202.09641 [gr-qc] .
Guerrero et al. (2022) M. Guerrero, G. J. Olmo, D. Rubiera-Garcia, and D. Sáez-Chillón Gómez, Phys. Rev. D 106, 044070 (2022), arXiv:2205.12147 [gr-qc] .
Singh and Ghosh (2018) B. P. Singh and S. G. Ghosh, Annals Phys. 395, 127 (2018), arXiv:1707.07125 [gr-qc] .
Filho et al. (2023d) A. A. A. Filho, J. A. A. S. Reis, and H. Hassanabadi, “Exploring antisymmetric tensor effects on black hole shadows and quasinormal frequencies,” (2023d), arXiv:2309.15778 [gr-qc] .
Heidari et al. (2023b) N. Heidari, H. Hassanabadi, A. A. A. Filho, and J. Kuríuz, “Exploring Non–commutativity as a Perturbation in the Schwarzschild Black Hole: Quasinormal Modes, Scattering, and Shadows,” (2023b), arXiv:2308.03284 [gr-qc] .
Filho et al. (2023e) A. A. A. Filho, H. Hassanabadi, N. Heidari, J. Kríz, P. J. Porfírio, and S. Zare, “Gravitational traces of bumblebee gravity in metric-affine formalism,” (2023e), arXiv:2305.18871 [gr-qc] .
Bambi (2013) C. Bambi, Phys. Rev. D 87, 107501 (2013), arXiv:1304.5691 [gr-qc] .
Jusufi et al. (2024) K. Jusufi, E. González, and G. Leon, “Addressing the Hubble tension in Yukawa cosmology?” (2024), arXiv:2402.02512 [astro-ph.CO] .
Jusufi and Sheykhi (2024) K. Jusufi and A. Sheykhi, “Apparent Dark Matter Inspired by Einstein Equation of State,” (2024), arXiv:2402.00785 [gr-qc] .
Allahyari et al. (2020) A. Allahyari, M. Khodadi, S. Vagnozzi, and D. F. Mota, JCAP 02, 003 (2020), arXiv:1912.08231 [gr-qc] .
Goebel (1972) C. J. Goebel, APJL 172, L95 (1972).
Cardoso et al. (2009) V. Cardoso, A. S. Miranda, E. Berti, H. Witek, and V. T. Zanchin, Phys. Rev. D 79, 064016 (2009), arXiv:0812.1806 [hep-th] .
Stefanov et al. (2010) I. Z. Stefanov, S. S. Yazadjiev, and G. G. Gyulchev, Phys. Rev. Lett. 104, 251103 (2010), arXiv:1003.1609 [gr-qc] .
Jusufi (2020) K. Jusufi, Phys. Rev. D 101, 084055 (2020), arXiv:1912.13320 [gr-qc] .
Cuadros-Melgar et al. (2020) B. Cuadros-Melgar, R. D. B. Fontana, and J. de Oliveira, Phys. Lett. B 811, 135966 (2020), arXiv:2005.09761 [gr-qc] .
Konoplya and Stuchlík (2017) R. A. Konoplya and Z. Stuchlík, Phys. Lett. B 771, 597 (2017), arXiv:1705.05928 [gr-qc] .
Konoplya (2023) R. A. Konoplya, Phys. Lett. B 838, 137674 (2023), arXiv:2210.08373 [gr-qc] .

Dark matter signatures of black holes with Yukawa potential

Abstract

keywords:

1 Introduction

2 Black hole solution with Yukawa potential

2.1 Correspondence with the black hole solution in dRGT gravity

2.2 Energy conditions

3 Accuracy of the approximate solution

3.1 Inner solution

3.2 Outer solution

3.3 Multiple scales - Derivative expansion method

4 Thermodynamic properties

4.1 Hawking temperature

4.2 Entropy

4.3 Heat capacity

5 Quasinormal frequencies

5.1 Scalar perturbations

5.2 Electromagnetic perturbations

5.3 Tensorial perturbations

6 Geodesic trajectories

6.1 The complete case

6.2 The approximated case

7 Photon spheres and Shadows: the case of optically thin radiating and infalling gas surrounding a black hole

7.1 Can we constrain α\alpha?

7.2 The galactic flat rotating curves

8 Eikonal quasinormal modes and shadow radius correspondence

9 Conclusion

CRediT authorship contribution statement

Declaration of competing interest

Acknowledgements

Data availability

References

7.1 Can we constrain $\alpha$ ?