Rosa Córdoba

Alfonso Ibarra

Dominique Mailly

Isabel Guillamón

Hermann Suderow

José María De

José María De Teresa
email: deteresa@unizar.es

3D superconducting hollow nanowires with tailored diameters grown by focused He+ beam direct writing

Keywords: electron tomography, focused ion beam induced deposition (FIBID), helium ion microscope, magneto-transport measurements, nano-superconductors, tungsten carbide (WC)

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

6.8 K and show high values of critical magnetic field and critical current density. Consequently, these 3D nano-objects could be implemented as components in the next generation of electronics, such as nano-antennas and sensors, based on 3D superconducting architectures.

Introduction

Superconductors are dissipationless carriers of electric current and provide macroscopic, and thus robust, quantum coherence. This allows for a wide range of applications, particularly at the nanometer-scale, where they can be easily integrated in circuits and used as ultrasensitive sensors of magnetic fields, temperature and as key elements for quantum computation. The behavior of nanosized superconductors as one-dimensional quantum oscillators [1], Josephson junction arrays [2], electronic transport devices [3][4][5][6][7], very small-scale devices [8,9], micrometerscale coolers [10], or thermal and spin sensors [11,12] has been studied in detail.

Nowadays, research on manufacturing highly energy-efficient three-dimensional (3D) structures [13] is critical for the development of future electronics. However, when approaching the nanometer-scale, the number of works on real 3D nano-superconductors [14][15][16][17][18][19] decreases dramatically, mostly due to the complex fabrication and characterization. A technique successfully utilized for fabricating 3D nano-objects is direct writing by a focused beam of positively charged particles, the so-called focused-ion-beam induced deposition (FIBID) [20]. A very promising development of FIBID is based on Ga + ions. Functional 3D nanomaterials have been grown by Ga + FIBID in the last decade [START_REF] Matsui | Three-Dimensional Nanostructure Fabrication by Focused-Ion-Beam Chemical Vapor Deposition[END_REF][START_REF] Morita | [END_REF][23][24][25][26]. In particular, Ga + FIBID in combination with W(CO) 6 as precursor material yielded 3D superconducting W-based wires with a critical temperature (T c ) below 5 K and a critical magnetic field (µ 0 H c2 (0)) up to 9.5 T [14][15][16]. Alternatively, in combination with Nb(NMe 2 ) 3 (N-t-Bu), Ga + FIBID yielded NbC wires with a broadened T c range from 4 to 11 K [18]. One significant limitation is that 3D elements below 100 nm in diameter cannot be obtained with Ga + FIBID, mainly due to the relatively large Ga + beam diameter (approx. 5 nm) and a high proximity effect generated by Ga + ion scattering.

