# Null geodesics and shadow of 4$D$ Einstein-Gauss-Bonnet black holes
surrounded by quintessence

Malihe Heydari-Fard1 and Mohaddese Heydari-Fard2
_1 Department of Physics, The University of Qom, 3716146611, Qom, Iran_
_2 Department of Physics, Shahid Beheshti University, G. C., Evin, Tehran
19839, Iran_ Electronic address<EMAIL_ADDRESS>address:
<EMAIL_ADDRESS>
###### Abstract

We study the properties of the null geodesics of four-dimensional Einstein-
Gauss-Bonnet black holes surrounded by quintessence matter. Due to the
quintessence correction, we discuss the radial and non-radial geodesics. For
the case of non-radial geodesics we obtain the effective potential, photon
sphere and impact parameter associated to the null geodesics in details. We
analyze the different possible orbits of massless particles (i.e. photons)
such as unstable circular orbits and unbounded orbits and investigate the role
of the Gauss-Bonnet coupling $\alpha$ and the quintessence parameter $q$ on
the null geodesic trajectories. Moreover, we study the effects of the model
parameters on the shadow cast by such black holes. The results are compared to
that of Schwarzschild black hole with and without quintessence matter and
four-dimensional Einstein-Gauss-Bonnet black holes without quintessence. As a
physical application of null geodesics, the deflection angle is calculated and
the effect of the quintessence matter on it is investigated.
PACS numbers: 97.60.Lf, 04.70.Bw, 95.36.+x
Keywords: Black holes, Radial and non-radial geodesics, Quintessence.

##  1 Introduction

Black holes are the most fascinating prediction of Einstein’s general
relativity theory. The first indirect evidence for an actual black hole in
nature came in 1971, from the X-ray astronomy . This astronomical X-rays
source was discovered in 1965 [1] so called Cygnus X-1 and was the first
source extensively accepted to be a stellar-mass black hole via dynamical
observations [2]–[3]. In recent years, observational evidence strongly
indicates for the existence of black holes in two mass ranges. The LIGO
detectors in USA and Virgo detectors in Italy discovered the gravitational
waves produced by the merger of a binary stellar-mass black hole in 1.4
billion years ago [4]. Also, strong astronomical observations in recent years
show the presence of a supermassive black hole, i.e Sgr A* at the center of
Milky Way galaxy [5]. The supermassive black holes of $10^{6}-10^{9}$ solar
masses are believed to be at the center of many galaxies in the universe.
Although the black holes have been considered as possible explanation for some
observed phenomena but not detected directly. Until 10 April 2019 when Event
Horizon Telescope (EHT) through an international collaboration provided the
first direct visual evidence for the existence of black holes. EHT revealed
the shadow images of M87* [6]–[8] an elliptical galaxy whose supermassive
black hole is 6 billion times heavier than the Sun. Therefore two
astrophysical tools the study of the black hole shadow and the gravitational
lensing could be used to understand the strong field regime of gravity and to
test the modified gravity theories. The shadow of a Schwarzschild black hole
was first discussed by Synge [9] and further explored by Luminet in [10].
Bardeen then studied the shadow cast by a Kerr black hole [11]. Also, the
study of the black hole shadow in modified theories of gravity has attracted
considerable attention in recent years [12]–[26]. For a short review on shadow
of black holes see [27] and [28].

The study of null geodesics around black holes is an important subject and
helps us to understand the properties of black holes and the geometric
structure of space-time. The exact analytical solutions of the geodesics
equations help us to have the best understanding of the particles motion along
the geodesics. But this is not always possible. Therefore, we can employ the
numerical solutions and the analytical approximation schemes. The first study
of the geodesic equations and their analytical solutions for the Schwarzschild
space-time in terms of the Weierstrass elliptic functions was carried out by
Hagihara in 1931 [29]. Then, Darwin solved the geodesic equations using the
Jacobian elliptic functions [30]. Moreover, the analytical solutions of the
geodesic equations in four-dimensional Schwarzschild space-time [31],
Schwarzschild-(anti) de Sitter and higher-dimensional Schwarzschild-(anti) de
Sitter [32]–[33], the Kerr space-time [34]–[35], Kerr-(anti) de Sitter
[36]–[37], Kerr-Newman [38], Reissner-Nordstrom [39]–[41], Reissner-
Nordstrom-(anti) de Sitter [42] and Myers-Perry space-times [43]–[44] have
been investigated. The study of geodesic motion in the context of modified
theories of gravity have been also extensively studied. The geodesic motion of
test particles in $f(R)$ modified gravity have been studied in [45]. In the
context of the brane world scenario, the null geodesics have been considered
in [46]–[47]. The geodesic analysis around black holes in Horava-Lifshitz
gravity, Einstein-Maxwell-dilaton gravity and conformal weyl gravity have been
discussed in [48]–[63], respectively. The null geodesics of Born-Infeld black
holes have been considered in [64]–[67]. In the space-time of the
Schwarzschild and Kerr black holes pierced by a cosmic string, test particles
motion has been studied in [68]–[69]. In [70]–[71], the null geodesics of
noncommutative geometry have been considered. For study of time-like and null-
like geodesics around the black holes surrounded by quintessence matter and
around hairy black holes, see [72]–[76] and [77]–[78], respectively. Also, in
the background of wormhole geometries the geodesic analysis have been
considered in [79]–[81].

In the past decades, the modified theories of gravity have been taken
considerable interest to address the theoretical and observational issues that
can not be explained in the framework of Einstein’s theory of general
relativity. A well-known example in this regard is the Gauss-Bonnet (GB)
gravitational theory which appears in the low energy limit of the string
theory [82]–[83]. As we know, in four-dimensional space-time, the integral of
the GB term is a topological invariant and does not contribute to the
gravitational field equations, unless a scalar field is coupled to the GB term
by a regular coupling function. In this context, the black hole solutions,
with various coupling functions, have been well studied since 1990s [84]–[89].

However, recently Glavan and Lin proposed a novel 4-dimensional Einstein-
Gauss-Bonnet (4$D$ EGB) gravity [90] that by re-scaling the GB coupling as
$\alpha\rightarrow\frac{\alpha}{(D-4)}$ and taking the limit $D\rightarrow 4$,
circumvent the conditions of Lovelock’s theorem. In such a way, the GB term
can make a non-trivial contribution to the gravitational field equations. This
theory preserves the number of degrees of freedom and avoids the Ostrogradsky
instability. But, it is found that taking the limit $D\rightarrow 4$ may not
be consistent and the vacua of the model are not also well-defined [91]–[97].
To solve these problems, some alternative approaches such as introducing a
counter term into the action and then take $D\rightarrow 4$ limit [98]–[99],
compactification of $D$-dimensional EGB gravity [100]–[101] and breaking the
temporal diffeomorphism invariance [102] have been suggested to obtain a
consistent EGB gravity. It is worth mentioning that these consistent theories
share the black hole solutions with the original 4$D$ EGB gravity that
attracted a great deal of attention. For instance, the black hole solutions
[103]–[111], quasinormal modes and stability analysis of black holes
[112]–[118], thermodynamics and weak cosmic censorship [119]–[122], thin
accretion disks [123]–[124], test particle motion and black hole shadow
[125]–[131], gravitational lensing [132]–[135] and also some cosmological
aspects of 4$D$ EGB gravity [136]–[141] have been extensively studied. In this
paper, we consider 4$D$ EGB black holes surrounded by quintessence matter
[142] to study the null geodesics around them and investigate the effects of
the model parameters on these trajectories.

The paper is structured as follows. In section 2, we present a brief review of
4$D$ EGB black holes surrounded by quintessence matter and some of their
properties. A detailed analysis of radial and non-radial null geodesics around
these black holes is done in section 3. In section 4, the shadow radius of
quintessence 4$D$ EGB black holes is calculated. Then, as an application of
null geodesics, we study the bending of light in the space-time of such black
holes in section 5 and finally we summarize the conclusions in section 6.

##  2 4$D$ EGB black holes with quintessence matter

In this section we introduce 4$D$ EGB black hole solutions surrounded by
quintessence matter obtained in [142]. The line element of this solution has
the following form

$ds^{2}=-f(r)dt^{2}+\frac{dr^{2}}{f(r)}+r^{2}\left(d\theta^{2}+\sin^{2}\theta
d\varphi^{2}\right),$ (1)

where

$f_{\pm}(r)=1+\frac{r^{2}}{2\alpha}\left[1\pm\sqrt{1+\frac{8\alpha
M}{r^{3}}+\frac{8\alpha q}{r^{3\omega_{q}+3}}}\right],$ (2)

that $M$ is the black hole mass, $\alpha$ is the GB coupling constant and
$\omega_{q}$ and $q$ are the state parameter of the quintessence matter and
the normalization factor (quintessence parameter), respectively. It is clear
that in the case of $q=0$, the above solution reduces to the 4$D$ EGB black
holes obtained by Glavan and Lin [90]. In the limit of $\alpha\rightarrow 0$,
the negative branch corresponds to the Schwarzschild solution surrounded by
quintessence matter, that was initially obtained by Kiselev in [143]. Thus, in
this paper we consider the negative branch of the above solution. Also, the
case of $q=0$ and $\alpha\rightarrow 0$ represents the ordinary Schwarzschild
solution. Moreover, by choosing the different values of the state parameter
$\omega_{q}$, the different black hole solutions can be obtained. For
instance, in the case of no quintessence $\omega_{q}=0$ the above metric
reduces to 4$D$ EGB black holes [90], and for $\omega_{q}=-1$ and
$q=\frac{\Lambda}{6}$ it reduces to 4$D$ EGB de Sitter solution[103]. When
$\omega_{q}=\frac{1}{3}$ and $q=-\frac{Q^{2}}{2}$, the solution (2)
corresponds to 4$D$ charged EGB black holes obtained in [103].

