# A Novel Analytical Model for LEO and MEO Satellite Networks based on Cox Point Processes

Chang-Sik Choi, François Baccelli

<figure><img src="image_1.png" /></figure>

## To cite this version:

Chang-Sik Choi, François Baccelli. A Novel Analytical Model for LEO and MEO Satellite Networks based on Cox Point Processes. IEEE Transactions on Communications, In press, pp.1-11. 10.1109/TCOMM.2024.3471978. [hal-04884682](https://hal.science/hal-04884682v1)

HAL Id: hal-04884682

https://hal.science/hal-04884682v1

Submitted on 3 Jun 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Novel Analytical Model for LEO and MEO
Satellite Networks based on Cox Point Processes

Chang-Sik Choi and François Baccelli

**Abstract**—This work develops an analytical framework for downlink low Earth orbit (LEO) or medium Earth orbit (MEO) satellite communications, leveraging tools from stochastic geometry. We propose a tractable approach to the analysis of such satellite communication systems, accounting for the fact that satellites are located on circular orbits. We accurately incorporate this geometric property of LEO or MEO satellite constellations by developing a Cox point process model that jointly produces orbits and satellites on these orbits. Our work contrasts with previous modeling studies that presumed satellite locations to be entirely random, thereby overlooking the fundamental fact that satellites are jointly positioned on orbits. Employing this Cox model, we analyze the network performance experienced by users located on Earth. Specifically, we evaluate the no-satellite probability of the proposed network and the Laplace transform of the interference created by such a network. Using it, we compute its SIR (signal-to-interference) distribution, namely its coverage probability. By presenting fundamental network performance as functions of key parameters, this model allows one to assess the statistical properties of downlink LEO or MEO satellite communications and can thus be used as a system-level design tool to operate and optimize forthcoming complex LEO or MEO satellite networks.

**Index Terms**—LEO satellite networks, MEO satellite networks, stochastic geometry, coverage probability, Cox point process, isotropic model.

## I. INTRODUCTION

### A. Motivation and Related Work

Satellite communications provide global-scale connectivity to users everywhere on earth without the need for deploying base stations and infrastructure on the ground [1], [2]. As orbiting the Earth at very high speeds, LEO and MEO satellites provide reliable and fast Internet connectivity to millions of devices [1], [2]. In the early stage of satellite communications, the number of deployed satellites was very small, and thus only a limited number of satellites were available for connections. For instance, the Iridium constellation [3]—6 orbits with 11 satellites on each orbit—provided a call coverage of up to 7 minutes. More precisely, because of the motions of satellites and their limited numbers, the call was either dropped or transferred after about 7 minutes. Preventing outages was one of the key design criteria for satellite communications of this early stage.

In modern LEO or MEO satellite communication systems with numerous satellites, the goal is not only to prevent

call outages but also to provide high-speed and low-latency Internet connections to millions of devices. Comprised of a large number of satellites [4], [5], often integrated with terrestrial network infrastructures [6], [7], satellite communication systems are envisioned to handle various demands of ground or even aerial devices around the globe. Recently, several companies intended to establish their satellite communication networks by building such large constellations for global connectivity [8]–[10]. It is not hard to imagine a large number of satellites from various companies will provide Internet connections to devices anywhere on earth. Since the spatial distribution of satellites determines the performance of communications on satellite networks, building an analytical framework is vital to the description and analysis of the satellite communications.

To provide an analytical framework for describing the locations of LEO or MEO satellites, several recent papers [11]–[16] used binomial or Poisson point process models in stochastic geometry [17]–[19]. The main advantage of employing these models lies in the fact that it not only provides geometric treatment to the interference and coverage in such networks but also identifies key network performance behaviors as functions of networks’ distributional parameters [20]–[23]. To take advantage of assessing typical network performance and of getting high-level insights, [11]–[16] modeled the distribution of LEO or MEO satellites as binomial or Poisson point processes, where the satellites’ locations are assumed to be uniformly distributed on a sphere. In [11]–[15], the impact of interference created by a binomial satellite point process was evaluated and accounted for the SINR coverage probability in [16]. It is important to note that the recent modeling technique based on binomial or Poisson point processes [11]–[16] overlooked the essential geometric facts that LEO or MEO satellite constellations are comprised of orbits and that they are always on those orbits. Although binomial or Poisson models provide new results on the performance metrics of downlink communications, the existence of an orbital structure and the impact of the structure to the network performance cannot be examined. This motivates us to develop an isotropic Cox point process model that incorporates these geometric characteristics. This process jointly produces orbits and satellites on these orbits, isotropically distributed in space. This Cox model allows one to analyze the network performance seen by users anywhere on Earth. To the best of our knowledge, this is the first work that explicitly models the orbits and satellites on them and analyzes the typical network performance of downlink communications from satellites to network users anywhere on

Chang-Sik Choi is with Hongik University, South Korea. (email: chang-sik.choi@hongik.ac.kr). François Baccelli is with Inria Paris and Telecom Paris, France. (email: francois.baccelli@inria.fr). The work of Chang-Sik Choi was supported by the NRF RS-2024-00334240. The work of Francois Baccelli was supported by the ERC NEMO grant 788851 to INRIA and by the French National Agency for Research via the project n°ANR-22-PEFT-0010 of the France 2030 program PEPR réseaux du Futur.
Earth. To demonstrate the use of the developed analytical framework, we derive the no-satellite probability and then the coverage probability of downlink communication from satellite to a typical user. It is important to note that analytical models such as the binomial and Poisson models approximate an existing or forthcoming LEO or MEO constellation by delineating the local geometry of the target constellation. Such a local approximation of the constellation is known to produce coverage probability approximating well that of the target constellation [11]–[14]. Similar to those studies, we also show that the Cox model effectively represents an existing or forthcoming constellation by better approximating the target network’s local geometry as seen by the typical user. We present numerical experiments demonstrating the use of the proposed framework.

## B. Theoretical Contributions

1) **Orbit-based Stochastic Geometry Model For LEO or MEO satellite constellations:** This work develops an analytical framework for LEO or MEO satellite networks using stochastic geometry. The proposed Cox model for LEO or MEO satellites stands in contrast to most existing frameworks that have employed binomial point processes to describe the locations of satellites as random points. While recent works show that the binomial or Poisson point process models locally portray some existing satellite constellations [11]–[16], this approach overlooked the fact that satellites are always located on orbits, which affects network geometry and the networks therein. This Cox framework serves as one among several modeling techniques for LEO or MEO satellite networks. To demonstrate its applicability, Section VI shows that how the Cox model approximates an upcoming constellation in terms of the coverage probability. It also demonstrates that how the Cox point process model can reproduce LEO or MEO satellite constellations by adjusting its two parameters: the mean number of orbits and the mean number of satellites per orbit.

2) **Performance Analysis of Downlink Communications::**

We first prove that our developed Cox point process is isotropic, namely invariant by all rotations. Leveraging this rotation invariance property, we derive a closed-form expression for the arc length of an orbit and provide an integral formula for the distribution of the distance from a typical user to its nearest satellite. These expressions are essential in the derivation of the network performance of LEO or MEO satellite downlink communications. For example, leveraging this, we derive the probability that a typical user has no visible satellite as a function of $\lambda$ the orbit density, $\mu$ the density of satellites per orbit, and $r_a$ the satellite altitude. Further, we focus on the interplay of network geometry and derive the SIR coverage probability of the typical user under Nakagami-m fading by assuming that the interference power dominates the noise. We compare the derived coverage formula to the results obtained by Monte Carlo simulations and validate the derived formulas.

3) **Design Insights for Practical LEO and MEO satellite Networks:** Our approach provides a comprehensive tool for

<figure><img src="image_3.png" /><figcaption>Fig. 1. Let A be the ascending/descending point. The longitude $\theta$ is the angle that $\overline{OA}$ makes with the x-axis. The inclination $\phi$ is the angle that the orbital plane makes with the reference plane. $\omega$ is the angle that the $\overline{OX}$ makes with $\overline{OA}$.</figcaption></figure>

designing and enhancing LEO or MEO satellite communication systems. For example, we derive an expression for the interference experienced by the typical user, which can be utilized to design interference management techniques in densely-deployed LEO or MEO satellite communication systems. Additionally, this paper obtains expressions for the no-satellite probability and the coverage probability as the functions of key parameters such as $\lambda$ and $\mu$, allowing us to assess each variable’s individual impact on the system’s large-scale performance. Leveraging the derived formula, network operators of satellite systems can explore and assess various deployment options by controlling the number of satellites per orbit or the number of orbits per altitude, without time-consuming system-level simulation. Lastly, the proposed model and analysis can serve as a basis for assessing the time-domain performance metrics of satellite networks, such as the time fraction of coverage or delay.

# II. SYSTEM MODEL

## A. Models for Orbits and Satellites

We denote by $r_e$ the radius of Earth ($r_e \approx 6400$ km) and the center of Earth is located at the origin of the Euclidean space. The reference plane is the $xy$ plane and the longitudinal zero point is the $x$-axis. We denote by $r_a$ the altitude of the satellites or equivalently the radius of orbits. All orbits are assumed to be circles centered at the origin. Let $r_s = r_e + r_a$, where $r_s$ is the radius of orbits.

To model circular orbits of LEO or MEO satellites in the simplest case, we first consider a Poisson point process $\Xi$ of density $\lambda \sin(\phi)/(2\pi)$ on the rectangle set $\mathcal{R}=[0, \pi) \times [0, \pi)$. We write $\Xi = \sum_i \delta_{(\theta_i, \phi_i)}$. Each point of the Poisson point process on $\mathcal{R}$, say $(\theta, \phi)$, is mapped to an undirected orbit $l(\theta, \phi)$ in the Euclidean space $\mathbb{R}^3$. Specifically, $\theta$ gives the longitude of the orbit $l(\theta, \phi)$ and $\phi$ is the inclination of the orbit $l(\theta, \phi)$. See Fig. 1 for the longitude and inclination. The orbit process $\mathcal{O}$ on $\mathbb{R}^3$ is

$$ \mathcal{O} = \bigcup_{i \in \mathbb{Z}} l(\theta_i, \phi_i) = \bigcup_{(\theta_i, \phi_i) \in \Xi} l(\theta_i, \phi_i). \quad (1) $$
<figure><img src="image_4.png" /><figcaption>Fig. 2. The proposed model with $\lambda = 30$, $\mu = 40$, and $r_s = 7000$ km.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. The proposed model with $\lambda = 28$, $\mu = 120$, and $r_s = 7000$ km.</figcaption></figure>

The orbit process is a union of circles located on the sphere $S_{r_s} = \{(x, y, z) | x^2 + y^2 + z^2 = r_s^2\}$. Since the density of the Poisson point process $\Xi$ is $\lambda \sin(\phi)/(2\pi)$, there are $\lambda$ points on $\mathcal{R}$ on average, or equivalently, there are $\lambda$ orbits on $\mathbb{R}^3$ on average. The orbit process $\mathcal{O}$ is isotropic, namely rotation invariant. We will shortly prove its rotation invariance property in Section III.

To represent the locations of satellites on orbits, we leverage the conditional structure based on a Cox point process. Specifically, conditionally on each orbit $l(\theta_i, \phi_i)$, the satellites on each orbit are modeled as a Poisson point process $\psi_i$ of mean $\mu$. Because of the conditional structure, our developed satellite model geometrically ensures that the satellites are exclusively on orbits given by $\mathcal{O}$. Based on the above construction, the orbital angles $\omega$ of satellites (See. Fig. 1.) are also characterized as a Poisson point process of intensity $\mu/(2\pi)$ on the finite interval $[0, 2\pi)$.

