# Nonhomogeneous Stochastic Geometry Analysis of Massive LEO Communication
Constellations

Niloofar Okati *  | and Taneli Riihonen *  | The work of N. Okati and T. Riihonen was supported by a Nokia University Donation.N. Okati and T. Riihonen are with Unit of Electrical Engineering, Faculty of Information Technology and Communication Sciences, Tampere University, FI-33720 Tampere, Finland (e-mail<EMAIL_ADDRESS>taneli.riihonen@tuni.fi).

###### Abstract

Providing truly ubiquitous connectivity requires development of low Earth
orbit (LEO) satellite Internet, whose theoretical study is lagging behind
network-specific simulations. In this paper, we derive analytical expressions
for the downlink coverage probability and average data rate of a massive
inclined LEO constellation in terms of total interference power’s Laplace
transform in the presence of fading and shadowing, ergo presenting a
stochastic geometry-based analysis. We assume the desired link to experience
Nakagami-$m$ fading, which serves to represent different fading scenarios by
varying integer $m$, while the interfering channels can follow any fading
model without an effect on analytical tractability. To take into account the
inherent non-uniform distribution of satellites across different latitudes, we
model the LEO network as a nonhomogeneous Poisson point process with its
intensity being a function of satellites’ actual distribution in terms of
constellation size, the altitude of the constellation, and the inclination of
orbital planes. From the numerical results, we observe optimum points for both
the constellation altitude and the number of orthogonal frequency channels;
interestingly, the optimum user’s latitude is greater than the inclination
angle due to the presence of fewer interfering satellites. Overall, the
presented study facilitates general stochastic evaluation and planning of
satellite Internet constellations without specific orbital simulations or
tracking data on satellites’ exact positions in space.

###### Index Terms:

Massive communication satellite networks, Low Earth orbit (LEO) Internet
constellations, interference, coverage probability, average achievable data
rate, stochastic geometry, Poisson point process.

##  I Introduction

Recent advances towards 6th generation (6G) wireless networks require
progression and development of non-terrestrial networks to provide seamless
connections with high transmission capacity[1, 2, 3, 4]. Among non-terrestrial
networks, low Earth orbit (LEO) satellite Internet constellations have gained
increasing popularity as they provide global connectivity for unserved or
underserved regions, where the deployment of terrestrial networks is not
feasible or economically reasonable[5, 6]. Deploying thousands of satellites
will ensure that every single person or appliance on Earth could be connected
and no location is left in outage.

While the performance of many LEO constellations (e.g., Starlink, OneWeb,
Kuiper) has been evaluated through network-specific simulations to put the
commercial plans forward, a general scientific understanding of their
performance is limited in the open literature. Conventional simulation-based
studies are restricted to few number of satellites with deterministic
locations which is not capable of evaluating the general performance of a
massive satellite network consisting of thousands of satellites. Moreover, in
most of the literature, the coverage regions are assumed to have fixed
circular footprints, while selecting smaller inclination angles and
simultaneous consolidated operation of several LEO networks render a not-so-
regular Voronoi tessellation.

In this paper, downlink coverage probability and average data rate of inclined
LEO constellations are analyzed under general shadowing and fading propagation
models. The satellites’ positions are assumed to be distributed as a
nonhomogeneous Poisson point process (NPPP), which models the satellites
distribution across varying latitudes precisely by setting the intensity
function to be the actual density of satellites in an actual constellation.

###  I-A Related Works

The literature around LEO networks is mostly limited to deterministic and
simulation-based analyses. In [7], the performance of two different LEO
constellations was simulated assuming specific constellation sizes. The
probability of average call drop and the distribution of the number of
handoffs were studied for the Iridium constellation in [8]. A deterministic
model to characterize the visibility time of one LEO satellite was presented
in [9]. Since the model in [9] is not valid for any arbitrarily located user,
authors in [10] contributed statistical analysis of coverage time in a mobile
LEO constellation. In [11], a LEO-based Internet-of-Things architecture was
presented so as to supply network access for devices distributed in remote
areas.

Stochastic geometry is an area of mathematics, which deals with the study of
random objects on Euclidean space. In the area of telecommunication,
stochastic geometry has been extensively utilized to model, evaluate, and
develop the wireless communication networks with irrigular topologies [12, 13,
14], especially for two-dimensional (planar) terrestrial networks [15, 12, 14,
16, 17, 18, 19, 20, 13]. Various studies in stochastic geometry modeling of
multi-tier and cognitive networks were reviewed in [15]. Observations in [16]
have shown that the Poisson point process (PPP) and a regular grid model
provide lower and upper bounds on the network performance metrics,
respectively, with the same deviation from the actual network performance. The
research in [17, 18], being an extension to [16], modeled a multi-tier network
considering the limitations for the achievable quality-of-service. The
coverage of uplink was studied in [19] assuming base stations and devices are
distributed as independent PPPs.

The application of stochastic geometry to three-dimensional wireless networks
has gained remarkable attention in the literature [21, 22]. In [21], a PPP
model was applied to model and analyze the coverage in three-dimensional
cellular networks. Since the PPP, despite providing tractable analysis, is not
accurate when applied to networks with limited nodes in a finite area [23], a
binomial point process (BPP) can be utilized instead to capture the
characteristics of such networks [24, 25]. A finite network of unmanned aerial
vehicles was modeled as a BPP in [22, 26]. In [24, 27], a planar network with
an arbitrary shape was studied assuming the transmitter is positioned at a
fixed distance. The results were then generalized in [20], by using two
protocols for selecting the transmitter.

On the literature around the LEO satellite networks, the analysis is limited
to few number of satellites with known locations and/or coverage spots. In
[28], with tools from stochastic geometry, authors have developed a method to
characterize the magnitude of Doppler shifts in a LEO network. Resource
control of a satellite–terrestrial network was investigated in [29], in order
to minimize the outage probability and maximize the data rate. Focusing on
only a single spotbeam, the hybrid satellite–terrestrial network supporting 5G
infrastructure has been presented in [30, 31]. In [32], the outage probability
of a satellite-based Internet-of-Things, in which LEO satellites relay uplink
data to ground stations, is derived in closed form by assuming a low number of
satellites at deterministic locations.

Recently, more research on LEO networks using stochastic geometry has started
emerging. In our study[33], generic performance of satellite networking has
been formulated assuming uniform distribution for satellites, and without
considering any explicit model of orbits and constellation geometry. However,
in practical constellations, the satellites are not evenly distributed across
different latitudes [34], i.e., as the user gets farther from the equator
towards the poles, more satellites are visible to it. Therefore, the
constellation densities are typically not uniform in practice which results in
an inconsistency between the performance of actual and uniformly distributed
constellations. In [33], we compensated for this mismatch by numerically
computing a parameter named as effective number of satellites, while in [34],
we derived a mathematical expression for it based on constellation geometry.

The work in [35] characterizes the distance distribution in two different
communication links in a LEO satellite network: link between a user on Earth
and the nearest satellite to it and the link between a satellite and its
nearest neighboring satellite. Unlike in [33], the satellites are assumed to
be placed at different altitudes, i.e., on multiple orbital shells. Stochastic
geometry and the results from [35] were then utilized in [36] to obtain the
downlink probability of coverage for a LEO network, where satellite gateways
act as relays between the satellites and users on Earth. An uplink
communication scenario was characterized by considering interfering
terrestrial transmitters in [37]. In [38], the coverage and rate of a noise-
limited interference-free LEO network were analyzed assuming the satellites’
positions are distributed as a nonhomogeneous Poisson point process (NPPP),
which models the actual distribution of satellites along different latitudes
more precisely by selecting proper intensity for it.

###  I-B Contributions and Paper Organization

We model the satellites’ positions in a LEO network as a nonhomogeneous PPP
which facilitates not only using the tools from stochastic geometry, but also
capturing the exact characteristics of the actual constellations, i.e., the
uneven distribution of satellites across different latitudes. Unlike in[33,
34, 36], by selecting the intensity of NPPP to fit the actual distribution of
satellites on an orbital shell, there is no mismatch between the performance
of theoretical stochastic constellations and actual deterministic LEO
networks. We derive the intensity of NPPP in closed form in terms of the
constellation parameters: the total number of satellites, altitude of the
constellation, latitude of the satellites, and the inclination of the orbits.

As the main contributions, we utilize stochastic geometry to formulate the
coverage and average achievable rate of a user served by a LEO constellation
in terms of the derivative of Laplace transform of interference power.111Thus,
the present study, unlike the preliminary results presented in [38] that are
limited to the special case of scheduling an orthogonal channel for every
satellite, includes the cumulative interference from all other satellites that
are visible to the user and share the same channel. Our derivations do not
rely on exact location of every single satellite and are applicable for
performance analysis of any given constellation as long as the constellation
parameters are known. Modeling the satellites’ locality as a NPPP, the
analytical expressions obtained from a stochastic constellation geometry can
be particularly used to analyze the actual deterministic constellations.

In propagation modelling, unlike most of literature on this topic, we take
into account the effect of shadowing caused by the presence of the obstacles
surrounding the user. To retain analytical tractability and still cover
different fading scenarios, we assume Nakagami-$m$ fading with integer $m$ as
well as shadowing with any desired distribution for the propagation model of
desired links.222Varying the value of $m$, we are able to control the multi-
path fading severity. For instance, $m=1$ corresponds to Rayleigh fading
environment while $m\to\infty$ represents non-fading channels. For interfering
signals, any arbitrarily distributed fading and shadowing can be considered,
since the analytical tractability is unaffected.

Finally, we evaluate two critical performance metrics, i.e., coverage and data
rate, in terms of several key design parameters, such as altitude of the
constellation and the number of frequency channels. From the numerical
results, we are able to observe optimal points for these parameters for some
specific network setup. Counter-intuitively, the user which resides in higher
latitudes, away from the constellation borders, has the best performance due
to existence of fewer interfering satellites in that region.

The remainder of this paper is organized as follows. Section II describes the
system model and the mathematical preliminaries for modeling a LEO network as
a NPPP. The main outcome of this study, which is the derivation of analytical
expressions for downlink coverage probability and average achievable rate of a
terrestrial user, is presented in Section III, which involves also the
analysis of the Laplace transform of interference power. We provide numerical
results in Section IV for the verification of our derivations and studying the
effect of key system parameters such as the size of the constellation and its
altitude as well as the channel parameters on the network performance.
Finally, we conclude the paper in Section V.

##  II System Model

In this section, first, we present the characteristics and geometries of
actual low Earth orbit satellite constellations. Next, we will introduce the
mathematical preliminaries for modeling the actual network as a stochastic
point process.

###  II-A Actual Inclined Constellations

As shown in Fig. 1, we consider a LEO communication satellite constellation
consisting of $N$ satellites launched uniformly on circular orbits with
inclination angle, $\iota$, and altitude that is denoted by $r_{\mathrm{min}}$
— the subscript indicates the minimum possible distance between a satellite
and a ground user (as measured at the zenith). Satellites’ spherical
coordinates in terms of their latitude and longitude are denoted by
$(\phi_{\mathrm{s}},\lambda_{\mathrm{s}})$.

A user terminal is located on any specific latitude, denoted by
$\phi_{\mathrm{u}}$, on the surface of Earth that is approximated as a perfect
sphere with radius $r_{\oplus}\approx 6371$ km. Satellites rising above the
horizon at an angle of $\theta_{\mathrm{s}}\geq\theta_{\mathrm{min}}$ are the
only ones capable of transmitting signals to the users. As such,
$r_{\mathrm{max}}$ refers to the maximum distance at which a satellite and a
user are able to communicate (and it occurs when
$\theta_{\mathrm{s}}=\theta_{\mathrm{min}}$), and

$\displaystyle\frac{r_{\mathrm{max}}}{r_{\oplus}}=\sqrt{\frac{r_{\mathrm{min}}}{r_{\oplus}}\left(\frac{r_{\mathrm{min}}}{r_{\oplus}}+2\right)+\sin^{2}(\theta_{\mathrm{min}})}-\sin(\theta_{\mathrm{min}}).$
(1)

Figure 1: A constellation in an example case of $N=400$ satellites flying on
$\iota=53^{\circ}$ inclined orbits. The borders of two spherical caps above
the user are shown: the outer one covers all visible satellites to the user
while the inner one is empty of satellites and the serving satellite is
located on its border.

In this paper, the serving satellite is the one with the shortest distance to
the user. Assuming that $K$, with $K\leq N$, orthogonal frequency channels are
available for the network, we will distribute $N/K$ satellites randomly among
the channels, which potentially causes co-channel interference to the user.
All satellites on the same channel that are elevated above the horizon to an
angle of $\theta_{\mathrm{s}}\geq\theta_{\mathrm{min}}$ cause interference to
reception of the user.