The behaviour of the $f(r)$ function and thus the horizons of black hole,
strongly depends on the values of $\omega_{q}$, $q$ and also the GB coupling
constant $\alpha$. The state parameter $\omega_{q}$ must be in the range
$-1<w_{q}<-\frac{1}{3}$, and the energy density of the quintessence matter is
as follows [142]

$\rho_{q}=-6w_{q}q\frac{1}{r^{3(1+w_{q})}},$ (3)

as is clear, in order to have a positive energy density the $q$ parameter must
be positive. In Figure 1 we have plotted the behaviour of the metric function
$f(r)$ for different values of $w_{q}$. For the case of $\alpha>0$, it is
clear that for central values of $w_{q}$ in the interval
$w_{q}\in(-1,-\frac{1}{3})$ there are three horizons, but when approaching to
the boundary value in the ‘right’ there is two horizon and approaching to the
‘left’ there is one horizon. For further study we refer the interested reader
to [142].

In the present work, we are interest to consider $\omega_{q}=-0.35$ and thus
the function of $f(r)$ in the equation (2) takes the following form

$f(r)=1+\frac{r^{2}}{2\alpha}\left[1-\sqrt{1+\frac{8\alpha
M}{r^{3}}+\frac{8\alpha q}{r^{1.95}}}\right].$ (4)

The behaviour of the above metric function $f(r)$ for $\alpha>0$ and
$\alpha<0$ with different values of the $q$ parameter has plotted in Figure 2.
The figure shows that for positive values of $\alpha$, the function of $f(r)$
has two real roots that are identified as the cosmological quintessential
horizon $r_{-}$ and the event horizon $r_{+}$, while for negative values of
$\alpha$ there is only one horizon at $r=r_{+}$.

As we mentioned before, the case of $q=0$ corresponds to 4$D$ EGB black holes
without quintessence. In this case, the inner and outer event horizons are
given by

$r_{\pm}={M\pm\sqrt{M^{2}-\alpha}},$ (5)

where for $\alpha>0$ there are two horizons while for $\alpha<0$ there is only
one horizon. However, in the case of $f(r)$ function in equation (4), we have
not an explicit expression for the event horizon and thus we have numerically
displayed the behaviour of $r_{+}$ and $r_{-}$ as a function of the GB
coupling in Figure 3. It is easy to see that with increasing $\alpha$ the
event horizon $r_{+}$ decreases. This is due to the fact that for $\alpha>0$,
the GB term play the role of dark energy, counteract gravity and thus the
event horizon radius takes smaller values. Also, we see that by increasing
$q$, the event horizon $r_{+}$ increases while, the quintessential horizon
$r_{-}$ decreases. Therefore, the presence of quintessence matter around 4$D$
EGB black holes increases the event horizon. The numerical results have also
summarized in Table 1. For $q=0$ the results of Table 1 are the same as the
values of Table 1 in ref [126].

Figure 1: The metric coefficient $f(r)$ as a function of the radial coordinate
$r$ for different values of $w_{q}$. Left panel: for positive GB coupling
$\alpha=0.3$. Right panel: for negative GB coupling $\alpha=-0.3$ with $q=0.1$
and $M=1$.

Figure 2: The metric coefficient $f(r)$ as a function of the radial coordinate
$r$ for different values of $q$. Left panel: for positive GB coupling
$\alpha=0.3$. Right panel: for negative GB coupling $\alpha=-0.3$ and $M=1$.

Figure 3:  Left panel: the behavior of the event horizon $r_{+}$ as a function of GB coupling. On the red curve, the point of $\alpha=q=0$ corresponds to the Schwarzschild event horizon at $r_{+}=2M$. Right panel: the cosmological quintessential horizon $r_{-}$ as a function of GB coupling for $M=1$. Table 1: The quintessential cosmological horizon $r_{-}$, the event horizon $r_{+}$, the radius of the photon sphere $r_{c}$ and the impact parameter $b_{c}$ of the photon sphere for different values of $\alpha$ and $q$ with $M=1$. $q$ | $\alpha$ | $r_{-}/M$ | $r_{+}/M$ | $r_{c}/M$ | $b_{c}$
---|---|---|---|---|---
0 | 0.7 | 0.45227 | 1.54772 | 2.61959 | 4.87959
| 0.5 | 0.29289 | 1.70711 | 2.74656 | 4.98198
| 0.3 | 0.16334 | 1.83666 | 2.85656 | 5.07351
| 0.1 | 0.05132 | 1.94868 | 2.95453 | 5.15692
| 0 | 0 | 2 | 3 | 5.19615
| -1.1 | NH | 2.44914 | 3.40373 | 5.55557
| -3.3 | NH | 3.07364 | 3.95844 | 6.07084
| -5.5 | NH | 3.54951 | 4.36744 | 6.64004
| -7.7 | NH | 3.94958 | 4.70134 | 6.78150
0.1 | 0.7 | 0.42197 | 2.10407 | 3.38368 | 7.01918
| 0.5 | 0.28175 | 2.24472 | 3.51361 | 7.15518
| 0.3 | 0.16027 | 2.36768 | 3.63108 | 7.28079
| 0.1 | 0.05104 | 2.47833 | 3.73897 | 7.39808
| 0 | 0 | 2.53005 | 3.78995 | 7.45401
| -1.1 | NH | 3.00061 | 4.26285 | 7.98741
| -3.3 | NH | 3.68224 | 4.95629 | 8.80054
| -5.5 | NH | 4.21107 | 5.49176 | 9.44392
| -7.7 | NH | 4.65903 | 5.94183 | 9.99154

##  3 Null geodesics of quintessence 4$D$ EGB black holes

The trajectory of photons in quintessence background of 4$D$ EGB black holes
describe by the null geodesics. The geodesic equations for this space-time can
be obtained from the following Lagrangian

${\cal
L}=\frac{1}{2}\left(-f(r)\dot{t}^{2}+\frac{1}{f(r)}\dot{r}^{2}+r^{2}\dot{\theta}^{2}+r^{2}\sin^{2}{\theta}\dot{\varphi}^{2}\right),$
(6)

where a dot denotes differentiation with respect to the proper time $\tau$.
Using the Euler-Lagrange equations for $t$ and $\varphi$ coordinates we have

$\dot{t}=\frac{E}{f(r)},$ (7) $\dot{\varphi}=\frac{L}{r^{2}\sin^{2}{\theta}},$
(8)

where $E$ and $L$ are two constants of motion, namely the energy and angular
momentum of test particles. We consider the motion in the equatorial plane,
$\theta=\frac{\pi}{2}$, and thus we have $\dot{\theta}=\ddot{\theta}=0$. By
substituting $\dot{t}$ and $\dot{\varphi}$ in the equation (6) the lagrangian
takes the following form

$\dot{r}^{2}+f(r)\left(\frac{L^{2}}{r^{2}}+h\right)=E^{2},$ (9)

where $2{\cal L}=h$ and $h=1$ and $h=0$ correspond to the time-like and null-
like geodesics, respectively. Comparing the above equations to the relation
$\dot{r}^{2}+V_{\rm eff}=E^{2}$, one can find the effective potential as

$V_{\rm eff}=\left(\frac{L^{2}}{r^{2}}+h\right)f(r).$ (10)

Since in the present work we focus on the null geodesics with $h=0$, the
effective potential is given by

$V_{\rm eff}=L^{2}\frac{f(r)}{r^{2}}.$ (11)

To describe the photon trajectories, one can find the following equation of
motion by eliminating the $\tau$ parameter from equations (8) and (9)

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

where the impact parameter is defined as $b\equiv\frac{L}{E}$.

In the next section we analyze both the radial and non-radial null geodesics
and investigate the role of the GB coupling constant $\alpha$ and quintessence
parameter $q$ on them.

###  3.1 Radial null geodesics $(L=0)$

The radial geodesics correspond to the motion of particles with vanishing
angular momentum, $L=0$. So, the effective potential for radial null geodesics
is zero, $V_{\rm eff}=0$, and the equations of motion for $\dot{r}$ and
$\dot{t}$ are now given by

$\dot{r}=\frac{dr}{d\tau}=\pm E,$ (13)
$\dot{t}=\frac{dt}{d\tau}=\frac{E}{f(r)}.$ (14)

From the above equations we have

$\frac{dt}{dr}=\pm\frac{1}{f(r)}=\pm\frac{1}{1+\frac{r^{2}}{2\alpha}\left[1-\sqrt{1+\frac{8\alpha
M}{r^{3}}+\frac{8\alpha q}{r^{1.95}}}\right]}.$ (15)

Also, the differential equation for the proper time is given by

$\frac{d\tau}{dr}=\pm\frac{1}{E},$ (16)

