Nicola Accettura

Balakrishna J Prabhu

Ubiquitous multi-gateway LoRa networks: models and performance evaluation

Keywords: ALOHA, duty cycle limitation, scalability, performance analysis

LoRaWAN is a low-power and wide-range communication technology that is well suited for the demands of Internet-of-Things. It uses ALOHA with a duty cycle limitation as a channel access method to avoid complexity and communication overheads at the end devices as well as to ensure a fair access to the channel. However, given the large number of devices that are expected to integrate this network, scalability or the network capacity will be one of the fundamental issues that will need to be investigated before any large network deployment. This paper analyses the throughput of multi-channel multi-gateway LoRaWAN with duty cycle limitation. Throughput formulae for diverse scenarios including regularly tiled multi-gateway networks are derived and validated through simulations. Also considered is the throughput for applications that require a packet to be received by multiple gateways -for example, radio localization requires reception by three gateways -simultaneously to be successful.

I. INTRODUCTION

The Internet of Things (IoT) is rapidly taking shape with a sharp growth in the number of objects connected to the network. The majority of these objects rely on wireless links to communicate and are constrained by their low energy budgets [START_REF] Warneke | Smart dust: communicating with a cubic-millimeter computer[END_REF]. Various generations of radio technologies saw the light of day to meet the growing demand and the specific needs of these applications [START_REF]IEEE Standard for Low-Rate Wireless Networks[END_REF]- [START_REF] Mackensen | Bluetooth low energy (ble) based wireless sensors[END_REF]. Within these set of communication solutions, Low Power Wide Area Networks (LPWAN) [START_REF] Sanchez-Iborra | State of the art in lp-wan solutions for industrial iot services[END_REF] specifically target low power devices that communicate sporadically over long distances. Most LPWANs operate on unlicensed spectra, thus reducing the cost of deployment and increasing accessibility. They are based on an hierarchical architecture similar to that of cellular networks with gateways acting as hubs for the end devices [START_REF] Palattella | Internet of things in the 5g era: Enablers, architecture, and business models[END_REF].

An emerging LPWAN technology is LoRaWAN [START_REF] Lora | 0.3 Specification[END_REF], which leverages Chirp Spread Spectrum techniques [START_REF] Semtech | LoRa TM Modulation Basics[END_REF], [START_REF] Petäjäjärvi | On the coverage of lpwans: range evaluation and channel attenuation model for lora technology[END_REF] to achieve long range communication. The channel access is based on an ALOHA-like protocol, which reduces the cost of signalling and other control messages required in centralized channel access methods. One of the drawbacks of ALOHA is that there is no admission control on the end devices which can transmit a new packet as soon as they have finished transmitting one. This induces unfairness in the channel use, can increase collision rate and reduce goodput. To improve fairness, LoRaWAN imposes a duty cycle limitation on the end devices. For each packet transmitted, an end device has to remain silent for certain amount of time (specified in the standards) before it can transmit another packet.

While standard pure ALOHA networks have been studied quite extensively, the addition of a fixed duty cycle limitation has a non-trivial impact on the performance and requires an in-depth investigation. Questions related to the throughput of nodes, optimal transmit rates, optimal density of nodes, etc. for a given duty cycle would help the national standardization bodies to choose the right duty cycle limitation for their environment.

A. Contributions

The main contribution of this paper is the throughput analysis of multi-channel multi-gateway LoRaWAN based on the ALOHA channel access with duty cycle limitation. The throughput formula is first obtained for a single gateway network for fixed as well as for random number of end devices sampled according to a Poisson distribution. One of the salient features of LoRaWAN is that, unlike standard ALOHA networks where a packet in successful when it is received by at least one gateway, certain applications can require a packet to be received simultaneously by more number of gateways to be successful. For example, for radio localization, a packet has to be received by at least three gateways for it to be successful. We obtain the throughput formula in the multigateway network with the additional requirement of a packet being received successfully by at least a certain number of gateways. The formulas are derived for the square, triangle and honeycomb architecture of the gateways. These formulas are then validated by extensive simulations carried out using an event-driven LoRaWAN simulator developed in Python.

B. Organization

In what follows, Sec. II overviews the related works, while III provides some relevant background on the LoRaWAN technology. Then, Sec. IV pictures a general probability model for single-gateway LoRaWAN networks, that extends the basic pure ALOHA one to account for the availability of multiple channels and duty cycle limitation enforcement; given the importance of designing a network to target an average density of sensor devices spread randomly on the ground, a model that depends on such a density value is also presented. To cope with the need of designing multi-gateway LoRaWAN deployments, a further step ahead is done in Sec. V: a general formulation targeting such a goal is presented together with a detailed description of the needed geometrical framework. Given that in many cases gateways are placed according to regular patterns, Sec. VI applies the findings of the previous sections to compute the theoretical performance in terms of expected throughput of LoRaWAN multi-gateway deployments where gateways are placed according to square or equilateral triangular lattices. Finally, Sec. VII concludes the paper and envisages future works.

II. RELATED WORKS

The ALOHA multi-access protocol was proposed and first analyzed in the early seventies for a single-channel [START_REF] Abramson | The aloha system: Another alternative for computer communications[END_REF], [START_REF] Roberts | Aloha packet system with and without slots and capture[END_REF] and later for multiple channels [START_REF] Zhang | Multichannel aloha data networks for personal communications services (pcs)[END_REF]. The basic protocol underwent improvements to meet the requirements of applications and network technologies. For satellite networks, reservation [START_REF] Crowther | A system for broadcast communication: Reservation-aloha[END_REF], [START_REF] Lam | Packet broadcast networks-a performance analysis of the r-aloha protocol[END_REF], diversity [START_REF] Choudhury | Diversity aloha -a random access scheme for satellite communications[END_REF], and contention resolution with diversity [START_REF] Casini | Contention resolution diversity slotted aloha (crdsa): An enhanced random access schemefor satellite access packet networks[END_REF], [START_REF] Liva | Graph-based analysis and optimization of contention resolution diversity slotted aloha[END_REF] were added. In [START_REF] Bajović | Slotted aloha for networked base stations[END_REF], slotted ALOHA is analysed for a network of base stations which can cooperatively or noncooperatively decode signals from the end devices.

Other variants include multichannel ALOHA [START_REF] Zhang | Multichannel aloha data networks for personal communications services (pcs)[END_REF], ALOHA with CDMA [START_REF] Makrakis | Spread slotted aloha techniques for mobile and personal satellite communication systems[END_REF], with OFDMA [START_REF] Choi | Multichannel random access in OFDMA wireless networks[END_REF], coded slotted ALOHA [START_REF] Paolini | High throughput random access via codes on graphs: Coded slotted aloha[END_REF], [START_REF] Choi | Throughput analysis for coded multichannel aloha random access[END_REF]. We note that multichannel slotted ALOHA has applications in Machine Type Communications (MTC) [START_REF] Arouk | Multi-channel slotted aloha optimization for machine-type-communication[END_REF] which is related to the Internet of Things.

ALOHA for LoRaWAN is quite specific in that it has a duty cycle limitation which restricts the rate at which an end device can inject packets into the networks. A duty cycle limitation of 1% implies that for one unit of on-air time, an end device has to remain silent for 99 units of time before it can transmit another packet. This is related to but not quite the same as backoff mechanisms [START_REF] Rivero-Angeles | Differentiated backoff strategies for prioritized random access delay in multiservice cellular networks[END_REF] which are only initiated when a transmission is unsuccessful. A backoff mechanism requires an end device to remain silent following a packet collision for a backoff period which is usually a random variable. For example, the backoff period in UMTS-LTE is a uniform distribution on a fixed window size whereas IEEE 802. [START_REF] Casini | Contention resolution diversity slotted aloha (crdsa): An enhanced random access schemefor satellite access packet networks[END_REF] doubles the window size after each collision [START_REF] Seo | Design and analysis of backoff algorithms for random access channels in umts-lte and ieee 802.16 systems[END_REF]. The duty cycle limitation in LoRaWAN, on the other hand, applies to all transmissions whether successful or not. After every transmission, an end device has to remain silent for the duration defined by the duty cycle.

LoRaWAN recommends a pure ALOHA based medium access though slotted versions have been proposed to improve throughput [START_REF] Polonelli | Slotted aloha overlay on lorawan -a distributed synchronization approach[END_REF]. Among the first works to study the scalability of LoRaWAN there is [START_REF] Mikhaylov | Analysis of capacity and scalability of the lora low power wide area network technology[END_REF]. A mathematical model of LoRaWAN in ACK mode with retransmissions as well as its performance analysis was done in [START_REF] Bankov | Mathematical model of lorawan channel access[END_REF] without capture and in [START_REF]Mathematical model of lorawan channel access with capture effect[END_REF] with capture. The impact of retransmissions on the network lifetime is analysed in [START_REF] Paul | A novel mathematical model to evaluate the impact of packet retransmissions in lorawan[END_REF]. A stochastic geometry based approach for the scalability analysis is followed in [START_REF] Błaszczyszyn | Analyzing lora long-range, lowpower, wide-area networks using stochastic geometry[END_REF] for a single gateway and in [START_REF] Aftab | Scalability analysis of multiple lora gateways using stochastic geometry[END_REF] for multiple gateways. These works do not model the duty cycle limitation which is an important aspect in LoRaWAN and do not focus on applications like radio localization, which require packets to be received by multiple gateways simultaneously. Other than analytical studies, there have also been simulation based approaches for the scalability of LoRaWAN [START_REF] Adelantado | Understanding the limits of lorawan[END_REF]- [START_REF] Magrin | Performance evaluation of lora networks in a smart city scenario[END_REF].