The variables $R_{0}$ and $R_{n}$, $n=1,2,\ldots,N$, represent the distances
from the user to the serving satellite and the other interfering satellites,
respectively, while $G_{0}$ and $G_{n}$ denote the corresponding channel gains
to model fading. Shadowing effect is modeled by random variables
$\mathcal{X}_{n}$, $n=1,2,\ldots,N$, correspondingly. Obviously,
$G_{n}=\mathcal{X}_{n}=0$ if $R_{n}>r_{\mathrm{max}}$ for some $n=1,\ldots,N$.
To simplify notation, when $N_{\mathrm{I}}>0$, we let indices
$n=1,2,\ldots,N_{\mathrm{I}}$ correspond to those $N_{\mathrm{I}}\leq N/K-1$
satellites with $\theta_{\mathrm{s}}\geq\theta_{\mathrm{min}}$ that cause co-
channel interference. To take into account the effect of directed
transmission, we assign serving and interfering satellites to different power
levels that are denoted by $p_{\mathrm{0}}$ and $p_{n}$, $n=1,2,\ldots,N$,
respectively, such that $p_{n}\leq p_{\mathrm{0}}$. Table I summarizes the
notation used in this paper.

According to the described model, the signal-to-interference-plus-noise ratio
(SINR) of the link is given by

$\displaystyle\mathrm{SINR}=\frac{p_{\mathrm{0}}G_{\mathrm{0}}\mathcal{X}_{\mathrm{0}}R_{\mathrm{0}}^{-\alpha}}{I+\sigma^{2}},$
(2)

where the constant $\sigma^{2}$ is the additive noise power, the parameter
$\alpha$ is a path loss exponent, and

$\displaystyle
I\triangleq\sum_{n=1}^{N_{\mathrm{I}}}p_{n}G_{n}\mathcal{X}_{n}R_{n}^{-\alpha}$
(3)

is the cumulative interference power from all $N_{\mathrm{I}}$ other
satellites above the user’s horizon that share the same frequency channel with
the serving satellite. The distance from the user to its nearest satellite is

$\displaystyle R_{\mathrm{0}}=\min_{n=1,2,\ldots,N_{\mathrm{vis}}}R_{n},$ (4)

where $N_{\mathrm{vis}}$ is a variable representing the number of visible
satellites to the user (cf. the outer cap in Fig. 1).

TABLE I: Summary of mathematical notation Notation | Description
---|---
$r_{\oplus};r_{\mathrm{min}};r_{\mathrm{max}}$ | Earth radius ($6371$ km); Constellation altitude; Maximum possible distance between a user and a visible satellite
$R_{\mathrm{0}}$; $R_{n}$ | Serving distance; Distance to the $n^{\mathrm{th}}$ satellite
$N;N_{\mathrm{I}};K$ | Constellation size; The number of interfering satellites; The number of frequency channels
$G_{\mathrm{0}};G_{n}$ | Channel fading gain of the serving link; Channel fading gain of the $n^{\mathrm{th}}$ link
$\mathcal{X}_{\mathrm{0}};\mathcal{X}_{n}$ | Shadowing component of the serving link; Shadowing component of the $n^{\mathrm{th}}$ link
$\alpha$ | Path loss exponent
$p_{\mathrm{0}};p_{n}$ | Transmission power of the serving satellite; Transmission power of the interfering satellites
$\sigma^{2}$ | Additive noise power
$T$ | SINR threshold for coverage probability
$P_{\mathrm{c}};\bar{C}$ | Coverage probability; Average achievable rate

###  II-B Nonhomogeneous PPP Model

In the constellation described earlier, the satellites appear unevenly along
the lines of latitudes, which means, e.g., that there are more visible
satellites for a user located close to inclination limits than for one on
equatorial region. In order to model the latitude-dependent distribution of
satellites, we assume that the satellites are distributed according to a
nonhomogeneous PPP, $\xi$, on a spherical surface with radius
$r_{\oplus}+r_{\mathrm{min}}$. The NPPP is characterized with its intensity,
$\delta(\phi_{\mathrm{s}},\lambda_{\mathrm{s}})$, which varies according to
the satellites’ latitude (and/or longitude).

By the definition of a NPPP, the number of points in some bounded region
${\cal A}$ of the orbital shell is a Poisson-distributed random variable
denoted by $\mathcal{N}$. Thereby, the probability to have $n$ satellites in
${\cal A}$ is given by

$\displaystyle P_{n}\left({\cal
A}\right)\triangleq{\mathbb{P}}\left(\mathcal{N}=n\right)$ (5)
$\displaystyle=\frac{1}{n!}\left(\iint_{\cal
A}\delta(\phi_{\mathrm{s}},\lambda_{\mathrm{s}})\,(r_{\mathrm{min}}+r_{\oplus})^{2}\cos(\phi_{\mathrm{s}})\,d\phi_{\mathrm{s}}d\lambda_{\mathrm{s}}\right)^{n}$
$\displaystyle\times\exp\left(-\iint_{\cal
A}\delta(\phi_{\mathrm{s}},\lambda_{\mathrm{s}})\,(r_{\mathrm{min}}+r_{\oplus})^{2}\cos(\phi_{\mathrm{s}})\,d\phi_{\mathrm{s}}d\lambda_{\mathrm{s}}\right),$

where $\delta(\phi_{\mathrm{s}},\lambda_{\mathrm{s}})$ is the intensity
function of nonhomogeneous PPP at latitude $\phi_{\mathrm{s}}$ and longitude
$\lambda_{\mathrm{s}}$. Based on the given system model, ${\cal A}$ is the
spherical cap where viewable satellites to the user exist (cf. the outer one
in Fig. 1), with surface area
$\left({\delta\pi\left(r_{\mathrm{max}}^{2}-r_{\mathrm{0}}^{2}\right)}\right)/({1-\frac{r_{\mathrm{min}}}{r_{\oplus}+r_{\mathrm{min}}}})$
(See [33, Appendix A]).

In order to precisely model a LEO network as a NPPP, we first need to
characterize the intensity function based on the actual physical network as
follows. The preliminaries obtained herein will be used shortly towards
contributing expressions for probability of coverage and average achievable
rate.

######  Lemma 1.

When satellites are distributed uniformly on low circular orbits with the
altitude, $r_{\mathrm{min}}$, and the inclination angle, $\iota$, the
intensity function of the nonhomogeneous PPP is a function of latitude,
$\phi_{\mathrm{s}}$, only and given by

$\delta(\phi_{\mathrm{s}})=\frac{N}{\sqrt{2}\pi^{2}(r_{\mathrm{min}}+r_{\oplus})^{2}}\cdot\frac{1}{\sqrt{\cos(2\phi_{\mathrm{s}})-\cos(2\iota)}},$
(6)

and we can denote
$\delta(\phi_{\mathrm{s}},\lambda_{\mathrm{s}})=\delta(\phi_{\mathrm{s}})$
since it does not depend on $\lambda_{\mathrm{s}}$, for
$\phi_{\mathrm{s}}\in[-\iota,\iota]$.

###### Proof.

The intensity function is equivalent to the actual density of the satellites
on an orbital shell element created by spanning the azimuthal angle from
$0^{\circ}$ to $360^{\circ}$ on the orbital shell at latitude
$\phi_{\mathrm{s}}$, that is calculated as

$\displaystyle\delta(\phi_{\mathrm{s}})=\frac{Nf_{\Phi_{\mathrm{s}}}(\phi_{\mathrm{s}})\,d\phi_{\mathrm{s}}}{2\pi(r_{\mathrm{min}}+r_{\oplus})^{2}\cos(\phi_{\mathrm{s}})\,d\phi_{\mathrm{s}}},$
(7)

which is the ratio of the number of satellites resided on the surface element
to the element’s area. Substituting the probability density
$f_{\Phi_{\mathrm{s}}}(\phi_{\mathrm{s}})$ of random latitude
$\Phi_{\mathrm{s}}$ [34, Lemma 2] completes the proof. ∎

If the intensity of satellites is simplistically presumed to be uniform all
over the orbital shell, it can be written as follows.

######  Lemma 2.

When satellites are uniformly distributed on a sphere, the point process turns
into a homogeneous Poisson point process with a constant intensity given by

$\delta=\frac{N}{4\pi(r_{\mathrm{min}}+r_{\oplus})^{2}},$ (8)

which does not depend on latitudinal/longitudinal parameters.

Thus, for the special case when satellites are distributed uniformly on the
orbital shell, by substitution from Lemma 2, the probability given in (5) can
be expressed in closed form as

$\displaystyle P_{n}\left({\cal
A}\right)=\frac{1}{n!}\left(\frac{N\left(r_{\mathrm{max}}^{2}-r_{\mathrm{min}}^{2}\right)}{4r_{\oplus}\left(r_{\oplus}+r_{\mathrm{min}}\right)}\right)^{n}\exp\left(-\frac{N\left(r_{\mathrm{max}}^{2}-r_{\mathrm{min}}^{2}\right)}{4r_{\oplus}\left(r_{\oplus}+r_{\mathrm{min}}\right)}\right).$
(9)

##  III Performance Analysis

In this section, we focus on the performance analysis of a LEO satellite
network in terms of coverage probability and data rate of a user in an
arbitrary location on Earth. We utilize stochastic geometry in order to
formulate coverage probability and rate as a function of the network and the
propagation parameters. Two main components of our analytical derivations are
the distribution of the distance to the nearest satellite and the Laplace
function of interference which will be presented throughout this section.

###  III-A The Distance to The Nearest Satellite

In this paper, the serving satellite is assumed to be the nearest one to the
user. Therefore, an important parameter for coverage and rate analysis is the
probability density function (PDF) of the distance to the nearest satellite,
$R_{\mathrm{0}}$, which is given as follows.

######  Lemma 3.

The PDF of the serving distance $R_{\mathrm{0}}$, when the satellites are
distributed according to a nonhomogeneous PPP with a latitude-dependent
intensity, $\delta(\phi_{\mathrm{s}})$, is

$\displaystyle f_{R_{\mathrm{0}}}\left(r_{\mathrm{0}}\right)$
$\displaystyle=2r_{\mathrm{0}}\left(\frac{r_{\mathrm{min}}}{r_{\oplus}}+1\right)\exp(-\gamma(r_{\mathrm{0}}))\int_{\max(\phi_{\mathrm{u}}-\theta_{\mathrm{s}},-\iota)}^{\min(\phi_{\mathrm{u}}+\theta_{\mathrm{s}},\iota)}\delta(\phi_{\mathrm{s}})$
$\displaystyle\hskip
85.35826pt\times\frac{\cos(\phi_{\mathrm{s}})}{\sqrt{\cos^{2}(\phi_{\mathrm{s}}-\phi_{\mathrm{u}})-\cos^{2}(\theta_{\mathrm{s}})}}d\phi_{\mathrm{s}},$
(10)

where

$\displaystyle\gamma(r_{\mathrm{0}})=2(r_{\mathrm{min}}+r_{\oplus})^{2}$
$\displaystyle\times\int_{\max(\phi_{\mathrm{u}}-\theta_{\mathrm{s}},-\iota)}^{\min(\phi_{\mathrm{u}}+\theta_{\mathrm{s}},\iota)}\delta(\phi_{\mathrm{s}})\cos(\phi_{\mathrm{s}})\cos^{-1}\left({\frac{\cos(\theta_{\mathrm{s}})}{\cos(\phi_{\mathrm{s}}-\phi_{\mathrm{u}})}}\right)d\phi_{\mathrm{s}}$
(11)

and $\theta_{\mathrm{s}}$ is the polar angle difference between the serving
satellite and the user which is given by
$\theta_{\mathrm{s}}=\cos^{-1}\left(1-\frac{r_{\mathrm{0}}^{2}-r_{\mathrm{min}}^{2}}{2(r_{\mathrm{min}}+r_{\oplus})r_{\oplus}}\right)$.
Equation (3) is valid for $\theta_{\mathrm{s}}\geq|{\phi_{\mathrm{u}}}|-\iota$
and $r_{\mathrm{0}}\in[r_{\mathrm{min}},2r_{\oplus}+r_{\mathrm{min}}]$ while
$f_{R_{\mathrm{0}}}\left(r_{\mathrm{0}}\right)=0$ otherwise.

###### Proof.

See Appendix -A. ∎

When the density of satellites is presumed to be uniform, i.e., it is not a
function of latitude, the PDF of the serving distance can be obtained in
closed form as follows.

######  Lemma 4.

The PDF of the serving distance $R_{\mathrm{0}}$ when the satellites are
distributed uniformly with constant intensity given in Lemma 2 is

