# Electrically charged regular black holes in nonlinear electrodynamics:
light rings, shadows and gravitational lensing

Marco A. A. de Paula<EMAIL_ADDRESS>Programa de Pós-Graduação em
Física, Universidade Federal do Pará, 66075-110, Belém, Pará, Brazil.
Consortium for Fundamental Physics, School of Mathematics and Statistics,
University of Sheffield, Hicks Building, Hounsfield Road, Sheffield S3 7RH,
United Kingdom.  Haroldo C. D. Lima Junior<EMAIL_ADDRESS>Programa de
Pós-Graduação em Física, Universidade Federal do Pará, 66075-110, Belém, Pará,
Brazil.  Departamento de Matemática da Universidade de Aveiro and Centre for
Research and Development in Mathematics and Applications (CIDMA), Campus de
Santiago, 3810-183 Aveiro, Portugal.  Pedro V. P. Cunha<EMAIL_ADDRESS>Departamento de Matemática da Universidade de Aveiro and Centre for Research
and Development in Mathematics and Applications (CIDMA), Campus de Santiago,
3810-183 Aveiro, Portugal.  Luís C. B. Crispino<EMAIL_ADDRESS>Programa de
Pós-Graduação em Física, Universidade Federal do Pará, 66075-110, Belém, Pará,
Brazil.  Departamento de Matemática da Universidade de Aveiro and Centre for
Research and Development in Mathematics and Applications (CIDMA), Campus de
Santiago, 3810-183 Aveiro, Portugal.

###### Abstract

Within nonlinear electrodynamics (NED), photons follow null geodesics of an
effective geometry, which is different from the geometry of the spacetime
itself. Over the last years, several works were dedicated to investigate the
motion of photons in the effective geometry of NED-based magnetically charged
regular black hole (RBH) solutions. However, there are few works considering
electrically charged RBHs. We study the light rings, shadows, and
gravitational lensing of the electrically charged RBH solution proposed by
Irina Dymnikova (ID), which is a static and spherically symmetric spacetime
with a NED source. We show that the shadow associated to the effective
geometry can be almost 10% bigger that the one associated to the standard
geometry. We also find that the ID solution may mimic the shadow properties of
the Reissner-Nordström (RN) BH, for low-to-extreme values of the electric
charge. Besides that, by using the backwards ray-tracing technique, we obtain
that ID and RN BH solutions can have a very similar gravitational lensing, for
some values of the correspondent electric charges. We also show that the
motion of photons in the effective geometry can be interpreted as a non-
geodesic curve submitted to a 4-force term, from the perspective of an
observer in the standard geometry.

##  I Introduction

General Relativity (GR) is a well-established classical gravitational theory
Will (2014); Abbott  _et al._ (2016); Akiyama  _et al._ (2019). Although it
has accumulated remarkable and numerous triumphs, it presents limitations,
specially at the core of the standard black hole (BH) solutions. GR predicts
the existence of singularities, which are pathologies where the laws of
physics break down Hawking and Ellis (1973), challenging the validity of
Einstein’s theory.

A possibility to overcome such pathologies is to consider appropriated
distributions of matter, leading to singularity-free BH solutions within GR.
The first line element for a non-singular BH geometry was proposed by James
Bardeen in 1968 Bardeen (1968). By minimally coupling GR and nonlinear
electrodynamics (NED), it was shown that it is possible to obtain various
exact charged RBH solutions (cf. Refs. Dymnikova (1992); Ayón-Beato and García
(1998, 1999a, 1999b)). In these theories, the Bardeen geometry can be
interpreted as a RBH sourced by a nonlinear magnetic Ayón-Beato and García
(2000) or electric monopole Rodrigues and de S. Silva (2018).

NED models can be seen as possible ultraviolet completions of linear
electrodynamics, i.e., for electromagnetic fields with magnitudes approaching
Greiner and Reinhardt (2008):

$E_{\rm{cri}}=1.3\times 10^{18}\,\dfrac{V}{m}\ \ \ \text{and}\ \ \
B_{\rm{cri}}=4.4\times 10^{9}\,T.$ (1)

One of the first covariant models of NED was proposed in 1934 (the so-called
Born-Infeld electrodynamics) as an attempt to obtain a finite self-energy
density for the electric charge Born (1934); Born and Infeld (1934). Another
influential model of NED is the Euler-Heisenberg theory Heisenberg and Euler
(1936), which is related with two important predictions of Quantum
Electrodynamics (QED): the light-by-light scattering Aaboud  _et al._ (2017);
Aad  _et al._ (2019) and the vacuum birefringence Mignani  _et al._ (2017);
Adam  _et al._ (2021). Beyond BH physics (see also Refs. Bronnikov (2001);
Matyjasek (2004); Dymnikova (2004); Balart and Vagenas (2014); Ma (2015);
Junior  _et al._ (2015); Fan and Wang (2016); Kruglov (2017); Toshmatov  _et
al._ (2018); de S. Silva and Rodrigues (2018)) and QED, NED has also
applications in string/M-theories Fradkin and Tseytlin (1985); Seiberg and
Witten (1999); Tseytlin (2000); Gibbons and Herdeiro (2001) and cosmology
García-Salcedo and Bretón (2000); De Lorenci  _et al._ (2002); Novello  _et
al._ (2004); Campanelli  _et al._ (2008). Among the applications of NED in BH
physics, one important result is that the motion of photons can be interpreted
as a null geodesic of an effective geometry Plebanski (1970); Boillat (1970);
Gutiérrez  _et al._ (1981); Novello  _et al._ (2000), which is different from
the geometry of the spacetime itself.

Since NED affects the motion of photons, the analysis of light rings (LRs),
shadows, and gravitational lensing – which are of utmost importance within the
context of BH physics – requires special attention. The LRs are circular
photon orbits that can be studied by analyzing the null geodesics in a given
(effective) geometry, as it was done for some NED-based magnetically charged
RBHs Vrba  _et al._ (2019); Habibina and Ramadhan (2020); Amaro and Macías
(2020); Habibina  _et al._ (2021). Noticeably, the analysis of the null
geodesics alone is not enough to distinguish the type of charge of a BH in the
same NED theory Toshmatov  _et al._ (2021). Besides that, in the
electromagnetic channel, the LRs are closely related to the BH shadow Cunha
and Herdeiro (2018), as seen by a distant observer. The BH shadow is related
to the dark region formed when a BH is illuminated by some source of light,
for instance, an accretion disk that surrounds the BH Akiyama  _et al._
(2019). Recently, some works studying the shadows of NED-based RBHs,
considering the effective geometry, were performed Stuchlík and Schee (2019);
Stuchlík  _et al._ (2019); Allahyari  _et al._ (2020); Kruglov (2020);
Rayimbaev  _et al._ (2020); Hu  _et al._ (2021), but focusing mainly on
magnetically charged solutions.

The study of the deflection of a light ray by a compact object due to the
gravitational interaction plays an important role in Einstein’s theory. For
instance, the first confirmed prediction of GR, the deflection of light by the
Sun Dyson  _et al._ (1920); Crispino and Kennefick (2019), is an example of
gravitational lensing effect. Over the last decades, several works on
gravitational lensing in standard BH spacetimes have been done (see, e.g.,
Refs. Virbhadra and Ellis (2000); Bozza (2002); Eiroa  _et al._ (2002); Keeton
and Petters (2005); Gibbons and Werner (2008); Bozza (2010); Arakida (2018);
Li and Zhou (2020); Pantig and Rodulfo (2020) and references therein). In the
background of NED-based RBH solutions, considering the effective geometry, the
gravitational lensing was studied for electric and magnetic models Eiroa
(2006); Liang (2017a, b); Ghaffarnejad  _et al._ (2018); Schee and Stuchlík
(2019). However, the computation of gravitational lensing using backwards ray-
tracing techniques Bohn  _et al._ (2015); Cunha and Herdeiro (2018) has not
been performed so far in the background of electrically charged NED-based
RBHs.

It is also important to emphasize that, within NED, electrically charged RBHs
are, in general, derived in the so-called $P$ framework Bronnikov (2001). In
this framework, the electric models could exist, from the theoretical point of
view, if they satisfy the weak energy condition Dymnikova (2004). The weak
energy condition leads to a de-Sitter behavior at the core of the central
object, providing a regular center, and the Maxwell limit can be satisfied at
infinity, which is the case, e.g., for the solutions in Refs. Ayón-Beato and
García (1998); Dymnikova (2004).

Although it is widely believed that astrophysical BHs are essentially neutral,
it has been argued that (at least) a small non-zero electric charge is
possible Zajaček  _et al._ (2019, 2018); Zajaček and Tursunov (2019), which
can affect the motion of charged particles. Therefore the study of
electrically charged BHs, in the spherically symmetric case, can be useful not
only to improve our theoretical understanding of BH physics, but also to gauge
the role of NED and its hypothetical impact in the context of astrophysical
BHs.

The aim of this work is to study the imprints of NED in the trajectories of
the photons by analyzing the LRs, shadows and gravitational lensing. For
concreteness, we focus on the static and spherically symmetric electrically
charged NED-based RBH solution proposed by Irina Dymnikova (ID) Dymnikova
(2004). Since the casual structure of the ID solution is similar to the
Reissner-Nordström (RN) one, we compare our results to those obtained in the
RN geometry. The remainder of this paper is organized as follows. In Sec. II
we review the ID geometry. The null geodesic equations, considering the
standard and effective geometries, are studied in Sec. III. Our main results
are presented in Sec. IV, and our final remarks in Sec. V. Throughout this
paper we use the natural units, for which $G=c=\hbar=1$, and the metric
signature ($+,-,-,-$).

##  II Background

In the $F$ framework, the action that describes NED minimally coupled with
gravity can be written as Bronnikov (2001)

$\mathrm{S}=\dfrac{1}{16\pi}\int
d^{4}x\sqrt{-g}\left[R-\mathcal{L}(F)\right],$ (2)

where $g$ is the determinant of the metric tensor $g_{\mu\nu}$, $R$ is the
corresponding Ricci scalar, and $\mathcal{L}(F)$ is a gauge-invariant
electromagnetic Lagrangian density. The function $F\equiv
F_{\mu\nu}F^{\mu\nu}$ is the Maxwell scalar, with $F_{\mu\nu}$ being the
standard electromagnetic field tensor. By introducing a structural function
$\mathcal{H}(P)$ through a Legendre transformation Salazar  _et al._ (1987),
namely

$\mathcal{H}(P)=2F\mathcal{L}_{F}-\mathcal{L}(F),$ (3)

one can obtain an alternative form for the NED theory in the so-called $P$
framework Bronnikov (2001). Within this context, the function $P\equiv
P_{\mu\nu}P^{\mu\nu}$ is a scalar obtained from the auxiliary anti-symmetric
tensor $P_{\mu\nu}$, defined as $P_{\mu\nu}\equiv\mathcal{L}_{F}F_{\mu\nu}$,
where $\mathcal{L}_{F}\equiv\partial\mathcal{L}/\partial F$. The relations
between the $F$ and $P$ frameworks are given by (see, for instance, Ref.
Dymnikova (2004)):

$P=(\mathcal{L}_{F})^{2}F,\ \ \ \mathcal{H}_{P}\mathcal{L}_{F}=1,\ \ \
\text{and}\ \ \ F_{\mu\nu}=\mathcal{H}_{P}P_{\mu\nu},$ (4)