The model considered here is closely related to the ones in [START_REF] Accettura | Low power radiolocation through long range wide area networks: A performance study[END_REF], [START_REF] Accettura | Simulating scalable long range wide area networks for very low power monitoring applications[END_REF]. In [START_REF] Accettura | Low power radiolocation through long range wide area networks: A performance study[END_REF], the focus was in finding the rate of successful transmissions without duty cycle limitation to at least 3 gateways in order to permit radio localization. There, a specific geometry for the gateway deployment was used, i.e., a honeycomb-like deployment. A small step ahead was done in [START_REF] Accettura | Simulating scalable long range wide area networks for very low power monitoring applications[END_REF], where the formulation for the throughput (i.e., rate of successful transmissions to at least 1 gateway) in the same honeycomb geometry was introduced. In this paper, we generalize the computation of the rate of successful transmissions to at least L gateways, where L can be any number, and all gateways are deployed at arbitrary locations that permit the coverage of any point on the deployment by at least L gateways. This generic formulation is then simplified for any regular deployment of gateways, i.e., according to equilateral triangular and square lattices. As matter of facts, the 2 formulas of [START_REF] Accettura | Low power radiolocation through long range wide area networks: A performance study[END_REF], [START_REF] Accettura | Simulating scalable long range wide area networks for very low power monitoring applications[END_REF] are just a few sub-models among the ones targeting equilateral triangular deployment of gateways. Another significant theoretical advancement presented in this contribution is a detailed modelling of duty cycle limitations on the frame transmission.

III. LORAWAN NETWORK ARCHITECTURE

A LoRa Wide Area Network [START_REF] Lora | 0.3 Specification[END_REF] structurally connects a number of very low power end devices to the Internet through wireless communications operating on distances of the order of kilometers. This is done by using a Chirp Spread Spectrum modulation, that enables the simultaneous transmission of several frames on the same bandwidth by using chirps, i.e., linear frequency sweeps. Indeed, a Spreading Factor (SF) defines the number of bits used to encode a symbol, i.e. the basic information carried by a chirp. A bigger SF corresponds to a bigger chirp time-on-air (ToA), that translates into a lower data rate and an increased transmission range. It can be quickly recognized that bigger SFs increase the frame sizes and the number of end devices accessing the same gateway, thus allowing very few chances for a frame to be transmitted without interference. Hence, the analysis performed throughout this paper just uses the smallest SF, since the expected corresponding capacity is the highest. In any case, similar reasoning can be replicated with higher SFs.

LoRaWAN end devices usually draw energy from batteries, whose lifetime is supposed to last around a decade. Such an energy saving requirement puts a constraint on the management of the most energy-expensive activity for end devices, i.e., on the wireless communication. In other words, the time for which the radio stays turned on for communication purposes should be kept as small as possible, while the efficiency in delivering any frame on the air should be taken as high as possible, i.e., frame collisions should be reduced and/or avoided. At the same time, the typical traffic handled by LoRaWANs is related to monitoring data flowing uplink from end devices towards any configured server throughout the Internet. In details, the applications targeted by LoRaWAN enable each end device to generate very low amount of frames sized in the order of tens to hundreds of bytes transmitted few times per day [START_REF] Thubert | Low-power wide-area networks at the IETf[END_REF].

Frames can be sent by end devices on any channel among a pool of available ones. The LoRaWAN specification for regional parameters [START_REF] Lora | RP002-1.0.1 LoRaWAN R Regional Parameters[END_REF] specifies that 3 channels must be always implemented by all end devices in the EU 863 -870 bandwidth for Europe. Any frame transmitted by a LoRaWAN end device is captured by all the LoRaWAN gateways in the coverage range of the end device itself. All captured replicas of the same frame by the gateways are forwarded over a reliable network to a LoRaWAN server. Such a server is in charge of discarding redundant copies and, if the frame requires a confirmation for its reception, it selects the gateway that will acknowledge the correct reception of the frame to the end device during a scheduled reception window. Therefore, the major energy consumption for wireless communications is confined on gateways, that are typically mains-powered and continuously listen on all the default channels to incoming LoRaWAN frames sent by end devices. As a major energy saving achievement, end devices preserve their own battery lifetime by consuming energy mostly for transmitting uplink.

The default wireless access in such network architecture is pure ALOHA, since any frame generated by the application running on any end device is immediately transmitted. Clearly, the low per-device data generation rate allows the pure ALOHA available capacity [START_REF] Abramson | The aloha system: Another alternative for computer communications[END_REF] to be shared by a higher number of end devices. In that, the objective of this contribution is to picture the effective capacity of a multigateway LoRaWAN, with the goal of providing a method able to identify how many end devices can be handled while insuring some minimal perceived Quality of Service. At the same time, an additional effort has been put to model possible duty cycle limitations [START_REF] Lora | RP002-1.0.1 LoRaWAN R Regional Parameters[END_REF] enforced by local radio access rules.

IV. SINGLE GATEWAY DEPLOYMENTS

The LoRaWAN framework holds some peculiar architectural features, whose joint performance has never been investigated from a detailed probabilistic point of view. Indeed, LoRaWAN enhances the ALOHA scheme by allowing the simultaneous exploitation of several channels at the same time. In addition, a frame transmitted by an end device can be received by several gateways: if at least one of those gateways is able to correctly receive the frame, the transmission is successful. A transmission is unsuccessful if all gateways experiment collisions. Intuitively, it is easy to figure out that the throughput is higher with respect to the case of radio coverage by a single end device.

One major difference with the standard ALOHA protocol is that LoRaWAN allows the regulatory authorities to impose a duty cycle (DC) limitation. A node that has just finished transmission has to wait a certain fixed amount of time before it can transmit another message. For example, the European Telecommunications Standards Institute (ETSI) has imposed a 1% duty cycle. One of the contributions of this work is the analysis of the ALOHA protocol with DC limitation. As mentioned in Sec. II, a very interesting work has investigated the throughput of Slotted ALOHA networks with coordinated gateways [START_REF] Bajović | Slotted aloha for networked base stations[END_REF]. The authors propose a model of the throughput based on the assumption that gateways are randomly spread on a given area according to a Poisson point process and that end devices can fall outside the coverage range of gateways. The throughput formulation is not straightforward and can be evaluated through approximations and numerical methods.

Contrariwise, in this paper it has been assumed that the gateways are non-randomly spread according to the following principle: wherever an end device is deployed, there is at least one gateway in its transmission range. This assumption comes from the common sense observation that realistic cellular networks deployments aim to get a complete coverage of an area. To the aim of evaluating scalable deployments, end devices are spread on the plane according to a Poisson point process.

Furthermore, it has been assumed that also the traffic generation process follows a Poisson process. In details, defining τ as the time-on-air of a LoRaWAN frame sent by a given end device to the closest gateways, any time interval will be expressed in units of τ throughout this contribution. The interarrival time between two consecutive frame generations will be assumed to follow an exponential distribution with average frame generation rate λ. Assuming a time interval t (expressed in units of τ ), the probability that the number N pkt (t) of frames generated in t equals a value k is:

P(N pkt (t) = k) = e -λt (λt) k k! . (1)
When a DC limitation applies, a single frame is transmitted at the beginning of a time interval long exactly ε units of time. This also means that, at the end of the transmission of a frame, the sending device shall remain silent for a duration of ε -1 units of time. As an example, a DC limitation of 1% will mean that ε = 100.

In the considered scenario, end devices are assumed to have no buffer. Any frame will be assumed to be lost if generated at an end device while the latter is either transmitting or staying silent due to a DC limitation. Frames that are not lost are transmitted as soon as they are generated. Let n be the number of channels. The end device selects a channel uniformly at random for transmitting a frame. Fig. 1 shows the evolution of the state of an end device. At time 0, a frame is transmitted (the gray rectangle), after which the end device has to wait for another ε -1 units before it can accept another frame for transmission. Since the inter-arrival times are exponentially distributed, this end device will wait for a time randomly picked according to an exponential distribution of rate λ before beginning a new cycle. The hatched rectangle is not needed for the illustrating the states of an end device. It will be needed later for the computation of the probability of success of a frame.

An end device is transmitting when it is in the gray area in Fig. 1, blocked during the ε -1 units when it cannot accept new frames for transmissions due to the duty cycle limitation. Otherwise, the end device is said to be idle. A cycle is thus a transmission state followed by a blocked state followed by an idle state. This cycle repeats itself ad infinitum.

Before computing the throughput for a multi-cell Lo-RaWAN network with Poisson number of end devices, it is instructive to compute the throughput for the simpler scenario of a single cell.

A. Performance measures: fixed number of nodes

First, we consider the case of a single cell with a fixed number of nodes. Theorem 1. A single cell LoRaWAN network operating over n frequency channels is assumed to be handling N end devices, each generating traffic according to an exponential distribution with average frame generation rate λ. Assume also that a duty cycle limitation applies, and after the transmission of a frame for 1 unit of time (i.e., the frame time-on-air), an end device remains idle for further ε -1 units of time. The throughput T (N ) is:

T (N ) = N gq N -1 .
(

) 2
where the average transmission rate g and the probability q that an end device does not interfere on the same channel with an ongoing transmission are respectively given by:

g = λ 1 + λε , (3)
and

