© 2022 IEEE. Personal use of this material is permitted. Permission from IEEE
must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes,
creating new collective works, for resale or redistribution to servers or
lists, or reuse of any copyrighted component of this work in other works.

# Optimal Friendly Jamming and Transmit Power Allocation in RIS-assisted
Secure Communication

Burhan Wafai1, Chinmoy Kundu2, Ankit Dubey3, and Mark F. Flanagan4,
{1burhan.wafai, 3ankit.dubey}<EMAIL_ADDRESS><EMAIL_ADDRESS>13Department of EE, Indian Institute of Technology
Jammu, Jammu & Kashmir, India  24School of Electrical and Electronic
Engineering, University College Dublin, Belfield, Ireland

###### Abstract

This paper analyzes the secrecy performance of a reconfigurable intelligent
surface (RIS) assisted wireless communication system with a friendly jammer in
the presence of an eavesdropper. The friendly jammer enhances the secrecy by
introducing artificial noise towards the eavesdropper without degrading the
reception at the destination. Approximate secrecy outage probability (SOP) is
derived in closed form. We also provide a simpler approximate closed-form
expression for the SOP in order to understand the effect of system parameters
on the performance and to find the optimal power allocation for the
transmitter and jammer. The optimal transmit and jamming power allocation
factor is derived by minimizing the SOP assuming a total power constraint. It
is shown that the SOP performance is significantly improved by the
introduction of the jammer and a gain of approximately $3$ dB is achieved at
an SOP of $10^{-4}$ by optimally allocating power compared to the case of
equal power allocation.

###### Index Terms:

Optimal power allocation, physical layer security, reconfigurable intelligent
surface, secrecy outage probability.

##  I Introduction

To address the demands of next-generation communication networks,
reconfigurable intelligent surface (RIS) technology is a promising candidate
for enabling smart control of the propagation of signals. An RIS consists of a
large number of passive meta-material elements, and can intelligently control
the phase shift and amplitude at each element and is also energy efficient due
to its passive elements [1, 2].

Despite RIS being a relatively new technology, its integration with other
wireless technologies has been studied extensively. Enhanced energy detection
using RIS was proposed in [3] for single-user spectrum sensing, cooperative
spectrum sensing; the performance of diversity reception and average
probability of detection was studied in each case. The authors in [4] studied
the ergodic capacity of user scheduling in RIS-assisted multi-user
communication. In [5], the authors analyzed the performance of an RIS-aided
communication system in terms of the outage probability, average symbol error
probability, and achievable rate bounds. Although the performance analysis in
an RIS-assisted communication is studied in the aforementioned papers, the
secrecy aspect has not been considered.

The broadcast nature of wireless transmission is always vulnerable to
eavesdropping. Hence the security of such networks is of utmost importance,
and consequently, physical layer security (PLS) has emerged as a viable
solution [6]. The authors in [7, 8, 9] have studied various optimization
techniques to maximize the secrecy rate of RIS enabled wireless systems and
designed the jointly optimal transmit beamformer and RIS phase shift matrix.
The above papers consider performance optimization rather than secrecy
performance analysis. The secrecy performance of an RIS-aided indoor
communication system via secrecy outage probability (SOP) and average secrecy
rate is considered in [10] assuming Rician fading model. The SOP is also
derived assuming a Rayleigh fading model in [11]. These papers do not consider
friendly jamming to enhance the secrecy performance or optimize the SOP
performance using power allocation.

Artificial noise (AN) can be added to the received signal by the transmitter
or a jammer to degrade the signal-to-noise ratio (SNR) at the eavesdropper
without affecting the SNR at the destination [12]. The legitimate user and the
transmitter/jammer cooperate in the system, and thus the user can cancel out
the effect of jamming. Furthermore, cooperative jamming (CJ), where a relay
node transmits a jamming signal to create interference at the eavesdropper, is
another way to improve the secrecy performance [13]. The authors in [14]
derived the SOP of a joint relay and friendly jammer selection scheme from a
set of multiple energy harvesting relays. In [15], the authors considered
optimal power allocation to cooperative jammers to maximize the secrecy rate
and obtained the corresponding SOP performance .

Though using AN and CJ to improve secrecy is common in the secrecy literature,
the application of AN or CJ to RIS-aided system is yet in the nascent stage.
The authors in [16] first investigated the joint transmit beamforming with AN
and RIS reflect beamforming in an RIS-aided secure communication system to
maximize the achievable secrecy rate. The authors in [17] obtained the optimal
secure beamforming and phase shift matrix for the RIS in a multiple-input
single-output (MISO) system with a friendly jammer by maximizing energy
efficiency. In [18], the authors deployed an RIS not to aid information
transmission, however but as a jammer to redirect the energy of the
confidential signal and deteriorate the reception at the eavesdropper. In this
context, we also mention that authors in [19] used jamming to aid
eavesdropping and cooperation between the jammer and eavesdropper was
considered.

The literature mentioned above with friendly jammers considered the joint
design of transmit beamforming and RIS phase shift matrix, to maximize the
secrecy rate or energy efficiency. The SOP performance was however not
evaluated. To the best of our knowledge, the SOP analysis using friendly
jamming and optimal power allocation to minimize the SOP has not been
considered yet in the literature. Motivated by this, we consider an RIS-aided
system in which the signals from the jammer and the source are impinging on
the RIS before being reflected to the destination, and obtain the approximate
closed-form SOP performance. Further, power is optimally allocated among the
source and the jammer in order to achieve minimum SOP assuming a total power
constraint.

###  I-A Key Contributions

The main contributions of this work are listed as follows:

  1. 1.

We find the approximate closed-form SOP of an RIS enabled system with a single
source, destination, eavesdropper, and a friendly jammer transmitting noise to
degrade eavesdropping.

  2. 2.

A simpler approximate closed-form SOP expression is further derived to
implement the optimal power allocation at the transmitter and the jammer, and
to provide a better insight regarding the performance depends on the system
parameters.

  3. 3.

Optimal power allocation for the transmitter and the jammer is carried out by
proving the convexity of the approximate SOP expression assuming a total power
constraint.

  4. 4.

The effect of the various system parameters on the SOP and the optimal power
allocation factor is shown.

Notation: For a random variable (r.v.) $X$, $\mathbb{E}_{X}[\cdot]$,
$F_{X}(\cdot)$, and $f_{X}(\cdot)$ denote the expectation, cumulative
distribution function (CDF), and probability density function (PDF),
respectively. $\mathbb{P}[\cdot]$ represents the probability of occurrence of
an event.

##  II System and channel model

We consider a wireless communication system wherein a source $S$ intends to
communicate with a legitimate destination $D$ with the aid of an RIS $R$. A
passive eavesdropper $E$ attempts to overhear the communication between $S$
and $D$. The RIS with $N$ elements enables the communication by intelligently
reflecting the signals from $S$ towards $D$. Since $E$ can overhear the
communication via RIS, a friendly jammer $J$ is introduced to transmit a
jamming signal via $R$ independent of the source signal, to degrade the signal
quality at $E$. A scenario is considered where an obstacle is present and the
blockage is so severe that the signal from $S$ and $J$ cannot reach $D$ and
$E$ directly. The jamming signal is assumed to be known to $D$ and hence, it
causes controlled interference at $E$ but not at $D$, thereby enhancing the
system’s secrecy performance. Nodes $S$, $J$, $D$, and $E$, are assumed to be
equipped with a single antenna. The signals from $S$ and $J$ are received by
$D$ and $E$ after being reflected from each element of the RIS having
amplitude reflection coefficient $\eta_{n}$ and adjustable phase
$\theta_{n}\in(0,2\pi]$, where $n\in\\{1,\ldots,N\\}$. Transmitting nodes, $S$
and $J$, have a total power budget of $P_{T}$, of which a fraction $\alpha$ is
used by $S$ and a fraction $(1-\alpha)$ fraction by $J$, where $0<\alpha<1$ is
the power allocation factor for $S$. The channel for node pair combinations
$XY\in\\{SR,RD,JR,RE\\}$ is denoted as
$h_{XY}^{(n)}=\sqrt{\zeta_{XY}}\tilde{h}_{XY}^{(n)}$, wherein
$\tilde{h}_{XY}^{(n)}$ is a complex Gaussian r.v. with zero mean and unit
variance and the path loss is
$\zeta_{XY}=\big{(}\frac{d_{0}}{d_{XY}}\big{)}^{\upsilon}$ where $d_{0}$ is
reference distance, $d_{XY}$ denotes the distance between the nodes, and the
path loss exponent $\upsilon$.

The received signal at $D$ through the RIS is given by

$\displaystyle y_{D}$ $\displaystyle=\sqrt{\alpha
P_{T}}\sum_{n=1}^{N}h_{RD}^{(n)}\eta_{n}e^{j\theta_{n}}h_{SR}^{(n)}s$
$\displaystyle+\sqrt{(1-\alpha)P_{T}}\sum_{n=1}^{N}h_{RD}^{(n)}\eta_{n}e^{j\theta_{n}}h_{JR}^{(n)}z+w_{D},$
(1)

where $s$ and $z$ denote the information and jamming signals transmitted by
$S$ and $J$, respectively, with
$\mathbb{E}\\{|s|^{2}\\}=\mathbb{E}\\{|z|^{2}\\}=1$, $w_{D}$ denotes the
complex additive white Gaussian noise (AWGN) at $D$ with zero-mean and
variance $N_{0}$. The received signal at $E$ through the RIS is given by an
expression similar to (II) but with $D$ replaced by $E$. Without loss of
generality, the noise at $E$ is also modelled as complex AWGN with the same
mean and a variance of $w_{D}$.

###  II-A SINR distribution of the legitimate link

Assuming RIS phase shift $\theta_{n}$ for each $n\in\\{1,\ldots,N\\}$ is
adjusted to maximize the SNR at $D$ by constructive signal combining with the
knowledge of the phases of $h_{RD}^{(n)}$ and $h_{SR}^{(n)}$; and $D$ has all
the knowledge to cancel the interference from $J$, the SINR at $D$ is written
as

$\displaystyle\Gamma_{D}$
$\displaystyle=\alpha\Gamma_{0}\Bigg{(}\sum_{n=1}^{N}\left|h_{RD}^{(n)}\right|\left|h_{SR}^{(n)}\right|\Bigg{)}^{2},$
(2)

where $\Gamma_{0}={P_{T}}/{N_{0}}$ and $\eta_{n}=1$ for each $n$ is assumed
throughout the paper. Since an RIS typically has large $N$, from the central
limit theorem, the distribution of the sum of r.v.s within the parenthesis in
(2) follows a Gaussian distribution, the mean $\mu$ and variance $\sigma^{2}$
of which is computed as $\mu=\frac{\pi N\sqrt{\zeta_{RD}\zeta_{SR}}}{4}$ and
$\sigma^{2}=N\zeta_{RD}\zeta_{SR}\left(\frac{16-\pi^{2}}{16}\right)$ [20].
Therefore, the CDF of $\Gamma_{D}$ is written as [5]

$\displaystyle F_{\Gamma_{D}}(x)$
$\displaystyle=1-\psi{Q}\Bigg{(}\frac{\sqrt{\frac{x}{\alpha\Gamma_{0}}}-\mu}{\sigma}\Bigg{)}\text{for}~{}x\geq
0,$ (3)

where ${Q}(\cdot)$ denote the Gaussian ${Q}$-function and
$\psi=1/Q\left(-{\mu/\sigma}\right)$ is the normalization factor.