Collectively, the satellite point process $\Psi$ on $\mathbb{R}^3$ is $\Psi = \sum_{i \in \mathbb{Z}} \psi_i$. The satellite point process is defined conditionally on the orbit process and it is hence a Cox point process [17], [18].

It is worth noting that the conditional structure of Cox point process was found to be useful in the modeling of vehicular networks on two-dimensional $\mathbb{R}^2$ plane since it jointly creates road systems and vehicles on them [22], [24], [25]. In the same way as vehicles are on roads in two-dimensional vehicular networks, satellites are on orbits in three-dimensional satellite networks.

Figs. 2–3 illustrate the proposed network model with various $\lambda$ and $\mu$ parameters. The proposed stochastic geometry framework of populating orbits and their corresponding satellites is designed to easily change the number of satellites, the number of orbits, or their topological characteristics. For instance, we may increase the total number of satellites by individually increasing $\lambda$ or $\mu$ or both. Fig. 3 shows $\mu = 120$ and since satellites are densely populated on the finite orbital planes, the clustering of the satellites on their orbits are more pronounced.

## B. User Location

In [11]–[16], users are assumed to be uniformly distributed on Earth. In the aforementioned work, the network users are specifically modeled as a Poisson point process of intensity $\lambda_u/(4\pi^2 r_e^2)$ on the sphere $\{(x, y, z) \in \mathbb{R}^3 | x^2 + y^2 + z^2 = r_e^2\}$ where $\lambda_u$ is the number of network users.

To reduce the computational complexity accompanied with the distribution of users, we also assume that users are uniformly distributed on Earth, independent of satellites. For downlink communications, users are assumed to get their downlink signals from their closest satellites [11]–[16], [26].

This nearest association assumption facilitates the analysis of SINR coverage by taking the most significant signal as the desired signal and treating the signals from other satellites as interference. In practice the nearest association can be implemented to each UE, first by detecting various signals from visible satellites and then by determining the closest satellites by time-averaging of those signals.

## C. Propagation Model

The downlink signals from satellites attenuate because of Doppler shift, weather, rain, reflection from objects, and path loss over the space. In practice, Doppler shift can be compensated by exploiting the existing data on satellites' orbits and speeds [27], [28]. A few large-scale uniform factors such as rain attenuation can be ignored under certain conditions [13]–[15], [28]. To emphasize the role of satellite geometry onto the network performance, we propose a simple propagation model where various attenuating and propagation factors are assumed to be aggregated to a single independent fading random variable. This approach was found useful in recent work [13]–[15], where analytical tractability was achieved by focusing on the impact of topology of satellite communication systems.

Suppose a transmitter and a user separated by a distance $d$. Assume the transmitter is visible to the user. We assume that the received signal power at the user is then given by

$$pG_1(\zeta)G_2(\zeta)Hd^{-\alpha}, \quad (2)$$
where $p$ is the receive signal power at 1 meter, $G_1(\zeta)$ is the transmit satellite antenna gain, $G_2(\zeta)$ is the receive antenna gain of network users, $H$ is a random variable representing small-scale fading, and $\alpha$ is the path loss exponent.

To make the analysis tractable, we assume that network users and satellites are able to direct their antennas toward their associated counterparts by using technologies such as phased antenna array [15], [29]. We assume that the gain $G_1(\zeta)$ of the LEO or MEO satellites is

$$G_1(\zeta) = \begin{cases} g & |\zeta| < \zeta_t, \\ 1 & |\zeta| > \zeta_t, \end{cases} \quad (3)$$

where $\zeta$ is the boresight angle from the antenna's maximum radiated power direction.

Similarly, network users are assumed to have isotropic antennas $G_2(\zeta) = g_r = 1$ for all $\zeta$, namely 0 dBi [26], [30].

Assuming general small-fading, we use a Nakagami-m fading [16]. The CCDF of the random variable $H$ is

$$\mathbf{P}(H > x) = e^{-mx} \sum_{k=0}^{m-1} \frac{(mx)^k}{k!} \quad \forall x \ge 0,$$

where $m \ge 1$. When $m = 1$, Nakagami-m fading becomes Rayleigh fading and $H$ becomes an exponential random variable with mean one.

### D. Performance Metrics

A downlink satellite communication system is built to cover all network users. Nevertheless, some users on the surface of the earth may not be well covered because of the lack of visible satellite or weak signals.

We define the no-satellite event as the event that there is no visible satellite from a network user. Since the Earth rotates and the relative locations of satellites vary over time, we characterize this through the no-satellite probability.

Then, to evaluate the basic performance of downlink satellite communication systems, we study the coverage probability, namely the CCDF of the SINR of network users. The coverage probability incorporates the path loss, the topological properties of the association satellite and of the interfering satellites, the small-scale fading, the antenna gains, and the background thermal noise. The coverage probability of a user at an arbitrary location $y$, or equivalently the CCDF of SINR of the user at $y$, is given by

$$\begin{aligned} \mathbf{P}(\text{SINR}_y > \tau) &= \mathbf{P}_y \left( \frac{pG_{X_y;y}H_i \|X_y - y\|^\alpha}{kTB_w + \sum_{X_j \in \bar{\Psi}_\backslash X_y} pG_{X_j;y}H \|X_j - y\|^\alpha} > \tau \right) \\ &= \mathbf{P}_y \left( \frac{pgH_i \|X_y - y\|^\alpha}{kTB_w + \sum_{X_j \in \bar{\Psi}_\backslash X_y} pH_j \|X_j - y\|^\alpha} > \tau \right), \quad (4) \end{aligned}$$

where $X_y$ is the location of the satellite that serves the user at $y$, $G_{X_i;y}$ is the aggregate antenna gain from the satellite at $X_i$ toward the user at $y$, $kTB_w$ is the noise power, and $\bar{\Psi}_y$ is the point process of satellites visible at $y$. The transmit or association satellite is assumed to direct its beam toward its intended receiver. Namely we let $G_{X_y;y} = gg_r = g$ and

$G_{X_y;y} = g_r = 1$ for all $X_j \neq X_y$. We use subscripts to distinguish satellites' locations and small-scale fading associated with them. The constant $\tau$ is the SINR threshold.

### III. STATISTICS OF ORBITS AND SATELLITES

This section provides statistical properties of the proposed satellite Cox point process that are essential to the analysis of the downlink communications from satellites to users.

#### A. Isotropy

**Theorem 1.** *O and Ψ are isotropic, namely invariant by all rotations.*

*Proof:* Below, the reference basis of $\mathbb{R}^3$ is denoted by $(e_x, e_y, e_z)$ and the unit sphere of center $o = (0, 0, 0)$ in $\mathbb{R}^3$ is denoted by $\mathbb{S}$. Let $U$ be a uniformly distributed random point on $\mathbb{S}$. For each $U$, there is a unique directed orbit $\mathcal{O}(U) \subset \mathbb{S}$, which is the orbit whose normal vector is $\vec{oU}$ with a direction that is the trigonometric (counterclockwise) direction with respect to (w.r.t.) this vector (namely seen from point $U$). Note that $\mathcal{O}(U)$ is a factor of $U$; namely, for all $\mathbb{R}^3$ rotations $R$ of center $o$, the orbit $\mathcal{O}(R(U))$ coincides with the orbit $R(\mathcal{O}(U))$. Since the law of $U$ is isotropic on $\mathbb{S}$ (namely left invariant by all $R$), it follows from the relation $\mathcal{O}(R(U)) = R(\mathcal{O}(U))$ that the law of $\mathcal{O}(U)$ is also isotropic.

It is well known [31] that the uniform random vector $U$ can be represented as

$$U = (\sqrt{1-V^2} \cos(\Theta), \sqrt{1-V^2} \sin(\Theta), V), \quad (5)$$

in the basis $(e_x, e_y, e_z)$ where $V \sim \text{Uniform}(-1, 1)$, $\Theta \sim \text{Uniform}(0, 2\pi)$, and $V \perp \Theta$.

For the directed orbit $\mathcal{O}(U)$, we define

* The longitude angle $\theta \in [0, 2\pi)$ to be the angle $\angle e_x \vec{oU}$, where $A$ is the ascending point of the orbit on the xy-plane;

* The inclination angle $\phi \in [0, \pi)$ to be the angle $\angle e_z \vec{oU}$.

It follows from Eq (5) that we can represent the longitude and inclination as follows:

$$\theta = \Theta + \pi/2 \quad \mod 2\pi, \quad (6)$$

$$\phi = \arccos(V). \quad (7)$$

Since $V$ and $\Theta$ are independent random variables, we have that $\phi$ and $\theta$ are independent. In other words, the longitude and inclination of the orbit are independent. Furthermore, since $\Theta \sim \text{Uniform}[0, 2\pi)$, we also have $\theta \sim \text{Uniform}[0, 2\pi)$. Then, based on the fact that $V \sim \text{Uniform}[-1, 1]$, we arrive at

$$\mathbb{P}(\phi < x) = \mathbb{P}(V < \cos(x)) = \frac{1 - \cos(x)}{2}, \quad 0 \le x < \pi.$$

Using the above CDF, we get the PDF of $\phi$ as follows:

$$f_{\phi}(x) = \begin{cases} \frac{\sin(x)}{2}, & \text{for } 0 \le x < \pi, \\ 0 & \text{otherwise.} \end{cases} \quad (8)$$

The isotropic directed orbit Poisson point process can hence be represented as a Poisson point process of density

$$\Lambda(\phi, \theta) = \frac{\lambda}{4\pi} \sin(\phi), \quad (9)$$
on the rectangle set $\mathcal{R} = [0, \pi) \times [0, 2\pi)$. Here, $\lambda$ corresponds to the mean number of directed orbits.

Furthermore, the directed orbit with angles $(\theta, \phi)$ and that with angles $(\theta + \pi, \phi + \pi/2 \mod \pi)$ reduce to the same orbit when forgetting the orbit direction. Therefore, for an undirected isotropic orbit, its longitude angle $\tilde{\theta}$ is defined as the angle that the orbital plane makes with the reference plane in $[0, \pi)$ and it is uniformly distributed in this interval.

Based on the same principle, the isotropic undirected orbit Poisson point process can hence be represented as a Poisson point process of density

$$ \tilde{\Lambda}(\tilde{\phi}, \tilde{\theta}) = \frac{\tilde{\lambda}}{2\pi} \sin(\tilde{\phi}), \quad (10) $$

on the rectangle set $\tilde{\mathcal{R}} = [0, \pi] \times [0, \pi]$. Here, $\tilde{\lambda}$ is the mean number of undirected orbits.

Since the density of the orbit process is given by Eq. (10), the proposed orbit process is isotropic. Furthermore, conditionally on each orbit, the satellite point process on each orbit is isotropic. As a result, the Satellite Cox point process is isotropic. ■

**Remark 1.** Note that this isotropy property also appears in recent work [11]–[16] where binomial or Poisson point processes were used to model the distribution of LEO or MEO satellites. Those isotropic models are capable of effectively depicting the local geometry of non isotropic satellite deployment scenarios, by using the mean parameters obtained locally from such scenarios at a given latitude. Similarly, our developed framework is also capable of reproducing the local geometry of satellite distribution and analyzing the downlink communications therein, by finding the mean number of orbits and the mean number of satellites per orbit which are based on a real-world scenario at a given latitude and then applying them in the analysis. See Section VI.