which has a solution as

$\tau=\pm\frac{r}{E}+\rm const_{\pm}.$ (17)

However, the equation (15) does not have an analytical solution, and thus we
need to solve it numerically. We have plotted the behaviour of the coordinate
time $t$ and also the proper time $\tau$ in Figure 4. Note that we study the
ingoing radial geodesics and thus consider the minus sign in these equations.
As can be seen, when particles approach the event horizon, $r\rightarrow
r_{+}$, we have $t\rightarrow\infty$, while it takes a finite period of the
proper time $\tau$, which is similar to the case of the Schwarzschild black
hole [31] and the Schwarzschild black hole surrounded by quintessence [73].

Figure 4: The coordinate time $t$ and the proper time $\tau$ as a function of
the radius $r$ for quintessence 4$D$ EGB black holes with $\alpha=q=0.1$ and
$M=1$.

###  3.2 Null geodesics with angular momentum ($L\neq 0$)

Now, we are going to study null geodesics with non-zero angular momentum.

####  3.2.1 Effective potential

In the case of non-radial geodesics the effective potential is given by

$V_{\rm
eff}=\frac{L^{2}}{r^{2}}f(r)=\frac{L^{2}}{r^{2}}\left[1+\frac{r^{2}}{2\alpha}\left(1-\sqrt{1+\frac{8\alpha
M}{r^{3}}+\frac{8\alpha q}{r^{1.95}}}\right)\right].$ (18)

The effective potential of massless particles in the space-time of 4$D$ EGB
black holes surrounded by quintessence has plotted in Figure 5. In the left
panel of the figure, we have shown the effective potential for different
values of $q$ parameter with $\alpha=0.3$ and for comparison we have also
plotted it for 4$D$ EGB black holes without quintessence, $q=0$. It is easy to
see that the potential for 4$D$ EGB black holes is higher than that for
quintessence 4$D$ EGB black holes; namely the effect of the quintessence
parameter is to decreases the effective potential. The effect of the GB
coupling $\alpha$ on the effective potential has been shown in the right panel
of the figure. We see that by increasing the GB coupling the effective
potential increases, too. For different values of $L$, the effective potential
is plotted in Figure 6. As one expected, the effective potential is larger for
larger values of $L$.

As we know, the motion of particles under the effective potential $V_{\rm
eff}(r)$ depends on their energy levels. So, according to the Figure 7, by
considering the photons coming from infinity and move radially inward
direction, the different motions of particle, can be described as follows:

  1. 1.

Scattering orbits: In region I, when the photons with $b>b_{c}$ start motion
from $r>r_{0}$, approach to the turning point at $r=r_{0}$ and thus pushed
away to the large $r$, again.

  2. 2.

Falling orbits: Another possibility in region I is the photons with $b>b_{c}$
starts in $r_{+}<r<r_{1}$. In such a case, photons cross the event horizon
$r_{+}$ and fall into the black hole.

  3. 3.

Circular orbits: For $b=b_{c}$, we have $\dot{r}=0$ at $r=r_{c}$. As is clear
from the shape of the potential, if photons start at $r>r_{c}$, they will have
unstable circular motion at $r=r_{c}$. The unstable circular orbits at
$r=r_{c}$ known as the photon sphere.

  4. 4.

Falling orbits: In region III, the photon with $b<b_{c}$ that coming from
infinity will continue its motion inward and finally enter the inside of the
black hole.

Figure 5: Left panel: The effective potential as a function of $r$ for
$\alpha=0.3$ and different values of $q$. Right panel: The effective potential
as a function of $r$ for $q=0.1$ and different values of GB coupling $\alpha$,
with $L=30$ and $M=1$. The solid curves correspond to 4$D$ EGB black holes
without quintessence matter ($q=0$) and the Schwarzschild solution surrounded
by quintessence ($\alpha=0$) in the left and right panels, respectively.
Figure 6: The effective potential as a function of the radial coordinate $r$
for different values of the angular momentum $L$, with $\alpha=0.3$, $q=0.1$
and $M=1$. Figure 7: The effective potential for quintessence 4$D$ EGB black
holes for $\alpha=0.3$, $q=0.1$ and $L=30$ with $M=1$. The different kinds of
motion are described with different values of the impact parameter $b$.

####  3.2.2 Circular orbits

As we mentioned in the previous section, the circular orbits occur at
$r=r_{c}$. In this point

$V_{\rm eff}=E_{c}^{2},$ (19)

and thus $\dot{r}=0$ and the photon circle the black hole at a constant radius
$r=r_{c}$. From the above equation we have

$\frac{dV_{\rm eff}}{dr}=0,$ (20)

which using the equation (11) leads to the following equation

$r\frac{df(r)}{dr}-2f(r)=0.$ (21)

By substituting the $f(r)$ function from equation (2), one can find the
following equation for the photon radius

$6q(1+w)-2r^{3w}\left(-3M+r\sqrt{1+\frac{8\alpha M}{r^{3}}+\frac{8\alpha
q}{r^{3(1+w)}}}\right)=0.$ (22)

Since the above equation does not have an exact solution for $\omega=-0.35$,
we numerically find the roots of it. In Figure 8, we have plotted the roots of
equation (22) which are identified with $r_{1c}$, and $r_{2c}$. Also, in left
panel of Figure 9, $r_{1c}$ and $r_{2c}$ have been plotted in one graph. The
solid red curve corresponds to the 4$D$ EGB black holes in the absence of
quintessence matter, $q=0$. The behaviour of $r_{1c}$, $r_{2c}$ and $r_{+}$
has been plotted in one graph, in the right panel of the figure and the photon
radius are compared to the radius of the event horizon. According to this
figure, we see that $r_{1c}>r_{+}$, while $r_{2c}<r_{+}$. Therefore,
$r_{1c}\equiv r_{c}$ is the radius of the unstable circular orbit since it is
always larger than the event horizon $r_{+}$. As is clear from the left panel
of Figure 8, with increasing the GB coupling constant the photon radius,
$r_{c}$, decreases while for a constant value of $\alpha$ by increasing the
$q$ parameter the photon radius increases, which is in agreement with Table 1.

The numerical results for the radius of the photon sphere for 4$D$ EGB black
hole with quintessence matter are compared to the case of 4$D$ EGB black hole,
the Schwarzschild black hole and also the Schwarzschild black hole with
quintessence matter, in Table 1. It shows that the presence of the positive GB
coupling constant decreases the radius of the photon sphere of 4$D$ EGB black
hole in comparison with the Schwarzschild black hole [126]. This is due to the
fact that the positive GB coupling $\alpha>0$, has the role of the dark
energy, counteract the gravity and decreases the instability area around the
black hole and thus the photon radius takes smaller values. However, the
quintessence matter increases the radius of the photon sphere in comparison
with 4$D$ EGB black holes as well as the Schwarzschild black hole. Also, we
see that in the limit of $\alpha\rightarrow 0$ and $q=0$, $r_{c}=3M$ and
$b_{c}=\sqrt{27}M$ which are respectively the radius of the photon sphere and
impact parameter for unstable circular orbits of the Schwarzschild space-time
[31].

The impact parameter at the radius of the photon sphere is also given by

$\frac{L_{c}}{E_{c}}=\frac{r_{c}}{\sqrt{f(r_{c})}}=\frac{r_{c}}{1+\frac{r_{c}^{2}}{2\alpha}\left[1-\sqrt{1+\frac{8\alpha
M}{r_{c}^{3}}+\frac{8\alpha q}{r_{c}^{1.95}}}\right]}\equiv b_{c}.$ (23)

The values of the impact parameter for different values of $q$ and $\alpha$
are presented in Table 1. As one can see, by increasing the GB coupling
$\alpha$ the impact parameter decreases, while with increasing $q$ the values
of impact parameter increase.

Figure 8: Left panel: The photon radius $r_{1c}$ as a function of GB coupling
constant for different values of $q$. Right panel: The photon radius $r_{2c}$
as a function of GB coupling constant for different values of $q$ and $L=30$
with $M=1$.

Figure 9: Left panel: The photon radii $r_{1c}$ and $r_{2c}$ as a function of
GB coupling constant for different values of $q$. Right panel: The comparison
of the photon radii $r_{1c}$ and $r_{2c}$ with the event horizon $r_{+}$ for
$q=0.3$, $L=30$ and $M=1$.

####  3.2.3 Unbounded orbits

By substituting the function $f(r)$ from equation (4) into the equation (12),
the equation of motion reads as

$\frac{dr}{d\varphi}=\sqrt{\frac{r^{4}}{b^{2}}-r^{2}-\frac{r^{4}}{2\alpha}\left[1-\sqrt{1+\frac{8\alpha
M}{r^{3}}+\frac{8\alpha q}{r^{1.95}}}\right]}.$ (24)

Since the above equation does not have an analytical solution, we numerically
plot this equation for different values of energy, as we explained in the
previous section. To this end, we first rewrite the above equation in terms of
the new variable $u=\frac{1}{r}$ as follows

$\frac{du}{d\varphi}=\sqrt{\frac{1}{b^{2}}-u^{2}-\frac{1}{2\alpha}\left[1-\sqrt{1+{8\alpha
Mu^{3}}+{8\alpha qu^{1.95}}}\right]}\equiv F(u).$ (25)