where $\mathcal{H}_{P}\equiv\partial\mathcal{H}/\partial P$. By using Eqs.
(2)-(4), we can write the corresponding action in the $P$ framework as

$\mathrm{S}=\dfrac{1}{16\pi}\int
d^{4}x\sqrt{-g}\left[R-\left(2P\mathcal{H}_{P}-\mathcal{H}(P)\right)\right].$
(5)

The corresponding field equations are given by

$G^{\mu}_{\ \ \nu}=-T^{\mu}_{\ \
\nu}=\dfrac{1}{2}\left[4\mathcal{H}_{P}P_{\nu\alpha}P^{\mu\alpha}-\delta^{\mu}_{\
\ \nu}\left(2P\mathcal{H}_{P}-\mathcal{H}\right)\right],$ (6)

which are the Einstein-NED (E-NED) field equations written in the $P$
framework. The conservation equation of $P_{\mu\nu}$ and the corresponding
Bianchi identities are given by

$\nabla_{\mu}P^{\mu\nu}=0\ \ \ \text{and}\ \ \
\nabla_{\mu}\left(\mathcal{H}_{P}\star P^{\mu\nu}\right)=0,$ (7)

respectively, where $\star$ is the Hodge symbol. A correspondence with
Maxwell’s theory is obtained if $\mathcal{H}(P)\rightarrow P$ and
$\mathcal{H}_{P}\rightarrow 1$, for small $P$. The $P$ framework is useful to
obtain exact solutions of Einstein field equations in the presence of NED
sources Salazar  _et al._ (1987) and it is equivalent to the $F$ framework
where the function $F(P)$ is a monotonic function of $P$ Bronnikov (2001).

Within the $P$ framework, NED-based RBHs may be found by specifying the NED
source $\mathcal{H}(P)$ and the appropriated function $P_{\mu\nu}$ 111BH
solutions obtained in the $P$ framework can also be formally derived in the
$F$ framework by using a suitable nonuniform variational method Dymnikova  _et
al._ (2015).. For the electrically charged RBH solution proposed by Irina
Dymnikova (ID) Dymnikova (2004), the NED source is specified by the following
structural function:

$\mathcal{H}(P)\equiv\dfrac{P}{(1+\alpha\sqrt{-P})^{2}},$ (8)

where $\alpha$ is a constant to be determined by the field equations.

To solve the field equations we need to take an ansatz for the line element
describing the spacetime. For the ID solution, it is considered a static and
spherically symmetric geometry, with the line element of the form

$ds^{2}=f(r)dt^{2}-f(r)^{-1}dr^{2}-r^{2}d\Omega^{2},$ (9)

in which $d\Omega^{2}=d\theta^{2}+\sin^{2}\theta d\varphi^{2}$ is the line
element of a unit 2-sphere and $f(r)$ is the metric function, given by

$f(r)=1-\dfrac{2\mathcal{M}(r)}{r}.$ (10)

The function $\mathcal{M}(r)$ is determined by the E-NED field equations. From
its asymptotic behavior it is possible to obtain the total mass
$\mathcal{M}(r\rightarrow\infty)=M$ of the (regular) BH Fan and Wang (2016);
Toshmatov  _et al._ (2018). Since we are considering a spherically symmetric
background and a purely electric NED source, the appropriated ansatz for
$P_{\mu\nu}$ can be written as

$P_{\mu\nu}=\left(\delta_{\mu}^{t}\delta_{\nu}^{r}-\delta_{\nu}^{t}\delta_{\mu}^{r}\right)D(r),$
(11)

where $D(r)$ is a function to be determined by the conservation equation of
$P_{\mu\nu}$ (7). Notice that the NED source (8) satisfies a correspondence
with Maxwell’s theory at infinity, which can be inferred by taking a series
expansion of the model around $P=0$. Taking this into account and integrating
Eq. (7), we obtain that $D(r)$ is given by

$D(r)=\dfrac{Q}{r^{2}},$ (12)

where $Q$ is the electric charge of the central object, then

$P=-\dfrac{2Q^{2}}{r^{4}}.$ (13)

The $G^{t}_{\ t}$ component of the E-NED field equations leads to

$\mathcal{M}(r)=-\dfrac{1}{4}\int_{0}^{r}\mathcal{H}(P)x^{2}dx,$ (14)

which, considering Eqs. (8) and (13), results in

$\mathcal{M}(r)=\dfrac{Q^{2}}{8}\left[\dfrac{2^{\frac{3}{4}}}{\sqrt{\alpha\,Q}}\arctan\left(\dfrac{r}{2^{\frac{1}{4}}\sqrt{\alpha\,Q}}\right)-\dfrac{2r}{\sqrt{2}\alpha\,Q+r^{2}}\right].$
(15)

The value of $\alpha$ can be fixed by recalling that the limit
$\mathcal{M}(r\rightarrow\infty)=M$ provides the unique mass of the BH, thus

$\alpha=\dfrac{\pi^{2}Q^{3}}{64\sqrt{2}M^{2}}.$ (16)

Since $\alpha$ is a model parameter, rather than an integration constant
(which is the case for $M,Q$) this means the choice of BH mass and charge,
fixes the model coupling $\alpha$. But the choice of charge to mass ratio
$Q/M$ does not fix $\alpha$.

Considering Eqs. (15) and (16), and defining

$z\equiv\dfrac{\pi Q^{2}}{8M},$ (17)

we obtain the metric function of the ID solution, given by Dymnikova (2004)

$f^{\rm{ID}}(r)=1-\dfrac{4M}{\pi
r}\left[\arctan\left(\dfrac{r}{z}\right)-\dfrac{rz}{r^{2}+z^{2}}\right].$ (18)

In the limit $r\rightarrow\infty$, $f^{\rm{ID}}(r)$ behaves as

$\left.f^{\rm{ID}}(r)\right|_{\infty}\approx
1-\dfrac{2M}{r}+\dfrac{Q^{2}}{r^{2}}-\dfrac{2Q^{2}z^{2}}{3r^{4}}+\mathcal{O}\left[\dfrac{1}{r^{5}}\right],$
(19)

which approaches the metric function of the RN spacetime, as expected since
the NED source associated with the ID solution satisfies a correspondence with
the linear electrodynamics in the weak field limit. Although the contributions
of order $r^{-n}$, with $n\geq 4$, in the metric function are negligible in
the weak field limit, they play an important role in the higher order
corrections of the weak deflection angle, as discussed in the Appendix A. On
the other hand, as we approach the core, the ID solution has a de Sitter
behavior, given by

$\left.f^{\rm{ID}}(r)\right|_{0}\approx
1-\dfrac{1}{3}\left(\dfrac{Q^{2}}{z^{4}}\right)r^{2}+\dfrac{2Q^{2}}{5z^{6}}r^{4}+\mathcal{O}\left[r^{5}\right],$
(20)

which is related with the finiteness of the self-energy density of the
electric NED source, with $\Lambda=Q^{2}/z^{4}$ being an effective
cosmological constant. In addition, the ID solution reduces to the
Schwarzschild solution in the chargeless limit ($Q\rightarrow 0$).

The event horizon of the BH solution (9) can be determined by $f(r)=0$. For
the ID solution, the equation $f(r)=0$ leads to a transcendental equation.
Hence we cannot obtain a closed expression for the event horizon radius as
functions of $M$ and $Q$, although we can obtain it numerically. The extreme
charge value, $Q_{\rm{ext}}$, can be obtained by solving $f(r)=0$ and
$f^{\prime}(r)=0$, simultaneously, where the prime denotes differentiation
with respect to the coordinate $r$. Therefore we can show that $Q_{\rm{ext}}$
for the ID solution is given by $Q_{\rm{ext}}^{\rm{ID}}\cong 1.07305M$ and the
corresponding extreme event horizon location by $r_{\rm{ext}}^{\rm{ID}}\cong
0.82532M$ (recall that for the RN BH solution,
$Q_{\rm{ext}}^{\rm{RN}}=r_{\rm{ext}}^{\rm{RN}}=M$).

In Fig. 1, we compare the metric functions of ID and RN BHs solutions, for a
given value of the normalized electric charge, defined as $q=Q/Q_{\rm{ext}}$,
which satisfies $0\leq q\leq 1$. We note that these solutions have a similar
causal structure. For $0<q<1$ we have a Cauchy horizon, $r_{-}$, and an event
horizon, $r_{+}$, while for $q=1$ the two horizons degenerate into a single
null hypersurface ($r_{\rm{ext}}$). The $q>1$ case is associated to
horizonless solutions. Here we will consider only BH solutions, which occur
when $0\leq q\leq 1$, with the case $q=0$ corresponding to the Schwarzschild
solution.

Figure 1: Comparison between the metric functions of ID
$\left[f^{\rm{ID}}(r)\right]$ and RN $\left[f^{\rm{RN}}(r)\right]$ BH
solutions, considering three distinct values of $q$ = $0.4$, $0.7$, and $1$,
as functions of $r/M$.

As a means to verify the regularity of the ID solution, we compute the
Kretschmann scalar, defined as $K\equiv
R_{\mu\nu\sigma\rho}R^{\mu\nu\sigma\rho}$. For the ID solution, $K$ is given
by

$\displaystyle K(r)=\ $
$\displaystyle\dfrac{64M^{2}}{\pi^{2}y^{6}A^{6}z^{6}}\big{(}3A^{6}B^{2}-2A^{3}By\left(9y^{4}+8y^{2}+3\right)$
$\displaystyle+y^{2}\left(7y^{4}+4y^{2}+1\right)\left(5y^{4}+4y^{2}+3\right)\big{)},$
(21)

where we defined the auxiliary functions:

$\displaystyle A(r)=\ $ $\displaystyle 1+y^{2},$ (22a) $\displaystyle B(r)=\ $
$\displaystyle\arctan y,$ (22b)

and $y\equiv r/z$. In Fig. 2, we display the behavior of the Kretschmann
scalar of the ID solution. We see that this scalar is finite for $r\geq 0$, as
long as $Q\neq 0$, which is enough to avoid the existence of curvature
singularities Bronnikov and Rubin (2021).

Figure 2: Kretschmann scalar of the ID RBH solution, normalized by
$\pi^{2}z^{6}$, as a function of $r/z$. At the core of the ID solution, the
Kretschmann scalar is finite and given by
$\left.K(r)\right|_{0}=512M^{2}/3\pi^{2}z^{6}$.

##  III Null geodesics

In this section we present the equations of motion for null geodesic in the
standard geometry (SG) [see Eq. (9)], as well as in the effective geometry
(EG) [see Eq. (33)], where photons in NED theory propagate along null
geodesics. Due to the spherical symmetry, we consider the motion in the
equatorial plane, i.e., $\theta=\pi/2$, without loss of generality.

###  III.1 Null rings in the standard geometry

The classical Hamiltonian $\mathrm{H}_{\rm{geo}}$ that provides the equations
of motion for massless particles is given by222In the remainder of this paper,
we use the term “massless particles” to refer to any particle that follow null
geodesics in the SG (9).

$\mathrm{H}_{\rm{geo}}=\dfrac{1}{2}g^{\mu\nu}p_{\mu}p_{\nu}=\dfrac{1}{2}\left(\dfrac{p_{t}^{2}}{f(r)}-f(r)p_{r}^{2}-\dfrac{p_{\varphi}^{2}}{r^{2}}\right),$
(23)

where $p_{\mu}$ are the components of the 4-momentum of massless particles. By
using the Hamilton’s equations, we obtain

