#
Hairy Kiselev Black Hole Solutions

Yaghoub<EMAIL_ADDRESS>.tr, Maxim Misyura
<EMAIL_ADDRESS>and Vitalii
<EMAIL_ADDRESS>

(a) Department of Mathematics, Faculty of Sciences, Bilkent University, 06800
Ankara, Turkey
(b) Department of High Energy and Elementary Particles Physics, Saint
Petersburg State University,
University Emb. 7/9, Saint Petersburg, 199034, Russia
(c) Physics department, Herzen state Pedagogical University of Russia, 48
Moika Emb.,
Saint Petersburg 191186, Russia
(d) SPB branch of SAO RAS, 65 Pulkovskoe Rd, Saint Petersburg 196140, Russia

In the realm of astrophysics, black holes exist within non-vacuum cosmological
backgrounds, making it crucial to investigate how these backgrounds influence
the properties of black holes. In this work, we first introduce a novel static
spherically-symmetric exact solution of Einstein field equations representing
a surrounded hairy black hole. This solution represents a generalization of
the hairy Schwarzschild solution recently derived using the extended
gravitational decoupling method. Then, we discuss how the new induced
modification terms attributed to the primary hairs and various background
fields affect the geodesic motion in comparison to the conventional
Schwarzschild case. Although these modifications may appear insignificant in
most cases, we identify specific conditions where they can be comparable to
the Schwarzschild case for some particular background fields.

Keyword: Gravitational decoupling, Kiselev black hole, hairy black hole,
cosmological fields, geodesics

##  1 Introduction

In 2019 the Event Horizon Telescope Collaboration unveiled the very first
image of a black hole located at the center of the massive elliptical galaxy
M87 [1, 2, 3]. More recently, scientists have successfully observed the shadow
of the supermassive black hole located in the center of our own galaxy [4].
These direct observations provide compelling evidence that black holes are not
merely abstract mathematical solutions of the Einstein field equations but
real astrophysical objects. Black holes possess a range of miraculous
properties. For instances, they allow for the extraction of energy from their
rotation and electric fields [5, 6, 7, 8]. In the vicinity of the black hole’s
event horizon, particles can possess negative energy [5, 9, 7, 10, 11, 12],
and black holes can even function as particle accelerators [14, 15, 16, 17,
18, 19, 20].

In the realm of astrophysics, black holes are not isolated objects and they
inhabit non-vacuum backgrounds. Some research has focused on investigating the
direct local effects of cosmic backgrounds on the known black hole solutions.
For instance, Babichev et al. [babi] have shown that in an expanding universe
by a phantom scalar field, the mass of a black hole decreases as a result of
the accretion of particles of the phantom field into the central black hole.
However, one notes that this is a global impact. To explore the local changes
in the spacetime geometry near the central black hole, one should consider a
modified metric that incorporates the surrounding spacetime. In this context,
an analytical static spherically symmetric solution to Einstein field
equations has been presented by Kiselev [21]. This solution generalizes the
usual Schwarzschild black hole to a nonvacuum background and is characterized
by an effective equation of state parameter of the surrounding field of the
black hole. Hence it can encompass a wide range of possibilities including
quintessence, cosmological constant, radiation and dust-like fields. Several
properties of the Kiselev black hole have been extensively investigated in the
literature [85-90]. Later, this solution has been generalized to the dynamical
Vaidya type solutions [22, 23, 24]. Such generalizations are well justified
due to the non-isolated nature of real-world black holes and their exitance in
non-vacuum backgrounds. Black hole solutions coupled to matter fields, such as
the Kiselev solution, are particularly relevant for the study of astrophysical
black holes with distortions [25, 26, 27, 28]. They also play a significant
role in investigating the no-hair theorem [29, 30, 31, 32]. This theorem
states that a black hole can be described only with three charges (i.e. mass
$M$, electric charge $Q$ and angular momentum $a$), and it relies on a crucial
assumption that the black hole is isolated, meaning that the spacetime is
asymptotically flat and free from other sources. However, real-world
astrophysical situations do not meet this assumption . For instances, one may
refer to black holes in binary systems, black holes surrounded by plasma, or
those accompanied by accretion disks or jets in their vicinity. Such
situations imply that a black hole may put on different types of wigs, and
hence the applicability of the standard no-hair theorem for isolated black
holes to these cases becomes questionable [30, 31, 32, 33, 34].

Recently, the minimal geometrical deformations [35, 36, 37] and the extended
gravitational decoupling methods [38, 39, 40] have been utilized to derive new
solutions from the known seed solutions of Einstein field equations. These
techniques have been particularly effective in investigating the violation of
the no-hair theorem, the emergence of novel types of hairy black holes, and
the exploration of alternative theories of gravity [65, 66, 67, 68, 69, 70,
71, 72] Using the extended gravitational decoupling method, Ovalle et all [41]
have introduced a generalization of a Schwarzschild black hole surrounded by
an anisotropic fluid and possesses primary hairs. This new solution has
motivated a substantial further research in generalizing this solution to
hairy Kerr [42], Vaidya and generalized Vaidya [43], regular hairy black holes
[44, 45] and many others. Indeed, the gravitational decoupling method
represents a novel and powerful tool for obtaining new solutions to the
Einstein equations.

In the present work, we introduce a novel class of exact solutions to the
Einstein field equations, which describe a surrounded hairy Schwarzschild
black hole. This solution serves as a generalization of the previously
obtained hairy Schwarzschild solution using the extended gravitational
decoupling method. Then, in order to analyze the properties of the solution,
we investigate the effect of the new modification terms, attributed to the
primary hairs and various surrounding fields, on the timelike geodesic motion.
Specifically, we compare the effects of modification terms to the conventional
Schwarzschild case. While these modifications may seem negligible in most
scenarios, we identify specific situations where they can be comparable to the
Schwarzschild case, particularly when specific surrounding fields are present.
This analysis sheds light on the significance of these modifications in
certain situations, providing insights into the behavior of geodesic motion
around real astrophysical black holes.

The structure of the present paper is as follows. In Section 2, we briefly
discuss the hairy Schwarzschild solution by the minimal geometrical
deformations and the extended gravitational decoupling method. In Section 3,
we solve the Einstein field equations in order to obtain the surrounded hairy
Schwarzschild black hole. In Section 4, we do analysis of the timelike
geodesic motion. In Section 5, we summarize the new findings and implications
of the study. The system of units $c=G=1$ will be used throughout the paper.

##  2 Gravitational decoupling and hairy Schwarzschild black hole

Gravitational decoupling method states that one can solve the Einstein field
equations with the matter source

$\tilde{T}_{ik}=T_{ik}+\Theta_{ik}\,,$ (1)

where $T_{ik}$ represents the energy-momentum tensor of a system for which the
Einstein field equations are

$G_{ik}=8\pi T_{ik}\,.$ (2)

The solution of the equations (2) is supposed to be known and represents the
seed solution. Then $\Theta_{ik}$ represents an extra matter sources which
causes additional geometrical deformations. The Einstein equations for this
new matter source are

$\bar{G}_{ik}=\alpha\Theta_{ik}\,,$ (3)

where $\alpha$ is a coupling constant and $\bar{G}_{ik}$ is the Einstein
tensor of deformed metric only. The gravitational decoupling method states
that despite of nonlinear nature of the Einstein equations, a straightforward
superposition of these two solutions (2) and (3)

$\tilde{G}_{ik}\equiv G_{ik}+\bar{G}_{ik}=8\pi
T_{ik}+\alpha\Theta_{ik}\equiv\tilde{T}_{ik}\,,$ (4)

is also the solution of the Einstein field equations.

Now, we briefly describe this method. Let us consider the Einstein field
equations,

$G_{ik}=R_{ik}-\frac{1}{2}g_{ik}R=8\pi T_{ik}\,.$ (5)

Let the solution of (5) is a static spherically-symmetric spacetime of the
form

$ds^{2}=-e^{\nu(r)}dt^{2}+e^{\lambda(r)}dr^{2}+r^{2}d\Omega^{2}\,.$ (6)

Here $d\Omega^{2}=d\theta^{2}+\sin^{2}\theta d\varphi^{2}$ is the metric on
unit two-sphere, $\nu(r)$ and $\lambda(r)$ are functions of $r$ coordinate
only, and they are supposed to be known. The metric (6) is termed as the seed
metric.

Now, we seek the geometrical deformation of (6) by introducing two new
functions $\xi=\xi(r)$ and $\eta=\eta(r)$ by:

$\displaystyle e^{\nu(r)}\rightarrow e^{\nu(r)+\alpha\xi(r)},$ $\displaystyle
e^{\lambda(r)}\rightarrow e^{\lambda(r)}+\alpha\eta(r).$ (7)

Here $\alpha$ is a coupling constant. Functions $\xi$ and $\eta$ are
associated with the geometrical deformations of $g_{00}$ and $g_{11}$ of the
metric (6), respectively. These deformations are caused by new matter source
$\Theta_{ik}$. If one puts $\xi(r)\equiv 0$ then the only $g_{11}$ component
is deformed, leaving $g_{00}$ unperturbed – this is known as the minimal
geometrical deformation. It has some drawbacks, for example, if one considers
the existence of a stable black hole possessing a well-defined event horizon
[37]. Deforming both $g_{00}$ and $g_{11}$ components is an arena of the
extended gravitational decoupling. One should note that gravitational
decoupling can lead to an energy exchange between two matter sources [64]. For
example, if one opts for gravitational decoupling of Vaidya spacetime, then
one can decouple the usual Vaidya spacetime without energy exchange. However,
in the generalized Vaidya spacetime, there is an energy exchange for arbitrary
mass function $M(v,r)$ [43].

Substituting (2) into (6), one obtains

$ds^{2}=-e^{\nu+\alpha\xi}dt^{2}+\left(e^{\lambda}+\alpha\eta\right)dr^{2}+r^{2}d\Omega^{2}\,.$
(8)

The Einstein equations for (8) as

$\tilde{G}_{ik}=8\pi\tilde{T}_{ik}=8\pi(T_{ik}+\Theta_{ik})\,,$ (9)

give

$\displaystyle
8\pi(T^{0}_{0}+\Theta^{0}_{0})=-\frac{1}{r^{2}}+e^{-\beta}\left(\frac{1}{r^{2}}-\frac{\beta^{\prime}}{r}\right),$
$\displaystyle
8\pi(T^{1}_{1}+\Theta^{1}_{1})=-\frac{1}{r^{2}}+e^{-\beta}\left(\frac{1}{r^{2}}+\frac{\nu^{\prime}+\alpha\xi^{\prime}}{r}\right),$
$\displaystyle
8\pi(T^{2}_{2}+\Theta^{2}_{2})=\frac{1}{4}e^{-\beta}\left(2(\nu^{\prime\prime}+\alpha\xi^{\prime\prime})+(\nu^{\prime}+\alpha\xi^{\prime})^{2}-\beta^{\prime}(\nu^{\prime}+\alpha\xi^{\prime})+2\frac{\nu^{\prime}+\alpha\xi^{\prime}-\beta^{\prime}}{r}\right),$
$\displaystyle e^{\beta}\equiv e^{\lambda}+\alpha\eta.$ (10)

Here the prime sign denotes the partial derivative with respect to the radial
coordinate $r$, and we have
$8\pi\left(T^{2}_{2}+\Theta^{2}_{2}\right)=8\pi\left(T^{3}_{3}+\Theta^{3}_{3}\right)$
due to the spherical symmetry.

From (2) one can define the effective energy density $\tilde{\rho}$, effective
radial and tangential $\tilde{P}_{r}$, $\tilde{P}_{t}$ pressures as

$\displaystyle\tilde{\rho}=-(T^{0}_{0}+\Theta^{0}_{0}),$
$\displaystyle\tilde{P}_{r}=T^{1}_{1}+\Theta^{1}_{1},$
$\displaystyle\tilde{P}_{t}=T^{2}_{2}+\Theta^{2}_{2}.$ (11)

