Downlink Performance of Cell-Free Massive MIMO for LEO Satellite Mega-Constellation

Many papers have studied the performance of LEO satellite networks. The coverage probability for downlink transmissions was investigated in [14], where satellite locations were distributed as a binomial point process (BPP). However, crucial instruments such as the distance distributions of the BPP are more challenging to analyze compared to those of a Poisson point process (PPP) in unbounded space and therefore, PPP-based satellite modeling was adopted in [15]. Many other works also used PPP to model the constellation. For example, ultra-dense LEO satellite constellations were studied in [16] to find the minimized number of satellites while satisfying the backhaul requirement of each UT. The uplink interference and performance of mega satellite constellations were examined in [17, 18] under the impacts of intra-constellation interference. The above work, although similar to our studied network, assumed that each UT was connected to a single satellite, thus overlooking the possibility of satellite cooperation.

The CF-mMIMO system, where each UT is coherently served in the same time-frequency resource block by a large number of geographically-located access points (APs) [19], has gained extensive research attention recently. Although there has been extensive research on CF-mMIMO in TNs such as [20, 21], relatively few studies have explored its application in LEO SatCom. An LEO satellite was initially considered an “add-on” to improve terrestrial CF networks in [22, 23]. However, they fail to introduce the CF-mMIMO architecture to large-scale LEO satellite networks. The authors in [24, 25] integrated CF-mMIMO into LEO satellite networks and developed an optimization framework to improve spectral efficiency. The benefits of CF-mMIMO in LEO satellite networks were further evaluated in [26] for broadband connectivity with multi-antenna satellites and handheld devices. In addition, the CF-mMIMO SatCom network was also studied from physical layer orthogonal time-frequency-space (OTFS) [27], uplink transmissions with CSI uncertainties [28, 29, 30], and dynamic clustering of satellite access points (SAPs) [31]. However, the above works focused on algorithm designs and assumed that the served UTs were confined in a small region on the earth, neglecting other coverage areas where inter-user interference could arise, as well as downlink interference from non-cooperative satellites. Our previous work [1] examined a basic model for coordinated multi-satellite joint transmissions. However, its reliance on Rayleigh fading and the omission of beamforming and channel estimation strategies made it misaligned with CF-mMIMO and impractical for channels with line-of-sight (LoS) components. Although satellite cooperation was considered in [32, 33], inter-satellite interference persisted, and the lack of beamforming and channel estimation made the structure incompatible. In addition, its tractability decreased as the number of satellites increased.

To thoroughly investigate the downlink performance of CF-mMIMO LEO SatCom networks from a system-level perspective, we incorporate essential factors of satellite networks for analysis using tools from stochastic geometry. The PPP model is adopted to model LEO satellite mega-constellation, providing analytical results with broad applicability and high tractability. Notably, the main contributions of this paper are summarized as follows:

• •

  CF-mMIMO LEO SatCom Network Design: We establish a CF-mMIMO system in the LEO satellite network using stochastic geometry. The locations of SAPs and UTs are modeled on two concentric sphere surfaces, one with the SAPs’ orbital radius and the other with the earth’s radius, respectively, based on two independent homogeneous PPPs. Simulation results demonstrate that the PPP constellation model can fit perfectly with practical constellations. Considering the shape of the earth and the service range, each UT is assigned a group of SAPs for service delivery according to its dome angle.

• •

  Modeling and Analysis: We introduce the related distance distribution for the service links of a typical UT and characterize the average number of UTs served by each SAP. Then, we provide expressions for the distribution of desired signal strength (DSS), the average inter-satellite interference, and the average multi-user interference. The coverage probability of a typical user is then presented accordingly, with both large-scale path-loss and small-scale Nakagami-$m$ fading. Furthermore, given the relatively minor impact of multi-path fading components against the long satellite-to-terrestrial distances, we derive the coverage probability for non-fading propagation environments.

