# On User Association in Large-Scale Heterogeneous LEO Satellite Network

Yuan Guo, Christodoulos Skouroumounis, Symeon Chatzinotas, and Ioannis
Krikidis This work received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(Grant agreement No. 819819). This work was also supported by the European
Union HORIZON programme under iSEE-6G (Grant agreement No. 101139291), and by
the Luxembourg National Research Fund (FNR) under the project MegaLEO
(C20/IS/14767486). Yuan Guo, Christodoulos Skouroumounis, and Ioannis Krikidis
are with the IRIDA Research Centre for Communication Technologies, Department
of Electrical and Computer Engineering, University of Cyprus, e-mail:
{yguo0001, cskour03<EMAIL_ADDRESS>Symeon Chatzinotas is with the
Interdisciplinary Centre for Security, Reliability and Trust (SnT), University
of Luxembourg, e-mail<EMAIL_ADDRESS>Part of this work was
presented at the IEEE International Conference on Communications, Rome, Italy,
May 2023 [1].

###### Abstract

In this paper, we investigate the performance of large-scale heterogeneous low
Earth orbit (LEO) satellite networks in the context of three association
schemes. In contrast to existing studies, where single-tier LEO satellite-
based network deployments are considered, the developed framework captures the
heterogeneous nature of real-world satellite network deployments. More
specifically, we propose an analytical framework to evaluate the performance
of multi-tier LEO satellite-based networks, where the locations of LEO
satellites are approximated as points of independent Poisson point processes,
with different density, transmit power, and altitude. We propose three
association schemes for the considered network topology based on: 1) the
Euclidean distance, 2) the average received power, and 3) a random selection.
By using stochastic geometry tools, analytical expressions for the association
probability, the downlink coverage probability, as well as the spectral
efficiency are derived for each association scheme, where the interference is
considered. Moreover, we assess the achieved network performance under several
different fading environments, including low, typical, and severe fading
conditions, namely non-fading, shadowed-Rician and Rayleigh fading channels,
respectively. Our results reveal the impact of fading channels on the coverage
probability, and illustrate that the average power-based association scheme
outperforms in terms of achieved coverage and spectral efficiency performance
against the other two association policies. Furthermore, we highlight the
impact of the proposed association schemes and the network topology on the
optimal number of LEO satellites, providing guidance for the planning of
multi-tier LEO satellite-based networks in order to enhance network
performance.

###### Index Terms:

Association scheme, coverage performance, heterogeneous networks, LEO
satellites, spectral efficiency, stochastic geometry.

##  I Introduction

The rapid evolution of massive Internet-of-Thing (IoT) devices demands the
emergence of sixth-generation (6G) networks [2]. As we transition from fifth-
generation (5G) to 6G, satellite communication networks are anticipated to
play a pivotal role in realizing the vision of a fully connected and
intelligent world. Satellite communications can significantly contribute to
achieving the goals of 6G by providing seamless global coverage and enhanced
connectivity in remote, rural, and underserved areas, as well as reinforcing
the resilience of communication networks in disaster-stricken regions [3, 4,
5]. Moreover, satellite communication systems hold an immense potential for
enabling advanced applications such as high-speed internet, ultra-low latency
communication, and massive IoT deployments, empowering sectors like smart
cities, agriculture, and transportation [6]. Hence, the concept of satellite
communication networks has emerged as a critical research paradigm, requiring
a comprehensive understanding of its potential impact, challenges, and
opportunities [7, 8, 9].

Satellite communication networks are characterized by a diverse range of
orbital configurations, such as geostationary (GEO), medium Earth orbit (MEO),
and low Earth orbit (LEO) satellites. Owing to the unique benefits of LEO
satellites, such as minimal latency, superior spatial flexibility, and low-
cost deployment, the employment of LEO satellites has emerged as an attractive
solution for 6G networks [10, 11, 12]. However, the proliferation of LEO
satellites and mega constellations pose significant challenges on network
modelling, design and analysis, necessitating the employment of sophisticated
mathematical tools that accurately capture the characteristics of large-scale
LEO satellite networks. Recently, tools from the field of stochastic geometry
have been leveraged to analyze the performance of large-scale LEO satellite-
based communication systems, highlighting its effectiveness as a powerful and
tractable mathematical tool for assessing the impact of key design parameters
on network performance [13, 14, 15, 16, 17]. Specifically, the authors in [13]
study the performance of large-scale LEO satellite communication networks,
illustrating the impact of satellites’ density and constellation altitude on
the coverage and communication latency. Moreover, the authors in [14]
investigate the coverage performance of LEO satellite networks, by assuming
that the satellite spatial deployment follows a homogeneous Poisson point
process (PPP), where a novel path-loss model is introduced for the satellite
communication networks. This work is further extended in [15], where the
optimal satellite altitude for attaining the highest coverage probability is
numerically investigated. Furthermore, the authors in [16] study the downlink
performance of LEO satellite-based networks by modelling the spatial
deployments of LEO satellites according to a homogeneous binomial point
process (BPP) and by investigating an iterative algorithm to maximize the
transmission rate and system throughput. The authors in [17] employ stochastic
geometry to analyze uplink performance of large-scale LEO satellite networks,
demonstrating the optimal LEO satellite density for achieving the maximum
ergodic capacity. Additionally, the authors in [18] adopt a non-homogeneous
PPP to model a massive LEO network and demonstrate the optimum altitude of LEO
satellites as well as number of orthogonal frequency channels for achieving
the highest network throughput. Nevertheless, these studies highlight the
capability of stochastic geometry in modelling large-scale LEO satellite-based
networks, under the assumption that LEO satellites are deployed at the same
altitude, opposed to real-world LEO satellite-based configurations [11, 13,
19]. Recently, the multi-altitude setup of LEO satellites has been discussed
in several works, by considering the same type of LEO satellites deployed at
various altitudes [19, 20, 21, 22]. In specific, the authors in [19] and [21]
model the multi-tier LEO satellites by using binomial point process and Cox
point process, respectively. In addition, a simulation-based study is
presented in [20] to evaluate the effect of interference on the downlink
systems. Moreover, the authors in [22] investigate the performance of multi-
hop routing in multi-tier LEO satellite networks. Although these works have
explored the altitude variations in LEO satellite-based networks, the real-
world multi-tier LEO satellites network exhibit more significant
heterogeneity, such as different transmit power and beamwidth of LEO
satellites, etc, which are overlooked from the current literature.

Given the heterogeneity presented in large-scale LEO satellite-based networks,
an appropriate association scheme, which methodically assigns each user to a
specific serving satellite, is of paramount importance for enhancing the end-
user performance. Various association schemes for conventional terrestrial
networks have been extensively explored in the literature, that are either
based on contact distance, cell area, received signal strength, or other
quality of service metrics [23, 24, 25]. In specific, the authors in [23]
propose two association schemes for millimeter wave networks, by taking into
consideration of the path-loss and the received signal power, respectively.
Additionally, a flexible association scheme is proposed in [24], where a
mobile user is associated with the strongest base stations (BSs) in terms of
long-term averaged biased power. The authors in [25] study a two-tier vertical
heterogeneous network consisting of terrestrial and aerial BSs, where an
average received power-based association is adopted. However, association
schemes for large-scale heterogeneous LEO satellite networks are missing from
the current literature; the majority of works that investigate massive LEO
satellite networks either ignore the heterogeneity or assume simplistic
distance-based association schemes [18, 14, 15, 16, 17]. This is mainly due to
the challenges arising from the different altitudes of multi-tier LEO
satellites that further escalate the heterogeneity of networks [1]. More
specifically, unlike conventional terrestrial networks, where the omission of
BSs’ heights does not deprive from the accuracy of analysis, the altitude
variation among multi-tier LEO satellites serves as a critical factor,
significantly affecting the coverage area and signal strength. The diverse
altitudes of LEO satellites not only lead to a higher level of complexity in
the design of appropriate association schemes but also pose substantial
challenges on analyzing the network performance.

Motivated by the above discussion, the aim of this work is to fill these gaps
by modelling and analyzing the performance of multi-tier LEO satellite-based
networks. Specifically, we extend our previous work [1] by developing
efficient association schemes to enhance the downlink performance of the
considered large-scale LEO satellite networks. The main contributions of this
paper are summarized as follows:

  * •

We develop an analytical framework based on stochastic geometry, which sheds
light on the modelling and analysis of large-scale heterogeneous LEO satellite
networks from a macroscopic point-of-view. The developed framework takes into
account the real-world heterogeneous nature of large-scale LEO satellite
networks, which involves the deployment of a general number of tiers of
satellites at various altitudes with different transmit power and density; the
locations of LEO satellites are modelled as the points of multiple independent
random point processes.

  * •

We propose three association schemes for the considered LEO satellite
networks, each with a distinct focus that collectively enhances network
performance. The first scheme aims at minimizing communication path-loss,
focusing on reducing the attenuation of signals to improve the quality of the
communication link. The second scheme is dedicated to maximizing received
signal strength, thereby ensuring robust communication by enhancing the power
of received signals. Finally, the third scheme addresses the computational
limitations of IoT devices by streamlining the association process, making
satellite network access more feasible for devices with restricted
computational capabilities.

  * •

We evaluate the performance of each association scheme, in terms of
association probability, coverage probability and spectral efficiency, by
considering different fading scenarios, namely, shadowed-Rician (SR), Rayleigh
and non-fading channels. By leveraging tools from stochastic geometry, we
initially evaluate the distance and the average signal power distributions,
based on which analytical expressions for the association probability with the
employment of each association schemes are derived. Additionally, we derive
the analytical expressions for the achieved coverage probability and spectral
efficiency in the context of proposed association schemes, with SR, Rayleigh
and non-fading channels.

  * •

Our results reveal the negative effect of fading channels on the coverage
probability for low signal-to-interference-plus-noise ratio (SINR) thresholds,
and their beneficial impact in high SINRs. In addition, the numerical results
illustrate that the average power-based association scheme outperforms the
distance-based association schemes in terms of coverage probability and
spectral efficiency, while the random selection-based scheme acts as a
performance lower bound. Finally, we highlight the trade-off imposed by the
number of LEO satellites between enhancing the coverage performance and
escalating the observed multi-user interference, while the optimal number of
LEO satellites is demonstrated for achieving the highest coverage probability.

The rest of the paper is organised as follows: Section II introduces the
considered system model. Section III presents our proposed association schemes
along with the detailed analysis of association probability. Section IV
presents the downlink network performance. Finally, numerical results are
presented in Section V, followed by our conclusions in Section VI.

Notation: $\Gamma(\cdot)$ and $\gamma(\cdot,\cdot)$ denote the Gamma and the
lower incomplete Gamma functions, respectively;$\
{}_{1}F_{1}\left(\cdot;\cdot;\cdot\right)$ is the confluent hypergeometric
function of the first kind; $(n)!$ denotes the factorial of $n$; $(m)_{z}$
denotes the Pochhammer symbol; $\Im\\{\cdot\\}$ denotes the imaginary
operator; $\mathcal{H}(\cdot)$ denotes the Heaviside function.

##  II System model

TABLE I: Summary of notations.

Notation | Description | Notation | Description
---|---|---|---
$\Phi_{k}$, $\lambda_{k}$ | PPP of the $k$-th tier LEO satellites with intensity $\lambda_{k}$ | $H_{k}$ | Altitude of the $k$-th tier LEO satellites
$x_{i,k}$ | Location of the $i$-th LEO satellite that belongs to the $k$-th tier | $\alpha$ | Path-loss exponent
$r_{i,k}$ | Distance between the typical gUE and the LEO satellite located at $x_{i,k}$ | $R_{\oplus}$ | Radius of the Earth
$G_{m,k}^{L},G_{s,k}^{L}$ | Main- and side-lobe antenna gain of the $k$-th tier LEO satellites | $G_{m}^{U},G_{s}^{U}$ | Main- and side-lobe antenna gain of the gUEs
$R_{\mathrm{max}_{k}}$ | Maximum distance from a gUE to the $k$-th tier satellites | $K$ | The number of tiers of LEO satellites
$(b,\Omega,m)$ | Parameters of the shadowed-Rician fading model | $P_{k}$ | Transmit power of the $k$-th tier LEO satellites
$N_{k}$ | Average number of the $k$-th tier LEO satellites | $\sigma^{2}$ | Thermal noise power

In this section, we provide details of the considered system model. The
network is studied from a large-scale perspective based on stochastic
geometry. The main mathematical notation related to the system model is
summarized in Table I.

###  II-A Network model

