A hybrid beamforming design for massive MIMO LEO satellite communications

5G and future cellular networks intend to incorporate low earth orbit (LEO) satellite communication systems (SatCom) to solve the coverage and availability problems that cannot be addressed by satellite-based or ground-based infrastructure alone. This integration of terrestrial and non terrestrial networks poses many technical challenges which need to be identified and addressed. To this aim, we design and simulate the downlink of a LEO SatCom compatible with 5G NR, with a special focus on the design of the beamforming codebook at the satellite side. The performance of this approach is evaluated for the link between a LEO satellite and a mobile terminal in the Ku band, assuming a realistic channel model and commercial antenna array designs, both at the satellite and the terminal. Simulation results provide insights on open research challenges related to analog codebook design and hybrid beamforming strategies, requirements of the antenna terminals to provide a given SNR, or required beam reconfiguration capabilities among others.


I. INTRODUCTION
Integrated satellite-terrestrial cellular networks are currently being pursued for long awaited commercial applications. The main expected benefit is an uninterrupted coverage, even in un-served/underserved areas. LEO constellations seem to be the most promising platforms for satellite-based non-terrestrial networks (NTNs), due to their relatively shorter propagation delay [1]. Although this propagation delay between ground terminals and LEO satellites (3 ms at 1000 km above the Earth surface) is non-negligible, LEO constellations can support continuity and ubiquity of the radio services, and also latency critical applications with requirements within tens of ms.
Beam configuration and reconfiguration in a link between a LEO satellite and a ground mobile terminal is challenging due to the long round trip delay, fast movement of the satellite and limited on board processing capabilities. Current LEO constellations [2] make use of fixed analog beams that illuminate a given area of the Earth's surface, without any capability to steer narrow beams in the user directions. Digital precoding stages are designed independently of the analog beam with the purpose of reducing inter-beam interference among the beams illuminating adjacent regions of interest (ROI) [3]. Moreover, these solutions assume that a single satellite is illuminating a large specific ROI, and handover procedures are put in place to enable satellite switching at the user side before the satellite movement causes loss of coverage.
Designing the LEO satellite footprint and the specific beam codebook at the satellite side in an integrated satellite-terrestrial network is challenging. The potential beampatterns need to be designed to guarantee coverage of the Earth's surface, limit the inter-beam interference, maximize system throughput, and, at the same time, exhibit some degree of compatibility with conventional codebooks defined in cellular standards. Commercial efforts on deploying LEO constellations, not necessarily for integration with the terrestrial network, are described in [2], [4]. Sizes of the beam footprints for the different deployments are provided, but the details of the beam codebook are not available. The massive MIMO LEO satellite communication system designed in [5] considers the set of all feasible beams generated by a uniform planar array at the satellite rather than a codebook, without trying to reduce the inter-beam interference or the amount of on-board processing.
In this paper, we propose and evaluate a massive MIMO LEO satellite communication system operating in the Ku band, based on a hybrid beamforming architecture. In particular, we design first the footprint associated to a given LEO satellite. Then, we design the beam codebook for the hybrid beamforming stage such that the derived footprint is covered by the beams in the codebook. We build our design using a 2-D DFT-based grid of beams as in Type I and Type II 5G New Radio CSI codebooks [6], [7], [8], [9], [10], [11]., considering an oversampling April [2], and also compliant with the 3GPP proposal for 5G broadband satellite services [12]. We denote as B DL the bandwidth allocated to the downlink. The system model we propose is general, but in the performance evaluation section we will focus on the Communications on the move (COOM) 5G use case. A given LEO satellite in the constellation can cover an elliptical region of interest (ROI), with semi-radius R x and R y , using a set of potential beams, as illustrated in Fig. 1. Thus, the system supports the direct transmission to a maximum number of simultaneous mobile user terminals (UTs) on the ground, denoted as N u , by using N b analog beams, with N u < N b . Time or frequency division multiplexing mechanisms will enable beam sharing when more than one user per beam is simultaneously served. The beam footprints are not necessarily identical, to account for the different elevation angles across beams. We assume that the beams generated at the LEO satellite belong to a fixed codebook, i.e. the beams move with the satellite, as considered in Scenario C1 and D2 in the 3GPP technical report that specifies solutions for 5G NR to support NTNs [13].
We also assume that the UTs know the position and trajectory of the satellites in the constellation at all times, as proposed in [14]. This information is periodically updated and broadcasted by the LEO satellite operations center to the cellular network. Knowledge of the planned trajectory allows the UTs to obtain an estimate of the satellite position even at locations where the terrestrial network is not available. It is out of the scope of this paper to consider the impact of the position or trajectory errors in the system model.