The geometry of the geodesics depends on the roots of the equation $F(u)=0$.
To understand the behaviour of light rays in the space-time of 4$D$ EGB black
holes surrounded by quintessence matter, we assume that light rays approach
the black hole from the infinity. By considering the different values of the
impact parameter in the above equation, we can display the different
trajectories. For $b=b_{c}$ the light rays that approach to the black hole
from $r>r_{c}$, revolve around the black hole at the radius of $r=r_{c}$. The
plot of this unstable circular orbit is shown with orange curve in Figure 10
which relates to the light rays in region II in Figure 7. For $b>b_{c}$, the
light rays approaching to the black hole from infinity, will be deflected at
the turning point $u_{0}=\frac{1}{r_{0}}$ which is the root of the function
$F(u)$, see Figure 7. This case corresponds to the unbounded null geodesics
which is shown with blue curve in the Figure 11. However, in region III with
$b<b_{c}$, the light rays approaching to the black hole from infinity cross
the event horizon and enter the black hole. We have shown these falling orbits
with the green curves in Figure 11.

Figure 10: The unstable circular orbit for the light ray approaches from
$r>r_{c}$ to the 4$D$ EGB black hole surrounded by quintessence matter. Left
panel: for positive GB coupling $\alpha=0.3$, $q=0.1$. Right panel: for
negative GB coupling $\alpha=-3.3$, $q=0.1$ and $M=1$. In each panel the event
horizon is at $r_{+}=2.36768$ and $r_{+}=3.68224$, respectively.

Figure 11: The trajectory of the photons in polar coordinate $(r,\varphi)$.
Left panel: for positive GB coupling $\alpha=0.3$, $q=0.1$. Right panel: for
negative GB coupling $\alpha=-3.3$, $q=0.1$ and $M=1$. The orange dashed curve
shows the unstable circular orbit, $b=b_{c}$, the green solid curves
correspond to the null geodesics with $b<b_{c}$ and the blue curves correspond
to $b>b{c}$. For all the null geodesics, we considered the spacing between
impact parameters as $\Delta b=0.25$. In each panel the solid disk shows the
black hole with the event horizon at $r_{+}=2.36768$ and $r_{+}=3.68224$,
respectively.

##  4 Shadow of quintessence 4$D$ EGB black holes

In order to determine the shape and the size of the black hole shadow one need
to obtain two celestial coordinates $X$ and $Y$ as follows [31]

$X=\lim_{r_{0}\rightarrow\infty}\left(-r_{0}^{2}\sin\theta_{0}\frac{d\varphi}{dr}\right),$
(26)
$Y=\lim_{r_{0}\rightarrow\infty}\left(r_{0}^{2}\frac{d\theta}{dr}\right),$
(27)

where $r_{0}$ is the distance from the black hole to the observer and
$\theta_{0}$ is the inclination angle between the line of sight of the
observer and the rotation axis of the black hole. Using the geodesic equations
and substituting the expressions $\frac{d\varphi}{dr}$ and
$\frac{d\theta}{dr}$ into the above equations one can find the equations which
relate the celestial coordinates to the constants of motion. However, in the
case of our study where the observer is situated in the equatorial plane with
the inclination angle $\theta_{0}=\frac{\pi}{2}$, it is easy to see that the
radius of the shadow is equivalent to the critical impact parameter given by
the following equation

$R_{s}\equiv\sqrt{X^{2}+Y^{2}}=b_{c}.$ (28)

In Figure 12, we have displayed the boundary of the shadow of quintessence
4$D$ EGB black holes for different values of quintessence parameter $q$ and
the GB coupling $\alpha$. In the top panels of figure we have plotted the
shadow of 4$D$ EGB black holes without quintessence, $q=0$. As is clear, for
positive values of the GB coupling $\alpha$ the shadow size of EGB black holes
is smaller than that for the Schwarzschild black hole, while these are larger
for negative $\alpha$. The size of the shadow for quintessence EGB black holes
with $q=0.1$ is shown in the bottom panels. We see that for positive $\alpha$
the shadow radius of quintessence EGB black holes is smaller than that of the
quintessence Schwarzschild black hole, while it is larger for the negative
ones. Moreover, it is easy to see that the shadow size of quintessence 4$D$
EGB black holes is larger than that of 4$D$ EGB black holes without
quintessence; namely the presence of the quintessence matter increases the
shadow radius of the black hole which is in agreement with Table 1.

Figure 12: Top-left panel and top-right panel the shadow cast by 4$D$ EGB
black holes without quintessence. Bottom-left panel and bottom-right panel the
shadow cast by 4$D$ EGB black holes with quintessence.

Figure 13:  The radius of shadow $R_{s}$. Left panel: for different values of
$q$ as a function of GB parameter $\alpha$. Right panel: for different values
of the GB coupling $\alpha$ as a function of $q$. The solid curves correspond
to 4$D$ EGB black holes without quintessence matter and quintessence
Schwarzschild black hole in the left and right panels, respectively.

Also, in the left panel of Figure 13 we have plotted the radius of shadow for
4$D$ EGB black holes surrounded by quintessence matter for different values of
quintessence parameter $q$. The solid red curve corresponds to 4$D$ EGB black
holes without quintessence, $q=0$. We see that for a fixed value of $\alpha$
by increasing the $q$ parameter the shadow radius increases. In the case of
$\alpha=0$ and $q=0$ we retain the radius of shadow for the Schwarzschild
black hole, i.e, $R_{s}=3\sqrt{3}$, the solid point on the red curve. The
right panel of the figure shows the radius of shadow for different values of
the GB coupling $\alpha$. The solid red curve corresponds to the quintessence
Schwarzschild black holes in the limit of $\alpha\rightarrow 0$. As can be
seen, for a fixed value of $q$ by increasing the $\alpha$ parameter the shadow
size decreases. Also, it is clear that the shadow radius for $\alpha<0$ is
larger than the quintessence Schwarzschild black hole, while it is smaller for
$\alpha>0$.

##  5 Deflection angle by quintessence 4$D$ EGB black holes

The problem of light bending in the space-time of 4$D$ EGB black holes has
been considered in [134]. Also, the deflection angle of charged massive
particles has been calculated in [144]. Now, in this section we are going to
calculate the bending angle of light for 4$D$ EGB black holes surrounded by
quintessence matter for $\omega_{q}=-0.35$. As we mentioned in the previous
section, when the light rays with $b>b_{c}$ approach to the black hole, will
be deflected at the turning point $r=r_{0}$, see Figure 7. So, the first step
is to find the closest approach $r_{0}$. To this end, we need to solve the
equation $\frac{dr}{d\varphi}=0$. From equation (24) we find

$\left(\frac{1}{r^{2}}\frac{d\varphi}{dr}\right)^{2}={\frac{1}{b^{2}}-\frac{1}{r^{2}}-\frac{1}{2\alpha}\left[1-\sqrt{1+\frac{8\alpha
M}{r^{3}}+\frac{8\alpha q}{r^{1.95}}}\right]}\equiv\xi(r).$ (29)

As is clear, the value of impact parameter has an important role in
determining the closest approach. In Figure 14, we have displayed $\xi(r)$ as
a function of $r$ for different values of the impact parameter. The figure
shows that the function $\xi(r)$ has two roots at $r=r_{0}$ and $r=r_{1}$, for
$b>b_{c}$. Also, according to the Figure 7 we see that $r_{0}>r_{c}$.
Therefore, the root at $r=r_{0}$ is considered as the closest approach.
Moreover, according to the figure, by increasing the impact parameter $b$ the
closest approach $r_{0}$ also increases.

Using equation (25) the bending angle is given by

$\delta=2\int_{0}^{u_{0}}\bigg{|}\frac{d\varphi}{du}\bigg{|}du-\pi=2\int_{0}^{u_{0}}\frac{1}{\sqrt{\frac{1}{b^{2}}-u^{2}-\frac{1}{2\alpha}\left[1-\sqrt{1+{8\alpha
Mu^{3}}+{8\alpha qu^{1.95}}}\right]}}du-\pi.$ (30)

In the right panel of Figure 15, the bending angle as a function of $u_{0}$ is
plotted for $q=0.1$ and different values of the GB coupling. It is easy to see
that in the limit $u_{0}\rightarrow 0$, the bending angle tends to zero, as
one would expect. As one decreases the $\alpha$ parameter, the effect of the
GB term as a candidate for dark energy also decreases, and thus the strength
of gravitational attraction increases leads to an increase of the bending
angle. Therefore, in the case of $\alpha\rightarrow 0$ which corresponds to
the Schwarzschild black hole surrounded by quintessence matter, we have the
largest values of the deflection angle. In order to see the effect of the $q$
parameter on the light bending, we have shown the bending angle as function of
$u_{0}$ for $\alpha=0.7$ and different values of $q$ in the left panel of
Figure 15. According to the figure, with decreasing $q$ the bending angle also
decreases, so that for $q\rightarrow 0$ which corresponds to the 4$D$ EGB
black holes, the bending angle has the minimum value.

Figure 14: The behavior of $\xi$ as a function of $r$ for different values of
impact parameter $b$ for $\alpha=0.3$, $q=0.1$ and $M=1$.