$\displaystyle\dot{t}$ $\displaystyle=\dfrac{p_{t}}{f(r)},$ (24)
$\displaystyle\dot{r}$ $\displaystyle=-f(r)p_{r},$ (25)
$\displaystyle\dot{\varphi}$ $\displaystyle=-\dfrac{p_{\varphi}}{r^{2}}.$ (26)

Since the Hamiltonian (23) does not depend explicitly on the coordinates $t$
and $\varphi$, $p_{t}\equiv E$ and $p_{\varphi}\equiv-L$ are constants of
motion, where $E$ and $L$ are the energy and angular momentum of the massless
particles, respectively. Recall also that for null geodesics in the SG
$\mathrm{H}_{\rm{geo}}=0$.

Using Eqs. (24)-(26), and $\mathrm{H}_{\rm{geo}}=0$, we may obtain a radial
equation for massless particles, given by

$\dot{r}^{2}+\mathrm{V}(r)=E^{2},$ (27)

where $\mathrm{V}(r)$ is the effective potential for the radial motion of
particles following null geodesics, defined as

$\mathrm{V}(r)\equiv L^{2}\dfrac{f(r)}{r^{2}}.$ (28)

In Fig. 3, we display the effective potential for massless particles on the ID
RBH background. Notice that the local maximum of the effective potential
increases as we consider higher values of the normalized electric charge. At
the local maximum of the effective potential we have unstable circular orbits
for massless particles.

Figure 3: Effective potential for null geodesics on the background of the ID
RBH solution, normalized by the angular momentum, as function of $r/r_{+}$,
and for distinct values of $q$.

Closed circular null orbits are described by $\dot{r}=0$ and $\ddot{r}=0$,
which implies that

$\mathrm{V}=E^{2}\ \ \ \text{and}\ \ \ \mathrm{V}^{\prime}=0,$ (29)

respectively. Moreover if $d^{2}\mathrm{V}/dr^{2}<0$, the closed circular
orbit is unstable. From Eqs. (29), we may find the radius $r_{c}$ of the
unstable circular orbit and the corresponding critical impact parameter
$b_{c}$ at this orbit, namely

$\displaystyle 2f_{c}-$ $\displaystyle r_{c}f^{\prime}_{c}=0,$ (30)
$\displaystyle b_{c}=\dfrac{L_{c}}{E_{c}}$
$\displaystyle=\dfrac{r_{c}}{\sqrt{f_{c}}}.$ (31)

The impact parameter is defined as $b\equiv L/E$ and the subscript “c” denotes
that the quantity under consideration is computed at the critical radius
$r_{c}$. In Fig. 4, we display the motion of massless particles on an ID RBH
spacetime with $q=0.8$, obtained by solving Eqs. (26) and (27) numerically.
For $b<b_{c}$, massless particles are absorbed by the central object, while
for $b>b_{c}$ they are scattered. At the threshold, when $b=b_{c}$, the
trajectories describe a circular orbit around the BH at $r=r_{c}$. Therefore,
we can interpret $b_{c}$ as the threshold between absorbed and scattered null
geodesics.

Figure 4: Null trajectories on the ID RBH with $q=0.8$, for distinct choices
of $b$. In this figure, $r_{c}=2.3929M$ and $b_{c}=4.4309M$, with the solid
and dashed black circles denoting the orbit at $r_{c}$ and the event horizon
location, respectively. Here the numerical infinity was placed at
$r_{\infty}=100M$.

In linear electrodynamics, photons follow null geodesics of the SG.
Consequently, the equations of motion for massless particles and photons
coincide. On the other hand, in NED theory, photons are interpreted to follow
the null geodesics of an effective geometry (EG) Gutiérrez  _et al._ (1981),
which is different of the SG. Accordingly, the null geodesics analysis in the
SG (9), considering NED-based RBHs, concern only to massless particles with a
nature other than electromagnetic. In other words, the trajectories examined
above do not describe photon motion. In Sec. III.2 we analyze the appropriate
equations that govern the photon trajectories in NED-based spacetimes.

###  III.2 Light rings (LRs) in the effective geometry

In NED, electromagnetic fluctuations propagate along an effective light cone,
that in general differs from the “light-cones” defined by the standard
geometry Plebanski (1970); Boillat (1970). In fact, for a general theory of
NED, depending on the two independent four dimensional relativistic
invariants, $F$ (defined above) and $F_{\mu\nu}\star F^{\mu\nu}$, there are
(in general) two effective light cones, one for each polarization. This
encodes the phenomenon of birefringence, which substantiates a medium
interpretation for electromagnetic fluctuations propagating on a NED
background (regardless of the coupling to gravity). For the particular case of
NED models depending solely on $F$ (no dependence on $F_{\mu\nu}\star
F^{\mu\nu}$), birefringence does not occur in general 333The birefringence
phenomena can take place for NED models that depends only on $F$ in the
presence of external magnetic fields Gaete and Helayël-Neto (2017, 2021)..
Then, the single effective light cone can be made geometric by considering
that photons propagate along null geodesics of an effective metric tensor
$g^{\mu\nu}_{\rm{eff}}$, which depends on the contributions of the NED source
to the energy-momentum tensor Plebanski (1970); Boillat (1970); Gutiérrez  _et
al._ (1981); Novello  _et al._ (2000).

The effective metric tensor of an electrically charged RBHs obtained in the
$P$ framework is given by Novello  _et al._ (2000)

$g^{\mu\nu}_{\rm{eff}}=\dfrac{1}{\mathcal{H}_{P}}g^{\mu\nu}+4\dfrac{\mathcal{H}_{PP}}{F_{P}}P^{\mu}_{\
\ \sigma}P^{\sigma\nu}.$ (32)

If the NED source is characterized by a purely electric field in the $P$
framework, with the SG given by (9), then the corresponding line element of
the EG can be written as Bronnikov (2001)

$ds^{2}_{\rm{eff}}=g_{\mu\nu}^{\rm{eff}}dx^{\mu}dx^{\nu}=\dfrac{1}{\Phi}\left(f(r)dt^{2}-\dfrac{dr^{2}}{f(r)}-\Phi\mathcal{H}_{P}r^{2}d\Omega^{2}\right),$
(33)

where $\Phi\equiv\mathcal{H}_{P}/F_{P}$. From Eq. (4), we notice that the
scalar $F$ can be written as $F=(\mathcal{H}_{P})^{2}P$, and if Maxwell’s weak
field limit is satisfied (as it is for the ID solution) then
$\mathcal{H}(P)\rightarrow P$ and $\mathcal{H}_{P}\rightarrow 1$, for small
$P$. Consequently, $F_{P}\rightarrow 1$, since $\mathcal{H}_{P}\rightarrow 1$,
and we see that the EG (33) reduces to the SG (9), in the weak field limit.

From Eq. (33) one concludes that apart from the overall $1/\Phi$ factor, that
is irrelevant for null geodesics (modulo possible singularities), the only
difference with the respect to the SG is the angular coefficient. Thus radial
photon orbits coincide with the null geodesics of the SG.

The classical Hamiltonian $\mathrm{H}_{\rm{geo}}^{\rm{eff}}$ for the effective
metric tensor $g_{\mu\nu}^{\rm{eff}}$ is given by

$\mathrm{H}_{\rm{geo}}^{\rm{eff}}=\dfrac{1}{2}g^{\mu\nu}_{\rm{eff}}\bar{p}_{\mu}\bar{p}_{\nu}=\dfrac{\Phi}{2}\left(\dfrac{\bar{p}_{t}^{2}}{f(r)}-f(r)\,\bar{p}_{r}^{2}-\dfrac{\bar{p}_{\varphi}^{2}}{\Phi\mathcal{H}_{P}r^{2}}\right),$
(34)

where $\bar{p}_{\mu}$ are the components of the 4-momentum of photons.
Following the same procedure presented in Sec. III.1, the equations of motion
can be written as

$\displaystyle\dot{t}$ $\displaystyle=\dfrac{E\,\Phi}{f(r)},$ (35)
$\displaystyle\dot{r}$ $\displaystyle=-f(r)\Phi\bar{p}_{r},$ (36)
$\displaystyle\dot{\varphi}$ $\displaystyle=\dfrac{L}{\mathcal{H}_{P}r^{2}}.$
(37)

Using Eqs. (34)-(37), and $\mathrm{H}_{\rm{geo}}^{\rm{eff}}=0$, we obtain a
radial equation for photons given by

$\left(\dfrac{F_{P}^{2}}{\mathcal{H}_{P}^{2}}\right)\,\dot{r}^{2}+\mathrm{U}(r)=E^{2},$
(38)

where $\mathrm{U}(r)$ is the effective potential for the radial motion of
photons, defined as

$\mathrm{U}(r)\equiv L^{2}\dfrac{F_{P}f(r)}{\mathcal{H}_{P}^{2}r^{2}}.$ (39)

In Fig. 5, we show the effective potential (28), for some values of $q$. One
observes that when $0\leq q<1$, there are no stable circular photon orbits for
$r>r_{+}$, but when $q=1$ we have a stable photon orbit exactly on the extreme
event horizon, i.e., $r_{\rm{ext}}=0.82532M$. These results are similar to
those obtained in the RN geometry Pradhan and Majumdar (2011). Moreover the
profile of the effective potential shown in Fig. 5 is similar to the profile
of other NED-based RBH solutions (see, e.g., Ref. Habibina and Ramadhan
(2020)).

Figure 5: Effective potential for photons on the background of the ID RBH
solution, as function of $r/r_{+}$.

Let us now quantitatively analyze the circular photon orbits, also known as
LRs. From $\dot{r}=0$, which implies $\mathrm{U}=E^{2}$, we obtain the
critical impact parameter associated to the LR, namely

$b_{l}=\dfrac{r_{l}\left(\mathcal{H}_{P}\right)_{l}}{\sqrt{\left(F_{P}\right)_{l}f_{l}}},$
(40)

whereas from $\ddot{r}=0$, which implies $\mathrm{U}^{\prime}=0$, we get the
corresponding radial coordinate of the LR $r_{l}$, given by

$f_{l}\left[2-\dfrac{r_{l}\left(F_{P}^{\prime}\right)_{l}}{\left(F_{P}\right)_{l}}+\dfrac{2r_{l}\left(\mathcal{H}_{P}^{\prime}\right)_{l}}{\left(\mathcal{H}_{P}\right)_{l}}\right]-r_{l}f_{l}^{\prime}=0.$
(41)

The subscript “$l$” denotes that the quantity under consideration is computed
at the LR coordinate $r_{l}$. Figure 6 compares the LR perimetral radius 444We
notice that the concept of distance is very subtle in curved spacetimes. In
particular, the radial coordinate $r$ is not a geometrical invariant measure
of distance. A meaningful geometrical quantity to compare distance in two
different geometries is the perimetral radius, defined by
$\overline{r}\equiv\sqrt{g_{\varphi\varphi}}|_{\theta=\frac{\pi}{2}}$. For the
SG, we have $\overline{r}=r$, whereas, for the EG, we obtain
$\overline{r}=\sqrt{\mathcal{H}_{P}}r$. In the remainder of this paper, we
shall plot the perimetral radius of the LR to compare radial distance in
different geometries. and the critical impact parameter of the ID and RN BHs
solutions. Generically, we see that these quantities diminish as we increase
the charge. The LR perimetral radius of the ID RBH solution is typically
smaller than the RN one, for the same values of $q$. For its turn, the
critical impact parameter of the ID RBH solution is smaller than the RN one
only up to $q\approx 0.8659\equiv q_{\rm{cri}}$.

