Sparse Satellite Constellation Design for Global and Regional Direct-to-Satellite IoT Services

In this article, we introduce and design sparse constellations for direct-to-satellite Internet of Things (DtS-IoT). DtS-IoT does not require a ground infrastructure, because the devices are directly connected to low earth orbit satellites acting as orbiting gateways. The key idea of sparse constellations is to significantly reduce the number of in-orbit DtS-IoT satellites by a proper dimensioning of the delivery delay anyway present in resource-constrained IoT services and an optimal positioning of the orbiting gateways. First, we analyze long-range modulation (LoRa)/LoRaWAN and narrowband Internet of Things (NB-IoT) standards and derive realistic constraints on the maximum gap time between two consecutive passing-by satellites. Then, we introduce and optimize an algorithm to design quasi-optimal topologies for sparse IoT constellations. Finally, we apply our design to both global and regional coverage and we analyze the tradeoff between latency, number of orbit planes, and total number of satellites. Results show that sparse constellations can provide world-wide IoT coverage with only 12.5 and 22.5% of the satellites required by traditional dense constellations considering 3 and 2-h gaps. Also, we show that region-specific coverage of Africa and Europe can be achieved with only four and three satellites for LoRa/LoRaWAN and NB-IoT, respectively.


I. INTRODUCTION
The availability of a global Internet of Things (IoT) service for connecting millions of devices is currently a hot topic for both the industrial and academic communities [1], and the satellite segment can play a key role to achieve this goal [2], [3]. In particular, the direct-to-satellite IoT (DtS-IoT) paradigm consists of connecting constrained devices directly to low earth orbit (LEO) satellites without relying on a terrestrial infrastructure [4], [5]. The recent in-orbit demonstration of nanosatellites, such as LacunaSat-1 (from Lacuna, U.K.) [6] and Enxaneta (from Sateliot, Spain) [7] proved the feasibility of this concept, so far only addressed by experiments and scientific studies [8]- [11].
DtS-IoT is well suited to work alone where a pervasive infrastructure is absent or temporarily out of service, or to complement a terrestrial network deployment. Clearly, the choice of the protocol is a fundamental step for the success of DtS-IoT services, as a seamless interoperability with gateways already operating on ground would strongly encourage its success. A number of short, medium, and longrange protocols, known as low-power wide area networks (LPWAN), is part of the IoT family [12]. Among them, the power-efficient long-range modulation (LoRa) [13] and the supporting LoRaWAN network specification [14] (LoRa/LoRaWAN) are emerging as appealing candidates for realizing the DtS-IoT vision. As another example, narrowband Internet of Things (NB-IoT) is becoming one of the most important IoT solution, thanks to its integration with existing cellular networks and its extremely good performance figures [15]. The LacunaSat-1 and Enxaneta satellites implement LoRa/LoRaWAN and NB-IoT demonstrators, respectively.
In general, an IoT service with persistent global or regional satellite coverage would require hundreds of orbiting gateways. For example, LacunaSat-1 is the first in a planned constellation of 240 nodes. Even with the important reduction of launch costs [16], the deployment of such dense constellations requires significant economic resources, making difficult the involvement of small/medium companies interested in developing the technology.
In this article, we study the new paradigm of sparse LEO constellations for DtS IoT services. The key idea is to balance the protocol latency introduced by the gap time between two consecutive passing-by satellites and the constellation cardinality. To exploit this tradeoff, we first derive realistic gap constraints for the protocols under study. Then, we develop an efficient algorithm able to achieve an optimal positioning of the orbiting gateways satisfying both the latency and the geographic constraints.
The rest of this article is organized as follows. An extended survey-class overview of state-of-the-art LoRa/LoRaWAN and NB-IoT capabilities is presented in Section II, from where we derive a novel vision of sparse DtS-IoT network constraints detailed in Section III. As an example, by allowing coverage gaps of 2 h, bounded by the Class-B mode defined in the LoRaWAN specification, the approach can significantly reduce the required in-orbit infrastructure at the expense of higher data delivery latency. This phenomenon can be further exploited in NB-IoT, where sophisticated network-negotiated power-saving modes (PSMs) provide unprecedented flexibility to accommodate coverage gaps up to 3 h. Thus, a properly designed sparse constellation could comply with standard protocol modes, indicating that space-terrestrial IoT integration can be achieved with fewer resources.
Section IV presents specific optimization techniques to design sparse DtS-IoT constellations for global and local service regions. With respect to [17], where sparse constellations were briefly introduced (but with reference to global LoRa/LoRaWAN coverage only), we describe in detail an extended gradient descent methodology to efficiently derive quasi-optimal orbital parameters for both regional and global sparse IoT constellations based on any of the key IoT technologies (LoRa/LoRaWAN and NB-IoT). The technique exploits specific ad hoc heuristics to accelerate the convergence toward the minimum satellite fleet size that provides a maximum coverage gap of arbitrary and configurable times. Furthermore, we present an extension of the model to deal with regional coverage definitions, so that specific areas of the planet are served with even smaller constellation sets. To this end, we study rectangular bounding-box and hexagonal meshing (tiling) models to represent the limited coverage area.
In Section V, we apply and compare these design algorithms with the global and regional coverage problem, while assessing the resulting constellations through key coverage and algorithmic complexity metrics. A number of case studies with varying time gaps adjusted to LoRa/LoRaWAN and NB-IoT is presented. Results on worldwide coverage show that 2-h gap LoRa/LoRaWAN DtS-IoT services can already be provided in a global scale with only nine satellites deployed in the right polar orbit, while 3-h gap NB-IoT requires only six satellites. For regional coverage, we show that only three-to-five satellites are needed for most cases. For example, covering Africa demands only three satellites and regions, such as Greenland can be served with just two LEO satellites. Conclusions based on these and other cases as well as future perspectives from this appealing research area are finally summarized in Section VI.