From (2) one can introduce the anisotropy parameter $\Pi$ as

$\Pi=\tilde{P}_{t}-\tilde{P}_{r}\,,$ (12)

where if $\Pi\neq 0$ then it indicates the anisotropic behaviour of fluid
$\tilde{T}_{ik}$.

The equations (2) can be decoupled into two parts444One should remember that
it always works for $T_{ik}\equiv 0$ i.e. the vacuum solution and for special
cases of $T_{ik}$ if one opts for Bianchi identities
$\nabla_{i}T^{ik}=\nabla_{i}\Theta^{ik}=0$ with respect to the metric (8)
otherwice there is an energy exchange i.e.
$\nabla_{i}\tilde{T}^{ik}=0\Rightarrow\nabla_{i}T^{ik}=-\nabla_{i}\Theta^{ik}\neq
0$.: the Einstein equations corresponding to the seed solution (6) and the one
corresponding to the geometrical deformations. If we consider the vacuum
solution i.e. $T_{ik}\equiv 0$ \- Schwarzschild solution, then by solving the
Einstein field equations which correspond the geometrical deformations, one
obtains the hairy Schwarzschild solution [41]

$ds^{2}=-\left(1-\frac{2M}{r}+\alpha e^{-\frac{r}{M-\frac{\alpha
l}{2}}}\right)dt^{2}+\left(1-\frac{2M}{r}+\alpha e^{-\frac{r}{M-\frac{\alpha
l}{2}}}\right)^{-1}dr^{2}+r^{2}d\Omega^{2}\,,$ (13)

where $\alpha$ is the coupling constant, $l$ is a new parameter with length
dimension and associated with a primary hair of a black hole. Here $M$ is the
mass of the black hole in relation with the Schwarzschild mass $\mathcal{M}$
as

$M=\mathcal{M}+\frac{\alpha l}{2}\,.$ (14)

The impact of $\alpha$ and $l$ on the geodesic motion, gravitational lensing,
energy extraction and the thermodynamics has been studied in Refs.[52, 53, 54,
55, 56], and the influence of primary hair on quasinormal frequencies for
scalar, vector and tensor perturbation fields has been investigated in [73].

##  3 Surrounded hairy Schwarzschild black hole

Recently, the hairy Schwarzschild black hole has been introduced in [41] by
using the gravitational decoupling method. This solution in the Eddington-
Finkelstein coordinates takes the form

$ds^{2}=-\left(1-\frac{2M}{r}+\alpha e^{-\frac{r}{M-\frac{\alpha
l}{2}}}\right)dv^{2}+2\varepsilon dvdr+r^{2}d\Omega^{2}\,.$ (15)

Here $v$ is the advanced $(\varepsilon=+1)$ or retarded $(\varepsilon=-1)$
Eddington time. In this section, using the approach in [21, 57, 22], we obtain
the generalization of this solution representing a hairy Schwarzschild
solution surrounded by some particular fields motivated by cosmology as in the
following theorem.

Theorem: Considering the extended gravitational decoupling [39] and the
principle of additivity and linearity in the energy-momentum tensor [21] which
allows one to get correct limits to the known solutions, the Einstein field
equations admit the following solution in the Eddington-Finkelstein
coordinates

$ds^{2}=-\left(1-\frac{2M}{r}-\frac{N}{r^{3\omega+1}}+\alpha
e^{-\frac{2r}{2M-\alpha l}}\right)dv^{2}+2\varepsilon
dvdr+r^{2}d\Omega^{2}\,,$ (16)

where $M=\mathcal{M}+\frac{\alpha l}{2}$ in which $\mathcal{M}$ and
$\mathcal{M}$ are integration constants. The metric represents a surrounded
hairy Schwarzschild solution or equivalently hairy Kiselev solution. We
summarize our proof as follows.

Let us consider the general spherical-symmetric spacetime in the form

$ds^{2}=-f(r)dv^{2}+2\varepsilon dvdr+r^{2}d\Omega^{2}\,.$ (17)

The Einstein tensor components for the metric (17) are given by

$\displaystyle
G^{0}_{0}=G^{1}_{1}=\frac{1}{r^{2}}\left(f^{\prime}r-1+f\right)\,,$
$\displaystyle
G^{2}_{2}=G^{3}_{3}=\frac{1}{r^{2}}\left(rf^{\prime}+\frac{1}{2}r^{2}f^{\prime\prime}\right)\,,$
(18)

where the prime sign represents the derivative with respect to the radial
coordinate $r$. The total energy-momentum tensor should be a combination of
$\Theta_{ik}$ associated to the minimal geometrical deformations and $T_{ik}$
associated to the surrounding fluid as

$\tilde{T}_{ik}=\alpha\Theta_{ik}+T_{ik}\,.$ (19)

One should note that here we do not demand the fulfilment of the condition
$\Theta^{ik}_{;k}=T^{ik}_{;k}=0$. Instead, we demand that
$\tilde{T}^{ik}_{;k}=0$ which follows the Bianchi identity. The total energy-
momentum tensor $\tilde{T}_{ik}$ follows the same symmetries of the Einstein
tensor (3) for (17), i.e., $\tilde{T}^{0}_{0}=\tilde{T}^{1}_{1}$ and
$\tilde{T}^{2}_{2}=\tilde{T}^{3}_{3}$.

An appropriate general expression for the energy-momentum tensor $T_{ik}$ of
the surrounding fluid can be [21]

$\displaystyle T^{0}_{0}=-\rho(r)\,,$ (20) $\displaystyle
T^{i}_{k}=-\rho(r)\left[-\xi\left(1+3\zeta\right)\frac{r^{i}r_{k}}{r^{n}r_{n}}+\zeta\delta^{i}_{k}\right]\,.$

From the form of the energy-momentum tensor (20), one can see that the spatial
profile is proportional to the time component, describing the energy density
$\rho$ with arbitrary constants $\xi$ and $\zeta$ depending on the internal
structure of the surrounding fields. The isotropic averaging over the angles
results in

$<T^{i}_{k}>=\frac{\xi}{3}\rho\delta^{i}_{k}=P\delta^{i}_{k}\,,$ (21)

since we considered $<r^{i}r_{k}>=\frac{1}{3}{\delta^{i}}_{k}r_{n}r^{n}$.
Then, we have a barotropic equation of state for the surrounding fluid as

$P(r)=\omega\rho(r)~{},\,~{}~{}\omega=\frac{\xi}{3}\,,$ (22)

where $P(r)$ and $\omega$ are the pressure and the constant equation of the
state parameter of the surrounding field, respectively. Here, one notes that
the source $T_{ik}$ associated to the surrounding fluid should possess the
same symmetries in $\tilde{T}_{ik}$ because $\Theta_{ij}$ associated to the
geometrical deformations has the same symmetries as 555One should note that
hairy Schwarzschild solution is supported with an anisotropic fluid
$\Theta^{i}_{k}$
$\Theta^{0}_{0}=-\bar{\rho}\,,\Theta^{1}_{1}=\bar{P}_{r}\,,\Theta^{2}_{2}=\Theta^{3}_{3}=\bar{P}_{t}\,.$
(23) Where the non-vanishing parameter $\Pi=\bar{P}_{t}-\bar{P}_{r}$ indicates
on the anisotropic nature of the energy momentum tensor. So, in order to
satisfy the condition $\Theta^{0}_{0}=\Theta^{1}_{1}$ the anisotropic fluid
should be satisfied with the equation of the state $P_{r}=-\bar{\rho}$.

$\displaystyle\Theta^{0}_{0}=\Theta^{1}_{1}=-\bar{\rho},$
$\displaystyle\Theta^{2}_{2}=\Theta^{3}_{3}=\bar{P}_{t}.$ (24)

It means that ${T}^{0}_{0}={T}^{1}_{1}$ and ${T}^{2}_{2}={T}^{3}_{3}$. These
exactly provide the so-called principle of additivity and linearity considered
in [21] in order to determine the free parameter $\zeta$ of the energy-
momentum tensor $T_{ik}$ of surrounding fluid as

$\zeta=-\frac{1+3\omega}{6\omega}\,.$ (25)

Now, substituting (22) and (25) into (20), the nonvanishing components of the
surrounding energy-momentum tensor $T_{ik}$ become

$\displaystyle T^{0}_{0}=T^{1}_{1}=-\rho,$ $\displaystyle
T^{2}_{2}=T^{3}_{3}=\frac{1}{2}\left(1+3\omega\right)\rho\,.$ (26)

Now, we know the Einstein tensor components (3) and the total energy-momentum
tensor (19). Putting all these equations together, the
$G^{0}_{0}=\tilde{T}^{0}_{0}$ and $G^{1}_{1}=\tilde{T}^{1}_{1}$ give us the
following equation:

$\frac{1}{r^{2}}\left(f^{\prime}r-1+f\right)=-\rho-\alpha\bar{\rho}\,.$ (27)

Similarly, the $G^{2}_{2}=\tilde{T}^{2}_{2}$ and $G^{3}_{3}=\tilde{T}^{3}_{3}$
components yields

$\frac{1}{r^{2}}\left(rf^{\prime}+\frac{1}{2}f^{\prime\prime}r^{2}\right)=\frac{1}{2}\left(1+3\omega\right)\rho+\bar{P}\,.$
(28)

Thus, there are four unknown functions $f(r),\,\rho(r),\,\bar{\rho}(r)$ and
$\bar{P}$ that can be determined analytically by the differential equations
(27) and (28) with the following ansatz:

$f(r)=g(r)-\frac{\alpha l}{r}+\alpha e^{-\frac{2r}{2M-\alpha l}}\,.$ (29)

Then, by substituting (29) into (27) and (28) and using (3) one obtains the
following system of linear differential equations 666Here we apply the
Einstein equation $\hat{G}^{i}_{k}=\alpha\Theta^{i}_{k}$ to eliminate
$\tilde{\rho}$ and $\tilde{P}$. $\hat{G}^{i}_{k}$ is the Einstein tensor for
the spacetime  $ds^{2}=-\left(1-\frac{\alpha l}{r}+\alpha
e^{-\frac{2r}{2M-\alpha l}}\right)dv^{2}+2\varepsilon dvdr+r^{2}d\Omega^{2}.$
for unknowns $\rho(r)$ and $g(r)$

$\displaystyle\frac{1}{r^{2}}\left(g^{\prime}r-1+g\right)=-\rho,$
$\displaystyle\frac{1}{r^{2}}\left(rg^{\prime}+\frac{1}{2}g^{\prime\prime}r^{2}\right)=\frac{1}{2}\left(1+3\omega\right)\rho\,.$
(30)

This second order linear system can be integrated to give the metric function
$g(r)$ as

$g(r)=1-\frac{2\mathcal{M}}{r}-\frac{N}{r^{3\omega+1}}\,,$ (31)

and the energy density $\rho(r)$ of the surrounding field as

$\rho(r)=-\frac{3\omega N}{r^{3(\omega+1)}}\,.$ (32)

Here $\mathcal{M}$ and $N$ are constants of integration representing the
Schwarzschild mass and the surrounding field structure parameter,
respectively. By putting all these solutions together, we arrive at the
surrounded hairy Schwarzschild solution or equivalently hairy Kiselev solution
as

$ds^{2}=-\left(1-\frac{2M}{r}-\frac{N}{r^{3\omega+1}}+\alpha
e^{-\frac{2r}{2M-\alpha l}}\right)dv^{2}+2\varepsilon
dvdr+r^{2}d\Omega^{2}\,,$ (33)

where $M=\mathcal{M}+\frac{\alpha l}{2}$. From (32), one can see that the weak
energy condition demands that parameters $\omega$ and $N$ have different
signs.

##  4 Timelike geodesics

Considering the geodesic motion in spherically-symmetric spacetime, without
loss of generality, one can consider the equatorial plane
$\theta=\frac{\pi}{2}$. The geodesic equations for the metric (17) can be
obtained by varying the following action:

$S=\int\mathcal{L}d\tau=\frac{1}{2}\int\left(-f\dot{v}^{2}+2\varepsilon\dot{v}\dot{r}+r^{2}\dot{\varphi}^{2}\right)d\tau\,,$
(34)

where the dot sign means the derivative with respect to the proper time
$\tau$. The spacetime (33) is spherically-symmetric and hence, in addition to
the time-translation Killing vector $\frac{\partial}{\partial t}$, there
exists another Killing vector $\varphi^{i}=\frac{\partial}{\partial\varphi}$
and the corresponding conserved quantity, the angular momentum per mass, is
given by

$\varphi^{i}u_{i}=\frac{\partial\mathcal{L}}{\partial\dot{\varphi}}=r^{2}\dot{\varphi}=L\,.$
(35)

Taking into account (34) and (35), one obtains the following three geodesic
equations

$\dot{\varphi}=\frac{L}{r^{2}}\,,$ (36)
$-\frac{1}{2}f^{\prime}\dot{v}^{2}+r\dot{\varphi}^{2}-\varepsilon\ddot{v}=0\,,$
(37) $\varepsilon\ddot{r}=f\ddot{v}+f^{\prime}\dot{v}\dot{r}\,,$ (38)

where the prime sign denotes the derivative with respect to the radial
coordinate $r$. Substituting (36) into (37), one obtains

$f\ddot{v}=\frac{\varepsilon fL^{2}}{r^{3}}-\frac{1}{2}\varepsilon
ff^{\prime}\dot{v}^{2}\,.$ (39)

Now, by applying the timelike geodesic condition $g_{ik}u^{i}u^{k}=-1$ into
the equation above, we find

$f^{\prime}\dot{v}\dot{r}=-\frac{1}{2}\varepsilon
f^{\prime}+\frac{1}{2}\varepsilon ff^{\prime}-\frac{1}{2}\varepsilon
f^{\prime}\frac{L^{2}}{r^{2}}\dot{v}^{2}\,.$ (40)

Substituting the equation (40) into (38) we arrive at the following general
equation of motion in terms of the metric function $f$ for the radial
coordinate

$\ddot{r}=-\frac{1}{2}\left(1+\frac{L^{2}}{r^{2}}\right)f^{\prime}+f\frac{L^{2}}{r^{3}}\,.$
(41)

Hence, using the obtained metric function (33), one obtains the geodesic
equation in the form

$\displaystyle\ddot{r}$ $\displaystyle=$
$\displaystyle\left(-\frac{M}{r^{2}}+\frac{L^{2}}{r^{3}}-\frac{3ML^{2}}{r^{4}}\right)_{sch}$
(42)
$\displaystyle+\left(-\gamma\frac{N}{2r^{\gamma+1}}-\left(\gamma+2\right)\frac{NL^{2}}{2r^{\gamma+3}}\right)_{s}$
$\displaystyle+\left(\frac{\alpha}{2M-\alpha l}e^{-\frac{2r}{2M-\alpha
l}}+\frac{\alpha L^{2}}{\left(2M-\alpha l\right)r^{2}}e^{-\frac{2r}{2M-\alpha
l}}-\frac{\alpha L^{2}}{r^{3}}e^{-\frac{2r}{2M-\alpha l}}\right)_{h}\,,$

where $\gamma=3\omega+1$. From (42), one can observe the following interesting
points.

  1. 1.

The three terms in the first line are the same as that of the standard
Schwarzschild black hole in which the first term represents the Newtonian
gravitational force, the second term represents the repulsive centrifugal
force, and the third term is the relativistic correction of Einstein’s general
relativity, which accounts for the perihelion precession.

  2. 2.

The terms in the second line are new correction terms due to the presence of
the background field, which surrounds the hairy Schwarzschild black hole, in
which its first term is similar to the term of the gravitational potential in
the first brackets, while its second term is similar to the relativistic
correction of general relativity. Then, regarding (42) one realizes that for
the more realistic nonempty backgrounds, the geodesic equation of any object
depends strictly not only on the mass of the central object of the system and
the conserved angular momentum of the orbiting body, but also on the
background field nature. The new correction terms may be small in general, in
comparison to their Schwarzschild counterparts (the first and the third term
in the first brackets). However, one can show that there are possibilities
that these terms are comparable to them. One also can observe, by using the
equation (32), that for $\omega\in(-\frac{1}{3},\,0)$ the Newtonian
gravitational force is strengthened by corrections caused by the surrounding
field, on the other hand, for other values of $\omega$ the force is weakened.
If we consider the same question regarding the second term, which corresponds
to the relativistic correction of Einstein’s general relativity, then for
values $\omega\in(-1,\,0)$ the force is strengthened and this is while this
force is weakened for other values $\omega$. The surrounding fluid does not
have any contributions to the repulsive centrifugal force.

  3. 3.

The terms in the third line represent modifications by the primary hairs
$\alpha$ and $l$. The second term here corresponds to the relativistic
correction of Einstein’s general relativity. The third term here represents a
new correction by the primary hairs to the repulsive centrifugal force. One
can define the effective distance $D$ to find out where this force disappears
by relation $\frac{A_{1}}{A_{r}}\approx 1$ where $A_{r}$ is the Schwarzschild
black hole repulsive centrifugal force, and $A_{1}$ is the correction to this
force caused by primary hairs. So the distance is given by

$D=\left(M-\frac{\alpha l}{2}\right)\ln\alpha\,.$ (43)

Considering a minimal geometrical deformations, $\alpha$ must be negligible,
i.e., $\alpha\ll 1$. So according to (43), the correction caused by primary
hairs can weaken the repulsive centrifugal force but it cannot cancel it, and
hence this correction is negligible, in general. The first term in (42)
contributes a correction to the Newtonian potential. This can be seen using
the effective potential $V_{eff}(r)$. One can write the geodesic equations in
the form

$V_{eff}(r)=\Phi(r)+\frac{L^{2}}{2r^{2}}+\Phi(r)\frac{L^{2}}{r^{2}}\,,$ (44)

where $\Phi(r)$ is related to $g_{00}$ metric component via relation

$g_{00}=-\left(1+2\Phi\right)\,.$ (45)

By comparing this with (33), we come to the conclusion that

$\Phi(r)=-\frac{M}{r}+\frac{N}{2r^{3\omega+1}}-\alpha e^{-\frac{r}{2M-\alpha
l}}\,.$ (46)

Now, taking the derivative of $V_{eff}$ in (44) with respect to $r$

$\frac{d^{2}r}{d\tau^{2}}=-\frac{dV_{eff}}{dr}\,,$ (47)

we arrive at the equation of motion (42).

In order to better understand the nature of the solution obtained in (33), one
can consider the following two groups of forces and investigate their
behaviour for various set of surrounding fields and primary hair parameters.

$G\equiv\frac{M}{r^{2}}+\gamma\frac{N}{2r^{\gamma+1}}-\frac{\alpha}{2M-\alpha
l}e^{-\frac{2r}{2M-\alpha l}}\\\ ,$ (48)
$H\equiv\frac{3ML^{2}}{r^{4}}+\left(\gamma+2\right)\frac{NL^{2}}{2r^{\gamma+3}}-\frac{\alpha
L^{2}}{(2M-\alpha l)r^{2}}e^{-\frac{2r}{2M-\alpha l}}\,,$ (49)

where $G$ group represents the Newtonian gravitational force with its
modifications and $H$ group corresponds to the relativistic corrections of the
general relativity. One can ask for the possibilities if the new modifications
caused by surrounding fields and primary hairs can cancel the original forces
or change their effect, i.e ., change their sign. Hence, we are interested in
possible cases in which for set of parameters $\omega$, $\alpha$ and $l$, the
$G$ and $H$ functions are getting negligible values or they change their
signs. In the following subsections, we consider some specific fields
possessing particular equations of state motivated by cosmology.

However, we can note the following facts which we can derive from (49). Let us
consider the first two terms: for $-1<\omega<0$ these two terms are always
positive. However, the second term is negative for positive $\omega$, and we
can expect the sign change of $H$. Let us consider two particular cases:

  * •

The radiation $\omega=\frac{1}{3}$. In this case, $|N|\leq M^{2}$ and the
first two terms become negative in the region $0\leq r\leq 2M/3$, which is
inside the event horizon. Because the third term in (49) is negligible we can
conclude that $H$ is always positive outside the event horizon region.

  * •

The stiff fluid $\omega=1$. In this case we can put $N=-M^{4}$ then
$f(r=M)>0$. Thus, in this case the event horizon location at the radius is
less than $M$. However, the first two terms in (49) become negative at $r=M$
and $H<0$ outside the event horizon region.

###  4.1 Stiff Fluid

We begin our analysis of timelike geodesics with the surrounding fluid having
the average equation of the state of a stiff fluid as

$P=\rho\,\Leftrightarrow\,\omega=1.$ (50)

As mentioned previously, the presence of the surrounding field has a weakening
effect on the forces given by (48) and (49). From (32), one observes that $N$
must be negative to maintain a positive energy density for the surrounding
fluid. Our objective is to determine whether the corrections by the
surrounding field and primary hairs can cancel out the initial Schwarzschild
forces or potentially can change their sign, and thereby, altering the
direction of the forces.

In Figure 1(a), we plotted three curves corresponding the usual Schwarzschild,
Kiselev and hairy Kiselev black holes. We observe that the function $G$ for
the hairy Kiselev black hole is negligible but positive near the event horizon
$r=2\mathcal{M}$ for the given specific set of parameters. However, in the
case of purely Kiselev black hole (i.e., $\alpha=0$), the function $G$ is
negative in the interval $2\leq r\leq 2.15$. One notes that in the purely
Kiselev case, we have a naked singularity (NS) (i.e., $g_{00}\neq 0$).777For
this set of parameters $g_{00}$ is always negative, i.e., there are not
positive roots of the equation $g_{00}=0$ for $r\in(0,+\infty)$. On this
reason, we have concluded that $r=0$ represents a NS because the Kretschmann
scalar diverges at $r=0$. By NS we mean that $r=0$ singularity is not covered
with the event horizon. The question about future-directed non-space-like
geodesics, which terminate at this singularity in the past, has not been
considered within this paper.

Figure 1(b) shows that the function $G$ becomes negative in the vicinity of
the event horizon (i.e., in the region $2\leq r\leq 2.02$) for the hairy
Kiselev black hole for the set of parameters $N=-5.186,\,l=1.567$. To have a
bigger distance from the event horizon, where the function $G$ can become
negative, one should increase $|N|$ and $l$, however, in this case,
$\mathcal{M}\sim\alpha l/2$ and it will not anymore be a minimal geometrical
deformation in (15). So we can conclude that $G$ might be negative outside the
event horizon but only in its vicinity.

Figure 1(c) compares the function $H$ for the Schwarzschild, Kiselev and hairy
Kiselev cases for the values considered in Figure 1(b).

Figure 1: Plot (a) shows the function $G$ versus the distance $r$ for
$N=-4.972,\,l=1.514,\,\alpha=0.5$ and $\mathcal{M}=1$. Plot(b) shows the
function $G$ versus the distance $r$ for $N=-5.186,\,l=1.567,\alpha=0.5$ and
$\mathcal{M}=1$, the small picture shows the function $G$ of hairy Kiselev
black hole in the horizon vicinity. Plot (c) shows the function $H$ versus the
distance $r$ for $N=-5.186,\,l=1.567,\alpha=0.5$ and $\mathcal{M}=1$. The red,
blue, and green curves represent the Schwarzschild, Kiselev, and hairy Kiselev
cases, respectively.

In order to understand better the influence of a primary hair on a geodesic
motion we put $\alpha=0.1$ in order to consider bigger values of $l$. Figures
2(a) and 2(b) show how $G$ changes with different values of $l$ and $N$. One
can see that there are regions where it becomes negative. However, from these
pictures one cannot realize if they deal with a black hole or a naked
singularity. For this purpose one should impose the condition of existence of
an event horizon . The Figure 2(c) shows how $G$ changes in this case.