We consider a multi-tier constellation setup for a downlink LEO satellite-
based network as illustrated in Fig. 1, where LEO satellites are deployed at
$K$ different spherical surfaces concentric with the Earth, of altitudes
$H_{k}$ above the mean sea level, with $H_{k}<H_{k+1}$, $\forall\ {}1\leq
k\leq K-1$. For facilitating the analytical tractability, we assume that at
any certain instance, the locations and the orbits of LEO satellites are
uncorrelated; while within each tier constellation, the LEO satellites’
locations are assumed to be distributed according to a homogeneous PPP
$\Phi_{k}$ with intensity
$\lambda_{k}=\frac{N_{k}}{4\pi(R_{\oplus}+H_{k})^{2}}$, where $N_{k}$ is the
average number of LEO satellites within the $k$-th tier; while $\Phi_{i}$ and
$\Phi_{j}$ for $i\neq j$ are assumed to be independent [14, 15, 13].
Therefore, the average number of LEO satellites within each tier constellation
is given as $N_{k}=4\pi(R_{\oplus}+H_{k})^{2}\lambda_{k}$, where $R_{\oplus}$
is the Earth radius; while $\Phi_{i}$ and $\Phi_{j}$ for $i\neq j$ are assumed
to be independent [14, 15, 16]. Furthermore, we assume that the locations of
ground user equipments (gUEs) follow a uniform distribution. Without loss of
generality and based on the Slivnyak’s theorem, we perform our analysis from
the perspective of the typical gUE located at $(0,0,R_{\oplus})$ as
illustrated in Fig. 1, while the results hold for any gUE in the network [26].
Let $x_{i,k}$ denote the location of the $i$-th LEO satellite that belongs to
the $k$-th tier constellation, and $r_{i,k}$ represent the Euclidean distance
from the $i$-th LEO satellite to the typical gUE. Additionally, let $x_{0,k}$
with $1\leq k\leq K$ depict the locations of the closest LEO satellites to the
typical gUE from each tier constellation. Since a gUE can only receive signals
from the LEO satellites above its local horizon, the maximum distance between
a gUE and a visible LEO satellite that belongs to the $k$-th tier
constellation is given by $R_{{\rm
max}_{k}}=\sqrt{H_{k}^{2}+2H_{k}R_{\oplus}}$, where $1\leq k\leq K$ as
illustrated in Fig. 1 [27, 18].

Figure 1: Network topology of a $K$-tier LEO satellites system.

###  II-B Satellite-to-ground channel model

We assume that the satellite-to-ground (StG) wireless channels experience both
small-scale block fading and large-scale path-loss [14, 15, 16, 17, 18].
Regarding the small-scale fading, we adopt the well-known SR fading model for
each link, which is typically a Rician fading channel with fluctuating line-
of-sight (LoS) components [16]. Let $h_{i,j}$ denote the channel power gain of
the link between the typical gUE and the $i$-th LEO satellite that belongs in
the $j$-th tier. Then, the cumulative distribution function (CDF) of the
channel power gain is given by

$\displaystyle F_{h}(x)=$
$\displaystyle\left(\frac{2bm}{2bm+\Omega}\right)^{m}\sum_{z=0}^{\infty}\frac{(m)_{z}}{z!\Gamma(z+1)}\left(\frac{\Omega}{2bm+\Omega}\right)^{z}\gamma\left(z+1,\frac{1}{2b}x\right),$
(1)

where $\Omega$ and $2b$ are the average power for the LoS and the multi-path
components, respectively, and $m$ is the fading parameters based on the
Nakagami-$m$ fading channels [28]. Regarding the large-scale fading, we assume
an unbounded singular path-loss model. More specifically, the path-loss effect
of the StG channels only depends on the spatial distance between the LEO
satellites and the gUEs; the path-loss between the typical gUE and a LEO
satellite located at $x_{i,k}$ is given by $\ell(r_{i,k})=r_{i,k}^{-\alpha}$,
where $\alpha>2$ is the path-loss exponent [19]. Finally, all StG wireless
links exhibit additive white Gaussian noise with variance $\sigma^{2}$.

We consider that both gUEs and LEO satellites are equipped with multiple
antennas. For the sake of simplicity and tractability, we assume a widely-
adopted simple sectorized model for the antenna directionality of LEO
satellites and gUEs [23, 29]. The adopted sectorized antenna model
approximates the actual radiation pattern by a step function, where a constant
main lobe gain is considered over the beamwidth as well as a constant side
lobe gain is considered outside the beamwidth. This approach provides a
suitable approximation of the actual beam pattern, while maintaining adequate
accuracy for analysis purposes and does not deprive from the intuitions [23,
29, 19]. Specifically, the antenna array gain of the LEO satellites that
belongs in the $k$-th tier is given by
$G_{k}^{L}=\\{G^{L}_{m,k},G^{L}_{s,k}\\}$, where $G^{L}_{m,k}$ and
$G^{L}_{s,k}$ represent the main- and the side-lobe gains of LEO satellites’
antenna, respectively. Similarly, the antenna array gain of the gUE is
characterized by two discrete values, i.e. the main-lobe gain $G^{U}_{m}$ and
the side-lobe gain $G^{U}_{s}$. Furthermore, we assume that perfect beam
alignment is achieved between each gUE and its serving LEO satellite by
estimating the angle of arrival at both gUE and LEO satellite, such that the
intended link between each gUE and its associated LEO satellite lies in the
boresight direction of the antennas of both terminals. Additionally, we assume
the beam misalignment between each gUE and interfering LEO satellites.
Therefore, the antenna array gain of the link between a gUE and the $k$-th
tier LEO satellite is given by [29]

$G=\begin{cases}G^{L}_{m,k}G^{U}_{m},&\text{intended link},\\\
G^{L}_{s,k}G^{U}_{s},&\text{interfering link}.\end{cases}$ (2)

It is worth mentioning that, although interfering satellites may also be in
the boresight direction of the gUE’s antenna, the employment of such
assumption offers a tractable analysis for the considered large-scale LEO
satellite networks.

###  II-C Power adjusting mechanism

The altitudes of LEO satellites have a significant impact on the downlink
network performance, e.g. higher altitudes of LEO satellites typically result
in longer propagation distances of wireless signals, degrading the received
signal quality at the gUEs. Moreover, according to International
Telecommunication Union’s regulations, the transmit power of LEO satellites is
strictly constrained, such that the interference inflicted at existing GEO
satellites systems is lower than a predefined threshold [30, 12]. Motivated by
the above discussion, we assume that a power adjusting mechanism is employed
by the LEO satellites, where the transmit power of the LEO satellites of each
tier is adjusted according to their altitudes, while the maximum achievable
received signal power at a gUE from each tier LEO satellites are identical,
i.e.
$\frac{P_{k}G_{m,k}^{L}}{H_{k}^{\alpha}}=\frac{P_{j}G_{m,j}^{L}}{H_{j}^{\alpha}},\
{}\ {}\forall\ {}1\leq j,k\leq K$, where $P_{k}$ is the transmit power of the
$k$-th tier LEO satellites. It is worth mentioning that the adopted power
adjusting mechanism not only captures the heterogeneity of multi-tier LEO
satellite networks but also facilitates the analytical tractability of the
developed framework.

##  III Association schemes for Heterogeneous LEO satellite networks

We propose three association schemes, which are associated with a two-stage
procedure. In the first stage, $K$ LEO satellites are selected into the set of
candidate LEO satellites; while in the second stage, the serving LEO satellite
is selected among the set of candidate LEO satellites, according to the rules
of the adopted association policy, i.e the Euclidean distance, the average
received power over the fading channels, or the random selection111It should
be noted that LEO satellites are operated on orbits with high velocity, and
thus the potential handover process is required to keep associating the gUEs
with specific serving LEO satellites. For simplicity, we assume that a soft-
handover policy is employed to ensure the seamless connection between each gUE
and its serving LEO satellite [31].. The three association schemes are
presented and the analytical expression for the association probability for
each scheme is derived.

###  III-A Distance-based association (DbA) scheme

The DbA scheme is based on the Euclidean distance between the typical gUE and
the LEO satellites. In specific, the DbA scheme enables each gUE to select and
communicate with its nearest LEO satellite to maintain an acceptable received
signal quality by minimizing the path-loss. Therefore, in the first stage of
the DbA scheme, the nearest LEO satellites from each tier are selected into
the set of candidate LEO satellites, i.e.
$\mathcal{C}=\\{x_{0,1},\cdots,x_{0,k},\cdots,x_{0,K}\\}$. Then in the second
stage, the typical gUE communicate with the closest LEO satellite from the set
of candidate LEO satellites, e.g.
$x_{0,k}=\operatorname*{arg\,min}\nolimits_{x\in\mathcal{C}}\|x\|$ denotes the
case where the typical gUE is associated with a LEO satellite belonging in the
$k$-th tier. It is worth mentioning that the DbA association scheme is a low-
complexity scheme, which requires an a prior knowledge of the locations of LEO
satellites, which can be obtained by monitoring the location of the LEO
satellites through a low-rate feedback channel or by a global positioning
system (GPS) [32, 33]. Hence, the DbA scheme holds significant utility for
gUEs equipped with GPS capabilities.

As the DbA relies on the Euclidean distance, we initially assess the
distribution of the distance between the typical gUE and candidate LEO
satellites, while the probability density function (PDF) of $r_{0,k}$ for
$1\leq k\leq K$, is evaluated in the following lemma.

######  Lemma 1.

The PDF of the distance between the typical gUE and its $k$-th candidate LEO
satellite is given by

$\displaystyle
f_{D}(r,\lambda_{k},H_{k})=\frac{2\pi\lambda_{k}r(R_{\oplus}+H_{k})}{R_{\oplus}}\exp\left(-\frac{\pi\lambda_{k}(R_{\oplus}+H_{k})(r^{2}-H_{k}^{2})}{R_{\oplus}}\right),$
(3)

where $H_{k}\leq r<R_{\mathrm{max}_{k}}$.

###### Proof.

See Appendix A. ∎

It is worth emphasizing that while several previous studies have focused on
deriving the Earth-centered zenith angle distribution to study single-tier LEO
satellite-based networks222For a single-tier LEO satellite network, the
distance between a gUE and a LEO satellite can be uniquely determined by the
Earth-centered zenith angle [14, 15, 13, 19]., the distance distribution
presented in Lemma 1 consists of a more general result and can be used for
both single- and multi-tier LEO satellite networks. Specifically, by letting
$H_{k}=H$ and $\lambda_{k}=\lambda$, $\forall\ {}1\leq k\leq K$, the
corresponding PDF for the special scenario, namely, the single-tier deployment
is obtained.

Since each gUE can only communicate with the LEO satellites that are located
above the local horizon, the probability that at least one $k$-th tier LEO
satellite is located above the local horizon of the typical gUE is computed as

$\displaystyle
P^{\circ}(N_{k},H_{k})=1-\exp\left(-\frac{\pi\lambda_{k}(R_{\oplus}+H_{k})(R_{\mathrm{max}_{k}}^{2}-H_{k}^{2})}{R_{\oplus}}\right)=$
$\displaystyle
1-\exp\left(-\frac{N_{k}}{2\left({R_{\oplus}}/{H_{k}}+1\right)}\right).$ (4)

Note that $R_{\oplus}\approx 6371$ km and the altitude of the LEO satellites,
i.e. $H_{k}$, is typically from $500$ km to $2000$ km. Moreover, for the
considered large-scale LEO satellite networks, each tier consists of large
number of LEO satellites, i.e. $N_{k}\gg 100$. Hence, the exponential term
approaches to zero, i.e.
$\exp\left(-\frac{N_{k}}{2\left({R_{\oplus}}/{H_{k}}+1\right)}\right)\rightarrow
0$, and $P^{\circ}(N_{k},H_{k})\rightarrow 1$. Therefore, we exclude the
extreme scenario where no LEO satellite exists above the local horizon of the
gUEs throughout our analysis [14]. The accuracy of the adopted assumption is
validated by our numerical and simulation results presented in Section V.
Then, by applying the result from Lemma 1, we present the association
probability for the DbA scheme, denoted as $\mathcal{A}_{D,k}$ in the
following lemma.

######  Lemma 2.

The probability of the typical gUE to communicate with a LEO satellite from
the $k$-th tier, for the case where DbA scheme is adopted, is given by

$\mathcal{A}_{D,k}=\sum\limits_{i=k}^{K}\int\nolimits_{H_{i}}^{H_{i+1}}\frac{2\pi\lambda_{k}r(R_{\oplus}+H_{k})}{R_{\oplus}}\prod\limits_{j=1}^{i}\exp\left(-\frac{\pi\lambda_{j}(R_{\oplus}+H_{j})(r^{2}-H_{j}^{2})}{R_{\oplus}}\right)\mathrm{d}r,$
(5)

where $H_{K+1}=r_{{\rm max}_{K}}$.

###### Proof.

See Appendix B. ∎

The results derived in Lemma 2 provide valuable insights into the association
probability of the DbA scheme and the interplay between various design
parameters that affect the network performance. Moreover, it offers a more
convenient and efficient way to evaluate the association probability than the
Monte Carlo simulation, especially in the considered 3D large-scale multi-tier
LEO satellite networks with a large number of network tiers. In specific,
based on the analytical expression presented in Lemma 2, the association
probability can be immediately determined for any given network parameters,
including the number of tiers, the height of each LEO satellite tier, and the
density of LEO satellites, while the exponential product term in (5) shows the
interdependence among LEO satellite tiers. Additionally, by substituting
$\lambda_{k}$ with $\frac{N_{k}}{4\pi(R_{\oplus}+H_{k})}$, we have,
$\mathcal{A}_{D,k}=\sum\nolimits_{i=k}^{K}\int\nolimits_{H_{i}}^{H_{i+1}}\frac{N_{k}r}{2R_{\oplus}(R_{\oplus}+H_{k})}\prod\nolimits_{j=1}^{i}\exp\left(-\frac{\pi\lambda_{j}(R_{\oplus}+H_{j})(r^{2}-H_{j}^{2})}{R_{\oplus}}\right)\mathrm{d}r$,
based on which we can easily observe that the association probability of each
tier increases with the average number of LEO satellites and decreases with
the altitude of this tier.

###  III-B Power-based association (PbA) scheme