A. Space-Based IoT-Like Technologies
In the space context, application-specific, low-power and low-throughput data transfer technologies targeting thousands of devices have existed for several years. However, these have not been prominently addressed under the IoT naming, a rather recent terminology. One of them is Argos [18]. Established in the late 1970s, Argos is a deviceto-satellite protocol designed to transport small amounts of environmental data. Seven satellites currently provide one-way (uplink) or two-way telemetry and telecontrol to more than 22 000 weather stations and buoys. Similarly, systems, such as Automatic Packet Reporting System [19] have been used for search and rescue among other applications. Furthermore, Satellite Automatic Identification System (S-AIS) [20] and Satellite Automatic Dependent Surveillance-Broadcast (ADS-B) [21] are popular space protocols fitting the IoT standards. Both systems use satellites to collect tracking data from vessels on the ocean or aircraft worldwide. AIS and ADS-B transceivers are currently on-board of the Orbcomm and Iridium satellite constellations, among other smaller providers.
Although these systems are currently operational, they are not prepared to be integrated with emerging state-ofthe-art terrestrial IoT technologies. Instead, we claim that a truly satellite-based IoT service will succeed (e.g., cheaper, power-efficient and feature-reach devices, and faster development) if it can profit from the economy of scale of terrestrial IoT technologies, currently led by LoRa/LoRaWAN and NB-IoT.

B. Direct-to-Satellite IoT
The IoT ecosystem is broad [see Fig. 1(a)]. It embraces short-and medium-range (1-10 km) cellular and sensor networks providing both high and low data rate services (0.1-100 Mbps) [22]. Moreover, and in accordance with the connectivity trend, low-power wide area (LPWA) technologies are also a crucial component of IoT capable of connecting applications that only need to transmit small amounts of information from long distances (100 km, at <50 Kbps) [23]. These applications range from agriculture to smart grid, environmental monitoring, emergency management, and others [24].
In order to further augment IoT coverage to provide a true global connectivity, satellite systems are natural candidates [2]. Compared with geostationary satellites, LEO satellites, orbiting below 1000 km height, can establish links with devices on the surface at reduced power budgets and round trip time delays [25]. However, the dynamics of near-earth orbits include velocities of ∼7 km/s, which results in overflights of <10 min over a given spot on the surface. These conditions demand constellations of several LEO satellites to achieve continuous coverage [26].
Existing LEO networks, such as Iridium [27] and Orbcomm [28], originally deployed to provide voice and data services, were already adapted to transport machine-tomachine data traffic [26]. Upcoming megaconstellations, such as OneWeb [29] and Starlink [30], will also include IoT services in their product portfolio [31], but the series of bankruptcies from the 1990s [32] suggest this is indeed a stressed business model. On the other hand, new startup companies, such as Kepler [33], Astrocast [34], and Lacuna [35], are deploying specific satellites constellations for IoT [36]. Because of the low-power and small IoT messages, nanosatellites are a great fit, enabling lower costs than traditional satellites and lower risk deployments. This concept has attracted the attention of both the industrial and academic community [37].
One approach to global IoT services is to use satellites as a backhaul to transport information from gateways deployed on ground that indirectly relay data from nearby IoT devices [38] [see Fig. 1(b), left-hand side part]. Indirect-tosatellite IoT is appealing outside highly accessible, dense, urban areas, where a reasonably high concentration of end devices justifies the deployment of dedicated IoT gateways to serve them. This configuration can thus leverage LEO or GEO satellites to backhaul the data to and from the core IoT network. On the other hand, applications on less accessible regions (i.e., oceans, mountains, and poles) might not justify or even hinder the deployment of IoT gateways. In this context, a more appealing but challenging architecture implies a DtS-IoT, where the IoT device directly transmits data to the passing-by satellite [4], [5] [see Fig. 1(b), right-hand side part]. Unlike the indirect approach, DtS-IoT motivates research on custom ground to space link protocols [39], [40]. However, leveraging features from existing LPWA standards would not only profit from the growing IoT market potential, but can also favor seamless interoperability with existing ground IoT infrastructure.

C. LoRa and LoRaWAN
Among LPWA technologies, LoRa is a chirp spreadspectrum modulation technique [13] particularly appealing for satellite as it provides very low power consumption and large link margin [23]. LoRa operates within the unlicensed industrial, scientific, and medical (ISM) frequency bands, which can also be an advantage in terms of meeting, otherwise, complex worldwide frequency regulations [41]. It is important to note that ISM bands are not allocated to space services (satellites are not defined as ISM applications), which still requires appropriate International Telecommunications Union (ITU) licensing and coordination [42].
The deployment of LacunaSat-1, the first LoRa/LoRaWAN nanosatellite by Lacuna in 2019, is compelling evidence of the feasibility of DtS-IoT, as expressed in related experiments [6] and scientific work, including long-range evaluations [8], modulation enhancements for LEO links [9], Doppler effect assessments [10], and adaptations for LoRa-based space links [11].
Running on top of the LoRa modulation, LoRaWAN is the specification responsible for the network layer service, enabling data handling over asynchronous bidirectional link protocols [14]. LoRaWAN is based on a star topology where gateways act as single-hop data concentrators bridging LoRa/LoRaWAN devices and a centralized network server in the Internet [43].
The specification render data rates from 0,3 kbps to 50 kbps via three classes of services, or device classes. 1) Class-A is the baseline, mandatory, and most powerefficient mode where the device turns the radio ON for the exact time needed to perform a frame transmission. This is an asynchronous pure ALOHA-based protocol [44]. Reception is supported by a receive window opened two times exactly 1 and 2 s after the end of the device transmission. 2) Class-B devices, typically equipped with batteries, allow gateway-triggered downlinks and uplink by periodically listening to so-called beacons (and pings in-between beacons) that synchronizes and coordinates reception and transmission episodes. This is a synchronous mode [45]. 3) Class-C devices, assumed with external power source, operate in continuous reception mode.