$\displaystyle
f_{R_{\mathrm{0}}}\left(r_{\mathrm{0}}\right)=\frac{Nr_{\mathrm{0}}}{2r_{\oplus}(r_{\mathrm{min}}+r_{\oplus})}\exp\left(-N\left(\frac{r_{\mathrm{0}}^{2}-r_{\mathrm{min}}^{2}}{4(r_{\mathrm{min}}+r_{\oplus})r_{\oplus}}\right)\right)$
(12)

for $r_{\mathrm{0}}\in[r_{\mathrm{min}},2r_{\oplus}+r_{\mathrm{min}}]$ while
$f_{R_{\mathrm{0}}}\left(r_{\mathrm{0}}\right)=0$ otherwise.

###### Proof.

In this proof, the same principles are used as in Lemma 3. However, the
integration from a constant density over the cap sphere will reduce to a
simple expression. Thus,
$F_{R_{\mathrm{0}}}\left(r_{\mathrm{0}}\right)=1-\exp\left(-N\left(\frac{r_{\mathrm{0}}^{2}-r_{\mathrm{min}}^{2}}{4(r_{\mathrm{min}}+r_{\oplus})r_{\oplus}}\right)\right)$.
Calculating the derivative of the CDF w.r.t. $r_{\mathrm{0}}$ completes the
proof. ∎

The PDF of the serving distance is also derived in [33] assuming satellites
are uniformly distributed as a BPP. It is worth noting that the Taylor series
expansion of Lemma 4 and the serving distance distribution given in [33, Lemma
2] are the same for the first two terms. The difference between the serving
distance in a uniformly distributed constellation given in [33] and the
homogeneous Poisson point process in Lemma 4 is insignificant since the
argument of exponential function in (12), i.e.,
$N\left(\frac{r_{\mathrm{0}}^{2}-r_{\mathrm{min}}^{2}}{4(r_{\mathrm{min}}+r_{\oplus})r_{\oplus}}\right)$,
is small.

###  III-B Coverage Probability and Average Data Rate

The coverage probability is the probability that the SINR at the receiver is
higher than the minimum SINR required to successfully transmit the data. The
coverage probability is defined as

$\displaystyle
P_{\mathrm{c}}\left(T\right)\triangleq{\mathbb{P}}\left(\mathrm{SINR}>T\right)={\mathbb{P}}\left(\frac{p_{\mathrm{0}}G_{\mathrm{0}}\mathcal{X}_{\mathrm{0}}R_{\mathrm{0}}^{-\alpha}}{I+\sigma^{2}}>T\right),$
(13)

where $T$ is the SINR threshold.

Using the above definition, we express the coverage probability of a user in
the following theorem.

######  Theorem 1.

The coverage probability for an arbitrarily located user under a Nakagami
fading serving channel while both shape parameter and rate parameter of gamma
distribution333Channel gain, being the square of Nakagami random variable,
follows a gamma distribution. are $m_{0}$, is

$\displaystyle P_{\mathrm{c}}\left(T\right)$
$\displaystyle\triangleq{\mathbb{P}}\left(\mathrm{SINR}>T\right)$
$\displaystyle=\int_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}\int_{0}^{\infty}f_{\mathcal{X}_{\mathrm{0}}}(\mathit{x}_{\mathrm{0}})f_{R_{\mathrm{0}}}\left(r_{\mathrm{0}}\right)\Bigg{[}e^{-s\sigma^{2}}$
$\displaystyle\sum_{k=0}^{m_{0}-1}\frac{\sum_{l=0}^{k}\binom{k}{l}\left(s\sigma^{2}\right)^{l}\left(-s\right)^{k-l}\frac{\partial^{k-l}}{\partial
s^{k-l}}\mathcal{L}_{I}(s)}{k!}\Bigg{]}_{s=\frac{m_{0}Tr_{\mathrm{0}}^{\alpha}}{p_{\mathrm{0}}\mathit{x}_{\mathrm{0}}}}\,d\mathit{x}_{\mathrm{0}}dr_{\mathrm{0}},$
(14)

where the PDF $f_{R_{\mathrm{0}}}\left(r_{\mathrm{0}}\right)$ is given in
Lemma 3 (or Lemma 4), $f_{\mathcal{X}_{\mathrm{0}}}(\mathit{x}_{\mathrm{0}})$
is the PDF of $\mathcal{X}_{\mathrm{0}}$ and $\mathcal{L}_{I}\left(s\right)$
is the Laplace transform of interference power $I$ calculated in the next
section.

###### Proof.

See Appendix -B. ∎

Let us then move on the average achievable data rate (in bit/s/Hz), which is
the ergodic capacity for a fading communication link derived from Shannon-
Hartley theorem normalized to unit bandwidth. The average achievable rate is
defined as

$\displaystyle\bar{C}\triangleq\frac{1}{K}~{}\mathbb{E}\left[\log_{2}\left(1+\mathrm{SINR}\right)\right].$
(15)

Unlike for the coverage probability, frequency reuse affects the average rate
in two opposite directions. One direction is the improvement in SINR due to
the reduction in the number of interfering satellites which use the same
channel. The other direction which results in lower data rate is induced by a
reduction in the availability of the frequency band shared among a group of
satellites.

In the following theorem, we calculate the expression for the average rate of
a user over Nakagami fading serving channel. The interfering channel gains may
follow any arbitrary distribution.

######  Theorem 2.

The average data rate of an arbitrarily located user under a Nakagami fading
serving channel and any fading or shadowing distribution for interfering
channels is given by

$\displaystyle\bar{C}$
$\displaystyle=\frac{1}{K}\int_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}\int_{0}^{\infty}\int_{0}^{\infty}f_{\mathcal{X}_{\mathrm{0}}}(\mathit{x}_{\mathrm{0}})f_{R_{\mathrm{0}}}\left(r_{\mathrm{0}}\right)\Bigg{[}e^{-s\sigma^{2}}$
$\displaystyle\sum_{k=0}^{m_{0}-1}\frac{\sum_{l=0}^{k}\binom{k}{l}\left(s\sigma^{2}\right)^{l}\left(-s\right)^{k-l}\frac{\partial^{k-l}}{\partial
s^{k-l}}\mathcal{L}_{I}(s)}{k!}\Bigg{]}_{s=\frac{m_{0}\left(2^{t}-1\right)r_{\mathrm{0}}^{\alpha}}{p_{\mathrm{0}}\mathit{x}_{\mathrm{0}}}}dtd\mathit{x}_{\mathrm{0}}dr_{\mathrm{0}},$
(16)

where $m_{0}$ is the parameter of Nakagami fading, and $\mathcal{L}_{I}(s)$
will be given in Lemma 5 and its corresponding corollaries, which cover some
special cases.

###### Proof.

See Appendix -C. ∎

###  III-C Interference Analysis

In this subsection, we derive the Laplace function of interference which is a
key element of the performance expressions in Theorems 1 and 2. We, first,
obtain the expression considering a general propagation model which means that
no assumption is made regarding the specific expressions of
$f_{\mathcal{X}_{n}}(x_{n})$ and $f_{G_{n}}(g_{n})$.

######  Lemma 5.

When the server is at distance $r_{\mathrm{0}}\geq r_{\mathrm{min}}$ from the
user and interfering channels experience an arbitrary distributed fading, the
Laplace transform of random variable $I$ is

$\displaystyle\mathcal{L}_{I}(s)$
$\displaystyle\triangleq\mathbb{E}_{I}\left[e^{-sI}\right]=\sum_{n=0}^{\infty}P_{n}\left(\mathcal{A}\left(r_{\mathrm{max}}\right)-\mathcal{A}\left(r_{\mathrm{0}}\right)\right)$
(17)
$\displaystyle\times\Bigg{(}\int_{r_{\mathrm{0}}}^{r_{\mathrm{max}}}\int_{0}^{\infty}\mathcal{L}_{G_{n}}(sp_{n}x_{n}r_{n}^{-\alpha})f_{\mathcal{X}_{n}}(x_{n})f_{R_{n}|R_{\mathrm{0}}}(r_{n}|r_{\mathrm{0}})d\mathit{x}_{\mathrm{0}}dr_{n}\Bigg{)}^{n},$

where

$f_{R_{n}|R_{\mathrm{0}}}(r_{n}|r_{\mathrm{0}})=\frac{{d\gamma(r_{n})}/{dr_{n}}}{\gamma(r_{\mathrm{max}})-\gamma(r_{\mathrm{0}})}$
(18)

is the probability density function of the distance from any visible satellite
to the user conditioned on the serving distance [33, Lemma 3]. The parameter
$\mathcal{A}\left(r_{\mathrm{max}}\right)$ represents the spherical cap where
all satellites that can be viewed by the user exist while
$\mathcal{A}\left(r_{\mathrm{0}}\right)$ is the cap above the user, empty of
satellites and with the serving satellite on its border (base of the cap). The
function $f_{\mathcal{X}_{n}}(x_{n})$ denotes the shadowing PDF for the $n$th
interfering channel.

###### Proof.

See Appendix -D. ∎

In the special cases, where channels experience Nakagami fading without any
shadowing, Lemma 5 can be reduced into the following corollary. The expression
thereof is obtained by substituting the Laplace function of a gamma random
variable, i.e.,
$\mathcal{L}_{G_{n}}(z)=\frac{m_{n}^{m_{n}}}{(m_{n}+z)^{m_{n}}}$, where
$m_{n}$ stands for both shape parameter and rate parameter of gamma
distribution for the $n$th link.

######  Corollary 1.

When the interfering channels experience only Nakagami fading (no shadowing),
the Laplace function of interference can be written as

$\displaystyle\mathcal{L}_{I}(s)$
$\displaystyle=\sum_{n=0}^{\infty}P_{n}\left(\mathcal{A}\left(r_{\mathrm{max}}\right)-\mathcal{A}\left(r_{\mathrm{0}}\right)\right)$
$\displaystyle\times\Bigg{(}\int_{r_{\mathrm{0}}}^{r_{\mathrm{max}}}\frac{m_{n}^{m_{n}}}{(m_{n}+sp_{n}r_{n}^{-\alpha})^{m_{n}}}f_{R_{n}|R_{\mathrm{0}}}(r_{n}|r_{\mathrm{0}})dr_{n}\Bigg{)}^{n},$
(19)

where $\mathcal{A}\left(r_{\mathrm{max}}\right)$ and
$\mathcal{A}\left(r_{\mathrm{0}}\right)$ are the visible cap and the null cap
above the user, respectively. The PDF
$f_{R_{n}|R_{\mathrm{0}}}(r_{n}|r_{\mathrm{0}})$ is given in (18), and $m_{n}$
is the Nakagami fading parameter for $n$th link.

When the intensity of the PPP is presumed to be constant (regardless of the
latitude), the Laplace function can be obtained from the following corollary
by simply substituting $\gamma(\cdot)$ in (18) by the product of the density
in Lemma 2 and the surface area of the spherical cap formed by the distance
between the user and the given interfering satellite.

######  Corollary 2.

The Laplace function of interference when the satellites are distributed
uniformly with constant intensity, and their channels experience Nakagami
fading, is given by

$\displaystyle\mathcal{L}_{I}(s)$
$\displaystyle=\sum_{n=0}^{\infty}\frac{1}{n!}\left(\frac{N\left(r_{\mathrm{max}}^{2}-r_{\mathrm{0}}^{2}\right)}{4r_{\oplus}\left(r_{\oplus}+r_{\mathrm{min}}\right)}\right)^{n}\exp\left(-\frac{N\left(r_{\mathrm{max}}^{2}-r_{\mathrm{0}}^{2}\right)}{4r_{\oplus}\left(r_{\oplus}+r_{\mathrm{min}}\right)}\right)$
$\displaystyle\times\Bigg{(}\int_{r_{\mathrm{0}}}^{r_{\mathrm{max}}}\int_{0}^{\infty}\frac{2r_{n}m_{n}^{m_{n}}f_{\mathcal{X}_{\mathrm{0}}}(\mathit{x}_{\mathrm{0}})}{\left(r_{\mathrm{max}}^{2}-r_{\mathrm{0}}^{2}\right)(m_{n}+sp_{n}r_{n}^{-\alpha})^{m_{n}}}\frac{}{}d\mathit{x}_{\mathrm{0}}dr_{n}\Bigg{)}^{n},$
(20)

where $f_{\mathcal{X}_{n}}(x_{n})$ is the PDF of the shadowing component and
$m_{n}$ is the fading parameter for Nakagami fading.

The following corollary presents the Laplace function of interference in
closed-form, under the assumptions given in Corollary 2 and additionally
excluding shadowing from the propagation model.

######  Corollary 3.