A. MIMO architecture and antenna models
We assume that the satellite is equipped with a large uniform planar array (UPA) with dimension N sat x × N sat y and N sat = N sat x N sat y denoting the total number of antennas. We consider x to be the axis in the direction of the satellite's movement and y to be the axis in the direction orthogonal to the movement. The possible beamforming architectures for the LEO satellite include: fully digital, reconfigurable analog, analog based on beam switching and hybrid [12]. Our choice is a partially connected hybrid MIMO architecture with a fixed codebook-based analog stage. The number of beamforming ports in each dimension is defined as N RF x and N RF y . The large array is organized into sub-arrays of size N sub x ×N sub y , with N sub Fig. 2. Each radio frequency (RF) chain drives one of the subarrays, generating an independent spot beam so the whole coverage area is illuminated. This approach allows the reduction of the number of beamforming ports and on board processing requirements, since the analog beams do not have to be constantly reconfigured.
Assuming a hybrid architecture, the precoding matrix F ∈ C Nsat×Nu can be written as F = F RF F BB , where F BB ∈ C N b ×Nu is the digital precoder and F RF ∈ C N sat ×N b is the analog precoder. . We envision a digital precoder at the satellite to manage inter-beam interference in the context of a full frequency reuse system. Finally, it is also assumed that all the power amplifiers at the satellite UPA are identical, and their impact can be modeled as a scaling factor √ P [12].
For the antenna arrays at the mobile terminal, we consider three possible solutions that enable beamforming in the satellite direction. The first one is a uniform planar array (UPA) of moderate size, with dimensions N UT x × N UT y , having a total of N UT = N UT x N UT y antenna elements. Recent advances on antenna technology enable the commercial applications of these planar arrays for SatCom terminals [15]. The second considered solution is a flat panel antenna based on leaky wave technology [16]. The final option is a commercial liquid crystal based metasurface antenna developed by Kymeta [17].
We define w ∈ C N UT as the vector containing the UT antenna weights over the set of feasible configurations. We consider first the UPA with antenna elements relative positions defined through k u ∈ C Nu×3 , being [k u ] n,: the n-th antenna element's relative position. In this case, the UT steering vector a u : ).
Note that there is a gain threshold 10 log 10 ) that cannot be exceeded by increasing the number of antenna elements in the leaky-wave antenna.
Finally, the metasurface antenna designed by Kymeta has a peak gain of 33 dB and a 1.2 cosine roll-off factor. This translates into a gain of 33 − 1.2 · 10 log 10 (sin(θ el )) dB towards a satellite with an elevation angle θ el relative to the user. The analytical expression of the array factor for this antenna is not available, so only system metrics that depend on the peak gain can be computed when the terminal is operating with this antenna.
Antenna sub-arrays

B. Channel model
We consider the channel between a single LEO satellite located at position x sat ∈ R 3 and a UT located at position x u ∈ R 3 on the Earth's surface, both measured from Earth's center. The satellite is orbiting Earth at a height h sat , moving with a linear speed of where R earth is the Earth's radius, m earth is the Earth's mass, and G is the gravitational constant.
This is equivalent to an angular speed of w sat = vsat r earth +hsat . The closest point to the satellite on the Earth surface, from now on the satellite's shadow, moves at a linear speed of v shadow = w sat r earth and circles the Earth with a period of 2π wsat . 1) Link budget model: The link power loss includes the free space path loss, denoted as LP fs , and the atmospheric absorption loss LP at . We neglect the impact of the atmospheric fading on the link budget. This way, the signal-to-noise ratio (SNR) can be written in terms of the the received signal strength (RSS) and the noise power σ 2 as [13] where P TX is the transmit power, G TX is the transmit antenna gain, LP cable is the cable loss between the antenna and the transmitter, G RX is the receiver antenna gain, T is the noise The free space path loss in dB is given by with x sat − x u being the traveled distance, f the frequency and c the speed of light.
Defining A(f, T (h), P (h), ρ(h)) as the loss per meter for temperature T (h), pressure P (h) and humidity ρ(h) at height h, the atmospheric path loss due to absorption is The integral is performed over the efficient limits of the atmosphere, that is between the user's height h min (h min = 0 for sea level) and the atmosphere thickness h max = 81km. The atmospheric temperature and pressure as a function of the height can be extracted using the 1976 Standard Atmospheric model depicted in Fig. 3, while the corresponding loss per meter can be computed using tabulated functions implemented in available software packages like gaspl in MATLAB.
The humidity density vertical profile ρ(h) can significantly vary depending on meteorological conditions, and its physical properties (condensation and freezing) can barely exist past its freezing point (< −0°C) without the formation of clouds. Therefore, we will assume that there are no clouds and that the water vapor is 7.5g/m 3 at sea level h = 0 and proportional to the air density in the first 2.3 km of the atmosphere, where the temperature is still over the 0 deg.