Figure 2: Plot(a) shows the function $G$ versus the parameters
$N\in[-7,-6.245],\,l\in[4,8]$ for $r=2.1,\,\alpha=0.1$ and $\mathcal{M}=1$.
Plot(b) shows the function $G$ versus the parameters
$N\in[-7,-4.367],\,l\in[4,8]$ for $r=2.5,\,\alpha=0.1$ and $\mathcal{M}=1$.
Plot(c) shows the function $G$ versus $l\in[4,8],\,r\in[2,3]$ for
$N\in[-6.183,-2.983],\,\alpha=0.1$ and $\mathcal{M}=1$. The event horizon,
located at $r=2\mathcal{M}$, follows the condition $N=-0.8l+1.6e^{-2}$.

###  4.2 Radiation

Here we consider the surrounding field having the average equation of state of
radiation field as

$P=\frac{\rho}{3}\,\Leftrightarrow\,\omega=\frac{1}{3}\,.$ (51)

In this case, the $N$ parameter must be negative, and akin to the previous
case, the surrounding radiation field and primary hairs weaken the forces in
(48) and (49).

Figure 3: Plot(a) shows the function $G$ versus the distance $r$ for
$N=-3.729,\,l=4,\,\alpha=0.5$ and $\mathcal{M}=1$. Plot(b) shows $G$ versus
$r$ for $N=-3.889,\,l=4.16,\,\alpha=0.5$ and $\mathcal{M}=1$, the small
picture shows the function $G$ of hairy Kiselev black hole in the horizon
vicinity. Plot(c) shows $H$ versus $r$ for $N=-3.889,\,l=4.16,\,\alpha=0.5$
and $\mathcal{M}=1$. The red, blue, and green curves correspond to the
Schwarzschild, Kiselev, and hairy Kiselev cases, respectively.

Figure 3(a) shows three curves in the pure Schwarzschild, Kiselev and hairy
Kiselev black holes for the parameter values $N=-3.729$ and $l=4$. For the
case of surrounding radiationlike field, one observes that the spacetime is
akin to the hairy Reissner-Nordstrom black hole such that the parameter $N$
plays the role of black hole’s electric charge, i.e. $N=-Q^{2}$. So, in purely
Reissner-Nordstrom case, the curve corresponds to the naked singularity
because $\mathcal{M}^{2}<Q^{2}$. In comparison to the stiff fluid case, one
notes that the parameters $l$ and $N$ are taken greater values to ensure that
the function $G$ is negligible.

In Figure 3(b), we plotted curves in order to show that hairs can affect the
geodesic motion and hence $G$ can become negative in the event horizon
vicinity (in the region $2\leq r\leq 2.042$). In this case, we set $N=-3.889$
and $l=4.16$. One can see that the smaller values of $\omega$ we take, the
bigger values of $l$ are required to ensure the negative values of $G$. For
example, if we take this value of $l$ (i.e., $l=4.16$), then in the case of
stiff fluid, we have $N=-15.557$ (we obtain this value by demanding that the
event horizon is located at $r=\mathcal{M}$), then the $G$ function is
negative in the region $2\leq r\leq 2.534$. Thus, one can see that the region,
where negative values of $G$ are possible, shrinks when $\omega$ tends to
zero. Figure 3(c) denotes the function $H$ with the values of $N$ and $l$ as
in the previous figure.

Similar to the stiff fluid case, we have several plots for $\alpha=0.1$.
Figures 4(a) and 4(b) show that $G$ becomes negative at the larger distances
in comparison to the stiff fluid case. This apparently contradicts our
previous statement that the smaller $\omega$ we consider, the region where $G$
becomes negative becomes smaller. However, one notes that this is a case of
the naked singularity because if one imposes an extra condition of the event
horizon existence, then for this case ($\alpha=0.1$) the $G$ function is
always positive outside the horizon as can be seen from Figure 4(c).

Figure 4: Plot(a) shows the function $G$ versus the parameters
$N\in[-7,-4],\,l\in[4,8]$ for $r=2.1,\,\alpha=0.1$ and $\mathcal{M}=1$.
Plot(b) shows the function $G$ versus the parameters $N\in[-7,-4],\,l\in[4,8]$
for $r=2.5,\,\alpha=0.1$ and $\mathcal{M}=1$. Plot(c) shows the function $G$
versus $r,\,l$ for $N\in[-1.546,-0.746],\,\alpha=0.1$ and $\mathcal{M}=1$. The
event horizon, located at $r=2\mathcal{M}$, must satisfy the condition
$N=-0.2l+0.4e^{-2}$.

###  4.3 Dust

For a dustlike field we have

$P=0\,\Leftrightarrow\,\omega=0\,,$ (52)

and we can show analytically that the function $G$ is positive near the event
horizon as follows. We have

$\frac{2M+N}{r}=1+\alpha e^{-\frac{r}{\mathcal{M}}}\,.$ (53)

Substituting this into (48) and considering the event horizon at
$r=2\mathcal{M}$, one obtains

$\frac{1}{4\mathcal{M}}-\frac{\alpha}{4\mathcal{M}e^{2}}>0\,.$ (54)

So, for physically relevant values of $\alpha,\,l$ and $N$, the function $G$
is positive outside the event horizon.

Figure 5(a) compares three curves of a hairy Kiselev black hole, purely
Kiselev when $\alpha=0$, and Schwarzschild case when $\alpha=0$ and $N=0$.
These curves are plotted for $l=0.5\,,N=-0.115$. Figure 5(b) is plotted for
the same values of black hole parameters and shows the behaviour of the
function $H$. For $\omega\geq 0$ the function $H$ is positive, and its
behaviour is shown in the Figure 5(c). For other values of $\omega$ we could
not find the condition (at small values of $\alpha$) where $H$ becomes
negative.

Figure 5: Plot(a) shows the function $G$ versus the distance $r$ for
$N=-0.115,\,l=0.5,\alpha=0.5$ and $\mathcal{M}=1$. Plot(b) shows the function
$H$ versus the distance $r$ for $N=-0.115,\,l=0.5,\,\alpha=0.5$ and
$\mathcal{M}=1$. The red, blue, and green curves correspond to the
Schwarzschild, Kiselev, and hairy Kiselev cases, respectively. Plot(c) shows
the function $H$ versus $r\,,l$ for the values
$N\in[-0.773,-0.373],\,l\in[4,8],\,r\in[2,5],\,\alpha=0.1$ and
$\mathcal{M}=1$. The event horizon, located at $r=2\mathcal{M}$, must satisfy
the condition $N=-0.1l+0.2e^{-2}$.

###  4.4 Quintessence

For a quintessencelike field, the equation of the state is

$P=-\frac{2}{3}\rho\,\Leftrightarrow\,\omega=-\frac{2}{3}\,.$ (55)

In this case, the parameter $N$ must be positive as one can see from (32). The
function $G$ can be negligible in the vicinity of the horizon only if either
$N$ or $L$ are negative. However, $G$ can take negative values but at large
distances from the event horizon. As can be shown from Figure 6(a) at values
$l=0.05\,,N=0.028$, the function $G$ for Kiselev black hole becomes negative
at $r>8.553$. The effect of $N$ and $\alpha$ on the function $H$ for these
values are negligible and they become considerable only at large distances, as
one can see from Figure 6(b) .

Figure 6: Plot(a) shows the function $G$ versus the distance $r$ for
$N=0.028,\,l=0.05,\alpha=0.5$ and $\mathcal{M}=1$. Plot(b) shows the function
$H$ versus the distance $r$ for $N=0.028,\,l=0.05,\,\alpha=0.5$ and
$\mathcal{M}=1$. The red, blue, and green curves correspond to the
Schwarzschild, Kiselev, and hairy Kiselev cases, respectively.

###  4.5 De Sitter background

In this case, the surrounded fluid has the effective equation of the state

$P=-\rho\,\Leftrightarrow\,\omega=-1\,.$ (56)

Like in the previous case, the parameter $N$ must be negative, and the
function $G$ must be positive near the event horizon.

Figure 7(a) shows that the function $G$ for $N=0.016\,,l=0.01$ becomes
negative for $r>3.841$. The function $H$ behaves very similar in all three
cases as can be seen in Figure 7(b).

Figure 7: Plot(a) shows the function $G$ versus the distance $r$ for
$N=0.016,\,l=0.01,\,\alpha=0.5$ and $\mathcal{M}=1$.Plot(b) shows the function
$H$ versus the distance $r$ for the same values of parameters. The red, blue,
and green curves correspond to the Schwarzschild, Kiselev, and hairy Kiselev
cases, respectively.

Figure 8: Plot (a): the dependence of the function $G$ on the $r\,,l$ for the
values $N=0.003..0.0008\,,l=0..0.2\,,r=2..25\,,\alpha=0.1$ and
$\mathcal{M}=1$. The event horizon, when it is located at $r=2\mathcal{M}$,
must satisfy the following condition $N=-0.0125l+0.025e^{-2}$. Plot (b): the
dependence of the function $G$ on the $r\,,l$ for the values
$N=-0.047..-0.097\,,l=4..8\,,r=2..16\,,\alpha=0.1$ and $\mathcal{M}=1$. The
event horizon, when it is located at $r=2\mathcal{M}$, must satisfy the
following condition $N=-0.0125l+0.025e^{-2}$.

Figure 8 shows the behaviour of $G$ at $\alpha=0.1$ and with an extra
condition of the event horizon existence. Here 8(a) is plotted for positive
cosmological constant as 8(b) for negative cosmological constant - anti-de-
Sitter case.

###  4.6 Phantom field

In general, the equation of the state a phantomlike field lies in the range
$\omega<-1$ [75, 76, 77, 78, 79, 80, 81]. In order to study the effect of a
phantom field, one can consider, as instance

$P=-\frac{4}{3}\rho\,\Leftrightarrow\,\omega=-\frac{4}{3}\,.$ (57)

The parameter $N$ must be positive, and as can be seen in Figure 9(a), the
function $G$ takes negative values at the region $r>3.056$ at
$l=0.05\,,N=0.007$. Figure 9(b) shows that for the same values of $l$ and $N$,
the function $H$ can be negative in the region $r>5.433$.

Figure 9: Plot(a) shows the function $G$ versus the distance $r$ for
$N=0.007\,,l=0.05,\alpha=0.5$ and $\mathcal{M}=1$. Plot(b) shows the function
$H$ versus the distance $r$ for $N=0.007\,,l=0.05,\alpha=0.5$ and
$\mathcal{M}=1$. The red, blue, and green curves correspond to the
Schwarzschild, Kiselev, and hairy Kiselev cases, respectively.

##  5 Conclusion

Inspired by the fact that black holes inhabit nonvacuum cosmological
backgrounds, we present a new solution to the Einstein field equations
representing a surrounded hairy Schwarzschild black hole. This solution takes
into account both the primary hair and surrounding fields (represented by an
energy-momentum tensor following the linearity and additivity condition [21]),
which affect the properties of the black hole. The effect of the corresponding
contributions on timelike geodesics are discussed. We find that the new
induced modifications can be considerable in certain cases. In particular, we
investigate how the specified surrounding fields and primary hairs affect the
Newtonian and perihelion precession terms. Our observations are as follows.

  * •

The surrounding fields with $-\frac{1}{3}<\omega<0$ contribute positively to
the Newtonian term, i.e ., strengthening the gravitational attraction.

  * •

The new corrections to the Newtonian term might be the same order or even
greater for all other cases if one considers a naked singularity
.888Considering the positive $\omega$, the weak energy condition demands
negative $N$ values. This restriction, for example, in the dust case requires
$|N|<2M$, otherwise, the metric function $f(r)$ is always positive for all
ranges of $r$ since all the being four terms are positive, and hence there is
no event horizon. In the case of the radiation, i.e., $\omega=\frac{1}{3}$,
the NS occurs if $M^{2}+N<0$ which requires large values of $|N|$. Hence one
observes that for bigger values of $|N|$, the function $|G|$ becomes bigger,
but this implies the violation of the condition required for the existence of
an event horizon.

  * •