Figure 6: Top panel: Comparison between the LR perimetral radius, given by
$\overline{r}=\sqrt{\mathcal{H}_{P}}r$, of the ID and RN BH solutions. Bottom
panel: Comparison between the critical impact parameter of the ID and RN BH
solutions, as functions of $q$.

In terms of $z$, the functions $\mathcal{H}_{P}$ and $F_{P}$ are given by:

$\mathcal{H}_{P}(r)=\dfrac{r^{6}}{\left(r^{2}+z^{2}\right)^{3}}\ \ \
\text{and}\ \ \
F_{P}(r)=\dfrac{r^{12}\left(r^{2}-2z^{2}\right)}{\left(r^{2}+z^{2}\right)^{7}},$
(42)

respectively. To ensure that the effective geometry does not flip its
signature along photon’s geodesic, we need to require that functions
$\mathcal{H}_{P}(r)$ and $F_{P}(r)$ must be positive. The function
$\mathcal{H}_{P}(r)$ is everywhere finite and positive for $r>0$. On the other
hand, the function $F_{P}(r)$ is zero at

$r=\sqrt{2}z\equiv r_{\rm{eff}}.$ (43)

For $r<r_{\rm{eff}}$, the signature of the metric changes. This also happens
to magnetically charged NED-based RBH solutions Allahyari  _et al._ (2020). In
Fig. 7, we compare the location of the event horizons and of the effective
radius $r_{\rm{eff}}$. We see that the region where the line element of the EG
changes its signature is always inside the event horizon. Hence, the motion of
photons outside the event horizon will not be affected by the sign flip of the
coordinates $t$ and $r$, which occurs only for $r<r_{+}$.

Figure 7: Comparison between the location of the event horizon with that of
the effective radius, as functions of $q$.

##  IV Shadows and gravitational lensing

###  IV.1 Observational setup consistent with NED

In this section, we discuss the observational setup consistent with a NED
model and the applications to the shadows and gravitational lensing on the
background of the ID RBH solution, considering the EG (33). We apply backwards
ray-tracing techniques Bohn  _et al._ (2015); Cunha and Herdeiro (2018), in
order to simulate the visual appearance of the ID RBH (33). We solve
numerically the following geodesic equations in the effective metric:

$\displaystyle\dot{t}=\frac{E\Phi}{f(r)},$ (44)
$\displaystyle\dot{\varphi}=\frac{L}{\mathcal{H}_{P}r^{2}\,\sin^{2}\theta},$
(45) $\displaystyle\ddot{r}+\bar{\Gamma}^{r}_{\
\mu\nu}\dot{x}^{\mu}\,\dot{x}^{\nu}=0,$ (46)
$\displaystyle\ddot{\theta}+\bar{\Gamma}^{\theta}_{\
\mu\nu}\dot{x}^{\mu}\,\dot{x}^{\nu}=0,$ (47)

where $\bar{\Gamma}^{\alpha}_{\ \mu\nu}$ are the components of the Christoffel
symbol computed with the EG (33). The initial conditions for Eqs. (44)-(47)
are obtained by projecting the 4-momentum of the photon into the vierbein of a
given observer. We assume that the observer follows a timelike world-line (of
the SG) and has no net charge. Hence the vierbein attached to the observer is
dictated by the SG (9). We consider a static observer in the ID geometry,
which is described by the following vierbein:

$\displaystyle\hat{\lambda}^{\hat{0}}_{\ \mu}=\left(\sqrt{f^{\rm{ID}}(r)},\
0,\ 0,\ 0\right),$ (48) $\displaystyle\hat{\lambda}^{\hat{1}}_{\
\mu}=\left(0,\ \frac{1}{\sqrt{f^{\rm{ID}}(r)}},\ 0,\ 0\right),$ (49)
$\displaystyle\hat{\lambda}^{\hat{2}}_{\ \mu}=\left(0,\ 0,\ r,\ 0\right),$
(50) $\displaystyle\hat{\lambda}^{\hat{3}}_{\ \mu}=\left(0,\ 0,\ 0,\
r\sin\theta\right),$ (51)

which is obtained by adopting $\hat{\lambda}^{\hat{0}}_{\ \mu}$ as the
4-velocity of the observer, and imposing orthonormality condition555The
orthonormality condition for the vierbein implies that
$g_{\mu\nu}\hat{\lambda}_{\hat{a}}^{\ \mu}\hat{\lambda}_{\hat{b}}^{\
\nu}=\eta_{\hat{a}\hat{b}}$, where $\eta_{\hat{a}\hat{b}}$ is the Minkowski
metric.  with $\hat{\lambda}^{\hat{1}}_{\ \mu}$, $\hat{\lambda}^{\hat{2}}_{\
\mu}$, $\hat{\lambda}^{\hat{3}}_{\ \mu}$. The components of the 4-momentum of
the photon projected into the vierbein are

$\displaystyle\bar{p}_{\hat{a}}=\hat{\lambda}^{\
\mu}_{\hat{a}}\,\bar{p}_{\mu}.$ (52)

We note that the components of $\hat{\lambda}_{\hat{a}}^{\ \mu}$ are computed
using the SG, since it is related to an observer following a timelike curve,
while the components of $\bar{p}_{\mu}$ are computed using the EG, since it is
related to the motion of photons. $\bar{p}_{\hat{a}}$ are the components of
the 4-momentum of the photon as measured by a static observer in the ID
spacetime. In particular, $\bar{p}^{\hat{t}}$ is the photon frequency and
$\bar{p}^{\hat{r}},\bar{p}^{\hat{\phi}},\bar{p}^{\hat{\theta}}$ are the
components of the spatial momentum measured by the static observer.

The 4-momentum of the photon $\bar{p}^{\mu}$ is null with respect to the
effective metric tensor $g^{\mu\nu}_{\rm{eff}}$. However it is, in general, a
non-null vector with respect to the standard metric tensor $g^{\mu\nu}$. In
particular, for a local static observer, the norm of the 4-momentum is given
by

$\displaystyle\bar{p}_{\hat{a}}\bar{p}^{\hat{a}}$
$\displaystyle=\eta^{\hat{a}\hat{b}}\hat{\lambda}_{\hat{a}}^{\
\mu}\hat{\lambda}_{\hat{b}}^{\
\nu}\bar{p}_{\mu}\bar{p}_{\nu}=g^{\mu\nu}\bar{p}_{\mu}\bar{p}_{\nu}\therefore$
$\displaystyle\bar{p}_{\hat{a}}\bar{p}^{\hat{a}}$
$\displaystyle=-4\frac{\mathcal{H}_{P}\,\mathcal{H}_{PP}}{F_{P}}P^{\mu}_{\
\sigma}\,P^{\sigma\nu}\,\bar{p}_{\mu}\,\bar{p}_{\nu},$ (53)

where we used Eq. (32) in the last equality and the fact that
$g^{\mu\nu}_{\rm{eff}}\bar{p}_{\mu}\bar{p}_{\nu}=0$. Moreover, using Eqs.
(11)-(12) and the geodesic equations (35)-(38), we obtain that

$\displaystyle\bar{p}_{\hat{a}}\bar{p}^{\hat{a}}=-4\frac{\mathcal{H}_{PP}L^{2}}{\mathcal{H}_{P}\,r^{2}}\left(P_{tr}\right)^{2}\leq
0,$ (54)

which is negative since $\mathcal{H}_{P}$ and $\mathcal{H}_{PP}$ are positive
outside the event horizon. Therefore, the 4-momentum of the photon is a space-
like or null-like vector with respect to the metric tensor $g^{\mu\nu}$,
namely

$\displaystyle g^{\mu\nu}\bar{p}_{\mu}\,\bar{p}_{\nu}\leq 0.$ (55)

We note that, outside the event horizon, $\bar{p}_{\hat{a}}$ is a null vector
only for radially moving photons, since $L=0$. For non-radial geodesics,
$\bar{p}_{\hat{a}}$ is a space-like vector. Thus, from the viewpoint of the
SG, a local static observer measures photons that travel with a speed greater
than that of massless particles (cf. Sec. III.1), except for radially moving
photons that are null also from the SG perspective.

We can parameterize the spatial components of the 4-momentum in terms of the
celestial coordinates $(\alpha,\beta)$:

$\displaystyle\bar{p}^{\hat{r}}=\textbf{p}\cos\alpha\cos\beta,$ (56)
$\displaystyle\bar{p}^{\hat{\theta}}=\textbf{p}\sin\alpha,$ (57)
$\displaystyle\bar{p}^{\hat{\phi}}=\textbf{p}\cos\alpha\sin\beta,$ (58)

where p is the norm of the photon’s spatial 3-momentum. Using Eqs. (52),
(56)-(58) and (35)-(37), we obtain that

$\displaystyle E=\bar{p}^{\hat{t}}\sqrt{f^{\rm{ID}}_{0}},$ (59)
$\displaystyle\dot{r}=\textbf{p}\sqrt{f^{\rm{ID}}_{0}}\left(\Phi\right)_{0}\cos\alpha\cos\beta,$
(60)
$\displaystyle\dot{\theta}=\frac{\textbf{p}\sin\alpha}{\left(\mathcal{H}_{P}\right)_{0}\,r_{0}},$
(61) $\displaystyle L=\textbf{p}\,r_{0}\,\sin\theta_{0}\cos\alpha\sin\beta,$
(62)

where $(r_{0},\theta_{0})$ is the location of the observer. The subscript “0”
denotes that the quantity under consideration is computed at the observer’s
radial coordinate $r_{0}$. We can explicitly compute the norm of the
4-momentum $\bar{p}_{\hat{a}}$. Using Eqs. (54) and (62) we find that

$\displaystyle\bar{p}_{\hat{a}}\,\bar{p}^{\hat{a}}=(\bar{p}^{\hat{t}})^{2}-\textbf{p}^{2}=\frac{-3\,z^{2}\,\sin^{2}\theta_{0}\cos^{2}\alpha\sin^{2}\beta\,\textbf{p}^{2}}{(r_{0}^{2}+z^{2})}.$
(63)

Thus, the relation between the norm of the spatial 3-momentum and the photon’s
frequency measured by the local observer $\bar{p}^{\hat{t}}$ is

$\displaystyle(\bar{p}^{\hat{t}})^{2}=\left(1-\frac{3\,z^{2}\,\sin^{2}\theta_{0}\cos^{2}\alpha\sin^{2}\beta}{(r_{0}^{2}+z^{2})}\right)\,\textbf{p}^{2}\,.$
(64)

We note that, for an observer located at the equatorial plane
$\theta_{0}=\pi/2$, the right side of Eq. (64) is positive for any direction
($\beta,\alpha$) if the observer is located at the region

$r_{0}>r_{\rm{eff}}.$ (65)

The trajectory of photons is independent of the specific value of the local
frequency $\bar{p}^{\hat{t}}$. A change in the local frequency simply implies
in a rescaling of the affine parameter along the geodesic. Thus, we can always
choose $\bar{p}^{\hat{t}}$, such that $\textbf{p}=1$, what simplifies the
initial conditions for the ray-tracing (59)-(62). This can be achieved by
choosing

$\displaystyle\bar{p}^{\hat{t}}=\left(1-\frac{3\,z^{2}\,\sin^{2}\theta_{0}\cos^{2}\alpha\sin^{2}\beta}{(r_{0}^{2}+z^{2})}\right)^{\frac{1}{2}}.$
(66)