The second association scheme, namely PbA scheme, takes into consideration of
the average received signal power at the gUEs. In specific, the PbA scheme
associates each gUE with the LEO satellite that provides the highest average
signal power. Note that LEO satellites which belong in the same tier have
identical transmit power and antenna array gain, and thus the nearest one to
the gUE provides the strongest average received signal. Hence, the set of
candidate LEO satellites is identical with that formulated for the DbA scheme.
In the second stage, the serving LEO satellite is selected among the set of
candidate LEO satellites, which provides the highest average received signal
power instead of shortest Euclidean distance.

We now focus on the association probability for the PbA scheme. Let
$\xi_{0,k}$ represent the average received signal power from the LEO satellite
located at $x_{0,k}$. Note that $\xi_{0,k}$ is computed by taking the
expectation over the channel power gain, i.e.
$\xi_{0,k}=\bar{h}P_{k}G_{m,k}^{L}r_{0,k}^{-\alpha}G_{m}^{U}$, where $\bar{h}$
is the first-order moment of the channel power gain, i.e.
$\bar{h}=(2b+\Omega)$ for the SR fading channel model [28]. In the following
lemma, we provide the PDF of $\xi_{0,k}$, which is essential to assess the
association probability for the PbA scheme.

######  Lemma 3.

The PDF of the average received signal power at the typical gUE from its
$k$-th candidate LEO satellite is given by

$\displaystyle
f_{P}(\xi,\lambda_{k},H_{k})=\frac{2\pi\lambda_{k}(R_{\oplus}+H_{k})}{\xi\alpha
R_{\oplus}}\left(\frac{\varrho_{k}}{\xi}\right)^{\frac{2}{\alpha}}\exp\left(-\frac{\pi\lambda_{k}(R_{\oplus}+H_{k})}{R_{\oplus}}\left(\left(\frac{\varrho_{k}}{\xi}\right)^{\frac{2}{\alpha}}-H_{k}^{2}\right)\right),$
(6)

where $\varrho_{k}=\bar{h}P_{k}G^{L}_{m,k}G_{m}^{U}$,
$\xi_{\mathrm{min}_{k}}<\xi\leq\xi_{\mathrm{max}_{k}}$,
${\xi}_{\mathrm{min}_{k}}={\varrho_{k}}{R_{\mathrm{max}_{k}}^{-\alpha}}$ and
${\xi}_{\mathrm{max}_{k}}={\varrho_{k}}{H_{k}^{-\alpha}}$.

###### Proof.

See Appendix C. ∎

It is worth mentioning that, while the distance distribution was studied to
analyze the performance of the DbA scheme, the power distribution presented in
Lemma 3 offers a more appropriate method for assessing the performance
achieved by the PbA scheme. In particular, the average received signal power
from the candidate LEO satellites exhibits identical lower and upper bounds.
The minimum value of the average received signal power from any candidate LEO
satellite approaches zero due to the extremely long communication distance,
i.e., $\xi_{\mathrm{min}_{k}}\approx 0$, $\forall~{}1\leq k\leq K$; whereas
the maximum average received signal power from all candidate LEO satellites is
also identical, i.e., $\xi_{\mathrm{max}_{k}}=\xi_{\mathrm{max}_{j}}$,
$\forall~{}1\leq j,k\leq K$, owing to the power adjusting mechanism introduced
in Section II-C. In following lemma, we provide the association probability
for the PbA scheme.

######  Lemma 4.

The probability of the typical gUE to communicate with a LEO satellite from
the $k$-th tier, for the case where PbA scheme is adopted, is given by

$\mathcal{A}_{P,k}=\frac{2\pi\lambda_{k}\varrho_{k}^{\frac{2}{\alpha}}(R_{\oplus}+H_{k})}{\alpha
R_{\oplus}}\int\nolimits_{{\xi}_{\mathrm{min}_{k}}}^{{\xi}_{\mathrm{max}_{k}}}\left(\frac{1}{\xi}\right)^{\frac{2+\alpha}{\alpha}}\prod_{j=1}^{K}\exp\left(-\frac{\pi\lambda_{j}(R_{\oplus}+H_{j})}{R_{\oplus}}\left(\left(\frac{\varrho_{k}}{\xi}\right)^{\frac{2}{\alpha}}-H_{j}^{2}\right)\right)\mathrm{d}\xi.$
(7)

###### Proof.

See Appendix D. ∎

From the expression derived in above lemma, we can observe that the PbA scheme
considers one additional parameter compared to the DbA scheme, i.e. the
received signal power; whereas the analytical expression for the association
probability of the PbA scheme has a more concise form than DbA, due to the
employment of the power adjusting mechanism.

###  III-C Random selection-based association (RbA) scheme

The RbA scheme is based on a random tier selection. In specific, the first
stage process of the RbA scheme is the same with the DbA and PbA schemes. In
the second stage, the typical gUE randomly selects one serving LEO satellite
among the set of candidate LEO satellites. The RbA scheme does not require any
operations for comparing distance and/or received power in the second stage
process, leading to low implementation/signaling complexity. On one hand, the
RbA scheme can be adopted in scenarios with a large number of network tiers,
where its simplicity offers significant advantages. On the other hand, for
scenarios with a small number of network tiers, the RbA scheme exhibits almost
the same complexity as the DbA and the PbA schemes. Furthermore, as the RbA
enables a gUE to randomly select one LEO satellite among $K$ candidate LEO
satellites, the probability of a gUE being associated with each tier is equal,
i.e. $\mathcal{A}_{R,k}=\mathcal{A}_{R,j}=\frac{1}{K},~{}\forall\ {}1\leq
k,j\leq K$.

##  IV Downlink Performance of LEO satellite networks

We investigate the downlink performance of the considered heterogeneous LEO
satellite networks in the context of proposed association schemes.
Specifically, the downlink performance is initially evaluated in terms of the
coverage probability, corresponding to the probability that the instantaneous
SINR exceeds a predefined threshold. Additionally, we discuss the performance
of the proposed association schemes under another two distinct extreme
scenarios, namely, Rayleigh fading and non-fading environments. Finally, we
investigate the spectral efficiency achieved at typical gUEs with the
employment of each association scheme.

###  IV-A Coverage performance

The analytical expressions for the coverage probability achieved by each
association scheme with SR, Rayleigh, and non-fading channels are presented in
the following subsections.

####  IV-A1 Coverage probability under the DbA scheme

In order to facilitate the analysis of the coverage probability achieved by
the DbA scheme, we evaluate the distribution of the contact distance, i.e. the
distance between the typical gUE and its serving LEO satellite. We present the
conditional PDF of the contact distance for the DbA scheme in the following
lemma, by considering the case where the typical gUE is associated with the
$k$-th tier network.

######  Lemma 5.

When the DbA scheme is employed, the PDF of the contact distance is given by

$\displaystyle f_{D,k}(r|\mathcal{A}_{D,k})=$
$\displaystyle\frac{2\pi\lambda_{k}(R_{\oplus}+H_{k})}{R_{\oplus}\mathcal{A}_{D,k}}\sum_{i=k}^{K}\mathbbm{I}(H_{i},H_{i+1})\exp\left(-\sum_{j=1}^{i}\frac{\pi\lambda_{j}(R_{\oplus}+H_{j})(r^{2}-H_{j}^{2})}{R_{\oplus}}\right)r,$
(8)

where $\mathbb{I}(H_{i},H_{i+1})$ is the indicator function, given by

$\mathbb{I}(H_{i},H_{i+1})=\begin{cases}1,&\text{ $H_{i}\leq r<H_{i+1}$},\\\
0,&\text{\rm otherwise}.\end{cases}$ (9)

###### Proof.

See Appendix E. ∎

It is worth noting that the results presented in Lemma 5 are applicable to
various LEO satellite network topologies. For instance, by setting $H_{k}=H$,
$G_{m,k}^{L}=G_{m}^{L}$, $P_{k}=P$ and $\lambda_{k}=\lambda$ $\forall~{}1\leq
k\leq K$, we capture the single-tier LEO satellite networks as investigated in
previous studies [14, 15, 16].

By aiming to derive the closed-form expression for the coverage probability,
we approximate the power of the aggregate interference signals with its mean,
and the coverage probability achieved by the DbA scheme is evaluated in the
following theorem.

######  Theorem 1.

The coverage probability achieved by the DbA scheme is given by

$\mathcal{P}_{D}(\beta)=\sum\limits_{k=1}^{K}\mathcal{A}_{D,k}\mathcal{P}_{D,k}(\beta),$
(10)

where

$\displaystyle\mathcal{P}_{D,k}(\beta)\approx\int\nolimits_{H_{k}}^{R_{\mathrm{max}_{k}}}\left(1-F_{h}\left(\frac{\beta
r^{\alpha}\big{(}{\widetilde{\mathcal{I}}_{D}}(r)+\sigma^{2}\big{)}}{P_{k}G^{L}_{m,k}G^{U}_{m}}\right)\right)f_{D,k}(r|\mathcal{A}_{D,k})\mathrm{d}r,$
(11)

and

$\displaystyle{\widetilde{\mathcal{I}}_{D}}(r)=\sum\limits_{j=1}^{K}\frac{2\pi\lambda_{j}\bar{h}P_{j}G_{s,j}^{L}G_{s}^{U}(R_{\oplus}+H_{j})(\max\\{\min\\{r^{2-\alpha},H_{j}^{2-\alpha}\\}-R_{\mathrm{max}_{j}}^{2-\alpha},0\\})}{(\alpha-2)R_{\oplus}}.$
(12)

###### Proof.

See Appendix F. ∎

Note that the term $\widetilde{\mathcal{I}}_{D}(r)$ represents the conditional
mean of the aggregate interference, which is a function of random variable
$r$. Such methodology takes into account the randomness of spatial deployments
and fading channels, and avoids the intractability that typically arises when
performing system-level analysis with the SR fading channel model. The
accuracy of the adopted approximation is validated by our numerical and
simulation results in Section V. Given that $F_{h}(\cdot)$ is a monotonically
increasing function, we can infer, based on Theorem 1, that the shorter
contact distance, i.e. $r$, the higher transmit power, i.e. $P_{k}$, the
smaller interference i.e. $\widetilde{\mathcal{I}}_{D}(r)$, lead to a higher
converge probability. Moreover, the closed-form expression of
$\widetilde{\mathcal{I}}_{D}(r)$ provides a convenient approach for analyzing
the average interference levels in the large-scale LEO satellite networks.
Specifically, we can easily observe that the average interference received by
the typical gUE is directly proportional to the density and transmit power of
LEO satellites, while a greater path-loss exponent mitigates the interference
power.

####  IV-A2 Coverage probability under the PbA scheme

We evaluate the distribution of the average received signal power at the
typical gUE from its serving LEO satellite. By considering the scenario where
the typical gUE is associated with the $k$-th tier, the conditional PDF of
$\xi_{0,k}$ is evaluated in the following lemma.

######  Lemma 6.

When the PbA scheme is employed, the PDF of the average received signal power
at the typical gUE from its serving LEO satellite is given by

$\displaystyle
f_{P,k}(\xi|\mathcal{A}_{P,k})=\frac{2\pi\lambda_{k}\varrho_{k}^{\frac{2}{\alpha}}(R_{\oplus}+H_{k})}{\alpha
R_{\oplus}\mathcal{A}_{P,k}}\left(\frac{1}{\xi}\right)^{\frac{2}{\alpha}+1}\prod_{j=1}^{K}\exp\left(-\frac{\pi\lambda_{j}(R_{\oplus}+H_{j})}{R_{\oplus}}\left(\left(\frac{\varrho_{k}}{\xi}\right)^{\frac{2}{\alpha}}-H_{j}^{2}\right)\right).$
(13)

###### Proof.

See Appendix G. ∎

By utilizing the aforementioned results, we provide the analytical expression
for the coverage probability achieved by the PbA scheme, in the following
theorem.

######  Theorem 2.

The coverage probability achieved by the PbA scheme is given by

$\mathcal{P}_{P}(\beta)=\sum\limits_{k=1}^{K}\mathcal{A}_{P,k}\mathcal{P}_{P,k}(\beta),$
(14)

where

$\displaystyle\mathcal{P}_{P,k}(\beta)\approx\int\nolimits_{\xi_{\mathrm{min}_{k}}}^{\xi_{\mathrm{max}_{k}}}\left(1-F_{h}\left(\frac{\beta\bar{h}\big{(}\widetilde{\mathcal{I}}_{P}(\xi_{0,k})+\sigma^{2}\big{)}}{\xi_{0,k}}\right)\right)f_{P,k}(\xi_{0,k}|\mathcal{A}_{P,k})\mathrm{d}\xi_{0,k},$
(15)

and

$\displaystyle\widetilde{\mathcal{I}}_{P}(\xi_{0,k})=\sum\limits_{j=1}^{K}\frac{2\pi\lambda_{j}P_{j}G_{s,j}^{L}G_{s}^{U}\bar{h}(R_{\oplus}+H_{j})\big{(}\min\big{\\{}(\frac{\varrho_{k}}{\xi_{0,k}})^{\frac{2-\alpha}{\alpha}},H_{j}^{2-\alpha}\big{\\}}-R_{\mathrm{max}_{j}}^{2-\alpha}\big{)}}{(\alpha-2)R_{\oplus}}.$
(16)

###### Proof.

See Appendix H. ∎