D. Narrowband Internet of Things
Another key LPWAN technology for low-power radio access is NB-IoT, part of 3GPP's Long Term Evolution (LTE) Release 13 [46] and 5G standards [47]. NB-IoT is capable of reusing and integrating with existing cellular network infrastructure and frequency allocations. Unlike LoRa's spread-spectrum modulation, it uses the same multiple-access technology than LTE: orthogonal frequency division multiple access in the downlink, and singlecarrier FDMA in the uplink, leveraging existing cellular channels.
To register with the NB-IoT network, the user equipment (UE) must first detect and synchronize with a g-NodeB (gNB) through the Narrowband primary and secondary synchronization signals. Then, after acquiring the narrowband physical broadcast channel, it receives the master information block. At this point, it can decode the narrowband physical downlink control channel (NPDCCH) and the narrowband physical downlink shared channel (NPDSCH) to obtain the system parameters.
Once this information has been successfully obtained, the device performs a random access procedure by contending for the cell's narrowband physical random-access channel (NPRACH). The access is based on slotted ALOHA in which users repeatedly transmit a random-access preamble within a time slot specified by the cell.
After accessing the NPRACH resources, it monitors the NPDCCH and NPDSCH, waiting to receive the random access response with uplink scheduling information. If the device successfully receives a reply, it will access the narrowband physical uplink shared channel to transmit its identity and a Radio Resource Control (RRC) Connection Request, among other information; otherwise, it will perform an exponential back-off and try again. When the user equipment (UE) receives a contention resolution message and an RRC connection setup from the gNB, it will acknowledge with an RRC Setup Complete, concluding the random access procedure. The UE can now proceed with authenticating and attaching itself to the network and, depending on the device, establishing a packet data network (PDN) connection to the network.
At the end of the network entry procedure, the UE will be registered and in the RRC Connected state, ready to transmit and receive data in coordination with the network. It is the core network that grants access to specific resource slots, specified and defined in time and frequency domains [48]. As a result, UE access to the radio channel is scheduled, mitigating interference at the cost of signaling traffic. This is, indeed, a tradeoff between complexity and efficiency (NB-IoT is based on stricter and more chatty device-tonetwork interaction than LoRa/LoRaWAN).
NB-IoT standards [46] define two classes of devices, NB1 and NB2, which can operate in 1) half-duplex frequency division duplex or 2) time division duplex modes. They use a 180-200 kHz bandwidth that can be split into multiple narrower tones of 3.75 and 15 kHz. Each device is associated with three possible power classes: +14 dBm (Class 6), and +20 dBm (Class 5), and +23 dBm (Class 3), enabling a maximum data rate of 150 kbps in the uplink (three times as much as LoRa/LoRaWAN).
To the best of our knowledge, the evaluation of LoRa/LoRaWAN service classes or NB-IoT features have not yet been properly conducted in the context of sparse DtS-IoT networks.

III. SPARSE DTS-IOT CONSTELLATIONS
The immediate approach toward a DtS-IoT is to deploy a LEO constellation with enough satellites to continuously cover the surface of interest, potentially the whole planet. This mimics the presence of at least one gateway to every device on ground. Any of the device classes from NB-IoT and LoRa/LoRaWAN can, in principle, be applied to such a dense constellation, whose orbital topology design does not differ from those currently providing voice and real-time data service (i.e., OneWeb and Starlink). In this work, we introduce an alternative vision: a sparse constellation based on LoRa/LoRaWAN and NB-IoT service classes.

A. Sparse LoRa/LoRaWAN Constellations
Being the basic and mandatory mode, the asynchronous Class-A LoRaWAN has received most of the community attention [45]. However, due to the device-initiated approach, there is no energy-efficient way for Class-A devices to effectively synchronize with satellites. Instead, we claim that devices implementing the synchronous Class-B mode can exploit the network-initiated downlink to become aware of the availability of a passing-by satellite. In practice, this occurs when the beacon can be successfully decoded [see Fig. 2(a)]. LoRa/LoRaWAN gateways would fly on the satellites, which can be provisioned with approximated location information of the IoT devices to coordinate the uplink access by means of the standard beacon and ping mechanisms. Uplink collisions, a notable issue for satellites illuminating hundreds of thousands of devices [49], can, thus, be conveniently controlled by the LoRa/LoRaWAN gateway in orbit. Strategies to achieve this include the proper estimation of the device set under coverage [50], [51] and a subsequent transmission probability moderation under control of the passing-by gateway [52]. Furthermore, broadcast messages can be timely transmitted to many LoRa/LoRaWAN devices, when on sight, via the so-called multicast groups [14].
When the beacons are no longer received (satellite no longer on sight), the LoRaWAN specification indicates a beaconless operation mode. In this mode, devices rely on their own clock to keep timing. However, due to the low-quality (low-cost) nature of IoT electronics, clocks are expected to quickly drift from the time reference. In this situation, the LoRaWAN specification states that compatible devices should be able to maintain beaconless mode up to 120 min (2 h) [14]. During this period, devices progressively increase their beacon listening window to compensate for the drift. Before the end of the period, another (or the same) satellite must make contact in order to adjust the device's internal clock, and reset the receive slots duration. As a result, the LoRa/LoRaWAN device can retain network-wide synchronization.