We note that the choice of $\bar{p}^{\hat{t}}$, such that $\textbf{p}=1$,
depends on the observation angles $(\alpha,\beta)$. Therefore, the initial
conditions for the ray-tracing with normalized p are given by

$\displaystyle
E=\left(1-\frac{3\,z^{2}\,\sin^{2}\theta_{0}\cos^{2}\alpha\sin^{2}\beta}{(r_{0}^{2}+z^{2})}\right)^{\frac{1}{2}}\sqrt{f^{\rm{ID}}_{0}},$
(67)
$\displaystyle\dot{r}=\sqrt{f^{\rm{ID}}_{0}}\left(\Phi\right)_{0}\cos\alpha\cos\beta,$
(68)
$\displaystyle\dot{\theta}=\frac{\sin\alpha}{\left(\mathcal{H}_{P}\right)_{0}\,r_{0}},$
(69) $\displaystyle L=\,r_{0}\,\sin\theta_{0}\cos\alpha\sin\beta.$ (70)

Dividing Eq. (70) by (67), we may obtain the relation between the critical
impact parameter $b_{l}$ and the observation angle of the shadow edge
$\beta_{l}$, measured in the observer frame666As far as we are aware, the
previous works about shadows in NED place the observer at spatial infinity
where the relation between the observation angle and the critical impact
parameter is rather simple, given by Eq. (72). In this work, we note that when
the observer is placed at a finite radial coordinate, the relation between the
observational angle and the critical impact parameter is non-trivial. The non-
triviality arises due to the fact that a local observer perceives the photon
as a spacelike particle [see Eq. (54)]. (considering $\alpha=0$ and
$\theta_{0}=\pi/2$)

$\displaystyle\sin\beta_{l}=\frac{b_{l}\,\sqrt{r_{0}^{2}+z^{2}}\,\sqrt{f_{0}^{\rm{ID}}}}{\sqrt{r_{0}^{2}\left(r_{0}^{2}+z^{2}\right)+3\,b_{l}^{2}\,z^{2}\,f_{0}^{\rm{ID}}}},$
(71)

while the shadow radius of the RBH in the observer’s screen is given by

$\displaystyle r_{s}=r_{0}\,\sin\beta_{l}.$ (72)

Notice that if we place the observer very far away from the RBH, i.e., for
$r_{0}\rightarrow\infty$, we have

$r_{s}=b_{l},$ (73)

with $b_{l}$ given by Eq. (40). Hence, as seen by a distant observer, the
impact parameter is the radius of the shadow. These features are in agreement
with linear electrodynamics.

The shadow boundary curve for a distant observer can be expressed in terms of
the so-called celestial coordinates $(x,y)$ as Bardeen (1973)

$r_{s}=\sqrt{x^{2}+y^{2}},$ (74)

where

$\displaystyle x$
$\displaystyle=\lim_{r_{0}\rightarrow\infty}\left(-r_{0}\,\frac{\bar{p}^{\hat{\phi}}}{\bar{p}^{\hat{t}}}\right),$
(75) $\displaystyle y$
$\displaystyle=\lim_{r_{0}\rightarrow\infty}\left(r_{0}\,\frac{\bar{p}^{\hat{\theta}}}{\bar{p}^{\hat{t}}}\right).$
(76)

The shape of the shadow can be obtained from a parametric plot of the circle
equation (74).

###  IV.2 Main results

In Fig. 8, we present some examples of shadows for ID RBH solutions, as seen
by an observer at spatial infinity. We note that the size of the shadows
decreases with the increase of $q$, as expected since the critical impact
parameter, which corresponds to the shadow radius [cf. Eq. (73)], diminishes
as we increase the charge [cf. Fig. 6].

Figure 8: Shadows of the ID RBH solution, as seen by an observer at spatial
infinity, for distinct values of $q$. We also consider the Schwarzschild case
$q=0$, for comparison.

It is possible to quantify the influence of the EG in the shadows size when
compared to the SG. To do this, we analyze the ratio between the shadow
radius, seen by an observer at infinity, obtained from Eq. (73)
($r_{\rm{s}}^{\rm{EG}}$) with the corresponding one obtained from Eq. (40)
($r_{\rm{s}}^{\rm{SG}}$), as showed in Fig. 9. We also consider the ratio
between $r_{\rm{s}}^{\rm{EG}}$ and $r_{\rm{s}}^{\rm{RN}}$. We see that
$r_{\rm{s}}^{\rm{EG}}$ is typically bigger than $r_{\rm{s}}^{\rm{SG}}$. In
particular, the highest difference between them, with
$r_{\rm{s}}^{\rm{EG}}>r_{\rm{s}}^{\rm{SG}}$, occurs for the extreme charge
case, for which $r_{\rm{s}}^{\rm{EG}}$ is $\approx 9.29\%$ bigger than
$r_{\rm{s}}^{\rm{SG}}$. Remarkably, for some $q=q_{\rm{crit}}$ (see Sec.
III.2, in particular, the bottom panel of Fig. 6) the shadow radius of ID and
RN BHs solutions coincide, with $r_{\rm{s}}^{\rm{EG}}<r_{\rm{s}}^{\rm{RN}}$
for $q<q_{\rm{crit}}$, while for $q>q_{\rm{crit}}$ one has
$r_{\rm{s}}^{\rm{EG}}>r_{\rm{s}}^{\rm{RN}}$.

Figure 9: Ratio between the shadows radius of the ID RBH solution, for two
different scenarios: (i) the effective and standard geometries (red curve);
and (ii) the effective geometry and the shadow radius of the RN BH solution
(blue-dashed curve).

At first sight, the result $r_{\rm{s}}^{\rm{EG}}>r_{\rm{s}}^{\rm{SG}}$ might
seem counter-intuitive, since null geodesics of the EG are in general
spacelike curves from the perspective of the SG. However, we notice that the
null geodesics of the EG can be interpreted as non-geodesic curves, from the
perspective of the SG, described by

$\ddot{x}^{\mu}+\,\Gamma^{\mu}_{\
\nu\beta}\,\dot{x}^{\nu}\dot{x}^{\beta}=\mathcal{F}^{\mu}_{\
\nu\beta}\,\dot{x}^{\nu}\,\dot{x}^{\beta},$ (77)

where $\Gamma^{\mu}_{\ \nu\beta}$ are the components of the Christoffel symbol
computed with the SG (33) and $\mathcal{F}^{\mu}_{\ \nu\beta}$ is a 4-force
term, whose analytical expression is given in Appendix B. In order to
substantiate the interpretation of photons following a non-geodesic curve,
submitted to a 4-force term, we show in the top panels of Fig. 10 the
trajectories of photons (continuous lines), compared to the trajectories of
massless particles in the SG (blue dashed lines) for the same observational
angle $\beta$. Along the photon’s trajectories we show, as a color map plot,
the absolute value of the 4-force along the radial direction. The regions in
red have a larger absolute value, while the regions in blue have a smaller
absolute value. We also show in the bottom panels of Fig. 10, the 4-force term
$\left(-\mathcal{F}^{r}_{\ \nu\beta}\dot{x}^{\nu}\dot{x}^{\beta}\right)$ as a
function of the radial coordinate. We notice that the 4-force term along the
radial direction is negative, meaning that photons experience an additional
attractive force, when compared to massless particles moving in the SG. Due to
this additional force, the photon is captured by the BH while the massless
particle (with the same observational angle $\beta$) is scattered to infinity.
Hence, such 4-force term in Eq. (77) explains why the shadows computed with
the EG are always larger than the shadows of massless particles in the SG.

Figure 10: Top panel: comparison between the trajectories of photons
(continuous lines) and the trajectories of massless particles in the SG (blue
dashed lines) for the same observational angle $\beta$. The color map plots
represent the contribution ($-\mathcal{F}^{r}_{\
\nu\beta}\dot{x}^{\nu}\dot{x}^{\beta}$) along the photon’s trajectory. The
black dashed circles represent the event horizon. Bottom panel: the
contribution ($-\mathcal{F}^{r}_{\ \nu\beta}\dot{x}^{\nu}\dot{x}^{\beta}$)
computed along the photon’s trajectory as a function of the radial coordinate.
We notice that the contribution is always negative for the ID spacetime,
representing an attractive force to the BH center. In the left panels, we have
chosen ($q=0.8$, $\beta=0.41$, $r_{0}=10M$), and O is the observer’s position;
while in the right panels, we have chosen ($q=1$, $\beta=0.04$, $r_{0}=100M$).

In Fig. 11, we show the shadows and gravitational lensing for the ID RBH
solution with different values of $q$. We have chosen the observer to be
located at $r_{0}=15M$ and $\theta=\pi/2$. This figure was obtained using
backwards ray-tracing techniques, which consists in evolving the light rays
from the observer position, and backwards in time, until it reaches a colored
celestial sphere with radius $r_{cs}=30M$ or falling to the event horizon. The
numerical code was written in C++ and it is a slightly modified version of the
code used in Refs. Junior  _et al._ (2021a, 2022). From Fig. 11 we notice that
the shadow decreases as we increase $q$, in agreement with the analytical
results presented in Fig. 8. We also notice that the gravitational lensing
varies with $q$. The major difference in the gravitational lensing arises
close to the shadow edge. Far from the shadow edge, the gravitational lensing
is essentially the same. In Appendix A, we derive an analytical approximation
for the scattering angle in the weak field limit, and we notice that the lower
contribution of the charge to the scattering is quadratic.

(a) Schwarzschild

(b) q = 0.3

(c) q = 0.6

(d) q = 1

Figure 11: Shadow and gravitational lensing of the ID RBH solution for
distinct values of $q$, considering the EG (33). In this figure, we have
chosen an observer located at $r_{0}=15M$ and $\theta=\pi/2$. We also
considered the Schwarzschild case ($q=0$) for comparison.

###  IV.3 Fine tuned degenerated shadows for asymptotic observers

The situation where $r_{s}^{\rm{EG}}=r_{s}^{\rm{RN}}$ suggests that the EG
(32) may mimic the shadow properties of singular BHs, such as the RN BH
solution, as seen by an observer at spatial infinity. This property, named as
shadow degeneracy, was investigated for static, as well as stationary BHs,
that are degenerated with respect to the Schwarzschild/Kerr BHs in Ref. Junior
_et al._ (2021b). To address the possibility of the ID RBH to be shadow
degenerated with respect to the RN BH, as seen by an observer at spatial
infinity, we may begin by searching for the values of the pairs
$(q^{\rm{ID}},q^{\rm{RN}})$, for which their corresponding $b_{l}$ coincide.
We name this property as fine tuned degenerated shadows, since we need to fine
tune the charges for the shadow to be degenerate. The fine tuned charge pairs
are shown in Fig. 12. We notice that it is possible to find fine tuned shadow
degenerated solutions for $q^{\rm{RN}}\lesssim 0.9728$.

Figure 12: Values of the normalized electric charges for which the ID and RN
BHs solutions are fine tuned shadow degenerated.

In Fig. 13, we compare the shadows of the RN BH with that of the ID RBHs for
some pairs $(q^{\rm{ID}},q^{\rm{RN}})$, for which the shadows are degenerated.
The ID RBH and the RN BH shadows can not be distinguished, as seen by a
distant observer, for low-to-extreme values of the normalized electric charge.