Figure 15: The bending angle versus the inverse of closet approach
$u_{0}=\frac{1}{r_{0}}$. Left panel: for $\alpha=0.3$ and different values of
$q$. Right panel: for $q=0.1$ and different values of the GB coupling
$\alpha$.

##  6 Conclusions

The exact $4D$ static spherically symmetric black hole solutions to EGB
gravity surrounded by quintessence matter have been derived in Ref. [142]. In
this paper, we considered these 4$D$ EGB black holes with quintessence matter
and studied the shadow and null geodesics around them and investigated the
effects of both the GB coupling constant and the quintessence parameter on the
motion of light rays.

It is shown that depending on the value of $\omega_{q}$, the metric of
quintessence 4$D$ EGB black holes could have one, two or three horizons. So,
we first presented the structure of the black hole event horizon for
$\omega_{q}=-0.35$ and then for this special case we studied both radial and
non-radial null geodesics in details. For non-radial geodesics, we have
numerically obtained the physical quantities such as the effective potential,
radius of the photon sphere and impact parameter for different values of the
GB coupling and the quintessence parameter, which the results summarized in
Table 1. We found that for a given value of $q$, with increasing the GB
coupling the values of the event horizon and the photon radius decrease. This
result is to be expected, because a positive GB coupling plays the role of
dark energy, counteract gravity and thus the event horizon and photon radius
take smaller values. However, it is shown that for a fixed value of $\alpha$,
by increasing the quintessence parameter $q$, both the event horizon and the
photon radius also increase. Therefore, the presence of the quintessence
matter around a 4$D$ EGB black hole increases the values of the event horizon
and the photon sphere. The achieved results are compared to that obtained for
4$D$ EGB black holes without quintessence [126] and to the Schwarzschild black
hole surrounded by quintessence matter [73].

Also, we have analyzed the null trajectories around quintessence 4$D$ EGB
black holes for different values of the impact parameter and we have
numerically plotted the unstable circular and unbounded orbits of photons in
details. Furthermore, we have studied the role of the GB coupling and the
quintessence parameter on the shadow size of quintessence 4$D$ EGB black holes
and found that the shape of the shadow is a perfect circle which its radius
increases with the quintessence parameter. However, for a given value of $q$,
the shadow radius decreases for positive GB coupling and increases for
negative ones. Finally, as an physical application of null geodesics we
calculated the deflection angle around 4$D$ EGB black holes with quintessence
and studied the effects of the GB coupling and the quintessence parameter on
it. The study of null geodesics of charged black holes in 4$D$ EGB gravity
surrounded by quintessence matter is the subject of the future investigations.