In the case that the solution represents a black holes, new corrections can be
of the same order or even greater than the Newtonian term in the event horizon
vicinity for $\omega>0$.

  * •

For $\omega<-\frac{1}{3}$, i.e., for effectively repulsive fluids akin to dark
energy models, the correction terms dominate far from the event horizon and
mainly near the cosmological horizon.

The Schwarzschild black hole is an idealized vacuum solution, and it is
important to consider how it gets deformed in the presence of matter fields.
Another crucial factor to consider is the impact of the surrounding
environment, particularly the shadow of a black hole in the cosmological
background, which serves as a potential cosmological ruler [58]. The solution
presented in this work can be further investigated to study the shadow of a
hairy Schwarzschild black hole in various cosmological backgrounds in order to
find out how anisotropic fluid can affect the observational properties [74],
which is a plan of our upcoming investigations. It is worthwhile to mention
that applying the Newman-Janis [59] and Azreg Ainou [60, 61] algorithms one
can obtain the rotating version of the solution presented here. Also,
investigation of quasinormal modes, thermodynamics properties, accretion
process, and gravitational lensing of these solutions can help us to
understand better the nature of these objects.

The obtained hairy Kiselev solution has many potential uses in various
cosmological and astrophysical scenarios. It can be an arena for high-energy
phenomena. If one considers the centre of mass energy $E_{c.m.}$ of two
colliding particles in usual Schwarzschild spacetime, than the value is quite
limited and small [14]. However, two extra terms here might lead to the
existence of the innermost stable equilibrium point in the horizon vicinity
[15], which can lead to unbound centre-of-mass energy $E_{c.m.}$ of two
colliding particles. Another tool to distinguish the hairy Kiselev black hole
from the usual Schwarzschild one is to study its shadow properties. The shape
of the shadow is the same as in the Schwarzschild case due to the spherical
symmetry. However, the existence of four extra parameters
$\omega\,,N\,,\alpha\,,l$ have, surely, impact on its size and intensity [62].
The study of the planet’s motion is the way to define if a primary hair can
have an impact on its trajectory . As we have shown, extra terms can
drastically change particle motion. However, in a realistic astrophysical
situation, one should consider this motion near the black hole where
$\alpha\,,l\,,N$ has a large impact on the particle motion. Based on the
parameters and variables considered in this model, it seems that attempting to
test it within the Solar system would be futile. This is because any
additional terms are essentially insignificant beyond the surface of the Sun,
resulting in a prediction that would be indistinguishable from that which is
already predicted in Schwarzschild spacetime. Therefore, it may be more
beneficial to focus on the study of black hole vicinity where more noticeable
results can be achieved. Although it has not yet been observed, the Hawking
temperature and radiation may get also influenced by a primary hair [55, 56].
The Schwarzschild black hole possesses a negative heat capacity. Cosmological
fields and primary hair might lead to positive specific heat capacity and
phase transition [63]. All these are the topics of our future investigations.

Acknowledgments: V. Vertogradov and M. Misyura say thanks to grant NUM.
22-22-00112 RSF for financial support.

## References

  * [1] The Event Horizon Telescope Collaboration, First M87 Event Horizon Telescope results. I. The shadow of the supermassive black hole, Astrophys. J. Lett. 875 (2019) L1
  * [2] The Event Horizon Telescope Collaboration, First M87 Event Horizon Telescope results. II. Array and instrumentation, Astrophys. J. Lett. 875 (2019) L2
  * [3] The Event Horizon Telescope Collaboration, First M87 Event Horizon Telescope results. III. Data processing and calibration, Astrophys. J. Lett. 875 (2019) L3
  * [4] Akiyama, K. et al. [Event Horizon Telescope Collaboration]. First Sagittarius A* Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole in the Center of the Milky Way. _Astrophys. J. Lett._ 2022, _930_ , L12.
  * [5] R. Penrose and R. M. Floyd, Extraction of Rotational Energy from a Black Hole Nature Physical Science 229, 177 (1971).
  * [6] O. B. Zaslavskii, Energy extraction from extremal charged black holes due to the Banados-Silk-West effect Phys. Rev. D 86, 124039 (2012) [arXiv:1207.5209]
  * [7] Lucas Timotheo Sanches, Mauricio Richartz. Energy extraction from non-coalescing black hole binaries Phys. Rev. D 104, 124025(2021), arXiv:2108.12408 (gr-qc)
  * [8] Vitalii Vertogradov, Extraction energy from charged Vaidya black hole via the Penrose process 2023 Commun. Theor. Phys. 75 045404
  * [9] O. B. Zaslavskii, Negative energy states in the Reissner-Nordstrom metric Mod. Phys. Lett. A 36, 2150120 (2021), arXiv:2006.02189 [gr-qc].
  * [10] Grib A.A., Pavlov Yu.V., Vertogradov V.D. Geodesics with negative energy in the ergosphere of rotating black holes. Modern Physics Letters A. Vol. 29,Iss. 20. 2014- P. 14501-14510. [arXiv:1304.7360]
  * [11] V. Vertogradov, The Negative Energy in Generalized Vaidya Spacetime Universe 6(9), 155 (2020), arXiv:2209.10976 [gr-qc]
  * [12] Vertogradov V.D. Geodesics with negative energy in the ergosphere of rotating black holes, Gravitation and Cosmology. Vol. 21, Iss. 2. 2015.- PP. 171-174. [arXiv:2210.04674]
  * [13] E. Babichev, V. Dokuchaev, Yu. Eroshenko, Black Hole Mass Decreasing due to Phantom Energy Accretion Phys. Rev. Lett. 93, 021102 (2004).
  * [14] M. Banados, J. Silk, S.M. West, Kerr Black Holes as Particle Accelerators to Arbitrarily High Energy Phys. Rev. Lett. 103, 111102 (2009) [arXiv:0909.0169].
  * [15] O. B. Zaslavskii, Acceleration of particles by black holes as a result of deceleration: Ultimate manifestation of kinematic nature of BSW effect Phys. Lett. B 712 (2012) 161 [arXiv:1202.0565]
  * [16] O. B. Zaslavskii, Circular orbits and acceleration of particles by near-extremal dirty rotating black holes: general approach. Class. Quantum Grav. 29 (2012) 205004 [arXiv:1201.5351 [gr-qc]]
  * [17] Grib, A. A., Pavlov Y. V. (2020) Rotating black holes as sources of high energy particles. Physics of Complex Systems, 1 (1), 40-49.DOI: 10.33910/2687-153X-2020-1-1-40-49
  * [18] Vertogradov, V.D. (2023) On particle collisions during gravitational collapse of Vaidya spacetimes. Physics of Complex Systems, 4 (1), 17-23.
  * [19] Mandar Patil, Tomohiro Harada, Ken-ichi Nakao, Pankaj S. Joshi, Masashi Kimura, Infinite efficiency of collisional Penrose process: Can over-spinning Kerr geometry be the source of ultra-high-energy cosmic rays and neutrinos? Phys. Rev. D 93, 104015 (2016) [arXiv:1510.08205 [gr-qc]]
  * [20] T. Harada, H. Nemoto, U. Miyamoto, Upper limits of particle emission from high-energy collision and reaction near a maximally rotating Kerr black hole Phys.Rev .D86:024027,2012 [arXiv:1205.7088]
  * [21] V. V. Kiselev, Class. Quintessence and black holes Quant. Grav. 20, 1187 (2003).
  * [22] Heydarzade, Y., Darabi, F. Surrounded Vaidya solution by cosmological fields. _Eur. Phys. J. C_ 2018, _78_ , 582.
  * [23] Heydarzade, Y.; Darabi, F. Surrounded Bonnor-Vaidya solution by cosmological fields. _Eur. Phys. J. C_ 2018, _78_ , 1004.
  * [24] Y. Heydarzade, F. Darabi Surrounded Vaidya black holes: apparent horizon properties, Eur. Phys. J. C (2018) 78:342.
  * [25] R. Geroch and J.B. Hartle, Distorted black holes J. Math. Phys. 23, 680 (1982).
  * [26] S. Fairhurst, B. Krishnan, Distorted Black Holes with Charge Int. J. Mod. Phys. D10, 691 (2001).
  * [27] S.R. Brandt, E. Seidel, Evolution of distorted rotating black holes. III. Initial data Phys. Rev. D54, 1403 (1996).
  * [28] M. Ansorg, D. Petroff, Black holes surrounded by uniformly rotating rings Phys. Rev. D 72.2, 024019 (2005).
  * [29] S. W. Hawking, Black holes in general relativity Commun. Math. Phys. 25, 152 (1972).
  * [30] J.D. Brown and V. Husian, Black holes with short hair Int. J. Mod. Phys. D6, 563 (1997).
  * [31] S.Droz, M. Heusler, N. Straumann, New black hole solutions with hair Phys. Lett. B, 268 (3-4), 371 (1991).
  * [32] J. Barranco, A. Bernal, J.C. Degollado, A. Diez-Tejedor, M. Megevand, M. Alcubierre, D. Ndnez, O. Sarbach, Schwarzschild Black Holes can Wear Scalar Wigs Phys. Rev. Lett. 109, 081102 (2012).
  * [33] J.D. Bekenstein, Novel ”no-scalar-hair” theorem for black holes Phys. Rev. D 51(12), R6608 (1995).
  * [34] Hawking, S.W.; Perry, M.J.; Strominger, A. Soft Hair on Black Holes. _Phys. Rev. Lett._ 2016, _116_ , 231301.
  * [35] J. Ovalle, Extending the geometric deformation: New black hole solutions. Int. J. Mod. Phys. Conf. Ser., 41, 1660132 (2016) [arXiv:1510.00855 [gr-qc]]
  * [36] Roberto Casadio, Jorge Ovalle, Roldao da Rocha, The Minimal Geometric Deformation Approach Extended. Class. Quantum Grav. 32 (2015) 215020 [arXiv:1503.02873 [gr-qc]]
  * [37] Ovalle, J.; Casadio, R.; Rocha, R.D.; Sotomayor, A.; Stuchlik, Z. Black holes by gravitational decoupling. _Eur. Phys. J. C_ 2018, _78_ , 960.
  * [38] Ovalle, J. Decoupling gravitational sources in general relativity: From perfect to anisotropic fluids. _Phys. Rev._ 2017, _D95_ , 104019.
  * [39] Ovalle, J. Decoupling gravitational sources in general relativity: The extended case. _Phys. Lett. B_ 2019, _788_ , 213.
  * [40] Contreras, E.; Ovalle, J.; Casadio, R. Gravitational decoupling for axially symmetric systems and rotating black holes. _Phys. Rev. D_ 2021, _103_ , 044020.
  * [41] Ovalle, J.; Casadio, R.; Contreras, E.; Sotomayor, A. Hairy black holes by gravitational decoupling. _Phys. Dark Universe_ 2021, _31_ , 100744.
  * [42] S. Mahapatra and I. Banerjee, Rotating hairy black holes and thermodynamics from gravitational decoupling, Phys. Dark Univ. 39 (2023) 101172 [arXiv: 2208.05796],
  * [43] Vitalii Vertogradov, Maxim Misyura ”Vaidya and Generalized Vaidya Solutions by Gravitational Decoupling”Universe 2022, 8(11), 567; doi:10.3390/universe8110567 [arXiv:2209.07441 [gr-qc]]
  * [44] Vitalii Vertogradov, Maxim Misyura, The Regular Black Hole by Gravitational Decoupling. Phys. Sci. Forum 2023, 7(1), 27
  * [45] Jorge Ovalle, Roberto Casadio, Andrea Giusti, Regular hairy black holes through Minkowski deformation. [arXiv:2304.03263 [gr-qc]]
  * [46] M. Jamil, S. Hussain, Dynamics of particles around a Schwarzschild-like black hole in the presence of quintessence and magnetic field B. Majeed, Eur. Phys. J. C 75. 1, 24 (2015).
  * [47] I. Hussain, S. Ali, Marginally stable circular orbits in the Schwarzschild black hole surrounded by quintessence matter Eur. Phys. J.l Plus 131. 8, 275 (2016).
  * [48] B. Malakolkalami, K. Ghaderi, The null geodesics of the Reissner-Nordstrom black hole surrounded by quintessence Mod. Phys. Lett . A 30 .10, 1550049 (2015).
  * [49] S. Fernando, Schwarzschild black hole surrounded by quintessence: Null geodesics Gen. Rel. Grav. 44.7, 1857 (2012).
  * [50] R. Uniyal, N.C. Devi, Geodesic Motion in Schwarzschild Spacetime Surrounded by Quintessence H. Nandan, K.D. Purohit, Gen. Rel. Grav. 47.2, 16 (2015).
  * [51] S. Fernando, S. Meadows, K. Reis, Null Trajectories and Bending of Light in Charged Black Holes with Quintessence Int. J. Theo. Phys. 54.10, 3634 (2015).
  * [52] Ramos, A.; Arias, C.; Avalos, R.; Contreras,E. Geodesic motion around hairy black holes. _Annals Phys._ 2021, _431_ , 168557\.
  * [53] Sohan Kumar Jha, Anisur Rahaman, Gravitational lensing by the hairy Schwarzschild black hole. arXiv:2205.06052[gr-qc]
  * [54] Zhen Li, Faqiang Yuan, Energy extraction via Comisso-Asenjo mechanism from rotating hairy black hole. [arXiv:2304 .12553 [gr-qc]]
  * [55] Cavalcanti, R.T.; Alves, K.d.S.; da Silva, J.M.H. Near horizon thermodynamics of hairy black holes from gravitational decoupling. _Universe_ 2022, _8_ , 363.
  * [56] Vertogradov V. D., Kudryavcev D. A. On the temperature of hairy black holes. Physics of Complex Systems. Vol. 4, no. 2, 2023
  * [57] Y. Heydarzade, F. Darabi, Black Hole Solutions Surrounded by Perfect Fluid in Rastall Theory Phys. Lett. B, 771, 365 (2017).
  * [58] Oleg Yu. Tsupko, Zuhui Fan, Gennady S. Bisnovatyi-Kogan, Black hole shadow as a standard ruler in cosmology. Classical and Quantum Gravity, 37, 065016 (2020) [arXiv:1905.10509 [gr-qc]]
  * [59] E. T. Newman and A. I. Janis, ”Note on the Kerr spinning particle metric,” 10.1063/1 J. Math. Phys. 6 (1965) 915-917.
  * [60] M. Azreg-Ainou, ”Regular and conformal regular cores for static and rotating solutions,” 10.1016/j.physletb.2014.01 Phys. Lett. B 730 (2014) 95-98, [arXiv:1401.0787 [gr-qc]].
  * [61] M. Azreg-Ai’nou, ”From static to rotating to conformal static solutions: Rotating imperfect fluid wormholes with(out) electric or magnetic field,” epjc/s10052-014-2865-8 Eur. Phys. J. C 74 no. 5, (2014) 2865, [arXiv:1401.4292 [gr-qc]].
  * [62] Rajibul Shaikh and others, Shadows of spherically symmetric black holes and naked singularities, Monthly Notices of the Royal Astronomical Society, Volume 482, Issue 1, January 2019, Pages 52-64
  * [63] Saleh, M., Thomas, B.B. Kofane, T.C. Thermodynamics and Phase Transition from Regular Bardeen Black Hole Surrounded by Quintessence. Int J Theor Phys 57, 2640-2647 (2018). https://doi.org/10 .1007/s10773-018-3784-5
  * [64] Ovalle, J., Contreras, E. Stuchlik, Z. Energy exchange between relativistic fluids: the polytropic case. Eur. Phys. J. C 82, 211 (2022). https://doi.org/10.1140/epjc/s10052-022-10168-5
  * [65] Gabbanelli, L., Rincón, Á. Rubio, C. Gravitational decoupled anisotropies in compact stars. Eur. Phys. J. C 78, 370 (2018). https://doi.org/10.1140/epjc/s10052-018-5865-2
  * [66] C. Deffayet, S. Deser, and G. Esposito-Farèse, Arbitrary $p$-form Galileons Phys. Rev. D 82, 061501 (R)
  * [67] Ángel Rincón and Grigoris Panotopoulos , Quasinormal modes of scale dependent black holes in $(1+2)$ -dimensional Einstein-power-Maxwell theory Phys. Rev. D 97, 024027
  * [68] Rincón, Á., Contreras, E., Bargueño, P. et al. Scale-dependent three-dimensional charged black holes in linear and non-linear electrodynamics. Eur. Phys. J. C 77, 494 (2017). https://doi.org/10 .1140/epjc/s10052-017-5045-9
  * [69] Benjamin Koch et al, A scale dependent black hole in three-dimensional space-time, 2016 Class. Quantum Grav. 33 225010
  * [70] Rebecca Briffa et al, Constraining teleparallel gravity through Gaussian processes, 2020 Class. Quantum Grav. 38 055007
  * [71] Grigoris Panotopoulos, Ángel Rincón, and Ilídio Lopes, Orbits of light rays in scale-dependent gravity: Exact analytical solutions to the null geodesic equations Phys. Rev. D 103, 104040
  * [72] Grigoris Panotopoulos, Ángel Rincón, Orbits of light rays in $(1+2)$-dimensional Einstein-power-Maxwell gravity: Exact analytical solution to the null geodesic equations, Annals of Physics, Volume 443, 2022, 168947 https://doi.org/10.1016/j.aop.2022.168947
  * [73] Avalos, R., Bargueño, P., Contreras, E., A Static and Spherically Symmetric Hairy Black Hole in the Framework of the Gravitational Decoupling. Fortschr. Phys. 2023, 71, 2200171. https://doi.org/10.1002/prop.202200171
  * [74] Sunny Vagnozzi et al, Horizon-scale tests of gravity theories and fundamental physics from the Event Horizon Telescope image of Sagittarius A, 2023 Class. Quantum Grav. 40–165007
  * [75] Parampreet Singh, M. Sami, and Naresh Dadhich, Cosmological dynamics of a phantom field, Phys. Rev. D 68, 023522
  * [76] Dzhunushaliev V., Folomeev V., Myrzakulov K., Myrzakulov R., Phantom fields: bounce solutions in the early Universe and s-branes, International Journal of Modern Physics DVol. 17, No. 12, pp. 2351–2358 (2008)
  * [77] Chakraborty, S., Debnath, U. The Effects of Tachyonic and Phantom Fields in the Intermediate and Logamediate Scenarios of the Anisotropic Universe. Int J Theor Phys 51, 1224-1238 (2012). https://doi.org/10.1007/s10773-011-0998-1.
  * [78] G.W. Gibbons, Phantom Matter and the Cosmological Constant, arXiv:hep-th/0302199.
  * [79] I. Maor, R. Brustein, J. McMahon, and P.J. Steinhardt, Measuring the equation of state of the universe: Pitfalls and prospects, Phys. Rev. D 65, 123003 (2002).
  * [80] A. E. Schulz, M. White, The Tensor to Scalar Ratio of Phantom Dark Energy Models, Phys. Rev. D 64, 043514 (2001).
  * [81] R.R. Caldwell, A Phantom Menace. Cosmological consequences of a dark energy component with super-negative equation of state, Phys. Lett. B 545, 23 (2002).