2) Satellite angular speed and Doppler effects:
The satellite movement at a speed v sat ∈ R 3 introduces a relative angular speed w rel with respect to the user and, for a signal at frequency f , a Doppler effect ∆f sat with expressions Note that the angular speed is computed in absolute terms and it is not directly related to a parametrization of choice. like for example polar coordinates. Therefore, it does not translate to variational speed of azimuth and elevation angles. We discourage the use of satellite tracking strategies in azimuth and elevation to avoid sharp changes in the azimuth direction. For example, The Doppler effect caused by the satellite's movement is increasing and anti-symmetric in the movement's direction, thus its maximum is located at the elliptical footprint boundary. By creating orbits with many satellites, we can reduce the maximum absolute Doppler effect, since the distance between satellites decreases, thus reducing the size of the covered area in the direction of movement. As for the relative angular speed, the effect is symmetric concave with a maximum in the satellite's shadow. Therefore, the maximum is not affected by the satellite's density. Using basic geometry, these maximum values are given by the following formulas The maximum Doppler caused by the satellite's movement can be simplified when The Doppler effect caused by the UT moving with the speed vector v u ∈ R 3 , for a random v ∈ S 2 which we are assuming to be uniformly distributed, is given by 3) Channel matrix: Due to the large distance and lack of medium interaction, a Rician geometric channel model with power β and Rician factor K r can be assumed [5]. We will also assume that the Doppler effect can be corrected and removed from the channel matrix. The UPA at the satellite can be described by steering vectors in the direction v ∈ S 2 , denoted as

A. Design of the coverage area
We consider a LEO constellation where satellites are distributed around N p orbital planes with an inclination θ op , each one with N s satellites. Each satellite illuminates an elliptical area as shown in Fig. 1. From a geometric perspective, to parametrize the area covered by the satellite, we are choosing a stereographic projection respect to the Earth's center as illustrated in Fig. 4(a).
This allows us to use Cartesian coordinates while still considering Earth's intrinsic geometry.
Mathematically, this parametrization of the covered area translates the Cartesian point (x, y) into the sphere point R earth for the unitary vector definition u sat = xsat xsat , u x = vsat vsat and u y = u x ∧ u sat using ∧ to denote the cross product.
To design the sizes of the semiradiuses of the elliptical ROI, we will force the LEO constellation to cover the entire Earth. Thus, the equivalent flat problem we want to solve is to cover the entire Earth surface with ellipses centered in the vertices of a rectangular grid with spacing πR earth /N s in the first dimension and πR earth sin(θ op )/N p in the second dimension.
April 23, 2021 DRAFT Due to ellipses being convex and symmetric in both axes, this is equivalent to design ellipses containing the center points of the rectangular grid, as illustrated in Fig. 4(b). This condition can be written as An example of a semi-radius pair that fulfills this condition is