## References

  * [1] S. Bowyer, E. T. Byrma, T. A. Chubb and H. Friedman, Science 147 (1965) 394.
  * [2] B. L. Webster and P. Murdin, Nature 235 (1972) 37.
  * [3] C. T. Bolton, Nature 235 (1972) 271.
  * [4] B. P. Abbott et al., (Virgo and LIGO Scientific Collaborations), Phys. Rev. Lett 116 (2016) 061102.
  * [5] Gravity Collaboration, A and A 602 (2017) A94.
  * [6] K. Akiyama et al., AstroPhys. J 875 (2019) L1.
  * [7] K. Akiyama et al., AstroPhys. J 875 (2019) L5.
  * [8] K. Akiyama et al., AstroPhys. J 875 (2019) L6.
  * [9] J. L. Synge, Mon. Not. R. Astron. Soc 131 (1966) 463.
  * [10] J. P. Luminet, Astron. Astrophys 75 (1979) 228.
  * [11] J. M. Bardeen, in Black holes, in Proceeding of the Les Houches Summer School, Session 215239 edited by C. De Witt and B.S. De Witt and B.S. De Witt (Gordon and Breach, NewYork, 1973).
  * [12] L. Amarilla, E. F. Eiroa and G. Giribet, Phys. Rev. D 81 (2010) 124045.
  * [13] S. W. Wei and Y. X. Liu, JCAP 11 (2013) 063.
  * [14] L. Amarilla and E. F. Eiroa, Phys. Rev. D 87 (2013) 044057.
  * [15] B. P. Singh and S. G. Ghosh, Ann. Phys 395 (2018) 127.
  * [16] M. Amir, A. Banerjee and S. D. Maharaj, Ann. Phys 400 (2019) 198.
  * [17] G. Z. Babar, A. Z. Babar and F. Atamurotov, Eur. Phys. J. C 80 (2020) 761.
  * [18] S. Kala, Saurabh, H. Nandan and P. Sharma, Int. J. Mod. Phys A 35 (2020) 2050177.
  * [19] I. Banerjee, S. Chakraborty and S. SenGupta, Phys. Rev. D 101 (2020) 041301.
  * [20] X. X. Zeng and H. Q. Zhang, Eur. Phys. J. C 80 (2020) 1058.
  * [21] M. Khodadi and E. N. Saridakis, Phys. Dark. Univ 32 (2021) 100835.
  * [22] M. Khodadi, G. Lambiase and D. F. Mota, arXiv:2107.00834 [gr-qc].
  * [23] J. A. V. Campos, M. A. Anacleto, F. A. Brito and E. Passos, arXiv:2103.10659 [hep-th].
  * [24] R. Shaikh, S. Paul, Pritam Banerjee and T. Sarkar, arXiv:2105.12057 [gr-qc].
  * [25] M. Okyay and A. Övgün, arXiv:2108.07766 [gr-qc].
  * [26] F. Rahaman, Ksh. Newton Singh, R. Shaikh, T. Manna and S. Aktar, arXiv:2108.09930 [gr-qc].
  * [27] J. W. Moffat, The Shadow of the Black Hole, (Oxford University Press, 2020).
  * [28] P. V. P. Cunha and C. A. R. Herdeiro , Gen. Rel. Grav 50 (2018) 42.
  * [29] Y. Hagihara, J. Astron. Geophys 8 (1931) 67.
  * [30] C. Darwin, Proc. Roy. Soc. London A 249 (1959) 180.
  * [31] S. Chandrasekhar, The Mathematical Theory of black holes, (Oxford University Press, 1998).
  * [32] E. Hackmann and C. Lammerzahl, Phys. Rev. D 78, (2008) 024035; ibid Phys. Rev. Lett 100 (2008) 171101.
  * [33] E. Hackmann, V. Kagramanova, J. Kunz and C. Lammerzahl, Phys. Rev. D 78 (2008) 124018; Erratum-ibid. 79 (2009) 029901.
  * [34] B. O’Neill, The Geometry of Kerr black holes, (A K Peters, Wellesley, 1995).
  * [35] D. Pugliese, H. Quevedo and R. Ruffini, Phys. Rev. D 84 (2011) 044030.
  * [36] G. V. Kraniotis, Class. Quant. Grav 21 (2004) 4743.
  * [37] E. Hackmann, C. Lammerzahl, V. Kagramanova and J. Kunz, Phys. Rev. D 81 (2010) 044020.
  * [38] E. Hackmann and H. Xu, Phys. Rev. D 87 (2013) 124030.
  * [39] D. Pugliese, H. Quevedo and R. Ruffini, Phys. Rev. D 83 (2011) 104052.
  * [40] S. Grunau and V. Kagramanova, Phys. Rev. D 83 (2011) 044009.
  * [41] G. Eskin, Math. Phys 31 (2019) 1950021.
  * [42] M. Olivares, J. Saavedra, C. Leiva and J. R. Villanueva, Mod. Phys. Lett. A 26 (2011) 2923.
  * [43] S. Grunau, H. Neumann and S. Reimers, Phys. Rev. D 97 (2018) 044011.
  * [44] H. Demirchian, A. Nersessian, S. Sadeghian and M. M. Sheikh-Jabbari, Phys. Rev. D 97 (2018) 104004.
  * [45] S. Soroushfar, R. Saffari, J. Kunz and C. Lmmerzahl, Phys. Rev. D 92 (2015) 044010.
  * [46] D. Youm, Mod. Phys. Lett. A 16 (2001) 2371.
  * [47] W. S. Klën and C. Molina, Phys. Rev. D 102 (2020) 104051.
  * [48] A. E. Mosaffa, Phys. Rev. D 83 (2011) 124006.
  * [49] J. Chen and Y. Wang, Int. J. Mod. Phys. A 25 (2010) 1439.
  * [50] V. Enolskii, B. Hartmann, V. Kagramanova, J. Kunz, C. Lammerzahl and P. Sirimachan, Phys. Rev. D 84 (2011) 84011.
  * [51] S. W. Wei, J. Yang and Y. Liu, Phys. Rev. D 99 (2019) 104016.
  * [52] A. Dashupta, H. Nandan and S. Kar, Phys. Rev. D 79 (2009) 124004.
  * [53] T. Maki and K. Shiraishi, Class. Quant. Grav 11 (1994) 227.
  * [54] S. Fernando, D. Krug and C. Curry, Gen. Rel. Grav. 35 (2003) 1243.
  * [55] S. Fernando, Phys. Rev. D 85 (2012) 024033.
  * [56] M. Olivares and J. R. Villanueva, Eur. Phys. J. C 73 (2013) 2659.
  * [57] P. Pradhan, Int. J. Mod. Phys. D 24 (2015) 1550086.
  * [58] C. Blaga, Applications. Math 22 (2013) 41.
  * [59] S. Soroushfar, R. Saffari and E. Sahami, Phys. Rev. D 94 (2016) 024010.
  * [60] K. Flathmann and S. Grunau, Phys. Rev. D 92 (2015) 104027.
  * [61] P. A. Gonzalez, M. Olivares, E. Papantonopoulos, J. Saavedra and Y. Vasquez, Phys. Rev. D 95 (2017) 104052.
  * [62] J. R. Villanueva, F. Tapia, M. Molina and M. Olivares, Eur. Phys. J. C 78 (2018) 853.
  * [63] B. Hoseini, R. Saffari and S. Soroushfar, Class. Quant. Grav 34 (2017) 055004.
  * [64] A. S. Habibina and H. S. Ramadhan, Phys. Rev. D 101 (2020) 124036.
  * [65] S. Fernando, ISRN Math. Phys 2012 (2012) 869069.
  * [66] N. Bretón, Class. Quant. Grav 19 (2002) 601.
  * [67] M. A. A. Paula, L. C. S. Leite and L. C. B. Crispino, Phys. Rev. D 102 (2020) 104033.
  * [68] E. Hackmann, B. Hartmann, C. Laemmerzahl and P. Sirimachan, Phys. Rev. D 81 (2010) 064016.
  * [69] E. Hackmann, B. Hartmann, C. Laemmerzahl and P. Sirimachan, Phys. Rev. D 82 (2010) 044024.
  * [70] R. S. Kuniyal, R. Uniyal, A. Biswas, H. Nandan and K. D. Purohit, Int. J. Mod. Phys. A 33 (2018) 1850098.
  * [71] A. Mohammadi, T. Golanbari, B. Malakolkalami and S. Jalalzadeh, arXiv:2101.08842 [gr-qc].
  * [72] S. Fernando, S. Meadows and K. Reis, Int. J. Theo. Phys 54 (2015) 3634.
  * [73] S. Fernando, Gen. Rel. Grav 44 (2012) 1857.
  * [74] B. Malakolkalami and K. Ghaderi, Mod. Phys. Lett. A 30 (2015) 1550049.
  * [75] R. Uniyal, D. C. Devi, H. Nandan and K. D. Purohit, Gen. Rel. Grav 47 (2015) 16.
  * [76] K. Nozari and M. Hajebrahimi, arXiv:2004.14775 [gr-qc].
  * [77] A. Ramos, C. Arias, R. Avalos and E. Contreras, Ann. Phys 431 (2021) 168557.
  * [78] S. Hod, Phys. Rev. D 84 (2011) 124030.
  * [79] A. Mishra and S. Chakraborty, Eur. Phys. J. C 78 (2018) 374.
  * [80] I. Potashov, J. Tchemarina and A. Tsirulev, Universe 6 (2020) 183.
  * [81] F. Willenborg, S. Grunau, B. Kleihaus and J. Kunz, Phys. Rev. D 97 (2018) 124002.
  * [82] B. Zwiebach, Phys. Lett. B 156 (1985) 315.
  * [83] R. I. Nepomechie, Phys. Rev. D 32 (1985) 3201.
  * [84] P. Kanti, N. Mavromatos, J. Rizos, K. Tamvakis, and E. Winstanley, Phys. Rev. D 54 (1996) 5049.
  * [85] B. Kleihaus, J. Kunz and E. Radu, Phys. Rev. Lett 106 (2011) 151104.
  * [86] T. P. Sotiriou and S.Y. Zhou, Phys. Rev. Lett 112 (2014) 251102.
  * [87] D. D. Doneva and S. S. Yazadjiev, Phys. Rev. Lett 120 (2018) 131103.
  * [88] G. Antoniou, A. Bakopoulos and P. Kanti, Phys. Rev. Lett 120 (2018) 131102.
  * [89] L. G. Collodel, B. Kleihaus, J. Kunz and E. Berti, Class. Quant. Grav 37 (2020) 075018.
  * [90] D. Glavan and C. Lin, Phys. Rev. Lett 124 (2020) 081301.
  * [91] W. Y. Ai, Commun. Theor. Phys 72 (2020) 095402.
  * [92] M. Gurses, T. C. Sisman and B. Tekin, Eur. Phys. J. C 80 (2020) 647.
  * [93] S. Mahapatra, Eur. Phys. J. C 80 (2020) 992.
  * [94] F. W. Shu, Phys. Lett. B 811 (2020) 135907.
  * [95] S. X. Tian and Z. H. Zhu, arXiv:2004.09954 [gr-qc].
  * [96] J. Bonifacio, K. Hinterbichler and L. A. Johnson, Phys. Rev. D 102 (2020) 024029.
  * [97] J. Arrechea, A. Delhom and A. Jiménez-Cano, Chin. Phys. C 45 (2021) 013107.
  * [98] P. G. S. Fernandes, P. Carrilho, T. Clifton and D. J. Mulryne, Phys. Rev. D 102 (2020) 024025.
  * [99] R. A. Hennigar, D. Kubizňák, R. B. Mann and C. Pollack, JHEP 2020 (2020) 27.
  * [100] H. Lü and Y. Pang, Phys. Lett. B 809 (2020) 135717.
  * [101] T. Kobayashi, JCAP 07 (2020) 013.
  * [102] K. Aoki, M. A. Gorji and S. Mukohyama, Phys. Lett. B 810 (2020) 135843.
  * [103] P. G. S. Fernandes, Phys. Lett. B 805 (2020) 135468.
  * [104] R. Kumar and S. G. Ghosh, JCAP 2020 (2020) 053.
  * [105] S. W. Wei and Y. X. Liu, Eur. Phys. J. Plus 136 (2021) 436.
  * [106] R. A. Konoplya and A. Zhidenko, Phys. Rev. D 101 (2020) 084038.
  * [107] A. Casalino, A. Colléaux, M. Rinaldi and S. Vicentini, Phys. Dark. Univ 31 (2021) 100770.
  * [108] S. G. Ghosh and S. D. Maharaj, Phys. Dark. Univ 30 (2020) 100687.
  * [109] S. G. Ghosh and R. Kumar, Class. Quant. Grav 37 (2020) 245008.
  * [110] A. Kumar and S. G. Ghosh, arXiv:2004.01131 [gr-qc].
  * [111] A. Kumar and R. Kumar, arXiv:2003.13104 [gr-qc].
  * [112] R. A. Konoplya and A. F. Zinhailo, Eur. Phys. J. C 80 (2020) 1049.
  * [113] A. K. Mishra, Gen. Rel. Grav 52 (2020) 106.
  * [114] A. Aragón, R. Bécar, P. A. González and Y. Vásquez, Eur. Phys. J. C 80 (2020) 773.
  * [115] S. Devi, R. Roy and S. Chakrabarti, Eur. Phys. J. C 80 (2020) 760.
  * [116] R. A. Konoplya and A. Zhidenko, Phys. Dark. Univ 30 (2020) 100697.
  * [117] M. A. Cuyubamba, Phys. Dark. Univ 31 (2021) 100789.
  * [118] P. Liu, C. Niu and C. Y. Zhang, Chin. Phys. C 45 (2021) 2.
  * [119] S. A. Hosseini Mansoori, Phys. Dark. Univ 31 (2021) 100776.
  * [120] S. Ying, Chin. Phys. C 44 (2020) 125101.
  * [121] D. V. Singh and S. Siwach, Phys. Lett. B 808 (2020) 135658.
  * [122] S. J. Yang, J. J. Wan, J. Chen, J. Yang and Y. Q. Wang, Eur. Phys. J. C 80 (2020) 937.
  * [123] C. Liu, T. Zhu and Q. Wu, Chin. Phys. C 45 (2021) 015105.
  * [124] M. Heydari-Fard, M. Heydari-Fard and H. R. Sepangi, Eur. Phys. J. C 81 (2021) 473.
  * [125] M. Guo and P. C. Li, Eur. Phys. J. C 80 (2020) 588.
  * [126] X. X. Zeng, H. Q. Zhang and H. Zhang, Eur. Phys. J. C 80 (2020) 872.
  * [127] B. Eslam Panah, Kh. Jafarzade, S. H. Hendi, Nucl. Phys. B 961 (2020) 115269.
  * [128] Y. P. Zhang, S. W. Wei and Y. X. Liu, Universe 6 (2020) 103.
  * [129] S. Shaymatov, J. Vrba, D. Malafarina, B. Ahmedov and Z. Stuchlik, Dark Universe 30 (2020) 100648.
  * [130] J. Rayimbaev, A. Abdujabbarov, B. Turimov and Farruh Atamurotov, arXiv:2004.10031 [gr-qc].
  * [131] F. Atamurotov, S. Shaymatov, P. Sheoran and S. Siwach, arXiv:2105.02214 [gr-qc].
  * [132] S. U. Islam, R. Kumar and S. G. Ghosh, JCAP 09 (2020) 030.
  * [133] R. Kumar, S. U. Islam and S. G. Ghosh, Eur. Phys. J. C 80 (2020) 1128.
  * [134] M. Heydari-Fard, M. Heydari-Fard and H. R. Sepangi, EPL 133 (2021) 50006.
  * [135] X. H. Jin, Y. X. Gao and D. J. Liu, Int. J. Mod. Phys. D 29 (2020) 2050065.
  * [136] S. L. Li, P. Wu and H. Yu, arXiv:2004.02080 [gr-qc].
  * [137] G. Narain and H. Q. Zhang, arXiv:2005.05183 [gr-qc].
  * [138] K. Aoki, M.A. Gorji and S. Mukohyama, JCAP 09 (2020) 014.
  * [139] K. Aoki, M.A. Gorji, S. Mizuno and S. Mukohyama, JCAP 01 (2021) 054.
  * [140] S. Shahidi and N. Khosravi, arXiv:2105.02372 [gr-qc].
  * [141] M. Ghodsi, A. Behnamfard, S. Fakhry and J. T. Firouzjaee, arXiv:2107.00562 [gr-qc].
  * [142] H. Shah, Z. Ahmad and H. H. Shah, Phys. Lett. B 818 (2021) 136383.
  * [143] V. V. Kiselev, Class. Quant. Grav 20 (2003) 1187.
  * [144] Z. Li, Y. Duan and J. Jia, arXiv:2012.14226[gr-qc].