###  II-B SINR distribution of the eavesdropping link

As $E$ cannot cancel out the jamming signal, the SINR at $E$ is expressed as

$\displaystyle\Gamma_{E}$ $\displaystyle=\frac{\left|\sqrt{\alpha
P_{T}}\sum_{n=1}^{N}h_{RE}^{(n)}\eta_{n}e^{j\theta_{n}}h_{SR}^{(n)}\right|^{2}}{N_{0}+\left|\sqrt{(1-\alpha)P_{T}}\sum_{n=1}^{N}h_{RE}^{(n)}\eta_{n}e^{j\theta_{n}}h_{JR}^{(n)}\right|^{2}}$
$\displaystyle=\frac{\left|H_{SE}\right|^{2}}{1+\left|H_{JE}\right|^{2}},$ (4)

where
$H_{SE}=\sqrt{\alpha\Gamma_{0}}\sum_{n=1}^{N}h_{RE}^{(n)}\eta_{n}e^{j\theta_{n}}h_{SR}^{(n)}$
and
$H_{JE}=\sqrt{(1-\alpha)\Gamma_{0}}\sum_{n=1}^{N}h_{RE}^{(n)}\eta_{n}e^{j\theta_{n}}h_{JR}^{(n)}$.

As the RIS phase shift $\theta_{n}$ is adjusted to benefit $D$, the phases of
channels $h_{SR}^{(n)}$ and $h_{RE}^{(n)}$ for each $n\in\\{0,\ldots,N\\}$
will be misaligned. Similarly, the phases of $h_{JR}^{(n)}$ and $h_{RE}^{(n)}$
for each $n\in\\{0,\ldots,N\\}$ will also be misaligned. Consequently,
$H_{SE}$ and $H_{JE}$ become complex Gaussian r.v.s, and $|H_{SE}|^{2}$ and
$|H_{JE}|^{2}$ become exponential r.v.s with parameters
${\lambda_{SE}}=\alpha\Gamma_{0}N\eta^{2}\zeta_{RE}\zeta_{SR}$ and
$\lambda_{JE}=(1-\alpha)\Gamma_{0}N\eta^{2}\zeta_{RE}\zeta_{JR}$ ,
respectively, following [11]. The PDF of $|H_{XY}|^{2}$ is thus expressed as

$\displaystyle f_{|H_{XY}|^{2}}(x)$
$\displaystyle=\frac{1}{\lambda_{XY}}\exp{\left(-\frac{x}{\lambda_{XY}}\right)}.$
(5)

Therefore, the distribution of $\Gamma_{E}$ is derived with the help of (5)
and its CDF as

$\displaystyle
F_{\Gamma_{E}}(x)=\int_{0}^{\infty}F_{|H_{SE}|^{2}}\left(x{\left(1+y\right)}\right)f_{|H_{JE}|^{2}}(y)dy$
$\displaystyle=\frac{1}{\lambda_{JE}}\int_{0}^{\infty}\left(1-\exp{\left(-\frac{(1+y)x}{\lambda_{SE}}\right)}\right)\exp{\left(-\frac{y}{\lambda_{JE}}\right)}dy$
$\displaystyle=1-\frac{\frac{\lambda_{SE}}{\lambda_{JE}}\exp{\left(-\frac{x}{\lambda_{SE}}\right)}}{\left({x+\frac{\lambda_{SE}}{\lambda_{JE}}}\right)}.$
(6)

By differentiating (II-B), the PDF of $\Gamma_{E}$ is expressed as

$\displaystyle
f_{\Gamma_{E}}(x)=\frac{\exp{\left(-\frac{x}{\lambda_{SE}}\right)}}{\lambda_{JE}\left({x+\frac{\lambda_{SE}}{\lambda_{JE}}}\right)}+\frac{{\lambda_{SE}}\exp{\left(-\frac{x}{\lambda_{SE}}\right)}}{{\lambda_{JE}}\left({x+\frac{\lambda_{SE}}{\lambda_{JE}}}\right)^{2}}.$
(7)

##  III Secrecy Outage Probability

This section evaluates the SOP of the system. The instantaneous secrecy rate
$C_{S}$ is defined as the difference between the legitimate channel capacity
and the wiretap channel capacity and is expressed as [21]

$\displaystyle
C_{S}=\max\left\\{\log_{2}\left(\frac{{1+{\Gamma}_{D}}}{1+{\Gamma}_{E}}\right),0\right\\}.$
(8)

The SOP is then defined as the probability that the secrecy rate $C_{S}$ falls
below a threshold rate $R_{th}$ and is expressed as [22]

$\displaystyle\mathcal{P}_{out}$
$\displaystyle=\mathbb{P}\left[C_{S}<R_{th}\right]$
$\displaystyle=\int_{0}^{\infty}F_{{\Gamma}_{D}}(\rho(x+1)-1)f_{{\Gamma}_{E}}(x)dx,$
(9)

where $\rho=2^{R_{th}}$. Using $F_{{\Gamma}_{D}}(x)$ and $f_{{\Gamma}_{E}}(x)$
derived in (3) and (7), respectively, in the previous section, the SOP is
expressed as

$\displaystyle\mathcal{P}_{out}$
$\displaystyle=1-\int_{0}^{\infty}\psi{Q}\left(\frac{\sqrt{\frac{\rho(x+1)-1}{\alpha\Gamma_{0}}}-\mu}{\sigma}\right)f_{{\Gamma}_{E}}(x)dx.$
(10)

To obtain the SOP in (10), we have to average the $Q$-function which is not
trivial. To obtain the integral solution in closed form, we use the following
approximation of $Q$-function [23]:

$\displaystyle
1-\sum_{i=1}^{3}\frac{w_{i}}{2}\exp\left(-\frac{p_{i}x^{2}}{2}\right)\textrm{when}~{}~{}x<0,$
(11a)
$\displaystyle\sum_{i=1}^{3}\frac{w_{i}}{2}\exp\left(-\frac{p_{i}x^{2}}{2}\right)~{}~{}~{}~{}\textrm{when}~{}~{}x\geq
0,$ (11b)

where
$w_{i}=\left\\{\frac{1}{6},\frac{1}{3},\frac{1}{3}\right\\}~{}\textrm{and}~{}p_{i}=\left\\{1,4,\frac{4}{3}\right\\}$
for $i\in\\{1,2,3\\}$. The integral in (10) is derived based on the positive
or negative argument of the $Q$-function following (11a) and (11b) as

$\displaystyle\mathcal{P}_{out}=1-\left(\mathcal{P}_{out}^{-}+\mathcal{P}_{out}^{+}\right),$
(12)

where

$\displaystyle\mathcal{P}_{out}^{-}$
$\displaystyle=\psi\int_{0}^{a}\left(1-\sum_{i=1}^{3}\frac{w_{i}}{2}\chi(x)\right)f_{{\Gamma}_{E}}(x)dx,$
(13) $\displaystyle\mathcal{P}_{out}^{+}$
$\displaystyle=\sum_{i=1}^{3}\frac{\psi
w_{i}}{2}\int_{a}^{\infty}\chi(x)f_{{\Gamma}_{E}}(x)dx,$ (14)

and
$\chi(x)=\exp{\Big{(}-\frac{p_{i}}{2\sigma^{2}}\big{(}\sqrt{\frac{\rho(x+1)-1}{\alpha\Gamma_{0}}}-\mu\big{)}^{2}\Big{)}}$.
The value $a$ which serves as the upper and lower limit of the integrals in
(13) and (14), respectively, is evaluated by finding limits of $x$ for which
$\frac{1}{\sigma}\Big{(}\sqrt{\frac{\rho(x+1)-1}{\alpha\Gamma_{0}}}-\mu\Big{)}$
is greater or less than zero. Thus, we obtain
$a=\frac{\alpha\mu^{2}\Gamma_{0}+1}{\rho}-1$.

Here we note that $\mathcal{P}_{out}^{+}$ tends to zero as $a$ becomes large
in (14) when $\Gamma_{0}$ or $N$ or both are large. Thus, the effect of
$\mathcal{P}_{out}^{+}$ is negligible as compared to $\mathcal{P}_{out}^{-}$
and the SOP approximates to $\mathcal{P}_{out}\approx
1-\mathcal{P}_{out}^{-}$. Then $\mathcal{P}_{out}^{-}$ is evaluated as

$\displaystyle\mathcal{P}_{out}^{-}=I_{0}-I_{1}-I_{2},$ (15)

where

$\displaystyle I_{0}$ $\displaystyle=\psi\int_{0}^{a}f_{\Gamma_{E}}(x)dx$
$\displaystyle=\psi\Bigg{[}\frac{e^{\frac{1}{\lambda_{JE}}}}{\lambda_{JE}}\left(\text{Ei}\left(-\frac{a}{\lambda_{SE}}-\frac{1}{\lambda_{JE}}\right)-\text{Ei}\left(-\frac{1}{\lambda_{JE}}\right)\right.\Bigg{.}$
$\displaystyle\Bigg{.}\left.+\Gamma\left(-1,\frac{1}{\lambda_{JE}}\right)-\Gamma\left(-1,\frac{a}{\lambda_{SE}}+\frac{1}{\lambda_{JE}}\right)\right)\Bigg{]},$
(16) $\displaystyle{I}_{1}$ $\displaystyle=\sum_{i=1}^{3}\frac{\psi
w_{i}}{2}\int_{0}^{a}\chi(x)\frac{\exp{\left(-\frac{x}{\lambda_{SE}}\right)}}{\lambda_{JE}\left({x+\frac{\lambda_{SE}}{\lambda_{JE}}}\right)}dx,$
(17) $\displaystyle{I}_{2}$ $\displaystyle=\sum_{i=1}^{3}\frac{\psi
w_{i}}{2}\int_{0}^{a}\chi(x)\frac{{\lambda_{SE}}\exp{\left(-\frac{x}{\lambda_{SE}}\right)}}{{\lambda_{JE}}\left({x+\frac{\lambda_{SE}}{\lambda_{JE}}}\right)^{2}}dx.$
(18)

The solutions of the integrals in (17) and (18) are similar. After some
algebraic manipulations, $I_{1}$ is expressed as

$\displaystyle{I}_{1}$ $\displaystyle=\sum_{i=1}^{3}\frac{\psi
w_{i}\xi}{2\lambda_{JE}}\Psi,$ (19)

where

$\displaystyle\xi$
$\displaystyle=\exp\left[-\frac{p_{i}\mu^{2}}{2\sigma^{2}}+\frac{p_{i}(\rho-1)\lambda_{SE}\Gamma_{0}}{\rho
p_{i}\lambda_{SE}+2\alpha\sigma^{2}\Gamma_{0}}+\frac{\rho
p_{i}\lambda_{SE}+2\alpha\sigma^{2}\Gamma_{0}}{2\alpha\sigma^{2}\lambda_{SE}\Gamma_{0}}\right.$
$\displaystyle\left.\times\left(\frac{\alpha\rho
p_{i}^{2}\mu^{2}\lambda_{SE}^{2}\Gamma_{0}}{\left(\rho
p_{i}\lambda_{SE}+2\alpha\sigma^{2}\Gamma_{0}\right)^{2}}+\frac{\rho-1}{\rho}\right)\right],$
(20)

and

$\displaystyle\Psi$
$\displaystyle=\int_{\sqrt{\rho-1}}^{\sqrt{a\rho+\rho-1}}\frac{x\exp{\left(-b\left(x-c\right)^{2}\right)}}{x^{2}+d}dx,$
(21)