We can easily observe from Theorem 2 that the interference level experienced
by the typical gUE using the PbA scheme is lower compared to that with the DbA
scheme. This observation can be justified by the fact that the PbA scheme
prioritizes the selection of LEO satellites based on the received signal
power. Therefore, the gUE is more likely to select a satellite with a stronger
signal, which in turn reduces the amount of interference from other LEO
satellites.

####  IV-A3 Coverage probability under the RbA scheme

By considering that the typical gUE is associated with a LEO satellite
belonging in the $k$-th tier, the PDF of the contact distance with the
employment of the RbA scheme is presented in the following lemma.

######  Lemma 7.

When the RbA scheme is employed, the PDF of contact distance is given by

$\displaystyle
f_{R,k}(r|\mathcal{A}_{R,k})=\frac{2\pi\lambda_{k}r(R_{\oplus}+H_{k})}{R_{\oplus}}\exp\left(-\frac{\pi\lambda_{k}(R_{\oplus}+H_{k})(r^{2}-H_{k}^{2})}{R_{\oplus}}\right).$
(17)

###### Proof.

The CDF of the contact distance by conditioning on the case that the typical
gUE is associated with the $k$-th tier LEO satellite is given by

$\displaystyle F_{R,k}(r|\mathcal{A}_{R,k})=1-\mathbb{P}\left\\{r_{0,k}\geq
r|\text{ the typical gUE is associated with the $k$-th tier}\right\\}.$ (18)

Since the tier selection is based on a random operation, i.e. it is
independent with $r_{0,k}$, we have
$F_{R,k}(r|\mathcal{A}_{R,k})=1-\mathbb{P}\left\\{r_{0,k}\geq r\right\\}$. By
directly following the results presented in Lemma 1, the final expression of
$f_{R,k}(r|\mathcal{A}_{R,k})$ is derived. ∎

We now present the analytical expression for the coverage probability achieved
by the RbA scheme in the following theorem.

######  Theorem 3.

The coverage probability achieved by the RbA scheme is given by

$\mathcal{P}_{R}(\beta)=\sum\limits_{k=1}^{K}\mathcal{A}_{R,k}\mathcal{P}_{R,k}(\beta),$
(19)

where

$\displaystyle\mathcal{P}_{R,k}(\beta)\approx\int\nolimits_{H_{k}}^{R_{\mathrm{max}_{k}}}\left(1-F_{h}\left(\frac{\beta
r_{0,k}^{\alpha}(\widetilde{\mathcal{I}}_{R}(r_{0,k})+\sigma^{2})}{P_{k}G^{L}_{m,k}G^{U}_{m}}\right)\right)f_{R,k}(r_{0,k}|\mathcal{A}_{R,k})\mathrm{d}r_{0,k},$
(20) and
$\displaystyle\widetilde{\mathcal{I}}_{R}(r_{0,k})=\frac{1}{(\alpha-2)R_{\oplus}}\bigg{(}{2\pi\lambda_{k}P_{k}G_{s,k}^{L}G_{s}^{U}\bar{h}(R_{\oplus}+H_{k})(r_{0,k}^{2-\alpha}-R_{\mathrm{max}_{k}}^{2-\alpha})}$
$\displaystyle\hskip 128.0374pt+\sum\limits_{j=1,j\neq
k}^{K}{2\pi\lambda_{j}P_{j}G_{s,j}^{L}G_{s}^{U}\bar{h}(R_{\oplus}+H_{j})(H_{j}^{2-\alpha}-R_{\mathrm{max}_{j}}^{2-\alpha})}\bigg{)}.$
(21)

###### Proof.

The proof follows similar steps presented in Appendix F, but with the average
interference power and association probability that correspond to the RbA
scheme. ∎

Some key insights can be obtained from the expression presented in the above
theorem. Although the RbA scheme has a low-complexity implementation, it
results in a higher level of interference experienced by the gUEs compared to
the other two schemes. Specifically, the random tier selection results in some
nearby LEO satellites from other tiers contributing significant interference
to the typical gUE.

####  IV-A4 Special fading scenarios

By aiming to evaluate the effects of fading channels on the heterogeneous LEO
satellites networks, we investigate the achieved performance in another two
fading scenarios, namely, Rayleigh and non-fading channels. Rayleigh fading
model is used to capture the worst scenario with weak or obstructed LoS links,
due to the near-ground blockages [14]. Although the Rayleigh fading channel
has been discussed in [27] for the single-tier LEO satellite-based networks
with simplistic distance-based association scheme, we now evaluate the
coverage probability with Rayleigh fading channels in terms of each proposed
association scheme. The exact analytical expression for the coverage
probability achieved by the each association scheme in Rayleigh fading
scenario is provided in following Proposition.

######  Proposition 1.

The coverage probabilities achieved by the DbA, PbA and RbA schemes with
Rayleigh fading channels are given by

$\mathcal{P}^{\dagger}_{\Delta}(\beta)=\sum_{k=1}^{K}\mathcal{A}_{\Delta,k}\mathcal{P}^{\dagger}_{\Delta,k}(\beta),$
(22)

where $\Delta\in\\{D,P,R\\}$ stands for the DbA, PbA and RbA schemes,
respectively, $\mathcal{P}^{\dagger}_{\Delta,k}(\beta)$ is given by

$\mathcal{P}^{\dagger}_{\Delta,k}(\beta)=\begin{cases}\int_{H_{k}}^{R_{\mathrm{max}_{k}}}\mathcal{L}^{\dagger}_{\mathcal{I}_{D}}(\varphi_{D})\exp(-\varphi_{D}\sigma^{2})f_{D,k}(r_{0,k}|\mathcal{A}_{D,k})\mathrm{d}r_{0,k},&\textit{\bf
if $\Delta=D$ },\\\
\int_{\xi_{\mathrm{min}_{k}}}^{\xi_{\mathrm{max}_{k}}}\mathcal{L}^{\dagger}_{\mathcal{I}_{P}}(\varphi_{P})\exp(-\varphi_{P}\sigma^{2})f_{P,k}(\xi_{0,k}|\mathcal{A}_{P,k})\mathrm{d}\xi_{0,k},&\textit{\bf
if $\Delta=P$ },\\\
\int_{H_{k}}^{R_{\mathrm{max}_{k}}}\mathcal{L}^{\dagger}_{\mathcal{I}_{R}}(\varphi_{R})\exp(-\varphi_{R}\sigma^{2})f_{R,k}(r_{0,k}|\mathcal{A}_{R,k})\mathrm{d}r_{0,k},&\textit{\bf
if $\Delta=R$ },\end{cases}$ (23)

$\mathcal{L}^{\dagger}_{\mathcal{I}_{\Delta}}(\varphi_{\Delta})$ is given by

$\mathcal{L}^{\dagger}_{\mathcal{I}_{\Delta}}(\varphi_{\Delta})=\begin{cases}\prod\limits_{j=1}^{K}\exp\left(-\frac{2\pi\lambda_{j}P_{j}G_{s,j}^{L}G_{s}^{U}\varphi_{D}(R_{\oplus}+H_{j})}{{R_{\oplus}}(\alpha-2)/(\varpi_{j}(\max\\{r_{0,k},H_{j}\\},\varphi_{D})-\varpi_{j}(R_{\mathrm{max}_{j}},\varphi_{D}))}\right),&\textit{\bf
if $\Delta=D$},\\\
\prod\limits_{j=1}^{K}\exp\Big{(}-\frac{2\pi\lambda_{j}P_{j}G_{s,j}^{L}G_{s}^{U}\varphi_{P}(R_{\oplus}+H_{j})}{{R_{\oplus}}(\alpha-2)/(\varpi_{j}(\max\\{({\varrho_{k}}/{\xi_{0,k}})^{1/\alpha},H_{j}\\},\varphi_{P})-\varpi_{j}(R_{\mathrm{max}_{j}},\varphi_{P}))}\Big{)},&\textit{\bf
if $\Delta=P$ },\\\
\exp\left(-\frac{2\pi\lambda_{k}P_{k}G_{s,k}^{L}G_{s}^{U}\varphi_{R}(R_{\oplus}+H_{k})}{{R_{\oplus}}(\alpha-2)/(\varpi_{k}(r_{0,k},\varphi_{R})-\varpi_{k}(R_{\mathrm{max}_{k}},\varphi_{R}))}\right)\times\\\
\hskip 56.9055pt\prod\limits_{j=1,j\neq
k}^{K}\exp\left(-\frac{2\pi\lambda_{j}P_{j}G_{s,j}^{L}G_{s}^{U}\varphi_{R}(R_{\oplus}+H_{j})}{{R_{\oplus}}(\alpha-2)/(\varpi_{j}(H_{j},\varphi_{R})-\varpi_{j}(R_{\mathrm{max}_{j}},\varphi_{R}))}\right),&\textit{\bf
if $\Delta=R$},\end{cases}$ (24)

$\varpi_{j}(r,\varphi)={r}^{2-\alpha}\
{}_{2}F_{1}\left[1,\frac{\alpha-2}{\alpha};2-\frac{2}{\alpha};-P_{j}G_{s,j}^{L}G_{s}^{U}r^{-\alpha}\varphi\right]$;
$\varphi_{D}=\varphi_{R}=\frac{\beta}{P_{k}G_{m,k}^{L}G_{m}^{U}r_{0,k}^{-\alpha}}$,
$\varphi_{P}=\frac{\beta}{\xi_{0,k}}$.

###### Proof.

See Appendix I. ∎

In contrast to the previous discussion with the Rayleigh fading channel, in
the following proposition, we investigate another extreme case, i.e. non-
fading scenario. In this case, the channel remains constant, i.e. $h=1$; this
scenario provides a useful benchmark to evaluate the impact of the fading
effect on the performance of LEO satellites networks.

######  Proposition 2.

The coverage probabilities achieved by the DbA, PbA and RbA schemes with non-
fading channels are given by

$\displaystyle\mathcal{P}^{\ddagger}_{\Delta}(\beta)=\sum_{k=1}^{K}\mathcal{A}_{\Delta,k}\mathcal{P}^{\ddagger}_{\Delta,k}(\beta),$
(25)

where

$\mathcal{P}^{\ddagger}_{\Delta,k}(\beta)=\begin{cases}\int\nolimits_{H_{k}}^{R_{\mathrm{max}_{k}}}\Big{(}\frac{1}{2}-\frac{1}{\pi}\int\nolimits_{0}^{\infty}\frac{\Im\left\\{\psi(t,\mathcal{S}_{D})\exp(it\beta\sigma^{2})\right\\}}{t}\mathrm{d}t\Big{)}f_{D,k}(r_{0,k}|\mathcal{A}_{D,k})\mathrm{d}r_{0,k},&\textit{\bf
if $\Delta=D$},\\\
\int\nolimits_{\xi_{\mathrm{min}_{k}}}^{\xi_{\mathrm{max}_{k}}}\Big{(}\frac{1}{2}-\frac{1}{\pi}\int\nolimits_{0}^{\infty}\frac{\Im\left\\{\psi(t,\mathcal{S}_{P})\exp(it\beta\sigma^{2})\right\\}}{t}\mathrm{d}t\Big{)}f_{P,k}(\xi_{0,k}|\mathcal{A}_{P,k})\mathrm{d}\xi_{0,k},&\textit{\bf
if $\Delta=P$},\\\
\int\nolimits_{H_{k}}^{R_{\mathrm{max}_{k}}}\Big{(}\frac{1}{2}-\frac{1}{\pi}\int\nolimits_{0}^{\infty}\frac{\Im\left\\{\psi(t,\mathcal{S}_{R})\exp(it\beta\sigma^{2})\right\\}}{t}\mathrm{d}t\Big{)}f_{R,k}(r_{0,k}|\mathcal{A}_{R,k})\mathrm{d}r_{0,k},&\textit{\bf
if $\Delta=R$},\end{cases}$ (26)

$i=\sqrt{-1}$,
$\psi(t,\mathcal{S}_{\Delta})=\exp(-it\mathcal{S}_{\Delta})\mathcal{L}^{\ddagger}_{\mathcal{I}_{\Delta}}(-it\beta)$,
$\mathcal{S}_{D}=\mathcal{S}_{R}=P_{k}G_{m,k}^{L}G_{m}^{U}r_{0,k}^{-\alpha}$,
$\mathcal{S}_{P}=\xi_{0,k}$,
$\mathcal{L}^{\ddagger}_{\mathcal{I}_{\Delta}}(z)$ is given by

$\mathcal{L}^{\ddagger}_{\mathcal{I}_{\Delta}}(z)=\begin{cases}\prod\limits_{j=1}^{K}\exp\Big{(}\frac{\pi\lambda_{j}(R_{\oplus}+H_{j})}{{R_{\oplus}}\alpha/(\vartheta_{j}(\max\\{r_{0,k},H_{j}\\},z)-\vartheta_{j}(R_{\mathrm{max}_{j}},z))}\Big{)},&\textit{\bf
if $\Delta=D$},\\\
\prod\limits_{j=1}^{K}\exp\Big{(}\frac{\pi\lambda_{j}(R_{\oplus}+H_{j})}{{R_{\oplus}}\alpha/(\vartheta_{j}\big{(}\max\\{({\varrho_{k}}/{\xi_{0,k}})^{1/\alpha},H_{j}\\},z\big{)}-\vartheta_{j}(R_{\mathrm{max}_{j}},z))}\Big{)},&\textit{\bf
if $\Delta=P$},\\\
\exp\left(\frac{\pi\lambda_{k}(R_{\oplus}+H_{k})}{{R_{\oplus}}\alpha/(\vartheta_{k}(r_{0,k},z)-\vartheta_{k}(R_{\mathrm{max}_{k}},z))}\right)\prod\limits_{{j=1,j\neq
k}}^{K}\exp\Big{(}\frac{\pi\lambda_{j}(R_{\oplus}+H_{j})}{{R_{\oplus}}\alpha/(\vartheta_{j}(H_{j},z)-\vartheta_{j}(R_{\mathrm{max}_{j}},z))}\Big{)},&\textit{\bf
if $\Delta=R$},\end{cases}$ (27)