q = 1 - 1 n λ min(ε, 2) + 1 -e λ min(ε-2,0) 1 + λε . (4)
Proof. Given that the overall network generated traffic is N λ, the throughput is related to the probability P succ of success as:

T (N ) = N λP succ .
To compute the probability of success of a frame, (called the tagged frame), the simultaneous occurrence of two independent events has to be considered: (a) the end device to which it arrives is idle; (b) the other nodes should not be transmitting and should not begin transmitting on the same channel where the tagged frame is being transmitted. In other words, P succ = P(a) • P(b).

With reference to the cycle pictured in Fig. 2a, the probability of event (a) can be computed from the Renewal Reward Theorem:

P(a) = λ -1 λ -1 + ε . 1 (TX) δ = 1 ε λ -1 Cycle Time (a) ε 2.
1 (TX) As matter of fact, a frame generated by an end device is transmitted with probability P(a), so that the average transmission rate g of an end device is obtained by multiplying the average frame generation rate λ by P(a). From this consideration, (3) follows.

δ = ε -1 ε λ -1 Cycle Time (b) ε < 2.
At the same time, the occurrence of event (b) is equivalent to saying that other nodes should not interfere during the transmission of the tagged frame. In details, the probability q that a single end device does not interfere with the transmission will be computed. Since the states of the end devices evolve independently, the probability of event (b) can be computed as the product of the individual probabilities q of the N -1 end devices, i.e., q N -1 . Hence, the throughput formulation of (2) is obtained.

The last step of this proof is related to the derivation of q. Consider an end device other than the one to which the tagged frame arrives. When the tagged frame arrives, this end device can be in 4 possible states (see Fig. 2a): (b.i) the end device is transmitting (TX); (b.ii) the end device is blocked but not in the hatched area; (b.iii) the end device is in the hatched area for a duration δ = 1; and (b.iv) the end device is idle for an average duration of λ -1 . It is worth noticing that such a general case refers to a DC limitation smaller than 50%, i.e., 2. For the sake of generality, bigger values of DC must be considered, i.e., 1 ε < 2. An example of such a case is pictured in Fig. 2b. Therein, an end device can be in one among the states (b.i), (b.iii), and (b.iv). Indeed, the state (b.ii) disappears, with the end device jumping to state (b.iii) at the end of state (b.i). Remarkably, the duration of (b.iii) is δ = ε -1. A more general formulation for the duration of state (b.iii) is δ = min(ε -1, 1).

In order to compute the probability q, it is worth defining some helping probability values. First, the probability that an end device selects the same channel adopted for the tagged frame delivery is 1/n. Define E I the event that an end device is either already transmitting when the tagged frame delivery starts or begins a transmission while the tagged frame is being transmitted. The probability q is then equal to:

q = 1 - 1 n P(E I ).
Herein, the probability of occurrence of E I can be expressed by the law of total probability:

P(E I ) = P(E I | b.i) P(b.i) + P(E I | b.ii) P(b.ii)+ + P(E I | b.iii) P(b.iii) + P(E I | b.iv) P(b.iv)
The probability of the events (b.i), (b.ii), (b.iii), and (b.iv) can be computed according to the Renewal Reward Theorem:

P(b.i) = 1 λ -1 + ε P(b.ii) = ε -1 -δ λ -1 + ε P(b.iii) = δ λ -1 + ε P(b.iv) = λ -1 λ -1 + ε
When an end device is in state (b.i), it is already transmitting, so that P(E I | b.i) = 1. Instead, when the end device is in state (b.ii), interference cannot happen, hence P(E I | b.ii) = 0. It is also very easy to deduce that if the end device is in state (b.iv) when the tagged frame delivery starts, then a possible interference can happen if at least one frame is generated on the considered end device during the tagged frame transmission, i.e., P(E I | b.iv) = 1 -e -λ . Finally, note that an end device in state (b.iii) will become idle during the transmission of the tagged frame: it could potentially need to transmit a generated frame which may collide with the tagged frame. This consideration leads to:

P(E I | b.iii) = 1 δ δ 0 1 -e -λ(x+1-δ) dx = = 1 - e -λ(1-δ) λδ + e -λ λδ .
Arranging together the computed probability values, the result of (4) is obtained.

Corollary 1. When a duty cycle limitation lower than 50% applies, the throughput over a single cell LoRaWAN network with N concurrent end devices is:

T LoRa (N ) = N λ 1 + λε 1 - 1 n • 2λ 1 + λε N -1 .
(5)

Proof. A DC 50% is equivalent to set ε 2 in (2), (3)
, and (4). Hence, (5) follows.

Corollary 2. The throughput over a single cell LoRaWAN network with N concurrent end devices and no duty cycle limitation is:

T P ure (N ) = N λ 1 + λ 1 - 1 n 1 - e -λ 1 + λ N -1 . (6)
Proof. The lack of DC limitation is equivalent to set ε = 1 in (2), (3), and (4). Hence, (6) follows.

Lemma 1. Assume a network scenario as the one pictured in Theorem 1. Then, the number N max of end devices allowing the throughput of (2) to achieve its maximum value is:

N max = 1 1 -q . ( 7
)
Proof. Such a result can be achieved by means of the ratio method [START_REF] Dahiya | An improved method of estimating an integer-parameter by maximum likelihood[END_REF] applied to (2).

Theorem 2. A duty cycle limitation lower than 50% (i.e., ε 2) allows the coexistence of a higher number of end devices with respect to a Pure ALOHA network with no duty cycle limitation. In other words:

N LoRa max N P ure max (8)
Proof. Thanks to Lemma 1, the number of end devices letting achieve the maximum throughput in a duty cycled LoRaWAN network is:

N LoRa max = n • 1 + λε 2λ , (9)
while in Pure ALOHA network such number would be:

N P ure max = n • 1 + λ 1 + λ -e -λ . (10)
Without loss of generality, proving ( 8) is equivalent to show that the argument of the floor function of ( 9) is major than the argument of the floor function of [START_REF] Abramson | The aloha system: Another alternative for computer communications[END_REF]. In other words, the following inequality must always hold true for any value of λ and ε 2:

1 + λε 2λ 1 + λ 1 + λ -e -λ . (11)
The previous inequality can be rearranged as:

ε -2 2e -λ 1 + λ -e -λ - 1 λ ,
where the first term must be greater or equal to 0. If the second term is lower or equal than 0, then the theorem is proved.

Rearranging again, it has to be proved that:

(1 + λ)e λ 1 + 2λ.

Since e λ 1 + λ for any value of λ, then:

(1 + λ)e λ (1 + λ) 2 = 1 + 2λ + λ 2 > 1 + 2λ,
thus proving the theorem.

The last Theorem shows that for a number of devices N LoRa max , the achievable throughput with a duty cycle limitation corresponding to ε 2 is maximum. In order to find the network scenarios where a duty cycle limitation helps in 12) and [START_REF] Crowther | A system for broadcast communication: Reservation-aloha[END_REF] better handling the available bandwidth, it must hold true the following inequality:

T LoRa (N LoRa max ) T P ure (N LoRa max ), or, in explicit terms and after a proper rearrangement:

n + λ(nε -2) (n -1) + λ(n -1) + e -λ 1- 1 n(1+λε) 2λ 1 + λε 1 + λ (12)
Given that x x -1 for any real number x, it is easy to recognize that if the following inequality holds true:

n + λ(nε -2) (n -1) + λ(n -1) + e -λ 1- 2λ n+λ(nε-2) 1 + λε 1 + λ , (13)
then [START_REF] Zhang | Multichannel aloha data networks for personal communications services (pcs)[END_REF] holds true as well.

The inequalities [START_REF] Zhang | Multichannel aloha data networks for personal communications services (pcs)[END_REF] and ( 13) depend on 3 parameters: the generation rate λ, the duty cycle parameter ε, and the number n of channels. Fixing the number of channels to n = 1, the region on the (λ, ε) plane where the 2 aforementioned inequalities are verified is plot in Fig. 3a. In details, consider the implicit function resulting by focusing on just the equality embedded in [START_REF] Zhang | Multichannel aloha data networks for personal communications services (pcs)[END_REF]. The sawtooth-shaped boundary of the light gray area corresponds to the locus of points (λ, ε) satisfying such an implicit function. As it can be seen, the non-linear behavior introduced by the floor function in the exponent of the leftmost side of ( 12) is reflected by the oscillatory trend of the aforementioned boundary. For the sake of an easy mathematical tractability, the relaxed inequality of ( 13) can be considered as the wanted design constraint: a duty cycle limitation can be beneficial for a range of traffic generation rates, and the suitability can be verified by checking the verification of (13), i.e., if a given (λ, ε) setting falls into the dark gray zone of Fig. 3a. For completeness, Fig. 3b shows similar representation when the the number of channels is fixed to 5.

As an example, a maximum size LoRaWAN frame on SF 7 is delivered for a time-on-air of τ = 368.896 ms [START_REF] Accettura | Low power radiolocation through long range wide area networks: A performance study[END_REF]. If a monitoring application requires the delivery of 1 LoRaWAN frame with maximum size per hour (i.e., 3600s), the corresponding frame generation rate expressed in τ units is slightly higher than λ = 10 -4 . By making a qualitative inspection of Fig. 3a, and 3b, it can be quickly verified that a duty cycle limitation lower than 1% is always beneficial with respect to the case when no duty cycle policy is adopted.

We end this subsection by looking at what happens when to [START_REF] Zhang | Multichannel aloha data networks for personal communications services (pcs)[END_REF] when the number of channels is large, that is when