with $b=\frac{\rho
p_{i}\zeta_{RE}\zeta_{SR}+2\sigma^{2}}{2\alpha\rho\zeta_{RE}\zeta_{SR}\sigma^{2}\Gamma_{0}N}$,
$c=\frac{\sqrt{\alpha\Gamma_{0}}\rho p_{i}\mu\zeta_{RE}\zeta_{SR}}{\rho
p_{i}\zeta_{RE}\zeta_{SR}+2\sigma^{2}}$, and
$d=1-\rho+\frac{\lambda_{SE}}{\lambda_{JE}}$. The solution of (21) is
difficult to achieve, hence, we find an approximation of the same. By noting
that $\exp{\big{(}-b\left(x-c\right)^{2}\big{)}}$ is symmetric around $x=c$
and if $c$ is large, which is typically the case when $N$ is large or
$\Gamma_{0}$ is large (high-SNR operation) or both, the product (21) tends to
be symmetric around $x=c$, we replace $\frac{x}{x^{2}+d}$ in (21) by its first
order Taylor series approximation around $x=c$, i.e.,
$\frac{x}{x^{2}+d}\approx\frac{c}{c^{2}+d}$. Then, $\Psi$ is easily
approximated as

$\displaystyle\Psi=\frac{c\sqrt{\pi}}{2\sqrt{b}(c^{2}+d)}\left[\text{erf}\left(\sqrt{b}\left(\sqrt{a\rho+\rho-1}-c\right)\right)\right.$
$\displaystyle\left.-\text{erf}\left(\sqrt{b}\left(\sqrt{\rho-1}-c\right)\right)\right].$
(22)

In the evaluation of $I_{2}$, we follow the same steps as were used while
evaluating $I_{1}$ and use the approximation
$\frac{x}{x^{2}+d^{2}}\approx\frac{c}{(c^{2}+d)^{2}}$. Therefore, $I_{2}$ is
evaluated as

$\displaystyle I_{2}=\sum_{i=1}^{3}\frac{\psi
c\sqrt{\pi}w_{i}\xi\lambda_{SE}}{2\sqrt{b}(c^{2}+d)^{2}\lambda_{JE}}\left[\text{erf}\left(\sqrt{b}\left(\sqrt{a\rho+\rho-1}-c\right)\right)\right.$
$\displaystyle-\left.\text{erf}\left(\sqrt{b}\left(\sqrt{\rho-1}-c\right)\right)\right].$
(23)

We observe from (15), (III), (19) and (III) that the closed-form approximate
SOP does not provide adequate information about how it behaves with the
different system parameters $\alpha$, $\rho$, $N$, $\Gamma_{0}$ and
$\zeta_{XY}$ where $XY\in\\{SR,RD,JR,RE\\}$. To obtain better insights into
the SOP performance, we use a further approximation in the following section.
This will also help us to obtain the optimal power allocation factor $\alpha$
in the next section.

##  IV Optimal Power Allocation

In this section, the optimal power allocation factor, $\alpha^{*}$, is derived
which minimizes the approximate SOP. To find $\alpha^{*}$, we first prove the
convexity of the approximate SOP with respect to $\alpha$. This is difficult
to prove from the approximate SOP derived in (12) in the previous section;
thus, we use further approximations on the SOP assuming $N$ is large. The
approximate SOP in this section also provides insights about the system
behaviour. The integrals in (III), (17), and (18) are respectively evaluated
by assuming $a\rightarrow\infty$ when $N$ is large following the same approach
considered in the previous section as

$\displaystyle I_{0}$
$\displaystyle=\psi,I_{1}=\psi\sum_{i=1}^{3}\frac{w_{i}c\rho\xi\sqrt{\pi}}{2\sqrt{b}(c^{2}+d)\lambda_{JE}},$
$\displaystyle I_{2}$ $\displaystyle=\sum_{i=1}^{3}\frac{\psi
w_{i}c\xi\sqrt{\pi}\lambda_{SE}}{2\sqrt{b}(c^{2}+d)^{2}\lambda_{JE}}.$ (24)

After using $I_{0}$, $I_{1}$, and $I_{2}$ from (IV) in (15) and substituting
for the values of $b$, $c$, $\lambda_{SE}$ and $\lambda_{JE}$, the approximate
SOP is written assuming $R_{th}>0$, i.e., $\rho>1$, as

$\displaystyle\mathcal{P}_{out}$
$\displaystyle=\sum_{i=1}^{3}\frac{\Omega_{1}\Omega_{2}}{1-\alpha}\left(1+\frac{\rho\lambda_{SE}}{c^{2}+d}\right)\exp\left(\frac{\rho-1}{\alpha\rho\zeta_{RE}\zeta_{SR}\Gamma_{0}N}\right),$
(25)

where

$\displaystyle\Omega_{1}=\frac{\psi
w_{i}\sqrt{\left({16-\pi^{2}}\right)}}{2\rho
p_{i}{N^{3/2}}\sqrt{\pi}\Gamma_{0}\zeta_{JR}\zeta_{RE}}\sqrt{{2\rho^{2}p_{i}+\frac{(16-\pi^{2})\rho\zeta_{RD}}{4\zeta_{RE}}}}$
(26)
$\displaystyle\Omega_{2}=\exp\left[-\frac{p_{i}\pi^{2}N}{2\left({16-\pi^{2}}\right)}\left(\frac{{1}}{\frac{8\rho
p_{i}\zeta_{RE}}{(16-\pi^{2})\zeta_{RD}}+{1}}\right)\right].$ (27)

To further simplify the equation, the term
$\frac{\rho\lambda_{SE}}{(c^{2}+d)}$ is neglected in (25) as
$\frac{\rho\lambda_{SE}}{(c^{2}+d)}<<1$ when $N$ is large. Finally, the
approximate SOP in expressed in a compact form as a function of $\alpha$ as

$\displaystyle\mathcal{P}_{out}=\sum_{i=1}^{3}\frac{\Omega_{1}\Omega_{2}}{1-\alpha}\exp\left(\frac{\rho-1}{\alpha\rho\zeta_{RE}\zeta_{SR}\Gamma_{0}N}\right).$
(28)

Now, it is easy to understand how the SOP behaves with the system parameters
$\alpha$, $N$, $\Gamma_{0}$ and $\zeta_{XY}$ where $XY\in\\{SR,RD,JR,RE\\}$.
We observe from (28) along with (26) and (27) that the SOP is inversely
proportional to $\Gamma_{0}$. This is because $\Omega_{1}$ is inversely
proportional to $\Gamma_{0}$ and the exponential term in (28) tends to unity
as $N$ becomes large. From (26) and (27) we observe that the SOP decreases
with $N$ at a rate of $\exp(-N)/N^{3/2}$ as the exponential term in (28) tends
to unity as $N$ increases. As the distance between the RIS and the
eavesdropper $d_{RE}$ increases, the SOP decreases exponentially. This is
because $\zeta_{RE}$ is inside the exponential in (27) and dominates over (26)
and (28) affecting the SOP. As the distance between the RIS and the
destination $d_{RD}$ increases, from (27) we observe that the SOP increases
exponentially and it dominates (26) which decreases slowly. As the distance
between the jammer and the RIS $d_{JR}$ increases, the SOP increases as we
observe from (26). As the distance between the transmitter and the RIS
$d_{SR}$ changes, the SOP does not change as $\Omega_{1}$ and $\Omega_{2}$ are
independent of $\zeta_{SR}$. Also, although $\zeta_{SR}$ is present in (28),
it does not affect the SOP as the exponential term in (28) tends to unity as
$N$ increases.

Finally, we observe from (28) that as $\alpha$ increases, the exponential term
in (28) decreases and hence, the SOP decreases. In contrast, the denominator
$(1-\alpha)$ decreases and hence, the SOP increases. This suggests that there
is an optimal value of $\alpha$ which can minimize the SOP. The optimal
$\alpha$ is derived in the next section.

###  IV-A Finding the optimal power allocation

To check the convexity of the approximate SOP with respect to $\alpha$, we
need to double differentiate the approximate SOP equation in (28). The first
derivative of $\mathcal{P}_{out}$ is given by

$\displaystyle\frac{d\mathcal{P}_{out}}{d\alpha}$
$\displaystyle=\sum_{i=1}^{3}\frac{\Omega_{1}\Omega_{2}\exp\left({\frac{\rho-1}{\alpha\rho\zeta_{RE}\zeta_{SR}\Gamma_{0}N}}\right)}{1-\alpha}$
$\displaystyle\times\left(\frac{1}{1-\alpha}-\frac{\rho-1}{\alpha^{2}\rho\zeta_{RE}\zeta_{SR}\Gamma_{0}N}\right),$
(29)

and the second derivative is evaluated as

$\displaystyle\frac{d^{2}\mathcal{P}_{out}}{d\alpha^{2}}=\sum_{i=1}^{3}\frac{\Omega_{1}\Omega_{2}\exp\left({\frac{\rho-1}{\alpha\rho\zeta_{RE}\zeta_{SR}\Gamma_{0}N}}\right)}{1-\alpha}\left(\frac{2-(1-\alpha)^{2}}{(1-\alpha)^{2}}\right.$
$\displaystyle\left.+\frac{(\rho-1)^{2}}{\alpha^{4}\rho^{2}\zeta_{RE}^{2}\zeta_{SR}^{2}N^{2}\Gamma_{0}^{2}}+\frac{2(\rho-1)}{(1-\alpha)\alpha^{3}\rho\zeta_{RE}\zeta_{SR}\Gamma_{0}N}\right).$
(30)

We observe that since the parameters $\rho$, $\zeta_{RE}$, $\zeta_{SR}$, $N$,
and $\Gamma_{0}$ are all positive and $0<\alpha<1$ and $\rho>1$, (IV-A) is
positive for all possible values of $\alpha$ and hence the function is convex
with respect to $\alpha$ and a global optimal minimum value exists.

Now, equating the first derivative of $\mathcal{P}_{out}$ in (IV-A) to zero,
we can find $\alpha^{*}$. The expression in (IV-A) will be equal to zero when
the term inside the brackets is equal to zero. After some mathematical
manipulations, this condition reduces to the quadratic equation as

$\displaystyle\alpha^{2}\rho\zeta_{RE}\zeta_{SR}\Gamma_{0}N+\alpha(\rho-1)-(\rho-1)=0.$
(31)

Since $\alpha^{*}$ cannot be a negative quantity, we take the positive root
for $\alpha$ and therefore the optimal power allocation factor is given by

$\displaystyle\alpha^{*}=\frac{-(\rho-1)+\sqrt{(\rho-1)^{2}+4\rho\zeta_{RE}\zeta_{SR}\Gamma_{0}N(\rho-1)}}{2\rho\zeta_{RE}\zeta_{SR}\Gamma_{0}N}.$
(32)

We observe from (32) that the optimal power allocation depends on parameters
$\rho$, $\zeta_{SR}$ and $\zeta_{RE}$, $N$ and $\Gamma_{0}$. However, it is
independent of $\zeta_{JR}$ and $\zeta_{RD}$. As $\rho$, $\zeta_{SR}$ and
$\zeta_{RE}$, $N$, and $\Gamma_{0}$ increase, the jamming power needs to be
increased in order to minimize the SOP.

##  V Numerical Results and Discussion

Figure 1: SOP versus $\Gamma_{0}$ for different values of $N$.

