{'fraction_non_alphanumeric': 0.029060250689508262, 'fraction_numerical': 0.008744982846379801, 'mean_word_length': 4.675489946551285, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Downlink Coverage and Rate Analysis of Low Earth Orbit Satellite Constellations Using Stochastic Geometry

As low Earth orbit (LEO) satellite communication systems are gaining increasing popularity, new theoretical methodologies are required to investigate such networks' performance at large. This is because deterministic and location-based models that have previously been applied to analyze satellite systems are typically restricted to support simulations only. In this paper, we derive analytical expressions for the downlink coverage probability and average data rate of generic LEO networks, regardless of the actual satellites' locality and their service area geometry. Our solution stems from stochastic geometry, which abstracts the generic networks into uniform binomial point processes. Applying the proposed model, we then study the performance of the networks as a function of key constellation design parameters. Finally, to fit the theoretical modeling more precisely to real deterministic constellations, we introduce the effective number of satellites as a parameter to compensate for the practical uneven distribution of satellites on different latitudes. In addition to deriving exact network performance metrics, the study reveals several guidelines for selecting the design parameters for future massive LEO constellations, e.g., the number of frequency channels and altitude.

simulated thoroughly before putting forward, general theoretical understanding on the performance of massive LEO communication constellations at large is still missing in the scientific literature.
In this paper, we apply stochastic geometry to acquire analytical tractable expressions for coverage probability and average data rate of downlink LEO networks. The approach formulated herein for the first time ever to the authors' knowledge paves the way to study the generic performance of satellite networking without relying on explicit orbit simulations and the actual geometry of any specific constellation.

A. Related Works
In [1], the uplink outage probability in the presence of interference was evaluated for two LEO constellations through time-domain simulations. A performance study of Iridium constellation was presented in [2] in terms of system capacity, the average number of beam-to-beam handoffs and satellite-tosatellite handoffs, the channel occupancy distribution and average call drop probability. The effect of traffic non-uniformity on signal-to-interference ratio was studied in [3] by assuming hexagonal service areas for satellites. A teletraffic analysis of a mobile satellite system based on a LEO constellation was performed in [4]. A general expression for a single LEO satellite's visibility time is provided in [5], but it is incapable of concluding the general distribution of visibility periods for any arbitrarily positioned user.
Stochastic geometry is a powerful mathematical and statistical tool for the modeling, analysis, and design of wireless networks with irregular topologies [6]- [8]. It has been mostly used in the literature to analyze two-dimensional (planar) terrestrial networks [6]- [14]. In [9], a comprehensive review of the literature related to the stochastic geometry modeling of multi-tier and cognitive cellular networks was presented. In [10], the Poisson point process (PPP) was observed to provide lower bounds on the coverage probability and the average transmission rate which are as tight as an upper bound provided by an idealized grid-based model. The work in [10] was extended to multi-tier networks under quality-of-service constraints in [11], [12]. Uplink coverage when base stations and users follow independent PPPs is evaluated in [13].
While stochastic geometry has been extensively applied for the analysis of planar scenarios, its application to threedimensional networks has started attracting significant attention recently [15], [16]. In [15], modeling and analysis of coverage in three-dimensional cellular networks have been investigated using a PPP model. Although PPP has provided tractable and insightful results, it is not valid for modeling a finite-area network with a limited number of nodes [17]. For such cases, binomial point process (BPP) is an appropriate model to capture the characteristics of the network [18], [19]. The performance for an arbitrarily shaped planar network is studied in [18], [20], but the reference transmitter is assumed to be located at a given distance, which can be overcome with two transmitter selection policies from [14]. A BPP modeling for a finite three-dimensional network of unmanned aerial vehicles (UAVs) was developed in [16], [21].
Tools from stochastic geometry have been applied to study satellite communications [22]- [25] only with a limited extent. For instance, in [22], stochastic geometry was used to model the locations of terrestrial users involving multi-beam satellite and terrestrial interference from cellular base stations. The authors in [23] derived coverage probability and data rate of a multi-UAV downlink network through assuming PPP distribution for users where, however, no interference among users and UAVs is considered due to their channel assignment policy. In [24], the performance of a cognitive satellite-terrestrial network is investigated; the secondary terrestrial network and users are modeled as independent point processes and share resources with a primary satellite system. An analysis of coverage times during LEO satellite visits has been conducted in [25] by inclusion of the distribution of users' positions.
We can conclude from the above that the utilization of stochastic geometry for satellite network analysis is limited to modeling the user's locality with a low number of satellites. Alternative methods to stochastic geometry are unable to provide a comprehensive analysis that fits to any arbitrary constellation due to being associated with some specific network design parameters. In these models, coverage footprints of satellites are assumed to be identical and typically form a regular circular or hexagonal grid, although uneven distribution of satellites along different latitudes, and differences in transmitted power and/or altitude, create irregular cells. Moreover, there is no analytical tool to model interference in a generic sense. Figure 1(a) demonstrates a random constellation in which a set of satellites is distributed on a sphere according to a uniform point process, while Figs. 1(b)-(d) show a class of regular deterministic Walker constellations with varying inclination angle, in which all satellites are evenly spaced and have the same period and inclination. Each cell represents the coverage area of a satellite wherein it is the closest to and, thus, serving all users located inside the polygon. The set of all coverage cells form a Voronoi tessellation which we can observe to be analogous in the zoomed regions for the random constellation and deterministic Walker constellations. While selecting the inclination angle in practice is based on the service area of interest, we see that smaller inclination results in a more irregular Voronoi tessellation. Furthermore, the inter-operation of multiple LEO constellations, as it might hold in the near future, will make satellites' mutual positions even more similar to those given by a random process.

B. Contributions and Organization of the Paper
Motivated by the above observations, in this paper, we first model the satellite constellation with an appropriate point process, which will then allow us to utilize the tools from stochastic geometry to analyze the performance of generic LEO networks in theory. We consider a network of a given number of satellites whose locations are modeled as a BPP on a sphere at a fixed altitude, which is justifiable due to the limited number of satellites covering a given finite region [26]. Users are located at some arbitrary locations on Earth and are associated with the nearest satellite while some other satellites above the horizon of a user can cause co-channel interference due to frequency reuse. However, since satellites in Walker constellations are distributed unevenly along different latitudes, i.e., the number of satellites is effectively larger on the inclination limit of the constellation than on the equatorial regions, the density of practical deterministic constellations is typically not uniform like with BPP modeling. Thus, we apply a new parameter in order to compensate for the uneven density w.r.t. practical Walker constellations and create a tight match between the results generated by BPP modeling and those from practical constellation simulations.
Based on the modeling summarized above, we present the following scientific contributions in this paper.
• We derive exact expressions for coverage probability and average achievable data rate of a user in terms of the Laplace transform of interference power distribution. • We validate our novel theoretical results with numerical simulations and also compare them with reference results from actual deterministic satellite constellations. • To suppress the performance mismatch between a random network and practical constellations, which mostly stems from uneven distribution of satellites along different latitudes, we define and calculate a new parameter, the effective number of satellites, for every user latitude. • We show that, with the above compensation, the generic performance of large deterministic constellations can be very accurately analyzed with theoretical expressions that are based on stochastic constellation geometry. • Finally, the two objectives of coverage and data rate are evaluated for different key design parameters, e.g., the number of frequency channels and satellite altitude.
As for propagation models, we consider two extreme cases, namely Rayleigh fading and static propagation, for the serving channels. The former corresponds to a more drastic fading environment when the received signal is subject to severe multi-path distortion due to the small elevation angle of the transmitting satellite; this case leads to simpler expressions for some specific path loss exponents. The latter represents the Fig. 2. A sketch of the considered system's stochastic geometry, where satellites are distributed randomly on a sphere with radius r ⊕ + r min and a user is located on the surface of another cocentric sphere with radius r ⊕ . typical cases where the serving satellite's elevation is large enough to provide line-of-sight to the user and, thus, weak components from multi-path propagation become insignificant. For interfering channels, any fading statistics can be adopted in general since this has no effect on analytical tractability.
We see that, although increasing the number of frequency bands improves the coverage probability, there is an optimal number of frequency channels that maximizes the data rate depending on the path loss exponent. Assuming a constant effective number of satellites for different altitudes, we observe that the optimum height which maximizes coverage probability or data rate is not within the practical altitude range of LEO networks (i.e., outside Earth's atmosphere), where the performance always declines with increasing the altitude.
The organization of the remainder of this paper is as follows. Section II describes the system model for a randomly distributed satellite network and characterizes some baseline probabilities stemming from the stochastic geometry of the system. As for the main results, we derive analytical expressions for downlink coverage probability and average achievable data rate for a terrestrial user in Sections III and IV, respectively. Numerical results are provided in Section V for studying the effect of key system parameters such as satellite altitude and the number of frequency bands allocated for the network. Finally, we conclude the paper in Section VI.

II. SYSTEM MODEL
Let us consider a downlink network of N satellite base stations that are uniformly distributed around Earth at the same altitude r min forming a BPP on a sphere with radius r ⊕ + r min , where r ⊕ denotes Earth's radius, as shown in Fig. 2, while user terminals are located on the surface of Earth. The altitude parameter r min specifies also the minimum possible distance from a satellite to a user (that is realized when it is directly above the user), hence the name. We assume that wireless transmissions propagate to a user from all and only the satellites that are above its horizon. Correspondingly, r max = 2r ⊕ r min + r 2 min denotes the maximum possible distance at which a satellite has any effect on the network service to a user (that is realized when the satellite is at the horizon). The notation followed in this paper is summarized in Table I.
Each user is associated to the nearest satellite that is referred to as the serving satellite in what follows, resulting in spherical Voronoi tessellation for satellites' coverage areas as illustrated in Fig. 1(a). The other satellites may cause cochannel interference to the user because we assume that there are only K, with K ≤ N, orthogonal frequency channels for the network and randomly assign a subset of N/K satellites to each channel. The scheduling performed in this way ensures that the nearest satellite in the constellation uses the channel that is assigned to a user. All the other satellites on the same channel, other than the serving satellite, cause interference to the user's reception whenever they are above the horizon.
We suppose that every satellite may be equipped with a directional antenna that radiates its main lobe towards the center of Earth. From the serving satellite, some lobe containing the higher power (but not necessarily the main lobe) is likely directed to the user, while the lobes with lower power levels radiate towards the user from the interfering satellites. In order to approximate the effect of directional transmission, we set differing power levels for serving and interfering satellites which are denoted by p s and p i , respectively, such that p i ≤ p s . If the user is located within the beamwidth of the serving satellite, the transmitted power from it will have the main lobe power level. Considering the main lobe transmitted power from the serving satellite and the largest sidelobe transmitted power from the interfering satellites, p s /p i corresponds to the sidelobe level of the antennas.
The distances from the user to the serving satellite and the other satellites are denoted by random variables R 0 and R n , n = 1, 2, . . . , N − 1, respectively, while G 0 and G n represent the corresponding channel gains. Obviously, G n = 0 if R n > r max for some n = 0, 1, . . . , N − 1, i.e., the satellite is below horizon. For notational convenience, when N I > 0, we let indices n = 1, 2, . . . , N I correspond to those N I ≤ N/K − 1 satellites (if any) that share the same frequency channel with the serving satellite and are above the user's horizon, so that they cause co-channel interference.
Based on the above modeling, the signal-to-interferenceplus-noise ratio (SINR) at the receiver can be expressed as where we assume that the user's receiver is subject to additive noise with constant power σ 2 , the parameter α is a path loss exponent, is the cumulative interference power from all other satellites above the user's horizon than the serving satellite, and N I is a random variable denoting the number of interfering satellites.
In the special case of having no interfering satellites, i.e., N I = 0, since they happen to be below the user's horizon, the SINR in (1) is reduced to signal-to-noise ratio (SNR) as and further SNR = 0 if R 0 > r max , i.e., also the nearest satellite that is supposed to be the serving satellite is below horizon. In order to contribute expressions for coverage probability and average achievable rate in the following sections, we first need to characterize some basic distance distributions that stem from the stochastic geometry of the considered system. In particular, we express the necessary cumulative distribution function (CDF) and probability density functions (PDFs) in the following three lemmas. Lemma 1. The CDF of the distance R from any specific one of the satellites in the constellation to the user is given by and the corresponding PDF is given by The above distribution of R will be a key tool in deriving other distance statistics (i.e., those of the serving and interfering satellites), which are required in turn to characterize the distribution of SINR in (1) using stochastic geometry.
Lemma 2. The PDF of the serving distance R 0 is given by Proof. Due to the channel assignment by which the serving satellite is the nearest one among all the N i.i.d. satellites, the CDF of R 0 can be expressed as ) N with the substitution of (4). The corresponding PDF is obtained by differentiation to complete the proof.
Lemma 3. When conditioned on the serving distance such that R 0 = r 0 , the PDF of the distance from any other satellite to the user is given by for r 0 < r n ≤ 2r ⊕ + r min while f R n |R 0 (r n |r 0 ) = 0 otherwise.
Proof. The CDF of R n | R 0 =r 0 is obtained by conditioning R on R 0 as follows:

Notation
Description r ⊕ ; r min Earth radius (6371 km); Altitude of satellites R 0 ; R n ; R Serving distance; Distance to the n th interfering satellite; Distance from the user to any satellite G 0 ; G n Channel fading gain of the serving link; Channel fading gain of the n th interfering link N; K; N I Total number of satellites; Total number of channels; Number of interfering satellites p s Transmission power from the serving satellite p i Transmission power from interfering satellites σ 2 Additive noise power α Path loss exponent T SINR threshold P c ;C Coverage probability; Average achievable rate The proof is then completed by taking the derivative of the CDF with respect to r n for obtaining the PDF in (7).
Furthermore, other auxiliary performance factors, which will be shortly used in analyzing coverage probability and average achievable data rate, are the number of interfering satellites N I and the probability P 0 of having an interferencefree situation (N I = 0) when R 0 = r 0 . They will be expressed in the following lemma and its corresponding corollary.
Lemma 4. When the serving satellite is at distance r 0 ≥ r min from the user, the number of interfering satellites, denoted by N I , is a binomial random variable with success probability Proof. The expression is directly given by the ratio of the surface area where visible interfering satellites can reside (the spherical cap formed by the intersection of the user's plane and the satellites' sphere) to the total surface area of the satellites' sphere excluding the shaded cap of Fig. 2.
It should be noted that there are N/K − 1 satellites sharing the same channel with the serving satellite and potentially causing co-channel interference. The probability of having no interference (N I = 0) then follows from the corollary below. Corollary 1. When the serving satellite is at distance r 0 ≥ r min from the user, the probability of having zero co-channel interference is given by for r 0 ≤ r max , and P 0 = 1 when r 0 > r max .
Especially, the above probability turns out to be a key factor for LEO constellation analysis since, for the derivation of the following performance metrics, we can separately consider the two complementary events of having either zero interference or non-zero interference.

III. COVERAGE PROBABILITY
In this section, we apply stochastic geometry to derive the downlink coverage probability of the LEO satellite network for a user in an arbitrary location on Earth. The performance measure of coverage probability is defined as where T represents the minimum SINR required for successful data transmission. In other words, whenever the SINR of the considered user from its nearest satellite is above the threshold level T, it is considered to be within the coverage of the satellite communication network.

A. Coverage Probability under Rayleigh Fading Channel
In this subsection, the propagation model takes into account the large-scale attenuation with path loss exponent α, as well as the small-scale fading. In particular, the serving channel is assumed to follow normalized Rayleigh fading so that the corresponding channel gain is an exponential random variable with unit mean, i.e., G 0 ∼ Exp (1), while interfering channels are considered to follow general fading (but Rayleigh fading for interfering channels is also analyzed as a special case). The Rayleigh fading model is applicable when the lineof-sight portion of signals received at the user's place is smaller than that of non-line-of-sight components. Due to lower altitudes in LEO constellations w.r.t. medium Earth or geostationary orbits (MEO/GEO), a strong LoS component is less likely available for the user if the number of satellites is limited. For instance, for constellations with a fewer number of satellites such as Iridium with only 66 satellites, the probability of having even medium elevation angle is rather low [27] that implies very low LoS probability. The probability of receiving LoS signals from interfering satellites is even lower due to their having smaller elevation angles. In addition, due to the shorter orbital period in LEO satellites, the channels vary rapidly with time and location which causes considerable Doppler variation at the user's location.
We express the coverage probability under the Rayleigh fading assumption for the serving channel as follows. In particular, we split the coverage probability into two terms, each of them corresponding to an important operational circumstance: The first term represents the case when there is no interference; and the second term corresponds to the case when there is at least one interfering satellite. Theorem 1. The probability of coverage for an arbitrarily located user under a Rayleigh-fading serving channel is with and where L I (s) is the Laplace transform of cumulative interference power I that is expressed in Lemma 5 and P 0 is given in (10) by Corollary 1.
Proof. To obtain (13), we start with the definition of coverage probability for the case of having zero interference: The upper limit for the integral is due to the fact that the satellites below the user's horizon are not visible to the user. Invoking the assumption regarding the distribution of G 0 , the expression in (13) is obtained. To prove (14), we first note that there is at least one interfering satellite, so that we have When G 0 ∼ Exp (1), we can express the first term in the integrand of (16) as Substituting (17) into (16) completes the derivation of (14).
It is worth highlighting that, in the case of having K = N orthogonal channels, the system becomes noise-limited and the coverage probability in (12) will reduce to its first term. In the following lemma, we will obtain the Laplace function of the random variable I to complete the derivation of (14).

Lemma 5.
When the serving satellite is at distance r 0 ≥ r min from the user, the Laplace transform of random variable I is where P I is given in (9) by Lemma 4 and L G n (·) is the Laplace transform of the random variable G n .
Proof. See Appendix B.
Consequently, by only specifying L G n (·) at the point sp i r −α n , i.e., assuming some specific fading model, the Laplace transform of interference for any fading distribution can be calculated using Lemma 5. For instance, when G n is exponentially distributed, L G n (sp i r −α n ) = 1 1+p i sr −α n . Thus, the Laplace function of interference when interfering channels are assumed to be Rayleigh is given by the following corollary.
Consequently, for Rayleigh-fading interfering channels and particular values of the path loss exponent, Lemma 5 will reduce to elementary functions, which will result in a simpler expression for (14). Using [28,Eq. 3.194.5], the integral in (19) can be rewritten as ∫ r n 1 + p i sr −α n dr n = − r 2 n 2 2 F 1 1, where 2 F 1 (·, ·; ·; ·) is the Gauss's hyper-geometric function.
Using the definition of the function [28,Eq. 9.100] and substituting with special arguments, it is reduced to elementary functions. If α = 2, we have [28, Eq. 9.121.5] so that (19) can be rewritten as On the other hand, if we assume α = 4, the Gauss hypergeometric function can be written as [28, Eq. 9.121.27] 2 F 1 1, by which (19)

B. Coverage Probability under Non-fading Channels
In this subsection, we derive the coverage probability for non-fading propagation environments. This propagation model is applicable (at least as an accurate approximation) when the number of satellites in a constellation is large enough, so that it is likely to have multiple satellites in LoS propagation range. Consequently, the serving satellite is likely high above the user and potential multi-path fading components are weak compared to the direct propagation path. Moreover, the LoS probability for interfering satellites increases as we decrease the number of frequency channels.
Similar to Theorem 1, the coverage probability is split into two terms. The first one is for the case of having no interference and, in the second term, we assume that there is at least one interfering satellite. As will be seen shortly, other than representing two important communication scenarios, such division is needed to make our approach tractable. For the non-fading serving channel, we denote G 0 = 1. Like in the previous section, any general fading statistics can be assumed for interfering channels, for which non-fading interference is considered as a special case towards the end of the section. Theorem 2. The coverage probability of a user with a nonfading serving channel is with and Proof. The derivation of (26) is straightforward by substituting G 0 → 1 in (3) The latter equality follows immediately from the Parseval-Plancherel property and from the fact that the Fourier transform of the square integrable function 1 [30].
The Laplace function L I (s) for general fading can be obtained using Lemma 5 given in the previous section. The following corollary presents the Laplace function when interfering signals experience non-fading propagation as well.
Corollary 3. When the serving satellite is at distance r 0 ≥ r min from the user, the Laplace transform of I for non-fading channels is where Γ(a, x) = ∫ ∞ x y a−1 e −y dy denotes the upper incomplete gamma function.
Proof. The first term of the integrand in Lemma 5 will be reduced to an exponential function. As a result, the integral can be expressed in form of the upper incomplete gamma function defined as Γ(a, x) = ∫ ∞ x y a−1 e −y dy by changing integration variable as sp i r −α n → y.

IV. AVERAGE ACHIEVABLE RATE
In this section, we focus on the average achievable data rate. The technical tools and expressions presented in Sections II and III are reused also herein. The average achievable rate (in bit/s/Hz) is defined as

A. Average Achievable Rate under Rayleigh Fading Channel
We can calculate the expression for the average rate of an arbitrary user under the assumption of a Rayleigh-fading serving channel as follows. It is worth noting that the average is taken over both the spatial BPP and fading distribution.
Theorem 3. The downlink average rate (in bits/s/Hz) of an arbitrarily located user and its serving satellite under Rayleigh fading assumption, i.e., G 0 ∼ Exp (1), is with E log 2 (1 + SNR) and Proof. See Appendix C.
The above applies to interfering channels with any fading statistics. However, analogous to the results in Section III, by considering G n ∼ Exp (1) for some specific α values, i.e., α = 2 and α = 4, and substituting (22) or (24) into (35) will lead to more simplified expressions for Theorem 3. Moreover, if we allocate K = N orthogonal channels to the system, the second term in (33) will be eliminated and the network's performance will become noise-limited.

B. Average Achievable Rate under Non-fading Channels
In this subsection, we derive the average data rate for nonfading serving channels and any fading statistics for interference. The rate expression is again split into two terms to consider zero and non-zero interference conditions separately.
Theorem 4. The downlink average rate of an arbitrary located mobile user and its serving satellite for a non-fading channel isC P 0 K E log 2 (1 + SNR) with Proof. See Appendix D.
In case of non-fading interfering channels, the Laplace transform in (38) can be calculated from Corollary 3.

V. NUMERICAL RESULTS
In this section, we verify our expressions for coverage probability and average data rate by comparing them with the results of Monte Carlo simulations. We obtain a new parameter-effective number of satellites (N eff )-to compensate for the effect of uneven distribution of satellites along different latitudes. We then illustrate the two performance metrics in terms of different network design parameters such as the number of frequency channels and the satellite altitude.
For numerical verification, we calculate the simulated coverage probability in Monte Carlo manner as follows. First, we randomly place N satellites with uniform distribution on a sphere centered at Earth's center and distanced r min from its surface. At the same time, we also model the channels between the satellites and the user by exponential random variables (to model Rayleigh fading) with unit mean value while, to model a non-fading environment, we set the channel gains to one. Then we calculate the SINR at the receiver and compare it with a pre-defined threshold value to evaluate coverage probability. We repeat this experiment for a large number of realizations to obtain the averaged performance metrics. Finally, we compare these numerical results with the analytical expressions in Theorems 1-4 to confirm our analysis.
For producing the numerical results, the transmitted power from all satellites is set to p s = p i = 10 W. The assumption corresponds to the case where all satellites are equipped with omni-directional antennas. The theory is applicable also to other antenna patterns by adjusting different power levels for the serving and interfering satellites accordingly to take into account that the serving satellite's transmission would be directed to the user while all interfering transmissions are sent to other, arbitrary and independent, directions. The noise power is assumed to be σ 2 = −98 dBm. The number of channels is considered to be K = 20 unless otherwise stated.

A. Corroboration of Theorems by Simulations
In this subsection, we validate the coverage probability and the data rate derived in Theorems 1-4 through Monte Carlo simulations. We considered static propagation or Rayleigh fading for serving channels and general fading for interference. Figure 3 demonstrates the coverage probability for different threshold values assuming α is 2 and 4. The parameters N and r min are set to 720 and 1200 km, respectively. We find that our theoretical results are perfectly in line with the simulations, which confirms the correctness of the derivations. The transition of coverage probability from one to zero occurs variously depending on path loss exponent and the type of fading. For α = 2, the coverage probability saturates when T > 20 dB which corresponds to the case when the number of interfering satellites is considerable. As a result, the SINR and, consequently, the second term in Theorem 1 approach zero. For larger path loss exponents, the effect of interference becomes insignificant and the transition from one to zero includes no saturation mode since the difference between SNR and SINR becomes less dominant. All curves tend to zero for thresholds greater than 80 dB although not shown in the figure. Figure 4 shows average achievable data rate versus the number of frequency bands for Rayleigh fading and nonfading environments. The behavior of the curves can be justified according to two contradictory effects of increasing the number of orthogonal channels on data rate. Allocating more orthogonal channels improves data rate by mitigating the interference but, at the same time, degrades it by making only a portion of the whole band available for each group of satellites. For α = 2, the former effect dominates for K < 45 as by increasing K, data rate increases accordingly, while for K > 45, the latter effects overcomes the interference elimination effect and data rate starts falling.
The effect of increasing the number of channels on the reduced available frequency band dominates the SINR enhancement for α = 4 since the largest average data rate for α = 4 corresponds to the lowest number of frequency channels. This is due to the fact that the larger path loss exponents make the interfering satellites less effective since they transmit from a farther distance. As it can be observed in the figure, for AWGN channel and α = 2, the behaviour of the data rate is slightly affected by the fading model.
The effect of frequency reuse on coverage probability can be determined from Fig. 5. The channel allocation policy is such that the channel corresponding to the nearest satellite to the user will be selected. As can be observed, the coverage probability rises with increasing the number of frequency bands due to interference mitigation, and all three curves present the same behaviour. However, the effect of the parameter K is more significant for the smaller path loss exponent as it affects the interference level more significantly. The coverage performance is superior for non-fading environment due to better quality of the serving signal. It can be interpreted from Figs. 4 and 5 that the coverage probability is an increasing function of the number of channels, i.e., the higher number of frequency channels will result in better coverage probability, while there is an optimum number of channels that maximizes the data rate. Therefore, the number of frequency channels should be compromised according to the performance demands of the intended constellation.

B. Effective Number of Satellites
In this subsection, we provide analysis of deterministic Walker satellite constellations with different inclination angles which are expected to provide global Internet broadband services to users by 720 LEO satellites. In Fig. 6, the coverage probability of the Walker constellation, with inclination angles, φ i , of 90 • (polar), 70 • and 40 • are plotted for comparison to a random constellation. The user's latitude, denoted by φ u , is assumed to be 30 • . The slight difference of coverage probability in random and Walker constellations is due to the fact that the latitudal density of satellites in deterministic constellations is not uniform, as there are more satellites around their inclination angle (i.e., the latitude limits) and less in the equatorial regions. Consequently, the coverage probability changes as we move from the equator to the latitude limits.
In order to compensate for uneven latitudal density, we introduce a new parameter, entitled as the effective number of satellites, and denote it by N eff (φ i , φ u ), which is the number of satellites for a random uniform constellation that corresponds to the same coverage probability as for a practical nonuniform constellation with inclination angle φ i and at the user latitude φ u . This parameter is approximated by mean absolute error minimization between the coverage probability given in Theorem 1 and the results from practical Walker constellation for a few threshold values on different latitudes. The approximated value for N eff is then refined to give the best matching between random and practical constellation based on the target performance metric. Although we used some of the values from simulated results to obtain N eff , the same N eff can be used for other system parameters, i.e., path loss exponent and/or SINR threshold values. Likewise, the acquired N eff can be also used for a constellation with different total number of satellites by linearly scaling it accordingly.
In Fig. 7, we eliminated the deviation between Walker and random constellations (cf. Fig. 6) by compensating the non-  uniform latitudal density. The solid markers on the curves were chosen from practical constellation in order to minimize the mean absolute error in coverage probability of a random constellation. As it can be seen, although only a few points are employed for fitting, the selected N eff also holds for other values that had no contribution on the minimization process. In other words, applying N eff in theorems provided in this paper will lead to more precise results for many system parameters and performance metrics by taking into account only a few input data from the actual network.
In Fig. 8, the effect of latitude of user on coverage probability is shown for three practical constellations. In this figure, we have two groups of data that demonstrate the coverage probability. The first group (depicted by markers) is obtained by simulation of a Walker constellation for polar,  70 • and 40 • inclined orbits. For a polar constellation, the number of satellites in the view increases monotonically from the equator to poles. As a result, the coverage probability declines as the user moves from the equator to poles due to intensifying interference. On the other hand, for the other two inclination angles, the number of satellites in the view varies non-uniformly due to the network geometry which results in coverage alternation for different latitudes. The zero coverage probability for user latitudes greater than 75 • in 40 • -inclined Walker results from the fact that there are no visible satellites for those latitudes.
The second group of data (depicted by lines) in Fig. 8 is generated using Theorem 1 with N eff (φ i , φ u ) satellites. As can be seen in the figure, coverage probability in Theorem 1 for N eff (φ i , φ u ) results in the same coverage from Walker constellation for the given φ u . Therefore, using N eff (φ i , φ u ) in Theorem 1, the effect of uneven distribution of satellites along different latitudes is compensated and the coverage probability can be obtained for any actual constellation. Most importantly, we can conclude that the approximation error in modeling deterministic constellations as random ones is actually rather minimal when the varying satellite density at different latitudes is taken into account.
The data rate of a Walker satellite constellation for different inclination angles to compare with a random constellation (viz. Theorem 3) is depicted in Fig. 9. As mentioned for Fig. 6, the difference between the random and Walker constellation is mainly due to the non-uniform density of satellites in Walker constellation. The optimal number of frequency channels that can be observed in this figure for each of the constellations increases by decreasing the inclination angle. Figure 10 depicts the effect of applying N eff to Theorem 3 in eliminating the observed deviation in Fig. 9. The selected N eff values are slightly different from those in Fig. 7 since they are refined to create the best possible matching between the data rate of the random and practical constellation. Figure 11 is a counterpart for Fig. 8 to illustrate the effect of latitude of the user on achievable data rate and N eff (φ i , φ u ). In this figure, N eff (φ i , φ u ) is the number of satellites in a random constellation that corresponds to the same data rate as a Walker constellation for the given φ u . Due to different values of N eff (φ i , φ u ) in Figs. 8 and 11, we can realize that this value is not only a function of geometrical characteristics of the constellation but also depends slightly on the performance metric. However, the compensation of latitudal density is still applicable to varying system parameters that are different from those used for fitting theoretical expressions to simulations.

C. The Effect of Satellite Altitude
Coverage probability for different altitudes is then depicted in Fig. 12. The figure is plotted using N eff (φ i , φ u ) for Theorem 1. In this figure, we assumed that the total number of satellites varies with altitude so that N eff remains constant irrespective of altitude. This balances between a practical constellation design in which the number of satellites (when aiming at covering economically the entire globe) typically decreases with increasing altitude and the direct relationship between altitude and the effective number of satellites. It can be observed from Fig. 12 that, for practical LEO altitudes, the coverage probability declines when altitude increases. Within Earth's atmosphere, there is an impractical optimal altitude, which moves even lower as the number of satellites increases.
Data rate for different altitudes is plotted in Fig. 13. The figure is created using Theorem 3 with N eff (φ i , φ u ) that corresponds to a specific latitude φ u for Walker constellations (polar and inclined). The same assumption as for Fig. 12 is also made regarding the invariance of N eff w.r.t. altitude. This serves to show that, having N eff (φ i , φ u ) at hand, we can use the expressions in Theorems 3 or 4 to analyze the data rate for any given satellite constellation and the approximation error in modeling a deterministic constellation as a binomial point process is insignificant. The same as for Fig. 12, it can be seen that the peak value of data rate does not exist within the LEO practical altitude range outside atmospheric drag.

VI. CONCLUSIONS
In this paper, we presented a tractable approach for downlink coverage and rate analysis of low Earth orbit satellite networks. Using the concept of stochastic geometry, we modeled the satellite network as a binomial point process.
We then applied this model to obtain exact expressions for coverage probability and data rate of an arbitrary user in terms of network parameters and the Laplace transform of interference power. The performance metrics of the random and real constellation match with each other almost perfectly; there is only a slight deviation between them which can be compensated by taking into account the effect of uneven satellite distribution along different latitudes. A frequency reuse approach was also applied to make the network spectral efficiency suitable for commercial operation. Its effect on data rate is more ambivalent than on coverage probability. The proposed framework in this paper paves the way for accurate analysis and design of the future dense satellite networks. The scalable results presented here are not only applicable to a single constellation but also to massive satellite networks that merge several different constellations.

A. Proof of Lemma 1
From basic geometry, the CDF of the surface area of the shaded spherical cap A cap , formed by any satellite at distance R from the user, in Fig. 2 To find the distribution of R, we need to find a relationship between A cap and R, when know that for which a and b are given in Fig. 2. Combining (39) and (40), we have Using the formula for the surface area of a spherical cap, i.e., A cap = 2π(r ⊕ + r min ) 2 (1 − cos θ), we obtain R 2 = r 2 min − 2r min (r ⊕ + r min ) The CDF can then be deduced as follows: Finally, the corresponding PDF can be derived by differentiating (43) with respect to r.

B. Proof of Lemma 5
In this Appendix, we provide the proof to the expression of Laplace function for general fading channels. Using the definition of the Laplace transform yields where (a) follows from the i.i.d. distribution of G n and its further independence from N I and R n , (b) is obtained using the interfering distance distribution from (7), (c) is the averaging over the binomial random variable N I with the success probability P I , which is given in Lemma 4,and (d) is the substitution from the definition of Laplace function.

C. Proof of Theorem 3
We only provide the proof for (35) herein, while the proof for (34) follows the same approach as well. In particular, where c 0 = N 2 ln(2)r ⊕ (r ⊕ +r min ) and (a) follows from the fact that for a positive random variable X, E [X] = ∫ t>0 P (X > t) dt.

D. Proof of Theorem 4
Assuming G n = 1, we have where c 0 = N 2 ln(2)r ⊕ (r ⊕ +r min ) . Following the same approach as in (28)-(30) and using Parseval-Plancherel formula will result in (38). The proof for (37) can be obtained easily using the same principles as the above and only by substituting I = 0.