Regarding a higher spatial resolution, the helium ion microscope (HIM) [27], based on a gas field-ionization source, has emerged as a tool for direct writing of complex 3D nano-objects taking advantage of its small beam diameter (approx. 0.3 nm) and low proximity effect [START_REF] Maas | Nanofabrication with a Helium Ion Microscope[END_REF]. When precursor molecules from the gas phase are adsorbed on a substrate surface, He + FIB dissociates them into non-volatile and volatile products. The non-volatile products attach to the surface, resulting in a deposit, whereas the volatile products ones are pumped out of the process chamber. Normally, the final deposit is a mixture of carbon, metallic elements and oxygen. As clearly described using analytical modelling [29] and Monte-Carlo simulations [30], the vertical growth of 3D nano-objects by He + FIBID is mainly caused by secondary electrons of the first order produced from the primary ion beam, whereas the lateral growth is induced by secondary electrons of the second order generated from scattered ions. Thus, the direct contribution of the primary He + ion beam and the scattered He + ions is almost negligible for the growth of these 3D nano-objects. Nevertheless, it is worth mentioning that its resolution, volume per dose and throughput are very sensitive to the selected growth conditions such as ion beam energy, ion beam current, precursor flux, surface interactions with the beam, and precursor molecules [29,30]. Hence, the He + FIBID technique is highly recommended for direct writing of 3D nano-objects with high resolution and aspect ratio [17,19,[31][32][33][34][35]. A successful example of tailored 3D nano-objects grown by He + FIBID has been reported by Kohama and co-workers [35]. The authors deposited W-based pillars with diameters down to approx. 40 nm and aspect ratios of approx. 50. The microstructure of the grown material consisted of fcc WC 1-x and W 2 (C,O) grains. Moreover, when the He + beam was well focused the authors observed columnar voids created at the center of the pillars with a diameter ranging from 1 to 15 nm, showing the path to build complex 3D nano-objects beyond standard nanowires (NWs). Recent breakthroughs in the growth of 3D WC superconducting nano-objects with extremely large aspect ratios using He + FIBID have been reported by some of the authors, such as hollow NW-like nanotubes as small as 32 nm in diameter [17] and nanohelices with controllable geometries, including the smallest and most densely packed nanohelix to date with a diameter of 100 nm [19].

In this work, we present the direct writing of 3D WC crystalline superconducting hollow NWs with tailored diameters grown using a HIM. The hollow NW geometry is successfully controlled by tuning the ion beam current and dose from 0.65 to 7 pA and from 0.1 to 0.4 nC, respectively, resulting in NWs with outer diameters from 36 to 142 nm and with inner diameters from 5 to 28 nm, and total length from 0.5 to 8.9 µm (aspect ratio ≈ 196). These values are significantly better than those reported in previous works [17,35]. The NWs microstructure consists of large grains of fcc WC 1-x , in good agreement with [17,35]. In addition, the NWs are hollow along the whole NW length, which could make them nonconventional nanopipettes, as demonstrated in 3D reconstructions of electron tomography experiments. Finally, these 3D hollow NWs exhibit superconductivity below 6.8 K (T c ), as well as high upper critical magnetic fields µ 0 H c2 ≈ 14.7 T, and large critical current densities J c ≈ 0.15 MA/cm 2 .

Results and Discussion

Growth of 3D hollow nanowires by He + FIBID

We use a HIM in combination with a W(CO) 6 precursor to grow individual, out-of-plane WC NWs in a single step, controlling inner and outer diameter and total length. The precursor gas is delivered to the process chamber and adsorbed onto the substrate surface, while the He + FIB spot remains fixed during the deposition favoring continuous vertical growth along [17].

Dimensional control for nanowires

We investigated the dimensional control for these NWs by optimizing in the deposition the following parameters: the ion beam current and ion dose. SEM images of typical NWs grown with ion beam currents ranging from 0.54 to 6.47 pA and doses from 0.1 to 1.4 nC are depicted in Figure 1a-d.

Varying these parameters enables us to fabricate 3D NWs with diameters ranging from 45 to 125 nm, lengths ranging from 0.5 to 8.9 µm, and with aspect ratios up to 198. Further details regarding the growth conditions are described in the Experimental section. We found a linear dependence of the NW volume (determined as π × (outer diameter/2 -inner diameter/2) 2 × height) as a function of the ion dose for the mentioned ion beam currents (Figure 1e). Moreover, we noted that the NW volume rapidly decreases as a function of the ion beam current for the same dose (Figure 2). When using high currents several effects can play a role in this dependence such as precursor depletion, local heating, which decreases the precursor molecule sticking coefficient, and low precursor diffusion from the substrate to the top of the pillar [36,37]. This shows the need for future systematic experiments varying the dwell time in pulsed growth or varying the flux of precursor gas.

(High-resolution) scanning transmission electron microscopy

Dependence of NW inner diameter on the ion beam current

To investigate the dependence of the NW inner diameter on the ion beam current, scanning transmission electron microscopy (STEM) experiments were performed. We found that inner diameter of the hollow NWs changes from 5 to 28 nm, whereas the outer diameter changes from 36 to 143 nm upon increasing the ion beam current from 1.3 to 7 pA. STEM images of these hollow NWs are shown in Figure 3a. The observed non-uniform shape of the cavity in the central nanowire could be explained by several reasons, such as He + FIB instability or irregular substrate surface. We find a linear dependence of the inner diameter on the ion beam current (Figure 3b), which indicates that the ion beam forms the cavity due to a milling effect. Thus via tuning the ion beam current and dose we have full control to tailor the diameters of the hollow 3D NWs. The specific deposition parameters and NW diameters are listed in Table 1.

Electron tomography

In order to further examine the NW diameters along their length, electron tomography experiments on two specific NWs were carried out. Figure 4 shows the tomographic reconstruction of hollow WC NWs grown at (a) 7 pA and (b) 2 pA. One can see from the images that the cavities are present up to the tip of the NW. On the left panel of Figure 4a, a STEM image of the NW with outer and inner diameter of 142 and 28 nm, respectively, is shown. On the right panel, a snapshot of the colored 3D tomographic reconstruction is depicted. Figure 4b shows a STEM image of the NW with outer and inner diameter of 77 and 8 nm, respectively, on the left panel, and a snapshot of the colored 3D tomographic reconstruction is displayed on the right panel. Three movies of the tomographic reconstruction for each hollow NW are added in Supporting Information File 1-7, including a transversal (x-y) and a longitudinal (y-z) section, and a colored three-dimensional reconstruction.

Microstructure

Concerning the microstructure of the NWs, high-resolution scanning transmission electron microscopy (HRSTEM) images have been acquired sequentially and processed to extract the crystallographic structure (Figure 5). We indexed the spots indicated in the fast Fourier transform (Figure 5b) of the image in Figure 5a with the planes (-11-1), (-200) and (-1-11), and the [011] zone axis of the WC 1-x fcc structure, with a lattice parameter of a = 0.4272 nm. A lower magnification STEM image of the NW grown at 1.3 pA is depicted in Figure 5c. These results are in good agreement with the previous work reported by some of the authors [17].

Magneto-electrical-transport study

To determine the critical superconducting parameters in NWs grown at 0.65, 1.3, and 2.18 pA (Figure 6 and Table 2), a magneto-electrical transport study using the typical four-pointprobe configuration has been performed. Following the procedure described in [17], first 3D NWs were placed flat on the SiO 2 layer of a Si/SiO 2 substrate by means of a nano-manipulator. Then, four Pt FIBID contacts were grown to connect the NWs to pre-patterned Ti pads. Finally, we made four-point-probe electrical measurements at low temperature (down to 0.5 K) and under a magnetic field perpendicular to the substrate plane (up to 9 T).

The NWs change from the normal to the superconducting state at T c (0.5R N ) values between 5.45 and 6.78 K (Figure 6a and Table 2). No clear trend was visible in T c values for NWs grown using different currents, although the identified T c range is in good agreement with the previously reported results [17]. Also, it is up to 1.5 times higher than that of Ga + FIBID nanostructures of similar dimensions [9]. The inset of Figure 6a shows the measured resistance as a function of the temperature in the full temperature range investigated for a NW grown at 0.65 pA. The value of µ 0 H c2 as a function of the temperature for NWs grown at 0.65 and 1.3 pA is depicted in Figure 6b. The values of µ 0 H c2 (0.9R N ) are extracted from the resistance-vs-temperature curves under perpendicular magnetic field (inset of Figure 6b). By fitting µ 0 H c2 (T) to a power-law equation µ 0 H c2 (T) ∝ (1 -T/T c ) n , µ 0 H c2 (0 K) is estimated to be approx. 14.5 T for the different NWs. The coherence length, ξ, at 0 K is around 4.77 nm and the estimated magnetic field penetration depth, λ, [38,39] ranges from 720 to 874 nm. Additionally, J c (0.6 K, 0 T) is approx. 0.15 MA/cm 2 .

Summarizing, the estimated superconducting parameters (T c , µ 0 H c2 , J c , ξ, λ) for these NWs (Table 2) are mostly compatible with those reported for He + FIBID out-of-plane WC nanotubes [17], nanohelices [19], and in-plane NWs used in hybrid microwave resonators [40]. They are potential building blocks for highly packed 3D nano-resonators, superconducting logic gates [41], quantum switches [42], and single-photon detectors [43][START_REF] Dauler | Superconducting nanowire single photon detectors[END_REF][START_REF] Natarajan | [END_REF].

Conclusion

We report a direct writing methodology to create 3D superconducting hollow NWs with tailored diameters using W(CO) 6 precursor with a highly focused He + beam. The resulting 3D hollow NWs have inner and outer diameters from 5 to 28 nm and from 36 to 142 nm, respectively, and aspect ratios above 196, which is unachievable by other additive manufacturing methods. The electron tomography study proved that the center hole is present along the whole length of the NWs.

As expected, the microstructure corresponds to the fcc WC 1-x phase. By studying their magnetotransport properties, we found T c ≈ 6.8 K, as well as µ 0 H c2 ≈ 14.7 T and J c ≈ 0.15 MA/cm 2 . The presented methodology yields an advanced bottom-up approach for the fabrication of innovative 3D nano-architectures, in which nano-superconductivity may provide an advantage, for future electronic components, particularly for sensors, energystorage components, and quantum computing.

Experimental Growth of 3D hollow WC nanowires

He + FIBID hollow WC NWs have been fabricated in a ZEISS ORION NanoFab instrument equipped with a helium ion beam column and a single-needle gas injection system (GIS) through which W(CO) 6 gas is delivered to the process chamber.

The NWs were deposited on top of the pre-patterned Ti pads (150 nm in thickness) to prevent charge effects on the insulator layer (250 nm thick of SiO 2 ) thermally grown on a silicon wafer [23]. These chips were fabricated following a routine recipe for UV optical lithography using a lift-off method. For the electron tomography and (HR)STEM experiments, NWs were directly grown on Cu TEM grids. Typical deposition conditions used for the He + FIBID process were as follows; precursor material: tungsten hexacarbonyl, W(CO) 6 ; T precursor = 55 °C; GIS needle diameter ≈ 500 µm; GIS z ≈ 500 µm; GIS x,y ≈ 60 µm; P base ≈ 3 × 10 -7 mbar; P process ≈ 4 × 10 -6 mbar; acceleration voltage = 30 kV; pattern shape: spot mode; ion beam current range of 0.54 to 6.47 pA and dose range of 0.1 to 1.4 nC.

Microstructure and tomography at the nanometer-scale

Scanning transmission electron microscopy (STEM) imaging and EDS were carried out in a probe-corrected FEI Titan 60-300 operated at 300 kV and equipped with a high-brightness X-FEG and a CETCOR C s corrector for the condenser system to provide sub-angstrom probe size.

STEM high-angle annular dark-field (HAADF) tomography was performed using a Thermo Fisher Tecnai field-emission gun operated at 300 kV. The angular range for the tilt series was ±70° with pictures taken every 1°. Image alignment and 3D reconstruction was carried out with FEI tomography acquisition software Inspect 3D after the acquisition of 140 images. The movies of the tomographic reconstruction for each hollow NW were performed using Amira 3D software.

Magneto-electrical transport study

The magneto-electrical transport measurements on the NWs were carried out using a ''Physical Property Measurement System'' (PPMS), from Quantum Design equipped with a helium-3 refrigerator insert.

regional Gobierno de Aragón (grants E13_20R and E28-20R) with European Social Fund (Construyendo Europa desde Aragón) and Comunidad de Madrid through project NANOMAGCOST-CM (Grant No. (S2018/NMT-4321)). The project that gave rise to these results received the support of a fellowship from "la Caixa" Foundation (ID 100010434). The fellowship code is LCF/BQ/PR19/11700008. The French RENATECH network (French national nanofabrication platform).

ORCID ® iDs

Rosa Córdoba -https://orcid.org/0000-0002-6180-8113 Alfonso Ibarra -https://orcid.org/0000-0002-4599-3013

Hermann Suderow -https://orcid.org/0000-0002-5902-1880

José María De Teresa -https://orcid.org/0000-0001-9566-0738

Figure 1 :

 1 Figure 1: (a-d) SEM images of hollow NWs grown by He + FIBID; (a) ion beam current: 0.54 pA, dose: 0.1-0.3 nC from right to left; (b) ion beam current: 1.17 pA, dose: 0.1-0.6 nC from left to right and upwards; (c) ion beam current: 2.91 pA, dose: 0.1-0.6 nC from left to right and upwards; (d) ion beam current: 6.47 pA, dose: 0.229, 0.456, 0.696, 0.935, 1.158, and 1.402 nC from left to right and upwards. (e) NW volume in (a) as a function of the ion beam dose.

Figure 2 :

 2 Figure 2: NW volume for a specific dose, as a function of the ion beam current.

Figure 3 :

 3 Figure 3: (a) STEM images of hollow NWs grown at 1.3, 2.3, and 7 pA, from left to right. (b) NW diameter in (a) as a function of the ion beam current.

Figure 4 :

 4 Figure 4: (a) Tomography of a hollow WC NW grown at 7 pA, with an outer diameter of 142 nm and an inner diameter of 28 nm; left panel: STEM image, right panel: snapshot of the 3D tomographic reconstruction. (b) Tomography of a hollow WC NW grown at 2 pA, with an outer diameter of 77 nm and inner diameter of 8 nm; left panel: STEM image, right panel: snapshot of the 3D tomographic reconstruction.

Figure 5 :

 5 Figure 5: (a) HRSTEM image of a typical hollow WC NW grown at 1.3 pA. (b) Fast Fourier transform of the image in (a), showing the crystalline nature of the material and indexed as the [011] zone axis of the fcc WC 1-x structure. (c) Lower magnification STEM image of the WC NW in (a).

Figure 6 :

 6 Figure 6: (a) Normalized resistance as a function of the temperature of hollow NWs grown using ion beam currents indicated in the legend. R N is the resistance for the normal state at T = 10 K. The inset shows the resistance as function of the temperature (full range) for a NW grown at 0.65 pA, I bias = 100 nA. (b) Upper critical magnetic field (µ 0 H c2 ) as a function of the temperature of NWs grown using ion beam currents indicated in the legend. Data is fitted to a power-law equation. The inset shows the resistance as a function of the temperature for a NW grown at 0.65 pA under perpendicular magnetic fields from 0 to 9 T.

Table 1 :

 1 Growth parameters and diameters of hollow WC NWs.

	hollow nanowire type	1	2	3
	ion beam current (pA)	1.3	2.3	7.0
	outer diameter (nm)	36	71	142
	inner diameter (nm)	5	9	28

Table 2 :

 2 Superconducting parameters of NWs estimated from experimental magneto-transport measurements.

	hollow
	nanowire type 1	2	3
	ion beam	0.65	1.30	2.18
	current (pA)
	outer diameter	36	41	72
	(nm)
	R N (Ω)	1959	2101	430
	T c (K)	6.78	6.16	5.45
	µ 0 H c2 (0 K) (T) 14.66	14.48	-
	(MA/cm 2 ) J c (0 T)	0.083 (0.6 K) 0.151 (0.6 K) 0.026 (1 K)
	ξ (0 K) (nm)	4.74	4.77	-
	λ (0 K) (nm)	839	720	874

Acknowledgements

The authors highly acknowledge Rubén Valero for the UV lithography process. The microscopy works have been conducted in the "Laboratorio de Microscopías Avanzadas" at "Instituto de Nanociencia de Aragón -Universidad de Zaragoza". The authors acknowledge the LMA-INA for offering access to their instruments and expertise. The authors would like to acknowledge the use of Servicio General de Apoyo a la Investigación-SAI, Universidad de Zaragoza, particularly the Servicio de Medidas Físicas.

Funding

This work was supported by the financial support from Spanish Ministry of Economy and Competitiveness through the projects MAT2017-82970-C2-2-R, PIE201760E027, including FEDER funding, FIS2017-84330-R, CEX2018-000805-M, EU ERC (Grant Agreement No. 679080), from the EU-H2020 research and innovation programme under grant agreement No 654360 NFFA-Europe, and Nanocohybri COST CA-16218, from

Supporting Information

Movies of 3D tomographic reconstruction.

Supporting Information File 1


René Lohmann

Beate Horn-Cosfeld

Jakob Schluck

Dominique Mailly

Nima H Siboni

Hans W Schumacher

Klaus Pierz

Thomas Heinzel

Jürgen Horbach

Rapid Communications

The magnetotransport of an electron gas in a two-dimensional, random arrangement of overlapping retroreflective crosses is studied using a combination of experiment and classical event-driven molecular dynamics simulation. The experimentally measured magnetoconductivity σ xx as a function of the magnetic field B displays an anomalous behavior at low B fields, accompanied with an increase of σ xx with increasing temperature at B = 0. The simulations show that at B = 0 the magnetoconductivity does not exist and is associated with anomalous diffusion in the asymptotic long-time limit that depends on the orientational order of the obstacles. At any finite B field, the motion of the tracer particle is asymptotically diffusive but the diffusion coefficient vanishes with a power law in the limit B → 0.

Introduction. Anomalous diffusion processes such as suband superdiffusion are ubiquitous in heterogeneous media [1][START_REF] Sahimi | Flow and Transport in Porous Media and Fractured Rock: From Classical Methods to Modern Approaches[END_REF][START_REF] Höfling | [END_REF]. However, there are only a few particle-based model systems that exhibit anomalous diffusion as an asymptotic process in the long-time limit. A celebrated example in this context is the disordered Lorentz gas (LG) [4] where, in its simplest version, a tracer particle probes the void space of a random arrangement of overlapping spherical objects. In this model, a conductor-insulator transition at a critical obstacle density n c coincides with the percolation transition of the obstacles, implying subdiffusion due to the fractal structure of the free space at n c [START_REF] Höfling | [END_REF][5][6][7][8][9][10][11][12][13][14][15][16].

Evidence for a completely different type of asymptotic subdiffusion has been found for Ehrenfest's wind-tree (EWT) model [START_REF] Ehrenfest | Collected Scientific Papers[END_REF]. Here, a point tracer particle moves through the void space of a two-dimensional (2D) random arrangement of overlapping or nonoverlapping squares. The flat boundaries of the squares induce retracing events (or retroreflections) where the tracer particle is reflected back such that even in short-time windows it returns several times to a region where it has been earlier. For overlapping squares, kinetic theory for low densities [START_REF] Hauge | [END_REF][19][20] predicts that the multiple backscattering events of the tracer particle result in a meansquared displacement (MSD) with an asymptotic sublinear growth in the long-time limit δr 2 (t ) ∝ t 1-4n /3 (with n = na 2 a dimensionless obstacle density, where n is the number density of the squares and 2a the length of their diagonals). An early simulation study [21] has qualitatively confirmed the existence of anomalous diffusion in the EWT model with overlapping squares.

These results for the EWT model raise many unsolved questions about anomalous transport in heterogeneous 2D media due to retroreflection. In previous studies, only the lowdensity limit of equally oriented squares has been considered, because this special case can be treated in the framework of kinetic theory. Using a combination of experiment and simulation, we go beyond this regime and consider intermediate and high densities of retroreflective obstacles. Moreover, we study the dependence of retroreflective transport on the obstacle's orientation. We consider 2D systems with random arrangements of overlapping retroreflective crosses (RCs) where effects due to retroreflection are expected to be more pronounced than in the EWT model [22][23][24].

As we shall see below, a key for the understanding of anomalous subdiffusive transport in retroreflective systems is the introduction of small perturbations from the straight trajectories by an external magnetic field B. The B field leads to a circular shape of the tracer particle's trajectory, with the cyclotron radius R cy ∝ B -1 acting as a characteristic length scale. We show that even at extremely large cyclotron radii, where the curvature of the trajectories is close to that of the straight lines at B = 0, the MSD grows always linearly in the long-time limit. As a consequence, we observe a power-law dependence of the magnetoconductivity at extremely low B fields. This anomaly comes along with an asymptotic subdiffusive behavior at B = 0.

A completely different issue is the investigation of the magnetotransport in an experimental realization of a retroreflective LG in terms of a two-dimensional electron gas (2DEG) with a lithographically imprinted obstacle structure [16,[25][26][27][28][29][30][31][32][33][34][35]. In this Rapid Communication, we present experiments on a 2DEG in the classical limit, exposed to a random arrangement of overlapping RCs, and compare them quantitatively to corresponding event-driven molecular dynamics (MD) simulations. In the experiments, background disorder prevents the occurrence of pure retroreflection and thus the presence of anomalous transport at B = 0. However, traces of anomalies due to retroreflection are seen in the behavior of the magnetoconductivity σ xx (B), such as an increase of σ xx (B = 0) with increasing temperature and a characteristic shoulder or peak of σ xx (B) at low B fields.

Experiments. A GaAs/Al 0.3 Ga 0.7 As heterostructure with a 2DEG 150 nm below the surface serves as the raw material. At liquid helium temperatures and after a brief illumination with infrared light, the pristine 2DEG has a density of n e = 2.5 × 10 15 m -2 and a typical electron mobility around 340 m 2 /V s, the exact value of which depends on the cooldown cycle as well as on the position of the particular sample on the wafer. This mobility corresponds to a mean free path of = 31 μm.

A quantum scattering length of q = 530 nm is extracted from the envelope of the Shubnikov-de Haas oscillations [36], indicating that remote donors are the primary source of scattering. A Hall bar of 100 μm width is prepared by optical lithography. Ni/AuGe Ohmic contacts allow electrical access to the 2DEG. Voltage probes of 10 μm width, with a spacing of 100 μm along the edge of the Hall bar, allow four-terminal measurements.

The Lorentz arrays are patterned by electron beam lithography and subsequent reactive ion etching (etch depth ≈ 150 nm). They are designed by first generating two sets of random numbers, with pairs of them defining the center coordinates of the obstacles. An RC, composed of two rectangles of size 300 nm × 1850 nm, is attached to each coordinate. Arrays with four number densities n = 0.0625, 0.125, 0.1875, and 0.25 μm -2 are defined by lithography and subsequent reactive argon ion etching with an etch depth of 150 nm, each with orientation order, where all RC edges form an angle of π/4 with respect to the current (x) direction, as well as without orientation order, i.e., these angles are randomly distributed. The lateral depletion of the 2DEG at the etched walls extends over a distance which is not known very well. We have estimated a depletion length of 75 nm [16,[37][38][39], corresponding to RCs of 2 μm size and an area of 1.6 μm 2 . In Fig. 1, a typical RC and an array thereof is shown. The dimensionless number density n , defined as the fraction of the area covered by the RCs in the absence of mutual overlaps, takes the values n = 0.1, 0.2, 0.3, and 0.4 for our arrays.

The transport measurements are carried out in the mixing chamber of a 3 He / 4 He dilution refrigerator with a base temperature of 25 mK as well as in a 4 He gas flow cryostat for temperatures above 1.2 K. Both cryostats are equipped with superconducting solenoids with maximum magnetic fields above B = 8 T. An AC current of amplitude 500 nA and a frequency of 17.7 Hz is applied, and the voltage drops across the arrays in longitudinal (x) and transverse (y) directions are recorded via the voltage probes using lock-in amplifiers.

Details of the simulation. For the event-driven MD simulations, randomly distributed overlapping RCs are placed in a rectangular box with dimension 200d o × 100d o , assuming periodic boundary conditions (here, the length scale d o corresponds to 1 μm in the experiments). The RCs have the same geometry as those in the experiment (Fig. 1). Both obstacle configurations with equally oriented and randomly oriented RCs are considered. Noninteracting tracer particles (electrons) with a velocity v are placed randomly in the free area between the crosses. For each B field and density n , considered below, we have simulated trajectories for 2000 independent host structures; in each structure we have placed 500 tracer particles. The initial velocity has a random direction and its magnitude is fixed to the Fermi velocity v F = | v|. The conversion of the Fermi velocity, the effective electron mass, and the electron charge to the experimental units are given by v F = 2.16 × 10 5 m/s, m = 6.097 × 10 -32 kg, and e = 1.6 × 10 -19 C, respectively (in simulation units, v F , m , and e are set to one). In the following, time is given in units of

t 0 = d o v F ≈ 4.63 ps. A dimensionless magnetic field is introduced via B = B/B 0 , with B 0 = m v F /(ed o ) ≈ 0.0823 T.
Results. Figure 1 displays trajectories for equally and randomly oriented RCs at B = 0 and n = 0.2 from the simulation. In both cases, one can identify retracing paths that, as we shall see below, lead to a subdiffusive transport in the long-time limit at zero magnetic field.

In the simulation, we study the diffusive transport of the tracer particles via the mean-squared displacement (MSD), defined as δr 2 (t ) = [ r(t )r(0)] 2 , with r(t ) the position of the particle at time t and • • • an ensemble average. These plots indicate that in all cases the MSDs do not grow linearly in the long-time limit, but as δr 2 (t ) ∝ t γ a with γ a < 1. In the case of the randomly oriented RCs, the simulation data suggest a slight increase of γ a towards the constant value γ a ≈ 0.73 at high density. For the equally oriented RCs, γ a decreases continuously with increasing obstacle density [see Fig. 2(d)].

To compare the behavior at finite magnetic fields to that at B = 0, we consider the derivative of the MSD (dMSD) dδr 2 (t )/dt. In the case of an asymptotic diffusive process, this quantity gives the diffusion coefficient D = lim t→∞ Figures 3(a) and 3(b) show the dMSD at n = 0.2 for the randomly and equally oriented RCs, respectively (note that we find a similar behavior for the other densities, too). Since for B = 0 the MSDs exhibit an asymptotic subdiffusive behavior, in this case the dMSD decays as t γ a -1 at long times. As we can infer from Figs. 3(a) and 3(b), the dMSDs for sufficiently low B fields first follow exactly the curve for B = 0 and then, above a timescale τ B they quickly reach a constant, corresponding to 4D . For B fields leading to a τ B smaller than the location of the maximum at τ m in the dMSD for B = 0, the behavior of the dMSDs changes qualitatively. This happens at about Bm = 0.09 for n = 0.2. Note that, for Bm , the time τ m corresponds to the time at which the tracer particle has made on average about 2.3 collisions with RCs. Thus, τ m can be considered as the elementary timescale for a retroreflection. The behavior for B > Bm is exemplified in the figures for B = 0.61. In the latter case, the dMSDs are only in agreement with the dMSD at B = 0 in the ballistic regime for t/t 0 < 1. Around t/t 0 = 1, there are some oscillations and the dMSD increases up to a timescale of the order of 10 3 , indicating a superlinear increase of the MSD in this time regime.

As Fig. 3(c) indicates, the timescale τ m is proportional to the inverse density, τ m ∝ 1/n . The magnetic field Bm is a linear function of n , Bm ∝ n [Fig. 3(d)].

The B dependence of the conductivity σ xx is plotted double logarithmically in Fig. 3(e). At small B, i.e., for B < 0.01, a power law, σ xx ∝ Bδ , is seen over about 3-4 orders of magnitude in B. While the data for the randomly oriented RCs can be fitted by power laws with a slightly decreasing exponent towards the value δ ≈ 0.09 at high densities, for the equally oriented RCs there is an increase of δ from about 0.15 to 0.21 for the considered densities [inset of Fig. 3(e)], reflecting the behavior of γ a at B = 0. Around Bm [marked by crosses in Fig. 3(e)], there is a shoulder in the data for the randomly oriented RCs and a local maximum in those for the equally oriented RCs. In this regime, the cyclotron radius is around 1/ Bm , i.e., of the order of 10d o . Here, the retroreflective character of the trajectories gets lost. Before the monotonic decay towards zero at higher B fields, another maximum can be seen in the curves for n 0.2. This feature is due to the change of the transport from a diffuse scatter by the obstacles to a directed motion along the contour that is also seen in systems with circular obstacles [16]. The question now arises as to what extent the transport anomalies, as seen in the simulations, are also visible in the experiment. Figures 4(a) and 4(b) compare the experimentally measured conductivities σ xx ( B) for equally oriented RCs at the different densities and the temperature T = 0.1 K to corresponding ones from the simulation (for experimental data of randomly oriented RCs, see the Supplemental Material [39]). With increasing RC density, the agreement between simulation and experiment becomes better. At n = 0.4, the discrepancies are always smaller than 15% for 0.1 < B < 2.4 [cf. Fig. 4(b)]. At very low B fields, however, the simulation exhibits the power-law behavior, as discussed above, while the experiment quickly saturates to a finite value in the limit B → 0. The fact that in the experimental system the conductivity at B = 0, σ xx (0), exists, indicates that there are background scattering processes that perturb the straight trajectories which are seen in the simulation at B = 0 (cf. Fig. 1). Therefore, the sequence of retracing events that lead to the nonexistence of σ xx (0) in the simulation are not present in the experimental system, and thus σ xx ( B) exhibits a regular behavior for B → 0 and at B = 0. Note that we can exclude effects due to weak localization in the experimental data, while the onset of Shubnikov-de Haas oscillation is shifted to magnetic fields above 0.2 T due to the strong scattering at the obstacles [39].

However, the increase of σ xx ( B) towards the shoulder in the interval 0.05 B 0.2 [cf. Fig. 3(d)] is a trace of the anomalous transport due to retroreflection in the experimental data at high densities. At n = 0.2 the experimental σ xx ( B) exhibits a maximum at B ≈ 0.1, in agreement with the simulation [cf. Fig. 4(a)]. With increasing temperature, background scattering events (e.g., due to phonons and local defects in the GaAs substrate) are expected to become more pronounced. In fact, as can be inferred from Fig. 4(c), at n = 0.3, with increasing temperature the shoulder around B = 0.4 disappears and at T = 32 K the conductivity decreases monotonically with increasing B. While the conductivity maximum caused by the directed motion along the RC's contours [16] is temperature stable up to 16 K, the shoulder around B = 0.4 is already suppressed at 4 K. This indicates the characteristic length scale of this process to be larger than the one corresponding to the process at B = 1. Moreover, the conductivity at B = 0 increases with increasing temperature [cf. Fig. 4(d) where σ xx (T ) at B = 0 is multiplied by n 2 to obtain a better visibility of the differences in the curves at different densities]. At the two highest densities, there is a significant decrease of σ xx towards low temperature, reflecting the fact that background scattering processes become less important with decreasing temperature as well as with increasing density. As a consequence, the traces of retroreflective transport become pronounced at low temperatures and high densities.

Discussion and conclusions. We have studied the magnetotransport of a tracer particle (electron) in 2D random arrangements of overlapping RCs, combining event-driven MD simulations with experiments on 2DEG. Our simulations reveal a system with anomalous transport due to retroreflective trajectories of the tracer particle, where at zero magnetic field the conductivity (or diffusion coefficient) does not exist at any density. Instead, subdiffusive transport is seen in the long-time limit. For randomly oriented RCs, the exponent for subdiffusion at B = 0 slightly increases towards the value γ a ≈ 0.73 at high density, while for the equally oriented RCs, γ a decreases with increasing density. Our results suggest that a value of γ a well below one is also approached for the equally oriented RCs in the limit of very low densities [cf. Fig. 2(d)]. Thus, one may see a very pronounced anomalous diffusion due to retroreflectivity at very low densities, i.e., in a regime that could be studied in the framework of kinetic theories. This is different in the EWT for equally oriented overlapping squares where the exponent γ a is only slightly different from one in the low-density limit. Thus, systems of disordered arrangements of RCs provide another paradigm for the study of anomalous transport processes via kinetic theory.

Traces of the anomalous transport due to retroreflections can be clearly seen in our experiments. Although background scattering processes cause a perturbation of the retracing events that lead to the nonexistence of the conductivity, several features can be seen in the experimental data that are also present in the classical simulation. Most remarkable is the complex behavior of the B dependence of σ xx for B < 2.0 with two underlying peaks, one being due to retroreflection and the other one due to the change of the transport from a diffusive scatter by the obstacles to a directed motion along their contour. Our study demonstrates that this complex behavior can be fully understood in terms of a classical LG model with RCs.

FIG. 1 .

 1 FIG. 1. Upper panel: Scanning electron microscopy image of the experimental system at the density n = 0.2. The inset shows a single RC and its geometry. Lower panel: Local configurations with equally oriented RCs (left) and randomly oriented RCs (right) at n = 0.2. The blue lines are trajectories, as obtained from the simulation at B = 0.

Figure 2 (

 2 a) shows the MSD at B = 0 and the different densities for the randomly and equally oriented RCs. The corresponding time dependence of the local exponent γ (t ) = d log[δr 2 (t )]/d log(t ) is displayed in Figs. 2(b) and 2(c).

FIG. 2 .

 2 FIG. 2. (a) MSDs, δr 2 (t ), at B = 0 and different densities for randomly oriented RCs ("ro") and equally oriented ones ("eo"), as obtained from the simulation. (b) and (c) show the exponent parameter γ (t ) corresponding respectively to the randomly oriented and equally oriented RCs. (d) Asymptotic exponent γ a = lim t→∞ γ (t ) as a function of n (see text).

1 4

 1 dδr 2 (t )/dt. A reduced, dimensionless diffusion coefficient D can be obtained from D = D d o v F . This diffusion coefficient is related to the conductivity via σ xx = n e e 2 d o m v F D = (0.486 × 10 -2 -1 )D .

FIG. 3 .

 3 FIG. 3. Derivative of MSD, dδr 2 (t )/dt, at n = 0.2 and different magnetic fields, including the case B = 0, for (a) the randomly oriented RCs and (b) the equally oriented RCs, as obtained from the simulation. (c) τ m as a function of 1/n . (d) Bm as a function of n . The dashed lines in (c) and (d) are fits to straight lines. (e) The conductivity σ xx ( B) at the different n for randomly oriented RCs (black solid circles) and equally oriented RCs (red open circles). The bold crosses mark the location of Bm at the different n . The solid lines are fits with power laws ∝ Bδ . The corresponding exponents δ as a function of n are shown in the inset.

FIG. 4 .

 4 FIG. 4. σ xx ( B), as obtained from experiment and simulation for equally oriented RCs ("eo"), for (a) n = 0.1, 0.2 (enlargement in the inset for n = 0.2) and (b) n = 0.3, 0.4. The experiments are done at the temperature T = 0.1 K. (c) Experimentally measured σ xx ( B) for different T and n = 0.3. (d) Temperature dependence of n 2 σ xx at B = 0, as obtained from the experiment, for the different densities.

We thank Dr. E. Neumann (HNF, Forschungszentrum Jülich) for providing scanning electron microscope pictures of the samples, as well as Dr. M. Cerchez and M. Hund for their support of the measurements at sub-Kelvin temperatures. Financial support from Heinrich Heine Universität Düsseldorf and the RENATECH program is gratefully acknowledged.


Martin Heusse
email: martin.heusse@imag.fr

Takwa Attia
email: takwa.attia@univ-grenoble-alpes.fr

Christelle Caillouet
email: christelle.caillouet@univ-cotedazur.fr

Franck Rousseau
email: franck.rousseau@imag.fr

Andrzej Duda
email: andrzej.duda@imag.fr

Capacity of a LoRaWAN Cell

Keywords: Networks → Network performance modeling, Very long-range networks, Sensor networks LoRaWAN, Packet Delivery Ratio, capacity, modeling

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

In this paper, we consider the problem of evaluating the capacity of a LoRaWAN cell. LoRaWAN deines a physical layer based on Chirp Spread Spectrum (CSS) for robust communication in the sub-GHz ISM band. The data rate depends on the Spreading Factor (SF) that characterizes the number of bits carried by a chirp. SF allows to trade data rate for range: a small value of SF (e.g., [START_REF] Waret | LoRa Throughput Analysis with Imperfect Spreading Factor Orthogonality[END_REF]) results in a higher data rate, short transmissions, and a small range, while SF12 gives the best sensitivity and the longest range, at the cost of long transmission times and worst energy consumption.

LoRa devices use a channel access method similar to unslotted ALOHA [START_REF] Kleinrock | Packet Switching in Radio Channels: Part I -Carrier Sense Multiple-Access Modes and Their Throughput-Delay Characteristics[END_REF]: a device may wake up at any instant and start transmitting a packet to the LoRa gateway without testing for on-going transmissions. The performance in terms of the Packet Delivery Ratio (PDR) and scalability strongly depend on the number of collisions arising when at least two end devices transmit packets overlapping in time. Transmissions using diferent spreading factors are quasi-orthogonalÐin case of a collision, both frame transmissions succeed if one is not several orders of magnitude stronger than the other. For collisions of transmissions with the same SF, a frame transmission succeeds if it is signiicantly stronger than the other one, a phenomenon called the capture efect 1 .

For an increasing distance to the gateway, a device needs to change SF to maintain the desired level of PDR: farther from the gateway, the nodes need to use larger SF to increase the probability of frame reception. We can thus deine SF boundaries around the gateway in form of concentric annuli in which devices share the same value of SF. In practice, the Adaptive Data Rate (ADR) commands issued by the LoRa network server give SF to use based on the measured SNR (Signal to Noise Ratio) and PDR at the gateway. Larger SFs result in more robust transmissions by using a lower data rate with similar transmission power. However, transmissions with larger SF take more time, which increases contention and collision probability.

Previous analytical studies investigated the LoRaWAN performance in terms of PDR given a number of devices around a gateway and its range [5ś10]. In this paper, we bring the following contributions:

• We propose a traic model in which all devices generate packets with the same inter-transmission times corresponding to the duty cycle limitation at the highest SF regardless of the distance to the gateway (Section 3). • We develop a simple yet precise model for PDR assuming that the following two events are dependent: successful capture during a collision and successful frame decoding despite ambient noise (Section 4). • We validate the analytical results with discrete event simulations and compare with models proposed in the literature in Section 5. • We evaluate, by simulation as well as analytically, the capacity of a LoRaWAN cell for a given density of nodes per km 2 and given SF boundaries. We irst study the case of boundaries located at distances corresponding to predetermined average SNR levels. Second, we optimize the SF boundaries so that the maximal number of nodes obtain given target PDR (Section 6). • We experimentally validate the hypothesis of Rayleigh fading for the LoRa channel model by analyzing the distribution of channel gains measured on a real-world testbed (Section 7). The important conclusion from our results is that a LoRa cell can handle a large number of devices, e.g., about 900 nodes can get at least PDR of 90% in a 10 km 2 cell or more than 1550 nodes can get 60% PDR in a 78 km 2 cell. We also show that the traic pattern that we use leads to practically no inter-SF interference (cross interference between transmissions with diferent SFs).

We start the paper with some background on LoRa in Section 2. We end it with the discussion of previous work in Section 8 and conclusions in Section 10.

LORAWAN BASICS AND ASSUMPTIONS

The physical layer of LoRaWAN has the following parameters [START_REF]SX1272/73 -860 MHz to 1020 MHz Low Power Long Range Transceiver[END_REF]: Bandwidth (BW) ś the range of the frequency sweep of CSS, between 7.8 kHz and 500 kHz. Larger bandwidth results in a higher data rate, but lowers the sensitivity. Spreading Factor (SF) SF varies between 6 (7 in practice) and 12, with SF12 resulting in the best sensitivity and range, at the cost of the lowest data rate and worst energy consumption. Decreasing SF by 1 unit roughly doubles the transmission rate and divides by 2 the transmission duration as well as energy consumption. Coding Rate (CR) ś the rate of Forward Error Correction (FEC) that improves packet error rate in presence of noise and interference. A lower coding rate results in better robustness, but increases the transmission time and energy consumption. The possible values are: 4/5, 4/6, 4/7, and 4/8. Transmission Power (P): LoRaWAN deines the following values of P for the EU 863-870 MHz band: 2, 5, 8, 11, and 14 dBm.

Table 1 presents SFj, data rate DRj, Signal-to-Noise Ratio (SNR) limit, airtime τ j , and s max , the maximum payload size. τ j denotes the transmission duration of the maximum size frame at data rate DRj. We consider in our analysis the basic LoRa CSS modulations with 125 kHz bandwidth and SF in 12, 11,. . . 7, which corresponds to data rates DRj, with j = 0, 1, . . . 5, and SF = 12 -j (see Table 1).

Assumptions and Notation

We consider a LoRaWAN cell with devices that conigure SF based on their transmission conditions, which corresponds to the annuli view presented in Figure 1. l j denotes the distance to the farthest device that uses SF{12 -j} and thus data rate DRj. Outside l 1 , there is no lower data rate available so all nodes need to use SF12. In this way, l 1 is the efective cell range to consider.

We use the Okumura-Hata model for path loss attenuation (also used in other work [12ś14]), the suburban environment variant with an antenna height of 15 m. This empirical model is slightly less favorable but more realistic than adopting an arbitrary path loss exponent as in most of previous work cited herein. We consider a gateway antenna gain of 6 dB, which compensates for a receiver noise factor of 6 dB. We neglect shadowing, since its net efect would be small in our case: it would modify the channel gain from each node without changing the general characteristics of the network [START_REF] Mahmood | Scalability Analysis of a LoRa Network under Imperfect Orthogonality[END_REF]. Finally, we consider a Rayleigh channel: a multiplicative random variable with an exponential distribution of unit mean afects the received signal power.

TRAFFIC MODEL

We propose a traic model in which all devices generate packets with the same inter-transmission times corresponding to the duty cycle limitation at the highest SF regardless of the distance to the gateway.

Such a model is more realistic than the assumptions of other work [5ś10] in which nodes saturate their share of bandwidth whatever their position in the network, which means that nodes closer to the gateway generate more traic than those farther away. In a LoRaWAN network, traic is application driven and should not depend on the position of the node in the network or channel quality. An application needs to send a given number of bytes (limited by the maximal frame size) and its operation should be transparent with respect to the position of devices that will use the same data generation pattern irrespectively of the conigured SF. Obviously, there are diferent IoT applications generating packets according to their needs, e.g., once a day, once an hour, under the constraint of the 1% duty cycle limit per band.

As devices operate with a given duty cycle and wake up to send data at constant intervals, the superposition of traic coming from a large number of devices tends to a homogeneous Poisson process and we denote its intensity with λ t that depends on the data generation pattern of IoT applications. The limiting case is when λ t is the traic intensity of devices operating at DR0 and SF12 at their maximal duty cycle and using 59 B packets, the maximum size at this rate. Our analysis will also hold for lower traic intensity coming from less demanding applications.

We assume that devices are homogeneously scattered with spatial density ρ. The number of nodes using DRj is a function of the radii deining the annulus between l j+1 and l j around the gateway (see Figure 1). We assume that l 6 = 0, so DR6 and above are not used.

The duty cycle limitation is for each frequency band: LoRa devices have to limit their occupation of each frequency band to 1% of time with 3 to 5 frequency channels in each band. In Europe, for instance, there are 3 channels in band h1.4 and 5 in h1.3, so if we assume 3 channels, devices can use each of them at 0.33% duty cycle. The airtime of maximum size packets at DR0 corresponds to τ 0 = 2.466 s, so they can be sent every 739.8 s on the same channel (or every 246.6 s on diferent channels) to achieve 0.33% duty cycle ( 1 300 ) per channel. Thus,

λ t = 1 739.8 s . (1)
For SF = 12 -j, and n j , the number of nodes contending with each other in annulus j, ofered traic is

v j = n j τ j λ t (2)
in Erlang, where τ j is the transmission duration at data rate DRj, and

n j = π ρ (l 2 j -l 2 j+1 ). (3)

MODEL FOR SUCCESSFUL PACKET DELIVERY

We present below the model to compute PDR based on the notation in Table 2. First, we derive the expression for PDR under the assumption that collisions and successful overcoming the ambient noise are independent. Then, we relax this assumption and derive the inal expression for PDR. In a LoRaWAN cell, a given frame is not correctly received for two reasons (and maybe both): i) SNR is below the reception threshold and/or ii) a collision occurs and the signal is not strong enough relatively to the interference.

Successful frame reception, no collisions

Provided that there is no collision, a frame transmission succeeds as long as SNR at the receiver for this transmission is above q j , the minimum SNR for the corresponding spreading factor [START_REF] Goursaud | Dedicated Networks for IoT: PHY/MAC State of the Art and Challenges[END_REF]. The received signal power depends on the distance and Rayleigh fading, whereas the noise power is the constant thermal noise for a 125 kHz-wide band: N = -123 dBm (-174 dBm per Hz).

Thus, the probability of successful transmission from distance d at data rate DRj is:

H = exp - Nq j Pд(d ) , (4)
where д(d ) is the average channel gain at distance d [START_REF] Georgiou | Low Power Wide Area Network Analysis: Can LoRa Scale[END_REF].

Successful frame reception, collisions

For a given transmission attempt, provided that it has already met the signal strength reception condition (with probability H ), we distinguish between three cases:

-Case 1 (0_collision) ś if the transmission does not overlap any other frame, it succeeds; -Case 2 (1_collision) ś if a single frame transmitted on the same SF interferes with the transmission, then the frame is captured, if it is received with a power margin of at least γ = 6 dB against the colliding frame [15]; -Case 3 (2+collisions) ś if two frames or more overlap the transmitted one, we simply deem it lost.

For a tagged transmission, the probability of Case 1 is:

Q 1 = exp -2v j , (5)
as there should be no other transmission event during time 2τ j to avoid overlap.

Case 2 (1_collision) happens when a single transmission occurs during time 2τ j with probability 2 v j exp(-2v j ) and the receiver successfully captures the frame. If we neglect the variability of д(d ) among the nodes using the same SF, then the received power in each annulus follows an exponential distribution. Successful frame capture occurs for a diference of γ in received power (typically, γ = 6 dB, which corresponds to a factor of 4). Since the probability that an exponential random variable is γ times greater than another one is (γ + 1) -1 , the probability of successful frame capture is just (γ + 1) -1 . Thus, the probability of success in Case 2 is the following:

Q 2 = 2 γ + 1 v j exp(-2v j ). (6)
For Case 3 (2+collisions), we can observe that collisions with successful packet capture against a single frame are already relatively rare, so meeting the capture condition for two interferers, or more, is efectively quite unlikely. Moreover, the co-occurrence of three transmissions is rare for manageable traic intensity.

Thus, concurrent traic impacts the successful packet reception with probability Q:

Q = Q 1 + Q 2 = (1 + 2 γ + 1 v j ) exp(-2v j ). (7)
Finally, combining Eqs. 4 and 7 yields the probability of successful packet reception under the independence assumption:

PDR I = H × Q. ( 8
)

Successful frame reception, dependent reception and capture

In the model above, and in all the cited papers, the probability of capture in case of a collision and the probability of successfully overcoming the ambient noise are considered independent. However, it is not the case, as a strong signal at the receiver is at the same time more likely to overcome thermal noise and the interference from colliding transmissions.

In case of a single collision, the probability of successful frame reception at DRj is thus:

PDR 1 = Pr д r > Nq j Pд(d ) ∩ д r > γд c (9) = Pr д r > max Nq j Pд(d ) , γд c , (10)
where д r and д c are the exponentially distributed Rayleigh channel gains for the considered frame and the colliding frame, д(d ) is the average path loss gain for both transmissions (thus assuming that both nodes are at distance d from the gateway), and γ is the required receiving power gap for successful capture (typically, 6 dB).

If we denote д t = N q j Pд(d ) , the threshold gain for successful reception against the thermal noise, the probability of successful frame reception becomes:

PDR 1 = +∞ 0 e -y +∞ max(д t ,γ y ) e -x dxdy (11) = e -д t γ + 1 1 + γ (1 -e -д t γ ) ( 12
)
by integral of the probability density function over the successful reception interval. Finally, the probability of successful packet reception under the dependence assumption is the following:

PDR D = HQ 1 + 2v j exp(-2v j ) PDR 1 . ( 13
)

MODEL VALIDATION

In this section, we compare our models and the model by Georgiou and Raza [START_REF] Georgiou | Low Power Wide Area Network Analysis: Can LoRa Scale[END_REF] with simulations.

Simulations

We have developed a discrete event simulator in Python for a single LoRa frequency channel used by a large number of nodes that takes into account Rayleigh fading, capture, as well as co-SF and inter-SF interference 2 . It simulates the transmission of a large number of frames (e.g., 10 7 ) generated by as many Poisson processes as SF zones, with intensities derived from the SF zone area, node density, and λ t . The transmission distance distribution in each SF zone matches the node distribution. When frame transmissions overlap, the simulator computes co-SF and inter-SF interference to evaluate the efect of collisions, capture, and inter-SF jamming, using the theoretical [START_REF] Goursaud | Dedicated Networks for IoT: PHY/MAC State of the Art and Challenges[END_REF] or experimental levels [START_REF] Croce | Impact of LoRa Imperfect Orthogonality: Analysis of Link-Level Performance[END_REF]. Based on the statistics of the number of correctly received frames, it computes the conditional PDR, given the distance to the gateway, with high precision (the conidence intervals are less than one percentage point). Compared to more complete simulators which include all protocol layers, focusing solely on contention allows to simulate hundred of thousand of frame transmission in a few minutes. The simulator and analytical results corroborate each other. Simulation brings to light the efect of signiicant path loss diferences between nodes in the same SF zone as well as the impact of inter-SF interference or the salience of the capture efect when multiple frames collide. Conversely, the analytical results are more suitable for SF boundary optimization that we develop in Section 6.

Comparisons

In this section, to compare our results with the familiar curves for the unslotted ALOHA model [START_REF] Kleinrock | Packet Switching in Radio Channels: Part I -Carrier Sense Multiple-Access Modes and Their Throughput-Delay Characteristics[END_REF], we consider a variable number of nodes n 0 , all at the same distance from the gateway. They send frames using SF12 with aggregate traic intensity v 0 = n 0 τ 0 λ t and we compute the channel utilization U = PDR × v 0 . We assume transmission power of P = 14 dBm. Figures 2 and3 present the channel utilization based on PDR I , PDR D , and PDR G (the model by Georgiou and Raza [START_REF] Georgiou | Low Power Wide Area Network Analysis: Can LoRa Scale[END_REF]) as well as the corresponding value for unslotted ALOHA at SF12, for the distances l 0 = 2500 m and l 0 = 7500 m, respectively.

As for Aloha, we consider that a collision happens if another frame starts transmitting within τ j before and after the start of a transmission. It is not the case in the model by Georgiou and Raza [START_REF] Georgiou | Low Power Wide Area Network Analysis: Can LoRa Scale[END_REF][START_REF] Georgiou | Coverage scalability analysis of multi-cell LoRa networks[END_REF]. Their model rather corresponds to slotted Aloha or to continuous transmissions as in CDMA networks. Interestingly, probability Q in this paper matches approximation Q 1 by Georgiou and Raza [START_REF] Georgiou | Low Power Wide Area Network Analysis: Can LoRa Scale[END_REF] with the precision of a few percent provided that we double the intensity of the interfering traic in their expressions to account for the 2τ j time during which at most one frame can start transmitting for no overlap. Consequently, in Figure 2 where Q is the dominant cause of losses, the corrected PDR G curve closely matches our simulations and analytical expressions.

Moreover, in this short range case, the curves for PDR I and PDR D overlap, which is not the case for the longer range. We observe that for light ofered traic, all models roughly give the same results whereas diferences become visible for higher ofered traic. Compared to all the models that take into account the capture efect, the unslotted ALOHA model [START_REF] Kleinrock | Packet Switching in Radio Channels: Part I -Carrier Sense Multiple-Access Modes and Their Throughput-Delay Characteristics[END_REF] clearly stands out: physical capture increases the capacity of the network and channel utilization.

Simulations allow to analyze the efect of capture in case of multiple colliding frames: the crosses show the utilization when three overlapping frames are all dropped, whereas the rounds correspond to a more ine-grained analysis in which any frame can be received as long as it is received with a γ margin over the sum of the power of the concurrent transmissions. For channel load below 1 Erlang, the assumption of a single interfering frame causes only a slight loss of precision.

NUMERICAL RESULTS

In this section, we present the results of the model and use it to optimize the SF boundaries. In this way, we obtain the number of nodes that a single gateway can serve, for various ranges and acceptable PDR. We compare the model with simulations that consider the efects of near-far co-SF and inter-SF interference. We highlight the necessity of power control in the SF7 zone and observe that inter-SF interference is moderate with our traic model.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • x x x x
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x o x Simulation w/o inter-SF Simulation w/ inter-SF Figure 6: Large cell radius (l 1 = 6.30 km), density of 5 nodes per km 2 , SNR-based SF boundaries. For lower density, nodes using small SF face little contention, but collisions still impact PDR for nodes located further away. 443 nodes within 5.3 km of the gateway experience PDR above 60%.

We consider three scenarios with three diferent radio ranges of 2.82, 5.23, and 7.36 km (small, medium, and large cell radius), respectively. For instance, in a medium range LoRaWAN cell, we would change SF each time the probability H drops below H target = 90%, which is the value of H at a range of 5.23 km. Solving Equation 4for each j gives the range for each data rate shown in Table 3. For completeness, the table includes the l 0 distances at which H drops below the target value at SF12, but the nodes beyond this radius cannot change to a lower SF and continue to use SF12. In the medium range case, the covered area is 86 km 2 , whereas the same gateway could cover an area of 170 km 2 only at the cost of setting H target = 70% (or a cell radius of 7.36 km).

Figures 4,5, and 6 present the values of H , v j , and PDR D for the three cell sizes. For the small cell size, we consider the density of 90 nodes per km 2 , 20 nodes per km 2 for the medium cell, and 5 per km 2 for the large one. The x-axis is the index of a nodeÐwe order the nodes by the distance to the gateway, e.g., the 1196 th node is at the distance of 2.07 km in the small cell.

The rounds and crosses represent the results obtained by simulation. The rounds give the average PDR for the transmissions emanating from the distance range on the x axis under the assumption of perfect inter-SF orthogonality. The crosses correspond to simulation results that take into account inter-SF interference [START_REF] Goursaud | Dedicated Networks for IoT: PHY/MAC State of the Art and Challenges[END_REF]. In the SF7 zone and especially for high node density, the simulation results diverge from the model because simulation takes into account the near-far problem: in case of collisions, nodes located close to the gateway always dominate the nodes further away.

All igures exhibit the same trend. Close to the gateway, for small SFs, there is little contention (v j ≪ 1) and PDR D is satisfactory. We can also observe that the smaller the cell, the more nodes it can handle, because more nodes can use a lower SF and have a high H probability. Shorter transmissions then reduce the channel load, which leaves room for incorporating many nodes before PDR D begins to drop due to collisions.

Conversely, in a larger cell, two adverse efects combine for the nodes located farther away. On the one hand, transmissions with higher SF take more airtime, resulting in more collisions. On the other hand, the annuli areas are larger going outward of the coverage disk, thus including more nodes at the same SF, which Node #: 0 1.22 1.52 1.67 1.78 l j (km):

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x Node #: 0 2.43 3.01 3.33 3.58 l j (km):

• • • • • • • • • • •• •• • •• •• •• ••• ••• •• • ••• ••• ••• ••• ••• ••••• ••• •• •••• •••• • ••••••••• x x
x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x also leads to more collisions, and the overall node capacity drops to several hundreds of nodes when we want to get 60% PDR. The next section gives numerical examples of the number of nodes that a gateway can serve with optimized boundaries.

PDR based boundaries and cell capacity

We start by inding the SF boundaries to meet a given target PDR D : we use the bisection method to position the SF7 boundary at the distance where PDR D matches the target value. In this way, the annulus includes the maximal number of nodes that obtain at least target PDR D . Then, we change to SF8 and increase its boundary to reach target PDR D , and so on for other SF boundaries.

• • • • • • • •• • • • • • •• • • • • •• • •• • • • • • •• •• •• •• •• • ••• •• •• •• •• • • ••• ••• ••••••••• • ••••••••• x x x x x x x
•••• •••• ••• ••• ••• •• •• • • ••• •••• ••• •••••• ••••••••••••••• •••••••• • •••
x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x Figures 7, 8, 9, 10, 11, and 12 present PDR D for various densities of nodes and two target PDR D values of 60% and 90%. In a small cell, probability H is high, so nodes can easily attain PDR D of 90% (see Figure 7). Compared to Figure 4, we can observe that the frontiers for SF7 and SF8 are farther away by 50 m and 100 m respectively, but those for SF9, SF10, SF11 are closer to the gateway (by as much as 0.62 km for SF11), which results in a smaller cell with the total of 908 nodes with PDR D of 90%. If we lower target PDR D to 60%, the resulting cell is much larger with a signiicantly larger number of nodes beneiting from target PDR D : 3648.

•• •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • •• • •• •• •• • ••• •• •• •• •• •• • ••• ••• •••• • ••• ••••• •• xx xx x x
The model leaves aside inter-SF interference and near-far unfairness. In most cases, these efects have little impact, except for heavily loaded small cells (see Figure 8) that we discuss in more detail in Sections 6.3 and 6.4. Table 4 summarizes the capacity limitsÐthe maximal number of nodes that beneit from a given level of target PDR D . For small cells, optimizing the SF allocation has more efect on capacity, as the optimization process has more freedom to spread nodes between the SF annuli. In this case, there is more than a threefold diference between the maximal number of nodes that a single gateway can handle, with vs. without SF boundary optimization (e.g., 3648 vs. 1196 for 60% PDR).

When the cell range increases, optimization of SF allocation is more constrained and the maximal number of nodes is slightly less sensitive to the tuning of the SF boundaries. Moreover, the maximal number of nodes considerably increases when we lower target PDR D (e.g., to 60%). This observation suggests that we can consider other ways of achieving high probability of successful packet delivery: nodes may repeat packets several times [START_REF] Birk | Judicious Use of Redundant Transmissions in Multichannel ALOHA Networks with Deadlines[END_REF]. The resulting packet delivery probability is: Pr(packet delivery) = 1 -(1 -PDR D ) n , where n is the number of packet repetitions. For instance, for PDR D = 60% and n = 3, Pr(packet delivery) = 0.94. Obviously, packet repetitions increase energy consumption; however, changing SF to larger values has a similar efectÐgoing from SFj to SF{j + 1} doubles the transmission duration. So, repeating 4 times a packet sent at SF7 has approximatively the same duration as a single transmission at SF9: the preamble overhead is the same because the preamble duration grows with SF and the only diference is that each symbol carries one more bit at each SF step. Using inter-frame error correction is even more efective, albeit more complex to implement [START_REF] Coutaud | Fragmentation and forward error correction for LoRaWAN small MTU networks[END_REF].

In these calculations, we have considered an application that would saturate a band with a 1% duty cycle limit and wide enough to include 3 channels. In this case, an application generating one packet every 246.6 s (4 mn and 7 s) on diferent channels reaches the duty cycle limit at SF12. Applications generating less traic can either:

• send less packets, which saves energy for the individual nodes and reduces contention or, alternatively, increases the network capacity in terms of the number of nodes; • repeat the data packets several times (packet repetition [START_REF] Birk | Judicious Use of Redundant Transmissions in Multichannel ALOHA Networks with Deadlines[END_REF]) to improve the desired data delivery probability signiicantly.

Power control in DR5 (SF7) zone

In our model, all nodes in the same annulus experience a similar average channel gain. As the diference in the coding gain when changing from SFj to SF{j + 1} is approximately 2.5 dB, the channel Node #: 0 2.43 3.01 3.33 3.58 l j (km):

• • • • •••••• • • •• • •• • •• •• •• •• ••• • •• •• ••• •• ••• ••• ••••• ••• •• ••••• ••• • ••••••••• x x x x
x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x gains for a given SF are spread over this interval (i.e., 2.5 dB), which is a quite reasonable assumption for all annuli except for the inner disk (DR5/SF7).

In this disk of radius l 5 , channel gains are heterogeneous, so power control is necessary for fair access, as some nodes may have a markedly shorter distance to cross compared to others. Thus, a proper ADR scheme would make a node close to the edge of the SF7 zone use P = 14 dBm, as for higher SFs. It would then attribute the next lower transmission power to the nodes enjoying between 2 dB and 4 dB of average SNR advantage compared to the edge, and so on by steps of 2 dB, until P = 2 dBm. This approach has two beneits, as illustrated in Figure 13 (to compare with Figure 8 presenting the results without power control): i) power-controlled nodes closer to the gateway are less harmful to farther nodes with respect to inter-SF interference and ii) energy consumption is marginally lower [START_REF] Ochoa | Evaluating LoRa Energy Eiciency for Adaptive Networks: From Star to Mesh Topologies[END_REF].