$\vartheta_{j}(r,z)=P_{j}G_{s,j}^{L}G_{s}^{U}\left(r^{2}\left(\alpha-2\mathcal{E}\left(\frac{2+\alpha}{\alpha},r^{-\alpha}z\right)\right)+2z^{\frac{2}{\alpha}}\Gamma\left[-\frac{2}{\alpha}\right]\right)$,
and $\mathcal{E}(a,b)=\int_{1}^{\infty}\frac{e^{-bx}}{x^{a}}\mathrm{d}x$.

###### Proof.

See Appendix J. ∎

###  IV-B Spectral efficiency

We assess the spectral efficiency achieved by each gUE. Specifically, the
spectral efficiency is measured in terms of the average ergodic rate per unit
bandwidth [24, 34]; for simplicity, the average ergodic rate is computed in
unit of nats/s/Hz. By following a similar approach with the derivation of the
coverage probability, the average ergodic rate achieved at the typical gUE is
given by

${\mathcal{R}}_{\Delta}=\sum_{k=1}^{K}\mathcal{R}_{\Delta,k}\mathcal{A}_{\Delta,k},$
(28)

where $\mathcal{R}_{\Delta,k}$ denotes the average ergodic rate per unit
bandwidth of a typical gUE associated with the $k$-tier. By considering that
the typical gUE is associated with the $k$-th tier, an analytical expression
for $\mathcal{R}_{\Delta,k}$ is presented in the following proposition.

######  Proposition 3.

The average ergodic rate per unit bandwidth of the typical gUE associated with
the $k$-th tier is given by

$\displaystyle\mathcal{R}_{\Delta,k}=\int_{0}^{\infty}\frac{\mathcal{P}_{\Delta,k}(\beta)}{1+\beta}\mathrm{d}\beta,$
(29)

where $\Delta\in\\{D,P,R\\}$ and $\mathcal{P}_{\Delta,k}(\beta)$ is given in
Theorem 1, 2& 3.

###### Proof.

The proof directly follows from the definition of the average ergodic rate and
the spectral efficiency [24, 34]. ∎

Hence, by combining (28) and (29), the spectral efficiency of the network
under each association scheme is obtained. Additionally, we can obtain the
spectral efficiency of the network under the Rayleigh fading and non-fading
scenarios, by replacing $\mathcal{P}_{\Delta,k}(\beta)$ with
$\mathcal{P}^{\dagger}_{\Delta,k}(\beta)$ and
$\mathcal{P}^{\ddagger}_{\Delta,k}(\beta)$, respectively.

TABLE II: SIMULATION PARAMETERS.

Parameters |  Values |  Parameters |  Values
---|---|---|---
Network tiers ($K$) |  $3$ |  Altitude of LEO satellites ($H_{1},H_{2},H_{3}$) |  $500$ km, $600$ km, $700$ km
Path-loss exponent ($\alpha$) |  $3$ |  Average number of LEO satellites ($N_{1},N_{2},N_{3}$) |  $500$, $1000$, $1500$
Radius of Earth ($R_{\oplus}$) |  $6371$ km |  Transmit power of LEO satellites ($P_{1},P_{2},P_{3}$) |  $32$ W, $55.3$ W, $87.8$ W
Thermal noise power ($\sigma^{2}$) |  $-81.4$ dBm |  Main-lobe gain of LEO satellites $(G_{m,1}^{L},G_{m,2}^{L},G_{m,3}^{L})$ |  $47$ dBi, $47$ dBi, $47$ dBi
Main-lobe gain of gUEs $\left(G_{m}^{U}\right)$ |  $10$ dBi |  Side-lobe gain of LEO satellites $(G_{s,1}^{L},G_{s,2}^{L},G_{s,3}^{L})$ |  $27$ dBi, $27$ dBi, $27$ dBi
Side-lobe gain of gUEs $\left(G_{s}^{U}\right)$ |  $0$ dBi |  Parameters of SR fading model $(b,\Omega,m)$ |  $(0.158,19.4,1.59)$

##  V Numerical and Simulation Results

We present analytical and simulation results to validate the accuracy of our
model and illustrate the performance of the proposed association schemes.
Unless otherwise stated, we use the parameters given in Table II. It is worth
mentioning that, the proposed analytical framework is generic, and the
selection for these parameter values is for the purpose of illustration [19].

(a) Coverage probability for DbA scheme versus the SINR threshold.

(b) Coverage probability for PbA scheme versus the SINR threshold.

(c) Coverage probability for RbA scheme versus the SINR threshold.

Figure 2: Coverage probability achieved by each association scheme, in the
Shadow Rician, the Rayleigh and the non-fading scenarios.

Fig. 2 illustrates the impact of fading channels on the coverage performance
of LEO satellite networks. In specific, Fig. 2 plots the coverage probability
achieved by the three proposed association schemes, i.e. DbA, PbA, and RbA, in
different fading scenarios that are denoted as “Non-fading”, “SR fading”, and
“Rayleigh fading” . For the low SINR regime, we observe that the absence of
fading results in a higher coverage probability in comparison to the cases
with fading, for all association schemes. However, as the SINR threshold
increases, the decrease in the performance of the non-fading scenario is
greater than the SR and Rayleigh fading, leading to the observation that non-
fading and SR fading achieves the worst and best performance at high SINRs,
respectively. This can be justified by the fact that fading channel leads to
the random attenuation of signals, including interference, which can
consequently enhance the SINR in certain operation scenarios. Additionally,
the non-fading scenario achieves a better performance at low SINR thresholds
due to its stable signal strength. Moreover, the SR fading channel with its
strong LoS component, corresponds to a higher signal quality than the Rayleigh
fading channel, resulting in a higher SINR. Finally, we can observe that for
each association scheme, the coverage performance achieved under the Rayleigh
and non-fading scenarios demonstrate a perfect match between simulation
(denoted as “Sim.”) and analysis (denoted as “Ana.”) results; the adopted
approximations for the analysis of the SR fading scenario only lead to minor
gaps between the simulation and analytical results at low SINR thresholds.

Figure 3: Coverage probability achieved by each association scheme versus the
SINR threshold for different LEO satellites’ altitudes.

Fig. 3 depicts the effect of LEO satellites’ altitudes on the coverage
performance. Specifically, Fig. 3 plots the coverage probability versus the
SINR threshold for different association schemes and different altitudes.
Initially, it can be observed that the coverage probability decreases with the
increase of the LEO satellites’ altitudes. This was expected since higher
altitudes of LEO satellites correspond to longer propagation distances between
each gUE and its serving satellites, which results in a lower SINR. Moreover,
we observe that the PbA scheme surpasses both the DbA and RbA schemes in terms
of the coverage probability. This can be justified by the fact that the PbA
scheme associates the gUE with the LEO satellite that delivers the highest
average signal power, thereby boosting the SINR. In contrast, the DbA scheme
associates the gUE with the nearest LEO satellite, that may lead to severe
interference due to the existence of nearby satellites with larger transmit
power, resulting in a degraded SINR at the gUE. Moreover, the RbA scheme
randomly associates the gUE with one of the candidate LEO satellites, leading
to the worst coverage probability. Finally, the difference in coverage
probability between the DbA and RbA schemes diminishes as the satellites’
altitude increases. This is due to the fact that as the altitude of the
satellites increases, the difference in path-loss between the nearest and
other candidate LEO satellites becomes less significant, and thus the
advantage of associating the gUE with the nearest LEO satellite in the DbA
scheme becomes negligible.

Figure 4: The coverage probability achieved by each association schemes versus
the average number of LEO satellites, where $\beta=10$ dB and
$N_{1}=N_{2}=N_{3}=N$.

Fig. 4 demonstrates the effect of the average number of LEO satellites on the
coverage performance. For the sake of simplicity, we consider that each tier
comprises an equal average number of LEO satellites, i.e.
$N_{1}=N_{2}=N_{3}=N$. Accordingly, Fig. 4 plots the coverage probability
versus the number of LEO satellites within each tier, for different
association schemes, with the SINR threshold set to $\beta=10$ dB. We can
clearly observe that for all association schemes, as the number of LEO
satellites increases, there is an initial increase in coverage probability,
which subsequently declines upon exceeding a specific threshold (denoted as
“Optimal number of LEO”). The observed enhancement in coverage probability can
be explained by the decreasing distance between the serving LEO satellite and
the gUE, which results in decreased path-loss and subsequently leads to an
elevated SINR at the gUE. However, when the number of LEO satellites exceeds a
certain limit, the increased number of interfering satellites leads to a
severe multi-user interference, and thus to a reduced coverage probability.
Another notable observation is that the PbA scheme consistently achieves the
highest coverage probability among the three association schemes.
Additionally, the DbA scheme outperforms the RbA scheme when the number of LEO
satellites is relatively small, while as the number of LEO satellites
increases, the RbA scheme surpasses the DbA scheme in terms of coverage
probability. This observation can be explained by the fact that, the RbA
scheme allows equal probability of association across tiers, increasing the
chance of associating a gUE to satellites with higher transmit power from
higher tiers, as opposed to DbA scheme which prioritizes proximity, resulting
in frequent associations with lower-tier satellites with lower transmit power.
Finally, we can observe that in the conventional single-tier network, the PbA
and the DbA result in the same coverage probability. While the multi-tier
network typically attains a higher coverage probability than the single-tier
network, it is worth noting that with a large number of LEO satellites, the
single-tier network can actually surpass it in coverage probability. This is
based on the fact that the multi-tier network offers more access options than
the single-tier network, which also leads to additional interference and thus
finally results in a lower coverage probability than the single-tier network
in scenario with dense deployments of LEO satellites.

Figure 5: The spectral efficiency versus the transmit power of LEO satellites,
where $P_{2}=\frac{H_{2}^{\alpha}G_{m,1}^{L}}{H_{1}^{\alpha}G_{m,2}^{L}}P_{1}$
and $P_{3}=\frac{H_{3}^{\alpha}G_{m,1}^{L}}{H_{1}^{\alpha}G_{m,3}^{L}}P_{1}$.

Finally, Fig. 5 shows the impact of LEO satellites’ transmit power and the
large-scale fading on the spectral efficiency of the considered networks. More
specifically, Fig. 5 plots the spectral efficiency versus the transmit power
of the first tier LEO satellites for different path-loss exponent. Firstly,
for all association schemes, i.e. the DbA, the PbA, and the RbA schemes, a
higher path-loss corresponds to a lower spectral efficiency, which can be
attributed to the reduced signal strength at the receiver. Secondly, the PbA
scheme consistently outperforms the DbA and RbA schemes, which is in line with
the remarks presented in Fig. 3 and Fig. 4. Finally, when the path-loss
exponent is small (e.g., $\alpha=2.5$), increasing the transmit power
initially leads to a higher data rate. However, the gUE experiences strong
interference from neighbouring satellites as a result of the reduced signal
attenuation, which counterbalances the benefits of higher transmit power,
causing the data rate to plateau and remain constant at a saturation point.

##  VI Conclusion

In this work, we present a comprehensive analysis of large-scale heterogeneous
LEO satellite networks based on stochastic geometry. We introduce three
association schemes and evaluate their performance in terms of association
probability, coverage probability, and spectral efficiency, while analytical
expressions are derived by leveraging tools from stochastic geometry. By
aiming to reveal the impact of fading channels on network performance, we
evaluate the coverage probability in different fading scenarios, i.e. SR,
Rayleigh and non-fading channel models. Our results highlight the critical
impact of fading channels, altitude, the number of satellites, and transmit
power on the coverage probability and spectral efficiency. Specifically, we
demonstrate that the PbA scheme consistently outperforms the other two
schemes, while the RbA scheme surpasses the DbA scheme as the number of
satellites increases; the optimal number of LEO satellites for achieving the
highest coverage probability is numerically demonstrated. We show that higher
path-loss results in reduced spectral efficiency, and increasing transmit
power initially enhances spectral efficiency which then remains constant due
to the interference from neighbouring satellites. Our work offers valuable
insights and guidance for the design and optimization of future large-scale
heterogeneous LEO satellite networks.

##  Appendix A Proof of Lemma 1

The proof follows a similar approach as [14]. We first compute the CDF of
$r_{0,k}$ based on the null probability of PPP [26], i.e.

$\displaystyle F_{D}(r,\lambda_{k},H_{k})=\mathbb{P}\\{r_{0,k}\leq
r\\}=1-\mathbb{P}\\{r_{0,k}>r\\}=1-\exp\left(-\lambda_{k}\Psi(r,H_{k})\right),$
(30)

where
$\Psi(r,H_{k})=\frac{\pi(r^{2}-H_{k}^{2})(R_{\oplus}+H_{k})}{R_{\oplus}}$
represents the area of the spherical dome as shown in Fig. 6. Finally, by
taking the derivative of $F_{D}(r,\lambda_{k},H_{k})$ with respect to $r$,
i.e. $f_{D}(r,\lambda_{k},H_{k})=\frac{\partial}{\partial
r}F_{D}(r,\lambda_{k},H_{k})$, the final expression in Lemma 1 can be derived.