n → ∞. Let O n = {(λ, ) : (12) is satisfied}. When n → ∞,
O n converges to the region given by the following proposition.

Proposition 1. Consider the inequality

1 + λ -e -λ 2λ • 1 + λε 1 + λ -1 log 1 + λε 1 + λ . ( 14
)
Then,

lim n→∞ O n = {(λ, ) : (14) is satisfied}, (15)
Proof. The proof follows by taking logarithm on both sides of (12) and writing the series expansion of both sides in n -k . The zeroth order term cancels out and for the inequality to hold in the limit, the first order coefficients have to be ordered the same way as [START_REF] Zhang | Multichannel aloha data networks for personal communications services (pcs)[END_REF]. Formally, taking the logarithm of both sides of ( 12) and focusing on the LHS in which we have ignored the floor function, we expand in powers of n -1 and get

1 - 2λ n(1 + λε) log n + λ(nε -2) (n -1) + λ(n -1) + e -λ = 1 - 2λ n(1 + λε) log 1 + λε 1 + λ 1 -2λ n(1+λε) 1 -1+λ-e -λ n(1+λ) = log 1 + λε 1 + λ + 1 n -2λ (1 + λε) log 1 + λε 1 + λ - 2λ (1 + λε) + 1 + λ -e -λ (1 + λ) + o(n -1 ).
Note that the coefficient of n 0 cancels out with the RHS. Further, for n → ∞, the n -2 terms become negligible and Ω Gateway Θ for [START_REF] Zhang | Multichannel aloha data networks for personal communications services (pcs)[END_REF] to hold, the coefficient of n -1 on both sides of [START_REF] Zhang | Multichannel aloha data networks for personal communications services (pcs)[END_REF] should also satisfy the same inequality. This gives

- 2λ (1 + λε) log 1 + λε 1 + λ - 2λ (1 + λε) + 1 + λ -e -λ (1 + λ) 0. ( 16
)
Rearranging the above inequality, we arrive at [START_REF] Lam | Packet broadcast networks-a performance analysis of the r-aloha protocol[END_REF].

It will seen in the next subsection that when the number of nodes is Poisson distributed, we surprisingly arrive at exactly the same inequality (see [START_REF] Mikhaylov | Analysis of capacity and scalability of the lora low power wide area network technology[END_REF]) as ( 14) irrespective of the number of channels.

B. Performance measures: Poisson number of nodes

With reference to Fig. 4, denote by Ω the compact set of points covered by the gateway of the cell, and assume that end devices are spread over Ω according to a Poisson Point Process of density µ. Then, the probability that the number N ed (Θ) of end devices lying on any compact subset Θ ⊆ Ω equals i is:

P(N ed (Θ) = i) = e -µA(Θ) [µA(Θ)] i i! , (17)
where A(Θ) is the measure of the area of Θ. The throughput of the cell Ω can be computed for this scenario by applying the result of the following Proposition.

Proposition 2. Given a random number of end devices deployed on the coverage area Ω around a gateway according to a Poisson Point Process of density µ, and that the throughput for a fixed number of end devices is computed according to (2), then the expected throughput S(Ω) is

S(Ω) = gµA(Ω) • e -(1-q)µA(Ω) . ( 18
)
Proof. By the law of total probability,

S(Ω) = ∞ i=1 T (i) P(N ed (Ω) = i) = = ∞ i=1 igq i-1 e -µA(Ω) [µA(Θ)] i i! = = ge -µA(Ω) µA(Θ) ∞ i=1 [qµA(Θ)] i-1 (i -1)! = = gµA(Ω) • e -(1-q)µA(Ω) .
As an example, when the coverage area Ω around the gateway is a disk of radius 1, so that A(Ω) = π, the throughput is expressed as follows:

S(Ω) = gµπ • e -(1-q)µπ . ( 19
)
We take a look at the set of parameters for which it is beneficial to have a duty cycle limitation. For n = 1 and A = 1, Figure 5 shows the region in the (λ, µ) plane where LoRaWAN has a better throughput than ALOHA for ε = 2 (Fig. 5a) and for ε = 100 (Fig. 5b).

Intuitively, for small frame generation rates, that is λ → 0, there will be very few collisions in ALOHA unless the the node density, µ, is sufficiently large to generate enough on-air traffic so as to adversely impact the throughput. Thus, for small λ, ALOHA is expected to give a better throughput for small µ whereas LoRaWAN would be better for larger µ. On the other hand, when λ → ∞, the duty cycle limitation is much more beneficial when ε > 2 since this leaves enough space an end devices to get a message through during the silence periods of of the others. When ε < 2, the silence periods are not long enough to allow a positive throughput in the limit, and both ALOHA and LoRaWAN have similar throughput. That is, ε = 2 defines a point of phase transition in the ε > 2.

The exact asymptotic behaviour of this boundary is stated in the following result. Theorem 3. Let μ = µA n be mean number of nodes per channel.

1) When λ → 0, for LoRaWAN to have a better throughput than ALOHA, we require

λμ ≥ 2 ε+1 , if ε ≤ 2; λμ ≥ 2( -1) 4 -5 , if ε > 2. ( 20
)
2) When λ → ∞, for LoRaWAN to have a better throughput than ALOHA, we require

e λ(ε-2) λ -1 μ ≥ ε log ε, if ε ≤ 2; μ ≥ ε ε-2 log ε, if ε > 2. ( 21
)
Proof. Starting from [START_REF] Bajović | Slotted aloha for networked base stations[END_REF] and doing some basic algebraic manipulations, for LoRaWAN to be better than ALOHA, we arrive at the following inequality.

µA(q Lora -q P ure ) ≥ log 1 + λε

1 + λ . ( 22
)
We give the proof for the case λ → ∞. The other case is similar.

Let ε > 2. For λ → ∞, from (4), observe that q P ure → 1 -1 n whereas q LoRa → 1 -2 nε . Substituting these values in [START_REF] Choi | Throughput analysis for coded multichannel aloha random access[END_REF], the subcase ε > 2 in (21) follows.

For ε ≤ 2, where the second line is a consequence of the fact that 1 < ε ≤ 2 which implies that e λ(ε-2) dominates e -λ .

q LoRa -q P ure ≈ 1 n e λ(ε-2) λε - e -λ λ , ≈ 1 n e λ(ε-2) λε ,
Theorem 3 given an indication of how the mean number of nodes per channel can scale with the arrival rate while keeping the benefits of duty cycle limits. Nest, we compare LoRaWAN and ALOHA for specific node densities which maximize the respective throughputs.

Lemma 2. Assume a network scenario as the one pictured in Proposition 2. Then, the density µ max of end devices allowing the throughput of (18) to achieve its maximum value is:

µ max = 1 A(Ω) • 1 1 -q . ( 23
)
Proof. The value of µ max is the one for which the derivative of ( 18) is equal to 0.

Lemma 3. A duty cycle limitation lower than 50% (i.e., ε 2) allows the coexistence of a higher density of end devices per area with respect to a Pure ALOHA network with no duty cycle limitation. In other words:

µ LoRa max µ P ure max (24)
Proof. Thanks to Lemma 2, the number of end devices letting achieve the maximum throughput in a duty cycled LoRaWAN network is:

µ LoRa max = 1 A(Ω) • n(1 + λε) 2λ , (25)
while in Pure ALOHA network such number would be:

µ P ure max = 1 A(Ω) • n(1 + λ) 1 + λ -e -λ . (26)
Plugging ( 25) and ( 26) into [START_REF] Rivero-Angeles | Differentiated backoff strategies for prioritized random access delay in multiservice cellular networks[END_REF] drives to solve the same inequality of [START_REF] Roberts | Aloha packet system with and without slots and capture[END_REF], thus this Lemma is proved by means of the same reasoning of Theorem 2.

The last lemma shows that for a density of devices µ LoRa max , the achievable throughput with a duty cycle limitation corresponding to ε 2 is maximum. In order to find the network scenarios where a duty cycle limitation helps in better handling the available bandwidth, it must hold true the following inequality:

S LoRa Ω, µ LoRa max S P ure Ω, µ LoRa max ,
where S LoRa Ω, µ LoRa max is obtained from (18) when ε 2, while setting ε = 1 in (18) leads to S P ure Ω, µ LoRa max . In explicit terms and after a proper rearrangement, it must hold true that

1 + λ -e -λ 2λ • 1 + λε 1 + λ -1 log 1 + λε 1 + λ . (27)
As it can be noticed, there is no dependence on the number n of channels. When a duty cycle limitation is enforced, [START_REF] Mikhaylov | Analysis of capacity and scalability of the lora low power wide area network technology[END_REF] provides a rule of thumb on the setting of the minimum and maximum data generation rate λ for each device. The resulting range of data rates better fits duty cycle limitations with respect to the case when no limitation is super-imposed. Remarkably, inequality [START_REF] Mikhaylov | Analysis of capacity and scalability of the lora low power wide area network technology[END_REF], which is derived for a Poisson distribution of nodes, turns out to be the same as [START_REF] Lam | Packet broadcast networks-a performance analysis of the r-aloha protocol[END_REF] which was obtained for large number of channels and deterministic number of nodes.

C. Model validation

The model of ( 18) has been validated through simulations. In details, the Python event-driven simulator LoRaWAN-sim has been used to mimic a single-gateway LoRaWAN network collecting information delivered by end devices in the gateway coverage range. All the configured settings reflect the regulatory policies [41] introduced by the European Telecommunications Standards Institute (ETSI) for the EU 863 -870 bandwidth [START_REF] Lora | RP002-1.0.1 LoRaWAN R Regional Parameters[END_REF].