# Physics-informed Gaussian process model for Euler-Bernoulli beam elements

Gledson Rodrigo Tondo⋆
Bauhaus-Universität Weimar
Weimar, Germany
&Sebastian Rau
Bauhaus-Universität Weimar
Weimar, Germany
&Igor Kavrakov
Bauhaus-Universität Weimar
Weimar, Germany
&Guido Morgenthal
Bauhaus-Universität Weimar
Weimar, Germany

###### Abstract

A physics-informed machine learning model, in the form of a multi-output
Gaussian process, is formulated using the Euler-Bernoulli beam equation. Given
appropriate datasets, the model can be used to regress the analytical value of
the structure’s bending stiffness, interpolate responses, and make
probabilistic inferences on latent physical quantities. The developed model is
applied on a numerically simulated cantilever beam, where the regressed
bending stiffness is evaluated and the influence measurement noise on the
prediction quality is investigated. Further, the regressed probabilistic
stiffness distribution is used in a structural health monitoring context,
where the Mahalanobis distance is employed to reason about the possible
location and extent of damage in the structural system. To validate the
developed framework, an experiment is conducted and measured heterogeneous
datasets are used to update the assumed analytical structural model.

_K_ eywords Gaussian process $\cdot$ physics-informed $\cdot$ machine learning
$\cdot$ stiffness regression $\cdot$ structural health monitoring $\cdot$
model updating

##  1 Introduction

Machine learning has been extensively applied in structural engineering,
especially in the structural health monitoring (SHM) field, as the
availability of data collected from sensors increases [1]. The heterogeneity
of the collected datasets along with the knowledge of physical relations
between them, usually represented as partial differential equations (PDEs),
have motivated the recent use of physics-informed machine learning models to
extract meaningful information from measured data. In these algorithms, the
governing PDE is built into the machine learning model, effectively
integrating measurements and mathematical models [2]. Gaussian processes (GP)
[3] have been extensively used for such a task, as they offer a non-parametric
probabilistic view of the modelling scheme, which is usually framed in a
Bayesian manner [4, 5]. Particular applications of the physics-informed GPs
have been observed as an analytical model of the latent curvature in a sleeper
beam and the inverse problem of identifying the Reynolds number of a CFD
simulation [6, 7]. In this paper, a physics-informed Gaussian process model
for the Euler-Bernoulli beam formulation is developed to simultaneously infer
physical quantities of interest, whilst considering the problem of identifying
the correct structural bending stiffness. The model is defined in section 2,
along with the optimization strategies, while section 3 is used to demonstrate
the inference capabilities and tolerances to measurement noise in a numerical
and controlled manner. Section 4 utilizes experimental measurements to update
an analytical model of a steel beam.

##  2 Physics-informed GP model of an Euler-Bernoulli beam

###  2.1 Model definition

Consider the general linear form of the Euler-Bernoulli beam equation:

$EI\frac{d^{4}u(\bm{x})}{d\bm{x}^{4}}=q(\bm{x}),$ (1)

where $u(\bm{x})$ are the beam deflections at positions $\bm{x}$ due to a
given input load $q(\bm{x})$ and $EI$ is the structural bending stiffness,
corresponding to Young’s modulus of elasticity $E$ and the second moment of
area $I$.

In the physics-informed model derived here, the deflections $u(\bm{x})$ are
represented as a zero-mean Gaussian process:

$u(\bm{x})\sim\mathcal{GP}(\bm{0},k_{uu}(\bm{x},\bm{x}^{\prime};\bm{\theta})),$
(2)

where $k_{uu}$ is a covariance kernel parametrized by the values in
$\bm{\theta}$. Several kernels are available in the literature, and without
loss of generality, the deflection model is here defined with the universal
squared exponential kernel

$k_{uu}(\bm{x},\bm{x}^{\prime};\bm{\theta}))=\sigma_{s}^{2}\mathrm{exp}\left(-\frac{1}{2}\left(\frac{\bm{x}-\bm{x}^{\prime}}{\ell}\right)\right),$
(3)

where $\theta=\\{\sigma_{s},\ell\\}$ contains the kernel standard deviation
$\sigma_{s}$ and the length scale $\ell$, which controls the covariance
smoothness [3].

By exploiting the linear aspect of GPs, a similar model can be derived for the
forces:

$q(\bm{x})\sim\mathcal{GP}(\bm{0},k_{qq}(\bm{x},\bm{x}^{\prime};\bm{\theta},EI)),$
(4)

where the kernel $k_{qq}$ is derived from $k_{u}u$ through the beam model
differential equation, such that