**Lemma 1.** The average number of all satellites is $\lambda\mu$.

*Proof:* The average number of satellites is given by

$$ \begin{align} \mathbf{E}[\Psi] &\stackrel{(a)}{=} \mathbf{E} \left[ \sum_{l(\theta_i, \phi_i) \in \mathcal{O}} \sum_{X_j \in \psi_i} 1 \right] \nonumber \\ &\stackrel{(b)}{=} \mathbf{E} \left[ \sum_{l(\theta_i, \phi_i) \in \mathcal{O}} \mathbf{E} \left[ \sum_{X_j \in \psi_i} 1 \middle| \mathcal{O} \right] \right] \nonumber \\ &\stackrel{(c)}{=} \mathbf{E} \left[ \sum_{l(\theta_i, \phi_i) \in \mathcal{O}} \frac{\mu}{2\pi} \int_0^{2\pi} d\omega \right] \nonumber \\ &\stackrel{(d)}{=} \mu \int_0^{\pi} \int_0^{\pi} \frac{\lambda \sin(\phi)}{2\pi} d\theta \, d\phi = \lambda\mu, \tag{11} \end{align} $$

where (a) follows from the definition of the number of points on all orbits. By conditioning on $\mathcal{O}$, we have (b). Since the satellite point process on the orbit $l(\theta_i, \phi_i)$ is created by the Poisson point process $\phi_i$ of intensity $\mu/(2\pi)$ on $\mathcal{I}$, we can use Campbell's mean value theorem [18] to get (c). We obtain (d) from Campbell's mean value theorem on the Poisson point process $\mathcal{I}$.

**Remark 2.** In this paper, we consider a typical user at $n = (0, 0, r_e)$. Since the satellite Cox point process is isotropic and the users are independent of the satellites, the typical user at the above location can represent all the users in the network. In other words, the probability law of a satellite network seen from the typical user is the same as that of the satellite network seen from any users at any locations on Earth. Therefore, the average number of satellites visible from the typical user is the mean number of satellites visible from any users on the Earth. Moreover the statistics of the SINR or the interference seen by the typical user represent the distributions of the SINR or interference of any user in the network.

Below, we evaluate the average number of satellites visible from the typical user. First, let us define a spherical cap as follows:

$$ C_d = \{(x, y, z) \in \mathbb{S}_{r_s} | \sqrt{x^2 + y^2 + (z - r_e)^2} \le d\}, \quad (12) $$

where $d$ is a positive constant. This spherical cap is the set of points on $\mathbb{S}_{r_s}$ the sphere of radius $r_s$ whose distance to the typical user is less than or equal to $d$. Fig. 15 shows the spherical caps $C_d$ and $C_{\bar{d}}$, where $r_s - r_e < d < \bar{d}$ and $\bar{d} = \sqrt{r_s^2 - r_e^2}$, the maximum distance to a visible satellite.

**Proposition 1.** From any user on the Earth, the average number of visible satellites is given by

$$ \frac{\lambda\mu}{\pi} \int_{0}^{\bar{\varphi}} \cos(\varphi) \arcsin \left( \sqrt{1 - \cos^2(\bar{\varphi}) \sec^2(\varphi)} \right) d\varphi, \quad (13) $$

where $\bar{\varphi} = \arccos(r_e/r_s)$.

*Proof:* Consider the spherical cap $C_{\bar{d}}$. Let $X_{j,i}$ be the location of the j-th satellite on the i-th orbit. Then, the average number of visible satellites from the typical user is given by

$$ \begin{align} \mathbf{E} \left[ \sum_{l(\theta_i, \phi_i) \in \mathcal{O}} \sum_{X_{j,i} \in \psi_i} \mathbf{1}_{X_{j,i} \in C_{\bar{d}}} \right] \nonumber \\ &= \mathbf{E} \left[ \sum_{l(\theta_i, \phi_i) \in \mathcal{O}} \mathbf{E} \left[ \sum_{X_{j,i} \in \psi_i} \mathbf{1}_{X_{j,i} \in C_{\bar{d}}} \middle| \mathcal{O} \right] \right] \nonumber \\ &\stackrel{(a)}{=} \mathbf{E} \left[ \sum_{l(\theta_i, \phi_i) \cap C_{\bar{d}} \neq \emptyset} \frac{\mu}{2\pi} \frac{\text{length}(l(\theta_i, \phi_i) \cap C_{\bar{d}})}{r_s} \right] \nonumber \\ &\stackrel{(b)}{=} \frac{\mu}{\pi} \mathbf{E} \left[ \sum_{l(\theta_i, \phi_i) \cap C_{\bar{d}} \neq \emptyset} \arcsin \left( \sqrt{1 - \cos^2(\bar{\varphi}) \csc^2(\phi_i)} \right) \right]. \end{align} $$

To obtain (a), we use the fact that the mean number of satellites on the arc $C_{\bar{d}} \cap l(\theta, \phi)$ is given by the product of the angle of the arc and of the intensity $\mu/(2\pi)$. To get (b), we use the result in Appendix A, which gives the arc length of the intersection of an orbit and a spherical cap.

Then, employing Campbell's averaging formula, the average number of visible satellites is given by

$$ \begin{align} \frac{\mu\lambda}{2\pi} \int_{\pi/2-\bar{\varphi}}^{\pi/2+\bar{\varphi}} \sin(\phi) \arcsin \left( \sqrt{1 - \cos^2(\bar{\varphi}) \csc^2(\phi)} \right) d\phi \nonumber \\ = \frac{\lambda\mu}{\pi} \int_{0}^{\bar{\varphi}} \cos(\varphi) \arcsin \left( \sqrt{1 - \cos^2(\bar{\varphi}) \sec^2(\varphi)} \right) d\varphi. \end{align} $$
To get the first expression, we use the fact that an orbit $l(\theta_i, \phi_i)$ meets $C_{\bar{d}}$ if and only if its inclination satisfies $|\pi/2 - \phi_i| \le \bar{\varphi} = \arccos(r_e/r_s)$. Then, we use the change of variables $\pi/2 - \phi = \varphi$ to obtain the final result. ■

In Eq. (13), the integration gives a constant. For instance, in a case of $r_a = 525$, the integration gives 0.038 and therefore, the mean number of visible satellites is given by $0.038\lambda\mu$. In a case of $r_a = 1100$ km, the integral formula gives 0.074 and there are about $0.074\lambda\mu$ visible satellites on average.

## IV. NO-SATELLITE PROBABILITY

### A. No-Satellite Probability

We leverage Theorem 1 to obtain performance metrics of network users by deriving those of the typical user. In this section, we derive (i) the no-satellite probability $\mathbf{P}(\text{no-satellite})$ and (ii) the distance distribution from the typical user to the association satellite, namely nearest visible satellite, $\mathbf{P}(D > d)$.

**Theorem 2.** The no-satellite probability of the typical user is

$$e^{-\lambda \int_{0}^{\bar{\varphi}} \cos(\varphi) \left(1 - e^{-\frac{\mu}{\pi} \arcsin\left(\sqrt{1 - \cos^2(\varphi) \sec^2(\varphi)/r_s^2}\right)}\right) d\varphi} \quad (14)$$

*Proof:* Let $D$ be the distance from the typical user to its nearest satellite. Then, the typical user observes no satellite if and only if the distance from the typical user to its nearest satellite $D$ is greater than $\bar{d}$. The no-satellite probability is given by

$$\begin{aligned}
\mathbf{P}(D > \bar{d}) \\
&= \mathbf{P}(\|X_i - n\| \ge \sqrt{r_s^2 - r_e^2}, \forall X_i \in \Psi) \\
&= \mathbf{P} \left( \bigcap_{(\theta_i, \phi_i) \in \Xi} \left( \bigcap_{X_{j,i} \in \psi_i} \|X_{j,i} - n\| \ge \bar{d} \right) \right) \\
&= \mathbf{E} \left[ \prod_{(\theta_i, \phi_i) \in \Xi} \mathbf{E} \left[ \prod_{X_{j,i} \in \psi_i} \mathbf{1}_{\|X_{j,i} - n\| > \bar{d}} \middle| \Xi \right] \right], \quad (15)
\end{aligned}$$

where we use the fact that, conditionally on orbits, the Poisson point processes on orbits $\{\psi_i\}_i$ are independent.

To evaluate Eq. (15), we consider the spherical cap $C_{\bar{d}}$. Then, we use the fact that the event that all satellites on orbit $l(\theta_i, \phi_i)$ are located at distances greater than $\bar{d}$ is equivalent to the event that the orbit $l(\theta_i, \phi_i)$ contains no satellite on $C_{\bar{d}}$. By using the density of the satellite Poisson point process $\psi_i$, we get

$$\mathbf{E} \left[ \prod_{X_{j,i} \in \psi_i} \mathbf{1}_{\|X_{j,i} - n\| > \bar{d}} \middle| \Xi \right] = e^{-\frac{\mu}{2\pi} \frac{1}{r_s} \nu(C_{\bar{d}} \cap l(\theta_i, \phi_i))},$$

where $\nu(C_{\bar{d}} \cap l(\theta_i, \phi_i))$ is the length of the arc $C_{\bar{d}} \cap l(\theta_i, \phi_i)$ given by

$$\nu(C_{\bar{d}} \cap l(\theta_i, \phi_i)) = 2r_s \arcsin(\sqrt{1 - \cos^2(\bar{\varphi}) \csc^2(\phi_i)}),$$

for $\pi/2 - \bar{\varphi} < \phi_i < \pi/2 + \bar{\varphi}$. Therefore, we have

$$\begin{aligned}
\mathbf{P}(D > \bar{d}) \\
&= \mathbf{E} \left[ \prod_{(\theta_i, \phi_i) \in \Xi} e^{-\mu \pi^{-1} \arcsin(\sqrt{1 - \cos^2(\bar{\varphi}) \csc^2(\varphi_i)})} \right] \\
&= e^{-\frac{\lambda}{2} \int_{\pi/2 - \bar{\varphi}}^{\pi/2 + \bar{\varphi}} \sin(\phi) \left(1 - e^{-\frac{\mu}{\pi} \arcsin(\sqrt{1 - \cos^2(\bar{\varphi}) \csc^2(\varphi)})}\right)} d\phi \\
&= e^{-\lambda \int_{0}^{\bar{\varphi}} \cos(\varphi) \left(1 - e^{-\frac{\mu}{\pi} \arcsin(\sqrt{1 - \cos^2(\bar{\varphi}) \sec^2(\varphi)})}\right)} d\varphi,
\end{aligned}$$

where we used the probability generating functional of the Poisson point process of density $\lambda \sin(\phi)/2\pi$ on the rectangle $\mathcal{R}$ and then the change of variables $\varphi = \pi/2 - \phi$. ■

The above no-satellite probability gives the probability that an arbitrarily located network user on the Earth is not able to find any satellite at a given time. The derived probability expression is a function of the network geometric parameters $\lambda, \mu$, and $r_s$.