B. Sparse NB-IoT Constellations
There are inherent challenges to provide NB-IoT services using DtS-IoT. The reason is that NB-IoT is built upon cellular technology that traditionally requires fixed cells (assuming there is only user mobility). Especially in the IoT case, this approach is predominantly nomadic, since devices are assumed to remain under the same cell coverage for prolonged periods of time [53].
When satellites and high-speed dynamics are introduced into NB-IoT, it becomes difficult to perform resource allocation, to compensate for differential Doppler and phase impairments, and to satisfy the link budget and battery constraints involved in communicating with low cost, low complexity, and limited power devices [54]. There are also several handover, cell selection, tracking update and paging, and random access capacity issues in providing a massive number of users with NB-IoT services through nonterrestrial networks (NTN) and nongeostationary satellites. These topics are currently being investigated by researchers, industry, and space agencies in the context of 5G NTN standards [55]- [59]. Moreover, NB-IoT supports a limited number of users per cell (>52 547); thus, careful cell management is required. This is particularly true for satellite-based infrastructures, which must serve large coverage footprints with tens to thousands of beams.
Since satellites are constantly moving, cell movement is supported by 5G-NTN [60]. However, a fixed cell scheme in which subsequent satellites take over gNB responsibility is specified in 5G-NTN. Under this approach, each cell is associated with a fixed region of the earth's surface under coverage of a given satellite's beam. As the satellite moves, it steers its antenna beams to compensate its own movement and ensure a constant cell footprint. When the satellite is moving away from the cell and reaches a certain elevation angle, it hands over all RRC Connected UEs to an incoming satellite [61], [62]. Instead, with moving cells, users are gradually handed over between satellites as they detect a new cell with better quality indicators.
In the case of dense constellations, if more than one satellite is available at a time, then a satellite selection policy (i.e., satellite with the lowest Doppler rate or highest signal power) should be applied. Furthermore, to ensure continuous coverage and seamless handover for singlebeam UEs, satellite beams can overlap. Emerging satellite operators, such as Sateliot [7], have already successfully deployed 5G NB-IoT demonstrator satellites in orbit and obtained the required permissions to provide IoT services from space [63], demonstrating that solving these issues is feasible.
In our alternative sparse constellation vision, we exploit two fundamental power-saving NB-IoT features available to registered UEs after releasing a connection: extended discontinuous reception (eDRX) and PSM, illustrated in Fig. 2(b). In eDRX, a UE cyclically switches its receiver ON and OFF (up to 175 min per cycle). At the start of the eDRX cycle, the UE receives and tracks the NPDCCH waiting to receive a paging message for the duration of a paging time window (PTW) before sleeping until the end of the cycle. During this sleep, it is unreachable by the network. When it wakes up, it will listen to the NPDCCH for a PTW and go back to sleep for N (configurable) cycles before entering in PSM [64]. In PSM, the UE switches OFF after the so-called Active Timer expires (up to 186 min) and is considered by the network as powered OFF. However, the UE is kept registered, allowing it to resume the PDN connections without reattaching or reestablishing a new connection (as discussed in Section II-D, a chatty and energy-demanding process). The UE may spend up to 413 d in PSM state without losing context and will not be paged by the network. The use of eDRX and PSM, and their duration (including that of the Active Timer) must be negotiated with the network and the UE may wake up, notify the network of its location through a Tracking Area Update (TAU), and trigger an RRC resume request to resume the suspended RRC Connected state at any time without having to reattach. The exact way on which eDRX and PSM modes are combined is implementation dependent (i.e., a UE may use eDRX until the PSM Active Timer expires and then immediately switch OFF [65]).
For eDRX-and PSM-enabled UEs, we propose sparse constellations with revisit times shorter than the maximum PSM Active Timer duration (186 min). As a result, the network can exploit this sophisticated coordination mechanism to provide delay-tolerant services (i.e., to high-latency sensor networks) to mobile-initiated traffic while ensuring device reachability for mobile-terminated traffic (network initiated). This assumption ensures the UE will be reachable for two satellite passes on average (the pass in which it negotiates the sleep and the next) to allow it to receive last-minute command abortions (i.e., aborting sleep) before entering in PSM mode.