• •

  Network Design Insights: We quantitatively analyze the performance of the CF-mMIMO LEO SatCom network concerning key network parameters, including the fading parameter, the service range of CF-mMIMO LEO SatCom network, the total number of SAPs, the number of UTs, and SAPs’ orbital altitudes. We observe that as the dome angle of the typical UT increases, the gain in desired signal power outweighs the increase in received multi-user interference, leading to improved coverage probability, albeit with increased signaling overhead. In addition, an increase in both the orbital altitude and the number of SAPs contributes positively to enhanced coverage performance. While system capacity increases as the number of UTs grows within a certain range, the per-user capacity continues to drop, exhibiting a tradeoff between maximizing system capacity and ensuring minimum throughput for each UT.

The rest of this paper is organized as follows. Section II describes the system model for the studied CF-mMIMO SatCom network, including spatial distributions, channel models, and downlink transmissions. In Section III, we present the statistical properties of the network and derive analytical expressions for the coverage probabilities. Section IV explores the studied system under non-fading channel conditions. Simulations and numerical results are provided in Section V. Finally, Section VI concludes the paper.

Notations: Bold lowercase letters denote column vectors. The superscripts $\left(\cdot\right)^{T}$, $\left(\cdot\right)^{*}$, and $\left(\cdot\right)^{H}$ are used to represent conjugate, transpose, and conjugate-transpose, respectively. The real-part operator, expectation operator, the statistical probability, and the Euclidean norm are represented as $\mathrm{Re}[\cdot]$, $\mathbb{E}\left\{\cdot\right\}$, $\mathbb{P}\left\{\cdot\right\}$ and $\|\cdot\|$, respectively. $y\sim\mathcal{CN}(\mu,\sigma^{2})$ denotes a complex circularly symmetric Gaussian random variable with mean $\mu$ and variance $\sigma^{2}$, while $\bm{y}\sim\mathcal{CN}(\mu,\mathbf{R})$ denotes that with mean $\mu$ and correlation matrix $\mathbf{R}$.

Denote the radius of the earth as $R_{\text{E}}$, and that of a sphere $S$ where the SAPs are located as $R_{\text{S}}$. An illustration of PPP-distributed satellites and PPP-distributed UTs are shown in Fig. 1(a). The distributions of UTs and SAPs are independent and each follows a homogeneous spherical Poisson point process (SPPPs). In a practical constellation such as Starlink, the distribution of satellites is correlated. For comparison, we simulate the Starlink constellations with random and fixed initial states, respectively. Note that in a random initial state shown in Fig. 1(b), the initial position and orbital parameters of each satellite are randomly distributed within a certain range, which can simulate the randomness of constellations at the time of deployment. This is very similar to the constellation in a practical constellation. In a fixed initial state shown in Fig. 1(c), the initial position and orbital parameters of all satellites are set in advance and strictly fixed to represent the ideal precise deployment situation. The satellites follow strict geometric symmetries in their orbital planes and phases in order to minimize relative deviations between satellites. Intuitively, the PPP model can perfectly simulate a practical constellation by setting a random initial state, which will be verified in Section V.

On this basis, we assume that there are $L$ SAPs in the constellation at the same orbital altitude $H_{\text{S}}=R_{\text{S}}-R_{\text{E}}$, and there are a total of $K$ UTs on the earth surface that need to be served by the CF-mMIMO SatCom network. Taking the mean of the Poisson random variable, the PPP density of UTs $\lambda_{\text{U}}$ and that of SAPs $\lambda_{\text{S}}$ is approximated by $\lambda_{\text{U}}=\frac{L}{4\pi{R_{\text{E}}}^{2}}$ and $\lambda_{\text{S}}=\frac{K}{4\pi{R_{\text{S}}}^{2}}$, respectively.

For the satellite-terrestrial signal transmissions, both large-scale path-loss attenuation and small-scale fading are taken into account. Specifically, a valid assumption is that the fading channel between every two links is assumed to be uncorrelated since SAPs are geographically distributed on the sphere $S$. We denote the channel coefficient between the $l$-th SAP and the $k$-th UT as

$$g_{lk}={\beta_{lk}}^{1/2}h_{lk},$$

(1)