The Laplace function of interference, when the satellites are distributed
uniformly with constant intensity and their channels experience only Nakagami
fading (no shadowing), is given by

$\displaystyle\mathcal{L}_{I}(s)$
$\displaystyle=\sum_{n=0}^{\infty}\frac{1}{n!}\left(\frac{N\left(r_{\mathrm{max}}^{2}-r_{\mathrm{0}}^{2}\right)}{4r_{\oplus}\left(r_{\oplus}+r_{\mathrm{min}}\right)}\right)^{n}\exp\left(-\frac{N\left(r_{\mathrm{max}}^{2}-r_{\mathrm{0}}^{2}\right)}{4r_{\oplus}\left(r_{\oplus}+r_{\mathrm{min}}\right)}\right)$
$\displaystyle\times\frac{1}{\left(r_{\mathrm{max}}^{2}-r_{\mathrm{0}}^{2}\right)}\Bigg{[}r_{\mathrm{max}}^{2}\,{}_{2}F_{1}\left(-\frac{2}{\alpha},m_{n};\frac{\alpha-2}{\alpha};-\,\frac{sp_{n}r_{\mathrm{max}}^{-\alpha}}{m_{n}}\right)^{n}$
$\displaystyle-{r_{\mathrm{0}}^{2}}\,{}_{2}F_{1}\left(-\frac{2}{\alpha},m_{n};\frac{\alpha-2}{\alpha};-\,\frac{sp_{n}r_{\mathrm{0}}^{-\alpha}}{m_{n}}\right)^{n}\Bigg{]},$
(21)

where ${}_{2}F_{1}\left(\cdot,\cdot;\cdot;\cdot\right)$ is the Gauss’s hyper-
geometric function and $m_{n}$ is the fading parameter.

Finally, using the function given in [39, Eq. 9.100] and substitution from
special parameter values, the above can be reduced to elementary functions.
For instance, when $m=1$ and $\alpha=2$, we have

$\displaystyle\mathcal{L}_{I}(s)$
$\displaystyle=\sum_{n=0}^{\infty}\frac{1}{n!}\left(\frac{N\left(r_{\mathrm{max}}^{2}-r_{\mathrm{0}}^{2}\right)}{4r_{\oplus}\left(r_{\oplus}+r_{\mathrm{min}}\right)}\right)^{n}\exp\left(-\frac{N\left(r_{\mathrm{max}}^{2}-r_{\mathrm{0}}^{2}\right)}{4r_{\oplus}\left(r_{\oplus}+r_{\mathrm{min}}\right)}\right)$
$\displaystyle\hskip
28.45274pt\times\left(1+\frac{sp_{n}}{\left(r_{\mathrm{max}}^{2}-r_{\mathrm{0}}^{2}\right)}\ln\left(\frac{k+r_{\mathrm{0}}^{2}}{k+r_{\mathrm{max}}^{2}}\right)\right).$
(22)

To perform frequency reuse, we assign each satellite randomly and
independently to a particular frequency channel. Therefore, the satellites
assigned to each of the orthogonal frequency channels form a thinned version
of the original PPP with intensity $\delta(\phi_{\mathrm{s}})/K$. Since
thinning preserves the Poisson point process according the thinning theorem of
PPP [13], we can take into account the effect of frequency reuse by
substitution $\delta(\phi_{\mathrm{s}})\to\delta(\phi_{\mathrm{s}})/K$ in
Laplace function of interference (in Lemma 5 or the corresponding
corollaries). Since the same frequency channel is used by the user and its
nearest satellite, the frequency reuse has no effect on the original value of
intensity that is used to obtain the PDF of the distance from the user to the
server in Lemma 3.

##  IV Numerical Results

In this section, we provide numerical results to study the effect of different
network parameters on coverage probability and average data rate using the
analytical expressions obtained in Section III. The performance of the network
is evaluated in terms of satellite altitude and the number of orthogonal
frequency channels, which provides important guidelines into the satellite
network design. Furthermore, our analytical derivations are all verified
through Monte Carlo simulations.

In the propagation model, we consider the large-scale attenuation with path
loss exponent $\alpha=2$, the small-scale Nakagami-$m$ fading with integer
$m\in\\{1,2,3\\}$, and lognormal shadowing. The lognormal shadowing is
represented as $\mathcal{X}_{\mathrm{0}}=10^{X_{0}/10}$ such that $X_{0}$ has
a normal distribution with mean $\mu_{\mathrm{s}}=0$ and standard deviation
$\sigma_{\mathrm{s}}=9$ decibels. Thus, the PDF of lognormal shadowing is

$\displaystyle
f_{\mathcal{X}_{\mathrm{0}}}(\mathit{x}_{\mathrm{0}})=\frac{10}{\ln(10)\sqrt{2\pi}\sigma_{\mathrm{s}}\mathit{x}_{\mathrm{0}}}\exp\left({-\frac{1}{2}\left(\frac{10\log_{10}(\mathit{x}_{\mathrm{0}})-\mu_{\mathrm{s}}}{\sigma_{\mathrm{s}}}\right)^{2}}\right).$
(23)

The number of orthogonal channels is set to $K=10$ in all the numerical
results unless otherwise stated. For the reference simulations, satellites are
placed uniformly on orbits centered at Earth’s center with radius
$r_{\oplus}+r_{\mathrm{min}}$.

Figure 2: Theorem 1 verification with simulations when
$\phi_{\mathrm{u}}=25^{\circ}$, $\iota=53^{\circ}$, $m\in\\{1,2,3\\}$, and
$\theta_{\mathrm{min}}=10^{\circ}$.

Figure 2 verifies our derivations given in Theorem 1 for $53^{\circ}$ inclined
orbits and a user located at $25^{\circ}$ latitude. The total number of
satellites and constellation altitude are chosen to be 2000 and 500 km,
respectively. As shown in the figure, the markers that depict the Monte Carlo
simulation results are completely matched with the lines that represent our
theoretical expressions. As the chance of the user being in outage increases,
shadowing affects the coverage probability more positively, the reason being
that shadowing randomness increases the chance of a user with poor average
$\mathrm{SINR}$ to be in coverage. Larger values of $m$ correspond to higher
elevation angles and, consequently, less multi-path distortion, which result
in slightly better coverage, but shadowing masks the effect of fading at
large.

Figure 3: Theoretical coverage probability on different users’ latitudes,
$T=5$ dB, $r_{\mathrm{min}}=500$ km, $\iota=53^{\circ}$, $m=2$, and
$\theta_{\mathrm{min}}=10^{\circ}$. The expression given in Theorem 1 is used
to plot this figure. Figure 4: Effect of altitude on coverage probability when
$T=5$ dB, $\phi_{\mathrm{u}}=\\{0^{\circ},25^{\circ},65^{\circ}\\}$,
$\iota=53^{\circ}$, $m=2$, and $\theta_{\mathrm{min}}=10^{\circ}$. The
expression given in Theorem 1 is used to plot this figure.

The effect of user’s latitude on coverage probability is depicted in Fig. 3.
The solid lines show the coverage probability when intensity of PPP is the
function given in Lemma 1 which forms a nonhomogeneous PPP with the intensity
increasing when the user moves from the equator to the polar regions towards
the inclination limits. As the user moves to higher latitudes, the performance
becomes more unreliable due to increase in the density of satellites that use
the same frequency channel as the user. As the spherical cap above the user
(where all visible satellites to the user reside) leaves the constellation
borders, the coverage probability starts increasing due to the reduction in
the number of visible interferers. The coverage reaches its maximum at about
$66^{\circ}$ where the only (if any) visible satellite to the user is the
serving one, i.e., the performance becomes noise-limited. For latitudes larger
than $66^{\circ}$, the coverage converges to zero quickly, since there are no
satellites above horizon to serve the user. When the intensity of satellites
is selected according to Lemma 2, the coverage probability remains constant
(dashed lines) all over the Earth for any latitude.

Figure 4 illustrates the probability of coverage at different latitudes of the
user. For all propagation environments, the coverage probability increases to
reach its maximum value as the altitude increases which is then followed by a
decline due to the rise in the number of visible interfering satellites. The
optimum altitude for $\phi_{\mathrm{u}}=0^{\circ}$ is slightly larger than
$\phi_{\mathrm{u}}=25^{\circ}$, the reason being that the intensity of
satellites, and consequently the number of interferers, is higher at upper
latitudes. When the user’s latitude is set to $\phi_{\mathrm{u}}=65^{\circ}$,
which means that the user is located out of the constellation borders
($\phi_{\mathrm{u}}>\iota=53^{\circ}$), a larger altitude is crucial for the
constellation so that the user can be served by a satellite within its visible
range. As a result, for altitudes lower than about $400$ km, no visible
satellite is available to serve the user which results in zero coverage
probability.

Figure 5: Theorem 2 verification with simulations when
$\phi_{\mathrm{u}}=25^{\circ}$, $\iota=53^{\circ}$, $m\in\\{1,2,3\\}$, and
$\theta_{\mathrm{min}}=10^{\circ}$.

Figure 5 verifies the derivations for average data rate of a user at the
latitude of $25^{\circ}$. As shown in the figure, the simulation results are
perfectly in line with theoretical expressions in Theorem 2. The disparate
behaviour of the curves is caused by the two opposite effects of frequency
reuse on the average achievable rate. As the number of frequency bands
increases, the total number of interfering satellites on the same frequency
band declines which results in an increase in the data rate. On the other
hand, by increasing the number of frequency channels, the bandwidth shared
among a group of satellites is reduced. An increase in the plot is observed at
first as a result of the decrease in interference received power, followed by
a drop which is due to comprising only $\frac{1}{K}$ of frequency band.

Figure 6, as a counterpart for Fig. 3, illustrates the variation of data rate
over different latitudes when the intensity of Poisson point process is chosen
to be according to Lemma 1 or, for comparison, Lemma 2. With intensity being
as in Lemma 1, data rate varies over the different user’s latitudes as shown
in Fig. 6. For $53^{\circ}$ inclined orbits, there is a minor decline in data
rate which is followed by a sharp rise due to a decrease in the number of
visible interfering satellites when the user leaves the inclination limits.

Figure 6: Data rate on different users’ latitudes, $K=10$,
$r_{\mathrm{min}}=500$ km, $\iota=53^{\circ}$, $m=2$, and
$\theta_{\mathrm{min}}=10^{\circ}$. The expression given in Theorem 2 is used
to plot this figure. Figure 7: Effect of altitude on data rate when $K=10$,
$\phi_{\mathrm{u}}=\\{0^{\circ},25^{\circ},65^{\circ}\\}$, $\iota=53^{\circ}$,
$m=2$, and $\theta_{\mathrm{min}}=10^{\circ}$. The expression given in Theorem
2 is used to plot this figure.

The effect of altitude on data rate is depicted in Fig. 7 for $K=10$. Similar
to Fig. 4, from the minimum altitude at which the user is able to visit at
least one satellite, the data rate increases rapidly until reaches a maximum
point. After the maximum point, the data rate decreases more slowly due to the
increase in the number of interfering satellites as well as satellites being
at a farther distance from the user. The altitude which maximizes the data
rate varies with the user’s latitude and, obviously, it has similar value
which results in the maximum coverage probability.

##  V Conclusions

In this paper, a generic approach to obtain the analytical performance of a
massive low Earth orbit network was presented by modeling the satellites’
locations as a nonhomogeneous Poisson point process and utilizing the tools
from stochastic geometry. The density of the nonhomogeneous Poisson point
process is derived from the actual geometry of the constellation which enables
us to take into account the non-uniform distribution of satellites across
different latitudes. Our next step was to apply this model to derive
analytical expressions for the coverage probability and average data rate of
an arbitrarily located user in terms of the distribution of fading and
shadowing as well as the Laplace function of interference. From numerical
results, we concluded that shadowing affects the coverage probability
positively as the probability of the user being in outage increases.
Furthermore, we showed how the analytical results allow one to find — without
involved orbital simulations — optimum values for altitude, the number of
orthogonal frequency channels, and user’s latitude which result in the largest
coverage and/or throughput, given constellation parameters.

###  -A Proof of Lemma 3

For an NPPP, the CDF of $R_{\mathrm{0}}$ can be written as

$\displaystyle
F_{R_{\mathrm{0}}}\left(r_{\mathrm{0}}\right)=1-{\mathbb{P}}(R_{\mathrm{0}}>r_{\mathrm{0}})=1-{\mathbb{P}}(\mathcal{N}=0),$
(24)

where ${\mathbb{P}}(\mathcal{N}=0)$ is the void probability of PPP in ${\cal
A}(r_{\mathrm{0}})$ that can be obtained from (5) by setting $n=0$. According
to (5), by integrating from the intensity over the spherical cap above the
user, we have