C. Sparse Constellation Design
The maximum beaconless (Class B in LoRaWAN) or power-saving (PSM in NB-IoT) coverage gap sits at the core of the proposed sparse IoT constellations paradigm. Instead of forcing continuous visibility to a dense satellite constellation, IoT devices could be left unserved and drifting for up to 2 h without detaching from the orbiting LoRa/LoRaWAN network, or up to 3 h in the case of NB-IoT services. Since most IoT applications are by definition asynchronous, they can tolerate delivery delays that would be prohibitive in traditional voice and data services. As a result, the number of satellites in a sparse constellation can be safely and significantly reduced, while still compliant with a flight-proven standardized DtS-IoT protocol, also widely adopted in ground systems.
In order to materialize the sparse DtS-IoT constellation vision, the adequate orbital parameters of the satellites need to be derived. This is a nontrivial constraint-programming optimization problem where 1) the maximum coverage gap G i, j must be kept under a given maximum value g max (maximum beaconless or power-saving times) in the T interval for a collection of ground latitude/longitude locations (i, j) in a set {L}, and 2) the size of the satellite fleet S arranged in P orbital planes carrying S/P satellites per plane shall be minimized. We are also interested in finding the inclination I of the orbits. Therefore, the search space is composed of tuples (S, P, I ), and we want to find those that lead to G i, j in which all points render revisit times lower than g max (e.g., g max = 120 min in LoRa/LoRaWAN Class-B, g max = 186 min for the Active Timer in NB-IoT). Table I summarizes the parameters and variables of the sparse constellation design problem. After presenting relevant considerations and assumptions, we introduce an algorithm to tackle the problem in Section IV. a) Constellation shape considerations: For the sake of simplicity, we initially assume that all satellites share the same orbital elements in terms of shape and size, as well as inclination and swath. Furthermore, we assume circular orbits (eccentricity e = 0), equally distributed orbital planes (in terms of ascending nodes), homogeneous distribution of satellites on each plane (satellite anomaly), and a plane phasing equal to the ascending node shift divided by the number of satellites per plane (known as Walker constellation pattern). We generalize this assumption to other constellation shapes in Section IV-A. b) Ground segment considerations: The sparse constellation design is governed by the maximum coverage gap g max , which is rooted in radio-access synchronization/power-saving aspects of LoRa/LoRaWAN and NB-IoT links between satellite and IoT devices. The connectivity from the LEO fleet to the ground segment (ground gateways and core IoT network infrastructure) is not involved in this issue, and thus, out of scope of this article. Nevertheless, it is worth mentioning that the ground segment in satellite IoT can also profit from relaxing the continuous visibility paradigm applied to dense constellations. IoT traffic is (in most cases) delay tolerant, which means data can rest stored on-board until delivered to the ground station (or device). Thus, although not part of the presented study, the reduction of the ground segment sizeat the expense of increased end-to-end latency-is another argument in favor of shifting toward a sparse paradigm in satellite IoT. c) Coverage region considerations: While most satellite IoT applications are designed for global coverage, a series of arguments favors regionalized/localized satellite IoT deployments. In particular, although remote regions without access to Internet infrastructure are distributed on a global scale, we envision satellite-based IoT networks launched into LEO orbits for the purpose of serving a subset of such areas. 1) The first argument is that we are living in a new space era in which the satellites (and ground supporting elements) are drifting from being expensive, long-lead items to cheap and quickly procured off-the-shelf assets. This enables architectural alternatives with degrees of scale and flexibility hitherto impossible [16]. Indeed, the cost reduction is already allowing governments and companies to deploy less ambitious satellite systems to serve specific regions of interest (even though the satellites' trajectories cover wider regions of the planet). 2) The second argument is based on region-specific regulations. At the time of writing, the definition of a suitable unique frequency band for the downlink of global-scale satellite IoT systems is still an open discussion topic among the LoRaWAN and NB-IoT satellite communities. As a result, early deployments will need to be designed to serve the areas for which the orbiting hardware was designed for. To accommodate regionalized satellite IoT use cases, the proposed design method should consider the maximum coverage gap only over arbitrary closed regions of the earth surface. These are typically expressed in shapefile format [66] from which {L} can be populated by comprising latitude/longitude points. d) Grid {L} considerations: In order to bound the (i, j) latitude/longitude grid resolution, the maximum coverage g max value can be leveraged. Let d P 1 g (l ) be distance that any point P 1 on the surface at latitude l shifts according to the rotation of the earth in g max time. There is no added information in measuring the coverage gap G P 2 of an intermediate point P 2 between the first and second location of P 1 . Thus, the longitudinal step of two points in {L} shall be no smaller than d P 1 g (l ) computed based on g max . 1 Furthermore, in the latitude dimension, the grid density can be adjusted by the cosine of the latitude to account for the increased density of points as they become closer to the poles [67]. e) Inclination I Considerations: Satellite coverage is known to be dependent of the inclination parameter [68]. In particular, Fig. 3 where L is the maximum latitude among all points in {L}, i is the orbital plane's inclination, and λ is the maximum offtrack ground angle (λ = s/2). This expression can be solved for i through numerical methods (we use the bisectional method), and it delivers a suitable lower bound inclination that can be used to initialize the following search heuristic.

IV. SPARSE DTS-IOT CONSTELLATION DESIGN
Based on the former parameters, assumptions, and considerations, we develop a gradient descent search algorithm that approximates the minimal S that can guarantee G i, j ≤ g max ∀i, j ∈ {L}.
To compute G i, j , each step of the algorithm propagates the orbit of all satellites in the constellation and subsequently performs an access evaluation for each point in the grid {L}. Since this can be a compute-demanding process, we care to preemptively discard noncompliant scenarios. We achieve this objective by making quick evaluations over shorter time horizons (T 0 T ), and progressively testing Fig. 4. Complexity-progressive gradient descent algorithm for sparse DtS-IoT constellation design [17].
each topology over {L} with increasing resolutions. In other words, simpler versions of the scenario are solved in an incremental fashion, and discarded as soon as the criterion is not met. The flow diagram of the proposed algorithm is presented in Fig. 4.
Step by step, the algorithm runs as follows. 1 Initialize: The algorithm takes as input (see Table I): antenna swath (s), satellite height (h), coverage region (grid {L}), maximum coverage gap (g max ), and the maximum time horizon (we find that T = 1 week is a suitable parameter). A solution space is then computed, comprised of 1) the possible number of satellites [s min , s max ] and planes [p min , p max ], where p max ≤ 360 • /s, and 2) an allowed inclination range [i min , i max ], where i min can be obtained by solving equation (1), 2 and i max = 90 • as we disregard retrograde orbits. Based on {L}, the steps N on which the input grid resolution will be approximated is determined. The shortened time horizon T 0 is computed as T 0 = g max + A max , where A max is the longest visibility period between a passing-by satellite and a point on ground, which can be computed by A max = (s × π a 3 /μ)/180, where a is the semi major axis of the circular orbit, and μ is the earth's gravitational parameter 3 (e.g., for 700 km height and s = 120 • , T 0 = 131.6 min).
2 Setup: The first scenario is generated from a new set L with a subset of the points in {L}, the coverage region. This region may be defined by its maximum and minimum latitude and longitudes or by shapes composed of coordinates, poly lines, or polygons, and bounded by a rectangular envelope. Results in Section V will discuss use cases for each of these region definitions.
Alternatively, to optimize the coverage gap for shapesensitive applications, such as regional surveillance and 2 i min can also be bounded by the launch site location.   Fig. 6 public safety, the envelope may be tessellated into hexagonal elements, constructing {L } set from the centroids of the elements that intersect the original shape.
After building the grid, the scenario is populated with the computed inclination I = i min , the minimum number of planes P = p min , and satellites S = s min . Furthermore, the first scenario is initially evaluated in a time horizon T 0 setting the so-called complexity level counter to zero. This counter indicates how many times the resolution of L is increased to reach the target resolution. We find that the reduced grid expression combined with the minimum time T 0 enables the direct elimination of most noncompliant scenarios (g n=1 > g max ) quite efficiently. Table II tabulates evidence on this regard based on the executions analyzed in Section V. 3 Simulate: The scenario is simulated for a time lapse T or T 0 for every first scenario for the tuple (S, P, I ). Orbits are propagated using the Simplified General Perturbations 4 (SGP4) algorithm [69] and accesses are computed for each of the points in the subset {L }. The new maximum gap G i, j is computed from the obtained gap intervals among all P points in the grid (G n = max(G P i, j |P ∈ L)). If the condition is met (G n > g max ), the resolution of the scenario is augmented; otherwise, the tuple (S, P, I ) is discarded and a new neighbor is allocated from the solution space as follows. 4 Obtain new neighbor: If the scenario was discarded, a new neighbor is obtained. The first criterion is to increment inclination on steps of one degree. If I = i max , then the next criterion is to increase the number of satellites per plane by one (S = S + P). Else, if S = s max , the final neighbor search criterion is to increment P by one. Else, if P = p max , the search over the complete search space is concluded. The increasing scenario complexity as the search advances avoids overshooting good quality solutions. If no more neighbors can be found, the algorithm terminates. 5 Increase resolution: If the scenario maximum gap was lower than G i, j > g max , the resolution and complexity of the scenario is increased. First, if the time horizon was set to T 0 , it is now extended to the full T simulation time.
Next, the resolution of {L } is increased, thus also moving to the next complexity level. If the envelope was divided into rings, the resolution is increased only in the latitude dimension: new rings are included by splitting the already evaluated intervals in two; otherwise, the resolution of the mesh elements is increased simultaneously in longitude and latitude, exponentially increasing their number and resulting in a complexity penalty, but allowing the constellation design to be finely tuned to a given shape.
If the target resolution of {L } is reached (i.e., {L } = {L}), it means the current scenario fulfills the required coverage gap and the corresponding generating tuple (S, P, I ) is stored. It is indeed possible that several tuples comply with the coverage gap condition with the same number of Center-dotted red-squares indicate device coordinates.
satellites S. In this case, we honor those with the least number of planes P (which are related with deployment cost), and then, in the last place, lower inclinations I (which reduces launch fuel requirements).
To conclude, Fig. 5 highlights the difference between the grids constructed for Brazil. Without tiling, points are sampled uniformly across the latitude/longitude range according to the selected resolutions. Instead, with tiling, the method considers only the centroids of the hexagons that belong to the regional shape (in gray).