$k_{qq}(\bm{x},\bm{x}^{\prime};\bm{\theta},EI)=EI\frac{d^{4}}{d\bm{x}^{4}}\left(EI\frac{d^{4}}{d\bm{x}^{\prime
4}}k_{uu}(\bm{x},\bm{x}^{\prime};\bm{\theta})\right).$ (5)

The cross-covariance functions between the applied loads and the resulting
deflections can be derived as

$k_{uq}(\bm{x},\bm{x}^{\prime};\bm{\theta},EI)=EI\frac{d^{4}}{d\bm{x}^{\prime
4}}k_{uu}(\bm{x},\bm{x}^{\prime};\bm{\theta}).$ (6)

In a similar manner, the Gaussian process models for other related quantities
of interest provided by the model, such as rotations $r$, strains $\epsilon$,
bending moments $m$ and shear forces $v$ can be derived from $k_{uu}$ through
the application of the respective linear differential operator. The complete
Euler-Bernoulli beam model is then given by the multi-output GP

$\left[\bm{u},\bm{r},\bm{\epsilon},\bm{m},\bm{v},\bm{q}\right]^{T}=\mathcal{GP}(\bm{\mu}_{p},\bm{K}_{p}),$
(7)

with the prior mean
$\mu_{p}=\left[\bm{0},\bm{0},\bm{0},\bm{0},\bm{0},\bm{0}\right]^{T}$ and the
prior covariance matrix $\bm{K}_{p}$ formed as

$\bm{K}_{p}=\begin{bmatrix}\bm{K}_{ww}^{n}\\!&\bm{K}_{wr}\\!&\bm{K}_{w\epsilon}\\!&\bm{K}_{w\\!m}\\!&\bm{K}_{wv}\\!&\bm{K}_{wq}\\!\\\
\bm{K}_{rw}\\!&\bm{K}_{rr}^{n}\\!&\bm{K}_{r\epsilon}\\!&\bm{K}_{r\\!m}\\!&\bm{K}_{rv}\\!&\bm{K}_{rq}\\!\\\
\bm{K}_{\epsilon w}\\!&\bm{K}_{\epsilon
r}\\!&\bm{K}_{\epsilon\epsilon}^{n}\\!&\bm{K}_{\epsilon\\!m}\\!&\bm{K}_{\epsilon
v}\\!&\bm{K}_{\epsilon q}\\!\\\
\bm{K}_{m\\!w}\\!&\bm{K}_{m\\!r}\\!&\bm{K}_{m\\!\epsilon}\\!&\bm{K}_{m\\!m}^{n}\\!&\bm{K}_{m\\!v}\\!&\bm{K}_{m\\!q}\\!\\\
\bm{K}_{vw}\\!&\bm{K}_{vr}\\!&\bm{K}_{v\epsilon}\\!&\bm{K}_{v\\!m}\\!&\bm{K}_{vv}^{n}\\!&\bm{K}_{vq}\\!\\\
\bm{K}_{qw}\\!&\bm{K}_{qr}\\!&\bm{K}_{q\epsilon}\\!&\bm{K}_{q\\!m}\\!&\bm{K}_{qv}\\!&\bm{K}_{qq}^{n}\\!\end{bmatrix}.$
(8)

where $K_{ab}=k_{ab}(\bm{x},\bm{x}^{\prime};\bm{\theta},EI)$ is a covariance
matrix generated with its respective kernel function. To account for
measurement noise in each dataset, an additional diagonal block matrix
$k_{a}^{n}$ is added to the covariance of each kernel matrix, calculated as:

$\bm{K}_{a}^{n}=\sigma_{a}^{2}\delta_{ij}(\bm{x},\bm{x}^{\prime}),$ (9)

where $\delta_{ij}$ is the Kronecker delta operator and $\sigma_{a}$ is the
noise standard deviation value for dataset a, such that
$\bm{K}_{aa}^{n}=\bm{K}_{aa}+\bm{K}_{a}^{n}$. The standard deviation for each
dataset, when unknown, is included in a parameter vector $\bm{\psi}$ and
becomes an optimizable variable within the model.

Accounting for boundary conditions (BCs) in the model is possible through the
application of a modified Green’s function to the kernel equation [8, 9],
which restricts the generated model to a single structural system. Instead, in
this work BCs are accounted for through the creation of an artificial dataset,
at the appropriate locations, that measures the BC values in a noise-less
manner.

###  2.2 Model selection

The selection of an appropriate model in face of the measured noisy data
involves the optimization of the parameter vector
$\bm{\psi}={\bm{\theta},EI,\bm{\psi}}$, that is, the kernel values, the
bending stiffness and the noise in each dataset, respectively. For that
matter, a probability distribution $p(\bm{\psi})$ can be created to represent
the prior beliefs on each of the model parameters. Following the assumption of
independence in the parameters $\bm{\psi}$, the log form of this distribution
is calculated as:

$\mathrm{log}p(\bm{\psi})=\sum_{i}\mathrm{log}p(\psi_{i}).$ (10)

Moreover, when data is presented to the model in the form of
$\bm{y}=[\bm{u},\bm{r},\bm{\epsilon},\bm{m},\bm{v},\bm{q}]^{T}$, measured at
locations
$\bm{x}=[\bm{x}_{u},\bm{x}_{r},\bm{x}_{\epsilon},\bm{x}_{m},\bm{x}_{v},\bm{x}_{q}]^{T}$,
for $\bm{x},\bm{y}\in\mathcal{R}^{N\times 1}$, the log-likelihood can be
analytically calculated by:

$\mathrm{log}p(\bm{y}|\bm{x},\bm{\psi})=-\frac{1}{2}\bm{y}^{T}\bm{K}_{p}^{-1}\bm{y}-\frac{1}{2}\mathrm{log}|\bm{K}_{p}|-\frac{N}{2}\mathrm{log}2\pi,$
(11)

where $\bm{K}_{p}$ is the global covariance matrix as shown in equation 8,
calculated for all datasets contained in $\bm{x}$ and $\bm{y}$ [3].

Following Bayes’ rule, the log-prior and log-likelihood can be combined to
generate a posterior distribution that reflects the influence of measured data
in the model parameters, such that

$\mathrm{log}(\bm{\psi}|\bm{y},\bm{x})\propto\mathrm{log}\
p(\bm{y}|\bm{\psi},\bm{x})+\mathrm{log}\ p(\bm{\psi}).$ (12)

The formulation above is proportional to the true log-posterior up to a value
defined by the marginal likelihood of the system, which is constant with
respect to the parameters $\bm{\psi}$. The marginal likelihood takes the form
of an intractable integral for most probabilistic systems, and therefore the
parameter selection is carried out as an optimization problem given by:

$\bm{\psi}=\mathrm{argmax}_{\bm{\psi}}\ \mathrm{log}\
p(\bm{\psi}|\bm{y},\bm{x}),$ (13)

which is popularly known as maximum a posteriori estimation. The maximization
of this posterior distribution is estimated numerically in a probabilistic
approach through the Metropolis-Hastings (MH) algorithm [10]. The algorithm is
a Markov chain Monte Carlo method that returns a sequence of auto-correlated
random samples of the parameters
$\bm{\Psi}=[\bm{\psi}_{0},\bm{\psi}_{1},...,\bm{\psi}_{k}]^{T}$, which are
representative of the true posterior distribution
$p(\bm{\psi}|\bm{y},\bm{x})$.

###  2.3 Model inference

Predictions of quantities $\bm{y}_{\star}$ at unobserved locations
$\bm{x}_{\star}$ can be made by conditioning the predictive distribution on
the noisy observations used during training, and integrating over the
identified parameter distribution,

$p(\bm{y}_{\star}|\bm{x}_{\star},\bm{y,\bm{x}})=\int
p(\bm{y}_{\star}|\bm{x}_{\star},\bm{y},\bm{x},\bm{\psi})p(\bm{\psi}|\bm{y},\bm{x})d\bm{\psi}.$
(14)

The inner probability distribution in the above formulation corresponds to the
standard GP predictive posterior when no uncertainty in the parameters
$\bm{\psi}$ exists. Due to the Gaussianity assumptions of the GP model, this
probability takes the closed form of
$p(\bm{y}_{\star}|\bm{x}_{\star},\bm{y},\bm{x},\bm{\psi})=\mathcal{N}(\bm{\mu}_{\bm{y}^{\star}},\bm{K}_{\bm{y}^{\star}})$,
with the mean $\bm{\mu}_{\bm{y}^{\star}}$ and the covariance matrix
$\bm{K}_{\bm{y}^{\star}}$ calculated by, respectively,

$\bm{\mu}_{\bm{y}^{\star}}=\bm{K}_{\star}^{T}\bm{K}_{p}^{-1}\bm{y},$ (15)

$\bm{K}_{\bm{y}^{\star}}=\bm{K}_{\star\star}-\bm{K}_{\star}^{T}\bm{K}_{p}^{-1}\bm{K}_{\star},$
(16)

where
$\bm{K}_{\star}=[k_{au}(\bm{x}_{\star},\bm{x}^{\prime}),...,k_{aq}(\bm{x}_{\star},\bm{x}^{\prime})]^{T}$
is the cross-covariance matrix between unobserved locations and training
positions for all training datasets, and
$\bm{K}_{(}{\star}{\star})=k_{aa}(\bm{x}_{\star},\bm{x}^{\prime}_{\star})$ is
the self-covariance matrix, for the particular quantity of interest $a$ [3,
11]. A closed-form solution of Equation 14 is generally intractable, and
therefore the mean and covariance parameters of the predictive distribution
are estimated as:

$p(\bm{y}_{\star}|\bm{x}_{\star},\bm{y},\bm{x})\approx\frac{1}{N}\sum_{i}^{N}p(\bm{y}_{\star}|\bm{x}_{\star},\bm{y},\bm{x},\bm{\psi}_{i}),$
(17)

where $\bm{\psi}_{i}\sim p(\bm{\psi}|\bm{y},\bm{x})$ are draws from the
parameter posterior approximated with the MH algorithm [11]. The predictive
posterior finally takes the form of a multivariate mixture of Gaussians,
following the assumption of Gaussian noise in the dataset $\bm{y}$.

##  3 Numerical experiments

In this section, a numerical investigation of the presented method is carried
out for a cantilever beam of length L with constant bending stiffness EI,
subjected to a uniformly distributed load of magnitude q, as shown in Figure
1. For such a case an analytical solution exists, such that

$u(x)=\frac{qx^{2}}{24EI}x^{2}-4Lx+6L^{2}.$ (18)

This solution is further used as a benchmark for the model’s predictions,
along with the true bending stiffness EI, which is compared to the value
regressed by the GP model. The software Matlab is used for the numerical
implementation of the method.

Figure 1: Benchmark structure for a numerical case: a cantilever beam with
constant stiffness $EI$ and length $L$, subjected to a uniformly distributed
load of magnitude $q$.

###  3.1 Stiffness regression and latent function inference

In this work, prior knowledge on the noise values and the kernel parameters is
modelled with a uniform distribution
$p(\bm{\psi}_{i})=\mathcal{U}(-\infty,\infty)$, while the bending stiffness is
treated in constrained manner as
$p(\psi_{E}I)=\mathcal{U}(0.1EI_{\mathrm{true}},2.0EI_{\mathrm{true}})$. It is
further assumed that four deflection sensors, equally spaced throughout the
length of the structure, monitor the system. These sensors normally operate
under the influence of environmental conditions and, due to the nature of the
measurement, present a variable level of uncertainty in their outputs. Noise
values are herein simulated assuming a Gaussian model
$u_{n}~{}\mathcal{N}(u_{\mathrm{ana}}(x),\sigma_{n}^{u})$, where
$u_{\mathrm{ana}}(x)$ is the analytical solution and $\sigma_{n}^{u}$ is a
standard deviation given by

$\sigma_{u}^{n}=\frac{|\hat{u}|}{\mathrm{SNR}}$ (19)