where $\beta_{lk}$ denotes the large-scale fading and $h_{lk}$ is the small-scale fading³³3Similar to [14, 15, 16, 17, 18], the perfect CSI is assumed. The derived expressions for the coverage probability provide upper-bound performance for the system. However, in this paper, more pilots will bring better CSI. We will evaluate how the number of pilots influences the DSS in Section V.. For large-scale fading, $\beta_{lk}$ can be represented using a distance-dependent model, i.e. $\beta_{lk}=\beta_{0}{d_{lk}}^{-\alpha}$, where $\beta_{0}$ is the path-loss at a reference distance, and $\alpha$ is the path-loss exponent. Then, the effective antenna gain is denoted as

$$G=G_{\text{t}}G_{\text{r}}\left(\frac{c}{4\pi f_{\text{c}}}\right)^{2},$$

(2)

where $G_{\text{t}}$ and $G_{\text{r}}$ are the transmit antenna gain of the SAP and the receive antenna gain of the UT, $c$ is the speed of light, and $f_{\text{c}}$ is the carrier frequency. The small-scale fading is modeled via the Nakagami-$m$ distribution, whose probability density function (PDF) is given by

$$f_{|h_{lk}|}(x)=\frac{2m^{m}}{\Omega^{m}\Gamma(m)}x^{2m-1}e^{-mx^{2}},$$

(3)

with $x\in[0,\infty]$, $m\in[0.5,\infty]$ and $\Gamma(m)=(m-1)!$ for integer $m>0$. The Nakagami-$m$ model is selected because of its adaptation to various small-scale fading conditions, e.g. Rayleigh channel when $m=1$ and Rician-$\mathcal{K}$ channel when $m=\frac{(\mathcal{K}+1)^{2}}{2\mathcal{K}+1}$ [15]. We also have $\mathbb{E}\left\{{{\left|{{h}_{lk}}\right|}^{2}}\right\}=\Omega$. Note that Shadowed-Rician (SR) fading can also be approximated by Nakagami-$m$ fading [18]. Moreover, changing this parameter $m$ makes it suitable for describing various channels, as it allows one to model the signal propagation conditions from severe to moderate.

In a CF-mMIMO system, the propagation channels are considered to be piecewise constant during each coherence time period and frequency coherence interval. Therefore, training must be conducted within each time-frequency coherence block. With channel reciprocity, channel estimation during uplink training can be used for uplink and downlink transmissions [19]. Herein, consider a coherence block of length $\tau_{\text{c}}$, where the durations of $\tau_{\text{p}}$ and $\tau_{\text{d}}=\tau_{\text{c}}-\tau_{\text{p}}$ are dedicated to uplink training and downlink transmission, respectively. In the uplink training stage, let $\sqrt{\rho_{\text{p}}}\bm{\varphi}_{k}\in\mathbb{C}^{\tau_{\text{p}}\times 1}$ be the pilot sequence transmitted by the $k$-th UT, where $||\bm{\varphi}_{k}||^{2}=\tau_{\text{p}}$, and $\rho_{\text{p}}$ represents the transmit power of the pilot symbols.

Denote the set of all served UTs of the $l$-th SAP as $\Phi_{l}^{\text{U}}$, the received signal at the $l$-th SAP can be given by

$$\bm{y}_{\text{p},l}=\sqrt{\tau_{\text{p}}\rho_{\text{p}}G}\sum_{k\in\Phi_{l}^{\text{U}}}{\beta_{lk}}^{1/2}h_{lk}\bm{\varphi}_{k}^{H}+\bm{w}_{\text{p},l},$$

(4)

where $\bm{w}_{\text{p},l}\in\mathbb{C}^{\tau_{p}\times 1}$ is the receiver additive white Gaussian noise (AWGN) vector whose elements are i.i.d. as $\mathcal{CN}(0,\sigma^{2})$. The post-processing signal is obtained by projecting $\bm{y}_{\text{p},l}$ onto $\bm{\varphi}_{k}$, which is represented as

	$\displaystyle\bar{y}_{\text{p},l,k}$	$\displaystyle=\bm{y}_{\text{p},l}\bm{\varphi}_{k}$		(5)
		$\displaystyle=\sqrt{\tau_{\text{p}}\rho_{\text{p}}G}\sum_{k^{\prime}\in\Phi_{l}^{\text{U}}}\beta_{lk^{\prime}}^{1/2}h_{lk^{\prime}}\bm{\varphi}_{k^{\prime}}^{H}\bm{\varphi}_{k}+\bm{w}_{\text{p},l}\bm{\varphi}_{k}.$