In this section, analytical results are plotted along with the numerical
results. Numerical results are denoted by a solid line with the symbol
‘$\times$’ having the same color of the corresponding analytical curve. The
system parameters, unless otherwise stated, are considered as $R_{th}=1$ bit
per channel use, $\\{d_{SR},d_{JR},d_{RD},d_{RE}\\}=\\{30,30,30,15\\}$ m, the
path loss in dB is given by $\zeta_{XY}(\text{dB})=z_{0}-10v\log_{10}(d_{XY})$
where $z_{0}=42$ dB and $v=3.5$ is the reference path loss and path loss
factor, respectively.

Fig. 1 plots the approximate SOP from (12) versus $\Gamma_{0}$ for different
number of RIS elements $N\in\\{16,32,64\\}$ with equal power allocation (EPA)
$\alpha=0.5$ and with the optimal power allocation (OPA) $\alpha^{*}$. We find
that the SOP decreases with the increase in the number of RIS elements. It is
observed that the OPA outperforms the EPA as expected. By using the OPA
factor, a gain of approximately $3$ dB is achieved at an SOP of $10^{-4}$ for
$N=64$.

Figure 2: SOP versus $\Gamma_{0}$ for different values of $d_{RD}/d_{RE}$.

In Fig. 2, the approximate SOP from (12) is plotted versus $\Gamma_{0}$, by
varying the ratio of distances $R$-$D$ and $R$-$E$, i.e.,
$\frac{d_{RD}}{d_{RE}}$ with $\\{d_{SR},d_{JR},d_{RE}\\}$ = $\\{30,30,15\\}$
m. It is observed that as the ratio increases, the secrecy performance
decreases. This signifies that as the $R$-$D$ distance increases, or as the
$R$-$E$ distance decreases, the SOP increases.

Figure 3: Optimal $\alpha^{*}$ for different numbers of RIS elements $N$.

Fig. 3 shows the effect of $N$ on $\alpha^{*}$ at different $\Gamma_{0}$. It
is observed that as $N$ increases, $\alpha^{*}$ decreases, i.e., more power is
allocated to the jammer. This is due to the fact that as the number of RIS
elements increases, the reflected signal strength increases for the
eavesdropper as well. The jamming power has to be increased in order to
minimize the SOP. We also notice that as $\Gamma_{0}$ increases for a given
$N$, $\alpha^{*}$ decreases. This means that as the SNR improves, less power
is required for signal transmission to minimize the SOP.

Figure 4: Optimal $\alpha^{*}$ versus $d_{RE}$ for different locations of $E$.

In Fig. 4, $\alpha^{*}$ versus $d_{RE}$ is plotted for different $d_{SR}$
values. With the increase in $d_{RE}$ for a given $d_{SR}$ or with the
increase in $d_{SR}$ for a given $d_{RE}$, $\alpha^{*}$ increases, i.e., more
power is allocated for signal transmission. This suggests that as distance of
$S$ increases from $R$, more signal power is required and as distance of $E$
increases from $R$, less jamming power is required to achieve the minimum SOP.

##  VI Conclusions

The secrecy performance of an RIS-assisted wireless communication has been
analyzed in the presence of an eavesdropper. A friendly jammer is introduced
to enhance the secrecy performance of the system. The approximate closed-form
SOP of the system is evaluated. A simpler approximate closed-form SOP
expression is also provided which shows the effect of the system parameters on
the performance. This also helps us to show the convexity of the approximate
SOP with the power allocation factor. The power allocation factor for the
transmitter and the jammer is derived in the closed-from. We find that the
optimal power allocation achieves approximately 3 dB performance gain as
compared to the case of equal power allocation at an SOP of $10^{-4}$ for
$N=64$.

## References

  * [1] Q. Wu and R. Zhang, “Towards smart and reconfigurable environment: intelligent reflecting surface aided wireless network,” _IEEE Communications Magazine_ , vol. 58, no. 1, pp. 106–112, Jan. 2020.
  * [2] E. Basar, M. Di Renzo, J. De Rosny, M. Debbah, M.-S. Alouini, and R. Zhang, “Wireless communications through reconfigurable intelligent surfaces,” _IEEE Access_ , vol. 7, pp. 116 753–116 773, Aug. 2019.
  * [3] W. Wu, Z. Wang, L. Yuan, F. Zhou, F. Lang, B. Wang, and Q. Wu, “IRS-enhanced energy detection for spectrum sensing in cognitive radio networks,” _IEEE Wireless Communications Letters_ , vol. 10, no. 10, pp. 2254–2258, Oct. 2021.
  * [4] X. Gan, C. Zhong, Y. Zhu, and Z. Zhong, “User selection in reconfigurable intelligent surface assisted communication systems,” _IEEE Communications Letters_ , vol. 25, no. 4, pp. 1353–1357, Apr. 2021.
  * [5] D. Kudathanthirige, D. Gunasinghe, and G. Amarasuriya, “Performance analysis of intelligent reflective surfaces for wireless communication,” in _Proc. IEEE International Conference on Communications_ , Jul. 2020, pp. 1–6.
  * [6] A. D. Wyner, “The wire-tap channel,” _Bell System Technical Journal_ , vol. 54, no. 8, pp. 1355–1387, Oct. 1975.
  * [7] X. Yu, D. Xu, and R. Schober, “Enabling secure wireless communications via intelligent reflecting surfaces,” in _Proc. IEEE Global Communications Conference_ , Dec. 2019, pp. 1–6.
  * [8] X. Yu, D. Xu, Y. Sun, D. W. K. Ng, and R. Schober, “Robust and secure wireless communications via intelligent reflecting surfaces,” _IEEE Journal on Selected Areas in Communications_ , vol. 38, no. 11, pp. 2637–2652, Nov. 2020\.
  * [9] M. Cui, G. Zhang, and R. Zhang, “Secure wireless communication via intelligent reflecting surface,” _IEEE Wireless Communications Letters_ , vol. 8, no. 5, pp. 1410–1414, Oct. 2019.
  * [10] V. P. Tuan and I. P. Hong, “Secrecy performance analysis and optimization of intelligent reflecting surface-aided indoor wireless communications,” _IEEE Access_ , vol. 8, pp. 109 440–109 452, Jun. 2020.
  * [11] L. Yang, J. Yang, W. Xie, M. O. Hasna, T. Tsiftsis, and M. D. Renzo, “Secrecy performance analysis of RIS-aided wireless communication systems,” _IEEE Transactions on Vehicular Technology_ , vol. 69, no. 10, pp. 12 296–12 300, Oct. 2020.
  * [12] R. Negi and S. Goel, “Secret communication using artificial noise,” in _Proc. IEEE Vehicular Technology Conference_ , vol. 3, Sep. 2005, pp. 1906–1910.
  * [13] L. Dong, Z. Han, A. P. Petropulu, and H. V. Poor, “Cooperative jamming for wireless physical layer security,” in _Proc. IEEE/SP Workshop on Statistical Signal Processing_ , Sep. 2009, pp. 417–420.
  * [14] T. M. Hoang, T. Q. Duong, N.-S. Vo, and C. Kundu, “Physical layer security in cooperative energy harvesting networks with a friendly jammer,” _IEEE Wireless Communications Letters_ , vol. 6, no. 2, pp. 174–177, Apr. 2017.
  * [15] K. Cumanan, G. C. Alexandropoulos, Z. Ding, and G. K. Karagiannidis, “Secure communications with cooperative jamming: optimal power allocation and secrecy outage analysis,” _IEEE Transactions on Vehicular Technology_ , vol. 66, no. 8, pp. 7495–7505, Aug. 2017.
  * [16] X. Guan, Q. Wu, and R. Zhang, “Intelligent reflecting surface assisted secrecy communication: Is artificial noise helpful or not?” _IEEE Wireless Communications Letters_ , vol. 9, no. 6, pp. 778–782, Jan. 2020.
  * [17] Q. Wang, F. Zhou, R. Q. Hu, and Y. Qian, “Energy efficient robust beamforming and cooperative jamming design for IRS-assisted MISO networks,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 4, pp. 2592–2607, Apr. 2021.
  * [18] S. Xu, J. Liu, and Y. Cao, “Intelligent reflecting surface empowered physical-layer security: Signal cancellation or jamming?” _IEEE Internet of Things Journal_ , vol. 9, no. 2, pp. 1265–1275, Jan. 2022.
  * [19] Y. Sun, K. An, J. Luo, Y. Zhu, G. Zheng, and S. Chatzinotas, “Intelligent reflecting surface enhanced secure transmission against both jamming and eavesdropping attacks,” _IEEE Transactions on Vehicular Technology_ , vol. 70, no. 10, pp. 11 017–11 022, Oct. 2021.
  * [20] J. Salo, H. El-Sallabi, and P. Vainikainen, “The distribution of the product of independent Rayleigh random variables,” _IEEE Transactions on Antennas and Propagation_ , vol. 54, no. 2, pp. 639–643, Feb. 2006.
  * [21] C. Kundu and M. F. Flanagan, “Ergodic secrecy rate of optimal source selection in a multi-source system with unreliable backhaul,” _IEEE Wireless Communications Letters_ , vol. 10, no. 5, pp. 1118–1122, May. 2021.
  * [22] C. Kundu, S. Ghose, and R. Bose, “Secrecy outage of dual-hop regenerative multi-relay system with relay selection,” _IEEE Transactions on Wireless Communications_ , vol. 14, no. 8, pp. 4614–4625, Aug. 2015.
  * [23] D. Sadhwani, R. N. Yadav, and S. Aggarwal, “Tighter bounds on the Gaussian $Q$ function and its application in Nakagami-${m}$ fading channel,” _IEEE Wireless Communications Letters_ , vol. 6, no. 5, pp. 574–577, Oct. 2017.


# Deviation of nonradial geodesics in a static spherically symmetric space-
time

V.P. Vandeev∗, A.N. Semenova+

###### Abstract

The article generalizes the description of tidal forces to the case of
geodesics with non-zero angular momentum in the metric of static spherically
symmetric black holes. We show that the geodesic deviation equation can be
diagonalized even with non-radial free motion of a test body in the
gravitational field. We present expressions for the spatial components of the
tidal force in a spherically symmetric metric. We also solve geodesic
deviation equation in the Schwarzschild metric and demonstrate how the
presence of angular momentum and its magnitude affect the solution.

∗BALTIC STATE TECHNICAL UNIVERSITY «VOENMEH» named after D.F. Ustinov,
1-Ya Krasnoarmeyskaya Ulitsa, 1, St Petersburg, 190005, Russia

+Petersburg Nuclear Physics Institute of National Research Centre ‘‘Kurchatov
Institute’’,
Gatchina, 188300, Russia

##  1 Introduction

The black hole properties are interesting for modern scientists, despite the
fact that the first solution to describe a black hole was found in 1916 by
Karl Schwarzschild [1], which has a single parameter – the mass of the black
hole, they are still relevant both for theoretical and experimental research.
Further, many solutions of Einstein’s equations appeared, which generalized
Schwarzschild’s solution: Reissner–Nordström (electrically charged) black hole
[2], [3], Kerr (rotating) black hole [4], Kerr–Newman (charged and rotating)
black hole [5]. Of particular interest are the properties of solutions to
Einstein’s equations in the presence of various matter. Examples are Kottler
metric [6], which solves the GR equations when the cosmological constant plays
the role of matter and Kiselev metric, where matter is represented by the
quintessence [7].