A. Constellation Shape
The presented complexity-progressive gradient descent algorithm assumes a Walker constellation shape-circular orbits, same period and inclination, and evenly spaced ascending nodes and satellites in each of the orbital planes. Thanks to the uniform orbital parameters, orbital perturbations affect each satellite in approximately the same way, reducing the dependency on station-keeping maneuvers (discussed and analyzed in Section V-A1) [70]. Walker constellations have been adopted for multiple missions, such as Iridium (deployed in the 1990s [71]) and modern megaconstellations (leveraging multiple shell configurations [72]). While we expect the main application scope of the introduced algorithm to rest in Walker formations, the logic can easily be adapted to accommodate other parameteric constellation shapes. To illustrate the flexibility of the approach, a case study on so-called flower constellations [73] with relaxed orbital circularity constraints is discussed in Section V-A2.

V. ANALYSIS
We implemented the sparse DtS-IoT constellation design algorithm in a specific Java-based application 4 that includes an SGP4 propagator released by the US Department of Defense (DoD), and improved as reported in Vallado's Revisiting Spacetrack Report number 3 [69]. Some tools of the Orekit open-source libraries were also leveraged [74]. Access intervals obtained by the implementation were contrasted against AGI's System Tool Kit software, rendering negligible differences in the order of few milliseconds in access times.

A. Global Coverage
Extensive simulations were performed for constellations devoted to global-scale IoT services with different values of elevation threshold in between 5 • and 40 • , which map to antenna swath s = 127.1 • and s = 87.3 • , respectively. We are also interested in studying maximum resolution {L} worldwide grids for which coverage is evaluated up to latitudes of 60 • and 80 • . All considered satellites are in a 700 km altitude, circular orbit, positioned in symmetrically distributed planes, RAAN, and mean anomaly in each plane. We let the solution space include one to five planes with one to eight satellites per plane (retrograde orbits not allowed). Obtained results are shown in Fig. 6 for LoRaWAN (with varying swaths angles from 127.1 • to 87.3 • ) and Fig. 7 for NB-IoT (with 127.1 • swath angle). In order to compare sparse and dense constellations, the former plots include information on the number of planes and satellites required to achieve a complete earth coverage using streets of coverage (a well-known pattern design for dense constellations deployed in polar orbits [68]).
In general, results from Fig. 6 confirm that the number of satellites increases with the reduction of the antenna swath, both for sparse and dense constellations. However, the difference between them is noticeable. With sparse LoRa/LoRaWAN constellations, up to nine satellites are needed, while 8 × 11 = 88 are required for a dense configuration (s = 127 • ), and when 40 are enough for sparse, dense constellations demands 28 × 31 = 868 satellites (s = 87 • ). To keep an adequate scale of the plot in Fig. 6, dense constellations are expressed in number of satellites per plane and number of planes, thus the multiplication to obtain the total satellites. Indeed, sparse constellations can operate with only 10% and 4% of the in-orbit infrastructure required for dense constellations at the expense of a worst-case latency of g max = 120 min for LoRaWAN and g max = 186 min for NB-IoT. This extra latency supported by NB-IoT reduces the constellation size by 33% to 45% (see  Fig. 7). For NB-IoT, five and six satellites (five and three planes) are sufficient for a sparse configuration, compared with the aforementioned nine satellites for LoRAWAN systems.
Results also show that the proposed approach is rather insensible to the maximum latitude. In particular, constellations 1 , 2 , and 3 (illustrated in Fig. 6 to give an intuition on the topology) are the best solution found for the same swaths in 60 • and 80 • maximum latitude values in {L}. Inclination, on the other hand, proved to be a more sensitive parameter to the maximum latitude, as already discussed in Section III-C-c). Specifically, the plots show that higher inclinations tend to be required to cope with smaller swaths, while some specific satellite fleets can meet the objective with exceptionally lower inclination values (e.g., s = 125 • with i = 78 • and i = 79 • ).
As an additional comment, plotted results correspond to orbital heights of 700 km, but observations from extended simulation campaigns showed that altitude is also quite insensitive. For example, the same S is obtained for all altitudes in the 400-1000 km range (s = 127 • ). At these altitudes, the orbital period is shorter than g max . Thus, the maximum contact gap is determined mainly by the time it takes for the orbital planes to rotate around earth. Also, for a fixed inclination and a given number of satellites, the constellation with the highest number of orbital planes will exhibit the lowest coverage gap.
1) Station Keeping: Station-keeping maneuvers are crucial to keep the formation of satellite constellations. Indeed, propellant consumption and budget are considered in traditional orbit design [75]. Maneuvers are also relevant in a LEO orbit with a growing debris count (i.e., space junk) that forces collision avoidance and end-of-life de-orbiting  operations. Moreover, efficient low-thrust electric propulsion systems are becoming increasingly popular as they enable maneuvers even for resource-and volume-constrained nanosatellite platforms [76]. With this in mind, an important aspect of the designed satellite IoT constellation is the time span that the resulting configuration can keep up with the g max constraint, without performing station-keeping maneuvers. In other words, we care about which is the maximum time period that the LEO satellites can drift due to orbital perturbations (external forces-earth oblateness, atmospheric drag, lunar gravity, and solar radiation pressure-not present in idealized orbital models). To study this phenomena, we leveraged the SGP4 drift model and applied it to three smallest global NB-IoT constellation candidates from Fig. 7.
Results from propagating the satellites until the g max constraint is no longer satisfied are tabulated in Table III. We observe that configurations with lower inclinations exhibit a higher sensibility to drift effects. In any case, we can expect 24-28 d of free flying without the need of station-keeping maneuvers to correct the constellation topology. Fig. 8 illustrates the impact of 28 d of uncorrected orbital perturbations on the ground track of six satellites in a Walker formation. Including the free-flying time as part of the optimization objective of the constellation design algorithm is not a trivial task and we leave it as a future research effort.
2) Non-Walker Constellations: As discussed in Section IV-A, the logic of the gradient descent algorithm can easily be accommodated for non-Walker (but parameteric) constellation topologies. In this section, we study an adaptation for Flower constellations, standing for the state-of-theart constellation patterns, which allow nonzero eccentricity (i.e., elliptical) orbits. An interesting property of Flower constellations is the repeating ground-track performance over a fixed number of days, which leads to a more uniform access distribution and simpler network operation. The core variables describing a Flower constellation are summarized  in Table IV. The interested reader is referred to the report in [73] for an in-depth discussion about the implications of these parameters. The gradient descent algorithm illustrated in Fig. 4 was adapted on step 4 (generation of new neighbors). Instead of performing the search on a space composed by (S, P, I ) tuples, we explore (F d, F n, N p, Nd, Ns, I ), a higher dimensional space. In order to quickly converge on a solution, we prune the search by applying symmetric phasing conditions, namely F n = F d and F d = Ns. Also, given that F h is an additional parameter introduced in later Flower constellations to change the satellite phasing, we set F h = 0 to further reduce the degree of freedom. Finally, the ascending node spread (RAAN) is fixed to 60 • . Also, following the Walker constellation case, we set at an antenna swath of 127.1 • . Regarding the initial parameters for this case study, we consider an eccentricity of 0.02 and a semimajor axis of 7150 km (apogee = 7293 km and perigee = 7000 km).
Results of the obtained constellation are summarized and illustrated in Fig. 9, together with the Walker topology that complies with the same coverage gap. In particular, the algorithm provided a Flower constellation comprising ten satellites (equal to the Walker case), with a slighter lower inclination value (60 • against 62.9 • ). However, the objective of this case study is to expose the flexibility of the design algorithm. An in-depth analysis of the applicability of Flower constellations for DtS-IoT services (and their corresponding optimization) is left for future work. B. Regional Coverage