The estimated channel $\hat{g}_{lk}$ using linear minimum mean-square-error (LMMSE) algorithm can then be written as

	$\displaystyle\hat{g}_{lk}$	$\displaystyle=\frac{\mathbb{E}\{\bar{y}_{\text{p},l,k}{g}_{lk}\}}{\mathbb{E}\{|\bar{y}_{\text{p},l,k}|^{2}\}}\bar{y}_{\text{p},l,k}$		(6)
		$\displaystyle=\begin{aligned} &\frac{\sqrt{\tau_{\text{p}}\rho_{\text{p}}G}\beta_{lk}\mathbb{E}\{h_{lk}h_{lk}^{H}\}}{\tau_{\text{p}}\rho_{\text{p}}G\sum_{{k^{\prime}}\in\Phi_{l}^{\text{U}}}\beta_{lk^{\prime}}|h_{lk^{\prime}}|^{2}|\bm{\varphi}_{k^{\prime}}^{H}\bm{\varphi}_{k}|^{2}+\sigma^{2}}\\ &\times\left(\sqrt{\tau_{\text{p}}\rho_{\text{p}}G}\sum_{k^{\prime}\in\Phi_{l}^{\text{U}}}\beta_{lk^{\prime}}^{1/2}h_{lk^{\prime}}\bm{\varphi}_{k^{\prime}}^{H}\bm{\varphi}_{k}+\bm{w}_{\text{p},l}\bm{\varphi}_{k}\right)\end{aligned}$
		$\displaystyle=\begin{aligned} &\beta_{lk}\mathbb{E}\{h_{lk}h_{lk}^{H}\}\\ &\times\begin{aligned} &\left(\frac{\sum_{k^{\prime}\in\Phi_{l}^{\text{U}}}\beta_{lk^{\prime}}^{1/2}h_{lk^{\prime}}\bm{\varphi}_{k^{\prime}}^{H}\bm{\varphi}_{k}}{\sum_{{k^{\prime}}\in\Phi_{l}^{\text{U}}}\beta_{lk^{\prime}}|h_{lk^{\prime}}|^{2}|\bm{\varphi}_{k^{\prime}}^{H}\bm{\varphi}_{k}|^{2}+\frac{\sigma^{2}}{\tau_{\text{p}}\rho_{\text{p}}G}}\right.\\ &\left.+\frac{\sqrt{\tau_{\text{p}}\rho_{\text{p}}G}\bm{w}_{\text{p},l}\bm{\varphi}_{k}}{\tau_{\text{p}}\rho_{\text{p}}G\sum_{{k^{\prime}}\in\Phi_{l}^{\text{U}}}\beta_{lk^{\prime}}|h_{lk^{\prime}}|^{2}|\bm{\varphi}_{k^{\prime}}^{H}\bm{\varphi}_{k}|^{2}+\sigma^{2}}\right).\end{aligned}\end{aligned}$

During the downlink transmission stage, SAPs treat the channel estimates as the actual channels to perform conjugate beamforming, allowing them to simultaneously send data symbols to UTs within their respective service areas. Focusing on the studied network, we employ equal power allocation at the transmitting SAP for each served UT by using conjugate beamforming.

Specifically, denoting the total transmit power of each SAP as $\rho_{\text{d}}$, and the number of UTs within the $l$-th SAP as $\left|\Phi_{l}^{\text{U}}\right|$, we have $\rho_{lk}^{\text{d}}=\frac{\rho_{\text{d}}}{\left|\Phi_{l}^{\text{U}}\right|}$ represent the transmit power from $l$-th SAP to $k$-th UT. Then, the transmitted symbols from the $l$-th SAP can be expressed as

	$\displaystyle x_{l}$	$\displaystyle=\sqrt{G}\sum_{k\in\Phi_{l}^{\text{U}}}\sqrt{\rho_{lk}^{\text{d}}}\frac{\hat{g}_{lk}^{*}}{|g_{lk}|}q_{k}$		(7)
		$\displaystyle=\sqrt{G}\sum_{k\in\Phi_{l}^{\text{U}}}\sqrt{\rho_{lk}^{\text{d}}}\frac{\hat{h}_{lk}^{*}}{|h_{lk}|}q_{k},$