The main subject of interest for us is geodesic deviation. The first work
devoted to solving geodesic deviation equation was [8], there the behavior of
the geodesic deviation vector in the Schwarzschild metric was obtained.
Despite the fact that modern studies of tidal forces are devoted to exotic
solutions of the Einstein’s equations: tidal effects in regular black holes
[9], tidal force in Schwarzschild spacetime in holographic massive gravity
[10], deviation vector behavior in 4D-Einstein–Gauss–Bonnet spacetime [11], in
this paper we concentrate on solving the problem of tidal forces in the metric
of a spherically symmetric spacetime. However unlike articles [12]-[14] we
consider non-radial geodesics.

Tidal effects play an important role in astrophysical phenomena. For example,
tidal disruption event (or sometimes called tidal disruption flare) is that
the star is destructed by the gravity of a supermassive black hole. Tidal
forces acting on a star in the Schwarzschild gravitational field are studied
in [15]. The study of tidal effects in axisymmetric spactimes in the vicinity
of rotation axis, for example, in the Kerr metric [16], can make it possible
to describe the formation of astrophysical jets. The description of
relativistic flow formation through tidal acceleration analysis is presented
in [17], [18], [19].

In our work the geodesic deviation equation depends on two integrals of
motion: energy $E$ and angular momentum $L$. We demonstrate how geodesic
deviation equation can be diogonalized for static spherically symmetric metric
defined by one function of radial variable. Such class of metrics include
Schwarzschild, Reissner–Nordström, Kottler, Kiselev spacetime etc. Geodesic
deviation equation with nonzero angular momentum have not previously been
considered in works devoted to the study of tidal forces in static metrics.
Probably the reason for this is the impossibility of expressing the solutions
of differential equations of geodesic deviation in elementary functions. We
found that these equations satisfy Fuchs’ theorem, which allows us to
construct solutions in the vicinity of the points of interest using the
Frobenius method. It is well known that in Schwarzschild spacetime a test body
falling towards the event horizon of BH experiences stretching along the
radial direction and compression along the both angular directions [20]. We
give a quantitative answer to the question of the angular momentum influence
on tidal force components and on the behavior of geodesic deviation equation
solution. The effect of angular momentum on tidal forces is especially
important for description of tidal destruction event, because the motion of
stars in the gravitational field of a black hole is non-radial. It should also
be noted that the relative tidal acceleration between the two gravitational
wave detectors was used to study the properties of gravitons [21], [22].

The article is organized as follows. Seс. 2 describes a system of geodesic
equations that helps to define the velocity 4-vector of a freely falling test
body. In Seс. 3 we consider the geodesic deviation equation and describe the
procedure for its diagonalization in a spherically symmetric spacetime. In
Sec. 4 we apply the expressions from the previous sections to the
Schwarzschild metric, then we solve geodesic deviation equation and compare
the solutions for radial geodesics ($L=0$) with solutions for non-radial
curves of free motion ($L\neq 0$). In the subsections of Sec. 4 we present the
numerical solutions of the geodesic deviation equation, local solutions in the
vicinity of zero and infinity in the form of generalized power series. Sec. 5
concludes our study, there we summarize and emphasize new results of our work
and describe the range of problems that can be solved within the framework of
the tidal forces study following our work. We use the metric signature
$(+,-,-,-)$ and set the speed of light $c$ and Newtonian gravitational
constant $G$ to $1$ throughout this paper.

##  2 Geodesics in static spherically symmetric spacetime

We consider spacetime with line element of static spherically symmetric black
hole given by

$ds^{2}=g_{\mu\nu}dx^{\mu}dx^{\nu}=f(r)dt^{2}-\frac{dr^{2}}{f(r)}-r^{2}\left(d\theta^{2}+\sin^{2}\theta
d\phi^{2}\right).$ (1)

At this stage we do not specify the form of the function $f(r)$. In this
metric there are four meaningful geodesic equations [23] for all coordinates:

$u^{0}=\frac{dt}{d\tau}=\frac{E}{f(r)},$ (2)
$\left(u^{1}\right)^{2}=\left(\frac{dr}{d\tau}\right)^{2}=E^{2}-f(r)\left(1+\frac{L^{2}}{r^{2}}\right),$
(3) $u^{2}=\frac{d\theta}{d\tau}=0,$ (4)
$u^{3}=\frac{d\phi}{d\tau}=\frac{L}{r^{2}},$ (5)

where $\tau$ is the proper time, $E$ and $L$ are energy and angular momentum
of a freely moving test body respectively. And it should be noted that
dynamics of the azimuthal variable is trivial due to the spherical symmetry of
spacetime, and the dynamics of the polar variable is nontrivial because there
is angular momentum $L$. The set of expressions (2)-(5) forms a unit covariant
4-velocity vector $u^{\mu}$ tangent to the geodesic. Also it should be noted
that work considers only timelike geodesics.

##  3 Geodesic deviation equation

Below we consider geodesic deviation equation. As it is well known [20], the
equation for the spacelike geodesic deviation vector $\tilde{\xi}^{\mu}$ is
given by

$\frac{D^{2}\tilde{\xi}^{\mu}}{d\tau^{2}}=R^{\mu}_{.\nu\alpha\beta}u^{\nu}u^{\alpha}\tilde{\xi}^{\beta},$
(6)

where $\frac{D^{2}}{d\tau^{2}}$ is covariant derivative along the geodesic,
$R^{\mu}_{.\nu\alpha\beta}$ is Riemann curvature tensor and $u^{\nu}$ is the
unit vector of 4-velocity tangent to the geodesic. $\tilde{\xi}^{\mu}$
describes the distance between two infinitesimally close free moving
particles.

The nonzero components of the Riemann tensor are calculated by the metric (1)
and have the form

$R^{0}_{.101}=-\frac{f^{\prime\prime}}{2f},\\\ $ (7a)
$R^{0}_{.202}=R^{1}_{.212}=-\frac{f^{\prime}r}{2},\\\ $ (7b)
$R^{0}_{.303}=R^{1}_{.313}=-\frac{f^{\prime}r\sin^{2}{\theta}}{2},\\\ $ (7c)
$R^{2}_{.323}=\left(1-f\right)\sin^{2}{\theta}.$ (7d)

It is seen that on the right-hand side of the Eq. (6) there is a matrix
$P^{\mu}_{.\beta}\equiv R^{\mu}_{.\nu\alpha\beta}u^{\nu}u^{\alpha}$ which is

$\begin{pmatrix}\frac{\dot{r}^{2}f^{\prime\prime}}{2f}+\frac{\chi
L^{2}\sin^{2}\theta}{r^{2}}&-\frac{E\dot{r}f^{\prime\prime}}{2f^{2}}&0&-\frac{\chi
EL\sin^{2}\theta}{f}\\\ \frac{E\dot{r}f^{\prime\prime}}{2}&\frac{\chi
L^{2}\sin^{2}\theta}{r^{2}}-\frac{E^{2}f^{\prime\prime}}{2f}&0&-\chi
L\dot{r}\sin^{2}\theta\\\
0&0&\frac{(f-1)L^{2}\sin^{2}\theta}{r^{4}}-\chi\omega&0\\\ \frac{\chi
EL}{r^{2}}&-\frac{\chi L\dot{r}}{fr^{2}}&0&-\chi\omega\end{pmatrix},$ (8)

where

$\chi=\frac{f^{\prime}}{2r},$ (9a) $\omega=1+\frac{L^{2}}{r^{2}},$ (9b)
$\dot{r}=u^{1}=\frac{dr}{d\tau}=\sqrt{E^{2}-f\left(1+\frac{L^{2}}{r^{2}}\right)}.$
(9c)

Matrix (8) has spectrum

$\lambda_{t}=0,$ (10a) $\lambda_{r}=\frac{\chi
L^{2}\sin^{2}\theta}{r^{2}}-\frac{f^{\prime\prime}\omega}{2},$ (10b)
$\lambda_{\theta}=\frac{(f-1)L^{2}\sin^{2}\theta}{r^{4}}-\chi\omega,$ (10c)
$\lambda_{\phi}=\frac{\chi L^{2}\sin^{2}\theta}{r^{2}}-\chi\omega.$ (10d)

Therefore, a tetrad basis for the free fall frame of reference can be
constructed. It has form

$e_{t}^{\mu}=\frac{1}{\sqrt{1+\frac{L^{2}\cos^{2}\theta}{r^{2}}}}\left(\frac{E}{f},\dot{r},0,\frac{L}{r^{2}}\right),$
(11a)
$e_{r}^{\mu}=\frac{1}{\sqrt{1+\frac{L^{2}}{r^{2}}}}\left(-\frac{\dot{r}}{f},-E,0,0\right),$
(11b) $e_{\theta}^{\mu}=\left(0,0,\frac{1}{r},0\right),$ (11c)
$e_{\phi}^{\mu}=\frac{L\sin\theta}{r\sqrt{1+\frac{L^{2}\cos^{2}\theta}{r^{2}}}\sqrt{1+\frac{L^{2}}{r^{2}}}}\left(\frac{E}{f},\dot{r},0,\frac{1+\frac{L^{2}}{r^{2}}}{L\sin^{2}{\theta}}\right).\\\
$ (11d)

This tetrad set $e^{\mu}_{\alpha}$ satisfy normalization condition
${e^{\mu}_{\alpha}\>e^{\nu}_{\beta}\>g_{\mu\nu}=\eta_{\alpha\beta}}$ with
Minkowski metric $\eta_{\alpha\beta}=\text{diag}\left(1,-1,-1,-1\right)$. The
geodesic deviation vector $\tilde{\xi}^{\mu}$ can be substituted as

$\tilde{\xi}^{\mu}=e^{\mu}_{\nu}\>\xi^{\nu}.$ (12)

Thus, the meaningful components of Eq. (6) in the new coordinate system are

$\ddot{\xi}^{r}=\left[-\frac{f^{\prime\prime}}{2}\left(1+\frac{L^{2}}{r^{2}}\right)+\frac{f^{\prime}}{2r}\frac{L^{2}}{r^{2}}\right]\xi^{r},$
(13)
$\ddot{\xi}^{\theta}=\left[-\frac{f^{\prime}}{2r}\left(1+\frac{L^{2}}{r^{2}}\right)+\frac{f-1}{r^{2}}\frac{L^{2}}{r^{2}}\right]\xi^{\theta},$
(14) $\ddot{\xi}^{\phi}=-\frac{f^{\prime}}{2r}\xi^{\phi}.$ (15)

These equations are the diagonal form of Eq. (6). It is worth noting that the
dependence on the $\theta$ disappeared from Eqs. (13)–(15) because there is no
dynamics along the direction of the azimuthal angle $\theta$ according to Eq.
(4) and geodesics lie in the equatorial plane $\theta=\frac{\pi}{2}$.

The diagonal form of Eq. (6) for a static spherically symmetric spacetime for
non-radial geodesics was obtained in [8] but a slightly more complex
mathematical apparatus was used there, although in fact it is sufficient to
use the methods of linear algebra. For $L=0$ Eqs. (13)–(15) coincide with Eqs.
(14) and (15) of Ref. [12].

To illustrate the result obtained, we can use superposition of charged black
hole and (Anti)-de Sitter spacetime, its metric is

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

where

$f(r)=1-\frac{2M}{r}+\frac{q^{2}}{r^{2}}-\frac{\Lambda}{3}r^{2}.$ (17)

Using (17) in Eqs. (13)–(15) we get