# Epitaxial Ferromagnetic Nanoislands of Cubic GdN in Hexagonal GaN

T. F. Kent  Department of Materials Science and Engineering, The Ohio State
University, Columbus, Ohio 43210, USA  J. Yang  Department of Materials
Science and Engineering, The Ohio State University, Columbus, Ohio 43210, USA
L. Yang  Department of Materials Science and Engineering, The Ohio State
University, Columbus, Ohio 43210, USA  M. J. Mills  Department of Materials
Science and Engineering, The Ohio State University, Columbus, Ohio 43210, USA
R. C. Myers  Department of Materials Science and Engineering, The Ohio State
University, Columbus, Ohio 43210, USA  Deparment of Electrical and Computer
Engineering, The Ohio State University, Columbus, Ohio 43210, USA

###### Abstract

Periodic structures of GdN particles encapsulated in a single crystalline GaN
matrix were prepared by plasma assisted molecular beam epitaxy. High
resolution X-ray diffractometery shows that GdN islands, with rock salt
structure are epitaxially oriented to the wurtzite GaN matrix. Scanning
transmission electron microscopy combined with in-situ reflection high energy
electron diffraction allows for the study of island formation dynamics, which
occurs after 1.2 monolayers of GdN coverage. Magnetometry reveals two
ferromagnetic phases, one due to GdN particles with Curie temperature of 70K
and a second, anomalous room temperature phase.

In this work, the epitaxial integration of discrete cubic GdN nanoparticles in
a continuous, high crystalline quality GaN matrix is reported. Although the
growth of coalesced epitaxial GdN films on III-nitrides by molecular beam
epitaxy (MBE)Scarpulla ; Natali  and by reactive ion sputtering on
AlNYoshitomi  has been previously reported, formation of discrete GdN islands
within a continuous, epitaxial III-nitride matrix has until now, yet to be
explored. Epitaxial growth of dissimilar crystal structures is in general met
with many challengesSands , but the rare earth pnictides (RE-Pn:EuO, ErAs.)
have been shown to grow well on zincblende III-V semiconductor
compoundsPalmstrom . Most widely studied has been the epitaxial integration of
semi-metallic ErAs in InGaAsKadow , which has resulted in high speed
photodetector and THz source applicationsGriebel . The mechanism of ErAs
embedded nanoparticle growth in GaAs has been proposedCrook  to proceed first
by incomplete layer formation of ErAs islands on the surface followed by
epitaxial lateral overgrowth of the uncovered GaAs.

---
Figure 1: (a) Cross sectional z-contrast HAADF STEM image of calibration
sample. (b) Atomic resolution STEM images for selected layers in the
calibration structure. The top image shows a GdN nanoparticle of clearly cubic
structure surrounded by a continuous GaN matrix. The bottom atomic resolution
data shows the normal wurtzite structure prior to island formation. ( c) RHEED
patterns taken during growth of GdN layer for different thicknesses showing
evolution of the surface reconstruction from 1$\times$2 to 2$\times$4
reconstruction and back to the 1$\times$2 with increasing layer thickness.
Figure 2: (a) Structural diagram of 50 period 10nm GaN/2.4ML GdN superlattice.
(b) Cross sectional z-contrast HAADF STEM image showing all periods of
superlattice heterostructure. (c) High resolution X-ray diffraction
$\omega$-2$\theta$ scan showing epitaxial orientation of the GdN (111) peak to
the wurtzite (0002) of GaN as well as superlattice fringes, indicating precise
layer thickness control.

The epitaxial integration of GdN with GaN is attractive for a number of
reasons. First, the dilute doping of III-Nitrides with Gd has attracted a
large amount of attention in recent years initially for its promise of
utilization of the intra-f-shell UV optical transitions of Gd in AlNGruber
and subsequently for the search of a room temperature dilute magnetic
semiconductor following the report of room temperature ferromagnetism in
Gd:GaNDhar ; Bedoya ; Davies . Gd is attractive for its magnetic properties,
possessing the most strongly correlated electronic structure of the
lanthanides, with 4f ground state of spin 7/2 . Devices for semiconductor
spintronics require efficient ferromagnetic spin injection and detection
layers, which currently are composed of either dilute magnetic
semicondutor(DMS) (GaMnAs) or metallic layers. No epitaxial spin injector is
currently available in the III-nitride materials system. For the preservation
of spin coherence through the device, interface and crystalline quality are
key considerations. Dilute Gd:GaN, though offering the promise of room
temperature ferromagnetism and realization of a nitride based DMS has proven
to be a difficult material to controlRoever  due to its poorly understood
defect mediated mechanism of ferromagnetism. GdN, in contrast is a well
understood classical ferromagnetGambino  with Tc around 70K. Furthermore,
unlike most other RE-Pn, which are well established to be semimetals, thin GdN
layers has been predictedMitra ; Duan ; Lambrecht  to be indirect gap
semiconductors, a claim which is consistent with recently reported absorption
features for thin GdN filmsYoshitomi . This leads to the possibility of a
controllable ferromagnetic semiconductor which can be epitaxially integrated
with GaN. In addition to intriguing magnetic properties, embedded GdN
nanoparticles in GaN could potentially function as carrier recombination
centers, giving rise to ultrafast photoconductivity in the same fashion as RE-
As particles in III-arsenides.

For the epitaxial structure of the matrix to remain single crystalline, the
layer coverage of the rock salt GdN must remain incomplete, allowing for
epitaxial laterall overgrowth of the surrounding matrix. This is due to the
lower symmetry of the rock salt structure (Fm$\overline{3}$m) than the host
(P63mc in the case of wurtzite). Complete films of GdN on wurtzite III-
Nitrides have been shown to epitaxial with the relationship
GdN[111]$||$GaN[0001] but containing two rotational variants due to crystal
symmetry considerations, which can be observed by an off-axis $\phi$ scan in
x-ray diffractometry and resulting in a polycrystalline overlayer of
GaNScarpulla .

Samples were prepared using the technique of plasma assisted molecular beam
epitaxy (PAMBE). In a Veeco GEN930 PAMBE system equipped with a Ga, Gd
effusion cell, and nitrogen plasma source. To study the GdN island formation
threshold, a calibration stack consisting of increasing effective thicknesses
of GdN are deposited from 0.2ML to 2.4ML in between 10nm GaN spacers on a GaN
buffer layer grown on an AlN on sapphire (KYMA) template at a substrate
temperature of 730C, beam equivalent pressure of 2$\times 10^{-5}$Torr and
III/V ratio of 2. During the period of GdN growth, the Ga shutter is closed,
meaning only Gd and N are being deposited. There is, however, a residual
amount of Ga present on the surface, due to growth of the GaN spacer under
metal rich conditions. To analyze the onset of GdN island formation, the
samples were characterized by cross-sectional, atomic resolution TEM using an
FEI Titan3 80-300 Probe-Corrected Monochromated (S)TEM, as can be seen Fig 1a.
Up to 1.2ML GdN, no change in the structure of the heavily Gd doped region is
observed, however at 1.2ML, discrete clusters of highly concentrated Gd atoms
appear. From the atomic resolution data shown in Fig.1b, cubic particles of
GdN are observed in the 1.4ML layer. From image analysis of the STEM
dataImageJ , the lattice parameter of cubic GdN in GaN is measured to be
4.8$\pm 0.1\text{\AA}$ and the nanoparticle size is roughly 2.6nm x 3.6nm.
After 1.2ML and up to 2.4ML, GdN particles with clearly cubic structure
surrounded by a hexagonal GaN matrix can be seen with the major change with
additional Gd deposition being increased lateral growth, suggesting that the
height of the nanoparticle is self limited and further growth will proceed by
lateral expansion of the GdN islands, which is similar to what has been
observed for Er-Pn in III-As nanoparticle structuresHanson .

Figure 3: Diamagnetic background corrected magnetization hysteresis loops for
low (a) and high (b) temperatures from SQUID magnetometry for in-plane and
out-of-plane film orientations relative to the applied field. Inserts show the
low field data and open nature of the loops. The low temperature scan clearly
shows the highly symmetric GdN magnetic phase with the expected saturation
magnetization of 7$\mu$B/Gd3+. (c). Diamagnetic background corrected, constant
field magnetization behavior with temperature after field cooling from 300K to
5K at 5T for in plane and out of plane orientation of the film with the
applied field.