where $q_{k}$ is the dummy or data symbol intended for the $k$-th UT with $\mathbb{E}\{|q_{k}|^{2}\}=1$, and $\hat{h}_{lk}$ is the estimated channel of $h_{lk}$ using LMMSE.

The signals received by the $k$-th UT can be written as

	$\displaystyle r_{k}^{\text{d}}$	$\displaystyle=\sum\limits_{l\in\Phi_{k}^{\text{Ser}}}{{{\beta}_{lk}}^{1/2}{{h}_{lk}}{{x}_{l}}}+\sum\limits_{l\in\Phi_{k}^{\text{Int}}}{\sqrt{{{\rho}_{\text{d}}}G}{{\beta}_{lk}}^{1/2}{{h}_{lk}}{{q}_{k}}}+w_{\text{d},k}$		(8)
		$\displaystyle=\begin{aligned} &\underbrace{\sum_{l\in\Phi_{k}^{\text{Ser}}}\sqrt{\rho_{lk}^{\text{d}}G}{\beta_{lk}}^{1/2}h_{lk}\frac{\hat{h}_{lk}^{*}}{|h_{lk}|}q_{k}}_{\mathrm{DS}_{k}}\\ &+\underbrace{\sum_{l\in\Phi_{k}^{\text{Ser}}}\sum_{\begin{subarray}{c}k^{\prime}\in\Phi_{l}^{\text{U}}\\ k^{\prime}\neq k\end{subarray}}\sqrt{\rho_{lk}^{\text{d}}G}{\beta_{lk}}^{1/2}h_{lk}\frac{\hat{h}_{lk^{\prime}}^{*}}{|h_{lk^{\prime}}|}q_{k^{\prime}}}_{\mathrm{MUI}_{k}}\\ &+\underbrace{\sum\limits_{l\in\Phi_{k}^{\text{Int}}}{\sqrt{{{\rho}_{\text{d}}}G}{{\beta}_{lk}}^{\frac{1}{2}}{{h}_{lk}}{{q}_{k}}}}_{\mathrm{ISI}_{k}}+w_{\text{d},k},\end{aligned}$

where $\Phi_{k}^{\text{Ser}}$ and $\Phi_{k}^{\text{Int}}$ are the set of service SAPs and interference SAPs, $\Phi_{l}^{\text{U}}$ is the set of UTs within $l$-th SAP’s service area, and $w_{\text{d},k}\sim\mathcal{CN}(0,\sigma^{2})$ is the downlink AWGN. Moreover, $\mathrm{DS}_{k}$, $\mathrm{MUI}_{k}$, and $\mathrm{ISI}_{k}$ represent the desired signal, the multi-user interference and the inter-satellite interference, received at the $k$-th UT, respectively.

Related parameters in the stochastic geometry sketch of UTs and SAPs are labeled in Fig. 2. According to the 3GPP Technical Standard [34], each ground UT has a maximum service range for satellites. Thus, the typical UT has a dome angle $\eta$ from its position and corresponds to the shortest vertical distance $H_{v}$ from the UT. We denote the area marked in light yellow as a “service range” corresponding to $\eta$. If an SAP is within the service range, that is, if its vertical distance from the UT is longer than $H_{v}$, it is considered a service SAP. A cluster of SAPs within the service range is coordinated by the CS to serve the UT coherently.

Let the distance from the typical UT to any one of the service SAPs be $D$. According to [1, Lemma 1], the cumulative distribution function (CDF) of $D$ is given in (9) at the bottom of the page,

where $\eta\in[0,\frac{\pi}{2}]$, and the corresponding PDF is given by

$\displaystyle f_{D}(d)=\frac{d}{R_{\text{E}}(R_{\text{S}}-R_{\text{E}}\sin^{2}\eta-\sqrt{{R_{\text{S}}}^{2}-{R_{\text{E}}}^{2}\sin^{2}\eta}\cos\eta)},$