$\ddot{\xi}^{r}=\left[\frac{2M}{r^{3}}-\frac{3q^{2}}{r^{4}}+\frac{\Lambda}{3}+\frac{L^{2}}{r^{2}}\left(\frac{3M}{r^{3}}-\frac{4q^{2}}{r^{4}}\right)\right]\xi^{r},$
(18a)
$\ddot{\xi}^{\theta}=\left[-\frac{M}{r^{3}}+\frac{q^{2}}{r^{4}}+\frac{\Lambda}{3}+\frac{L^{2}}{r^{2}}\left(-\frac{3M}{r^{3}}+\frac{2q^{2}}{r^{4}}\right)\right]\xi^{\theta},$
(18b)
$\ddot{\xi}^{\phi}=\left[-\frac{M}{r^{3}}+\frac{q^{2}}{r^{4}}+\frac{\Lambda}{3}\right]\xi^{\phi}.$
(18c)

##  4 Tidal force

Expressions (18) are rather complicated due to the large number of parameters,
therefore, below we restrict ourselves to an electrically neutral black hole
with no cosmological constant ($q=0$ and $\Lambda=0$). In other words, we
consider the Schwarzschild metric

$\ddot{\xi}^{r}=\left[\frac{2M}{r^{3}}+\frac{3ML^{2}}{r^{5}}\right]\xi^{r},$
(19)
$\ddot{\xi}^{\theta}=\left[-\frac{M}{r^{3}}-\frac{3ML^{2}}{r^{5}}\right]\xi^{\theta},$
(20) $\ddot{\xi}^{\phi}=-\frac{M}{r^{3}}\xi^{\phi}.$ (21)

This case for $L=0$ was considered in [20]. Note that the presence of angular
momentum makes the expressions for the angular components (20) and (21) of the
tidal force different. The presence of terms with angular momentum in the
radial (19) and azimuthal (20) components of the tidal force does not change
the monotonicity of these components depending on the radius $r$. This can be
seen in Figs. 1 and 2. Also it is obvious that all components of the tidal
force are constant sign and do not vanish for finite $r$.

Figure 1: Radial tidal force component for Schwarzschild BH with different
choices of $L$. Figure 2: Azimuthal tidal force component for Schwarzschild BH
with different choices of $L$.

There is no need to depict the polar component of the tidal force (21) because
it is identical to the azimuthal one at $L=0$. It looks the same as the black
solid curve in Fig. 2. Also it should be noted that in the case of
Schwarzschild metric we do not consider very large values of the angular
momentum $L$ because Eq. (3) for $L\to\infty$ has form of a power series

$\frac{dr}{d\tau}=\frac{L}{r}\sqrt{\frac{2M-r}{r}}+O\left(L^{-1}\right),$ (22)

which is defined only under the black hole horizon for $r\leq 2M$. The absence
of geodesics above the horizon $r>2M$ naturally leads to the absence of their
deviation.

###  4.1 Solutions of the geodesic deviation equations

Using Eq. (3) we pass from differentiation with respect to proper time $\tau$
to differentiation with respect to dimensionless radial variable
$\rho=\frac{r}{M}$ in Eqs. (19)–(21)

${\xi^{r}}^{\prime\prime}+Q(\rho){\xi^{r}}^{\prime}+P_{r}(\rho)\xi^{r}=0,$
(23)
${\xi^{\theta}}^{\prime\prime}+Q(\rho){\xi^{\theta}}^{\prime}+P_{\theta}(\rho)\xi^{\theta}=0,$
(24)
${\xi^{\phi}}^{\prime\prime}+Q(\rho){\xi^{\phi}}^{\prime}+P_{\phi}(\rho)\xi^{\phi}=0.$
(25)

Where $l=\frac{L}{M}$ is dimensionless angular momentum and the following
notation is used for abbreviation

$Q(\rho)=-\frac{\rho^{2}-l^{2}\rho+3l^{2}}{(E^{2}-1)\rho^{4}+2\rho^{3}-l^{2}\rho^{2}+2l^{2}\rho},$
(26a)
$P_{r}(\rho)=-\frac{2\rho^{2}+3l^{2}}{(E^{2}-1)\rho^{5}+2\rho^{4}-l^{2}\rho^{3}+2l^{2}\rho^{2}},$
(26b)
$P_{\theta}(\rho)=\frac{\rho^{2}+3l^{2}}{(E^{2}-1)\rho^{5}+2\rho^{4}-l^{2}\rho^{3}+2l^{2}\rho^{2}},$
(26c) $P_{\phi}(\rho)=\frac{1}{(E^{2}-1)\rho^{3}+2\rho^{2}-l^{2}\rho+2l^{2}}.$
(26d)

It should be noted that for $l=0$ the solutions of the Eqs. (23)–(25) are
expressed in terms of elementary functions

$\xi^{r}=A\alpha(\rho)+\frac{B}{\gamma^{4}}\left[6+\gamma^{2}\rho+\frac{3\alpha(\rho)}{\gamma}\ln\left(\frac{\alpha(\rho)-\gamma}{\alpha(\rho)+\gamma}\right)\right],$
(27) $\xi^{\theta}=\xi^{\phi}=\rho\bigg{(}C+D\alpha(\rho)\bigg{)},$ (28)

where

$\gamma=\sqrt{E^{2}-1},\;\alpha(\rho)=\sqrt{E^{2}-1+\frac{2}{\rho}},$ (29)

and $A,B,C$ and $D$ are integration constants. These solutions were firstly
obtained in Ref. [8].

These solutions in the vicinity of $\rho=0$ are

$\xi^{r}(\rho)=\frac{A\sqrt{2}}{\sqrt{\rho}}+O\left(\sqrt{\rho}\right)\;\text{for}\;\rho\to
0,$ (30a)
$\xi^{\theta}(\rho)=\sqrt{2}D_{\theta}\sqrt{\rho}+O\left(\rho\right)\;\text{for}\;\rho\to
0,$ (30b)
$\xi^{\phi}(\rho)=\sqrt{2}D_{\phi}\sqrt{\rho}+O\left(\rho\right)\;\text{for}\;\rho\to
0.$ (30c)

The same solutions in the vicinity of $\rho=\infty$ behave

$\xi^{r}(\rho)=\frac{B}{E^{2}-1}\rho+O\left(1\right)\;\text{for}\;\rho\to\infty,$
(31a)
$\xi^{\theta}(\rho)=\left(C_{\theta}+D_{\theta}\sqrt{E^{2}-1}\right)\rho+O\left(1\right)\;\text{for}\;\rho\to\infty,$
(31b)
$\xi^{\phi}(\rho)=\left(C_{\phi}+D_{\phi}\sqrt{E^{2}-1}\right)\rho+O\left(1\right)\;\text{for}\;\rho\to\infty.$
(31c)

Solutions of (23)–(25) can not be expressed in terms of elementary functions
for the case $l\neq 0$ so we will find them in the form of generalized series
in the vicinity of points $\rho=0$ and $\rho=\infty$. We will present a
comparison of such solutions for non-zero angular momentum and Eqs. (30) and
(31) in the following sections.

###  4.2 Numerical solutions

Nevertheless, the solutions of Eqs. (23)–(25) can be found numerically. In
Fig. 3 we construct a numerical solution of the radial Eq. (23) for various
values of the angular momentum $l$. It can be seen that the black solid curve
with $l=0$ in the vicinity of $\rho=0$ corresponds to (30a).

Figure 3: Dependence of $\xi^{r}$ on $\rho$ at $E=10$. Initial data are
$\xi^{r}(10)=1$ and ${\xi^{r}}^{\prime}(10)=0$.

In Fig. 4 we construct a numerical solution of the azimuthal Eq. (24) for
various values of the angular momentum $l$. It can be seen that the black
solid curve with $l=0$ in the vicinity of $\rho=0$ corresponds to (30b).

Figure 4: Dependence of $\xi^{\theta}$ on $\rho$ at $E=10$. Initial data are
$\xi^{\theta}(10)=1$ and ${\xi^{\theta}}^{\prime}(10)=0$.

In Fig. 5 we construct a numerical solution of the polar Eq. (25) for various
values of the angular momentum $l$. It can be seen that the black solid curve
with $l=0$ in the vicinity of $\rho=0$ corresponds to (30c).

Figure 5: Dependence of $\xi^{\phi}$ on $\rho$ at $E=10$. Initial data are
$\xi^{\phi}(10)=1$ and ${\xi^{\phi}}^{\prime}(10)=0$.

###  4.3 Local behavior of solutions

Below we solve Eqs. (23)–(25) in the vicinity of points $\rho=0$,
$\rho=\infty$. Since these equations cannot be solved in terms of elementary
functions, we will construct solutions in the form of power series.

####  4.3.1 In the BH singularity vicinity

Here we represent functions (26) in the form of a power series near $\rho=0$

$Q(\rho)=-\frac{3}{2\rho}-\frac{1}{4}+\frac{8-l^{2}}{8l^{2}}\rho+\frac{12E^{2}-l^{2}}{16l^{2}}\rho^{2}+O\left(\rho^{3}\right),$
(32a)
$P_{r}(\rho)=-\frac{3}{2\rho^{2}}-\frac{3}{4\rho}+\frac{4-3l^{2}}{8l^{2}}+\frac{12E^{2}-3l^{2}+4}{16l^{2}}\rho+O\left(\rho^{2}\right),$
(32b)
$P_{\theta}(\rho)=\frac{3}{2\rho^{2}}+\frac{3}{4\rho}-\frac{8-3l^{2}}{8l^{2}}-\frac{12E^{2}-3l^{2}+8}{16l^{2}}\rho+O\left(\rho^{2}\right),$
(32c)
$P_{\phi}(\rho)=\frac{1}{2l^{2}}+\frac{\rho}{4l^{2}}+\frac{l^{2}-4}{8l^{4}}\rho^{2}-\frac{4E^{2}-l^{2}+4}{16l^{4}}\rho^{3}+O\left(\rho^{4}\right).$
(32d)

It follows from the form of these series that the Eqs. (23)–(25) satisfy
Fuchs’ theorem. This means that their solutions can be found using generalized
power series by the Frobenius method. Therefore, each of them has two linearly
independent solutions

$\xi^{j}_{1}=\sum_{n=0}^{\infty}c^{j}_{n}\rho^{n+\zeta_{1}},$ (33a)
$\xi^{j}_{2}=\sum_{n=0}^{\infty}d^{j}_{n}\rho^{n+\zeta_{2}},$ (33b)

where $j=r,\theta,\phi$. And $\zeta_{1},\zeta_{2}$ are roots of the quadratic
equation $\zeta(\zeta-1)+a^{j}_{-1}\zeta+b^{j}_{-2}=0$ (where $a^{j}_{-1}$ is
coefficient of $z^{-1}$ in series (32a) and $b^{j}_{-2}$ is coefficient of
$z^{-2}$ in series (32b), (32c) or (32d) correspondingly). Consider each of
the Eqs. (23)–(25) separately.

  1. 1.

For radial Eq. (23) there is a quadratic equation for powers
$\zeta(\zeta-1)-\frac{3}{2}\zeta-\frac{3}{2}=0$. They are
$\zeta_{1}=3,\zeta_{2}=-\frac{1}{2}$. This allows us to substitute the series
(33) with the leading powers into Eq. (23), equate the coefficients at all
powers of $\rho$ to zero. We have expressed all unknown coefficients in terms
of $c^{r}_{0}$ and $d^{r}_{0}$, which are assumed to be equal to one,
therefore the solutions to the Eq. (23) are