Figure 6: The closest LEO satellite of the $k$-th tier constellation.

##  Appendix B Proof of Lemma 2

By employing the DbA scheme, the probability that the typical gUE is
associated with the $k$-th tier LEO satellite can be expressed as
$\mathcal{A}_{D,k}=\mathbb{P}\left\\{r_{0,k}<\min_{j,j\neq
k}r_{0,j}\right\\}$. Note that for the considered heterogeneous LEO satellites
networks, the orbit altitude is much smaller than the Earth radius, i.e.
$H_{k}\ll R_{\oplus}$, and thus it is easy to verify that $H_{1}\leq
H_{2}\leq\cdots\leq H_{K}\leq R_{\mathrm{max}_{1}}\leq
R_{\mathrm{max}_{2}}\leq\cdots\leq R_{\mathrm{max}_{K}}$. Hence, by applying
the law of total probability and by denoting $R_{{\rm max}_{K}}$ with
$H_{K+1}$, we can have

$\mathcal{A}_{D,k}=\sum_{i=k}^{K}\mathbb{P}\left\\{\min\limits_{{1\leq j\leq
K,j\neq k}}r_{0,j}>r_{0,k}~{}\&~{}H_{i}\leq r_{0,k}\leq H_{i+1}\right\\}.$
(31)

Upon observing that $r_{0,k}\leq\min\limits_{1\leq j\leq K-i}r_{0,i+j}$ is
valid for $H_{i}\leq r_{0,k}\leq H_{i+1}$, we can further simplify
$\mathcal{A}_{D,k}$ as

$\mathcal{A}_{D,k}=\sum_{i=k}^{K}\mathbb{P}\left\\{\min_{{1\leq j\leq i,j\neq
k}}r_{0,j}>r_{0,k}~{}\&~{}H_{i}\leq r_{0,k}\leq H_{i+1}\right\\}.$ (32)

Then, by noticing that $r_{0,j}$ and $r_{0,k}$, $\forall\ {}j\neq k$ are
mutually independent, each term of the probability within $\mathcal{A}_{D,k}$
can be evaluated as

$\displaystyle\mathbb{P}\left\\{\min_{{1\leq j\leq i,j\neq
k}}r_{0,j}>r_{0,k}~{}\&~{}H_{i}\leq r_{0,k}\leq H_{i+1}\right\\}$ (33)
$\displaystyle=$
$\displaystyle\int\nolimits_{H_{i}}^{H_{i+1}}{\prod_{j=1,j\neq
k}^{i}\exp\left(-\frac{\pi\lambda_{j}(R_{\oplus}+H_{j})(r^{2}-H_{j}^{2})}{R_{\oplus}}\right)}f_{D}(r,\lambda_{k},H_{k})\mathrm{d}r$
$\displaystyle=$
$\displaystyle\int\nolimits_{H_{i}}^{H_{i+1}}\prod_{j=1}^{i}\exp\left(-\frac{\pi\lambda_{j}(R_{\oplus}+H_{j})(r^{2}-H_{j}^{2})}{R_{\oplus}}\right)\frac{2\pi\lambda_{k}r(R_{\oplus}+H_{k})}{R_{\oplus}}\mathrm{d}r,$

where the last step is obtained by substituting the expression of
$f_{D}(r,\lambda_{k},H_{k})$ given in (3) and by cancelling out the common
factors of the numerator and denominator.

##  Appendix C Proof of Lemma 3

To begin with, we evaluate the CDF of $\xi_{0,k}$ as following

$\displaystyle F_{P}(\xi,\lambda_{k},H_{k})=$
$\displaystyle\mathbb{P}\\{\xi_{0,k}\leq\xi\\}=1-\mathbb{P}\left\\{r_{0,k}<(\bar{h}P_{k}G^{L}_{m,k}G^{U}_{m}/\xi)^{\frac{1}{\alpha}}\right\\}$
(34) $\displaystyle=$ $\displaystyle
1-F_{D}\big{(}(\bar{h}P_{k}G^{L}_{m,k}G^{U}_{m}/\xi)^{\frac{1}{\alpha}},\lambda_{k},H_{k}\big{)}.$

The PDF of $\xi_{0,k}$ is obtained by taking the derivative of
$F_{P}(\xi,\lambda_{k},H_{k})$ with respect to $\xi$, i.e.
$f_{P}(\xi,\lambda_{k},H_{k})=\frac{\partial}{\partial\xi}F_{P}(\xi,\lambda_{k},H_{k})$.

##  Appendix D Proof of Lemma 4

According to the procedure of the PbA scheme, the probability that the typical
gUE is associated with the $k$-th tier LEO satellite is formulated as
$\mathcal{A}_{P,k}=\mathbb{P}\left\\{\xi_{0,k}>\max\limits_{{1\leq j\leq
K,j\neq k}}\xi_{0,j}\right\\}$. The deployments of $K$ tiers of LEO satellites
follow $K$ independent PPPs, which indicates that $\xi_{0,j}$ for $1\leq j\leq
K$ are mutually independent between each other. Hence, based on the order
statistics [35], we have

$\displaystyle\mathcal{A}_{P,k}=$ $\displaystyle\prod_{j=1,j\neq
k}^{K}\mathbb{P}\left\\{\xi_{0,k}>\xi_{0,j}\right\\}=\mathbb{E}_{\xi_{0,k}}\left\\{\prod_{j=1,j\neq
k}^{K}\Big{(}1-F_{P}\left(\xi_{0,k},\lambda_{j},H_{j}\right)\Big{)}\right\\}$
$\displaystyle=$ $\displaystyle\mathbb{E}_{\xi_{0,k}}\left\\{\prod_{j=1,j\neq
k}^{K}\exp\left(-\pi\lambda_{j}\frac{R_{\oplus}+H_{j}}{R_{\oplus}}\left(\left(\frac{\bar{h}P_{j}G^{L}_{m,j}G_{m}^{U}}{\xi_{0,k}}\right)^{\frac{2}{\alpha}}-H_{j}^{2}\right)\right)\right\\}.$
(35)

Finally, by evaluating the expectation over $\xi_{0,k}$, the result in Lemma 4
is derived.

##  Appendix E Proof of Lemma 5

Without loss of generality, we consider that the typical gUE is associated
with the $k$-th tier network for $1\leq k\leq K$. Then, the CDF of the contact
distance is computed as

$\displaystyle F_{D,k}(r|\mathcal{A}_{D,k})=$
$\displaystyle\mathbb{P}\left\\{r_{0,k}\leq r\ {}\Big{|}\
{}r_{0,k}<\min_{\begin{subarray}{c}1\leq j\leq K,\\\ j\neq
k\end{subarray}}r_{0,j}\right\\}=1-\frac{\mathbb{P}\left\\{r_{0,k}>r~{}\&~{}r_{0,k}<\min\limits_{{1\leq
j\leq K,j\neq
k}}r_{0,j}\right\\}}{\mathbb{P}\left\\{r_{0,k}<\min\limits_{{1\leq j\leq
K,j\neq k}}r_{0,j}\right\\}}.$ (36)

Note that $\mathbb{P}\left\\{r_{0,k}<\min_{1\leq j\leq K,j\neq
k}r_{0,j}\right\\}=\mathcal{A}_{D,k}$. Then, for the case where $H_{i}\leq
r<H_{i+1}$, the probability term in (36) is computed as

$\displaystyle\mathbb{P}\left\\{r<r_{0,k}<\min\limits_{{1\leq j\leq K,j\neq
k}}r_{0,j}~{}\&~{}H_{i}\leq r<H_{i+1}\right\\}$
$\displaystyle=\int\nolimits_{r}^{H_{i+1}}\prod_{j=1}^{i}\exp\left(-\frac{\pi\lambda_{j}(R_{\oplus}+H_{j})(r^{2}-H_{j}^{2})}{R_{\oplus}}\right)\frac{2\pi\lambda_{k}r(R_{\oplus}+H_{k})}{R_{\oplus}}\mathrm{d}r.$
(37)

By taking into account all possible ranges of $r$ within $H_{k}\leq r\leq
R_{\mathrm{max}_{k}}$ and by substituting (E) into (36), we obtain the final
expression of $F_{D}(r|\mathcal{A}_{D,k})$, i.e.

$\displaystyle F_{D,k}(r|\mathcal{A}_{D,k})$ $\displaystyle=$ $\displaystyle
1-\frac{\sum_{i=k}^{K}\mathbbm{I}(H_{i}\leq
r<H_{i+1})\int\nolimits_{r}^{H_{i+1}}\prod_{j=1}^{i}\exp\left(-\frac{\pi\lambda_{j}(R_{\oplus}+H_{j})(r^{2}-H_{j}^{2})}{R_{\oplus}}\right)\frac{2\pi\lambda_{k}r(R_{\oplus}+H_{k})}{R_{\oplus}}\mathrm{d}r}{\mathcal{A}_{D,k}}.$
(38)

Finally, the PDF of the contact distance, i.e. $f_{D,k}(r|\mathcal{A}_{k})$,
can be derived by taking the derivative of $F_{D,k}(r|\mathcal{A}_{D,k})$ with
respect to $r$, i.e. $f_{D,k}(r|\mathcal{A}_{k})=\frac{\partial}{\partial
r}F_{D,k}(r|\mathcal{A}_{D,k})$.

##  Appendix F Proof of Theorem 1

By the law of the total probability and according to the definition of the
coverage probability, we have
$\mathcal{P}_{D}(\beta)=\sum_{k=1}^{K}\mathcal{P}_{D,k}(\beta)\mathcal{A}_{D,k}$,
where $\mathcal{P}_{D,k}(\beta)$ represents the coverage probability when the
typical gUE is associated with the $k$-th tier, which is given by

$\displaystyle\mathcal{P}_{D,k}(\beta)=$
$\displaystyle\mathbb{P}\left\\{\frac{P_{k}h_{0,k}G_{m,k}^{L}G^{U}_{m}r_{0,k}^{-\alpha}}{\sum\nolimits_{x_{i,j}\in\cup_{j=1}^{K}\Phi_{j}^{\backslash
x_{0,k}}}P_{j}h_{i,j}r_{i,j}^{-\alpha}G^{L}_{s,j}G^{U}_{s}+\sigma^{2}}\geq\beta\right\\}$
(39) $\displaystyle=$
$\displaystyle\mathbb{E}_{r_{0,k},h_{i,j}}\left\\{\mathbb{P}\left\\{h_{0,k}\geq\frac{\beta\big{(}\sum\nolimits_{x_{i,j}\in\cup_{j=1}^{K}\Phi_{j}^{\backslash
x_{0,k}}}P_{j}h_{i,j}r_{i,j}^{-\alpha}G^{L}_{s,j}G^{U}_{s}+\sigma^{2}\big{)}}{P_{k}G^{L}_{m,k}G^{U}_{m}r_{0,k}^{-\alpha}}\right\\}\right\\}$
$\displaystyle\approx$
$\displaystyle\mathbb{E}_{r_{0,k}}\left\\{\mathbb{P}\left\\{h_{0,k}\geq\frac{\beta\big{(}\widetilde{\mathcal{I}}_{D}(r_{0,k})+\sigma^{2}\big{)}}{P_{k}G^{L}_{m,k}G^{U}_{m}r_{0,k}^{-\alpha}}\right\\}\right\\}=\mathbb{E}_{r_{0,k}}\left\\{1-F_{h}\left(\frac{\beta\big{(}\widetilde{\mathcal{I}}_{D}(r_{0,k})+\sigma^{2}\big{)}}{P_{k}G^{L}_{m,k}G^{U}_{m}r_{0,k}^{-\alpha}}\right)\right\\},$

where $F_{h}(\cdot)$ is the CDF of channel power gain;
$\widetilde{\mathcal{I}}_{D}(r_{0,k})$ is the mean of the aggregate
interference power, and is computed as following

$\displaystyle\widetilde{\mathcal{I}}_{D}(r_{0,k})=$
$\displaystyle\mathbb{E}\left\\{\sum\nolimits_{x_{i,j}\in\cup_{j=1}^{K}\Phi_{j}^{\backslash
x_{0,k}}}P_{j}G^{L}_{s,j}h_{i,j}r_{i,j}^{-\alpha}G^{U}_{s}\right\\}$ (40)
$\displaystyle\overset{(a)}{=}$
$\displaystyle\sum_{j=1}^{K}\int\nolimits_{0}^{2\pi}\int\nolimits_{\theta}\bar{h}\lambda_{j}P_{j}G^{L}_{s,j}r_{i,j}^{-\alpha}G^{U}_{s}(R_{\oplus}+H_{j})^{2}\sin\theta\mathrm{d}\theta\mathrm{d}\phi$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{K}\int\nolimits_{\theta}2\pi\bar{h}\lambda_{j}P_{j}G^{L}_{s,j}r_{i,j}^{-\alpha}G^{U}_{s}(R_{\oplus}+H_{j})^{2}\sin\theta\mathrm{d}\theta,$

where $(a)$ follows from the Campbell’s Theorem of PPP with the spherical
coordinates [26, 14]. Then, by replacing $\theta$ with a function of $r$ based
on the law of cosines, i.e.

$\theta=\arccos\left(\frac{R_{\oplus}^{2}+(R_{\oplus}+H_{j})^{2}-r^{2}}{2R_{\oplus}(R_{\oplus}+H_{j})}\right),$
(41)