where $|\hat{u}|$ is the maximum absolute displacement at the tip of the
cantilever beam, and SNR is a signal-to-noise ratio parameter.

By providing the load value q, along with a set of sensor readings obtained
from different sensor locations along the structure, and contaminating the
measurements with a noise defined by a SNR=10, the GP model can be trained and
used for further inference. A sample of the predicted displacement field is
shown in Figure 2. Despite the noisy measurements, the GP model is able to
accurately predict the displacement values along the length of the beam. In
addition, the model uncertainty, represented by the prediction’s standard
deviation $\sigma_{u}$, increases along the length of the structure. This
effect reflects the boundary condition assumption, introduced as an additional
artificial noise-less sensor located at $x/L=0$.

Figure 2: Normalized mean $\mu_{u}$ and standard deviation $\sigma_{u}$ of the
displacement field predictions based on 5 data points collected by 4 noisy
sensors along the structure.

Once the model is optimized, a set of auto-correlated parameters is returned
from the MH algorithm, which approximates the posterior parameter distribution
$p(\bm{\psi}|\bm{y},\bm{x})$. The results, shown in Figure 3, reflect the GP
model uncertainty on each of the parameters in question.

Figure 3: Correlation matrix of the probabilistic model parameters optimized
using the MH algorithm. The kernel values $\sigma_{s}$ and $\ell$ show a
positive degree of correlation, whereas the remaining parameters are
uncorrelated.

Of particular interest in the parameter model is the stiffness distribution,
shown in Figure 4. The mean value of $p(\psi_{E}I)$ approximates the true
stiffness $EI_{\mathrm{true}}$ with a 0.63% error. The standard deviation,
calculated as $0.023EI_{\mathrm{true}}$, reflects the model uncertainty on the
stiffness and is a function of the training dataset. Further, a normalized
version of the stiffness distribution is defined
as$p(\psi_{EI}/EI_{\mathrm{true}})=\mathcal{N}(\mu_{EI},\sigma_{EI})$, and the
Mahalanobis distance, given by

$d_{M}=\sqrt{\frac{(\mu_{EI}-1)^{2}}{\sigma_{EI}^{2}}},$ (20)

is used to evaluate the quality of the predicted stiffness distribution. This
formulation doesn’t only reflect the accuracy of the mean value but also
penalizes uncertain stiffness models. The Mahalanobis distance of the
stiffness distribution in Figure 4 was calculated as $d_{M}=0.40$.

Figure 4: The regressed probabilistic stiffness value of the structure. The
distribution’s mean value approaches the correct stiffness
$EI_{\mathrm{true}}$, while the standard deviation represents the model
uncertainty.

In addition to the regressed parameters, the fully trained GP model can also
be used to infer physical quantities that were not directly measured, as shown
in Figure 5.

Figure 5: Mean $\mu$ and standard deviation $\sigma$ of the inferences on the
latent rotations $r$, strains $\epsilon$, moments $m$ and shear forces $v$,
normalized by the maximum absolute analytical result of the respective
physical quantity.

Similarly to the displacement case, the mean values predicted for rotations,
strains, bending moments and shear forces are in accordance with the
analytical results, with a root mean squared error (RMSE) calculated for each
of the normalized predictions in the order of $10^{(}-3)$. The standard
deviations of the rotation predictions follow the trend of the displacement
values, increasing along the length of the structure, as a result of the
boundary condition enforced as $r(x=0)=0$ rad. Similarly, the strains, bending
moments and shear values show an increase in uncertainty closer to the support
location, once again, due to the boundary conditions $m(x=L)=0$ Nm and
$v(x=L)=0$ N.

###  3.2 Influence of noise and number of data points

The training of this physics-informed GP model takes place iteratively, by
consecutive evaluations of the probabilistic a posteriori distribution in the
Markov chain of the Metropolis-Hastings algorithm. A drawback of a standard
Gaussian process model is that, for every model evaluation, a computational
cost of $\mathcal{O}(N^{3})$ is involved.

Considering typical sampling rates for displacement sensors, a large amount of
data can be quickly generated, and the usage of the GP model becomes
unfeasible. Therefore, it is important to reduce the size of the dataset
provided to the GP algorithm during training, while retaining a good
prediction accuracy. An obvious correlation exists between the number of
provided data points $N_{dp}$ and the level of noise of the measurements,
which is numerically controlled by the signal-to-noise ratio (SNR). In Figure
6 the interaction between noise and data size is observed. Considering four
equidistant sensors (cf. Figure 2), the GP model has a clear improvement in
accuracy when higher quality data, and more data points, are used as inputs.
The trade-off exists as a model response to correctly identify the mean and
standard deviation values, given the noisy points provided.

Figure 6: Influence of measurement noise and the number of data points
provided to the algorithm. A total of 4 equidistant sensors were used in this
study, as shown in Figure 2

###  3.3 Damage identification capabilities

The parameter study presented in section 3.2 allows for an informed decision
on sensor quality, as well as the computational costs of the GP model, given a
minimum desired stiffness prediction accuracy. If a good prediction is
guaranteed, and continuous monitoring of a structure is available, then the
generated model can be used to detect, locate and estimate the severity of
eventual damages. Damage is estimated through the deviation of the original
numerical stiffness $EI_{\mathrm{true}}$ using the Mahalanobis distance,
calculated as per Equation 20. To simulate damage, a finite element model of
the cantilever beam is generated, with a total of 20 elements. The damage is
simulated as a reduction in the bending stiffness of up to 40% of the original
structure, in one element at a time, for all the finite elements in the model.
In Figure 7 the deviations of the original structural stiffness are shown, in
terms of $d_{M}$, as a function of damage location and extent.

Figure 7: Damage study: influence of stiffness reduction and corresponding
location along the structure on the model’s predictions.

The results indicate an increase in the Mahalanobis distance value for
progressive reductions of remaining stiffness and damage locations closer to
the support. Given that the structural system is a cantilever, the system’s
displacement response is more sensitive to changes in close to the support.
For damage cases at locations $x>0.60L$, the indication of damage may be
obscured by the model’s prediction error, as a function of the noise levels.
Nevertheless, given a damage condition closer to the support, the Mahalanobis
distance allows for an informative decision of the location of the damage.
Once the location is known, the damage extent identification can be achieved,
given the linearity condition of the system.

##  4 Experimental model updating

To validate the proposed model, an experimental study conducted on a simply
supported steel beam is now presented. The true bending stiffness is initially
unknown but estimated through material and geometrical parameters to be
$EI_{\mathrm{init}}$. The goal is to use measurement data and provide an
optimized version of the bending stiffness value, in a model updating sense.
Figure 8 displays the structure and measurement setup. Different deflection
sensors, such as laser deflectometers, displacement transducers, draw wire
transducers and dial gauges were used to monitor the deflections due to the
applied load.

Figure 8: Experimental set-up: deflection, inclination and strain sensors are
installed in a simply-supported steel beam. A static load uniformly
distributed is applied with bags of steel spheres, resting on wooden plates.

The degree of noise that contaminates the measurements varies with each
measurement device and is estimated by the GP through the standard deviation
$\sigma_{a}$ in equation 9, for each dataset a. In that manner, the developed
framework is able to differentiate high and low-quality measurements and
prioritize high-fidelity datasets. In addition, inclinometers were used at
both support locations to measure the structural rotation, and a strain gauge
was installed at mid-span, on the bottom side of the beam. The SNRs calculated
for the datasets are all higher than 15, and based on Figure 6, the number of
data points provided to the GP algorithm can be defined. In this model, a
total of $N_{dp}=7$ is used for all of the available measurement sets, as
shown in Figure 9. The results indicate that the original calculated stiffness
value $EI_{\mathrm{init}}$ is smaller than the real structural stiffness, as
the structural responses, in the form of the measurements, tend to be smaller
than predicted. The GP framework is able to correctly identify a model that
explains the heterogeneous datasets composed of deflections, rotations and
strains. In addition, it also ignores the deflection information at
approximately $x=0.6L$, as it appears to originate from a malfunctioning
sensor. The uncertainty levels, similar to the numerical example, reflect the
boundary conditions informed by the GP model in the form of noise-less
datasets.

Figure 9: Normalized inference on the trained model using the experimental
results, as compared to the initial, analytical model using
$EI_{\mathrm{init}}$, shown in blue.

To quantify the updated stiffness value, the probabilistic model of
$\psi_{EI}$ is shown in Figure 10. The normalized $EI_{\mathrm{init}}$ is used
as an initial point for the optimization algorithm, and a stable distribution
is obtained with a mean $\mu_{EI}$ that is 19.08% larger than initially
assumed. The model uncertainty, measured in the form of the standard deviation
of $\psi_{EI}$, amounts to $\sigma_{EI}=0.11$%.

Figure 10: Updated probabilistic bending stiffness model and the initial
assumption $EI_{\mathrm{init}}$.

##  5 Conclusions

A physics-informed Gaussian process model has been developed on the basis of
the Euler-Bernoulli beam theory. The model is of a hybrid nature, being driven
simultaneously by the heterogeneous data provided to it, as well as the
mathematical theory in the form of differential equations. The probabilistic
model selection involves the optimization of the GP parameters, the
identification of noise in the datasets and, most importantly, the regression
of the bending stiffness value in the form of a probability distribution. The
identified model can further be used to make probabilistic inferences in any
physical quantity linked through the differential equation, even when no data
from them is directly provided. The identified stiffness distribution is
directly related to the noise levels of the data provided to the model. A
highly noisy dataset requires, therefore, more information to provide a
reasonable stiffness prediction. In consequence, a trade-off exists between
the data quality and the number of data points provided to the GP model, which
leads to higher computation costs. The regressed stiffness parameter, once
correctly identified, can be used in a damage identification scenario to
reason about the location and the severity of such a damage case. Finally, an
experimental case has been used to showcase that the developed physics-
informed GP in a model updating framework. For that matter, the model is also
able to integrate heterogeneous datasets, with different measured physical
quantities, and different sets of the same quantity with varying levels of
noise.

## References

  * [1] Keith Worden and Graeme Manson.  The application of machine learning to structural health monitoring.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 365(1851):515–537, 2007.
  * [2] George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang.  Physics-informed machine learning.  Nature Reviews Physics, 3(6):422–440, 2021.
  * [3] Christopher KI Williams and Carl Edward Rasmussen.  Gaussian processes for machine learning, volume 2.  MIT press Cambridge, MA, 2006.
  * [4] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis.  Inferring solutions of differential equations using noisy multi-fidelity data.  Journal of Computational Physics, 335:736–746, 2017.
  * [5] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis.  Numerical gaussian processes for time-dependent and nonlinear partial differential equations.  SIAM Journal on Scientific Computing, 40(1):A172–A198, 2018.
  * [6] Alastair Gregory, F Din-Houn Lau, Mark Girolami, Liam J Butler, and Mohammed ZEB Elshafie.  The synthesis of data from instrumented structures and physics-based models via gaussian processes.  Journal of Computational Physics, 392:248–265, 2019.
  * [7] Maziar Raissi and George Em Karniadakis.  Hidden physics models: Machine learning of nonlinear partial differential equations.  Journal of Computational Physics, 357:125–141, 2018.
  * [8] Lars Lau Raket.  Differential equations, splines and gaussian processes.  arXiv preprint arXiv:2102.03306, 2021.
  * [9] Simo Särkkä.  Linear operators and stochastic partial differential equations in gaussian process regression.  In Artificial Neural Networks and Machine Learning–ICANN 2011: 21st International Conference on Artificial Neural Networks, Espoo, Finland, June 14-17, 2011, Proceedings, Part II 21, pages 151–158. Springer, 2011.
  * [10] W Keith Hastings.  Monte carlo sampling methods using markov chains and their applications.  1970\.
  * [11] Vidhi Lalchand and Carl Edward Rasmussen.  Approximate inference for fully bayesian gaussian process regression.  In Symposium on Advances in Approximate Bayesian Inference, pages 1–12. PMLR, 2020.