Since $r_s, \lambda$ and $\mu$ are system parameters, large-scale impact of modifying geometry of the satellite network, such as increasing the number of orbits or decreasing the satellite density, can easily be evaluated. For instance, we observe that the no-satellite probability expression of Theorem 2 is exponential in the number of orbits. Figs. 4 and 5 show the no-satellite probability when the satellite altitude is 550 km. The simulation results are produced by Monte Carlo simulation with the sample size $N = 10^6$. For each simulation instance, we create orbits and their satellites based on the definition of the Cox point process. Then, we count the event that there is no visible satellite at the typical user and combine such events to get the probability. We confirm that the simulation results validate the accuracy of the derived formula in Theorem 2.

For the moment, suppose $\mu$ is very high. In this case, each orbit has a very large number of satellites and thus the no-satellite probability is strongly dictated by the distribution of orbits. From the above formula, with $\mu \gg 1$, we have $e^{-2r_s\mu} \approx 0$ and thus the above derived no-satellite probability is asymptotically approximated as $\exp(-\lambda \sin(\bar{\varphi}))$. For instance, if $r_s = 7000$ km, we have $\exp(-0.4\lambda)$. If $\lambda = 52$, the no-satellite probability is approximately $10^{-9}$.

### B. Distribution of the Distance to the Nearest Satellite

Let $D$ be the distance from the typical user to its nearest satellite. Note that when the nearest satellite is visible, we have $r_s - r_e \le D \le r_s^2 - r_e^2$. We set $D = \infty$ if there is no visible satellite from the typical user.

**Theorem 3.** $\mathbf{P}(D > d)$ the CCDF of the distance to the nearest satellite is given by $1$ for $d < r_s - r_e$. Then, for $r_s - r_e \le d < \sqrt{r_s^2 - r_e^2}$,

$$e^{-\lambda \int_{0}^{\xi(d)} \cos(\varphi) \left(1 - e^{-\frac{\mu}{\pi} \arcsin(\sqrt{1 - \cos^2(\xi(d)) \sec^2(\varphi)})}\right)} d\varphi.$$
where $\xi(d) = \arccos((r_s^2 + r_e^2 - d^2)/(2r_e r_s))$. Then, for $d \ge \sqrt{r_s^2 - r_e^2}$, it is equal to the no-satellite probability.

*Proof:* The distance to the nearest visible satellite is greater than $d$ if and only if all the satellites are at distances greater than $d$. Since the orbit radius is $r_s$, the CCDF is 1 for $d < r_s - r_e$.

For $r_s - r_e \le d < \sqrt{r_s^2 - r_e^2}$, we exploit Appendix A and the proof of Theorem 2 to obtain