we have

$\displaystyle\mathrm{d}\theta=$
$\displaystyle\frac{r}{R_{\oplus}(H_{j}+R_{\oplus})\sqrt{1-\frac{\left((H_{j}+R_{\oplus})^{2}-r^{2}+R_{\oplus}^{2}\right)^{2}}{4R_{\oplus}^{2}(H_{j}+R_{\oplus})^{2}}}}\mathrm{d}r,$
(42) $\displaystyle\sin\theta=$
$\displaystyle\sqrt{1-\frac{\left((H_{j}+R_{\oplus})^{2}-r^{2}+R_{\oplus}^{2}\right)^{2}}{4R_{\oplus}^{2}(H_{j}+R_{\oplus})^{2}}}.$

Hence, we can further rewrite $\widetilde{\mathcal{I}}_{D}(r_{0,k})$ as
following

$\displaystyle\widetilde{\mathcal{I}}_{D}(r_{0,k})=$
$\displaystyle\sum_{j=1}^{K}\int\nolimits_{\max\\{r_{0,k},H_{j}\\}}^{R_{\mathrm{max}_{j}}}2\pi\bar{h}\lambda_{j}\frac{R_{\oplus}+H_{j}}{R_{\oplus}}P_{j}G^{L}_{s,j}r^{1-\alpha}G_{s}^{U}\mathrm{d}r.$
(43)

Finally, by solving the above integral and by evaluating the expectation over
$r_{0,k}$, the final results in Theorem 1 are derived.

##  Appendix G Proof of Lemma 6

Consider that the typical gUE is associated with the $k$-th tier, the CDF of
$\xi_{0,k}$ is computed as

$\displaystyle F_{P}(\xi|\mathcal{A}_{P,k})=$
$\displaystyle\mathbb{P}\left\\{\xi_{0,k}\leq\xi\Big{|}\xi_{0,k}>\max_{{1\leq
j\leq K,j\neq
k}}\xi_{0,j}\right\\}=\frac{\mathbb{P}\left\\{\max\limits_{{1\leq j\leq
K,j\neq
k}}\xi_{0,j}<\xi_{0,k}\leq\xi\right\\}}{\mathbb{P}\left\\{\xi_{0,k}>\max\limits_{{1\leq
j\leq K,j\neq k}}\xi_{0,j}\right\\}}.$ (44)

Note that the denominator of (44) is exactly the association probability
presented in Lemma 4, i.e.
$\mathbb{P}\left\\{\xi_{0,k}>\max\limits_{\begin{subarray}{c}1\leq j\leq
K,~{}j\neq k\end{subarray}}\xi_{0,j}\right\\}=\mathcal{A}_{P,k}$; while the
numerator can be further computed as

$\displaystyle\mathbb{P}\left\\{\max\limits_{\begin{subarray}{c}1\leq j\leq
K\\\ j\neq
k\end{subarray}}\xi_{0,j}<\xi_{0,k}\leq\xi\right\\}=\int_{\xi_{\mathrm{min}_{k}}}^{\xi}\prod_{\begin{subarray}{c}j=1,j\neq
k\end{subarray}}^{K}\mathbb{P}\left\\{\xi_{0,k}>\xi_{0,j}\right\\}f_{P}(\xi_{0,k},\lambda_{k},H_{k})\mathrm{d}\xi_{0,k}$
$\displaystyle=$
$\displaystyle\int_{\xi_{\mathrm{min}_{k}}}^{\xi}\prod_{\begin{subarray}{c}j=1,j\neq
k\end{subarray}}^{K}\big{(}1-F_{P}(\xi_{0,k},\lambda_{j},H_{j})\big{)}f_{P}(\xi_{0,k},\lambda_{k},H_{k})\mathrm{d}\xi_{0,k}$
(45) $\displaystyle=$
$\displaystyle\int_{\xi_{\mathrm{min}_{k}}}^{\xi}\prod_{j=1,j\neq
k}^{K}\exp\left(-\pi\lambda_{j}\frac{R_{\oplus}+H_{j}}{R_{\oplus}}\left(\left(\frac{\bar{h}P_{j}G^{L}_{m,j}G_{m}^{U}}{\xi_{0,k}}\right)^{\frac{2}{\alpha}}-H_{j}^{2}\right)\right)f_{P}(\xi_{0,k},\lambda_{k},H_{k})\mathrm{d}\xi_{0,k}.$

Then by substituting (G) into (44) and by taking the partial derivative of
$F_{P}(\xi|\mathcal{A}_{P,k})$ with respect to $\xi$, the conditional PDF of
$\xi_{0,k}$ is derived.

##  Appendix H Proof of Theorem 2

The proof follows a similar methodology presented in Appendix F. In specific,
based on the definition of the PbA scheme and by the law of the total
probability, the coverage probability is formulated as
$\mathcal{P}_{P}(\beta)=\sum_{k=1}^{K}\mathcal{P}_{P,k}(\beta)\mathcal{A}_{P,k}$,
where $\mathcal{P}_{P,k}(\beta)$ is given by

$\displaystyle\mathcal{P}_{P,k}(\beta)=$
$\displaystyle\mathbb{P}\left\\{\frac{P_{k}h_{0,k}G_{m,k}^{L}G^{U}_{m}r_{0,k}^{-\alpha}}{\sum\nolimits_{x_{i,j}\in\cup_{j=1}^{K}\Phi_{j}^{\backslash
x_{0,k}}}P_{j}h_{i,j}r_{i,j}^{-\alpha}G^{L}_{s,j}G^{U}_{s}+\sigma^{2}}\geq\beta\right\\}$
(46) $\displaystyle\approx$
$\displaystyle\mathbb{E}_{\xi_{0,k}}\left\\{\mathbb{P}\left\\{h_{0,k}\geq\frac{\beta\bar{h}\big{(}\widetilde{\mathcal{I}}_{P}(\xi_{0,k})+\sigma^{2}\big{)}}{\xi_{0,k}}\right\\}\right\\},$

where $\widetilde{\mathcal{I}}_{P}(\xi_{0,k})$ is the mean of the aggregate
interference signal power observed at the gUE with the employment of the PbA
scheme, which is computed as

$\displaystyle{\widetilde{\mathcal{I}}_{P}}(\xi_{0,k})=$
$\displaystyle\mathbb{E}\left\\{\sum\nolimits_{x_{i,j}\in\cup_{j=1}^{K}\Phi_{j}^{\backslash
x_{0,k}}}P_{j}h_{i,j}r_{i,j}^{-\alpha}G^{L}_{s,j}G^{U}_{s}\right\\}=\bar{h}\mathbb{E}\left\\{\sum\nolimits_{x_{i,j}\in\cup_{j=1}^{K}\Phi_{j}^{\backslash
x_{0,k}}}P_{j}r_{i,j}^{-\alpha}G^{L}_{s,j}G^{U}_{s}\right\\}$ (47)
$\displaystyle=$
$\displaystyle\sum_{j=1}^{K}\int\nolimits_{0}^{2\pi}\int\nolimits_{\theta}\bar{h}\lambda_{j}P_{j}G^{L}_{s,j}r_{i,j}^{-\alpha}G^{U}_{s}(R_{\oplus}+H_{j})^{2}\sin\theta\mathrm{d}\theta\mathrm{d}\phi$
$\displaystyle=$
$\displaystyle\sum\limits_{j=1}^{K}\int_{\max\left\\{\left({\varrho_{k}}/{\xi_{0,k}}\right)^{1/\alpha},H_{j}\right\\}}^{R_{\mathrm{max}_{j}}}2\pi\bar{h}\lambda_{j}\frac{R_{\oplus}+H_{j}}{R_{\oplus}}P_{j}G^{L}_{s,j}r^{1-\alpha}G_{s}^{U}\mathrm{d}r,$

where the last two steps directly follow from the methodology used for the
derivation of $\widetilde{\mathcal{I}}_{D}(r_{0,k})$ presented in Appendix F,
but with a different lower limit for the integral. Finally, by applying the
CDF of $h_{0,k}$, and by evaluating the average over $\xi_{0,k}$, the results
in Theorem 2 are proven.

##  Appendix I Proof of Proposition 1

We initially evaluate the coverage probability for DbA scheme under the
Rayleigh fading. Considering that the typical gUE is associated with the
$k$-th tier, then the coverage probability achieved with the DbA scheme is
computed as

$\displaystyle\mathcal{P}^{\dagger}_{D,k}(\beta)=$
$\displaystyle\mathbb{E}\left\\{\mathbb{P}\left\\{h_{0,k}\geq\frac{\beta\big{(}\sum\nolimits_{x_{i,j}\in\cup_{j=1}^{K}\Phi_{j}^{\backslash
x_{0,k}}}P_{j}h_{i,j}r_{i,j}^{-\alpha}G^{L}_{s,j}G^{U}_{s}+\sigma^{2}\big{)}}{P_{k}G^{L}_{m,k}G^{U}_{m}r_{0,k}^{-\alpha}}\right\\}\right\\}$
$\displaystyle=$
$\displaystyle\mathbb{E}\left\\{\exp\left(-\frac{\beta\big{(}\sum\nolimits_{x_{i,j}\in\cup_{j=1}^{K}\Phi_{j}^{\backslash
x_{0,k}}}P_{j}h_{i,j}r_{i,j}^{-\alpha}G^{L}_{s,j}G^{U}_{s}+\sigma^{2}\big{)}}{P_{k}G^{L}_{m,k}G^{U}_{m}r_{0,k}^{-\alpha}}\right)\right\\}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{r_{0,k}}\left\\{\mathbb{E}\left\\{\exp\left(-{\varphi_{D}\sum\nolimits_{x_{i,j}\in\cup_{j=1}^{K}\Phi_{j}^{\backslash
x_{0,k}}}P_{j}h_{i,j}r_{i,j}^{-\alpha}G^{L}_{s,j}G^{U}_{s}}\right)\right\\}\exp(-\varphi_{D}\sigma^{2})\right\\}$
$\displaystyle=$
$\displaystyle\int_{H_{k}}^{R_{\mathrm{max}_{k}}}\mathcal{L}^{\dagger}_{\mathcal{I}_{D}}(\varphi_{D})\exp(-\varphi_{D}\sigma^{2})f_{D,k}(r_{0,k}|\mathcal{A}_{D,k})\mathrm{d}r_{0,k},$
(48)

where the first step is based on the fact that for the Rayleigh fading model,
the channel power gain is an exponential random variable with mean one, i.e.
$h\sim\exp(1)$;
$\varphi_{D}=\frac{\beta}{P_{k}G_{m,k}^{L}G_{m}^{U}r_{0,k}^{-\alpha}}$ and
$\mathcal{L}^{\dagger}_{\mathcal{I}_{D}}(z)$ is the Laplace transform of the
interference, which is evaluated as following, i.e.

$\displaystyle\mathcal{L}_{\mathcal{I}_{D}}^{\dagger}(\varphi_{D})=$
$\displaystyle\mathbb{E}\left\\{\prod\nolimits_{x_{i,j}\in\cup_{j=1}^{K}\Phi_{j}^{\backslash
x_{0,k}}}\mathbb{E}_{h_{i,j}}\bigg{\\{}\exp\left(-\varphi_{D}P_{j}h_{i,j}r_{i,j}^{-\alpha}G^{L}_{s,j}G^{U}_{s}\right)\bigg{\\}}\right\\}$
(49) $\displaystyle=$
$\displaystyle\mathbb{E}\left\\{\prod\nolimits_{x_{i,j}\in\cup_{j=1}^{K}\Phi_{j}^{\backslash
x_{0,k}}}\exp\left(\frac{1}{1+\varphi_{D}P_{j}r_{i,j}^{-\alpha}G^{L}_{s,j}G^{U}_{s}}\right)\right\\}$
$\displaystyle=$
$\displaystyle\mathbb{E}\left\\{\prod_{j=1}^{K}\exp\left(-2\pi\lambda_{j}\frac{R_{\oplus}+H_{j}}{R_{\oplus}}\int_{\max\\{r_{0,k},H_{j}\\}}^{R_{\mathrm{max}_{j}}}\left(1-\frac{1}{1+\varphi_{D}P_{j}r^{-\alpha}G^{L}_{s,j}G^{U}_{s}}\right)r\mathrm{d}r\right)\right\\},$

where the last step is obtained by using the probability generating function
of the PPP [26], where the above integral can be solved based on [36,
3.194.5]. Finally, by applying the law of total probability, the coverage
probability of DbA scheme in the Rayleigh fading scenario is derived. Note
that the proof of the coverage probability of PbA and RbA schemes in the
Rayleigh fading scenario, follows a similar methodology and thus is omitted.

##  Appendix J Proof of Proposition 2

We provide the detailed proof of the coverage probability for the DbA scheme
in the non-fading scenario, while a similar approach can be easily applied to
the PbA and RbA schemes. Hence, by considering that the DbA scheme is employed
and the typical gUE is associated with the $k$-th tier, the coverage
probability is formulated as