$\displaystyle F_{R_{0}}\left(r_{0}\right)$
$\displaystyle=1-\exp\Bigg{(}-\int_{\max(\phi_{\mathrm{u}}-\theta_{\mathrm{s}},-\iota)}^{\min(\phi_{\mathrm{u}}+\theta_{\mathrm{s}},\iota)}\beta(\phi_{\mathrm{s}})\delta(\phi_{\mathrm{s}})\left(r_{\mathrm{min}}+r_{\oplus}\right)^{2}\cos(\phi_{\mathrm{s}})d{\phi_{\mathrm{s}}}\Bigg{)}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}1-\exp\Bigg{(}-2\left(r_{\mathrm{min}}+r_{\oplus}\right)^{2}\int_{\max(\phi_{\mathrm{u}}-\theta_{\mathrm{s}},-\iota)}^{\min(\phi_{\mathrm{u}}+\theta_{\mathrm{s}},\iota)}\delta(\phi_{\mathrm{s}})\cos(\phi_{\mathrm{s}})$
$\displaystyle\hskip
85.35826pt\times\cos^{-1}\left({\frac{\cos(\theta_{\mathrm{s}})}{\cos(\phi_{\mathrm{s}}-\phi_{\mathrm{u}})}}\right)d\phi_{\mathrm{s}}\Bigg{)},$
(25)

where $\beta(\phi_{\mathrm{s}})$ is the longitude range inside the spherical
cap above the user at latitude $\phi_{\mathrm{s}}$. Equality (a) follows from
substitution of $\beta(\phi_{\mathrm{s}})$ using the basic geometry. Taking
the derivative of (-A) with respect to $r_{\mathrm{0}}$ completes the proof of
Lemma 3. Note that for $\theta_{\mathrm{s}}\leq|{\phi_{\mathrm{u}}}|-\iota$
the CDF given in (-A) is zero since the spherical cap formed by polar angle
$\theta_{\mathrm{s}}$ above the latitude $\phi_{\mathrm{u}}$ is much farther
from the constellation’s borders to contain any satellite.

###  -B Proof of Theorem 1

To obtain the expression given (1), let us begin with the definition of
coverage probability:

$\displaystyle P_{\mathrm{c}}\left(T\right)$
$\displaystyle=\mathbb{E}_{R_{\mathrm{0}}}\left[{\mathbb{P}}\left(\mathrm{SINR}>T|R_{\mathrm{0}}\right)\right]$
(26)
$\displaystyle=\int_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}{\mathbb{P}}\left(\mathrm{SINR}>T|R_{\mathrm{0}}=r_{\mathrm{0}}\right)f_{R_{\mathrm{0}}}\left(r_{\mathrm{0}}\right)dr_{\mathrm{0}}$
$\displaystyle=\int_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}{\mathbb{P}}\left(\frac{p_{\mathrm{0}}G_{\mathrm{0}}\mathcal{X}_{\mathrm{0}}r_{\mathrm{0}}^{-\alpha}}{I+\sigma^{2}}>T\right)f_{R_{\mathrm{0}}}\left(r_{\mathrm{0}}\right)\,dr_{\mathrm{0}}$
$\displaystyle=\int_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}\mathbb{E}_{I}\left[{\mathbb{P}}\left(G_{\mathrm{0}}\mathcal{X}_{\mathrm{0}}>\frac{Tr_{\mathrm{0}}^{\alpha}\left(I+\sigma^{2}\right)}{p_{\mathrm{0}}}\right)\bigg{|}I>0\right]$
$\displaystyle\hskip 42.67912pt\times
f_{R_{\mathrm{0}}}\left(r_{\mathrm{0}}\right)dr_{\mathrm{0}}.$

Since satellites with elevation angles smaller than $\theta_{\mathrm{min}}$
are not visible to the user, the integral has an upper limit. Then

$\displaystyle P_{\mathrm{c}}\left(T\right)=$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\int_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}\mathbb{E}_{I}\left[\int_{0}^{\infty}f_{\mathcal{X}_{\mathrm{0}}}(\mathit{x}_{\mathrm{0}})\left(1-F_{G_{\mathrm{0}}}\left(\frac{Tr_{\mathrm{0}}^{\alpha}\left(I+\sigma^{2}\right)}{p_{\mathrm{0}}\mathit{x}_{\mathrm{0}}}\right)\right)\right]$
$\displaystyle\hskip 42.67912pt\times
f_{R_{\mathrm{0}}}\left(r_{\mathrm{0}}\right)\,d\mathit{x}_{\mathrm{0}}dr_{\mathrm{0}}$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}\int_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}\mathbb{E}_{I}\left[\int_{0}^{\infty}f_{\mathcal{X}_{\mathrm{0}}}(\mathit{x}_{\mathrm{0}})\left(\frac{\Gamma\left(m_{0},m_{0}\frac{Tr_{\mathrm{0}}^{\alpha}\left(I+\sigma^{2}\right)}{p_{\mathrm{0}}\mathit{x}_{\mathrm{0}}}\right)}{\Gamma\left(m_{0}\right)}\right)\right]$
$\displaystyle\hskip 42.67912pt\times
f_{R_{\mathrm{0}}}\left(r_{\mathrm{0}}\right)\,d\mathit{x}_{\mathrm{0}}dr_{\mathrm{0}}$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{=}}\int_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}\int_{0}^{\infty}f_{\mathcal{X}_{\mathrm{0}}}(\mathit{x}_{\mathrm{0}})f_{R_{\mathrm{0}}}\left(r_{\mathrm{0}}\right)e^{-\frac{m_{0}Tr_{\mathrm{0}}^{\alpha}\sigma^{2}}{p_{\mathrm{0}}\mathit{x}_{\mathrm{0}}}}\mathbb{E}_{I}\Bigg{[}e^{-\frac{m_{0}Tr_{\mathrm{0}}^{\alpha}I}{p_{\mathrm{0}}\mathit{x}_{\mathrm{0}}}}$
$\displaystyle\sum_{k=0}^{m_{0}-1}\frac{\sum_{l=0}^{k}\binom{k}{l}\left(m_{0}\frac{Tr_{\mathrm{0}}^{\alpha}\sigma^{2}}{p_{\mathrm{0}}\mathit{x}_{\mathrm{0}}}\right)^{l}\left(m_{0}\frac{Tr_{\mathrm{0}}^{\alpha}I}{p_{\mathrm{0}}\mathit{x}_{\mathrm{0}}}\right)^{k-l}}{k!}\Bigg{]}\,d\mathit{x}_{\mathrm{0}}dr_{\mathrm{0}}$
$\displaystyle=\int_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}\int_{0}^{\infty}f_{\mathcal{X}_{\mathrm{0}}}(\mathit{x}_{\mathrm{0}})f_{R_{\mathrm{0}}}\left(r_{\mathrm{0}}\right)\Bigg{[}e^{-s\sigma^{2}}$
$\displaystyle\sum_{k=0}^{m_{0}-1}\frac{\sum_{l=0}^{k}\binom{k}{l}\left(s\sigma^{2}\right)^{l}\left(-s\right)^{k-l}\frac{\partial^{k-l}}{\partial
s^{k-l}}\mathcal{L}_{I}(s)}{k!}\Bigg{]}_{s=\frac{m_{0}Tr_{\mathrm{0}}^{\alpha}}{p_{\mathrm{0}}\mathit{x}_{\mathrm{0}}}}\,d\mathit{x}_{\mathrm{0}}dr_{\mathrm{0}}.$
(27)

The substitution in (a) follows from the product distribution of two
independent random variables, (b) follows from the distribution of gamma
random variable $G_{\mathrm{0}}$ (being the square of the Nakagami random
variable), and (c) is calculated by applying the incomplete gamma function for
integer values of $m_{0}$ to (b).

###  -C Proof of Theorem 2

Most of the steps in derivation of the data rate are similar to those given in
Appendix B. According to the definition of the average data rate given in
(15), we have

$\displaystyle\bar{C}=\mathbb{E}_{I,G_{\mathrm{0}},\mathcal{X}_{\mathrm{0}},R_{\mathrm{0}}}\left[\log_{2}\left(1+\mathrm{SINR}\right)\right]$
$\displaystyle=\int_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}\mathbb{E}_{I,G_{\mathrm{0}},\mathcal{X}_{\mathrm{0}}}\left[\log_{2}\left(1+\frac{p_{\mathrm{0}}G_{\mathrm{0}}\mathcal{X}_{\mathrm{0}}r_{\mathrm{0}}^{-\alpha}}{\sigma^{2}+I}\right)\right]f_{R_{\mathrm{0}}}(r_{\mathrm{0}})\,dr_{\mathrm{0}}$
$\displaystyle\stackrel{{\scriptstyle}}{{=}}\int_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}\int_{0}^{\infty}\mathbb{E}_{I,G_{\mathrm{0}},\mathcal{X}_{\mathrm{0}}}\left[{\mathbb{P}}\left(\log_{2}\left(1+\frac{p_{\mathrm{0}}G_{\mathrm{0}}\mathcal{X}_{\mathrm{0}}r_{\mathrm{0}}^{-\alpha}}{\sigma^{2}+I}\right)>t\right)\right]$
$\displaystyle\hskip 85.35826pt\times
f_{R_{\mathrm{0}}}(r_{\mathrm{0}})\,dtdr_{\mathrm{0}},$ (28)

where the last step follows form the fact that for a positive random variable
$X$, $\mathbb{E}\left[X\right]=\int_{t>0}{\mathbb{P}}\left(X>t\right)dt$.
Thus, we have

$\displaystyle\bar{C}=\int_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}\int_{0}^{\infty}\mathbb{E}_{I,G_{\mathrm{0}},\mathcal{X}_{\mathrm{0}}}\left[{\mathbb{P}}\left(G_{\mathrm{0}}\mathcal{X}_{\mathrm{0}}>\frac{r_{\mathrm{0}}^{\alpha}\left(\sigma^{2}+I\right)}{p_{\mathrm{0}}}\left(2^{t}-1\right)\right)\right]$
$\displaystyle\hskip 85.35826pt\times
f_{R_{\mathrm{0}}}(r_{\mathrm{0}})\,dtdr_{\mathrm{0}}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\int_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}\int_{0}^{\infty}\mathbb{E}_{I}\Bigg{[}\int_{0}^{\infty}f_{\mathcal{X}_{\mathrm{0}}}(\mathit{x}_{\mathrm{0}})$
$\displaystyle\times\left(1-F_{G_{\mathrm{0}}}\left(\frac{r_{\mathrm{0}}^{\alpha}\left(I+\sigma^{2}\right)\left(2^{t}-1\right)}{p_{\mathrm{0}}\mathit{x}_{\mathrm{0}}}\right)\right)d\mathit{x}_{\mathrm{0}}\Bigg{]}f_{R_{\mathrm{0}}}(r_{\mathrm{0}})\,dtdr_{\mathrm{0}}$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}\int_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}\int_{0}^{\infty}\mathbb{E}_{I}\Bigg{[}\int_{0}^{\infty}f_{\mathcal{X}_{\mathrm{0}}}(\mathit{x}_{\mathrm{0}})$
$\displaystyle\times\left(\frac{\Gamma\left(m_{0},m_{0}\frac{r_{\mathrm{0}}^{\alpha}\left(I+\sigma^{2}\right)\left(2^{t}-1\right)}{p_{\mathrm{0}}\mathit{x}_{\mathrm{0}}}\right)}{\Gamma\left(m_{0}\right)}\right)d\mathit{x}_{\mathrm{0}}\Bigg{]}f_{R_{\mathrm{0}}}(r_{\mathrm{0}})\,dtdr_{\mathrm{0}}$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{=}}\int_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}\int_{0}^{\infty}\int_{0}^{\infty}f_{\mathcal{X}_{\mathrm{0}}}(\mathit{x}_{\mathrm{0}})f_{R_{\mathrm{0}}}(r_{\mathrm{0}})e^{-\frac{m_{0}\left(2^{t}-1\right)r_{\mathrm{0}}^{\alpha}\sigma^{2}}{p_{\mathrm{0}}\mathit{x}_{\mathrm{0}}}}$
$\displaystyle\times\mathbb{E}_{I}\Bigg{[}e^{-\frac{m_{0}\left(2^{t}-1\right)r_{\mathrm{0}}^{\alpha}I}{p_{\mathrm{0}}\mathit{x}_{\mathrm{0}}}}$
$\displaystyle\sum_{k=0}^{m_{0}-1}\frac{\sum_{l=0}^{k}\binom{k}{l}\left(\frac{m_{0}\left(2^{t}-1\right)r_{\mathrm{0}}^{\alpha}\sigma^{2}}{p_{\mathrm{0}}\mathit{x}_{\mathrm{0}}}\right)^{l}\left(\frac{m_{0}\left(2^{t}-1\right)r_{\mathrm{0}}^{\alpha}I}{p_{\mathrm{0}}\mathit{x}_{\mathrm{0}}}\right)^{k-l}}{k!}\Bigg{]}dtd\mathit{x}_{\mathrm{0}}dr_{\mathrm{0}}$
$\displaystyle=\int_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}\int_{0}^{\infty}\int_{0}^{\infty}f_{\mathcal{X}_{\mathrm{0}}}(\mathit{x}_{\mathrm{0}})f_{R_{\mathrm{0}}}\left(r_{\mathrm{0}}\right)\Bigg{[}e^{-s\sigma^{2}}$
$\displaystyle\sum_{k=0}^{m_{0}-1}\frac{\sum_{l=0}^{k}\binom{k}{l}\left(s\sigma^{2}\right)^{l}\left(-s\right)^{k-l}\frac{\partial^{k-l}}{\partial
s^{k-l}}\mathcal{L}_{I}(s)}{k!}\Bigg{]}_{s=\frac{m_{0}\left(2^{t}-1\right)r_{\mathrm{0}}^{\alpha}}{p_{\mathrm{0}}\mathit{x}_{\mathrm{0}}}}dtd\mathit{x}_{\mathrm{0}}dr_{\mathrm{0}}.$
(29)