For the sake of validation, both gateways and end devices are configured to use spreading factor SF 7, which allows the fastest datarate in the currently most used hardware platforms, and to work on a single channel in the pool of available (and default) 3 in the sub-band between 868.0MHz and 868.6MHz. Each end device generates LoRaWAN frames of maximum size on SF 7, that corresponds to a time-on-air of τ = 368.896ms [START_REF] Accettura | Low power radiolocation through long range wide area networks: A performance study[END_REF]. Each frame is randomly generated according to an exponential distribution featured by an average interarrival time of 60 seconds. In other words, the average frame generation rate is configured to be λ = 0.368896 60 ∼ 0.006148267. Frames do not require an acknowledgement, so that retransmissions are disabled.

Two scenarios are considered, on one hand the lack of any duty cycle limitation, on the other hand the enforcement of a 1% duty cycle, as stated by the ETSI regulatory specification for the considered sub-band. In facts, the first case shows a typical pure ALOHA network, while the second scenario refers to a realistic LoRaWAN network. Without loss of generality, it has been assumed that distances are measured as multiples of the transmission range. A unit area is the area corresponding to a square having side equal to the transmission range. In such a context and for each of the duty cycle scenarios considered above, the density of end devices per unit area is increased 5 by 5 starting from 0 until 80. Fixed the duty cycle option and the density of end devices, 20 different network sizes are randomly generated by varying the seed feeding the random number generator.

Each of the considered networks is evaluated for a simulated time of 1 day. The resulting estimated throughput is computed as the portion of time used by the gateway to receive LoRaWAN frames not corrupted by collisions. The 20 throughput values related to the same density of end devices and duty cycle option are then aggregated to produce an average estimated throughput. All these results are plot in Fig. [START_REF] Palattella | Internet of things in the 5g era: Enablers, architecture, and business models[END_REF].

In details, the dark red lines refer to the lack of any duty cycle policy, while the light red lines correspond to the 1% duty cycle expected in LoRaWAN networks. For any configured duty cycle option, there are 2 plot lines: the dotted lined comes from simulation results, and the resulting confidence interval defined by a 95% level is also plot for each evaluation point; instead, the continuous line is obtained by properly substituting the needed parameters in [START_REF] Bajović | Slotted aloha for networked base stations[END_REF]. As it can be quickly verified by inspection, the model of ( 18) is validated by simulation results.

D. Preliminary discussion for extending the model to multigateway deployments

Since the remaining part of this paper details the throughput computation in multi-gateway deployments, it is worth providing some preliminary result that will be widely exploited in that theoretical landscape. Proposition 3. Consider the whole compact set Ω falling under the coverage range of 1 gateway (or more). Assume that a random number of end devices is deployed on Ω according to a Poisson Point Process of density µ. Assume also that q is the probability that an end device is not transmitting on a given channel. Then, the probability Q(Ω) that none of the end devices ∈ Ω is transmitting on a given channel is:

Q(Ω) = e -(1-q)µA(Ω) . ( 28
)
Proof. If the number of end devices on a given area is i, then the joint probability that none of them is transmitting equals q i . By the law of total probability,

Q(Ω) = ∞ i=0 q i P(N ed (Ω) = i) = = ∞ i=0 q i e -µA(Ω) [µA(Ω)] i i! = = e -(1-q)µA(Ω) .
Suppose to be interested to the throughput of end devices located in Θ ⊆ Ω as illustrated in Fig. 4. Since the number of end devices is Poisson-distributed, the throughput S Ω (Θ) has the following decomposition. Proposition 4. Let Θ be a compact set falling under the coverage range of 1 gateway (or multiple gateways), and the superset Ω ⊇ Θ that results as the whole coverage area around the gateway (or the considered multiple gateways). Given a random number of end devices deployed on Ω according to a Poisson Point Process of density µ, and that the throughput for a fixed number of end devices is computed according to (2), then the expected throughput S Ω (Θ) related to Θ is

S Ω (Θ) = gµA(Θ) • Q(Ω). (29)
Proof. Consider the sets Θ and Ω \ Θ. Since these sets are disjoint, the number of nodes in each of them is an independent Poisson random variable with density µ. Conditioning on the number of nodes in each of these two subsets and then using the law of total probability, it results:

S Ω (Θ) = ∞ i=1 ∞ j=0 igq i-1+j • P(N ed (Θ) = i)• • P(N ed (Ω \ Θ) = j) = = ∞ i=1 igq i-1 e -µA(Θ) [µA(Θ)] i i! • • ∞ j=0 q j e -µA(Ω\Θ) [µA(Ω \ Θ)] j j! = = [gµA(Θ) • e -(1-q)µA(Θ) ] • [e -(1-q)µA(Ω\Θ) ] = = gµA(Θ) • e -(1-q)µA(Ω) = = gµA(Θ) • Q(Ω).
The decomposition property says that the throughput is a product of two terms: (i) gµA(Θ), which is the traffic generated from Θ, and (ii) e -(1-q))µA(Ω) , which is the probability that the other nodes in Ω are not interfering. The first term depends only on Θ, whereas the second term depends on the whole cell Ω.

V. MULTI-GATEWAY DEPLOYMENTS

Assuming that end devices are randomly spread on the Euclidean plane according to a homogeneous point Poisson process and that they transmit link-layer frames according to a Poisson process as well, the purpose of this section is to compute the average throughput of ALOHA multi-gateway networks, as the rate of successful receptions by at least L gateways. When L = 1, the searched formulation expresses the common meaning of throughput used in the literature until now, i.e., the rate of successful reception by at least one gateway. When L > 1, frames are expected to be received successfully by many gateways in order to enable redundancy for smart applications as geo-fencing, radiolocation, low power tracking, etc.

As matter of fact, the searched formulation is related to the portion of the Euclidean plane covered by at least L gateways, and is general enough to account for any compliant deployment of gateways. The arguments needed for such a computation are not much straightforward, but some interesting simplifying considerations can be made. Indeed, the first sub-problem to be overcome is the computation of the rate of successful transmissions from end devices positioned within a compact subset Θ of the Euclidean space R 2 to all the gateways in the subset Γ of the set of gateways Γ * = {γ 1 , γ 2 , . . . } deployed on R 2 as well (i.e., γ i ∈ R 2 ). In order to gently introduce the needed procedure, the definitions of some "interesting" compact subsets of the plane are presented. Definition 1. Given a device positioned at α ∈ R 2 , its associated coverage area is a disk centered at α with radius equal to R:

C(α) := {β ∈ R 2 | ||β -α|| R}. (30)
In other words, positioning a gateway on the plane is equivalent to draw a circle with radius R (equal to the transmission range) around the gateway itself: the resulting disk defines the portion of the plane where any deployed end device would be able to wirelessly interact with the gateway, i.e., the gateway coverage area. With similar arguments, an end device positioned on the plane defines a coverage area around itself: any gateway placed not further than R around the end device will be able to communicate with it. Definition 2. Given a set Γ of gateways positioned on the Euclidean space R 2 , an end device positioned at α ∈ R 2 can communicate with all the gateways ∈ Γ if it belongs to the intersection I Γ among all the disks centered at the considered gateways:

I Γ = γ∈Γ C(γ). (31)
Definition 3. Given a set Γ of gateways positioned on the Euclidean space R 2 , an end device positioned at α ∈ I Γ is able to communicate with all the gateways ∈ Γ and it can suffer from interfering communications started by end devices in the union U Γ among all the disks centered at the considered gateways:

U Γ = γ∈Γ C(γ). (32)
The last two definitions are used to setup a specific topological environment, that in turn eases the computation of the rate of successful transmissions from end devices ∈ Θ to all the gateways ∈ Γ. Indeed, to perform such a computation, it is required that Θ ⊆ I Γ , i.e., any end device positioned within Θ must be able to communicate with all gateways belonging to Γ. Here, frames delivered by end devices in Θ can be successfully received at the gateways if none of the other end devices in the same area interfere; in addition, to insure that frames are successfully received by all gateways ∈ Γ, all the end devices ∈ U Γ \ Θ must not interfere too. Summing up these considerations and applying Proposition 4 to such a multi-gateway deployment, the rate of successful transmissions from end devices ∈ Θ to all the gateways ∈ Γ is:

S Γ (Θ) = gµA(Θ) • Q(U Γ ). ( 33
)
The throughput quantity of ( 33) is shaped as the product of two terms: (i) the overall generated traffic, represented by gµA(Θ); and (ii) the probability of successful reception, that is defined, in the case under consideration, as the probability that none of the other end devices interferes and causes a collision on any of the gateways in Γ. However, ( 33) is subject to some specific and stringent conditions: a frame sent by an end device in Θ must be successfully received by all gateways in Γ; the compact subset on which the calculation is done, Θ, must be a subset of all coverage areas defined by gateways ∈ Γ. The target of this section is instead to picture the rate of successful transmissions S L (Ω) from end devices on any compact subset Ω of the plane to at least L gateways among the ones reachable by each end device ∈ Ω. To approach such a goal, it is important to figure out how the coverage areas around each deployed gateway partition the plane in smaller compact subsets like, for instance, E 3 and E 4 in Fig. 7. Indeed, the property of each of these compact subsets is that an end device positioned in it can communicate exclusively with all the gateways belonging to a "well-defined" subset of the ones deployed on the plane. To this end, some further definitions are presented in what follows.