$\displaystyle\mathcal{P}^{\ddagger}_{D,k}(\beta)=$
$\displaystyle\mathbb{P}\left\\{\frac{P_{k}G_{m,k}^{L}G^{U}_{m}r_{0,k}^{-\alpha}}{\sum\nolimits_{x_{i,j}\in\cup_{j=1}^{K}\Phi_{j}^{\backslash
x_{0,k}}}P_{j}r_{i,j}^{-\alpha}G^{L}_{s,j}G^{U}_{s}+\sigma^{2}}\geq\beta\right\\}$
(50) $\displaystyle=$
$\displaystyle\mathbb{P}\left\\{P_{k}G_{m,k}^{L}G^{U}_{m}r_{0,k}^{-\alpha}-\beta\sum\nolimits_{x_{i,j}\in\cup_{j=1}^{K}\Phi_{j}^{\backslash
x_{0,k}}}P_{j}r_{i,j}^{-\alpha}G^{L}_{s,j}G^{U}_{s}\geq\beta\sigma^{2}\right\\}$
$\displaystyle=$
$\displaystyle\int_{H_{k}}^{R_{\mathrm{max}_{k}}}\mathbb{P}\left\\{\chi(r_{0,k})\geq\beta\sigma^{2}\right\\}f_{D,k}(r_{0,k}|\mathcal{A}_{D,k})\mathrm{d}r_{0,k},$

where $\chi(r_{0,k})=\mathcal{S}(r_{0,k})-\beta\mathcal{I}_{D}(r_{0,k})$,
$\mathcal{S}(r_{0,k})=P_{k}G_{m,k}^{L}G^{U}_{m}r_{0,k}^{-\alpha}$, and

$\mathcal{I}_{D}(r_{0,k})=\sum\nolimits_{x_{i,j}\in\cup_{j=1}^{K}\Phi_{j}^{\backslash
x_{0,k}}}P_{j}r_{i,j}^{-\alpha}G^{L}_{s,j}G^{U}_{s}.$ (51)

Then, by applying the Gil-Pelaez inversion theorem [37, 38, 39], the
probability term in (50) can be further computed as

$\displaystyle\mathbb{P}\left\\{\chi(r_{0,k})\geq\beta\sigma^{2}\right\\}=\frac{1}{2}-\frac{1}{\pi}\int_{0}^{\infty}\frac{\Im\left\\{\psi_{\chi(r_{0,k})}(t)\exp(it\beta\sigma^{2})\right\\}}{t}\mathrm{d}t,$
(52)

where $\Im\\{\cdot\\}$ is the imaginary operator and $\psi_{\chi(r_{0,k})}$ is
the characteristic function of $\chi(r_{0,k})$ which is given by

$\displaystyle\psi_{\chi(r_{0,k})}(t)=\mathbb{E}\left\\{\exp\bigg{(}-it\big{(}\mathcal{S}(r_{0,k})-\beta\mathcal{I}_{D}(r_{0,k})\big{)}\bigg{)}\right\\}=e^{-it\mathcal{S}(r_{0,k})}\mathcal{L}^{\ddagger}_{\mathcal{I}_{D}}(-it\beta),$
(53)

where
$\mathcal{L}^{\ddagger}_{\mathcal{I}_{D}}(z)=\mathbb{E}\left\\{e^{-z\mathcal{I}_{D}(r_{0,k})}\right\\}$
is the Laplace transform of the interference, which can be evaluated as
following

$\displaystyle\mathcal{L}^{\ddagger}_{\mathcal{I}_{D}}(z)$
$\displaystyle=\mathbb{E}\left\\{\prod\nolimits_{x_{i,j}\in\cup_{j=1}^{K}\Phi_{j}^{\backslash
x_{0,k}}}\exp\left(-{zP_{j}r_{i,j}^{-\alpha}G^{L}_{s,j}G^{U}_{s}}\right)\right\\}$
(54) $\displaystyle=$
$\displaystyle\mathbb{E}\left\\{\prod_{j=1}^{K}\exp\left(-2\pi\lambda_{j}\frac{R_{\oplus}+H_{j}}{R_{\oplus}}\int_{\max\\{r_{0,k},H_{j}\\}}^{R_{\mathrm{max}_{j}}}\big{(}1-\exp\left(-{zP_{j}r^{-\alpha}G^{L}_{s,j}G^{U}_{s}}\right)\big{)}r\mathrm{d}r\right)\right\\}.$

Hence, by solving the above integral and by substituting (52)-(54) into
$\eqref{AppF1}$, the final expression of $\mathcal{P}_{D,k}^{\ddagger}(\beta)$
is obtained. Finally, by the law of total probability, the coverage
probability of the DbA scheme for the non-fading scenario is derived.

## References

  * [1] Y. Guo, C. Skouroumounis, and I. Krikidis, “Large-scale heterogeneous ultra-dense leo satellite-based cellular networks,” in _Proc. IEEE Int. Conf. Commun. (ICC)_ , Rome, Italy, May. 2023.
  * [2] D. C. Nguyen, M. Ding, P. N. Pathirana, A. Seneviratne, J. Li, D. Niyato, O. Dobre, and H. V. Poor, “6G internet of things: A comprehensive survey,” _IEEE Internet Things J._ , vol. 9, no. 1, pp. 359–383, Jan. 2022\.
  * [3] X. Fang, W. Feng, T. Wei, Y. Chen, N. Ge, and C.-X. Wang, “5G embraces satellites for 6G ubiquitous IoT: Basic models for integrated satellite terrestrial networks,” _IEEE Internet Things J._ , vol. 8, no. 18, pp. 14 399–14 417, Sep. 2021.
  * [4] E. Yaacoub and M.-S. Alouini, “A key 6G challenge and opportunity—connecting the base of the pyramid: A survey on rural connectivity,” _Proc. IEEE_ , vol. 108, no. 4, pp. 533–582, Apr. 2020.
  * [5] X. Zhu and C. Jiang, “Integrated satellite-terrestrial networks toward 6G: Architectures, applications, and challenges,” _IEEE Internet Things J._ , vol. 9, no. 1, pp. 437–461, Jan. 2022.
  * [6] M. Centenaro, C. E. Costa, F. Granelli, C. Sacchi, and L. Vangelista, “A survey on technologies, standards and open challenges in satellite IoT,” _IEEE Commun. Surv. Tutor._ , vol. 23, no. 3, pp. 1693–1720, May 2021.
  * [7] 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, Aug. 2017\.
  * [8] B. Di, L. Song, Y. Li, and H. V. Poor, “Ultra-dense LEO: Integration of satellite access networks into 5G and beyond,” _IEEE Wirel. Commun._ , vol. 26, no. 2, pp. 62–69, Apr. 2019.
  * [9] A. Mohammed, A. Mehmood, F.-N. Pavlidou, and M. Mohorcic, “The role of high-altitude platforms (HAPs) in the global wireless connectivity,” _Proc. IEEE_ , vol. 99, no. 11, pp. 1939–1953, Nov. 2011.
  * [10] S. Liu, Z. Gao, Y. Wu, D. W. Kwan Ng, X. Gao, K.-K. Wong, S. Chatzinotas, and B. Ottersten, “LEO satellite constellations for 5G and beyond: How will they reshape vertical domains?” _IEEE Commun. Mag._ , vol. 59, no. 7, pp. 30–36, Jul. 2021.
  * [11] H. Xie, Y. Zhan, G. Zeng, and X. Pan, “LEO mega-constellations for 6G global coverage: Challenges and opportunities,” _IEEE Access_ , vol. 9, pp. 164 223–164 244, Dec. 2021.
  * [12] H. Al-Hraishawi, H. Chougrani, S. Kisseleff, E. Lagunas, and S. Chatzinotas, “A survey on nongeostationary satellite systems: The communication perspective,” _IEEE Commun. Surv. Tutor._ , vol. 25, no. 1, pp. 101–132, Aug. 2022.
  * [13] R. Wang, M. A. Kishk, and M.-S. Alouini, “Ultra-dense LEO satellite-based communication systems: A novel modeling technique,” _IEEE Commun. Mag._ , vol. 60, no. 4, pp. 25–31, Apr. 2022.
  * [14] A. Al-Hourani, “An analytic approach for modeling the coverage performance of dense satellite networks,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 4, pp. 897–901, Apr. 2021.
  * [15] ——, “Optimal satellite constellation altitude for maximal coverage,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 7, pp. 1444–1448, Jul. 2021\.
  * [16] D.-H. Jung, J.-G. Ryu, W.-J. Byun, and J. Choi, “Performance analysis of satellite communication system under the shadowed-rician fading: A stochastic geometry approach,” _IEEE Trans. Commun._ , vol. 70, no. 4, pp. 2707–2721, Apr. 2022.
  * [17] H. Jia, C. Jiang, L. Kuang, and J. Lu, “An analytic approach for modeling uplink performance of mega constellations,” _IEEE Trans. Veh. Technol._ , vol. 72, no. 2, pp. 2258–2268, Feb. 2023.
  * [18] N. Okati and T. Riihonen, “Nonhomogeneous stochastic geometry analysis of massive LEO communication constellations,” _IEEE Trans. Commun._ , vol. 70, no. 3, pp. 1848–1860, Mar. 2022.
  * [19] A. Talgat, M. A. Kishk, and M.-S. Alouini, “Stochastic geometry-based analysis of LEO satellite communication systems,” _IEEE Commun. Lett._ , vol. 25, no. 8, pp. 2458–2462, Aug. 2021.
  * [20] Z. Qiu, W. Wang, J. Geng, and Y. Liu, “Interference analysis of multi-tier NGSO based on stochastic geometry,” in _2023 IEEE Wireless Communications and Networking Conference (WCNC)_ , 2023, pp. 1–6.
  * [21] C.-S. Choi and F. Baccelli, “Cox point processes for multi-altitude leo satellite networks,” _arXiv preprint arXiv:2301.02469_ , 2023.
  * [22] R. Wang, M. A. Kishk, and M.-S. Alouini, “Reliability analysis of multi-hop routing in multi-tier leo satellite networks,” _IEEE Transactions on Wireless Communications_ , pp. 1–1, 2023.
  * [23] M. Di Renzo and W. Lu, “System-level analysis and optimization of cellular networks with simultaneous wireless information and power transfer: Stochastic geometry modeling,” _IEEE Trans. Veh. Technol._ , vol. 66, no. 3, pp. 2251–2275, Mar. 2017.
  * [24] H.-S. Jo, Y. J. Sang, P. Xia, and J. G. Andrews, “Heterogeneous cellular networks with flexible cell association: A comprehensive downlink SINR analysis,” _IEEE Trans. Wireless Commun._ , vol. 11, no. 10, pp. 3484–3495, Oct. 2012.
  * [25] M. Alzenad and H. Yanikomeroglu, “Coverage and rate analysis for vertical heterogeneous networks (VHetNets),” _IEEE Trans. Wireless Commun._ , vol. 18, no. 12, pp. 5643–5657, Sep. 2019.
  * [26] M. Haenggi, _Stochastic Geometry for Wireless Networks._ Cambridge University Press, 2012.
  * [27] 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 Trans. Commun._ , vol. 68, no. 8, pp. 5120–5134, Aug. 2020.
  * [28] A. Abdi, W. Lau, M.-S. Alouini, and M. Kaveh, “A new simple model for land mobile satellite channels: first- and second-order statistics,” _IEEE Trans. Wireless Commun._ , vol. 2, no. 3, pp. 519–528, May 2003.
  * [29] J. Park, J. Choi, and N. Lee, “A tractable approach to coverage analysis in downlink satellite networks,” _IEEE Trans. Wireless Commun._ , vol. 22, no. 2, pp. 793–807, Feb. 2023.
  * [30] ITU-R, “Simulation methodologies for determining statistics of short-term interference between co-frequency, codirectional non-geostationary-satellite orbit fixed-satellite service systems in circular orbits and other non-geostationary fixed-satellite service systems in circular orbits or geostationary-satellite orbit fixed-satellite service networks,” 2003.
  * [31] A. Al-Hourani, “Session duration between handovers in dense LEO satellite networks,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 12, pp. 2810–2814, Dec. 2021.
  * [32] Y. Yang, X. Yue, and A. G. Dempster, “GPS-based onboard real-time orbit determination for leo satellites using consider Kalman filter,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 52, no. 2, pp. 769–777, Apr. 2016.
  * [33] A. Al-Hourani and I. Guvenc, “On modeling satellite-to-ground path-loss in urban environments,” _IEEE Commun. Lett._ , vol. 25, no. 3, pp. 696–700, Mar. 2021.
  * [34] S. S. Kalamkar, F. Baccelli, F. M. Abinader, A. S. M. Fani, and L. G. U. Garcia, “Beam management in 5G: A stochastic geometry analysis,” _IEEE Trans. Wireless Commun._ , vol. 21, no. 4, pp. 2275–2290, Apr. 2022\.
  * [35] H.-C. Yang and M.-S. Alouini, _Order statistics in wireless communications: diversity, adaptation, and scheduling in MIMO and OFDM systems_. Cambridge University Press, 2011\.
  * [36] I. S. Gradshteyn and I. M. Ryzhik, _Table of Integrals, Series, and Products_. New York, NY, USA: Elsevier, 2007.
  * [37] J. Gil-Pelaez, “Note on the inversion theorem,” _Biometrika_ , vol. 38, no. 3-4, pp. 481–482, Dec. 1951.
  * [38] J. Sayehvand and H. Tabassum, “Interference and coverage analysis in coexisting RF and dense TeraHertz wireless networks,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 10, pp. 1738–1742, Oct. 2020.
  * [39] Y. Guo, C. Skouroumounis, and I. Krikidis, “Joint information and energy transfer of SWIPT-enabled mobile users in wireless networks,” _IEEE Trans. Green Commun. Netw._ , vol. 6, no. 2, pp. 1141–1156, Jun. 2022.