Similar to the proof of Theorem 1, (a) follows from the product distribution
of two independent random variables, (b) follows from the gamma distribution
of serving channel gain $G_{\mathrm{0}}$, and (c) is calculated by applying
the incomplete gamma function for integer values of $m_{0}$ to (b).

###  -D Proof of Lemma 5

In this appendix, we derive the expression for Laplace function of
interference assuming arbitrary distributions for fading and shadowing. Let us
start with the definition of Laplace function for random variable $I$ which is

$\displaystyle\mathcal{L}_{I}(s)\triangleq\mathbb{E}_{I}\left[e^{-sI}\right]$
$\displaystyle=\mathbb{E}_{\mathcal{N},R_{n},\mathcal{X}_{n},G_{n}}\left[\exp{\left(-s\sum_{n\in\xi\backslash\\{\mathrm{s}\\}}p_{n}G_{n}\mathcal{X}_{n}R_{n}^{-\alpha}\right)}\right]$
$\displaystyle=\mathbb{E}_{\mathcal{N},R_{n},\mathcal{X}_{n},G_{n}}\left[\prod_{n\in\xi\backslash\\{\mathrm{s}\\}}\exp{\left(-sp_{n}G_{n}\mathcal{X}_{n}R_{n}^{-\alpha}\right)}\right]$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\mathbb{E}_{\mathcal{N},R_{n}}\left[\prod_{n\in\xi\backslash\\{\mathrm{s}\\}}\mathbb{E}_{\mathcal{X}_{n},G_{n}}\left[\exp{\left(-sp_{n}G_{n}\mathcal{X}_{n}R_{n}^{-\alpha}\right)}\right]\right]$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}\mathbb{E}_{\mathcal{N}}\Bigg{[}\prod_{n\in\xi\backslash\\{\mathrm{s}\\}}\int_{r_{\mathrm{0}}}^{r_{\mathrm{max}}}\mathbb{E}_{\mathcal{X}_{n},G_{n}}\left[\exp{\left(-sp_{n}G_{n}\mathcal{X}_{n}r_{n}^{-\alpha}\right)}\right]$
$\displaystyle\hskip 85.35826pt\times
f_{R_{n}|R_{\mathrm{0}}}(r_{n}|r_{\mathrm{0}})dr_{n}\Bigg{]},$ (30)

where (a) follows from the i.i.d. distribution of $G_{n}$ and
$\mathcal{X}_{n}$ as well as their independence from $\mathcal{N}$ and $R_{n}$
and (b) is obtained by taking the expectation over the random variable $R_{n}$
conditioned on $R_{\mathrm{0}}$. Then

$\displaystyle\mathcal{L}_{I}(s)\stackrel{{\scriptstyle(c)}}{{=}}\sum_{n=0}^{\infty}P_{n}\left(\mathcal{A}\left(r_{\mathrm{max}}\right)-\mathcal{A}\left(r_{\mathrm{0}}\right)\right)$
$\displaystyle\hskip
2.84526pt\times\Bigg{(}\int_{r_{\mathrm{0}}}^{r_{\mathrm{max}}}\mathbb{E}_{\mathcal{X}_{n},G_{n}}\left[\exp{\left(-sp_{n}G_{n}\mathcal{X}_{n}r_{n}^{-\alpha}\right)}\right]f_{R_{n}|R_{\mathrm{0}}}(r_{n}|r_{\mathrm{0}})dr_{n}\Bigg{)}^{n}$
$\displaystyle\stackrel{{\scriptstyle(d)}}{{=}}\sum_{n=0}^{\infty}P_{n}\left(\mathcal{A}\left(r_{\mathrm{max}}\right)-\mathcal{A}\left(r_{\mathrm{0}}\right)\right)$
$\displaystyle\times\Bigg{(}\int_{r_{\mathrm{0}}}^{r_{\mathrm{max}}}\int_{0}^{\infty}\mathcal{L}_{G_{n}}(sp_{n}x_{n}r_{n}^{-\alpha})f_{\mathcal{X}_{n}}(x_{n})f_{R_{n}|R_{\mathrm{0}}}(r_{n}|r_{\mathrm{0}})d\mathit{x}_{\mathrm{0}}dr_{n}\Bigg{)}^{n},$
(31)

where
$\mathcal{A}\left(r_{\mathrm{max}}\right)-\mathcal{A}\left(r_{\mathrm{0}}\right)$
indicates the region above the user where satellites which are more distanced
from the user than the serving satellite exist, (c) is obtained by taking the
average over the Poisson random variable $\mathcal{N}$, and applying the law
of total expectation on independent random variables $G_{n}$ and
$\mathcal{X}_{n}$ results in (d).

## References

  * [1] M. Giordani and M. Zorzi, “Non-terrestrial networks in the 6G era: Challenges and opportunities,” _IEEE Network_ , vol. 35, no. 2, pp. 244–251, Mar. 2021.
  * [2] Z. Jia, M. Sheng, J. Li, D. Zhou, and Z. Han, “Joint HAP access and LEO satellite backhaul in 6G: Matching game-based approaches,” _IEEE Journal on Selected Areas in Communications_ , vol. 39, no. 4, pp. 1147–1159, Apr. 2021.
  * [3] M. Giordani and M. Zorzi, “Satellite communication at millimeter waves: a key enabler of the 6G era,” in _Proc. International Conference on Computing, Networking and Communications_ , Feb. 2020.
  * [4] M. Giordani, M. Polese, M. Mezzavilla, S. Rangan, and M. Zorzi, “Toward 6G networks: Use cases and technologies,” _IEEE Communications Magazine_ , vol. 58, no. 3, pp. 55–61, Mar. 2020.
  * [5] J. L. Grubb, “The traveler’s dream come true,” _IEEE Communications Magazine_ , vol. 29, no. 11, pp. 48–51, Nov. 1991.
  * [6] A. Mokhtar and M. Azizoglu, “On the downlink throughput of a broadband LEO satellite network with hopping beams,” _IEEE Communications Letters_ , vol. 4, no. 12, pp. 390–393, Dec. 2000.
  * [7] F. Vatalaro, G. E. Corazza, C. Caini, and C. Ferrarelli, “Analysis of LEO, MEO, and GEO global mobile satellite systems in the presence of interference and fading,” _IEEE Journal on Selected Areas in Communications_ , vol. 13, no. 2, pp. 291–300, Feb. 1995.
  * [8] A. Ganz, Y. Gong, and B. Li, “Performance study of low Earth-orbit satellite systems,” _IEEE Transactions on Communications_ , vol. 42, no. 234, pp. 1866–1871, Feb. 1994.
  * [9] I. Ali, N. Al-Dhahir, and J. E. Hershey, “Predicting the visibility of LEO satellites,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 35, no. 4, pp. 1183–1190, Oct. 1999.
  * [10] Y. Seyedi and S. M. Safavi, “On the analysis of random coverage time in mobile LEO satellite communications,” _IEEE Communications Letters_ , vol. 16, no. 5, pp. 612–615, May 2012.
  * [11] Z. Qu, G. Zhang, H. Cao, and J. Xie, “LEO satellite constellation for Internet of Things,” _IEEE Access_ , vol. 5, pp. 18 391–18 401, Sep. 2017.
  * [12] M. Haenggi, _Stochastic geometry for wireless networks_. Cambridge University Press, 2012.
  * [13] B. Błaszczyszyn, M. Haenggi, P. Keeler, and S. Mukherjee, _Stochastic geometry analysis of cellular networks_. Cambridge University Press, 2018.
  * [14] M. Haenggi, J. G. Andrews, F. Baccelli, O. Dousse, and M. Franceschetti, “Stochastic geometry and random graphs for the analysis and design of wireless networks,” _IEEE Journal on Selected Areas in Communications_ , vol. 27, no. 7, pp. 1029–1046, Sep. 2009.
  * [15] H. ElSawy, E. Hossain, and M. Haenggi, “Stochastic geometry for modeling, analysis, and design of multi-tier and cognitive cellular wireless networks: A survey,” _IEEE Communications Surveys and Tutorials_ , vol. 15, no. 3, pp. 996–1019, Jun. 2013.
  * [16] J. G. Andrews, F. Baccelli, and R. K. Ganti, “A tractable approach to coverage and rate in cellular networks,” _IEEE Transactions on Communications_ , vol. 59, no. 11, pp. 3122–3134, Nov. 2011.
  * [17] H. S. Dhillon, R. K. Ganti, F. Baccelli, and J. G. Andrews, “Modeling and analysis of k-tier downlink heterogeneous cellular networks,” _IEEE Journal on Selected Areas in Communications_ , vol. 30, no. 3, pp. 550–560, Apr. 2012.
  * [18] D. Cao, S. Zhou, and Z. Niu, “Optimal combination of base station densities for energy-efficient two-tier heterogeneous cellular networks,” _IEEE Transactions on Wireless Communications_ , vol. 12, no. 9, pp. 4350–4362, Sep. 2013.
  * [19] H. S. Dhillon, T. D. Novlan, and J. G. Andrews, “Coverage probability of uplink cellular networks,” in _Proc. IEEE Global Communications Conference_ , Dec. 2012.
  * [20] M. Afshang and H. S. Dhillon, “Fundamentals of modeling finite wireless networks using binomial point process,” _IEEE Transactions on Wireless Communications_ , vol. 16, no. 5, pp. 3355–3370, May 2017.
  * [21] Z. Pan and Q. Zhu, “Modeling and analysis of coverage in 3-D cellular networks,” _IEEE Communications Letters_ , vol. 19, no. 5, pp. 831–834, May 2015.
  * [22] V. V. Chetlur and H. S. Dhillon, “Downlink coverage analysis for a finite 3-D wireless network of unmanned aerial vehicles,” _IEEE Transactions on Communications_ , vol. 65, no. 10, pp. 4543–4558, Oct. 2017.
  * [23] S. Srinivasa and M. Haenggi, “Distance distributions in finite uniformly random networks: Theory and applications,” _IEEE Transactions on Vehicular Technology_ , vol. 59, no. 2, pp. 940–949, Feb. 2010.
  * [24] J. Guo, S. Durrani, and X. Zhou, “Outage probability in arbitrarily-shaped finite wireless networks,” _IEEE Transactions on Communications_ , vol. 62, no. 2, pp. 699–712, Feb. 2014.
  * [25] Z. Khalid and S. Durrani, “Distance distributions in regular polygons,” _IEEE Transactions on Vehicular Technology_ , vol. 62, no. 5, pp. 2363–2368, Jun. 2013.
  * [26] X. Wang, H. Zhang, Y. Tian, and V. C. M. Leung, “Modeling and analysis of aerial base station-assisted cellular networks in finite areas under LoS and NLoS propagation,” _IEEE Transactions on Wireless Communications_ , vol. 17, no. 10, pp. 6985–7000, Oct. 2018.
  * [27] J. Guo, S. Durrani, and X. Zhou, “Performance analysis of arbitrarily-shaped underlay cognitive networks: Effects of secondary user activity protocols,” _IEEE Transactions on Communications_ , vol. 63, no. 2, pp. 376–389, Feb. 2015.
  * [28] T. A. Khan and M. Afshang, “A stochastic geometry approach to Doppler characterization in a LEO satellite network,” in _Proc. IEEE International Conference on Communications_ , 2020.
  * [29] J. Hu, G. Li, D. Bian, L. Gou, and C. Wang, “Optimal power control for cognitive LEO constellation with terrestrial networks,” _IEEE Communications Letters_ , vol. 24, no. 3, pp. 622–625, Mar. 2020.
  * [30] J. Zhang, B. Evans, M. A. Imran, X. Zhang, and W. Wang, “Performance analysis of C/U split hybrid satellite terrestrial network for 5G systems,” in _Proc. 20th IEEE International Workshop on Computer Aided Modelling and Design of Communication Links and Networks_ , Sep. 2015.
  * [31] A. J. Roumeliotis, C. I. Kourogiorgas, and A. D. Panagopoulos, “Optimal dynamic capacity allocation for high throughput satellite communications systems,” _IEEE Wireless Communications Letters_ , vol. 8, no. 2, pp. 596–599, Apr. 2019.
  * [32] A. K. Dwivedi, S. P. Chokkarapu, S. Chaudhari, and N. Varshney, “Performance analysis of novel direct access schemes for LEO satellites based IoT network,” in _Proc. IEEE 31st Annual International Symposium on Personal, Indoor and Mobile Radio Communications_ , Sep. 2020.
  * [33] N. Okati, T. Riihonen, D. Korpi, I. Angervuori, and R. Wichman, “Downlink coverage and rate analysis of low Earth orbit satellite constellations using stochastic geometry,” _IEEE Transactions on Communications_ , vol. 68, no. 8, pp. 5120–5134, Aug. 2020.
  * [34] N. Okati and T. Riihonen, “Stochastic analysis of satellite broadband by mega-constellations with inclined LEOs,” in _Proc. IEEE 31st Annual International Symposium on Personal, Indoor and Mobile Radio Communications_ , Sep. 2020.
  * [35] A. Talgat, M. A. Kishk, and M.-S. Alouini, “Nearest neighbor and contact distance distribution for binomial point process on spherical surfaces,” _IEEE Communications Letters_ , vol. 24, no. 12, pp. 2659–2663, Dec. 2020\.
  * [36] ——, “Stochastic geometry-based analysis of LEO satellite communication systems,” _IEEE Communications Letters_ , pp. 1–1, Oct. 2020.
  * [37] A. Yastrebova, I. Angervuori, N. Okati, M. Vehkaperä, M. Höyhtyä, R. Wichman, and T. Riihonen, “Theoretical and simulation-based analysis of terrestrial interference to LEO satellite uplinks,” in _Proc. IEEE Global Communications Conference_ , Dec. 2020.
  * [38] N. Okati and T. Riihonen, “Modeling and analysis of LEO mega-constellations as nonhomogeneous Poisson point processes,” in _Proc. IEEE 93rd Vehicular Technology Conference_ , Apr. 2021.
  * [39] I. S. Gradshteyn and I. M. Ryzhik, _Table of Integrals, Series, and Products_. Academic Press, 2014.