(10)

for $r_{\text{S,min}}\leq d\leq r_{\text{S,max}}$, while $f_{D}(d)=0$, otherwise. Note that $r_{\text{S,min}}={R_{\text{S}}}-{R_{\text{E}}}$ and $r_{\text{S,max}}=\sqrt{{R_{\text{S}}}^{2}-{R_{\text{E}}}^{2}{{\sin}^{2}}\eta}-{R_{\text{E}}}\cos\eta$ are the minimum and maximum distances from the typical UT to a service SAP, respectively, and ${{r}_{\max}}=\sqrt{R_{\text{S}}^{2}-R_{\text{E}}^{2}}$ is the maximum distance from the typical UT to any SAP [15].

Next, the shortest vertical distance from service SAPs to the typical UT is given by [1, Lemma 2]

$$H_{v}=\begin{cases}\cos^{2}\eta\left(\sqrt{{R_{\text{E}}}^{2}+\frac{{R_{\text{S}}}^{2}-{R_{\text{E}}}^{2}}{\cos^{2}\eta}}-R_{\text{E}}\right),&0\leq\eta<\frac{\pi}{2},\\ 0,&\eta=\frac{\pi}{2}.\end{cases}$$

(11)

Then, the maximum distance from a service SAP to the typical UT is $\frac{H_{v}}{\cos\eta}$ for $0\leq\eta<\frac{\pi}{2}$, and $\sqrt{{R_{\text{S}}}^{2}-{R_{\text{E}}}^{2}}$ for $\eta=\frac{\pi}{2}$.

Finally, due to the change of $\eta$, each SAP has a service area on the earth’s surface. The average number of UTs within this service area is written as [1, Lemma 3]

$$\left|\Phi_{l}^{\text{U}}\right|_{\text{avg}}=2\pi R_{\text{E}}\lambda_{\text{U}}\begin{aligned} &\left(R_{\text{E}}-\frac{{R_{\text{E}}}^{2}}{R_{\text{S}}}\sin^{2}\eta\right.\\ &\left.-\frac{R_{\text{E}}\sqrt{{R_{\text{S}}}^{2}-{R_{\text{E}}}^{2}\sin^{2}\eta}\cos\eta}{R_{\text{S}}}\right).\end{aligned}$$

(12)

The coverage probability is determined by the signal-to-interference-plus-noise ratio (SINR) of the message received by a typical UT and a threshold $\gamma_{\text{th}}$. When the SINR is above the threshold $\gamma_{\text{th}}$, this UT can successfully decode the data. Following the expression for the signals received by the $k$-th UT in (8), the corresponding SINR is

$$\mathrm{SINR}_{k}=\frac{{{\left|\mathrm{DS}_{k}\right|}^{2}}}{\sum\limits_{k^{\prime}\neq k}{\left|\mathrm{UI}_{kk^{\prime}}\right|}^{2}+{\left|\mathrm{SI}_{k}\right|}^{2}+{{\sigma}^{2}}},$$

(13)

where $\mathrm{DS}_{k}=\sum\limits_{l\in\Phi_{k}^{\text{Ser}}}{\sqrt{\rho_{lk}^{\text{d}}G}{{\beta}_{lk}}^{\frac{1}{2}}{{h}_{lk}}\frac{\hat{h}_{lk}^{*}}{\left|{{\hat{h}}_{lk}}\right|}}q_{k}$, $\mathrm{UI}_{kk^{\prime}}=\sum\limits_{l\in\Phi_{k}^{\text{Ser}}}{\sqrt{\rho_{lk^{\prime}}^{\text{d}}G}{{\beta}_{lk}}^{\frac{1}{2}}{{h}_{lk}}\frac{\hat{h}_{lk^{\prime}}^{*}}{\left|{{\hat{h}}_{lk^{\prime}}}\right|}}q_{k^{\prime}}$, $\mathrm{ISI}_{k}=\sum\limits_{l\in\Phi_{k}^{\text{Int}}}{\sqrt{{{\rho}_{d}}G}{{\beta}_{lk}}^{\frac{1}{2}}{{h}_{lk}}}q_{k}$.