$$ 
\begin{aligned}
\mathbf{P}(D > d) &= e^{-\lambda \int_{\pi/2 - \xi(d)}^{\pi/2 + \xi(d)} \frac{\sin(\phi)}{2} \left( 1 - e^{-\frac{\mu}{\pi} \sin^{-1}\left(\sqrt{1 - \cos^2(\xi(d)) \csc^2(\phi)}\right)} \right) d\phi \\
&= e^{-\lambda \int_{0}^{\xi} \cos(\varphi) \left( 1 - e^{-\frac{\mu}{\pi} \sin^{-1}\left(\sqrt{1 - \cos^2(\xi(d)) \sec^2(\varphi)}\right)} \right) d\varphi,
\end{aligned}
$$

where we use the fact that only the orbits having inclinations between $\pi/2 - \xi(d)$ and $\pi/2 + \xi(d)$ will meet the spherical cap $C_d$. For $r_s - r_e \le d < \sqrt{r_s^2 - r_e^2}$, we have

$$ \xi(d) = \arccos((r_s^2 + r_e^2 - d^2)/(2r_s r_e)). $$

Note that $d$ denotes the distance to the nearest satellite and it is the value that $D$ will take.

For $d \ge \sqrt{r_s^2 - r_e^2}$, the CCDF evaluated at such $d$ is the probability that there is no visible satellite, namely the no-satellite probability of Theorem 2. ■

<figure><img src="image_6.png" /><figcaption>Fig. 4. The no-satellite probability. The satellite altitude is 550 km. Simulation results validates the derived formula.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. The no-satellite probability. The satellite altitude is 550 km. Simulation results validates the derived formula.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. The no-satellite probability. The satellite altitude is 1100 km. Simulation results validates the derived formula.</figcaption></figure>

## V. COVERAGE PROBABILITY

This section evaluates the coverage probability of the typical user. We start with the Rayleigh fading case. Table I summarizes the system-level simulation parameters used for the SINR computation.

**Theorem 4.** In the interference-limited regime with $m = 1$, namely in Rayleigh fading environment, the coverage probability of the typical user is given by Eq. (16).

*Proof:* Let $X_*$ denote the location of the satellite nearest to the typical user at $n$, $I$ the power of the interference created by the visible satellites at distances greater than $\|X_* - n\|$, and $g$ the aggregate antenna gain when the transmit and receiver antennas are aligned to provide the maximum antenna gain, namely $g = g_1 + g_2$. We have

$$ \mathbf{P}(\text{SINR} > \tau) = \mathbf{E} \left[ \mathcal{L}_I \left( \frac{\tau \|X_* - n\|^\alpha}{pg} \right) \right], \quad (17) $$

where $\mathcal{L}_I(s)$ is the Laplace transform of $I$. We used here the fact that, for $m = 1$, $H$ follows an exponential random variable with mean one independent of $I$ and $X_*$. The difficulty with Eq. (17) is that the random variables $I$ and $\|X_* - n\|$ are not independent. In order to compute their joint distribution, we proceed in 3 main steps.
$$
\begin{equation}
\begin{split}
& \int_{r_s - r_e}^{\sqrt{r_s^2 - r_e^2}} \frac{\lambda \mu z}{\pi r_s r_e} \exp \left( -\lambda \int_{\xi(z)}^{\bar{\varphi}} \left( 1 - e^{-\frac{\mu}{\pi} \int_0^{\omega_2(\nu)} \left(1 - \mathcal{L}_H\left(\frac{\tau z^\alpha}{gK_{\nu,\omega}^\alpha}\right)\right) d\omega} \right) \cos(\nu) d\nu \right) \\
& \quad \times \exp \left( -\lambda \int_0^{\xi(z)} \left( 1 - e^{-\frac{\mu}{\pi} \operatorname{arcsin}\left(\sqrt{1-\cos^2(\xi(z)) \sec^2(\nu)}\right) - \frac{\mu}{\pi} \int_{\omega_1(\nu,z)}^{\omega_2(\nu)} \left(1 - \mathcal{L}_H\left(\frac{\tau z^\alpha}{gK_{\nu,\omega}^\alpha}\right)\right) d\omega \right) \cos(\nu) d\nu \right) \\
& \quad \times \left( \int_0^{\xi(z)} \frac{e^{-\frac{\mu}{\pi} \operatorname{arcsin}\left(\sqrt{1-\cos^2(\xi(z)) \sec^2(\varphi)}\right) - \frac{\mu}{\pi} \int_{\omega_1(\varphi,z)}^{\omega_2(\varphi)} \left(1 - \mathcal{L}_H\left(\frac{\tau z^\alpha}{gK_{\varphi,\omega}^\alpha}\right)\right) d\omega}{\sqrt{1-\cos^2(\xi(z)) \sec^2(\varphi)}} d\varphi \right) dz.
\end{split}
\tag{16}
\end{equation}
$$

The first step consists in introducing a partition of all
possibilities concerning D and the orbit that contains X. Let
i denote the index of the latter. We have

$$
\begin{align}
& \mathbf{E} \left[ \mathcal{L}_I \left( \frac{\tau \|X_* - n\|^\alpha}{pg} \right) \right] \tag{18} \\
&= \int_z \mathbf{E} \left[ \sum_i \mathcal{L}_I \left( \frac{\tau z^\alpha}{pg} \right) 1_{i_* = i} 1_{D \in [z, z+dz]} \right] \nonumber \\
&\stackrel{\text{(a)}}{=} \int_z \int_{(\theta, \phi)} \mathbf{E}_{\theta, \phi}^0 \left[ \mathcal{L}_I \left( \frac{\tau z^\alpha}{pg} \right) 1_{l_{i_*} = l(\theta, \phi)} 1_{D \in [z, z+dz]} \right] \nonumber \\
&\stackrel{\text{(b)}}{=} \lambda \int_z \int_{\phi} \mathbf{E}_{\theta, \phi}^0 \left[ \mathcal{L}_I \left( \frac{\tau z^\alpha}{pg} \right) 1_{\phi_{i_*} = \phi} 1_{D \in [z, z+dz]} \right] \sin(\phi) \, d\phi, \nonumber
\end{align}
$$

where $\mathbf{E}_{\theta, \phi}^0$ denotes the Palm probability of $\Xi$. In (a), we used
Campbell's formula. In (b), we used that fact that the function
of interest are invariant by a change of $\theta$. Note that the law
of $\Xi$ under the Palm probability in question is the law of
$\Xi + \delta_{(\theta, \phi)}$, with $\Xi$ distributed as above. Hence

$$
\begin{align}
& \mathbf{E}_{\theta, \phi}^{0} \left[ \mathcal{L}_I \left( \frac{\tau z^{\alpha}}{pg} \right) 1_{\phi_{i_*} = \phi} 1_{D \in [z, z+dz]} \right] \nonumber \\
&= \mathbf{E} \left[ \mathcal{L}_{\hat{I}} \left( \frac{\tau z^{\alpha}}{pg} \right) 1_{\phi_{\hat{i}_*} = \phi} 1_{\hat{D} \in [z, z+dz]} \right], \tag{19}
\end{align}
$$

with $\hat{I}$, $\hat{i}_*$ and $\hat{D}$ for $\Xi + \delta_{(\theta, \phi)}$ in place of $\Xi$.
Let $\mathcal{F}$ be the sigma-algebra generated by $\Xi$. We have

$$
\begin{align}
& \mathbf{E} \left[ \mathcal{L}_{\hat{I}} \left( \frac{\tau z^{\alpha}}{pg} \right) 1_{\phi_{\hat{i}_*} = \phi} 1_{\hat{D} \in [z, z+dz]} \right] \nonumber \\
&= \mathbf{E} \left[ \mathbf{E} \left[ \mathcal{L}_{\hat{I}} \left( \frac{\tau z^{\alpha}}{pg} \right) 1_{\phi_{\hat{i}_*} = \phi} 1_{\hat{D} \in [z, z+dz]} \mid \mathcal{F} \right] \right] \tag{20}
\end{align}
$$

and, by first principles,

$$
\begin{align}
& \mathbf{E} \left[ \mathcal{L}_{\hat{I}} \left( \frac{\tau z^{\alpha}}{pg} \right) 1_{\phi_{\hat{i}_*} = \phi} 1_{\hat{D} \in [z, z+dz]} \mid \mathcal{F} \right] \tag{21} \\
&= \mathbf{E} \left[ \mathcal{L}_{\hat{I}} \left( \frac{\tau z^{\alpha}}{pg} \right) \mid \phi_{\hat{i}_*} = \phi, \hat{D} \in [z, z+dz], \mathcal{F} \right] \nonumber \\
& \qquad \mathbf{P} \left[ \phi_{\hat{i}_*} = \phi, \hat{D} \in [z, z+dz] \mid \mathcal{F} \right]. \nonumber
\end{align}
$$

The next steps consist in computing the two terms in the last
expression. We start with the conditional Laplace transform
of the interference $\hat{I}$. As we now show, conditionally on $\mathcal{F}$
and the event $\hat{G}(z)$ that $\hat{\phi}_{\hat{i}_*} = \phi$ and $\hat{D} \in [z, z+dz]$, the

interference created by the satellites on different orbits are
conditionally independent and there is a closed form for the
Laplace transform of the interference on each orbit. To derive
this closed form, we use

$$
\xi(z) = \arccos \left( \frac{r_s^2 + r_e^2 - z^2}{2r_s r_e} \right). \quad (22)
$$

The angle $\xi(z)$ is that between the north direction and any
point on the rim of the spherical cap $C_z$ for $r_s - r_e < z <
\bar{d} = \sqrt{r_s^2 - r_e^2}$. Then we need the following

$$
K_{i,\omega} = \sqrt{r_s^2 - 2r_s r_e \sin(\omega) \sin(\phi_i) + r_e^2}, \quad (23)
$$

$$
K_{\omega} = \sqrt{r_s^2 - 2r_s r_e \sin(\omega) \sin(\phi) + r_e^2}, \quad (24)
$$

$$
\kappa_{i,1} = \pi/2 - \arcsin(\sqrt{1 - (r_e/r_s)^2 \csc^2(\phi_i)}), \quad (25)
$$

$$
\kappa_1 = \pi/2 - \arcsin(\sqrt{1 - (r_e/r_s)^2 \csc^2(\phi)}), \quad (26)
$$

$$
\kappa_{i,2}(z) = \pi/2 - \arcsin(\sqrt{1 - \cos^2(\xi(z))\csc^2(\phi_i)}), \quad (27)
$$

$$
\kappa_2(z) = \pi/2 - \arcsin(\sqrt{1 - \cos^2(\xi(z))\csc^2(\phi)}). \quad (28)
$$

The function $K_{i,\omega}$ gives the distance between the typical user and the point of orbit $l(\theta_i, \phi_i)$ with orbital angle $\omega$. The function $\kappa_{i,2}(z)$ gives the angle (in the orbital plane of orbit $i$, i.e., $l(\theta_i, \phi_i)$) between the apex of the orbit and the point of the orbit where the distance to $n$ is $z$. The parameter $\kappa_{i,1}$ is the angle (in the orbital plane of orbit $i$) between the apex of the orbit and the point of the orbit beyond which a satellite is no more visible from $n$.

Let $\mathcal{E}(z) = \mathcal{F} \cap \hat{\mathcal{G}}(z)$. For short we use the notation $\mathcal{E}$ for $\mathcal{E}(z)$. Let $N$ denote the number of points of $\xi$ and let $(\theta_i, \phi_i)$, $i = 1, \dots, N$, denote their coordinates. Let $\psi_i$ denote the linear Poisson point process on $l(\theta_i, \phi_i)$. Note that these random variables are $\mathcal{F}$-measurable. Using the conditional independence alluded to above, and the independence of the number of points of the linear Poisson point processes on disjoint subsets of an orbit, we get

$$
\mathcal{L}_{\hat{I}|\mathcal{E}}(s) = \mathbf{E}_{\psi} \left[ \prod_{X_j \in \psi}^{\bar{d} > K_{j,z}} \mathcal{L}_H(spK_j^{-\alpha}) \right] \\
\prod_{i=1,N}^{\bar{d} > K_{i,j} > z} \mathbf{E}_{\psi_i} \left[ \prod_{X_{j,i} \in \psi_i}^{\bar{d} > K_{i,j} > z} \mathcal{L}_H(spK_{i,j}^{-\alpha}) \mid \mathcal{E} \right].
$$
Then, we have

$$
\begin{align*}
\mathcal{L}_{\hat{I}|\mathcal{E}}(s) ={}& \exp \left( -\frac{\mu}{\pi} \int_{\kappa_1}^{\kappa_2(z)} (1 - \mathcal{L}_H(spK_{\omega}^{-\alpha})) \, d\omega \right) \\
& \prod_{|\phi_i - \pi/2| \le \xi(z)} \exp \left( -\frac{\mu}{\pi} \int_{\kappa_{i,1}}^{\kappa_{i,2}(z)} (1 - \mathcal{L}_H(spK_{i,\omega}^{-\alpha})) \, d\omega \right) \\
& \prod_{\xi(z) < |\phi_i - \pi/2| < \bar{\varphi}} \exp \left( -\frac{\mu}{\pi} \int_{\kappa_{i,1}}^{\pi/2} (1 - \mathcal{L}_H(spK_{i,\omega}^{-\alpha})) \, d\omega \right),
\end{align*}
$$

where $\mathcal{L}_H(s)$ denotes the Laplace transform of the mean 1 exponential random variable. Using the change of variables $\omega \to \pi/2 - \omega$, we get

$$
\begin{align}
\mathcal{L}_{\hat{I}|\mathcal{E}}(s) ={}& \exp \left( -\frac{\mu}{\pi} \int_{\omega_1(z)}^{\omega_2} (1 - \mathcal{L}_H(sp\widetilde{K}_{\omega}^{-\alpha})) \, d\omega \right) \nonumber \\
& \prod_{|\phi_i - \pi/2| < \xi(z)} \exp \left( -\frac{\mu}{\pi} \int_{\omega_{i,1}(z)}^{\omega_{i,2}} (1 - \mathcal{L}_H(sp\widetilde{K}_{i,\omega}^{-\alpha})) \, d\omega \right) \nonumber \\
& \prod_{\xi(z) < |\phi_i - \pi/2| < \bar{\varphi}} \exp \left( -\frac{\mu}{\pi} \int_{0}^{\omega_{i,2}} (1 - \mathcal{L}_H(sp\widetilde{K}_{i,\omega}^{-\alpha})) \, d\omega \right), \tag{29}
\end{align}
$$

where

$$
\omega_{i,1}(z) = \arcsin(\sqrt{1 - \cos^2(\xi(z)) \sec^2(\phi_i)}), \quad (30)
$$

$$
\omega_1(z) = \arcsin(\sqrt{1 - \cos^2(\xi(z)) \sec^2(\phi)}), \quad (31)
$$

$$
\omega_{i,2} = \arcsin(\sqrt{1 - (r_e/r_s)^2 \sec^2(\phi_i)}), \quad (32)
$$

$$
\omega_2 = \arcsin(\sqrt{1 - (r_e/r_s)^2 \sec^2(\phi)}), \quad (33)
$$

$$
\widetilde{K}_{i, \omega} = \sqrt{r_s^2 - 2r_s r_e \cos(\omega) \sin(\phi_i) + r_e^2}. \quad (34)
$$

Our third step consists in computing the following condi-
tional probability

$$
\mathbf{P}(\hat{G}(z) | \mathcal{F}) = \mathbf{P}(\exists j_{\star} \in \psi : \|X_{j_{\star}} - n\| \in [z, z + dz], \\
\forall j \neq j_{\star}, \|X_j - n\| > z, \forall i = 1, N, \forall j, \|X_{i,j} - n\| > z | \mathcal{F}).
$$

Using now the fact that the length of the arc $C_z \cap l(\theta_i, \phi_i)$ is
$2r_s \arcsin(\sqrt{1 - \cos^2(\xi(z)) \csc^2(\phi_i)})$, we get

$$
\begin{align}
\mathbf{P}(\hat{D} \in [z, z + dz), \phi_{\hat{i}_*} = \phi | \mathcal{F}) \tag{35} \\
&= \frac{\partial}{\partial z} \left( 1 - e^{-\frac{\mu}{\pi} \arcsin\left(\sqrt{1-\cos^2(\xi(z)) \csc^2(\phi)}\right)} \right) dz \nonumber \\
&\quad \prod_{\pi/2-\xi(z) < \phi_i < \pi/2+\xi(z)} \nonumber \\
&\qquad \qquad \qquad e^{-\frac{\mu}{\pi} \arcsin\left(\sqrt{1-\cos^2(\xi(z)) \csc^2(\phi_i)}\right)} \nonumber \\
&= \frac{\mu z |\csc(\phi)| e^{-\frac{\mu}{\pi} \arcsin\left(\sqrt{1-\cos^2(\xi(z)) \csc^2(\phi)}\right)}}{\pi r_e r_s \sqrt{1-\cos^2(\xi(z)) \csc^2(\phi)}} \nonumber \\
&\quad \prod_{\pi/2-\xi(z) < \phi_i < \pi/2+\xi(z)} \prod_{i=1,N} \nonumber \\
&\qquad \qquad \qquad e^{-\frac{\mu}{\pi} \arcsin\left(\sqrt{1-\cos^2(\xi(z)) \csc^2(\phi_i)}\right)}, \nonumber
\end{align}
$$

where cos(ξ(z)) is given by Eq. (22). Here, we used the fact
that conditionally on F, the satellite point processes on orbits
are independent Poisson point processes.

<table><caption>TABLE I
SIMULATION PARAMETERS</caption>
  <thead>
    <tr>
      <th>System parameter</th>
      <th>Value</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>Mean number of orbits &lambda;</td>
      <td>30</td>
    </tr>
    <tr>
      <td>Mean number of satellites per orbit &mu;</td>
      <td>30</td>
    </tr>
    <tr>
      <td>Mean number of total satellites &lambda;&mu;</td>
      <td>900</td>
    </tr>
    <tr>
      <td>Transmit antenna gain g</td>
      <td>20 dB</td>
    </tr>
    <tr>
      <td>Transmit antenna power p</td>
      <td>30 dB</td>
    </tr>
    <tr>
      <td>Transmitter EIRP (pg)</td>
      <td>80 dBm</td>
    </tr>
    <tr>
      <td>Satellite orbit radius r<sub>s</sub></td>
      <td>6950 km</td>
    </tr>
    <tr>
      <td>Boltzmann constant k</td>
      <td>-228.6 dBW/K/Hz</td>
    </tr>
    <tr>
      <td>Noise temperature T</td>
      <td>290 K</td>
    </tr>
    <tr>
      <td>Receive antenna gain g<sub>r</sub></td>
      <td>0 dB</td>
    </tr>
    <tr>
      <td>System bandwidth B<sub>w</sub></td>
      <td>30 MHz</td>
    </tr>
    <tr>
      <td>Carrier frequency f<sub>c</sub></td>
      <td>2 GHz</td>
    </tr>
  </tbody>
</table>

We now plug in Eq. (29) and (35) (together with (19), (20), and (21)) in (18) to get (16). This expression is obtained by using the Laplace functional of the Poisson point process Ξ and by the change of variable $\varphi = \pi/2 - \phi$. It uses the functions

$$
\begin{align*}
\overline{K_{\varphi, \omega}} &= \sqrt{r_s^2 - 2r_s r_e \cos(\omega) \cos(\varphi) + r_e^2}, \\
\omega_1(\varphi, z) &= \arcsin(\sqrt{1 - \cos^2(\xi(z)) \sec^2(\varphi)}), \\
\omega_2(\varphi) &= \arcsin(\sqrt{1 - (r_e/r_s)^2 \sec^2(\varphi)}).
\end{align*}
$$

This completes the proof.

<figure><img src="image_9.png" /></figure>

In a case that thermal noise cannot be neglected, the coverage probability of the typical user can be assessed by scaling the formula in Theorem 5 with a constant.

Theorem 5 gives the coverage probability of the typical user in terms of the mean number of orbits, the mean number of satellites per orbit, the radius of orbits, and the transmit antenna gain. As our Cox-distributed satellite model exhibits rotation invariance, the coverage probability of the typical user located at *n* statistically represent the coverage probabilities of all users in the network. Consequently, the derived coverage probability inherently characterizes the fundamental performance of the LEO or MEO satellite network. By utilizing the formula, one can readily forecast and estimate the behavior of the coverage probability by varying the geometric parameters.

In the following, we use the system-level parameter values in Table I, directly obtained and derived from values in [26], [30]. Specifically, we use a transmit antenna gain of $g = 20$ dB based on the satellite parameters provided in [26, Section 6] and [30, Section 6], where the maximum transmit antenna gains are equal to 30 dB or 24 dB.

Fig. 7 and 8 illustrate the coverage probability of the typical user for various values of $\lambda$ and $\mu$. The simulation results validate the accuracy of the derived formula in Theorem 5. While obtaining simulation results necessitates significant computation time due to the construction of large-scale network layout, the derived formula enables the generation of the entire coverage probability graph in much less time. This efficiency facilitates the execution of complex multivariate analyses of satellite networks, that we will see shortly. In the following, we explore the coverage behavior of the satellite network by varying various geometric parameters.

Fig. 9 displays the impact of satellite altitude on coverage probability. For the given $\lambda$ and $\mu$, it is observed that higher
$$ \int_{0}^{\infty} \int_{r_s-r_e}^{\sqrt{r_s^2-r_e^2}} \frac{\lambda \mu z}{\pi r_s r_e} \exp \left( -\lambda \int_{\xi(z)}^{\bar{\varphi}} \left( 1 - e^{-\frac{\mu}{\pi} \int_{0}^{\omega_2(\varphi)} \left( 1 - \mathcal{L}_H\left(\left(2^u-1\right)z^\alpha/(g\bar{K}_{\varphi,\omega}^\alpha)\right) \right)} d\omega \right) \cos(\varphi) d\varphi \right) \\ \times \exp \left( -\lambda \int_{0}^{\xi(z)} \left( 1 - e^{-\frac{\mu}{\pi} \arcsin\left(\sqrt{1-\cos(\xi(z))^2} \sec(\varphi)^2\right) - \frac{\mu}{\pi} \int_{\omega_1(\varphi,z)}^{\omega_2(\varphi)} \left( 1 - \mathcal{L}_H\left(\left(2^u-1\right)z^\alpha/(g\bar{K}_{\varphi,\omega}^\alpha)\right) \right)} d\omega \right) \cos(\varphi) d\varphi \right) \\ \times \left( \int_{0}^{\xi(z)} \frac{e^{-\frac{\mu}{\pi} \arcsin\left(\sqrt{1-\cos(\xi(z))^2} \sec(v)^2\right) - \frac{\mu}{\pi} \int_{\omega_1(v,z)}^{\omega_2(v)} \left( 1 - \mathcal{L}_H\left(\left(2^u-1\right)z^\alpha/(g\bar{K}_{v,\omega}^\alpha)\right) \right)}{\sqrt{1-\cos^2(\xi(z)) \sec^2(v)}} dv \right) dz du. \quad (36) $$

<figure><img src="image_10.png" /><figcaption>Fig. 7. The coverage probability for $\lambda$. We use $r_s = 7050$ km.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9. The coverage probability for various satellite altitudes.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 8. The coverage probability for various $\mu$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10. The 0-dB coverage probability for various $\lambda$ and $\mu$.</figcaption></figure>

satellite altitudes typically lead to lower coverage probabilities. It is worth noting that the difference in coverage probability from various satellite altitudes is noticeable for higher thresholds.

Fig. 10 shows the 0-dB coverage probability, namely the probability that the typical user has the SINR is greater than 0 dB. The 0-dB coverage probability of the typical user also indicates the total fraction of network users having their SINR greater than 0 dB. Such an interpretation leads to important insights for satellite network design or its optimization. For instance, for $\lambda = 50$ and $\mu = 50$, about 85% of users have

SINR greater than 0 dB. For $\lambda = 100$ and $\mu = 50$ or $\lambda = 50$ and $\mu = 100$, where the total average numbers of satellites are the same, their 0-dB coverage probabilities are almost identical, at 0.74. A deployment plan of forthcoming LEO or MEO satellite constellation—such as adding orbits or adding satellites on each orbit—can be effectively evaluated by adapting the proposed framework, offering a comprehensive tool to design LEO or MEO satellite networks for network operators.

Fig. 11 gives the 0-dB coverage probability for various $\mu$ and $g$. Here, we consider $g \in (15, 25)$ dB. In practice, LEO
<figure><img src="image_14.png" /><figcaption>Fig. 11. The 0-dB coverage probability for various $\mu$ and $g$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12. The 0-dB coverage probability for various altitude. We use $\mu = 40$.</figcaption></figure>

or MEO satellite network operators deploy more satellites to cope with an increased number of ground users. Yet, this may lead to increased interference or decreased coverage due to the additional satellites. Occupying additional orbital planes may not be a feasible option, and in this case, a viable solution may be to deploy new LEO or MEO satellites with better antenna gain. Using the derived formula, Fig. 11 shows how much antenna gain from the LEO or MEO satellite is required to achieve the same coverage probability level for various densities of $\mu$. For instance, if the satellite density $\mu$ changes from $\mu_1$ to $\mu_2$, the network operator may also need additional antenna gain from $g$ to $g'$ to ensure the same quality of service for network users. In the same vein, we conduct a similar experiment in Fig. 12 where we examine the coverage benefit of having additional antenna gain for satellite altitudes.

**Remark 3.** In the evaluation of the coverage probability in Theorem 5, we assume that the interference power dominates the noise power at the typical user. For the values in Table I, we have $\mathbf{E}[I] \gg kTB_w$ namely $\sum_{X_i \in \Psi \setminus X_*} pgr_i \|X_i\|^{-\alpha} \gg kTB_w$. For the values of $\lambda$ and $\mu$ in Table I, the average number of visible satellites is 35. The typical user is associated

with the nearest satellite whereas the remaining 34 satellites are creating interference to the typical user.

**Corollary 1.** Accounting for the noise at the user, the coverage probability is given by

$$ \mathbf{P}(\text{SINR} > \tau) = e^{-\frac{kTB_w\tau}{pg}} \mathbf{E} \left[ \mathcal{L}_I \left( \frac{\tau \|X_* - n\|^\alpha}{pg} \right) \right] \quad (37) $$

where the second expectation of Eq. (37) is given by Eq. (16). In other words, the coverage probability with thermal noise is equal to the expression of the SIR coverage probability of Eq. (37) multiplied by the following constant: $e^{-kTB_w\tau/(pg)}$.

The achievable rate or the ergodic capacity of downlink LEO or MEO satellite communications is given by Eq. (36).

*Proof:* Since the noise is a constant, we obtain the result by using Eq. (??). The achievable rate is

$$ \begin{aligned} \text{Throughput} &= \int_{0}^{\infty} \mathbf{P}(\log_{2}(1 + \text{SINR}) > r) \, dr \\ &= \int_{0}^{\infty} \mathbf{P}(\text{SINR} > 2^{r} - 1) \, dr, \end{aligned} \quad (38) $$

where we use the fact that $\log_2(1+SINR)$ is a positive random variable. ■

In below, we derive the coverage probability of the typical user for a general fading with $m = 1, 2, ...$

**Theorem 5.** In the interference-limited regime with general $m$, the coverage probability of the typical user is given by Eq. (39) where $\xi$, $\omega_{i,1}(z)$, $\omega_{i,2}$, $\tilde{K}_{i,\omega}$ are given by Eqs. (22), (30), (32), and (34), respectively. We also have

$$ \begin{aligned} h_1(z, \phi_i, s) &= \exp \left( -\frac{\mu}{\pi} \int_{\omega_{i,1}(z)}^{\omega_{i,2}} \left( 1 - \mathcal{L}_H \left( \frac{sp}{\tilde{K}_{i,\omega}} \right)^{\alpha} \right) d\omega \right), \\ h_2(z, \phi_i, s) &= \exp \left( -\frac{\mu}{\pi} \int_{0}^{\omega_{i,2}} \left( 1 - \mathcal{L}_H \left( \frac{sp}{\tilde{K}_{i,\omega}} \right)^{\alpha} \right) d\omega \right), \\ g(z, \phi) &= \exp \left( -\frac{\mu}{\pi} \arcsin \left( \sqrt{1 - \cos^2(\xi(z)) \csc^2(\phi)} \right) \right). \end{aligned} $$

*Proof:* As in the proof of Theorem 4, in the interference-limited regime, the coverage probability of the network is

$$ \begin{aligned} \mathbf{P}(\text{SINR} > \tau) &= \mathbf{E}_{\Xi} \left[ \mathbf{P} \left( H > \frac{\tau I \|X_* - n\|^\alpha}{pg} \middle| \Xi \right) \right] \\ &= \mathbf{E}_{\Xi, l_*, D, I} \left[ e^{-\frac{m \tau I z^\alpha}{pg}} \sum_{k=0}^{m-1} \frac{\frac{m^k \tau^k I^k z^{k\alpha}}{p^k g^k}}{k!} \middle| D, l_* \Xi \right] \right] \\ &= \sum_{k=0}^{m-1} \frac{1}{k!} \mathbf{E}_{\Xi, l_*, D} \left[ s^k (-1)^k \frac{\partial^k}{\partial s^k} \mathcal{L}_{I|D, l_*, \Xi}(s) \bigg|_{s=\frac{m\tau z^\alpha}{pg}} \right], \end{aligned} $$

where we used the probability distribution function of the Nakagami-$m$ random variable. We also used the fact that the $k$-th moment of the interference can be obtained by taking the $k$-order derivative of the Laplace transform of the interference. Hence, we have

$$ \mathbf{E}_I \left[ e^{-\frac{m \tau I z^\alpha}{pg}} \sum_{k=0}^{m-1} \frac{\frac{m^k \tau^k I^k z^{k\alpha}}{p^k g^k}}{k!} \middle| \mathcal{E} \right] = \sum_{k=0}^{m-1} \frac{(-s)^k}{k!} \frac{\partial^k}{\partial s^k} \mathcal{L}_{I|\mathcal{E}}(s), $$
$$
\sum_{k=0}^{m-1} \frac{\lambda}{2k!} \int_{r_s-r_e}^{\sqrt{r_s^2-r_e^2}} \int_{\pi/2-\xi}^{\pi/2+\xi} \mathbf{E}_{\Xi} \left[ \left\{ (-s)^k \frac{\partial^k}{\partial s^k} \left( h_1(z, v, s) \prod_{(\theta_i, \phi_i) \in \Xi} h_1(z, \phi_i, s) \prod_{(\theta_i, \phi_i) \in \Xi} h_2(z, \phi_i, s) \right) \right\} \frac{\mu ug(z, v) \prod_{(\theta_i, \phi_i) \in \Xi} g(z, \phi_i)}{\pi r_e r_s \sqrt{1 - \cos^2(\xi) \csc^2(v)}} \right] dv dz. \quad (39)
$$

where $s = \frac{m\tau z^\alpha}{pg}$. As in the proof of Theorem 4, we have

$$
\mathcal{L}_{I|\mathcal{E}}(s) = \exp\left(-\frac{\mu}{\pi} \int_{\omega_1(z)}^{\omega_2} \left(1 - \mathcal{L}_H(sp\tilde{K}_{\omega}^{-\alpha}) d\omega\right)\right) \\
\prod_{\substack{|\phi_i - \pi/2| < \xi(z) \\ i=1,N}} \exp\left(-\frac{\mu}{\pi} \int_{\omega_{i,1}(z)}^{\omega_{i,2}} \left(1 - \mathcal{L}_H(sp/\tilde{K}_{i,\omega}^{\alpha}) d\omega\right)\right) \\
\prod_{\substack{\xi(z) < |\phi_i - \pi/2| < \bar{\varphi} \\ i=1,N}} \exp\left(-\frac{\mu}{\pi} \int_{0}^{\omega_{i,2}} \left(1 - \mathcal{L}_H(sp/\tilde{K}_{i,\omega}^{\alpha}) d\omega\right)\right), \quad (40)
$$

where $\mathcal{L}_H(s)$ is the Laplace transform of the random variable $H$. The variables $\xi, \omega_{i,1}(z), \omega_{i,2}$, and $\tilde{K}_{i,\omega}$ are given by Eqs. (22), (30), (32), and (34).

Conditionally on $\Xi$ and $l_*$, the PDF of $D$ is given by taking
the derivative of the CDF of $D$ as follows:

$$
f_D(z) = \frac{\mu z |\csc(\phi)| e^{-\frac{\mu}{\pi} \arcsin(\sqrt{1-\cos^2(\xi(z)) \csc^2(\phi)})}}{\pi r_e r_s \sqrt{1 - \cos^2(\xi(z)) \csc^2(\phi)}} \prod_{\substack{|\phi_i - \pi/2| < \xi(z) \\ i=1,N}} e^{-\frac{\mu}{\pi} \arcsin(\sqrt{1-\cos^2(\xi(z)) \csc^2(\phi_i)})} . \quad (41)
$$

Finally, we get the result by combining all expressions and
following the same steps as in the proof of Theorem 4. ■

VI. DISCUSSION

A. Cox Model to a Forthcoming Constellation

In this section, we focus on describing the Starlink constel-
lation second-generation (2A) leveraging the proposed Cox
point process by comparing the SINR coverage probability
of the Starlink and the Cox model. Note that the Starlink
constellation 2A is a forthcoming plan and it is an example
among many future satellite constellations. Therefore, the aim
here is to demonstrate that the proposed model can be used to
model a planned constellation. Since the local distribution of
the Starlink constellation varies with user latitudes, we develop
a local approximation technique based on a moment matching.
Specifically, we adjust the orbital and satellite parameters,
denoted by $\lambda$ and $\mu$, to ensure that the Cox constellation and
the Starlink constellation have, on average, the same number
of satellites.

For the upcoming Starlink constellation, we refer to the
Starlink deployment plan available at the FCC [8], which
involves 28 orbital planes at an altitude of 525 km and an
inclination of 43 degrees, 28 orbital planes at 530 km with

an inclination of 53 degrees, and 28 planes at an altitude of
535 km with an inclination of 33 degrees. Each plane will
accommodate 120 satellites. To numerically derive the network
performance of this constellation, we consider a frequency
reuse factor of 4 where 30 satellites on each orbit use the
same spectrum resource.

Fig. 13 illustrates the coverage probability based on the upcoming Starlink constellation and the proposed Cox point process. Note that in order to numerically evaluate the SIR coverage probability of the typical user in the Starlink constellation, we simulate 10⁶ satellite layouts according to the available information above. Then, the SIR of the typical user is obtained for each layout and the collected SIR values are finally combined together as the CCDF of the SIR random variable. Based on the moment matching method specified above, we find the parameter settings (λ, μ) = (38, 80) for the latitude of 30 degrees and (λ, μ) = (100, 21) for the latitude of 0 degree. From Fig. 13, we observe that at the latitude of 30 degrees, the proposed Cox point process accurately emulates the constellation. Similarly, we observe that at a latitude of 0 degree, the proposed Cox point process replicates a coverage probability comparable to that of the Starlink constellation, with differences of less than 0.5 dB. The marginal difference arises from the geometric distinction that while the developed Cox point process is isotropic, the satellites of the forthcoming Starlink 2A constellation are regularly separated on each orbit and the orbits' inclinations are having only three values at 43, 53, 33 degrees.

<figure><img src="image_16.png" /><figcaption>Fig. 13. The coverage probability of Starlink and Cox constellations for users at latitudes of 0 and 30 degrees. $r_s = 7000$.</figcaption></figure>
## B. Cox model and Binomial Model

Here, we highlight the difference between our Cox point process and the binomial point process. Analytical stochastic geometry models, such as the binomial model [11]–[16] or our Cox model, approximate existing or forthcoming target constellations. For example, the binomial point process features a geometric parameter, denoted as $N$, representing the number of satellites. This parameter $N$ is then fine-tuned to match the mean number of visible satellites of the target constellation. In contrast, the proposed Cox point process features two parameters, $\lambda$ and $\mu$, and these two parameters are adjusted jointly to match the mean number of orbits and the mean number of satellites at the same time. Consequently, the presence of two parameters allows one to approximate the higher-order geometric characteristics of the target constellation, offering a better representation of orbital planes that LEO or MEO satellite constellations feature in space.

It is important to note that a proper comparison must be based on various mathematical properties of those point processes such as the repulsion between points often examined by the Ripley-K function of the point process. The detailed comparison between binomial and Cox is left for future work. Here, we simply study them by underlining their differences in performance metrics such as the no-satellite probability or the coverage probability.

First, we compare the no-satellite probability of the proposed Cox point process to the no-satellite probability of the binomial point process. Consider the total numbers of satellites to be $N = 100, 200, 400$, and $1000$, all at the altitudes of $550$ km. For the binomial model, we empirically compute the no-satellite probability by first creating the binomial satellite point process and then counting each event that there is no satellite point visible from the typical observer $(0, 0, r_e)$. We obtain the no-satellite probabilities for $N = 100, 200, 400$, and $1000$ as follows: 0.0176 for $N = 100$ and 0.0003 for $N = 200$. For $N = 400$ and $N = 1000$, the no-satellite probabilities are negligible.

On the other hand, the no-satellite probability of the proposed Cox point process is given by Eq. (14). For a mean of 100 satellites, we basically have numerous ways to distribute the mean of 100 satellites by varying the parameters $\lambda$ and $\mu$ together, namely $100 = \lambda\mu$. For instance, we can use the following 10 cases: $(\lambda, \mu) = (10, 100/10), (20, 100/20), \dots$, and $(100, 100/100)$. The first case indicates that there are 100 satellites distributed over 10 distinct orbits whereas the last case indicates the 100 satellites are distributed over 100 distinct orbits. For the above 10 cases, we numerically obtain the following no-satellite probabilities: $0.0886, 0.0471, \dots$, and $0.0235$. It is worth noting that the first and last cases created by the Cox model are topologically very different because of the existence of the orbit parameter $\lambda$. We observe that the proposed Cox model can reproduce an accurate no-satellite probability relevant to each deployment scenario in practice, by delineating various deployment scenarios such as a dense orbit scenario with few orbits and many satellites per orbit or a sparse orbit scenario with many orbits and few satellites per orbit. Note, on the other hand, that the binomial

model produces a single value of 0.0176.

Similarly, for a mean of 200 satellites, the no-satellite probability of the binomial model is very low at 0.0003. On the other hand, for the proposed Cox model, we can vary $(\lambda, \mu) = (10, 200/10), (20, 200/20), \dots$, and $(200, 200/200)$. The no-satellite probabilities are given by 0.033, 0.0079, \dots, and 0.0006. Consider a deployment scenario where 200 satellites are uniformly distributed over 10 isotropic orbits. The proposed Cox model produces an accurate estimate of the no-satellite probability as 0.033, whereas the binomial model produces an estimate of the no-satellite probability value as 0.0003, which is twice larger than the extreme case that the Cox model produces with even 200 orbits. This low no-satellite probability of the binomial model stems from the fact that the points of the binomial point process are independent and uniformly distributed over the sphere. Therefore, it is statistically enforced to produce an extremely low no-satellite probability. In conclusion, unlike the existing binomial point process where only a single layout is created for a given number of satellites, the developed Cox model can reproduce various deployment scenarios.

Secondly, in Fig. 14 we numerically compare the SIR coverage probabilities of both the binomial and Cox satellite point processes. We consider a single target constellation having 300 satellites at an altitude of 550 km and the rest of simulation parameters are given in Table I. The modeling technique based on the binomial point process results in a single curve for the coverage probability of the typical user. On the other hand, the modeling based on Cox point process has additional flexibility through the number of orbits. This allows for a wide range of achievable coverage probabilities by varying $\lambda$ and $\mu$. Thanks to the flexibility, we observe that even when the two point processes have the same mean number of satellites, the achievable region of the coverage probabilities derived under the Cox point process may include the coverage probability derived under the binomial point process. In Fig. 14, the shaded area is the coverage probabilities produced by the Cox model for a single mean number of satellites. For this mean number, the Cox point process produces a continuum of orbital planes by varying the orbit parameter $\lambda$ from a very small number to a very large number. From the wide range of coverage probabilities that the Cox model produces, Fig. 14 suggests that the Cox model may include the binomial model in terms of the coverage probability.

## C. Direct Orbit Representation

In this paper, we consider a non-homogeneous density function $\lambda \sin(\phi)/(2\pi)$ to create the orbits that are isotropic and invariant by rotations. Some forthcoming constellations may exhibit regularly spaced orbits where orbits' longitudes are regularly spaced while their inclinations are fixed to some values. For such cases, we modify the density function to directly represent the constellation. For instance, the satellite constellation with regularly spaced $N$ orbits with an inclination angle of 43 degrees (0.75 rad) are represented as follows: $\Xi = \sum_{i=1}^{N} \delta_{\pi/N(i-1)+\bar{U}, 0.75}$ where $\bar{U}$ is a uniform random variable in the interval $(0, \pi/N)$. Each point of $\Xi$
<figure><img src="image_17.png" /><figcaption>Fig. 14. Our Cox model features a wide range of coverage probabilities whereas a single coverage probability curve appears for the binomial model.</figcaption></figure>

represents an orbit of inclination 0.75 rad and of longitude $\{(i-1)\pi/N\}_{i=1,...,N}$. The uniform random shift is required in the above expression since it ensures that the longitudinally rotated version of $\Xi$ and its original $\Xi$ have the same distribution. In other words, $\Xi$ is longitudinally rotation invariant. The above representation exactly captures the regularly spaced single constellation, and therefore, leveraging the machinery that we explored in this paper, one can also derive the no-satellite probability or the SIR coverage probability of the regularly spaced satellite network.

Another example of regularly spaced LEO or MEO satellite constellation is a Walker constellation. A simple Walker constellation exhibits $N$ orbits each with the inclination of $\pi/2$ degrees; $\tilde{\Xi} = \sum_{i=1}^{N} \delta_{\frac{\pi}{N}(i-1)+\tilde{U}, \pi/2}$ where $\tilde{U}$ is a uniform random variable in the interval $(0, \pi/N)$.

Nevertheless, it is essential to note that since the above orbit processes are longitudinally rotation invariant, the probability law of the orbit process will alter as we rotate the orbit process along any axes other than the z-axis. In other words, observers at different latitudes will see different satellite distributions, and therefore, the no-satellite probability, the interference, and the SIR coverage probability must depend on the latitude of the typical observer. A more detailed analysis based on the modified orbit process is left for future work.

## VII. CONCLUSION

We have constructed an analytical framework for modeling and analyzing downlink LEO or MEO satellite communications. By developing an isotropic Cox point process for LEO or MEO satellites, we introduced a novel method for incorporating the essential geometry of LEO or MEO satellite constellations, where satellites are always located on orbits. We calculate the no-satellite probability for the typical user. Assuming that network users receive their downlink signals from their closest satellites, we evaluate the Laplace transform of the interference and assess the coverage probability for the typical user. The derived metrics represent network performance as functions of key network parameters. Thanks

<figure><img src="image_18.png" /><figcaption>Fig. 15. Illustration of spherical caps $C_d$ and $C_{\bar{d}}$ where $r_s - r_e < d < \bar{d}$ where $\bar{d}$.</figcaption></figure>

to isotropy, our analysis of the typical user's network performance can be interpreted as the network performance spatially averaged over all users across the network. Therefore, our analysis can be applied to the systematic design of large-scale downlink satellite communication networks. We demonstrate that our Cox point process effectively replicates a forthcoming LEO or MEO satellite constellation through moment matching approximation. Additionally, our Cox model provides numerous ways to model real or forthcoming constellations, in contrast to models based on binomial or Poisson point processes.

Future work will concentrate on addressing the limitations of the proposed model. For example, our framework currently assumes isotropic orbits, and a similar approach can be used to explore and investigate non-isotropic orbits. Similarly, the developed isotropic orbit may feature regularly-spaced satellites exclusively located on each orbit. Moreover, the proposed model assumes a fixed altitude for satellites. A similar approach can also be used to analyze satellite networks at various altitudes.

## APPENDIX A

Let us treat $d$ as a constant between $r_s - r_e$ and $\sqrt{r_s^2 - r_e^2}$ for the moment. Fig. 15 shows the bottom of spherical cap $C_d$ where the arc $\overline{BF}$ is defined as $C_d \cap l(\theta, \varphi)$.

For a triangle $\triangle AOU$, in Fig. 15 $\overline{AO} = r_s$ and let $\xi$ be the angle $\angle AOU$. Then, we have $\cos(\xi) = \frac{r_s^2+r_e^2-d^2}{2r_e r_s}$.

For $\triangle BOC$, $\angle BOC = \xi$ and we have $\overline{OB} = r_s$, $\overline{OC} = r_s \cos(\xi)$, $\overline{BC} = r_s \sin(\xi)$. Then, since the orbital plane of $l(\theta, \phi)$ has the inclination of $\phi$, $\angle DOC = \pi/2 - \phi$, $\overline{OC} = r_s \cos(\xi)$, and $\overline{CD} = r_s \cos(\xi) \cot(\phi)$.

For $\triangle BCD$, we have $\overline{BC} = r_s \sin(\xi)$ and $\overline{CD} = r_s \cos(\xi) \cot(\phi)$. Since $\angle BDC = \pi/2$, we have $\overline{BD} = \sqrt{r_s^2 \sin^2(\xi) - r_s^2 \cos^2(\xi) \cot^2(\phi)}$.

Now consider a triangle $\triangle OBD$ in the orbital plane. We have $\overline{OB} = r_s$ and $\overline{BD} = \sqrt{r_s^2 \sin^2(\xi) - r_s^2 \cos^2(\xi) \cot^2(\phi)}$. Let $\kappa' = \angle BOD$. Since $\angle BDO = \pi/2$, we get

$$ \sin(\kappa') = \frac{\overline{BD}}{r_s} = \sqrt{\sin^2(\xi) - \cos^2(\xi) \cot^2(\phi)}. \quad (42) $$

As a result, for $\pi/2 - \xi < \phi < \pi/2 + \xi$ the length of the arc $\overline{BF}$ is $2r_s \arcsin(\sqrt{1 - \cos^2(\xi) \csc^2(\phi)})$.
APPENDIX B

Based on the simple geometry, the (x, y, z) coordinates of
the satellite with orbital angle ωⱼ on the orbit l(θᵢ, φᵢ) are

$$
\begin{align*}
x &= \sqrt{r_s^2 \cos^2(\omega_j) + r_s^2 \sin^2(\omega_j) \cos^2(\phi_i)} \cos(\tilde{\theta}_i + \theta), \\
y &= \sqrt{r_s^2 \cos^2(\omega_j) + r_s^2 \sin^2(\omega_j) \cos^2(\phi_i)} \sin(\tilde{\theta}_i + \theta), \\
z &= r_s \sin(\omega_j) \sin(\phi_i), \\
\tilde{\theta} &= \arctan(\tan(\omega_j) \cos(\phi_i)).
\end{align*}
$$

Moreover, the distance from X_{j,i} to n is given by

$$
\|X_{j,i} - n\| = \sqrt{r_s^2 - 2r_s r_e \sin(\omega_j) \sin(\phi_i) + r_e^2} = K_{i,j}. \quad (43)
$$

ACKNOWLEDGMENT

This paper is a joint study with the support of Korea-France
Cooperative Development Program (STAR) from the NRF,
funded by the Korean government in the year 2024 (RS-2024-
00247682). The work of Chang-Sik Choi was supported by
the NRF-RS-2024-00334240. The work of Francois Baccelli
was supported by ERC NEMO grant 788851 to INRIA and
by the French National Agency for Research (ANR) via the
project n°ANR-22-PEFT-0010 of the France 2030 program
PEPR. Réseaux du Futur.

REFERENCES

[1] Z. Qu, G. Zhang, H. Cao, and J. Xie, "LEO satellite constellation for Internet of Things," *IEEE Access*, vol. 5, pp. 18391–18401, 2017.

[2] Y. Su, Y. Liu, Y. Zhou, J. Yuan, H. Cao, and J. Shi, "Broadband LEO satellite communications: Architectures and key technologies," *IEEE Wireless Commun.*, vol. 26, no. 2, pp. 55–61, 2019.

[3] Iridium Communications, Inc, "History of Iridium," June 14, 2024. [Online]. Available: https://www.iridiummuseum.com/timeline/

[4] Q. Chen, G. Giambene, L. Yang, C. Fan, and X. Chen, "Analysis of inter-satellite link paths for LEO mega-constellation networks," *IEEE Trans. Veh. Technol.*, vol. 70, no. 3, pp. 2743–2755, 2021.

[5] Q. Chen, L. Yang, Y. Zhao, Y. Wang, H. Zhou, and X. Chen, "Shortest path in LEO satellite constellation networks: An explicit analytic approach," *IEEE J. Sel. Areas Commun.*, vol. 42, no. 5, pp. 1175–1187, 2024.

[6] A. Guidotti, A. Vanelli-Coralli, M. Conti, S. Andrenacci, S. Chatzinotas, N. Maturo, B. Evans, A. Awoseyila, A. Ugolini, T. Foggi, L. Gaudio, N. Alagha, and S. Cioni, "Architectures and key technical challenges for 5G systems incorporating satellites," *IEEE Trans. Veh. Technol.*, vol. 68, no. 3, pp. 2624–2639, 2019.

[7] 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, 2021.

[8] Federal Communications Commission, "Kuiper systems, LLC application for authority to deploy and operate a Ka-band non-geostationary satellite orbit system," no. SAT-LOA-20190704-00057, Jul 30 2020. [Online]. Available: https://www.fcc.gov/document/fcc-authorizes-kuiper-satellite-constellation

[9] ——, "FCC authorizes Boeing broadband satellite constellation," Nov. 3 2021. [Online]. Available: https://www.fcc.gov/document/fcc-authorizes-boeing-broadband-satellite-constellation

[10] J. Fomon. Starlink slowed in Q2, competitors mounting challenges. [Online]. Available: https://www.ookla.com/articles/starlink-hughesnet-viasat-performance-q2-2022

[11] 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, 2020.

[12] A. Talgat, M. A. Kishk, and M.-S. Alouini, "Nearest neighbor and contact distance distribution for binomial point process on spherical surfaces," *IEEE Commun. Lett.*, vol. 24, no. 12, pp. 2659–2663, 2020.

[13] —, "Stochastic geometry-based analysis of LEO satellite communication systems," *IEEE Commun. Lett.*, vol. 25, no. 8, pp. 2458–2462, 2021.

[14] D.-H. Na, K.-H. Park, Y.-C. Ko, and M.-S. Alouini, "Performance analysis of satellite communication systems with randomly located ground users," *IEEE Trans. Wireless Commun.*, vol. 21, no. 1, pp. 621–634, 2022.

[15] 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, 2022.

[16] 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, 2023.

[17] S. N. Chiu, D. Stoyan, W. S. Kendall, and J. Mecke, *Stochastic geometry and its applications*. John Wiley & Sons, 2013.

[18] F. Baccelli and B. Blaszczyszyn, "Stochastic geometry and wireless networks: volume I theory," *Foundations and Trends in Networking*, vol. 3, no. 3–4, pp. 249–449, 2010.

[19] —, "Stochastic geometry and wireless networks: Volume II applications," *Foundations and Trends in Networking*, vol. 4, no. 1–2, pp. 1–312, 2010.

[20] J. G. Andrews, F. Baccelli, and R. K. Ganti, "A tractable approach to coverage and rate in cellular networks," *IEEE Trans. Commun.*, vol. 59, no. 11, pp. 3122–3134, 2011.

[21] M. Haenggi, *Stochastic geometry for wireless networks*. Cambridge University Press, 2012.

[22] C.-S. Choi and F. Baccelli, "Poisson Cox point processes for vehicular networks," *IEEE Trans. Veh. Technol.*, vol. 67, no. 10, pp. 10160–10165, Oct 2018.

[23] C.-S. Choi, F. Baccelli, and G. de Veciana, "Densification leveraging mobility: An IoT architecture based on mesh networking and vehicles," in *Proc. IEEE/ACM MobiHoc*, 2018, p. 71–80.

[24] V. V. Chetlur and H. S. Dhillon, "Coverage analysis of a vehicular network modeled as Cox process driven by Poisson line process," *IEEE Trans. Wireless Commun.*, vol. 17, no. 7, pp. 4401–4416, July 2018.

[25] C.-S. Choi and F. Baccelli, "An analytical framework for coverage in cellular networks leveraging vehicles," *IEEE Trans. Commun.*, vol. 66, no. 10, pp. 4950–4964, Oct 2018.

[26] 3GPP TR 38.821, "Solutions for NR to support non-terrestrial networks (NTN)," 3GPP TR 38.821.

[27] M. Arti, "Two-way satellite relaying with estimated channel gains," *IEEE Trans. Commun.*, vol. 64, no. 7, pp. 2808–2820, 2016.

[28] G. Zheng, S. Chatzinotas, and B. Ottersten, "Generic optimization of linear precoding in multibeam satellite systems," *IEEE Trans. Wireless Commun.*, vol. 11, no. 6, pp. 2308–2320, 2012.

[29] B. Pavan Kumar, C. Kumar, V. Senthil Kumar, and V. V. Srinivasan, "Active spherical phased array design for satellite payload data transmission," *IEEE Trans. Antennas and Propagation*, vol. 63, no. 11, pp. 4783–4791, 2015.

[30] 3GPP TR 38.821, "Study on new radio (NR) to support non-terrestrial networks)," 3GPP TR 38.811.

[31] M. E. Muller, "A note on a method for generating points uniformly on n-dimensional spheres," *Commun. of the ACM*, vol. 2, no. 4, pp. 19–20, 1959.