We assume that the satellite is equipped with a large uniform planar array (UPA) with dimension N x sat × N y sat and N x sat × N y sat 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 (Angeletti and De Gaudenzi, 2020). Our choice is a partially connected hybrid MIMO (Multi-Input Multi-Output) architecture with a fixed codebook-based analog stage, as a practical trade-off between complexity and reconfigurability. In this regard, a fully analogue reconfigurable beamforming network entails a large mass and power consumption, whereas a fully digital solution is costly and power hungry due to all the required connections. The advantages of the hybrid approach become especially relevant for very large arrays, with major efforts currently underway to develop massive hybrid arrays based on multiple subarrays operating with analog beam-forming (De Gaudenzi et al., 2020). The number of beamforming ports in each dimension is defined as N x RF and N y RF making a total of N RF = N x RF N y RF RF-chains. The large array is organized into sub-arrays of size N x sub × N y sub , with N x sub = N x sat / N x RF , N y sub = N y sat / N y RF as illustrated in Figure 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 consisting on N b beam-patterns F ∈ ℂ N sat × N b can be written as F = F RF F BB , where F BB ∈ ℂ N RF × N b is the digital precoder and F RF ∈ ℂ N sat × N RF is the analog precoder. Due to power constraints in the mth RF-chain, we have that | | [ F BB ] m , : | | 2 ≤ P RF . We define the index function I RF : { 1 , … , N sat } → { 1 , … , N RX } , where I RF ( n ) is the index of the RF-chain to which the nth antenna is connected. With this definition, the RF precoder F RF can be written as [ F RF ] n , m = { e ϕ n i if  m = I RF ( n ) 0 otherwise . 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 (Angeletti and De Gaudenzi, 2020).

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 x UT × N y UT , having a total of N UT = N x UT N y UT antenna elements. Recent advances on antenna technology enable the commercial applications of these planar arrays for SatCom terminals (Sikri and Jayasuriya, 2019). The second considered solution is a flat panel antenna based on leaky wave technology (Mehdipour and Eleftheriades, 2014). The final option is a commercial liquid crystal based metasurface antenna developed by Kymeta (Mehdipour et al., 2019).

We define w ∈ ℂ 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 ∈ ℝ N RF × 3 , being [ k u ] n , : the nth antenna element’s relative position. In this case, the UT steering vector a u : S 2 → ℂ N UT is defined over the unitary sphere S 2 as [ a u ( v ) ] n = e − i 2 π λ 〈 [ k u ] n , : , v 〉 . The maximum gain in a direction v ∈ S 2 can be computed as max w | 〈 a u ( v ) , w 〉 | | | w | | . The set of feasible beam-patterns for a uniform planar array is { z ∈ ℂ : | z | = 1 } N UT , and its maximum achievable gain in any direction is 10 log 10 ( N UT ) when each one of the antenna elements is assumed to provide a unitary gain.

For a leaky wave antenna with leakage factor α, a feasible beam-patterns w ∈ ℂ N UT follows the expression [ w ] n = e − 1 λ ( [ k u ] n , 1 α + [ k u ] n , 2 α ) − i λ ( [ k u ] n , 1 θ 1 + [ k u ] n , 2 θ 2 ) for θ 1 , θ 2 ∈ [ − π , π ] (Mehdipour and Eleftheriades, 2014). When α ≠ 0 , the gain of the leaky wave antenna can be computed as 10 log 10 ( ( 1 − e − α N x UT ) ( 1 − e − α N y UT ) ( 1 + e − α ) ( 1 + e − α ) ( 1 + e − α N x UT ) ( 1 + e − α N y UT ) ( 1 − e − α ) ( 1 − e − α ) ) . Note that there is a gain threshold 10 log 10 ( ( 1 + e − α ) ( 1 + e − α ) ( 1 − e − α ) ( 1 − e − α ) ) 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.

We consider the channel between a single LEO satellite located at position x sat ∈ ℝ 3 and a UT located at position x u ∈ ℝ 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 v sat = G m earth R earth + h sat , (1) 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 = v sat r earth + h sat . 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 π w sat .

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 Figure 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 Figure 4A. 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 R earth u sat + x u x + y u y R earth 2 + x 2 + y 2 for the unitary vector definition u sat = x sat | | x sat | | , u x = v sat | | v sat | | 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. 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 Figure 4B. This condition can be written as ( π R earth / N s R x ) 2 + ( π R earth sin ( θ op ) / N p R y ) 2 ≤ 1. (12) An example of a semiradius pair that fulfills this condition is R x = 2 π R earth / N s ,     R y = 2 π R earth sin ( θ op ) / N p . (13)

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 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 Figure 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 (Miao et al., 2018; Ericsson, 2017a, b; Samsung, 2017; Technical Specification Group Radio Access Network, 2017a, b). 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. Figure 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 ∈ ℂ N sat × N b can be written in terms of this oversampled DFT codebook as [ F RF ] n , I RF ( n ) = e − i 2 λ ( [ k sat ] n , 1 θ I RF ( n ) , 1 + [ k sat ] n , 2 θ I RF ( n ) , 2 ) , (14) for { ( θ m , 1 , θ m , 2 ) } m ⊆ { 2 π ( n − 1 ) O N x sub } n O N x sub × { 2 π ( n − 1 ) O N y sub } n O N y sub . Note that the m-th element of the codebook is contained in the non-zero entries of the mth column of F RF , that is [ F RF ] n , m = e − i 2 λ ( [ k sat ] n , 1 θ m , 1 + [ k sat ] n , 2 θ m , 2 ) . This expression is also proportional to the steering vector in the direction v ∈ S 2 with [ v ] 1 = θ m , 1 π and [ v ] 2 = θ m , 2 π . To create the codebook, we prune the complete oversampled DFT dictionary, keeping only those beam-patterns pointing to the elliptical ROI. 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 (Angeletti and De Gaudenzi, 2020).

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 5 includes examples of the beam footprints when this type of codebook is assumed on board.

When operating with a transmitting beam-pattern f ∈ ℂ N sat in the satellite beam codebook, and a receiving beam-pattern w ∈ ℂ N UT , the RSS is given by R S S [ W ] = | w H Hf | 2 = | γ | 2 | < a u ( x u − x sat | | x u − x sat | | ) + 1 K r a R , w > | 2 | < a sat H ( x sat − x u ‖ x sat − x u ‖ ) , f > | 2 . (15)

Note that this expression depends on the user and satellite positions in x u and x sat respectively. The different terms in Eq. 15 can be identified with the terms in the link budget equation. Thus G RX [ dB ] = 10 log 10 ( | < a u ( x u − x sat | | x u − x sat | | ) + 1 K r a R , w > | 2 ) , the transmitter gain G TX [ dB ] = 10 log 10 ( | < a sat H ( x sat − x u | | x sat − x u | | ) , f > | 2 ) , and 10 log 10 ( | γ | 2 ) comprises all the losses L P cable [ d B ] + L P at [ d B ] + L P fs [ d B ] .

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 | 〈 a u ( 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 ( | 〈 a u ( v ) + 1 K r a R , w u 〉 | 2 | | w u | | 2 ) . Since a R is Gaussian and independent of any other variable, we can write the average gain as E ( | 〈 a u ( v ) + 1 K r a R , w u 〉 | 2 | | w u | | 2 ) = E ( | 〈 a u ( v ) , w u 〉 | 2 ‖ w u ‖ 2 ) + 1 K r E ( | 〈 a R , w u 〉 | 2 | | w u | | 2 ) . (16) 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 E ( max w | 〈 a u ( v ) , w 〉 | 2 ‖ w 2 ‖ ) + 1 K r . (17) Since the LEO satellite is transmitting simultaneously through multiple beams with a gain G n TX , 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 I [ d W ] = P TX [ d B W ] − L P cable [ d B ] + G TX interf [ d B ] − L P at [ d B ] − L P fs [ d B ] + G RX [ d B ] (18) with G n selected TX the gain corresponding to the desired beampattern and G TX interf = ∑ n ≠ n selected G n TX . With these definitions, the signal to interference plus noise ratio (SINR) has the expression S I N R = R S S [ W ] I [ W ] + σ 2 [ W ] . (19)

We consider the throughput expression given by the capacity B DL log 2 ( 1 + S I N R ) in bits per second. We also consider the throughput without interference B DL log 2 ( 1 + S N 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 + S I N R ) and without interference log 2 ( 1 + S N R ) as performance metric.