Definition 4. Given the set Γ * of gateways deployed on the Euclidean space R 2 , the subset of gateways associated to a point α ∈ R 2 is defined as:

σ(α) := {γ ∈ Γ * | γ ∈ C(α)}. ( 34
)
As pictured in Fig. 7, the presence of multiple gateways partitions the surface into areas, each of them featured by the number of gateways that "cover" that same area. In the example of Fig. 7, any point ∈ E 3 (whose border is drawn with dashed blue lines) simultaneously falls into the coverage area of the gateways placed at P , Q, and S. The computation of the rate of successful transmissions by end devices deployed in E 3 suffers from the interference of end devices deployed on the union area bounded by the blue line. Instead, any end device positioned at a point ∈ E 4 (whose border is drawn with dashed red lines) can deliver a frame to the gateways positioned at P , Q, R, and S. The computation of the rate of successful transmissions by end devices deployed in E 4 suffers from the interference of end devices deployed on the union area bounded by a red line. These cases show that there are compact subsets of the Euclidean plane, where any point simultaneously falls into the coverage area of a specific subset of gateways. The following definition describes such an evidence. Definition 5. A subset Γ of the set Γ * of gateways deployed on the Euclidean space R 2 and a compact subset Θ ⊂ R 2 are exclusively associated if the following condition is valid:

∀α ∈ Θ : σ(α) = Γ. ( 35
)
As matter of facts, the compact subsets E 3 and E 4 in Fig. 7 have specific properties. E 3 is the compact subset that contains all the points that simultaneously belongs to all and only the coverage areas centered on P , Q, and S. Instead, E 4 is the compact subset that contains all the points that simultaneously belongs to all and only the coverage areas centered on P , Q, R, and S. The following definition helps in shedding some light on the previous finding. Definition 6. The compact subset ∆ Γ ⊂ R 2 defined as:

∆ Γ := {α ∈ R 2 | σ(α) = Γ}, ( 36
)
contains all points having Γ as associated set of gateways. This subset is assigned to Γ.

As a consequence, E 3 is assigned to {P, Q, S}, while E 4 is assigned to {P, Q, R, S}. The computation of the throughput on E 3 (and on any of its subsets) is not correlated with the computation of the throughput on E 4 (and on any of its subsets) and elsewhere, since E 3 and E 4 are assigned to different sets of gateways. It is then important to identify such separate areas: the next Lemma and Theorem prove that different subsets of gateways always determine separate assigned compact subsets on the Euclidean space R 2 . Lemma 4. A compact subset Θ ⊂ R 2 associated to a subset Γ of the set Γ * of gateways deployed on the Euclidean space R 2 , is a subset of the subset ∆ Γ assigned to Γ:

∀α ∈ Θ : σ(α) = Γ ⇐⇒ Θ ⊆ ∆ Γ . ( 37
)
Proof. Assuming that Θ is exclusively associated to ∆ Γ , the proof of Θ ⊆ ∆ Γ can be done by "reductio ad absurdum". As matter of facts, assuming that Θ is not subset of ∆ Γ is equivalent to state that Θ \ ∆ Γ = Θ 0 = ∅. Hence, ∃α ∈ Θ 0 so that α / ∈ ∆ Γ , thus implying that σ(α) = Γ according to Definition 6. This leads to a contradiction with the hypothesis, so that the implication is true, because it can not be false.

Instead, proving the opposite implication is straightforward: given that any point α ∈ ∆ Γ has got the property that σ(α) = Γ, any subset Θ ⊆ ∆ Γ contains points having the same property.

Γ = B =⇒ ∆ Γ ∩ ∆ B = ∅. ( 38
)
Proof. Assuming by "reductio ad absurdum" that ∆ Γ ∩ ∆ B = Θ 0 = ∅, then, according to Lemma 4, both the following implications hold true:

• Θ 0 ⊆ ∆ Γ , so that ∀α ∈ Θ 0 : σ(α) = Γ; • Θ 0 ⊆ ∆ B , so that ∀α ∈ Θ 0 : σ(α) = B.
In other words, ∀α ∈ Θ 0 ⊂ R 2 : σ(α) = Γ = B, which contradicts with the assumption Γ = B. So, [START_REF] Thubert | Low-power wide-area networks at the IETf[END_REF] is true, because it can not be false.

Given a set of gateways Γ, a throughput formulation involving all and only the gateways ∈ Γ can be computed only on the assigned compact subset ∆ Γ (and its subsets too). The following proposition provides a method for identifying the assigned compact subset ∆ Γ .

Proposition 5. Given a subset Γ of the set Γ * of gateways deployed on the Euclidean space R 2 , its assigned compact subset ∆ Γ ⊂ R 2 is expressed by:

∆ Γ = I Γ \ U Γ C , (39)
where

Γ C = Γ * \ Γ.
Proof. According to Definition 6, the subset ∆ Γ contains all the points whose associated subset of gateways is Γ. Therefore, any point α ∈ ∆ Γ belongs simultaneously to the coverage area C(γ) around each gateway γ ∈ Γ, so it must be true that α ∈ I Γ .

At the same time, α must not fall into the coverage areas centered at other gateways on the plane, i.e, α / ∈ C(γ), with γ ∈ Γ * \Γ = Γ C . In other words, it must be true that α / ∈ U Γ C . All in all, in order to have σ(α) = Γ, it must be true that α ∈ I Γ \ U Γ C , from which [START_REF] Lora | RP002-1.0.1 LoRaWAN R Regional Parameters[END_REF] follows.

The rate of successful transmissions S L (Ω) from end devices on any compact subset Ω of the plane to at least L gateways is formally derived later in Theorem 5. For the sake of clearness, two preliminary needed lemmas are separately presented in what follows.

Lemma 5. Given two integer values F and L, with L F , then:

F l=L F l l -1 L -1 (-1) l-L = 1 (40)
Proof. Using Theorem 1 in [START_REF] Call | Pascal's matrices[END_REF], it is straightforward to prove that for any couple of integer values F and L, with L F , the following equation is true:

F l=L F l l L (-1) l-L = 0. (41)
Fixing the values F and L, it is by now assumed that the LHS expression in ( 40) is equal to a value x. Subtracting the LHS expression in (41) from the LHS expression in [START_REF] Dahiya | An improved method of estimating an integer-parameter by maximum likelihood[END_REF], it is found:

F l=L F l l -1 L -1 - l L (-1) l-L = = F l=L F l l -1 L (-1) l-L-1 = = F l=L+1 F l l -1 L (-1) l-L-1 = L =L+1 = F l=L F l l -1 L -1 (-1) l-L = x. (42)
Recursively, the leftmost expression in ( 40) is then equal to x for any couple of values L and F . This is true also when L = F , so that x is found to be equal to 1, and the Lemma is proved.

Lemma 6. Given a set Γ and a positive number L |Γ|, the system of linear equations defined by:

Ψ⊆Φ |Ψ| L c Ψ,L = 1, ∀Φ ⊆ Γ, |Φ| L, (43)
has got a unique solution:

c Φ,L = (-1) |Φ|-L |Φ| -1 L -1 , ∀Φ ⊆ Γ, |Φ| L. (44)
Proof. It is first proved that each value c Φ,L is only dependent on the the size of the subset Φ. Indeed, the leftmost term in (43) can be expanded as:

Ψ⊆Φ |Ψ| L c Ψ,L = |Φ| l=L Ψ⊆Φ |Ψ|=l c Ψ,L = c Φ,L + |Φ|-1 l=L Ψ⊂Φ |Ψ|=l c Ψ,L ,
so that, any variable c Φ,L associated to a set Φ can be expressed as function of the variables associated to all the subsets of Φ:

c Φ,L = 1 - |Φ|-1 l=L Ψ⊂Φ |Ψ|=l c Ψ,L . (45)
By means of the induction principle, the statement c Φ,L = c |Φ|,L is now proved, i.e., the variable c Φ,L associated to a set Φ is only dependent on its size. The base case in the induction is represented by any subset Φ ⊂ Γ with minimal acceptable size, i.e., |Φ| = L. The resulting coefficients are:

c Φ,L = 1 = c L,L ∀Φ ⊆ Γ, |Φ| = L, (46)
c Φ,L = 1 - F -1 l=L Ψ⊂Φ |Ψ|=l c Ψ,L = 1 - F -1 l=L c l,L Ψ⊂Φ |Ψ|=l 1 = = 1 - F -1 l=L F l c l,L = c F,L , (47)
hence, the coefficient depends only on the set size, and the assertion is proved.

The previous result is used to simplify the system of (43) into the following one:

F l=L F l c l,L = 1, ∀F L ∧ F |Γ|. (48)
By using Lemma 5, such a system of equations is solved when the coefficients are shaped as in (44).

Theorem 5. Given a compact subset Ω ⊂ R 2 , if the size of the set of gateways σ(α) associated to any point α ∈ Ω is greater or equal to L, it is possible to compute the rate S L (Ω) of successful transmissions from end devices ∈ Ω to at least L gateways:

S L (Ω) = gµ Γ⊂Γ * |Γ| L ∆Γ∩Ω =∅ A(∆ Γ ∩ Ω) • P Γ L , with ∀α ∈ Ω : |σ(α)| L, (49)
where

P Γ L = |Γ| l=L (-1) l-L l -1 L -1 Φ⊆Γ |Φ|=l Q(U Φ ) (50)
Proof. By using the result of Theorem 4, a deployment of gateways Γ * partitions the Euclidean plane in non overlapping compact subsets, each of them assigned to a specific subset Γ ⊂ Γ * . Therein, Ω is also partitioned by Γ * : given a subset