Co-SF and inter-SF interference

Experimental measurements by Mikhaylov et al. [START_REF] Mikhaylov | On LoRaWAN Scalability: Empirical Evaluation of Susceptibility to Inter-Network Interference[END_REF] showed that physical capture in presence of co-SF interference may happen well below 6 dB of the power margin. Compared to our model, this efect would tend to increase the actual packet delivery ratio. Some authors [START_REF] Croce | Impact of LoRa Imperfect Orthogonality: Analysis of Link-Level Performance[END_REF][START_REF] Mikhaylov | On LoRaWAN Scalability: Empirical Evaluation of Susceptibility to Inter-Network Interference[END_REF] also question the theoretical levels of inter-SF interference: in particular, low SFs may be jammed by other signals even when the power margin is of the order of -6 dB, whereas the expected values would be at least 10 dB lower [START_REF] Goursaud | Dedicated Networks for IoT: PHY/MAC State of the Art and Challenges[END_REF]. This aspect is not overly critical: irst, interference from higher SFs comes from nodes located farther away, so they face greater attenuation and thus, they are relatively harmless. Second, for the optimized SF boundaries, channel utilization is much reduced for higher SFs, so that nodes in the lower SF annuli face lower levels of accumulated interference from farther nodes compared to the SNR-based SF allocation for instance.