During growth, the surface reconstruction was monitored using reflection high
energy electron diffraction (RHEED) operating at 10kV and cathode current of
1.4A, results are shown in Fig. 1c. For the first 0.4ML of GdN coverage, the
pattern is representative of the wurtzite Ga-face 1$\times$2 reconstruction.
After 0.5ML and until 1.2ML of GdN the pattern changes to a 2$\times$4
reconstruction. Past 1.2ML of coverage, the wurtzite 1$\times$2 pattern again
is visible, indicating a temporary change in the surface structure during
growth of the GdN layer.

After calibration of the GdN precipitation threshold, the heterostructure show
in Fig. 2a. consisting of a GaN buffer on an AlN template on sapphire (KYMA)
and a 50 period superlattice of alternating 10nm uid-GaN and 2.4ML GdN layers
was prepared under identical growth conditions as the calibration structure.
Cross sectional STEM images, shown in Fig.2b, using a Technai F20 operating in
HAADF imaging mode, which provides atomic number contrast, shows expected
discrete GdN particles in a GaN matrix with GaN spacer thickness of
11.8$\pm$0.4nm and GdN layer thickness of 5.6$\pm$0.3nm obtained from image
analysis.The structure of the sample was further characterized by high
resolution x-ray diffractometry using a Bruker D8 triple axis system.
Diffraction data shown in Fig. 2c. exhibits clear epitaxial orientation of the
GdN [111] to the wurtzite [0001]. Also visible are superlattice fringes,
indicating precise layer thickness control of the GaN spacing layers between
the GdN regions. From analysis of the superlattice fringes, the GaN spacer
thickness can be determined to be 10.98nm which is close to the value
determined from STEM. Furthermore, from the diffraction angle, we can
determine the lattice parameter of GdN to be 4.97Å, which is in very good
agreement with the value obtained from STEM of the nanoparticles and with
values for bulk GdNNatali .

The magnetic properties of the sample were analyzed by superconducting quantum
interference device (SQUID) magnetometry using a Quantum Design MPMS XL.
Results, depicted in Fig. 3 clearly show evidence of two distinct
ferromagnetic phases in the sample. The dominant phase at low temperature
(Fig. 3a) can be identified as rocksalt GdN due to a saturation magnetization
of nearly exactly 7$\mu_{B}$/Gd3+ (158.2 emu/cm3), the expected configuration
for GdN. The low remanent magnetization but correct saturation of
7$\mu_{\text{B}}$ is indicative that a large fraction of Gd is paramagnetic,
which is further supported by a temperature dependence containing both a mean
field like behavior with a T${}_{\text{c}}$ 70K but an additional 1/T
contribution. Samples were characterized in both the in-plane
($\vec{B}$$||$GaN [0001]) and out-of-plane ($\vec{B}$$\perp$GaN [0001])
configuration. The low temperature phase shows very little anisotropy which is
consistent with small particles of a cubic structure, which should be free
from the shape anisotropy of a fully coalesced film. The coercive field,
H${}_{\text{c}}$ is measured to be 363Oe for the in-plane configuration and
170Oe in the out-of-plane configuration, respectively.

Past the curie point of GdN, a second, weaker and anisotropic ferromagnetic
phase persists to room temperature. This anomalous phase is hypothesized to be
the result of interaction of the Gd with local point defects in the GaN matrix
and is of the same type as observed by Dhar, et. al.Dhar . It was previously
reportedDavies  that ferromagnetic films of dilutely doped Gd:GaN exhibit
anisotropy in their saturation magnetization between the in-plane and out of
plane orientations of the film. As observed in Fig 3b, the room temperature
phase exhibits anisotropy with M${}_{\text{s}}$=6.84emu/cm3 and
H${}_{\text{c}}$ = 100Oe for the out-of-plane configuration and
M${}_{\text{s}}$ =2.8emu/cm3, H${}_{\text{c}}$ = 30.6Oe for the in-plane
orientation of the film.

Measurement of the magnetization behavior with temperature, after cooling in a
5T field, is shown in Fig 3c. These data reveal a sharp decrease in the
magnetization with temperature up to 70K, the reported Curie point of
GdNScarpulla . After 70K and up to the highest temperature measured, 350K, a
residual amount of magnetization persists, again pointing to the possibility
of an anomalous room temperature phase. In the M vs. T data, anisotropy is
observed at low temperatures, which is consistent with the low field
anisotropy present in the 5K hysteresis scan. For the out-of-plane
configuration, a distinct knee is visible in the M vs. T scan which could be
due to error in the orientation of the film as mounted in the magnetometer,
causing signal from the in-plane configuration to contribute slightly to the
measured magnetization. Due to the sample mounting technique employed, the
out-of-plane configuration has a larger uncertainty in the absolute
orientation of the film.

In summary, we have extended the growth of embedded rare earth pnictide
nanoparticles in III-V semiconductors to the family of the III-nitrides.
Samples show clear rocksalt structure in cross sectional TEM above a threshold
value of 1.2ML GdN and x-ray diffractometry indicates epitaxial orientation of
the [111] direction of GdN to the wurtzite c-axis. Magnetic characterization
shows evidence of two magnetic phases, one due to the rocksalt ferromagnet GdN
with Curie temperature of 70K and a second, anisotropic phase whose
magnetization persists past room temperature. The room temperature phase is
hypothesized to be the same type of defect mediated ferromagnetism reported in
Gd:GaN and shows a prominent out-of plane easy axisDhar ; Bedoya ; Davies .

Funding provided by the Center for Emergent Materials under NSF Award Number
DMR-0820414 and the Institute for Materials Research at OSU under the
Interdisciplinary Materials Research Grants (IMRG) program. Jim O’Brien of
Quantum Design is thanked for valuable discussions about advanced SQUID
measurement techniques.

## References

  * (1) M.A. Scarpulla, C. S. Gallinat, S. Mack, J. S. Speck, A. C. Gossard, J. Crys. Grow. 311, 1239 (2009)
  * (2) F. Natali, N. O. V. Plank, J. Galipaud, B. J. Ruck, J. J. Trodahl, F. Semond, S. Sorieul, L. Hirsch, J. Crys. Growth 312 3583
  * (3) H. Yoshitomi, S. Kitayama, T. Kita, O. Wada, M. Fujisawa, H. Ohta, T. Sakurai, Phys. Rev. B 83, 155202 (2011)
  * (4) T. Sands, C. J. Palmstr$\o$m, J. P. Harbison, V. G. Keramidas, N. Tabatabaie, T. L. Cheeks, R. Ramesh, Y. Silberberg, Mat. Sci. Rep. v5, 3, 99-170 (1990)
  * (5) C. J. Palmstr$\o$m, Annu. Rev. Mater. Sci. 25: 389-415 (1995)
  * (6) C. Kadow, S. B. Fleischer, J. P. Ibbetson, J. E. Bowers, A. C. Gossard, J. W. Dong, C. J. Palmstr$\o$m, Appl. Phys. lett. 75, 3548 (1999)
  * (7) M. Griebel, J. H. Smet, D. C. Driscoll, J. Kuhl, C. A. Diez, N. Freytag, C. Kadow, A. C. Gossard, K. von Kiltzing, Nature Materials 2, 122 - 126 (2003)
  * (8) A. M. Crook, H. P. Nair. D. A. Ferrer. S. R. Bank, Appl. Phys. Lett. 99, 072120 (2011)
  * (9) J. B. Gruber, U. Vetter, H. Hofsäss, B. Zandi, M. F. Reid, Phys. Rev. B 69, 195202 (2004)
  * (10) C. Mitra and W. R. L. Lambrecht, Phys. Rev. B 78, 195203 (2008)
  * (11) W. R. Lambrecht, Phys. Rev. B. 62, 13538 (2000)
  * (12) R. P. Davis, B. P. Gila, C. R. Abernathy, S. J. Pearton, C. J. Stanton, Appl. Phys. Lett. 96, 212502 (2010)
  * (13) C. Duan, R. F. Sabiryanov, J. Liu, W. N. Mei, Pa. A. Dowben, J. R. Hardy, Phys. Rev. Lett. 94 237201 (2005)
  * (14) S. Dhar, O. Brandt, M. Ramsteiner, V. F. Sapega, K. H. Ploog, Phys. Rev. Lett. 94 037205 (2005)
  * (15) A. Bedoya-Pinto, J. Malindretos, M. Roever, D. D. Mai, A. Rizzi, Phys. Rev. B. 80, 195208 (2009)
  * (16) M. Roever, J. Malindretos, A. Bedoya-Pinto, A. Rizzi, C. Rauch, F. Tuomisto, Phys. Rev. B. 84, 081201(R) (2011)
  * (17) M. Hanson, Ph.D. thesis, University of California at Santa Barbara (2007).
  * (18) R. J. Gambino, T. R. McGuire, H. A. Alperin, S. J. Pickart, J. Appl. Phys. 41, 933 (1970)
  * (19) M. D. Abramoff, Biophotonics International 11 36, (2004)
  * (20) L. F. Schneemeyer, R. B. van Dover, E. M. Gyorgy, J. Appl. Phys. 61, 3543 (1987)