Γ ⊂ Γ * with |Γ| L and ∆ Γ = ∅, if the intersection Θ = Ω ∩ ∆ Γ = ∅,
then Θ is one of the partitions created by Γ * on Ω. According to Lemma 4, Θ and Γ are exclusively associated. Hence, computing the rate of successful transmissions throughout Ω is equivalent to sum up the rate of successful transmissions computed on each partition created by the deployment of gateways:

S L (Ω) = Γ⊂Γ * |Γ| L ∆Γ∩Ω =∅ S L (∆ Γ ∩ Ω), with ∀α ∈ Ω : |σ(α)| L. (51)
In the previous formula, each term S L (Θ) in the summation is the rate of successful transmissions from end devices ∈ Θ = ∆ Γ ∩ Ω, and can be expressed as:

S L (Θ) = gµA(Θ) • P Γ L , with ∀α ∈ Θ : σ(α) = Γ. (52)
Similarly to [START_REF] Adelantado | Understanding the limits of lorawan[END_REF], S L (Θ) in ( 52) is the product of (i) the overall generated traffic, represented by gµA(Θ), and (ii) the probability P Γ L of successful reception by at least L gateways in Γ, i.e., the probability that no collision happens on at least L gateways. Combining (51) and (52), the first part of the theorem, i.e., equation (49), is proved. Herein, the formulation of ( 50) is proved too.

To this purpose, it is worth remarking that the value Q(U Φ ), defined by [START_REF] Bankov | Mathematical model of lorawan channel access[END_REF], is associated to the union of the disks centered at the gateways ∈ Φ and represents the probability that none of those gateways detect collisions. In general, given a set Γ of gateways, with ∆ Γ = ∅, the probability of successful reception by at least L |Γ| gateways measures how likely is to occur that all the gateways in any subset Φ ⊆ Γ with |Φ| L do not detect collisions. As a consequence, such a probability is a linear combination of the values Q(U Φ ), with L |Φ| |Γ|:

P Γ L = Φ⊆Γ |Φ| L c Φ,L Q(U Φ ), (53)
and the computation of all c Φ,L coefficients is the only missing building block needed for drawing the final formulation.

It is very important to notice that the probability value Q(U Φ ), with Φ ⊆ Γ, is computed regardless of whether gateways ∈ Γ \ Φ detect collisions or not. Interestingly, consider the simultaneous occurrence of the following conditions: (i) all gateways ∈ Φ do not detect collisions, and (ii) all gateways ∈ Γ\Φ detect collisions. Such a case is not only accounted by Q(U Φ ), but also by the probability Q(U Ψ ) associated to any subset Ψ of Φ. Clearly, it must be taken into account just once in the linear combination under analysis. This consideration results in the set of

|Γ| l=L |Γ| l
constraints of (43). Thus, the searched solution can be found by using Lemma 5 and Lemma 6, and the resulting coefficients can be plugged into (53) to obtain (50).

For the sake of comparison with the case of single gateway presented in [START_REF] Makrakis | Spread slotted aloha techniques for mobile and personal satellite communication systems[END_REF], it is worth considering an average value of the aforementioned quantity that relates to any compact subset of the plane Ω 0 featured by an area A(Ω 0 ) = π. In facts, consider any compact subset of the plane Ω, whose measure is A(Ω) π, and that is featured by the property that any point α ∈ Ω falls under the coverage area of at least L gateways. Then, the average rate of successful transmissions S Ω L from end devices deployed on any compact subset Ω 0 (A(Ω 0 ) = π, and Ω 0 ⊆ Ω) is given by the following formula: In the next sections, some special cases of the general model are presented.

S Ω L = π Ω • S L (Ω),

VI. GATEWAYS DEPLOYED ACCORDING TO REGULAR

PATTERNS By assuming that gateways are deployed according to 2D lattices, it is possible to easily assess the rate quantity defined in the previous section for surfaces of any size. Such regular arrangements of both the gateways and their associated coverage ranges shape the Euclidean plane as tessellated with "tiles" of the same form and area. The two cases pictured in Fig. 8 will help the following discussion.

Focusing of the equilateral triangular lattice of Fig. 8a (as an example, yet similar reasoning can be done for the square lattice of Fig. 8b), the nearby gateways A, B, and C define a tile Ω T , i.e., the blue equilateral triangle, that can be repeated indefinitely to cover the Euclidean plane. Once computed the value of (49) on the considered triangular tile, it is possible to easily compute the same index for any compact subset Ω k T formed by a set of k tiles. Indeed, the measure of the compact subset Ω k T is the k-th multiple of the measure of the compact subset Ω T , i.e., A(Ω k T ) = k•A(Ω T ). As well, it can be quickly verified that the rate S L (Ω k T ) of successful transmission to at least L gateways related to Ω k T is the k-th multiple of the rate S L (Ω T ) of successful transmission to at least L gateways related to Ω T , i.e., S L (Ω k T ) = k • S L (Ω T ). As a consequence, considering a sufficiently large compact subset S L (Ω T ), the quantity defined by (54) can be derived as:

S Ω T L = π Ω k T • S L (Ω k T ) = = π k • Ω T • k • S L (Ω T ) = = π Ω T • S L (Ω T ). (55)
As already anticipated, a similar rationale can be used for the deployment pictured in Fig. 8b, where the significant tile is the square having as vertices A, B, C, and D.

As matter of facts, the two deployments pictured in Fig. 8 are the only ones that allow tiling the Euclidean plane with regular polygons. 1 In the remaining part of this section, the considered honeycomb and square deployments are featured in many details, respectively in Sec. VI-A and Sec. VI-B.

A. Honeycomb deployment

When gateways are arranged according to an equilateral triangular lattice, the resulting honeycomb deployments are organized in the forms pictured in Fig. 9a and10a. Assuming that the transmission range R is an immutable property of the underlying transmission technology, it is possible to design the arrangement of gateways in a way that any end device of the LoRa deployment falls into the coverage area of at least L gateways. This is done by properly setting the distance d between two nearby gateways.

The honeycomb structure for gateways of Fig. 9a is obtained by setting the distance between two nearby gateways to the minimum value insuring that each point on the Euclidean plane is covered by at least L = 1 gateway. The equilateral triangle ABC is the "tile" used for computing the rate of successful transmissions of (55). To achieve such a coverage, it is first worth noting that the farthest point from A, B and C within the equilateral triangle ABC is its centroid P . In facts, the distance d between two nearby gateways, e.g., A and B, must be set in a way that permits the coverage of P by at least one gateway. This design constraint sets up the distance d to √ 3 times the transmission range R. The resulting tile is partitioned into compact subsets covered by just 1 gateway (colored with red), and compact subsets covered by 2 gateways (colored with green). In such a deployment, referred to with H1, the average rate of successful transmissions to at least 1 gateway is indicated as S H1

1 . The analytical expression of the latter can be found by using (55) combined with (49), and it is reported in Tab. I. Similarly to the validation of Sec. IV-C, such a model has also been validated through simulations, as shown in Fig. 9b. It has to be specified that the simulated deployment area is a square whose side is 6 times the transmission range.

The measurements are instead related to the inner centered square with side equal to 2 times the transmission range. Such a restriction on the evaluation is done to avoid some border effect: in the area closer to the borders of the whole simulated area, the end devices would have more chances to see their frames being received by gateways, since there are less interfering devices around the gateways close to the same borders [START_REF] Accettura | Simulating scalable long range wide area networks for very low power monitoring applications[END_REF].

If the design requirement on the honeycomb deployment is to get a coverage of the plane to allow the reachability of any point by at least L = 2 gateways, the distance between any pair of nearby gateways must be further reduced. Referring to Fig. 9a, to achieve a complete coverage of the plane with at least 2 gateways, the positions A and B where two nearby gateways are placed must come closer. This requirement is translated in enforcing a design constraint: the position of a gateway falls under the coverage range of any nearby gateway in the equilateral triangle lattice. The minimum distance satisfying such requirement is d = R, as pictured in Fig. 10a. As matter of facts, setting up d = R is also the minimum distance allowing the coverage of any point on the

= 4 √ 3π 9 gµ π 2 Q (π) - π 2 - 3 √ 3 4 Q 5π 3 + √ 3 2 R 1 S H3 1 = 4 √ 3π 3 gµ π 2 Q (π) -π - 3 √ 3 4 Q 4π 3 + √ 3 2 - π 2 - 3 √ 3 4 Q 5π 3 + √ 3 2 + [37] + π 2 - √ 3 2 Q 3π 2 + √ 3 + π - 3 √ 3 2 Q 5π 3 + √ 3 - π 2 - 3 √ 3 4 Q 5π 3 + 3 √ 3 2 2 S H3 2 = 4 √ 3π 3 gµ π - 3 √ 3 4 Q 4π 3 + √ 3 2 + π 2 - 3 √ 3 4 Q 5π 3 + √ 3 2 -2 π 2 - √ 3 2 Q 3π 2 + √ 3 + -2 π - 3 √ 3 2 Q 5π 3 + √ 3 + 3 π 2 - 3 √ 3 4 Q 5π 3 + 3 √ 3 2 3 S H3 3 = 4 √ 3π 3 gµ π 2 - √ 3 2 Q 3π 2 + √ 3 + π - 3 √ 3 2 Q 5π 3 + √ 3 -3 π 2 - 3 √ 3 4 Q 5π 3 + 3 √ 3 2
[36], [START_REF] Accettura | Simulating scalable long range wide area networks for very low power monitoring applications[END_REF] Euclidean plane by at least L = 3 gateways. As shown in Fig. 10a, such a considered honeycomb deployment, referred to with the notation H3, partitions the resulting ABC tile into compact subsets covered by 3 gateways (colored with blue), and compact subsets covered by 4 gateways (colored with yellow). In this deployment, it is possible to compute the average rate of successful transmissions to at least 1, 2, or 3 gateways. Such quantities are referred to respectively with S H3 1 , S H3 2 , S H3 3 , they are reported in Tab. I, and they are respectively plot and validated in red, green and blue in Fig. 10b. The throughput formula for this case but without duty cycle limitation was derived in [START_REF] Accettura | Simulating scalable long range wide area networks for very low power monitoring applications[END_REF] when packets have to be received by at least one gateway and for the case of at least three gateways in [START_REF] Accettura | Low power radiolocation through long range wide area networks: A performance study[END_REF]. Their formulas are a special case of (49) applied to this specific configuration with d = R. For the sake of completeness, they have been mentioned in Tab. I.