In fact, using the less favorable inter-SF rejection levels found by Croce et al. [START_REF] Croce | Impact of LoRa Imperfect Orthogonality: Analysis of Link-Level Performance[END_REF], inter-SF interference is non-negligible in a dense cell as we can see in Figure 14. However, since inter-SF interference mostly comes from transmissions of nodes located near the gateway, power control greatly helps (see Figure 15). In all cases, inter-SF interference is less of a problem with our homogeneous Node #: 0 2.43 3.01 3.33 3.58 l j (km):

• • • • • • • • • •• •• •• •• •• •• ••• ••• •• • •• ••• •• •• ••• ••• ••••• ••• •• •••• •••• • ••••••••• x x
x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x

• • • •••••• • • •• •• •• •• •• •• •• •• • • ••• •• ••• ••• •• •• ••••• •••• • •••• •••• • ••••••••• x x x x
x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x o x Simulation w/o inter-SF Simulation w/ inter-SF Figure 15: Small cell radius, density of 90 nodes per km 2 , SF boundaries for target PDR of 60%. This is for the case of low inter-SF rejection and using power control in SF7 zone. Compared to Figure 14, PDR is more even among SF7 nodes and there is less (although still noticeable) inter-SF interference for SF8 and above.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x
Figure 16: Medium cell radius, density of 20 nodes per km 2 , SF boundaries for target SNR of 90%. In this simulation, all nodes use 1% duty cycle, which corresponds to the traic model in most of the litterature [START_REF] Georgiou | Low Power Wide Area Network Analysis: Can LoRa Scale[END_REF][START_REF] Li | 2D Time-Frequency Interference Modelling Using Stochastic Geometry for Performance Evaluation in Low-Power Wide-Area Networks[END_REF][START_REF] Mahmood | Scalability Analysis of a LoRa Network under Imperfect Orthogonality[END_REF][START_REF] Georgiou | Coverage scalability analysis of multi-cell LoRa networks[END_REF][START_REF] Hoeller | Analysis and Performance Optimization of LoRa Networks With Time and Antenna Diversity[END_REF]. Compared to Figure 5, a lot more packets are lost due to collisions and inter-SF interference still erase many transmissions.