Figure 13: Comparison of the shadows of the RN BH with the shadows of the ID
RBH, as seen by an observer at infinity, considering the pairs
$(q^{\rm{ID}},q^{\rm{RN}})=(0.2,0.2138)$, $(0.4,0.4231)$, $(0.6,0.6231)$,
$(0.8,0.8089)$, $(1,0.9728)$. Recall that the smaller the value of the
normalized electric charges, the bigger the shadows.

##  V Final Remarks

With the recent experimental tests of NED Aaboud  _et al._ (2017); Aad  _et
al._ (2019); Mignani  _et al._ (2017); Adam  _et al._ (2021), it is clear that
a full comprehension of the nature of electromagnetic fields requires the
consideration of nonlinear effects in the appropriated field regime. However,
the imprints of these nonlinearities in the astrophysical environment of BHs
still need to be better understood. By studying LRs, shadows, and
gravitational lensing of the ID RBH solution, and considering the effective
geometry, which describes the motion of photons, we revealed some imprints of
NED in BH physics. Our main results can be summarized as follows:

(i) We performed the shadow analysis for an observer at a finite radial
coordinate, as well as for an observer at spatial infinity. We noticed that
the observer in a finite radial coordinate perceives the photon as a spacelike
particle (one may interpret that the NED “accelerates” the photons to a
superluminal speed). The fact that photons are perceived as spacelike
particles by a local observer has implications to the shadow’s observational
angle, as shown in Eq. (71). For an observer placed at spatial infinity, we
recover results available in the literature, namely that the critical impact
parameter is the radius of the shadow.

(ii) We obtained that the shadow size decreases as we consider higher values
of the electric charge, in agreement with linear electrodynamics. In addition,
the effective geometry can increase the shadow radius in more than $9\%$, in
comparison with the standard geometry. We explained the fact that the shadow
size computed with the effective geometry is larger than the standard geometry
by writing the photon’s equation of motion as a non-geodesic curve submitted
to a 4-force term, from the standard geometry perspective. We obtained an
analytic expression for the 4-force term and showed that it acts as a radially
attractive force, thus increasing the shadow size for the effective geometry.

(iii) For $0<q<q_{\rm{cri}}$, the shadow radius of the ID solution is smaller
than the RN one, while for $q>q_{\rm{cri}}$, it is bigger, since the shadow
radius corresponds to the critical impact parameter [cf. Eq. (73)]. At the
threshold value, i.e., $q=q_{\rm{cri}}$, the shadows of ID and RN BHs solution
are the same. Noticeably, it is possible to find other configurations for
which the shadows of these BHs coincide. We named these configurations as fine
tuned degenerated shadows, since it is necessary to fine tune the electric
charges in order to obtain two geometries with degenerated shadows.

(iv) We also observed that the main difference in the gravitational lensing
appears close to the shadow edge. In the weak field limit, the gravitational
lensing is essentially the same, since the contributions of the charge are
very small.

As an extension of this work, the study of other optical phenomena, such as
the birefringence, can be performed in future work, aiming to reveal more
signatures of NED within BH physics. We also plan to consider rotating NED-
based RBH spacetimes due to their relevance in astrophysical scenarios.

###### Acknowledgements.

The authors are grateful to Carlos Herdeiro for his important contributions to
this work. We acknowledge Fundação Amazônia de Amparo a Estudos e Pesquisas
(FAPESPA), Conselho Nacional de Desenvolvimento Científico e Tecnológico
(CNPq) and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES)
– Finance Code 001, from Brazil, for partial financial support. MP thanks the
University of Sheffield, in England, while LC and HLJ thank University of
Aveiro, in Portugal, for the kind hospitality during the completion of this
work. This work is supported by the Center for Research and Development in
Mathematics and Applications (CIDMA) through the Portuguese Foundation for
Science and Technology (FCT - Fundacão para a Ciência e a Tecnologia),
references UIDB/04106/2020 and UIDP/04106/2020. PC is supported by the
Individual CEEC program 2020 funded by the FCT. There has been further support
from the projects CERN/FISPAR/0027/2019, PTDC/FIS-AST/3041/2020,
CERN/FISPAR/0024/2021 and 2022.04560.PTDC. This work has further been
supported by the European Union’s Horizon 2020 research and innovation (RISE)
programme H2020-MSCARISE-2017 Grant No. FunFiCO-777740 and by the European
Horizon Europe staff exchange (SE) programme HORIZONMSCA-2021-SE-01 Grant No.
NewFunFiCO-101086251.

##  Appendix A Weak deflection angle in the ID metric by using the geodesic
method

In Sec. IV.2, we numerically computed the gravitational lensing of the ID RBH
solution. In this appendix, we derive an expression for the deflection angle
of the ID metric in the weak field limit by using the geodesic method and
considering the standard and effective geometries.

Considering the SG, we can rewrite Eq. (27) as

$\left(\dfrac{dr}{d\varphi}\right)^{2}=\dfrac{r^{4}}{b^{2}}-r^{2}f(r),$ (78)

and, similarly for the EG, Eq. (38) as

$\left(\dfrac{dr}{d\varphi}\right)^{2}=\dfrac{r^{4}\mathcal{H}_{P}^{4}}{F_{P}^{2}b^{2}}-\dfrac{r^{2}\mathcal{H}_{P}^{2}f(r)}{F_{P}}.$
(79)

The impact parameter associated with the radius of maximum approximation of
the particle $r_{m}$ is obtained by solving
$\left(dr/d\varphi\right)|_{r=r_{m}}=0$ for $b$. For the SG, we get

$b=\dfrac{r_{m}}{\sqrt{f_{m}}},$ (80)

while for the EG we have

$b=\dfrac{r_{m}\left(\mathcal{H}_{P}\right)_{r_{m}}}{\sqrt{\left(F_{P}\right)_{r_{m}}f_{m}}}.$
(81)

The deflection angle of the scattered massless particle is Newton (1982)

$\Theta(b)=\gamma(b)-\pi,$ (82)

where

$\gamma(b)=2\int_{r_{m}}^{\infty}\left(\dfrac{dr}{d\varphi}\right)^{-1}dr.$
(83)

Since we are interested in the weak field limit, we can expand the integrand
of Eq. (83) in powers of $1/r$, considering Eqs. (78) or (79), up to the
fourth order. The radius $r_{m}$ as function of $b$ is obtained by solving
Eqs. (80) or (81) and expanding the results in powers of $2M/b$ up to the
fourth order. Following these steps, we obtain that the weak deflection angle
for the SG is given by

$\displaystyle\Theta(b)=\ $
$\displaystyle\dfrac{4M}{b}+\dfrac{3\pi\left(5M^{2}-Q^{2}\right)}{4b^{2}}+\dfrac{16M\left(8M^{2}-3Q^{2}\right)}{3b^{3}}+$
$\displaystyle\dfrac{5Mz\left(21\left(33M^{4}-18M^{2}Q^{2}+Q^{4}\right)+8Q^{2}z^{2}\right)}{8b^{4}Q^{2}}$
$\displaystyle+\mathcal{O}\left[\dfrac{1}{b^{5}}\right].$ (84)

Note that up to the third order in $1/b$, the results for the ID RBH solution,
considering the SG, coincide with the RN result Keeton and Petters (2005);
Crispino  _et al._ (2009), with $4M/b$ being Einstein’s deflection angle Wald
(2010). However for $b^{-n}$, with $n\geq 4$, the results are different, due
to the higher order contributions of the ID metric function (18) in the far
field.

In the case of the EG we obtain

$\displaystyle\Theta(b)=\ $
$\displaystyle\dfrac{4M}{b}+\dfrac{3\pi\left(5M^{2}-Q^{2}\right)}{4b^{2}}+\dfrac{6Mz^{3}}{b^{2}Q^{2}}+$
$\displaystyle\dfrac{16M\left(8M^{2}-3Q^{2}\right)}{3b^{3}}+\dfrac{16Mz^{2}}{b^{3}}+\mathcal{O}\left[\dfrac{1}{b^{4}}\right].$
(85)

We see that the weak deflection angle, computed considering the EG, reproduces
the results of the SG with corrections for $b^{-n}$, with $n\geq 2$. These
corrections can be related with the nonlinearity of the NED source, since the
EG is a direct consequence of the nonlinearities of the electromagnetic field.
Besides that, in the chargeless limit $(Q\rightarrow 0)$, we obtain the
Schwarzschild deflection angle, as expected.

##  Appendix B The description of photon’s motion from the SG perspective

In Sec. III, the photons followed null geodesics of an effective geometry,
which is different from the standard spacetime geometry. This is the standard
approach, adopted by several authors, concerning the motion of photons in NED
geometries. In this Appendix, we propose an alternative (but equivalent)
interpretation for the motion of photons in NED spacetimes. Namely, we show
that, from the perspective of the SG, the motion of photons can be interpreted
as a non-geodesic curve submitted to a 4-force term $\mathcal{F}^{\mu}_{\
\alpha\beta}$.

In order to establish this result, we rewrite Eq. (32) as

$g^{\mu\nu}_{\rm{eff}}=\dfrac{1}{\mathcal{H}_{P}}g^{\mu\nu}+h^{\mu\nu},\qquad
h^{\mu\nu}=4\dfrac{\mathcal{H}_{PP}}{F_{P}}P^{\mu}_{\ \ \sigma}P^{\sigma\nu}.$
(86)

Using $g^{\mu\nu}_{\rm{eff}}\,g_{\nu\lambda}^{\rm{eff}}=\delta^{\mu}_{\
\lambda}$, we obtain an analytical expression for the covariant components of
$g_{\nu\lambda}^{\rm{eff}}$, given by

$\displaystyle
g_{\mu\nu}^{\rm{eff}}=\mathcal{H}_{P}\,g_{\mu\nu}+\Sigma_{P}\,h_{\mu\nu},$
(87)

where

$\displaystyle\Sigma_{P}\equiv-\frac{\mathcal{H}_{P}^{2}}{1-2\,P\,\mathcal{H}_{PP}\Phi}.$
(88)

We notice that the geodesic equation for the effective geometry is written as

$\displaystyle\ddot{x}^{\mu}+\bar{\Gamma}^{\mu}_{\
\nu\alpha}\dot{x}^{\nu}\,\dot{x}^{\alpha}=0.$ (89)

Using Eqs. (86) and (87) into the geodesic equation (89), we obtain that:

$\displaystyle\ddot{x}^{\mu}+\Gamma^{\mu}_{\
\nu\alpha}\dot{x}^{\nu}\,\dot{x}^{\alpha}=\mathcal{F}^{\mu}_{\
\nu\alpha}\dot{x}^{\nu}\,\dot{x}^{\alpha},$ (90)

where $\Gamma^{\mu}_{\
\nu\alpha}\equiv\frac{1}{2}g^{\mu\beta}\left(\partial_{\nu}g_{\beta\alpha}+\partial_{\alpha}g_{\nu\beta}-\partial_{\beta}g_{\nu\alpha}\right)$
are the components of the Christoffel symbol computed with the SG (9), and

$\displaystyle\mathcal{F}^{\mu}_{\ \nu\alpha}=$
$\displaystyle\frac{1}{2}\left(\dfrac{1}{\mathcal{H}_{P}}g^{\mu\nu}+h^{\mu\nu}\right)\left[\left(\partial_{\beta}\mathcal{H}_{P}\right)\,g_{\nu\alpha}-\left(\partial_{\alpha}\mathcal{H}_{P}\right)\,g_{\nu\beta}-\left(\partial_{\nu}\mathcal{H}_{P}\right)\,g_{\alpha\beta}+\partial_{\beta}\left(\Sigma_{P}h_{\nu\alpha}\right)-\partial_{\nu}\left(\Sigma_{P}h_{\beta\alpha}\right)-\partial_{\alpha}\left(\Sigma_{P}h_{\nu\beta}\right)\right]$
$\displaystyle-\mathcal{H}_{P}\,h^{\mu\beta}\Gamma_{\beta\nu\alpha},$ (91)