$\xi^{r}_{1}(\rho)=\rho^{3}+\frac{\rho^{4}}{3}+\frac{\left(8l^{2}-21\right)\rho^{5}}{66l^{2}}-\frac{\left(132E^{2}-40l^{2}+105\right)\rho^{6}}{858l^{2}}+O\left(\rho^{7}\right).$
(34a)
$\xi^{r}_{2}(\rho)=\frac{1}{\sqrt{\rho}}-\frac{\rho^{\frac{1}{2}}}{4}-\frac{\rho^{\frac{3}{2}}}{32}+\frac{\left(32E^{2}-l^{2}\right)\rho^{\frac{5}{2}}}{128l^{2}}+O\left(\rho^{\frac{7}{2}}\right),$
(34b)

And the general solution to Eq. (23) is a linear combination

$\xi^{r}(\rho)=A_{r}\xi^{r}_{1}(\rho)+B_{r}\xi^{r}_{2}(\rho),$ (35)

where $A_{r}$ and $B_{r}$ are arbitrary constants.

  2. 2.

For azimuthal Eq. (24) there is a quadratic equation for powers
$\zeta(\zeta-1)-\frac{3}{2}\zeta+\frac{3}{2}=0$. They are easy to find
$\zeta_{1}=\frac{3}{2},\zeta_{2}=1$. This allows us to substitute the series
(33) with the leading powers into Eq. (24), equate the coefficients at all
powers of $\rho$ to zero. We have expressed all unknown coefficients in terms
of $c^{\theta}_{0}$ and $d^{\theta}_{0}$, which are assumed to be equal to
one, therefore the solutions to the Eq. (24) are

$\xi^{\theta}_{1}(\rho)=\rho^{\frac{3}{2}}-\frac{\rho^{\frac{5}{2}}}{4}-\frac{\left(5l^{2}+16\right)\rho^{\frac{7}{2}}}{160l^{2}}-\frac{\left(160E^{2}+35l^{2}-368\right)\rho^{\frac{9}{2}}}{4480l^{2}}+\left(\rho^{\frac{11}{2}}\right),$
(36a)
$\xi^{\theta}_{2}(\rho)=\rho-\rho^{2}+\frac{\rho^{4}}{5l^{2}}+\frac{\left(5E^{2}-1\right)\rho^{5}}{70l^{2}}+O\left(\rho^{6}\right).$
(36b)

And the general solution to Eq. (24) is a linear combination

$\xi^{\theta}(\rho)=A_{\theta}\xi^{\theta}_{1}(\rho)+B_{\theta}\xi^{\theta}_{2}(\rho),$
(37)

where $A_{\theta}$ and $B_{\theta}$ are arbitrary constants.

  3. 3.

For polar Eq. (25) there is a quadratic equation for powers
$\zeta(\zeta-1)-\frac{3}{2}\zeta=0$. They are
$\zeta_{1}=\frac{5}{2},\zeta_{2}=0$. This allows us to substitute the series
(33) with the leading powers into Eq. (25), equate the coefficients at all
powers of $\rho$ to zero. We have expressed all unknown coefficients in terms
of $c^{\phi}_{0}$ and $d^{\phi}_{0}$, which are assumed to be equal to one,
therefore the solutions to the Eq. (25) are

$\xi^{\phi}_{1}(\rho)=\rho^{\frac{5}{2}}+\frac{5\rho^{\frac{7}{2}}}{28}+\frac{\left(5l^{2}-32\right)\rho^{\frac{9}{2}}}{96l^{2}}-\frac{5\left(224E^{2}-35l^{2}+160\right)\rho^{\frac{11}{2}}}{9856l^{2}}+O\left(\rho^{\frac{11}{2}}\right),$
(38a)
$\xi^{\phi}_{2}(\rho)=1+\frac{\rho^{2}}{2l^{2}}-\frac{\rho^{4}}{8l^{4}}-\frac{E^{2}\rho^{5}}{25l^{4}}+\left(\rho^{6}\right).$
(38b)

And the general solution to Eq. (25) is a linear combination

$\xi^{\theta}(\rho)=A_{\phi}\xi^{\phi}_{1}(\rho)+B_{\phi}\xi^{\phi}_{2}(\rho),$
(39)

where $A_{\phi}$ and $B_{\phi}$ are arbitrary constants.

Leaving only the leading members in (35), (37) and (39) we have

$\xi^{r}(\rho)=\frac{B_{r}}{\sqrt{\rho}}+O\left(\rho^{\frac{1}{2}}\right)\;\text{for}\;\rho\to
0,$ (40a)
$\xi^{\theta}(\rho)=B_{\theta}\rho+O\left(\rho^{\frac{3}{2}}\right)\;\text{for}\;\rho\to
0,$ (40b)
$\xi^{\phi}(\rho)=B_{\phi}+O\left(\rho^{2}\right)\;\text{for}\;\rho\to 0.$
(40c)

Comparing expressions (30) and (40) it is seen that angular momentum presence
does not affect the behavior of the radial geodesic deviation vector component
in the vicinity of zero. At small $\rho$ (30a) and (40a) behave in the leading
order as inverse square root $\xi^{r}=\frac{C}{\sqrt{\rho}}$. So near physical
singularity the test free-falling body experiences tidal stretching along the
radial direction regardless of angular momentum presence. Note that Eq. (40a)
corresponds to numerical solutions in Fig. 3 for $l=3,10,15$.

However, the behavior of the angular components changed with the appearance of
the angular momentum $l$. The azimuthal component (30b) of geodesic deviation
vector without $l$ behaves as $\xi^{\theta}=C\sqrt{\rho}$, but with $l\neq 0$
it behaves linearly $\xi^{\theta}=C\rho$ (40b) at small $\rho$. Note that Eq.
(40b) corresponds to numerical solutions in Fig. 4 for $l=3,10,15$.

Also the polar component of geodesic deviation vector changes its dependence
on the radial variable $\rho$ when angular momentum appears. Without angular
momentum the $\phi$-component behaves like $\theta$-component (30c) as
$\xi^{\phi}=C\sqrt{\rho}$. But when there is angular momentum, the polar
component near $\rho=0$ becomes constant (40c). Note that Eq. (40c)
corresponds to numerical solutions in Fig. 5 for $l=3,10,15$.

####  4.3.2 Infinitely far from BH

To determine the asymptotic behavior of equations solutions (23)–(25) at large
$\rho$ we need to replace the variable $\rho=\frac{1}{s}$ and consider the
solution in the vicinity of the point $s=0$. New form of Eqs. (23)–(25) is

${\xi^{r}}^{\prime\prime}+W(s){\xi^{r}}^{\prime}+G_{r}(s)\xi^{r}=0,$ (41)
${\xi^{\theta}}^{\prime\prime}+W(s){\xi^{\theta}}^{\prime}+G_{\theta}(s)\xi^{\theta}=0,$
(42)
${\xi^{\phi}}^{\prime\prime}+W(s){\xi^{\phi}}^{\prime}+G_{\phi}(s)\xi^{\phi}=0,$
(43)

where

$W(s)=\frac{2}{s}+\frac{1-l^{2}s+3l^{2}s^{2}}{(E^{2}-1)+2s-l^{2}s^{2}+2l^{2}s^{3}},$
(44a)
$G_{r}(s)=-\frac{2+3l^{2}s^{2}}{(E^{2}-1)s+2s^{2}-l^{2}s^{3}+2l^{2}s^{4}},$
(44b)
$G_{\theta}(s)=\frac{1+3l^{2}s^{2}}{(E^{2}-1)s+2s^{2}-l^{2}s^{3}+2l^{2}s^{4}},$
(44c) $G_{\phi}(s)=\frac{1}{(E^{2}-1)s+2s^{2}-l^{2}s^{3}+2l^{2}s^{4}}.$ (44d)

Now we solve these equations in the vicinity of the point $s=0$, which
corresponds to $\rho=\infty$. Therefore, we expand the functions (44) in power
series for $s\to 0$. These series are

$W(s)=\frac{2}{s}+\frac{1}{E^{2}-1}-\left[\frac{l^{2}(E^{2}-1)+2}{(E^{2}-1)^{2}}\right]s+O(s^{2}),$
(45a)
$G_{r}(s)=-\frac{2}{(E^{2}-1)s}+\frac{4}{(E^{2}-1)^{2}}-\left[\frac{3l^{2}(E^{2}-1)^{2}+2l^{2}(E^{2}-1)+8}{(E^{2}-1)^{3}}\right]s+O\left(s^{2}\right),$
(45b)
$G_{\theta}(s)=\frac{1}{(E^{2}-1)s}-\frac{2}{(E^{2}-1)^{2}}+\left[\frac{3l^{2}(E^{2}-1)^{2}+l^{2}(E^{2}-1)+4}{(E^{2}-1)^{3}}\right]s+O\left(s^{2}\right),$
(45c)
$G_{\phi}(s)=\frac{1}{(E^{2}-1)s}-\frac{2}{(E^{2}-1)^{2}}+\left[\frac{l^{2}(E^{2}-1)+4}{(E^{2}-1)^{3}}\right]s+O\left(s^{2}\right).$
(45d)

According to the Fuchs’ theorem Frobenius method is applicable again and the
point $s=0$ is a regular singular point for Eqs. (41)–(43). For all series
(45b)–(45d) coefficient ${b_{-2}=0}$ since the functions (44b)–(44d) in $s=0$
have first order pole and coefficient $a_{-1}=2$ of (45a). Thus, the quadratic
equation for the leading powers has form $\zeta(\zeta-1)+2\zeta=0$. Leading
powers are $\zeta_{1}=0$ and $\zeta_{2}=-1$, their difference is equal to an
integer. Therefore, two linearly independent solutions have the form

$\xi^{j}_{1}(s)=\sum_{n=0}^{\infty}c^{j}_{n}s^{n+\zeta_{1}},$ (46a)
$\xi^{j}_{2}(s)=\sum_{n=0}^{\infty}d^{j}_{n}s^{n+\zeta_{2}}+A\xi^{j}_{1}(s)\ln(s),$
(46b)

where $j=r,\theta,\phi$ corresponds to different equations. To find the
coefficients $\\{c^{j}_{n},d^{j}_{n}\\}_{n=1}^{\infty}$ and $A$, one need to
substitute these series into the Eqs. (41)–(43) and after simplifying like
terms equate the coefficients at all powers of $s$ to zero.

  1. 1.

By substituting expression (46a) in Eq. (41), the coefficients
$c^{r}_{j},\;j\geq 1$ can be expressed in terms of $c^{r}_{0}$, equating which
to one we will find the first solution

$\xi^{r}_{1}(s)=1+\frac{s}{E^{2}-1}-\frac{s^{2}}{2\left(E^{2}-1\right)^{2}}+\frac{\left(E^{4}l^{2}-E^{2}l^{2}+2\right)s^{3}}{4\left(E^{2}-1\right)^{3}}+O\left(s^{4}\right).$
(47)

To find the second linearly independent solution, substitute expression (46b)
into Eq. (41), assuming $\xi^{r}_{1}(s)$ as a solution. Equating the
coefficients at all powers of $s$ to zero and setting $d^{r}_{0}=1$ and
$d^{r}_{1}=0$ since $d^{r}_{0}$ is free coefficient, and $d^{r}_{1}$ is a
factor at $s^{\zeta_{1}}$ in second solution, we get

$\xi^{r}_{2}(s)=\ln(s)\left[\frac{3}{E^{2}-1}+\frac{3s}{\left(E^{2}-1\right)^{2}}-\frac{3s^{2}}{2(E^{2}-1)^{3}}\right]+\frac{1}{s}-\\\
-\frac{\left(E^{2}l^{2}-l^{2}+18\right)s}{2\left(E^{2}-1\right)^{2}}+\frac{\left(4E^{4}l^{2}-3E^{2}l^{2}-l^{2}+9\right)s^{2}}{4\left(E^{2}-1\right)^{3}}+O\left(s^{3}\right).$
(48)