B. Beam codebook design
Our goal is to design a set of analog beampatterns or beam codebook, such that the set of beam footprints covers an ellipsoidal ROI as that described in the previous section. Depending on the location of the active users, some active beams will be selected from the set available in the codebook. We will consider a DFT-based analog precoder codebook structure, a popular solution when operating with planar arrays. In this case, the array response vectors in each dimension have a linearly increasing phase, like the columns in the DFT matrix. Therefore, the MIMO channel can be written in terms of array response vectors with this structure. DFT-based precoders match in this case the channel structure providing a good beamforming gain. Note April 23, 2021 DRAFT that the column vectors in a DFT codebook generate a grid of beams (GoB) that can be used to span all angular directions, as illustrated in Fig. 5 for a planar array with 16 panels (subarrays), and 7 × 5 antenna elements per subarray.
In particular, 5G New Radio defines what is called a Type II DFT codebook, where an oversampling factor is introduced in the definition of array steering vectors for azimuth and elevation [6], [7], [8], [9], [10], [11]. The goal of the oversampling factor, which is fixed and equal to 4 in 5G NR, is to mitigate the loss of beamforming gain in directions between orthogonal beams (straddling loss). This way, beam orthogonality is lost in general, and only beams spaced the oversampling factor are orthogonal. Fig. 5 shows the set of orthogonal beams marked in red when an oversampling factor equal to 2 has been used in each dimension. In our design, we will also consider a 2D oversampled DFT codebook for each subarray, although the oversampling factor O is not fixed, allowing the adjustment of the number of beams that illuminate footprints of different sizes for some given RF resources. The elements of the block diagonal RF precoding matrix F RF ∈ C N sat ×N b can be written in terms of this oversampled DFT codebook as for {(θ m,1 , θ m,2 )} m ⊆ { 2π(n−1) The oversampling factor can be used to adapt the RF resources to the size of the ROI and is the key parameter that drives the tesselation of the coverage area. By increasing O, the beam footprint is reduced, and the number of beampatterns pointing to the desired ROI increases.
The optimal value for O, in terms of provided SNR, allows the creation of as many beams pointing to the ROI as available RF chains. Note that when increasing the oversampling factor, the dictionary provides a higher SNR over the ROI because the straddling loss decreases, but the interference between adjacent beams increases, leading to a reduction of the SINR. Therefore, the oversampling factor has to be carefully selected to find the appropriate tradeoff between SNR and SINR. The inter-beam interference problem can be addressed by designing a suitable digital precoding matrix F BB that enables full frequency reuse, or by allocating different carriers to users served by adjacent beams, as proposed in [12].
Once the codebook is created, the cells are defined as the areas inside the elliptical ROI where a given beam provides a higher gain than any other beam in the available codebooks for all the subarrays. Note that the DFT nature of the designed codebook creates a rectangular tiling of the ROI. Section V includes examples of the beam footprints when this type of codebook is assumed on board.

A. Coverage
When operating with a transmitting beam-pattern f ∈ C N sat in the satellite beam codebook, and a receiving beam-pattern w ∈ C N UT , the RSS is given by Note that this expression depends on the user and satellite positions in x u and x sat respectively.
The different terms in (15)  To compute the combining gain at the UT side, we will consider the LoS direction to be known (in the system model we have assumed that the satellite position is known) and the Rician component of the channel to be unknown. Thus, the combiner w u will be computed as the one maximizing |<au(v),w>| 2 w 2 for the set of feasible w as discussed in Section 2.1. The average gain associated to this combiner can be computed as E( ). Since a R is Gaussian and independent of any other variable, we can write the average gain as Now, if we make use of the definition of w u and consider the Rician component to be uniform, i.e. Σ = I, we obtain the final expression for the receive gain as Since the LEO satellite is transmitting simultaneously through multiple beams with a gain G TX n , n = 1, . . . , N b , besides considering the RSS or the SNR, we also need to understand the impact of interference when computing the coverage. To this aim, we will compute the interference as the RSS of the signals transmitted through undesired beam-patterns with G TX n selected the gain corresponding to the desired beampattern and G interf RX = n =n selected G TX n . With these definitions, the signal to interference plus noise ratio (SINR) has the expression .

B. Throughput and Spectral efficiency
We consider the throughput expression given by the capacity B DL log 2 (1 + SIN R) in bits per second. We also consider the throughput without interference B DL log 2 (1 + SN R) as a bound to what can be achieved with inter-beam interference mitigation techniques. Assuming that frequency division duplexing mechanisms enable beam sharing among close users illuminated by the same beam, we will also consider the spectral efficiency with interference log 2 (1 + SIN R) and without interference log 2 (1 + SN R) as performance metric.

V. SIMULATION RESULTS
Unless otherwise specified, we consider a LEO satellite orbiting at a 1300 km height, covering an ellipsoidal ROI with semi-radiuses 534.1 km and 170.5 km. This area has been designed x-axis (m) y-axis (m) 10 5 Rectangular mesh design x-axis (m) For the next simulations we will assume that the satellite phased arrays consist of 5 × 3 sub-arrays of size 12 × 24, with a total of 60 × 72 antennas. Note that while we are arranging the sub-arrays into a larger plannar antenna array structure, each subarray will be transmitting an independent signal. We found out that, for this number of antennas and RF-chains, with an oversampling factor of 1.2 we can adjust all the available beams to the size of the ROI. The remaining system parameters are summarized in Table I.
To understand the impact of different antennas at the user terminal, we will consider both the UPA and the leaky wave antenna (denoted as CRLH-LW) to have a number of elements equal to 24 × 24. Beamforming with the planar antenna array at the user side introduces a gain of 27.6 dB. The gain provided by the leaky wave antenna is a function of the leakage factor, as illustrated in Fig. 7(a). The dependance of the antenna developed by Kymeta with the elevation angle [18] is also shown in Fig. 7(a). We consider now the effect of beam miss-alignment on the gain through Monte Carlo simulation for the phased antenna array and the leaky wave antenna. SINR when using the best beam for each cell.
shown in Fig. 7(c). It can be seen that a hybrid architecture offers a much better performance than an analog one under phase quantization, but both perform very similarly when the phase shifters have a high resolution.
Next terms of SNR and SINR can be seen in Fig. 8. While the SNR shows a good coverage of the ROI for an oversampling factor of 1.2, when visualizing the SINR we can observe areas with high interference. This is a consequence of using the DFT-type codebook employed in 5G NR, since its associated rectangular mesh creates points that are covered with the same gain by 4 different satellite beams. When selecting one of the beams, the other three will create a strong interference, leading to a SINR not larger than −4.77 [dB]. Furthermore, the design of static codebooks will always lead to points in which at least 3 beam-patterns have equal gain, thus leading to a SINR never larger than −3[dB].
In Fig. 10 we see the SNR variation over each axis for the lines with maximum and minimum SNR. The curves with minimum SNR corresponds to the boundary of the elliptical ROI, while the curve for maximum SNR corresponds to the line y = 0 in the coverage map. It can be observed that the gain difference between the center of a cell and the edge is in the order of 3dB for both axes if the boundary of the ROI is excluded. This SNR and SINR translate into the spectral efficiency shown in Fig. 9. Note that there is a significant inter-beam interference when operating with DFT-type codebooks as proposed in LTE and 5G NR. Future work needs to be developed to design alternative analog codebooks and hybrid beamforming strategies at the satellite side to reduce inter-beam interference.  Finally, we analyze the SNR degradation over time due to the satellite's movement for a user initially located at the center of the elliptical ROI. We have considered a beam selection strategy based on maximizing SNR. Fig. 12(a) shows an example of different satellite positions and associated selected beams. Fig. 12(b) shows the SNR corresponding to the different satellite positions illustrated in (a). It can be concluded that, for this codebook design and system parameters, the user is illuminated by the same beam for a period of around 30 seconds. After that, another beam has to be selected to keep a proper gain. Note that even after updating the best beam, there is a non-negligible gain loss due to the satellite movement, as also illustrated in Fig. 12(b).

VI. CONCLUSION
Research on integrating LEO-based satellite networks with terrestrial networks is still in its infancy. We implemented a complete simulator of the downlink of a massive MIMO LEO SatCom system to understand the physical layer challenges and evaluate a hybrid beamforming strategy based on an analog codebook. We designed the satellite footprint as an elliptical ROI given the number of orbital planes and number of satellites per orbit. We also designed a static DFT-type codebook, as commonly used in current cellular standards, to adjust the beamwidth and RF resources to the size of the satellite footprint by means of an interpolation factor. Simulation results show the relevance of the antenna terminal design to achieve an acceptable SNR. Our numerical experiments also show that DFT codebooks introduce a high inter-beam interference which reduces the ability of the system to provide a high spectral efficiency at all locations inside the ROI. Designing fully or partially connected hybrid beamforming architectures with per antenna power constraints and a digital precoder that mitigates inter-beam interference is an open research challenge that needs to be addressed to overcome the SINR performance of current approaches. The design of low complexity dynamic codebooks that better adapt the resources to the specific location of the users and reduce the SNR loss due to the satellite movement is a critical challenge. Our simulator also shows that the variation of the relative angular speed between the satellite and the UT is smooth when working with spherical coordinates. This leads to the need of designing beam tracking methods for the UT that track the satellite position rather than its variation in elevation and azimuth. Our numerical experiments also show the effectiveness of exploiting satellite position information at the UT side to avoid the need of channel estimation and tracking, which will provide old estimates due to the long transmission delay. Finally, an accurate model for the error in the satellite's position and trajectory is also necessary to better understand its impact in the performance of the satellite beam tracking strategies.