1) Simultaneous LoRa/LoRaWAN and NB-IoT Services:
We first consider the 2-h maximum coverage gap and study the constellation design for two key regions: Europe and Africa (see Fig. 10). We also concentrate on selected countries that suffer from poor land-based network coverage in remote and rural areas [77], where DtS-IoT can play a major role.
Africa By applying the algorithm to the African continent, the advantage of using hexagonal tiling can be clearly perceived: it provides the shortest contact gap metrics (G i, j ∀(i, j)) and lowest inclinations (I) for all constellation sizes. Already with S = 4 satellites (P = 4 and I = 20 • ), the g max = 120 min gap constraint is satisfied, while at least S = 6 satellites (P = 6 and I = 35.82 • ) are required at the highest inclination. Increasing the constellation to S = 9 satellites (which renders an inclination equivalent to 25 • , a third of that of the global constellation) can cut the maximum coverage gap in half, to slightly over an hour. Finally, with S = 36 satellites, the maximum G i, j among all points on the tiled {L} is shorter than 15 min, suitable for emergency-like services over the whole continent.
Europe On the other hand, for the European region, the tiling effect on {L} is not so noticeable. Since Europe includes several islands and land masses at high latitudes, distributed across the rectangular envelope, the constellation design is pushed toward medium-to-high inclinations (I ∈ (60.67 • , 85.7 • )). However, thanks to the higher accuracy of the tiling, the algorithm was able to identify solutions that were discarded when using the simpler annular approach. In particular, since the heuristic prioritizes low -inclination solutions, solutions were biased toward candidates with a lower inclination but higher G i, j for a given constellation size. Nevertheless, as more satellites are added, there is a convergence toward a solution with 10 min MCG using 36 satellites ((S, P ) = (6, 6)). Table V shows the (S, P, I ) and MCG (max(G i, j )) for some countries and for a large but sparse region corresponding to the world's mountain ranges, including Antarctica (see Fig. 11). For comparison, we have selected solutions that minimize the number of satellites (the main optimization metric), and then configurations with minimum inclination. For these countries, satellites distributed across a number of orbital P planes equal to the number of satellites (i.e., {(S, P ) = (2, 2), (3, 3), . . ., ∀ S ∈ (2, 5)}) proves to be sufficient to comply with g max = 120 min. This shows that connectivity of sparse LEO constellations with 2-3 h revisit time constraint is mainly determined by the interplane spacing, and Algorithm resource consumption metrics in terms of compute time and memory consumption for the above corresponding cases (bottom). not the number of satellites in each plane. In other words, once a sufficient number of planes has been reached, the constraint is satisfied. Then, adding more satellites allows either to reduce the MCG coverage metric G i, j or the orbital inclination I. Indeed, results from Fig. 10 and Table V indicate that up to a 15-20 • reduction in inclination can be obtained, at the expense of 50 to 300% increase on the satellite fleet size.