is interpreted as a 4-force term that acts on photons along their world-line.
Hence we conclude that, from the SG perspective, the motion of photons are
described as a non-geodesic curves subjected to a 4-force term
$\mathcal{F}^{\mu}_{\ \alpha\nu}$. In Sec. IV, we show the contribution from
the 4-force term along the photon’s motion. We notice that the contribution is
negative along the radial direction. Therefore the photons experience an
additional inward force in the radial direction, arising due to the NED. This
explains why the shadows computed with the EG are always larger when compared
to the SG shadows.

## References

  * Will (2014) C. M. Will, “The Confrontation between General Relativity and Experiment,” Living Rev. Relativ. 17, 4 (2014), arXiv:1403.7377 [gr-qc] .
  * Abbott  _et al._ (2016) B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), “Observation of Gravitational Waves from a Binary Black Hole Merger,” Phys. Rev. Lett. 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Akiyama  _et al._ (2019) K. Akiyama _et al._ (Event Horizon Telescope Collaboration), “First M87 Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole,” ApJL 875, L1 (2019), arXiv:1906.11238 [astro-ph.GA] .
  * Hawking and Ellis (1973) S. W. Hawking and G. F. R. Ellis,  _The Large Scale Structure of Space-Time_ (Cambridge University Press, Cambridge, England, 1973).
  * Bardeen (1968) J. M. Bardeen, “Non-singular general-relativistic gravitational collapse,” in  _Proc. Int. Conf. GR5, Tbilisi_ , Vol. 174 (1968) p. 174.
  * Dymnikova (1992) I. Dymnikova, “Vacuum nonsingular black hole,” Gen. Relat. Gravit. 24, 235–242 (1992).
  * Ayón-Beato and García (1998) E. Ayón-Beato and A. García, “Regular black hole in general relativity coupled to nonlinear electrodynamics,” Phys. Rev. Lett. 80, 5056–5059 (1998), arXiv:gr-qc/9911046 [gr-qc] .
  * Ayón-Beato and García (1999a) E. Ayón-Beato and A. García, “Nonsingular charged black hole solution for nonlinear source,” Gen. Relat. Gravit. 31, 629–633 (1999a), arXiv:gr-qc/9911084 [gr-qc] .
  * Ayón-Beato and García (1999b) E. Ayón-Beato and A. García, “New regular black hole solution from nonlinear electrodynamics,” Phys. Lett. B 464, 25–29 (1999b), arXiv:hep-th/9911174 [hep-th] .
  * Ayón-Beato and García (2000) E. Ayón-Beato and A. García, “The Bardeen model as a nonlinear magnetic monopole,” Phys. Lett. B 493, 149–152 (2000), arXiv:gr-qc/0009077 [gr-qc] .
  * Rodrigues and de S. Silva (2018) M. E. Rodrigues and M. V. de S. Silva, “Bardeen Regular Black Hole With an Electric Source,” JCAP 06, 025 (2018), arXiv:1802.05095 [gr-qc] .
  * Greiner and Reinhardt (2008) W. Greiner and J. Reinhardt,  _Quantum electrodynamics_ (Springer Science & Business Media, Berlin, Germany, 2008).
  * Born (1934) M. Born, “On the quantum theory of the electromagnetic field,” Proc. Roy. Soc. Lond. A 143, 410–437 (1934).
  * Born and Infeld (1934) M. Born and L. Infeld, “Foundations of the new field theory,” Proc. Roy. Soc. Lond. A 144, 425–451 (1934).
  * Heisenberg and Euler (1936) W. Heisenberg and H. Euler, “Folgerungen aus der Diracschen Theorie des Positrons,” Z. Physik 98, 714–732 (1936).
  * Aaboud  _et al._ (2017) M. Aaboud _et al._ (ATLAS Collaboration), “Evidence for light-by-light scattering in heavy-ion collisions with the ATLAS detector at the LHC,” Nature Phys. 13, 852–858 (2017), arXiv:1702.01625 [hep-ex] .
  * Aad  _et al._ (2019) G. Aad _et al._ (ATLAS Collaboration), “Observation of light-by-light scattering in ultraperipheral Pb+Pb collisions with the ATLAS detector,” Phys. Rev. Lett. 123, 052001 (2019), arXiv:1904.03536 [hep-ex] .
  * Mignani  _et al._ (2017) R. P. Mignani, V. Testa, D. G. Caniulef, R. Taverna, R. Turolla, S. Zane, and K. Wu, “Evidence for vacuum birefringence from the first optical-polarimetry measurement of the isolated neutron star RX J1856.5$-$3754,” Mon. Not. Roy. Astron. Soc. 465, 492–500 (2017), arXiv:1610.08323 [astro-ph.HE] .
  * Adam  _et al._ (2021) J. Adam _et al._ (STAR Collaboration), “Measurement of $e^{+}e^{-}$ Momentum and Angular Distributions from Linearly Polarized Photon Collisions,” Phys. Rev. Lett. 127, 052302 (2021), arXiv:1910.12400 [nucl-ex] .
  * Bronnikov (2001) K. A. Bronnikov, “Regular magnetic black holes and monopoles from nonlinear electrodynamics,” Phys. Rev. D 63, 044005 (2001), arXiv:gr-qc/0006014 [gr-qc] .
  * Matyjasek (2004) J. Matyjasek, “Extremal limit of the regular charged black holes in nonlinear electrodynamics,” Phys. Rev. D 70, 047504 (2004), arXiv:gr-qc/0403109 [gr-qc] .
  * Dymnikova (2004) Irina Dymnikova, “Regular electrically charged vacuum structures with de Sitter centre in nonlinear electrodynamics coupled to general relativity,” Class. Quantum Grav. 21, 4417–4429 (2004), arXiv:gr-qc/0407072 [gr-qc] .
  * Balart and Vagenas (2014) L. Balart and E. C. Vagenas, “Regular black holes with a nonlinear electrodynamics source,” Phys. Rev. D 90, 124045 (2014), arXiv:1408.0306 [gr-qc] .
  * Ma (2015) M.-S. Ma, “Magnetically charged regular black hole in a model of nonlinear electrodynamics,” Annals Phys. 362, 529–537 (2015), arXiv:1509.05580 [gr-qc] .
  * Junior  _et al._ (2015) E. L. B. Junior, M. E. Rodrigues, and M. J. S. Houndjo, “Regular black holes in $f(T)$ Gravity through a nonlinear electrodynamics source,” JCAP 10, 060 (2015), arXiv:1503.07857 [gr-qc] .
  * Fan and Wang (2016) Z.-Y. Fan and X. Wang, “Construction of Regular Black Holes in General Relativity,” Phys. Rev. D 94, 124027 (2016), arXiv:1610.02636 [gr-qc] .
  * Kruglov (2017) S. I. Kruglov, “Black hole as a magnetic monopole within exponential nonlinear electrodynamics,” Annals Phys. 378, 59–70 (2017), arXiv:1703.02029 [gr-qc] .
  * Toshmatov  _et al._ (2018) B. Toshmatov, Z. Stuchlík, and B. Ahmedov, “Comment on “Construction of regular black holes in general relativity”,” Phys. Rev. D 98, 028501 (2018), arXiv:1807.09502 [gr-qc] .
  * de S. Silva and Rodrigues (2018) M. V. de S. Silva and M. E. Rodrigues, “Regular black holes in $f(G)$ gravity,” Eur. Phys. J. C 78, 638 (2018), arXiv:1808.05861 [gr-qc] .
  * Fradkin and Tseytlin (1985) E. S. Fradkin and A. A. Tseytlin, “Nonlinear Electrodynamics from Quantized Strings,” Phys. Lett. B 163, 123–130 (1985).
  * Seiberg and Witten (1999) N. Seiberg and E. Witten, “String theory and noncommutative geometry,” JHEP 09, 032 (1999), arXiv:hep-th/9908142 [hep-th] .
  * Tseytlin (2000) A. A. Tseytlin, “Born-infeld action, supersymmetry and string theory,” The many faces of the superworld 2, 417–452 (2000), arXiv:hep-th/9908105 [hep-th] .
  * Gibbons and Herdeiro (2001) G. W. Gibbons and C. A. R. Herdeiro, “Born-Infeld theory and stringy causality,” Phys. Rev. D 63, 064006 (2001), arXiv:hep-th/0008052 [hep-th] .
  * García-Salcedo and Bretón (2000) R. García-Salcedo and N. Bretón, “Born-Infeld cosmologies,” Int. J. Mod. Phys. A 15, 4341–4354 (2000), arXiv:gr-qc/0004017 [gr-qc] .
  * De Lorenci  _et al._ (2002) V. A. De Lorenci, R. Klippert, M. Novello, and J. M. Salim, “Nonlinear electrodynamics and FRW cosmology,” Phys. Rev. D 65, 063501 (2002), arXiv:arXiv:gr-qc/9806076 [gr-qc] .
  * Novello  _et al._ (2004) M. Novello, S. E. P. Bergliaffa, and J. Salim, “Non-linear electrodynamics and the acceleration of the universe,” Phys. Rev. D 69, 127301 (2004), arXiv:astro-ph/0312093 [astro-ph] .
  * Campanelli  _et al._ (2008) L. Campanelli, P. Cea, G. L. Fogli, and L. Tedesco, “Inflation-Produced Magnetic Fields in Nonlinear Electrodynamics,” Phys. Rev. D 77, 043001 (2008), arXiv:0710.2993 [astro-ph] .
  * Plebanski (1970) Jerzy Plebanski,  _Lectures on non linear electrodynamics_ (Nordita, Copenhagen, Denmark, 1970).
  * Boillat (1970) G. Boillat, “Nonlinear electrodynamics - Lagrangians and equations of motion,” J. Math. Phys. 11, 941–951 (1970).
  * Gutiérrez  _et al._ (1981) S. A. Gutiérrez, A. L. Dudley, and J. F. Plebanski, “Signals and Discontinuities in General Relativistic Nonlinear Electrodynamics,” J. Math. Phys. 22, 2835–2848 (1981).
  * Novello  _et al._ (2000) M. Novello, V. A. De Lorenci, J. M. Salim, and R. Klippert, “Geometrical aspects of light propagation in nonlinear electrodynamics,” Phys. Rev. D 61, 045001 (2000), arXiv:gr-qc/9911085 [gr-qc] .
  * Vrba  _et al._ (2019) J. Vrba, A. Abdujabbarov, A. Tursunov, B. Ahmedov, and Z. Stuchlík, “Particle motion around generic black holes coupled to non-linear electrodynamics,” Eur. Phys. J. C 79, 778 (2019), arXiv:1909.12026 [gr-qc] .
  * Habibina and Ramadhan (2020) A. S. Habibina and H. S. Ramadhan, “Geodesic of nonlinear electrodynamics and stable photon orbits,” Phys. Rev. D 101, 124036 (2020), arXiv:2007.03211 [gr-qc] .
  * Amaro and Macías (2020) D. Amaro and A. Macías, “Geodesic structure of the Euler-Heisenberg static black hole,” Phys. Rev. D 102, 104054 (2020).
  * Habibina  _et al._ (2021) A. S. Habibina, B. N. Jayawiguna, and H. S. Ramadhan, “Bound orbits around charged black holes with exponential and logarithmic electrodynamics,” Gen. Relat. Gravit. 53, 113 (2021), arXiv:2104.12071 [gr-qc] .
  * Toshmatov  _et al._ (2021) B. Toshmatov, B. Ahmedov, and D. Malafarina, “Can a light ray distinguish charge of a black hole in nonlinear electrodynamics?” Phys. Rev. D 103, 024026 (2021), arXiv:2101.05496 [gr-qc] .
  * Cunha and Herdeiro (2018) P. V. P. Cunha and C. A. R. Herdeiro, “Shadows and strong gravitational lensing: a brief review,” Gen. Relat. Gravit. 50, 42 (2018), arXiv:1801.00860 [gr-qc] .
  * Stuchlík and Schee (2019) Z. Stuchlík and J. Schee, “Shadow of the regular Bardeen black holes and comparison of the motion of photons and neutrinos,” Eur. Phys. J. C 79, 44 (2019).
  * Stuchlík  _et al._ (2019) Z. Stuchlík, J. Schee, and D. Ovchinnikov, “Generic Regular Black Holes Related to Nonlinear Electrodynamics with Maxwellian Weak-field Limit: Shadows and Images of Keplerian Disks,” ApJ 887, 145 (2019).
  * Allahyari  _et al._ (2020) A. Allahyari, M. Khodadi, S. Vagnozzi, and D. F. Mota, “Magnetically charged black holes from non-linear electrodynamics and the Event Horizon Telescope,” JCAP 02, 003 (2020), arXiv:1912.08231 [gr-qc] .
  * Kruglov (2020) S. I. Kruglov, “The shadow of M87* black hole within rational nonlinear electrodynamics,” Mod. Phys. Lett. A 35, 2050291 (2020), arXiv:2009.07657 [gr-qc] .
  * Rayimbaev  _et al._ (2020) J. Rayimbaev, M. Figueroa, Z. Stuchlík, and B. Juraev, “Test particle orbits around regular black holes in general relativity combined with nonlinear electrodynamics,” Phys. Rev. D 101, 104045 (2020).
  * Hu  _et al._ (2021) Z. Hu, Z. Zhong, P.-C. Li, M. Guo, and B. Chen, “QED effect on a black hole shadow,” Phys. Rev. D 103, 044057 (2021), arXiv:2012.07022 [gr-qc] .
  * Dyson  _et al._ (1920) F. W. Dyson, A. S. Eddington, and C. Davidson, “A Determination of the Deflection of Light by the Sun’s Gravitational Field, from Observations Made at the Total Eclipse of May 29, 1919,” Phil. Trans. Roy. Soc. Lond. A 220, 291–333 (1920).
  * Crispino and Kennefick (2019) L. C. B. Crispino and D. Kennefick, “100 years of the first experimental test of General Relativity,” Nature Phys. 15, 416 (2019), arXiv:1907.10687 [physics.hist-ph] .
  * Shapiro (1964) I. I. Shapiro, “Fourth Test of General Relativity,” Phys. Rev. Lett. 13, 789–791 (1964).
  * Shapiro  _et al._ (1968) I. I. Shapiro, G. H. Pettengill, M. E. Ash, M. L. Stone, W. B. Smith, R. P. Ingalls, and R. A. Brockelman, “Fourth Test of General Relativity: Preliminary Results,” Phys. Rev. Lett. 20, 1265–1269 (1968).
  * Virbhadra and Ellis (2000) K. S. Virbhadra and G. F. R. Ellis, “Schwarzschild black hole lensing,” Phys. Rev. D 62, 084003 (2000), arXiv:astro-ph/9904193 [astro-ph] .
  * Bozza (2002) V. Bozza, “Gravitational lensing in the strong field limit,” Phys. Rev. D 66, 103001 (2002), arXiv:gr-qc/0208075 [gr-qc] .
  * Eiroa  _et al._ (2002) E. F. Eiroa, G. E. Romero, and D. F. Torres, “Reissner-Nordstrom black hole lensing,” Phys. Rev. D 66, 024010 (2002), arXiv:gr-qc/0203049 [gr-qc] .
  * Keeton and Petters (2005) C. R. Keeton and A. O. Petters, “Formalism for testing theories of gravity using lensing by compact objects. I. Static, spherically symmetric case,” Phys. Rev. D 72, 104006 (2005), arXiv:gr-qc/0511019 [gr-qc] .
  * Gibbons and Werner (2008) G. W. Gibbons and M. C. Werner, “Applications of the Gauss-Bonnet theorem to gravitational lensing,” Class. Quantum Grav. 25, 235009 (2008), arXiv:0807.0854 [gr-qc] .
  * Bozza (2010) V. Bozza, “Gravitational Lensing by Black Holes,” Gen. Relativ. Gravit. 42, 2269–2300 (2010), arXiv:0911.2187 [gr-qc] .
  * Arakida (2018) H. Arakida, “Light deflection and Gauss–Bonnet theorem: definition of total deflection angle and its applications,” Gen. Relativ. Gravit. 50, 48 (2018), arXiv:1708.04011 [gr-qc] .
  * Li and Zhou (2020) Z. Li and T. Zhou, “Equivalence of Gibbons-Werner method to geodesics method in the study of gravitational lensing,” Phys. Rev. D 101, 044043 (2020), arXiv:1908.05592 [gr-qc] .
  * Pantig and Rodulfo (2020) R. C. Pantig and E. T. Rodulfo, “Weak deflection angle of a dirty black hole,” Chin. J. Phys. 66, 691–702 (2020), arXiv:2003.00764 [gr-qc] .
  * Eiroa (2006) E. F. Eiroa, “Gravitational lensing by Einstein-Born-Infeld black holes,” Phys. Rev. D 73, 043002 (2006), arXiv:gr-qc/0511065 [gr-qc] .
  * Liang (2017a) J. Liang, “Strong gravitational lensing by regular electrically charged black holes,” Gen. Relativ. Gravit. 49, 137 (2017a).
  * Liang (2017b) J. Liang, “Regular Magnetic Black Hole Gravitational Lensing,” Chinese Phys. Lett. 34, 050401 (2017b).
  * Ghaffarnejad  _et al._ (2018) H. Ghaffarnejad, M. Amirmojahedi, and H. Niad, “Gravitational Lensing of Charged Ayon-Beato-Garcia Black Holes and Nonlinear Effects of Maxwell Fields,” Adv. High Energy Phys. 2018, 3067272 (2018), arXiv:1601.05749 [physics.gen-ph] .
  * Schee and Stuchlík (2019) J. Schee and Z. Stuchlík, “Effective Geometry of the Bardeen Spacetimes: Gravitational Lensing and Frequency Mapping of Keplerian Disks,” ApJ 874, 12 (2019).
  * Bohn  _et al._ (2015) A. Bohn, W. Throwe, F. Hébert, K. Henriksson, D. Bunandar, M. A. Scheel, and N. W. Taylor, “What does a binary black hole merger look like?” Class. Quantum Grav. 32, 065002 (2015), arXiv:1410.7775 [gr-qc] .
  * Zajaček  _et al._ (2019) M. Zajaček, A. Tursunov, A. Eckart, S. Britzen, E. Hackmann, V. Karas, Z. Stuchlík, B. Czerny, and J. A. Zensus, “Constraining the charge of the Galactic centre black hole,” J. Phys. Conf. Ser. 1258, 012031 (2019), arXiv:1812.03574 [astro-ph.GA] .
  * Zajaček  _et al._ (2018) M. Zajaček, A. Tursunov, A. Eckart, and S. Britzen, “On the charge of the Galactic centre black hole,” Mon. Not. Roy. Astron. Soc. 480, 4408–4423 (2018), arXiv:1808.07327 [astro-ph.GA] .
  * Zajaček and Tursunov (2019) M. Zajaček and A. Tursunov, “The Electric Charge of Black Holes: Is It Really Always Negligible,” The Observatory 139, 231–236 (2019), arXiv:1904.04654 [astro-ph.GA] .
  * Salazar  _et al._ (1987) I. H. Salazar, A. Garcia, and J. Plebanski, “Duality Rotations and Type $D$ Solutions to Einstein Equations With Nonlinear Electromagnetic Sources,” J. Math. Phys. 28, 2171–2181 (1987).
  * Dymnikova  _et al._ (2015) I. Dymnikova, E. Galaktionov, and E. Tropp, “Existence of Electrically Charged Structures with Regular Center in Nonlinear Electrodynamics Minimally Coupled to Gravity,” Adv. Math. Phys. 2015, 496475 (2015).
  * Bronnikov and Rubin (2021) K. A Bronnikov and S. G. Rubin,  _Black Holes, Cosmology and Extra Dimensions_, 2nd ed. (World Scientific Publishing, Singapore, 2021) https://www.worldscientific.com/doi/pdf/10.1142/12186 .
  * Gaete and Helayël-Neto (2017) P. Gaete and J. A. Helayël-Neto, “A note on nonlinear electrodynamics,” EPL 119, 51001 (2017), arXiv:1709.03869 [physics.gen-ph] .
  * Gaete and Helayël-Neto (2021) P. Gaete and J. A. Helayël-Neto, “Remarks on inverse electrodynamics,” Eur. Phys. J. C 81, 899 (2021), arXiv:2108.07929 [hep-ph] .
  * Pradhan and Majumdar (2011) P. Pradhan and P. Majumdar, “Circular orbits in extremal Reissner-Nordstrom spacetime,” Phys. Lett. A 375, 474–479 (2011), arXiv:1001.0359 [gr-qc] .
  * Bardeen (1973) J. M. Bardeen, “Timelike and null geodesics in the Kerr metric,” Proceedings, Ecole d’Eté de Physique Théorique: Les Astres Occlus : Les Houches, France, August, 1972, 215-240 , 215–240 (1973).
  * Junior  _et al._ (2021a) H. C. D. Lima Junior, P. V. P. Cunha, C. A. R. Herdeiro, and L. C. B. Crispino, “Shadows and lensing of black holes immersed in strong magnetic fields,” Phys. Rev. D 104, 044018 (2021a), arXiv:2104.09577 [gr-qc] .
  * Junior  _et al._ (2022) H. C. D. Lima Junior, J.-Z. Yang, L. C. B. Crispino, P. V. P. Cunha, and C. A. R. Herdeiro, “Einstein-Maxwell-dilaton neutral black holes in strong magnetic fields: Topological charge, shadows, and lensing,” Phys. Rev. D 105, 064070 (2022), arXiv:2112.10802 [gr-qc] .
  * Junior  _et al._ (2021b) H. C. D. Lima Junior, L. C. B. Crispino, P. V. P. Cunha, and C. A. R. Herdeiro, “Can different black holes cast the same shadow?” Phys. Rev. D 103, 084040 (2021b), arXiv:2102.07034 [gr-qc] .
  * Newton (1982) R. G. Newton,  _Scattering theory of waves and particles_ (Dover Publications, New York, United States, 1982).
  * Crispino  _et al._ (2009) L. C. B. Crispino, S. R. Dolan, and E. S. Oliveira, “Scattering of massless scalar waves by Reissner-Nordström black holes,” Phys. Rev. D 79, 064022 (2009), arXiv:0904.0999 [gr-qc] .
  * Wald (2010) Robert M Wald,  _General relativity_ (University of Chicago Press, Chicago, United States, 2010).