traic pattern, whereas it is exacerbatedÐas well as collisionsÐ when all nodes transmit using the same duty cycle as in the model by Georgiou and Raza presented in Figure 16. Figure 18: Measured RSSI distribution for a 2.4 km Lo-RaWAN link, from an indoor sender, transmission power 14dBm, SF11. The distribution is truncated to the left because the reception power is then below the sensibility.

EXPERIMENTAL ASSESSMENT OF THE CHANNEL MODEL

The theoretical model presented in this paper as well as those in related work assume a Rayleigh fading channel. This assumption is quite usual, except in rare cases [START_REF] Markkula | Simulating LoRaWAN: On Importance of Inter Spreading Factor Interference and Collision Efect[END_REF], but we draw attention to the fact that it was originally a model for mobile nodes. In the case of LPWAN networks, we consider mostly static nodes so that it would also be reasonable to assume less variability of the channel gain.

To get insight on how the channel gain behaves in the real-world, we have run long term experiments using an STMicroelectronics mote [START_REF]B-L072Z-LRWAN1[END_REF] and several gateways of The Things Network (TTN), a public LoRa network [START_REF] Giezeman | The Things Network[END_REF]. We have generated traic and gathered statistics over 3 months to obtain the RSSI distributions of which two examples are shown in Figures 17 and18. We send at SF11 with transmission power of 14 dBm. The transmissions alternate between 6 diferent frequency channels. Both histograms are for non-line of sight links in a suburban environment. The RSSI distribution of the long link (see Figure 17a) follows well the expected distribution of gains for a Rayleigh channel described in Section 2.1, which appears in Figure 17bÐdeep fades of -10 dB below the average gain (0 dB in Figure 17b) are common.