# Light quark energy loss in a soft-wall AdS/QCD model

Xiangrong Zhu<EMAIL_ADDRESS>School of Science, Huzhou University,
Huzhou 313000, China  Zi-qiang Zhang<EMAIL_ADDRESS>School of Mathematics
and Physics, China University of Geosciences, Wuhan 430074, China

###### Abstract

We investigate the energy loss of light quarks in a holographic QCD model with
conformal invariance broken by a background dilaton. We perform the analysis
within falling string and shooting string, respectively. It turns out that the
two methods give the same result: the presence of chemical potential and
confining scale tends to enhance the energy loss, in accord with previous
findings of drag force and jet quenching parameter.

###### pacs:

11.25.Tq, 11.15.Tk, 11.25-w

##  I Introduction

One of the main purposes of the heavy-ion collisions experiments is to explore
the QCD phase diagram and the properties of new state of matter produced
through collisions at high energy density. It is believed that the
experimental program at Relativistic Heavy-Ion Collider (RHIC) and Large
Hadron Collider (LHC) have created a new state of matter so-called quark gluon
plasma (QGP) ev ; mg . One of the striking features of such substance is jet
quenching: the energy loss of high energy partons produced through collisions
as they interact with the plasma before they fragment into hadrons (for recent
reviews see mc ; gy ). On the other hand, it has been observed that QGP
behaves as a strongly coupled fluid es ; kad ; ja , which involves
nonperturbative physics suitable for the application of the anti-de
Sitter/conformal field theory (AdS/CFT) correspondence.

AdS/CFT Maldacena:1997re ; Gubser:1998bc ; MadalcenaReview , the duality
between the type IIB superstring theory formulated on AdS${}_{5}\times S^{5}$
and $\mathcal{N}=4$ super Yang-Mills theory (SYM) in four dimensions, provides
a useful tool for understanding the strong interaction, i.e. quantum
chromodynamics (QCD). Although SYM differs from QCD in many properties at zero
temperature, it reveals some qualitative features of QCD in the strongly
coupled regime at non-zero temperature. During the last two decades, the
AdS/CFT correspondence has yielded many important insights for studying
various aspects of QGP (see JCA ; OD  for recent reviews with many
phenomenological applications). An interesting example of such applications is
jet quenching. For instance, the drag force which describes the energy loss
for heavy quarks moving through $\mathcal{N}=4$ SYM plasma was studied in cp ;
ss0 . Moreover, the jet quenching of light quarks moving through
$\mathcal{N}=4$ SYM plasma has been addressed based on various approaches,
e.g., jet quenching parameter hl ; hl1 , falling string ss ; pm ; pm1 ; pa ;
pa1 , shooting string ss1 ; ss2 , etc.

Here we present an alternative holographic approach to study the energy loss
of light quarks. The motivation is that AdS/QCD models, e.g., hard wall H1 ;
H2 , soft wall AKE  and some improved holographic models JP ; AST ; DL ; DL1 ;
SH ; SH1 ; RRO  could provide a nice phenomenological description of quark-
antiquark interaction and some hadronic properties. In particular, we will
employ the SWT,μ model PCO  which is defined by the AdS with a charged black
hole to describe finite temperature and density multiplied by a warp factor to
generate confinement. It turns out that this model has achieved considerable
success in describing various aspects of hadron physics PCO ; CPA ; PCO1 ;
PCO2 ; YH ; zq0 ; XCH ; zq . Another motivation for this paper is that the
drag force and jet quenching parameter have been studied in the SWT,μ model YH
; zq0  and the results show that the presence of chemical potential and
confining scale increases the two parameters thus enhancing the energy loss.
Inspired by this, we wonder whether chemical potential and confining scale
have the same effect on the energy loss of light quarks using falling string
and shooting string.

The outline of the paper is as follows. In the next section, we introduce the
SWT,μ model presented in PCO . In section 3, we study the energy loss of light
quarks in this model within falling string and shooting string, in turn. The
last part is devoted to conclusion and discussion.

##  II Setup

In this section, we briefly review the SWT,μ model given in PCO . In the
string frame, the SWT,μ model has the following metric

$ds^{2}=\frac{R^{2}}{z^{2}}h(z)(-f(z)dt^{2}+d\vec{x}^{2}+\frac{dz^{2}}{f(z)}),$
(1)

with

$f(z)=1-(1+Q^{2})(\frac{z}{z_{h}})^{4}+Q^{2}(\frac{z}{z_{h}})^{6},\qquad
h(z)=e^{c^{2}z^{2}},$ (2)

where $R$ is the AdS radius. $Q$ is the charge of black hole. $z$ represents
the fifth coordinate with $z=z_{h}$ the horizon, defined by $f(z_{h})=0$. The
warp factor $h(z)$, characterizing the soft wall model, distorts the
background metric and brings the confining scale $c$ (see JP  for a
anatlytical way to introduce the warp factor within potential reconstruction
approach).

The temperature of the black hole reads

$T=\frac{1}{\pi z_{h}}(1-\frac{Q^{2}}{2}),\qquad 0\leq Q\leq\sqrt{2}.$ (3)

The chemical potential reads

$\mu=\sqrt{3}Q/z_{h}.$ (4)

Notice that for $Q=0$, the SWT,μ model becomes the Andreev model OA . For
$c=0$, it reduces to the AdS-Reissner Nordstrom black hole CV:1999 ; DT:2006 .
For $Q=c=0$, it restores to the AdS black hole.

##  III Energy loss of light quarks in the SWT,μ model

###  III.1 Falling string

In previous research, various authors ss ; pm ; pm1 ; pa ; pa1  have studied
the jet quenching of light quark in $\mathcal{N}=4$ SYM plasma using falling
string in different ways, e.g, with or without the addition of fundamental-
charge matter, specify different initial conditions. Though each approaches
the subject from different vantage points, the analysis of the stopping
distance traveled by the falling string in $AdS_{5}$-Schwarzschild are
consistent.

Next, we will follow the argument in pa ; pa1  to study the jet quenching of
light quark in the SWT,μ model by analyzing the stopping distance of an image
jet induced by a massless source field, characterized by a massless particle
falling along the null geodesic in the WKB approximation. According to this
scenario, the R-charged current is generated by a massless gauge field in the
gravity dual and the induced current is regarded as an energetic jet passing
through the medium. When the wave packet of the massless gauge field falls
into the horizon of the dual geometry, the image jet on the boundary
dissipates and then thermalizes in the medium. The stopping distance is
therefore defined as the (maximum) distance for a jet passing through the
medium before it thermalizes.

In the WKB approximation, the wave packet of the massless gauge field in the
gravity dual is supposed to be localized in the momentum space such that the
wave function of the gauge field could be factorized as

$A_{j}(t,z)=exp[\frac{i}{\hbar}(q_{k}x_{k}+\int dzq_{z})]\tilde{A}_{j}(t,z),$
(5)

where $q_{k}$ denotes the 4-momentum, conserved as the metric preserves the
translational symmetry along the 4-dimensional spacetime. $q_{z}$ represents
the momentum along the bulk direction. $\tilde{A}_{j}(t,z)$ refers to the
slow-varying with respect to $t$ and $z$. $j,k$ are the 4-dimensional
spacetime coordinates.

In the classical limit, i.e., $\hbar\rightarrow 0$, the equation of motion of
the wave pack reduces to a null geodesic,

$0=(ds^{2})=dx^{i}g_{ij}dx^{j}+dzg_{zz}dz,$ (6)

yielding

$\frac{dz}{d\zeta}=\frac{1}{\sqrt{g_{zz}}}[-g_{ij}\frac{dx^{i}}{d\zeta}\frac{dx^{j}}{d\zeta}]^{1/2},$
(7)

where $\zeta$ is an affine parameter for the trajectory. As the 4-dimensional
translation invariance,

$g_{ij}\frac{dx^{j}}{d\zeta},$ (8)

is conserved and proportional to $q_{i}$, leading to

$\frac{dx^{i}}{d\zeta}\propto g^{ij}q_{j}.$ (9)

Then, dividing (9) by (7) gives

$\frac{dx^{i}}{dz}=\sqrt{g_{zz}}\frac{g^{ij}q_{j}}{(-q_{k}g^{kl}q_{l})^{1/2}},$
(10)

one can check that the null geodesic in the above equation remains unchanged
even when one uses the Einstein frame.

To proceed, we calculate the stooping distance. Supposing that the 3-momentum
$\vec{q}$ to point in one of $\vec{x}$ directions, e.g., the $x_{3}$
direction, implying $q_{i}=(-\omega,0,0,|\vec{q}|)$, where $\omega$ and
$\vec{q}$ are the energy and spacial momentum of the light quark,
respectively. Then plugging (1) into (10), the stopping distance for the SWT,μ
model can be obtained,

$x=\int_{0}^{z_{h}}\frac{dz}{e^{c^{2}z^{2}}\sqrt{\frac{\omega^{2}}{|\vec{q}|^{2}}-\frac{1-(1+\frac{\mu^{2}z_{h}^{2}}{3})\frac{z^{4}}{z_{h}^{4}}+\frac{\mu^{2}z^{6}}{3z_{h}^{4}}}{e^{c^{4}z^{4}}}}},$
(11)

note that when one turns off the chemical potential and confining scale
effects by setting $\mu=c=0$, the above equation recovers the result of SYM pa
; pa1 .

Before going further, we need to turn to numerics. First, we determine the
value range of $c$. Here we tend to study the light quark energy loss in a
class of models parametrized by $c$ rather than in a specific model with fixed
$c$. To that end, we make $c$ dimensionless by normalizing it at fixed $T$ and
express $\mu$ in unit of it as well. The knowledge of lattice calculations
suggests HLL  that the range of $0\leq c/T\leq 2.5$ is most relevant for a
comparison with QCD. We take that range here.