With a similar reasoning, it is possible to compute the minimum distance that allows the coverage of any point on the plane by at least 4 gateways. With reference to Fig. 10a, the lattice of gateway positions must be shrunk to allow the disk centered at D to "touch" the centroid of the ABC triangle. This is obtained by setting the distance d to be equal to √ 3/2 times R. The graphical result is shown in Fig. 11. However, the analytical results related to the average rate of successful transmissions to at least L = 1, 2, 3, 4 gateways are not reported, since their computation is toilsome and does not add any significant value from a scientific point of view.

B. Square deployment

In square deployments, the basic tile used for tessellation is in facts a square, as pictured in Fig. 12a, 13a, and 14a. In the remaining part of this subsection, two nearby gateways are those lying on the same side of a square tile (e.g., A and B in Fig. 14a), and whose distance is indicated as d.

The minimum distance d among nearby gateways allowing coverage by at least L = 1 gateway in a square deployment is obtained by imposing that the coverage area centered at a gateway must reach the farthest point in the square tile from its vertices, i.e., its centroid. This is reported in Fig. 12a, where the square ABCD is the tile used for tessellation and P is its centroid. As it can be seen, setting the distance d to √ 2 times the transmission range R permits the coverage of any point on the plane by at least 1 gateway. The resulting deployment is indicated as S1 and the related average rate of successful transmission to at least 1 gateway is referred to with S S1

1 . The exact expression of such value is reported in Tab. II, and it is plot and validated in Fig. 12b.

To allow a coverage of the plane with at least L = 2 gateways for any point, the minimum distance d between nearby gateways is obtained by imposing that the position of gateway on the plane falls under the coverage range of its own nearby gateways. Fig. 13a shows the resulting S2 deployment, where d = R. On such a deployment, it is possible to compute the average rate of successful transmission to at least 1 gateway (referred to with S S2 1 ), and 2 gateways (referred to with S S2

2 ). The exact expression of such models is reported in Tab. II, and S S2 1 and S S2 2 are respectively plot and validated in red and green in Fig. 13b.

Reducing the distance of nearby gateways to 2 √ 5/5 times the transmission range R, the resulting S3 square deployment permits each point of the plane to fall under the coverage range of at least 3 gateways, as shown in Fig. 14a. The average rate of successful transmissions to at least L = 1, 2, 3 gateways are respectively indicated S S3

1 , S S3 2 , and S S3 3 , they are reported in Tab. II, and they are respectively plot and validated in red, green and blue in Fig. 14b. In the S2 deployment of Fig. 13a, the green zones are those covered by just 2 gateways. By bringing gateways closer each others, green zones becomes smaller and smaller until they disappear. This is the case pictured in Fig. 14a. From a geometrical point of view, the minimum distance d between nearby gateways that allows a coverage of any point in the deployment by at least L = 3 gateways is achieved by imposing that the coverage area around a vertex of a square tile touches the middle point of the non adjacent sides in the same square tile. In other words, when the disk centered at A, touches the middle point P of the segment BC.

In order to permit a complete coverage of the deployment with at least 4 gateways, the distance between nearby gateways must be further reduced. In particular, the minimum distance d that fits this requirement is that resulting from the accomplishment of the graphical condition shown in Fig. 15: the circles that bound the coverage areas centered at A, E, and F must intersect at the same point P . In other words, the blue areas (those covered by 3 gateways) of 2, 3, and 4 gateways. The exact computation of the rate of successful transmissions for S4 deployments is not reported, since this would be toilsome and, in any case, it would be not relevant from a scientific point of view.

VII. CONCLUDING REMARKS AND FUTURE WORKS

Throughput formulas for multi-channel multi-gateway Lo-RaWAN using ALOHA with duty cycle limitation were derived in this paper. They give insights into the scalability of LoRaWAN and can serve as guidelines for determining the right density of end devices the network can support given a constraint on the minimum throughput.

The analysis was done under the assumptions of nobuffering at the end devices as well as with ACK mode disabled. It would be interesting to study the effects of these features on the throughput. Other possible lines of research are to include some of the specific characteristics of wireless channels such as fading and capture effects.

√ 2R 1 S S1 1 = π 2 gµ πQ (π) -(π -2) Q 3π 2 + 1 R 1 S S2 1 = πgµ πQ (π) - 4π 3 - √ 3 Q 4π 3 + √ 3 2 -(π -2) Q 3π 2 + 1 + 5π 3 -2 √ 3 Q 19π 12 + √ 3 2 + 1 + - π 3 - √ 3 + 1 Q 5π 3 + √ 3 + 1 2 S S2 2 = πgµ 4π 3 - √ 3 Q 4π 3 + √ 3 2 + (π -2) Q 3π 2 + 1 -2 5π 3 -2 √ 3 Q 19π 12 + √ 3 2 + 1 + +3 π 3 - √ 3 + 1 Q 5π 3 + √ 3 + 1 2 √ 5 5 R 1

Fig. 1 .

 1 Fig. 1. Evolution of the state of a node

Fig. 2 .

 2 Fig. 2. Cycle states.

  is matched[START_REF] Zhang | Multichannel aloha data networks for personal communications services (pcs)[END_REF] and (13) are matched (14) is matched(b) n = 3.

Fig. 3 .

 3 Fig.3. Acceptance regions for[START_REF] Zhang | Multichannel aloha data networks for personal communications services (pcs)[END_REF] and[START_REF] Crowther | A system for broadcast communication: Reservation-aloha[END_REF]

Fig. 4 .

 4 Fig. 4. Single gateway deployment.

Fig. 5 .

 5 Fig. 5. Region where LoRaWAN for a given ε has a better throughput than ALOHA. n = 1. A = 1.

Fig. 6 .

 6 Fig.6. Throughput in single-gateway deployments.

Fig. 7 .

 7 Fig. 7. Example of deployment of gateways.

Theorem 4 .

 4 Given two subsets Γ and B of the set Γ * of gateways deployed on the Euclidean space R 2 , if Γ = B, then the related assigned compact subsets ∆ Γ ⊂ R 2 and ∆ B ⊂ R 2 are disjoint, i.e., ∆ Γ ∩ ∆ B = ∅:

  thus proving the statement. If L = |Γ|, the proof is finished. Otherwise, consider the value F > L. Assume that c Φ,L = c |Φ|,L for any set Φ with |Φ| < F . The coefficient c Φ,L with |Φ| = F is then:

  with A(Ω) π ∧ ∀α ∈ Ω : |σ(α)| L. (54)

Fig. 8 .

 8 Fig. 8. Regular grids.

  th, L=1 1% DC, th, L=1 No DC, sim, L=1 1% DC, sim, L=1 (b) Simulation results validation.

Fig. 9 .

 9 Fig. 9. Honeycomb deployments with d = √ 3R, L = 1.

Fig. 10 .

 10 Fig. 10. Honeycomb deployments with d = R, L = 1, 2, 3.

  th, L=1 1% DC, th, L=1 No DC, sim, L=1 1% DC, sim, L=1 (b) Simulation results validation.

Fig. 12 .

 12 Fig. 12. Square deployments with d = √ 2R, L = 1.

  Fig. 14a collapse to single points (like P ) when reducing the distance d to equal 3 √ 2/5 times the transmission range R. Such deployment, referred to as S4, allows the computation of S S4 1 , S S4 2 , S S4 3 , and S S4 4 as the rate of successful transmission to at least (respectively) 1, th, L=1 1% DC, th, L=1 No DC, th, L=2 1% DC, th, L=2 No DC, sim, L=1 1% DC, sim, L=1 No DC, sim, L=2 1% DC, sim, L=2 (b) Simulation results validation.

Fig. 13 .

 13 Fig. 13. Square deployments with d = R, L = 1, 2.

  transmissions to at least L gateways

  th, L=1 1% DC, th, L=1 No DC, th, L=2 1% DC, th, L=2 No DC, th, L=3 1% DC, th, L=3 No DC, sim, L=1 1% DC, sim, L=1 No DC, sim, L=2 1% DC, sim, L=2 No DC, sim, L=3 1% DC, sim, L=3 (b) Simulation results validation.

Fig. 14 .

 14 Fig. 14. Square deployments with d = 2 √ 5 5 R, L = 1, 2, 3.

Another viable tiling with regular polygons would be the one that uses hexagons, yet an hexagonal tiling is a special case of the tiling with equilateral triangles.

, L = 1, 2, 3, 4.