For the longer link, we found that the use of several frequency channels explains the observed variability: the nodes are not mobile but frequency diversity has a similar efect, which is reminiscent of the reason why GSM networks use frequency hopping whereas GSM and LoRaWAN have similar transmission range, frequency channel separation, and carrier frequency.

The link in Figure 18 shows more attenuation because the sender is indoor, although the link in shorter. in this case, the variability would also randomize the outcome of the capture efect should a collision occur. Note that the sensitivity of approximately -140 dBm for SF11 truncates the distribution because the frames experiencing high fading are lost.

Finally, the measures conirm that the channel gains are always variable and that even for nodes with similar average channel gains, it would be very common to witness reception powers with a difference of 6 dB. In this case, the capture efect allows the gateway to receive one of the colliding frames.

The Rayleigh channels assumption discussed in this section is important because a channel with less variability would dramatically change the results, with a longer range and also little inter-SF impact [START_REF] Markkula | Simulating LoRaWAN: On Importance of Inter Spreading Factor Interference and Collision Efect[END_REF]. Experimentally, our inding is that the Rayleigh channel assumption is reasonable at least for links longer than a couple of km.

DISCUSSION AND RELATED WORK 8.1 Collision model

Regarding collisions, we point out that solely meeting the capture condition (with probability Q in Equation 7) is not suicient to ensure successful reception, which is in contradiction with the model by Waret et al. [START_REF] Waret | LoRa Throughput Analysis with Imperfect Spreading Factor Orthogonality[END_REF]. We think that attenuation must always be taken into account, so PDR always integrates probability H (Equation 4) in this paper.