In fig.1, we compare the stopping distance of a light quark moving in the
SWT,μ model with the same one moving in $\mathcal{N}=4$ SYM plasma. The left
panel corresponds to $x/x_{SYM}$ versus $\mu/T$ with fixed $c/T$ while the
right one represents $x/x_{SYM}$ versus $c/T$ with fixed $\mu/T$, where
$x_{SYM}$ denotes the stopping distance of SYM. From these figures, one finds
that the stopping distance in the SWT,μ model is smaller than that of SYM. In
particular, the left panel tells us that with fixed $c/T$, increasing $\mu/T$
leads to decreasing $x/x_{SYM}$. Namely, the inclusion of chemical potential
decreases the stopping distance thus enhancing the energy loss, in accord with
that found in EL . On the other hand, one can see from the right panel that
the confining scale has similar effect. Given the above, one could reach the
following conclusion: the presence of chemical potential and confining scale
both decrease the stopping distance thus enhancing the energy loss,
consistently with the findings of the drag force YH  and jet quenching
parameter zq0 .

Figure 1: Left: $x/x_{SYM}$ versus $\mu/T$. Right: $x/x_{SYM}$ versus $c/T$.
Here we take $|\vec{q}|=0.99\omega$.

###  III.2 Shooting string

This subsection is devoted to the analysis of the light quarks energy loss
using shooting string. According to ss1 ; ss2 , one considers a particular
type of classical string motion: the string endpoint is close to horizon
initially and move towards the boundary, carrying some energy and momentum
which are gradually bled off into the rest of the string during its rise,
hence, this motion is called finite-endpoint-momentum shooting string, or
shooting string for short.

Next, we will follow the approach in ss1 ; ss2  to study the light quark
energy loss within shooting string for the background metric (1). The
instantaneous energy loss of light quarks takes the form

$\frac{dE}{dx}=-\frac{|L|}{2\pi\alpha^{\prime}}\frac{1}{z^{2}},$ (12)

where $L$ represents the null geodesics that the endpoint follows. From the
above equation, it appears that small $z$ (meaning the endpoint starts near
the boundary) will result in large energy loss, indicating the jets will be
quenched quickly and cannot be seen. So to hedge against this, one assumes the
strings start close to the horizon.

As before, one supposes the quark moving along the $x_{3}$ direction. Given
that, the energy and momentum of the quark become

$E=-\frac{h(z)}{\eta}\frac{f(z)}{z^{2}},$ (13)

and

$p_{x_{3}}=\frac{h(z)}{\eta z^{2}}\frac{dx}{dt},$ (14)

with $\eta$ the auxiliary field. Then the null geodesics reads

$L=\frac{E}{p_{x_{3}}}=-f(z)\frac{dt}{dx}.$ (15)

The finite momentum endpoints will move along $ds^{2}=0$, which gives

$(\frac{dx}{dz})^{2}=\frac{1}{L^{2}-f(z)}.$ (16)

The denominator of (16) will vanish at $z=z_{*}$, yielding

$L=-\sqrt{f(z_{*})}.$ (17)

So the null geodesics equation becomes

$\frac{dx}{dz}=\frac{1}{\sqrt{f(z_{*})-f(z)}}.$ (18)

Finally, using (1), (12) and (17), one obtains the energy loss for the SWT,μ
model as

$\frac{dE}{dx}=-\frac{e^{c^{2}z^{2}}}{2\pi\alpha^{\prime}}\frac{\sqrt{1-(1+\frac{\mu^{2}z_{h}^{2}}{3})\frac{z_{*}^{4}}{z_{h}^{4}}+\frac{\mu^{2}z_{*}^{6}}{3z_{h}^{4}}}}{z^{2}},$
(19)

note that for $\mu=c=0$, the above equation recovers the result of SYM ss1 ;
ss2 .

Before going further, one needs to solve the null geodesics equation (18).
However, it is difficult to solve it analytically, but it is possible
numerically. The procedures of the numerical evaluation are as follows: First,
sending $z_{*}\rightarrow 0$, then for a given value of $\mu/T$, one can
integrate (18) and invert to get $z(x)$. Next, substituting $z(x)$ into (19)
one can obtain $dE/dx$ as a function of $c/T$, $\mu/T$ as well as $x$.

In fig.2, we compare the instantaneous energy loss of a light quark moving in
the SWT,μ model with its counterpart in $\mathcal{N}=4$ SYM plasma. The left
panel corresponds to $(dE/dx)/(dE/dx)_{SYM}$ versus $\mu/T$ with fixed $c/T$
while the right one denotes $(dE/dx)/(dE/dx)_{SYM}$ versus $c/T$ with fixed
$\mu/T$, where $(dE/dx)_{SYM}$ represents the energy loss of SYM. For both
panels, one finds the energy loss in the SWT,μ model is larger than that of
SYM. In particular, increasing $\mu/T$ or $c/T$ both increase
$(dE/dx)/(dE/dx)_{SYM}$. Namely, the inclusion of chemical potential and
confining scale both increase the energy loss, consistently with the analysis
of the stopping distance in the last subsection. The physical significance of
the results will be discussed in the next section.

Figure 2: Left: $(dE/dx)/(dE/dx)_{SYM}$ versus $\mu/T$. Right:
$(dE/dx)/(dE/dx)_{SYM}$ versus $c/T$. Here we take $x=0.01$.

##  IV conclusion and discussion

Jet quenching in high-energy heavy-ion collisions can be used to probe
properties of QGP. Studying the jet quenching in strongly coupled nonconformal
plasma with finite quark density may shed qualitative insights into analogous
questions in QCD. In this paper, we investigated the energy loss of light
quarks in the SWT,μ model using falling string and shooting string, in turn.
We discussed how chemical potential and confining scale modify the energy loss
for both cases, respectively. For the former, we calculated the stopping
distance of a massless particle moving along the null geodesic and found the
presence of chemical potential and confining scale both decrease the stopping
distance thus increasing the energy loss. For the latter, we computed the
instantaneous energy loss based on the finite endpoint momentum framework and
found the same result: increase chemical potential and confining scale both
enhance the energy loss. These results agree with those obtained by the drag
force and jet quenching parameter in the SWT,μ model YH ; zq0 . Taking all of
this together, one may draw a conclusion that the effects of chemical
potential and confining scale on the energy loss of heavy quarks and light
quarks are consistent.

We would like to make some comment on the possible significance of our results
to the current study on heavy-ion physics although such comparison may be
qualitative. Recent results of the experiments at RHIC and LHC indicate that
after collision, QGP expands and the jet quenching decreases with a decrease
in temperature. Since chemical potential and confining scale both have the
effect of increasing the energy loss, one may infer that including the two
effects may lower the required temperature region in order to observe the
obtained RHIC or LHC values for jet quenching.

However, it must be admitted that there are some limitations to the present
study. The major disadvantage is that the SWT,μ model is not a bonafide
solution of the classical equations of motion. Considering jet quenching in
some consistent models, e.g. JP ; AST ; DL ; DL1 ; SH ; SH1 ; RRO  would be
instructive. We hope to report our progress in this regard in the near future.

##  V Acknowledgments

This work is supported by Zhejiang Provincial Natural Science Foundation of
China Nos. LY19A050001, LY18A050002 and the NSFC under Grant No. 11947410.

## References

  * (1) E. V. Shuryak, Nucl. Phys. A 750, 64 (2005).
  * (2) M. Gyulassy and L. McLerran, Nucl. Phys. A750, 30 (2005).
  * (3) M. Connors, C. Nattrass, R. Reed and S. Salur, Rev. Mod. Phys. 90, 025005 (2018).
  * (4) G. Y. Qin, X.-N. Wang, Int. J. Mod. Phys. E 24 (2015) no.11, 1530014\.
  * (5) E. Shuryak, Prog. Part. Nucl. Phys. 53, 273 (2004).
  * (6) K. Adcox et al. (PHENIX Collaboration), Nucl. Phys. A757, 184 (2005).
  * (7) J. Adams et al. (STAR Collaboration), Nucl. Phys. A757, 102 (2005).
  * (8) J. M. Maldacena, Adv. Theor. Math. Phys. 2, 231 (1998).
  * (9) S. S. Gubser, I. R. Klebanov and A. M. Polyakov, Phys. Lett. B 428, 105 (1998).
  * (10) O. Aharony, S. S. Gubser, J. Maldacena, H. Ooguri and Y. Oz, Phys. Rept. 323, 183 (2000).
  * (11) J. C. Solana, H. Liu, D. Mateos, K. Rajagopal, and U. A. Wiedemann, arXiv:1101.0618.
  * (12) O. DeWolfe, S. S. Gubser, C. Rosen and D. Teaney, Prog. Part. Nucl. Phys. 75, 86 (2014).
  * (13) C. P. Herzog, A. Karch, P. Kovtun, C. Kozcaz and L. G. Yafe, JHEP 07 (2006) 013.
  * (14) S. S. Gubser, Phys. Rev. D 74, 126005 (2006).
  * (15) H. Liu, K. Rajagopal and U. A. Wiedemann, Phys. Rev. Lett. 97, 182301 (2006).
  * (16) H. Liu, K. Rajagopal and U. A. Wiedemann, JHEP 03 (2007) 066.
  * (17) S. S. Gubser, D. R. Gulotta, S. S. Pufu and F. D. Rocha, JEHP 10 (2008) 052.
  * (18) P. M. Chesler, K. Jensen and A. Karch, Phys. Rev. D 79, 025021 (2009).
  * (19) P. M. Chesler, K. Jensen, A. Karch and L. G. Yaffe, Phys. Rev. D 79, 125015 (2009).
  * (20) P. Arnold and D. Vaman, JHEP 10 (2010) 099.
  * (21) P. Arnold and D. Vaman, JHEP 04 (2011) 027.
  * (22) A. Ficnar and S. S. Gubser, Phys. Rev. D 89, 026002 (2014).
  * (23) A. Ficnar, S. S. Gubser and M. Gyulassy, Phys. Lett. B 738, 464 (2014).
  * (24) J. Erlich, E. Katz, D. T. Son and M. A. Stephanov, Phys. Rev. Lett. 95, 261602 (2005).
  * (25) J. Polchinski, M. J. Strassler, JHEP 05 (2003) 012.
  * (26) A. Karch, E. Katz, D. T. Son and M. A. Stephanov, Phys. Rev. D 74, 015005 (2006).
  * (27) J. P. Shock, F. Wu, Y-L. Wu and Z-F. Xie, JHEP 03 (2007) 064.
  * (28) A. Stoffers and I. Zahed, Phys. Rev. D 83, (2011) 055016.
  * (29) D. n. Li and M. Huang, JHEP 11 (2013) 088.
  * (30) D. n. Li, S. He, M. Huang and Q. S. Yan, JHEP 09 (2011) 041.
  * (31) S. He, M. Huang and Q. S. Yan, Phys. Rev. D 83, 045034 (2011).
  * (32) S. He, S. Y. Wu, Y. Yang and P. H. Yuan, JHEP 04 (2013) 093.
  * (33) R. Rougemont, A. Ficnar, S. Finazzo and J. Noronha, JHEP 04 (2016) 102\.
  * (34) P. Colangelo, F. Giannuzzi and S. Nicotri, Phys. Rev. D, 83, 035015 (2011).
  * (35) C. Park, D.-Y. Gwak, B.-H. Lee, Y. Ko and S. Shin, Phys. Rev. D 84, (2011) 046007.
  * (36) P. Colangelo, F. Giannuzzi and S. Nicotri, JHEP 05 (2012) 076.
  * (37) P. Colangelo, F. Giannuzzi, S. Nicotri and F. Zuo, Phys. Rev. D 88 (2013) 115011.
  * (38) Y. Xiong, X. Tang, and Z. Luo, Chin. Phys. C 43, 113103 (2019).
  * (39) X. R. Zhu and Z.-q. Zhang, Chin. Phys. C 44, 105105 (2020).
  * (40) X. Chen, S.-Q. Feng, Y.-F. Shi, Y. Zhong, Phys. Rev. D 97, 066015 (2018).
  * (41) Z.-q. Zhang and X. R. Zhu, Phys. Lett. B 793 (2019) 200.
  * (42) O. Andreev and V. I. Zakharov, Phys. Rev. D 74, 025023 (2006).
  * (43) M. Cvetic, et al., Nucl. Phys. B 558 (1999) 96.
  * (44) D. T. Son and A. O. Starinets, JHEP 03 (2006) 052.
  * (45) H. Liu, K. Rajagopal and Y. Shi, JHEP 08 (2008) 048.
  * (46) E. Caceres, A. Kundu and D. L. Yang, JHEP 03 (2014) 073