Therefore, the general solution to Eq. (41) is a linear combination of (47)
and (48)

$\xi^{r}(s)=A_{r}\xi^{r}_{1}(s)+B_{r}\xi^{r}_{2}(s),$ (49)

where $A_{r}$ and $B_{r}$ are arbitrary constants.

  2. 2.

By substituting expression (46a) in Eq. (42), the coefficients
$c^{\theta}_{j},\;j\geq 1$ can be expressed in terms of $c^{\theta}_{0}$,
equating which to one we will find the first solution

$\xi^{\theta}_{1}(s)=1-\frac{s}{2\left(E^{2}-1\right)}+\frac{s^{2}}{2\left(E^{2}-1\right)^{2}}-\frac{\left(2E^{4}l^{2}-3E^{2}l^{2}+l^{2}+5\right)s^{3}}{8\left(E^{2}-1\right)^{3}}+O\left(s^{4}\right).$
(50)

To find the second linearly independent solution, substitute expression (46b)
into Eq. (42), assuming $\xi^{\theta}_{1}(s)$ as a solution. Equating the
coefficients at all powers of $s$ to zero and setting $d^{\theta}_{0}=1$ and
$d^{\theta}_{1}=0$ since $d^{\theta}_{0}$ is free coefficient, and
$d^{\theta}_{1}$ is a factor at $s^{\zeta_{1}}$ in second solution, we get

$\xi^{\theta}_{2}(s)=\frac{1}{s}-\frac{l^{2}s}{2\left(E^{2}-1\right)}+\frac{l^{2}s^{2}}{2\left(E^{2}-1\right)^{2}}-\frac{\left(E^{2}l^{4}-l^{4}+5l^{2}\right)s^{3}}{8\left(E^{2}-1\right)^{3}}+O\left(s^{4}\right).$
(51)

Therefore, the general solution to Eq. (42) is a linear combination of (50)
and (51)

$\xi^{r}(s)=A_{\theta}\xi^{\theta}_{1}(s)+B_{\theta}\xi^{\theta}_{2}(s),$ (52)

where $A_{\theta}$ and $B_{\theta}$ are arbitrary constants.

  3. 3.

By substituting expression (46a) in Eq. (43), the coefficients
$c^{\phi}_{j},\;j\geq 1$ can be expressed in terms of $c^{\phi}_{0}$, equating
which to one we will find the first solution

$\xi^{\phi}_{1}(s)=1-\frac{s}{2\left(E^{2}-1\right)}+\frac{s^{2}}{2\left(E^{2}-1\right)^{2}}-\frac{\left(E^{2}l^{2}-l^{2}+5\right)s^{3}}{8\left(E^{2}-1\right)^{3}}+O\left(s^{4}\right).$
(53)

To find the second linearly independent solution, substitute expression (46b)
into Eq. (43), assuming $\xi^{\phi}_{1}(s)$ as a solution. Equating the
coefficients at all powers of $s$ to zero and setting $d^{\phi}_{0}=1$ and
$d^{\phi}_{1}=0$ since $d^{\phi}_{0}$ is free coefficient, and $d^{r}_{1}$ is
a factor at $s^{\zeta_{1}}$ in second solution, we get

$\xi^{\phi}_{2}(s)=\frac{1}{s}-\frac{l^{2}s}{2\left(E^{2}-1\right)}+\frac{E^{2}l^{2}s^{2}}{2\left(E^{2}-1\right)^{2}}-\frac{\left(E^{2}l^{4}+5E^{2}l^{2}-l^{4}\right)s^{3}}{8\left(E^{2}-1\right)^{3}}+O\left(s^{4}\right).$
(54)

Therefore, the general solution to Eq. (43) is a linear combination of (53)
and (54)

$\xi^{r}(s)=A_{\phi}\xi^{\phi}_{1}(s)+B_{\phi}\xi^{\phi}_{2}(s),$ (55)

where $A_{\phi}$ and $B_{\phi}$ are arbitrary constants.

Leaving only the leading members in (49), (52) and (55) we have

$\xi^{r}(s)=B_{r}\left(\frac{1}{s}+\frac{3\ln(s)}{E^{2}-1}\right)+A_{r}+O\left(s\right),$
(56a) $\xi^{\theta}(s)=\frac{B_{\theta}}{s}+A_{\theta}+O\left(s\right),$ (56b)
$\xi^{\phi}(s)=\frac{B_{\phi}}{s}+A_{\phi}+O\left(s\right).$ (56c)

This means that in terms of the radial variable $\rho=\frac{1}{s}$, all
deviation vector components grow linearly for $\rho\to\infty$

$\xi^{r}(r)=B_{r}\left(\rho-\frac{3\ln(\rho)}{E^{2}-1}\right)+O\left(1\right),$
(57a) $\xi^{\theta}(\rho)=B_{\theta}\rho+O\left(1\right),$ (57b)
$\xi^{\phi}(\rho)=B_{\phi}\rho+O\left(1\right).$ (57c)

Result (57) qualitatively coincides with (31). This is to be expected because
the Schwarzschild metric is asymptotically flat, and the angular momentum $l$
in flat space does not matter because in Minkowski space geodesics are
straight lines and they can deviate from each other only linearly in the
radial variable $\rho=\frac{r}{M}$.

##  5 Conclusion

In this article we consider tidal properties of spherical symmetric
spacetimes. In Sec. 3 it was possible to show an alternative way to
diagonalize the geodesic deviation equation in the case of geodesics with non-
zero angular momentum. We perform specific calculations in the Schwarzschild
metric, which is determined only by the black hole mass. Sec. 4 demonstrates
that the presence of angular momentum does not lead to the appearance of
extremum points for all the tidal force components and leaves them constant in
sign.

Further, in Sec. 4.1 we investigate the solutions of geodesic deviation
equations in Schwarzschild spacetime. First we present the previously known
[8] solutions (27) and (28) for radial geodesics, in order to see how the
angular momentum affects geodesic deviation vector behavior. Eqs. (23)–(25)
for $l\neq 0$ have no solutions in terms of elementary functions, therefore,
we found numerical solutions and presented them at different angular momentum
values in Figs. 3–5 in Sec. 4.2. These figures demonstrate that for all Eqs.
(23)–(25) point $\rho=0$ is special, which is expected because it is physical
singularity. Thus, in Sec. 4.3 we consider solutions near $\rho=0$,
$\rho=\infty$. In a vicinity of zero and infinity we found solutions in the
form of generalized power series using the Frobenius method because Equations
(23)–(25) satisfy Fuchs’ theorem.

From a comparison of (30a) and (40a) we found that angular momentum presence
does not change the behavior of radial geodesic deviation vector component in
the vicinity of zero $\rho=0$. Therefore, regardless of fall trajectory on the
Schwarzschild black hole, the test body will experience stretching, which is
described by the leading order of radial geodesic deviation vector component
$\xi^{r}=\frac{C}{\sqrt{r}}$. Which corresponds to Fig. 3. However, the
behavior of the angular geodesic deviation vector components changed as
expected. So the behavior of the azimuthal geodesic deviation vector component
with angular momentum appearance changed from (30b) $\xi^{\theta}=C\sqrt{r}$
to (40b) $\xi^{\theta}=Cr$. Which corresponds to Fig. 4. The polar geodesic
deviation vector component also changes, the appearance of angular momentum
leads to the fact that square root growth (30c) $\xi^{\phi}=C\sqrt{r}$ becomes
constant (40c). This can be seen in Fig. 5.

At very large distances from the black hole, the angular momentum presence, as
expected, does not change the behavior of geodesic deviation vector
components, which can be seen from the comparison (31) and (57). This is due
to the fact that the Schwarzschild metric is asymptotically flat and at large
distances from the gravitating mass the second order differential Eq. (6) has
a zero right-hand side, which means that its solution will be a linear
function of $\tau$. And far from a black hole proper time $\tau$ is
proportional to radial variable $\rho=\frac{r}{M}$ according to Eq. (3). So
the linearity of solutions of Eq. (6) in vacuum and expressions (31) and (57)
have the same origin.

The study of tidal effects for non-radial geodesics in general relativity can
evolve towards considering more general metrics, for example,
Reissner–Nordström or Kottler metric. Also quite interesting is the question
of geodesic deviation when both a test free falling body and a black hole have
angular momentum. In addition, it is interesting to investigate the geodesic
deviation equation in multidimensional spacetimes and, in particular, in the
metrics of black rings. As we can see, despite the fact that the classical
effects in the general relativity have been studied for a very long time, this
issue still contains many unsolved problems.

## Acknowledgement

We would like to express our gratitude to Yuri Viktorovich Pavlov for
meaningful discussions and useful advice.

## References

  * [1] Schwarzschild K. Berliner Sitzungsbesichte (Phys. Math. Klasse), 189 – 196 3 Feb. (1916).
  * [2] Reissner H. Ann. Physik, 50, 106 – 120 (1916).
  * [3] Nordström G. Proc. Kon. Ned. Akad. Wet., 20, 1238 – 1245 (1918).
  * [4] Kerr R. P. Phys. Rev. Lett., 11, 237 – 238 September 1. (1963).
  * [5] Newman E. T., Couch E., Chindapared K. et al. J. Math. Phys., 6, 918 – 919 (1965).
  * [6] Kottler F. Ann. Physik, 56, 401 – 462 (1918).
  * [7] Kiselev, V. V.: Class. Quant. Grav. 20, 1187 (2003).
  * [8] Fuchs H. Ann. Physik, 495, 231 – 233 (1983).
  * [9] Sharif M., Sadiq S., JETP 153, 232 – 239 (2018).
  * [10] Soon-Tae Hong, et al. Phys. Lett. B 881 135967 (2020).
  * [11] Jing Li, et al. Eur. Phys. J. C 81, 590 (2021).
  * [12] Crispino L.C.B., et al.: Eur. Phys. J. C 76, 168 (2016).
  * [13] Shahzad M. U., Jawad A., Eur. Phys. J. C 77, 372 (2017).
  * [14] Vandeev V.P., Semenova A.N.: Eur. Phys. J. C 81: 610 (2021).
  * [15] Luminet J. P., Marck J .A., Mon. Not. R. Astron. Soc. 212, 57 (1985).
  * [16] Haroldo C. D. Lima Junior, et al.: Eur. Phys. J. Plus 135: 334 (2020).
  * [17] Chicone C., Mashhoon B., Punsly B., Int. J. Mod. Phys. D 13, 945 (2004).
  * [18] Chicone C., Mashhoon B., Annalen Phys. 14, 290 (2005).
  * [19] Bini D., Chicone C., Mashhoon B., Phys. Rev. D 95, 104029 (2017).
  * [20] Misner C. W., Thorne K. S., Wheeler J. A., Gravitation (San Francisco: W.H. Freeman, 1973).
  * [21] B.P. Abbott, et al., LIGO Scientific, Virgo, Phys. Rev. Lett. 116 061102 (2016).
  * [22] B.P. Abbott, et al., LIGO Scientific, VIRGO, Phys. Rev. Lett. 118 221101 (2017), erratum: Phys. Rev. Lett. 121 129901 (2018).
  * [23] S. Chandrasekhar, The Mathematical Theory of Black Holes (Clarendon Press, Oxford, 1983).