Also, our model leads to a more conservative estimate of the probability of packet capture than those proposing to spread the energy of the overlapping interferer over the entire considered frame [START_REF] Li | 2D Time-Frequency Interference Modelling Using Stochastic Geometry for Performance Evaluation in Low-Power Wide-Area Networks[END_REF][START_REF] Magrin | Performance Evaluation of LoRa Networks in a Smart City Scenario[END_REF][START_REF] Markkula | Simulating LoRaWAN: On Importance of Inter Spreading Factor Interference and Collision Efect[END_REF].

We do not take into account the collision timing, although it has some efect [START_REF] Rahmadhani | When LoRaWAN Frames Collide[END_REF] on the frame capture. We leave this reinement for future work.

Traic models

Many recent papers focusing on the LoRaWAN radio capacity assume that all nodes completely use their allowed duty cycle regardless of their SF [START_REF] Georgiou | Low Power Wide Area Network Analysis: Can LoRa Scale[END_REF][START_REF] Li | 2D Time-Frequency Interference Modelling Using Stochastic Geometry for Performance Evaluation in Low-Power Wide-Area Networks[END_REF][START_REF] Mahmood | Scalability Analysis of a LoRa Network under Imperfect Orthogonality[END_REF][START_REF] Georgiou | Coverage scalability analysis of multi-cell LoRa networks[END_REF][START_REF] Hoeller | Analysis and Performance Optimization of LoRa Networks With Time and Antenna Diversity[END_REF]. This assumption means sending much more packets when a node is closer to the gateway, which is unrealistic: let us assume that an IoT application needs to send a packet every hour and the node is at the distance resulting in SF8. If the transmission conditions change and it can switch to SF7, it will not start sending packets every half an hour because of the possibility of doing so.