Selected Regions
From these results, we conclude that the proposed algorithm can be quite sensitive to the regional shape, but the presence of tiling may be of limited value for some localized cases (i.e., Europe). Nevertheless, if designers also aim at minimising the coverage gap, then it is best to apply tiling. Without tiling, the algorithm may produce larger constellations at higher inclinations. This is an intuitive result that can be understood from Fig. 5: uniformly sampling points in the rectangular envelope results in the coverage grid being populated with locations outside the target region (in gray) and this effect is amplified for triangular regions because of the greater mismatch between the rectangular envelope and the region itself.
Focusing on the algorithm complexity metrics, defined in terms of compute time and memory usage, we observe that applying tiling results in a 41.18%-158.46% increase in execution time and 74.95-% 258.43% greater memory usage, especially at high complexity levels (> 2). As the tiling resolution of {L} doubles in latitude and longitude, there is a quick exponential growth in execution time due to limited parallelism. There is also a considerable penalty in computing the tiles that intersect the coverage region. The greater its detail, the more computationally intensive this step is; thus, region contours should be defined carefully. Memory usage, on the other hand, barely increases across complexity levels since the coverage gap is computed one point at a time. The main memory requirement is storing all propagated satellite trajectories over the simulation interval. With tiling, there are also complex, nested, data structures being used to keep track of the regional shape and hexagonal tiles, resulting in a higher memory consumption.
Although not detailed, it is important to highlight that with tiling, the algorithm converges to viable solutions at lower complexity levels (≥ 2) because it only samples points that intersect the regional shape. Therefore, the MCG metric does not vary abruptly as the resolution of {L} is increased. As a consequence, tiling requires more memory but may be comparable in overall execution time.
2) Exclusively NB-IoT Services: Having assessed the combined LoRa/LoRaWAN and NB-IoT services, we study and compare the flexibility that the eDRX and PSM powersaving modes of NB-IoT brings to the constellation design problem. In particular, we assess the impact of the approximately 3-h maximum coverage gap (g max = 186 min) for Europe and Africa.
To begin, the satellite fleet is reduced to three satellites for both cases, but only for Europe there is an increase of 20 • in inclination, as shown in Fig. 12. An additional result from simulations is that combined Africa-Europe coverage can be achieved with S = 5 satellites (P = 5 and I = 60.7 • ), which is one more satellite than required by LoRa/LoRaWAN for Africa coverage, and three fewer satellites than required for combined Africa-Europe Coverage, highlighting the advantages of supporting higher latencies. Finally, for services with a latency constraint shorter than 1 h, at least 12 satellites in four orbital planes are required to service Europe and Africa.

C. Discussion
LoRa/LoRaWAN vs NB-IoT constellation size The obtained results for LoRa/LoRaWAN maximum coverage gap of g max = 120 min showed that a total of S = 9 satellites are required to provide global coverage (for h = 700 km up to latitudes of 80 • ) [17]. However, this constellation configuration varies when studying LoRa/LoRaWAN coverage on specific regions. For those regions extending beyond 80 • latitudes, the size of the constellation can rise up to S = 10, such as the world mountain case. Indeed, pole-topole coverage for sparse mountain regions is a challenging setup for an appealing use case we intend to address in the future. We expect that even smaller S values could be found when relaxing the uniform inclination assumption, a promising future research emerging from this article. But for more specific cases, such as Greenland, we found the fleet size can be reduced down to just S = 2 satellites. On the other hand, thanks to NB-IoT's PSM 3-h gap, the same global coverage mentioned previously can be achieved with S = 5 satellites. Most notably, Europe and Africa can be covered with P = 3 equally spaced orbital planes, each with one satellite (S = 3). Three planes and satellites are the minimum that allows us to cover the nearly 120 • arc with high inclination to serve targets at 80 • latitudes (as planes rotates around the earth) within constraints.

Cellular versus noncellular
The possible infrastructure reduction based on NB-IoT coverage gap flexibility is, of course, in tradeoff with the enhanced complexity of the enabling technology. In particular, while LoRa/LoRaWAN leverages a simpler ALOHA-based access where devices are not bonded to a specific gateway, NB-IoT implements cellular-based techniques supported by strong and strict interaction and negotiation with the network. Moreover, it is possible that a multibeam scheme with hundreds of steerable beams may be required to provide effective cell coverage over the earth surface, which might hinder the utilization of cost-effective Cubesats [78]. This flexibility renders a higher technological complexity for NB-IoT that might be detrimental on other system parameters, such as overall cost, cell-level scalability (i.e., NB-IoT supports up to 52k UE per cell), and power consumption [79]. Nevertheless, an efficient beam-management solution may enable the required system-level scalability to reach millions of users with guaranteed quality of service at higher rates, while reducing interference and contention probability. There is not a one-size-fits-all solution.
Traffic considerations Another crucial aspect of constellation design is the traffic load the fleet is expected to move between the user's devices and the ground segment. Since there is only one satellite in visibility of each user terminal in sparse constellations, under NB-IoT cell limits, thousands of cells per satellite might be needed in heavily demanded regions. But being remote areas the core focus of sparse constellations, accounting for low user density with reduced data demand could also relax the coverage need, and save valuable in-orbit infrastructure. Thus, even if data rate is low, and even under constrained duty-cycling/sleep time limits, the expected traffic load shall be considered. Including traffic model abstractions in the constellation design process is indeed an appealing research line emerging from this work.

VI. CONCLUSION
In this article, we have explored the reduction of the required in-orbit infrastructure to provide massive machine-to-machine DtS IoT services on a global and regional level, while leveraging standard LPWAN protocols. By exploiting key LoRa/LoRaWAN and NB-IoT features, such as beaconless and power saving mechanisms, we have designed and implemented an appealing optimal sparse constellation design algorithm to realize the vision of DtS-IoT at a fraction of the size of traditional satellite fleets. Extensive results cover a large set of realistic case studies. In particular, we have achieved global coverage with only six to nine satellites, and continental coverage with only three to four satellites. Such a significant reduction with respect to classical (mega-)constellations is possible thanks to the delay-tolerant nature of IoT communications. The provided low-cost sparse satellite configurations represent a very cost-efficient approach to extend existing LPWAN terrestrial networks. As future work, we identify the following two key lines in the algorithmic and protocol domains. 1) Algorithmic: We foresee the following three clear improvements and extensions to the gradient descent algorithm.
a) Constellation stability: We analyzed the stationkeeping implications, but formally including maneuvers as part of the overall goal of constellation design is left as future work. b) Constellation patterns: While we focused on Walker and on a case study of Flower constellations, an appealing research line is to explore and compare other non-uniform shapes for different service areas. c) Constellation traffic: The final goal of the satellite fleet is to transport data, and the geographical distribution of traffic demands [80] as well as more sophisticated application models (e.g., public safety, emergency services) can also be considered as a design constraint (possibly involving simulations on the loop).
2) Protocol: Authors are already working on the detailed evaluation of LoRa/LoRaWAN and NB-IoT focusing on collision probabilities, contention rates, and power consumption. Even though the obtained constellations provide the required maximum connectivity gap for both, the objective comparison of LoRa/LoRaWAN and NB-IoT in the satellite IoT context remains an open research question.