Under the assumption that all nodes use their allowed duty cycle like in Figure 16, if the metric for the capacity is the overall throughput, then the models give optimistic results (i.e., high throughput) because all nodes are able to ill up all available quasi-orthogonal sub-channels [START_REF] Li | 2D Time-Frequency Interference Modelling Using Stochastic Geometry for Performance Evaluation in Low-Power Wide-Area Networks[END_REF]. On the contrary, if the metric to express the capacity is the number of nodes [START_REF] Georgiou | Low Power Wide Area Network Analysis: Can LoRa Scale[END_REF][START_REF] Mahmood | Scalability Analysis of a LoRa Network under Imperfect Orthogonality[END_REF][START_REF] Georgiou | Coverage scalability analysis of multi-cell LoRa networks[END_REF][START_REF] Hoeller | Analysis and Performance Optimization of LoRa Networks With Time and Antenna Diversity[END_REF], the models result in dramatically small node capacity (e.g., 20 nodes for PDR of 80%) because all nodes saturate their duty cycle, so there is little beneit in moving a node to a smaller SF, since it immediately transmits much more traic and thus creates similar collision levels.

ACKNOWLEDGEMENTS

This work has been partially supported by the ANR projects: PERSY-VAL-Lab under contract ANR-11-LABX-0025-01 and DiNS under contract ANR-19-CE25-0009-01.

CONCLUSION

In this paper, we propose a model for the capacity of a LoRaWAN cell that takes into account collisions and the capture efect. It allows to assess the number of nodes that can share a single LoRaWAN cell for a given node density, cell range, traic generation rate, and target PDR. We have validated the model with discrete event simulations.

Our results show that a LoRaWAN cell can handle thousands of devices. We also consider the efect of inter-SF interference, which turns out to be relatively small: interference from higher SFs comes from nodes located farther away, so they face greater attenuation and thus, they do not interfere with lower SF nodes. Moreover, power control can greatly help to avoid the jamming of distant nodes by the near ones.

Our hypothesis of similar traic intensity in the cell puts under the spotlight the problem of suitable SF allocation. In particular, for devices close to the gateway, it is critical to distribute as much traic as possible on lower SFs. For farther devices, there is a balance to ind: they need to use larger SFs to reduce the impact of fading and at the same time, use the smallest possible SFs to avoid collisions.

Finally, a salient question left for future work is how to practically allocate SFs: gateways have limited downlink transmission opportunities [START_REF] Vincenzo | Improving Downlink Scalability in LoRaWAN[END_REF], so that a practical ADR control method would have to take this limitation into account. Besides, a real network also has to consider the presence of multiple gateways and hence, macro-diversity reception, a central feature of LoRaWAN, which we plan to analyze in the future.

Figure 2 :Figure 3 :

 23 Figure 2: Channel utilization based on the values of PDR I , PDR D , and PDR G (the model by Georgiou and Raza[START_REF] Georgiou | Low Power Wide Area Network Analysis: Can LoRa Scale[END_REF], with or without correction) as well as the corresponding value for unslotted ALOHA, when all nodes are 2.5 km away from the gateway. Points correspond to the simulation values: o -success for γ margin over the sum of the powers of all colliding frames; x -success conditioned by at most one colliding frame and γ margin. 0.0 0.5 1.0 1.5 2.0 0.00 0.05 0.10 0.15 0.20 0.25 Offered traffic (Erlang)

Figure 5 :

 5 Figure 5: Medium cell radius (l 1 = 4.54 km), density of 20 nodes per km 2 , SNR-based SF boundaries. Losses due to attenuation and fading start to weigh on PDR. 950 nodes within 3.9 km of the gateway experience PDR above 60%

Figure 7 :

 7 Figure 7: Small cell radius, density of 90 nodes per km 2 , SF boundaries for target PDR of 90%.

x x x x x x x x x x o x Figure 8 :

 x8 Figure 8: Small cell radius, density of 90 nodes per km 2 , SF boundaries for target PDR of 60%.

Figure 9 :

 9 Figure 9: Medium cell radius, density of 20 nodes per km 2 , SF boundaries for target PDR of 90%.

Figure 11 :

 11 Figure 11: Large cell radius, density of 5 nodes per km 2 , SF boundaries for target PDR of 90%.

Figure 12 :

 12 Figure 12: Large cell radius, density of 5 nodes per km 2 , SF boundaries for target PDR of 60%.

Figure 13 :

 13 Figure13: Small cell radius, density of 90 nodes per km 2 , SF boundaries for target PDR of 60%, power control in DR5 zone.

Figure 14 :

 14 Figure14: Small cell radius, density of 90 nodes per km 2 , SF boundaries for target PDR of 60%. Low inter-SF rejection[START_REF] Croce | Impact of LoRa Imperfect Orthogonality: Analysis of Link-Level Performance[END_REF].

Figure 17 :

 17 Figure 17: (a) Measured RSSI distribution for a 4.4 km LoRaWAN link from an outdoor sender, transmission power 14dBm, SF11. (b) Expected fading gain (in dB) distribution for a Rayleigh channel.

Table 1 :

 1 LoRa parameters for BW of 125 kHz.

	SFj SNR limit (q j ) Airtime (τ j )	Bit rate DRj s max
	7	-7.5 dB	102.7 ms DR5: 5469 b/s 230 B
	8	-10 dB	184.8 ms DR4: 3125 b/s 230 B
	9	-12.5 dB	328.7 ms DR3: 1758 b/s 123 B
	10	-15 dB	616.5 ms DR2: 977 b/s	59 B
	11	-17.5 dB	1315 ms DR1: 537 b/s	59 B
	12	-20 dB	2466 ms DR0: 293 b/s	59 B
			SF7 SF8	l 5 l 4	l 3 l 2	l 1
			SF10	SF9	l 0
				SF11

SF12

Figure 1: Annuli of SF allocation around the gateway

Table 2 :

 2 Notation

	Spatial density of nodes	ρ
	Packet generation intensity	λ t
	Frame transmission duration at data rate DRj	τ j
	Distance of farthest node using DRj	l j
	Ofered traic (in Erlang) at DRj	v j
	Average channel gain at distance d	д(d )
	SNR threshold for DRj	q j
	Power gap for successful capture (typically 6 dB)	γ
	Transmission power, In-band noise power	P, N
	Transmission success probability taking into account
	attenuation, fading, thermal noise	H
	Transmission success probability due to collisions	Q

Table 3 :

 3 SNR-based SF boundaries (km)

				S F 7	SF8	SF9 SF10 SF11 SF12
				l 5	l 4	l 3	l 2	l 1	l 0
		H target = 99%	1.18	1.43	1.72	2.07 2.41	2.82
		H target = 90%	2.23	2.68	3.23	3.89 4.54	5.23
		H target = 70%	3.09	3.72	4.48	5.40 6.30	7.36
		τ j (ms)	102.7 184.8 328.7 616.5 1315 2466
	%	20 40 60 80 100	H (attenuation, fading) PDR D v j (Erlang)
		0	SF7	SF8	SF9	SF10	SF11
	Node #:	0	500		1000	1500
	l j (km):	0	1.18 1.43		1.72	2.07	2.41

  2 , SNR-based SF boundaries. Losses due to attenuation and fading start to weigh on PDR. 950 nodes within 3.9 km of the gateway experience PDR above 60%6.1 SNR-based SF boundariesWe start with the SF allocation and the annuli boundaries based on the measured SNR at the gateway, which corresponds to the value of probability H . The allocation relects the operation of Adaptive Data Rate (ADR) in which the network server issues ADR commands to nodes indicating which SF to use.

	%	20 40 60 80 100		H (attenuation, fading) PDR D v j (Erlang)
		0		SF7	SF8	SF9	SF10	SF11
	Node #:	0	100	200	300	400	500	600
	l j (km):	0		3.09 3.72	4.48	5.4		6.3

  2 , SF boundaries for target PDR of 90%.

	%	20 40 60 80 100		H (attenuation, fading) PDR D v j (Erlang)
		0		SF7	SF8	SF9	SF10	SF11
	Node #:	0	100	200	300	400	500	600
	l j (km):	0		3.34	4.04	4.76	5.46	5.93

Table 4 :

 4 Maximal number of nodes with guaranteed PDR values of 60% and 90% for various node densities

	Nodes per km 2 Target PDR	Number of served nodes	Coverage radius l 1 (km)
	90	90% 60%	908 3648	1.79 3.59
	20	90% 60%	510 1563	2.85 4.99
	5	90% 60%	198 553	3.56 5.94

The capture efect refers to the capacity of correctly receiving a signiicant fraction of colliding frames [2ś4].

https://gricad-gitlab.univ-grenoble-alpes.fr/heussem/lora-simple-discrete-eventsimulator