Alessandro Ceci
email: aceci@kth.se

Romain Gojon

Mihai Mihaescu

Large Eddy Simulations for Indirect Combustion Noise Assessment in a Nozzle Guide Vane Passage

Keywords: CFD, LES, Aeroacoustics, Indirect combustion noise Nomenclature

The combustion noise in aero-engines is known to originate from two different sources. First, the unsteady heat release in the combustion chamber generates the direct combustion noise. Second, hot and cold spots of air generated by the combustion process are convected and accelerated by the turbine stages and give rise to the so-called indirect combustion noise. The present work targets, by using a numerical approach, the generation mechanism of indirect combustion noise for a simplified geometry of a turbine stator passage. Periodic temperature fluctuations are imposed at the inlet, permitting to simulate hot and cold packets of air coming from the unsteady combustion. Three-dimensional Large Eddy Simulation (LES) calculations are conducted for transonic operating conditions to evaluate the blade acoustic response to the forced temperature perturbations at the inlet plane. Transonic conditions are characterized by trailing edge expansion waves and shocks. It is notably shown that their movement can be excited if disturbances with a particular frequency are injected in the domain.

A

amplitude of the entropy wave f 1 forcing frequency of the entropy wave z cell width in z-direction L z total mesh depth in z-direction

C

Introduction

Jet noise suppression technologies for aero-engines are currently available and implemented. Thus, the research attention was directed recently towards the relative importance of other noise sources in jet engines, e.g. indirect combustion (or entropy) noise. Its source was early investigated in the works of Candel [START_REF] Candel | Analytical Studies of Some Acoustic Problems of Jet PhD thesis[END_REF] and Marble & Candel [START_REF] Marble | Acoustic disturbances from gas non uniformities convected through a nozzle[END_REF]. They proposed an analytical model for the convection of non-uniform temperature regions through a nozzle in several configurations. In 1977, Cumpsty and Marble [START_REF] Cumpsty | The interaction of entropy fluctuations with turbine blade rows; a mechanism of turbojet engine noise[END_REF] proposed an analytical method, based on the actuator disk theory, for the evaluation of entropy noise through several turbine stages. Leyko et al. [START_REF] Leyko | Simulation and modeling of the waves transmission and generation in a stator blade row in a combustion-noise framework[END_REF] analysed the wave generation and transmission mechanisms from a numerical perspective. In their work, they assessed the range of validity of the compact assumption for a stator blade row, comparing the analytical results with simulation data. In this paper, a transonic off-design condition of a nozzle guide vane (NGV) is analysed. The aim of the present work is to assess the indirect combustion noise generation mechanism using 3D LES calculations. A simplified topology will be adopted in the simulations in order to have a deeper physical comprehension of the phenomenon. The derived acoustic reflection and transmission coefficients will be compared with the ones predicted by the analytical model of Cumpsty and Marble [START_REF] Cumpsty | The interaction of entropy fluctuations with turbine blade rows; a mechanism of turbojet engine noise[END_REF]. Furthermore, singular flow features arising from transonic operating conditions are captured by the analysis of the frequency content of two-dimensional slices sampled during the simulation. The first part of the analysis will focus on the validation of the numerical solver against experimental results available for the baseline case without incoming disturbances. The second part, instead, will move towards the indirect combustion noise mechanism assessment.

Geometry and Mesh

The geometry adopted in the present study consists in a simplified topology of a real NGV cascade, for which experimental data were acquired by Yasa et al. [START_REF] Yasa | Performance of a nozzle guide vane in subsonic and transonic regimes tested in an annular sector[END_REF] in a transonic test tunnel facility. The annular sector cascade was designed by Siemens Industrial Turbomachinery AB in 1998 and the guide vanes have identical geometry as the NGVs used in the first stage High Pressure Turbine (HPT) of a Siemens gas turbine. The complete 3D data of the studied NGV can be found in the work of Yasa et al. [START_REF] Yasa | Performance of a nozzle guide vane in subsonic and transonic regimes tested in an annular sector[END_REF]. The simplification consists in adapting a linear cascade model as a representation of the 3D annular NGV cascade at mid-span. Afterwards, the mesh is built around the central unwrapped airfoil. It consists in a structured hexahedral-mesh with 18 blocks adapted manually to the high curvature stator blade. These hexa-cells are aligned along the local flow direction with the aim of capturing flow gradients in an accurate manner. The computational domain together with the mesh blocking strategy can be visualized in Fig. 1. For the LES calculations this 2D topology is extruded in the spanwise flow direction for capturing 3D effects. Periodic boundary conditions at the upper and lower boundaries are used in order to simulate an infinite linear cascade. Periodic boundary conditions are also applied at the extruded faces to model a channel of infinite width. Note that the height of the annulus enclosing the guide vanes is one order of magnitude lower than the inner radius of the cascade [START_REF] Atassi | Interaction of acoustic and vortical disturbances with an annular cascade in a swirling flow[END_REF], thus a linear cascade model can be adopted. The 2D mesh is hence extruded with 75 layers of chosen width z = 0.0001 [m]. The width of the cell in the spanwise direction is based on the transverse Taylor micro-scale ≈ 0.0001 [m]. The total mesh depth in the spanwise direction L z is equal to 0.1095C ax , where C ax represents the axial cord length, which is equal to 0.0685 [m]. The blade to blade distance H is, instead, equal to 0.08345[m], which is also the distance between the upper and lower periodic boundaries. Sponge zones with gradual grid stretching are used at the inlet and at the outlet boundaries in order to damp

Numerical Setup

The numerical simulations are performed using the fully coupled density based solver dbnsTurbFoam, developed within the foam-extend framework 1 . It is an open-source, cellcentred, finite volume solver for compressible flows. The compressible Navier-Stokes equations are solved in their conservative form, which is necessary for capturing phenomena characteristic to transonic flows [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF][START_REF] Lax | Systems of conservation laws[END_REF]. The solver adopts the so-called method of lines to decouple spatial and temporal discretizations, reducing the governing equations to a system of ordinary differential equations (ODEs) [9]. The spatial convective discretization is performed with the Rusanov flux scheme [START_REF] Rusanov | The calculation of the interaction of non-stationary shock waves and obstacles[END_REF], which is an approximation of the exact Rienman fluxes. This scheme is first order in space but the solver dbnsTurbFoam interpolates linearly the state from the cell center to the face center which makes the method of second order. Moreover, in order to avoid spurious oscillations near regions with sharp gradients, the Barth Jespersen limiter [START_REF] Barth | The design and application of upwind schemes on unstructured meshes[END_REF] is applied. The time integration is performed using an explicit 4-stage low-storage Runge-Kutta algorithm. The solver has been previously validated on two canonical compressible flow simulations by Chandramouli et al. [START_REF] Chandramouli | Numerical characterization of entropy noise with a density based solver[END_REF]. The LES simulations are initialized with the solution of the 2D Reynolds Averaged Navier-Stokes (RANS) equations using Favré averaging. The k -ωSST compressible turbulence model was adopted for the 2D RANS calculations. The LES turbulence model relies, instead, on the one-equation model for the subgrid-scale kinetic energy [START_REF] Yoshizawa | Bridging between eddy-viscosity-type and second-order turbulence models through a two-scale turbulence theory[END_REF].

Baseline Case

The baseline case of the present study refers to the numerical simulation matching the off-design transonic condition of the experiment [START_REF] Yasa | Performance of a nozzle guide vane in subsonic and transonic regimes tested in an annular sector[END_REF], without any inflow disturbances. The boundary conditions at the inlet are specified in terms of total pressure p 01 = 187100[Pa] and total temperature T 01 = 304.51[K]. The inflow velocity direction is perpendicular to the inlet boundary patch.

At the outlet, a face averaged static pressure of p 2 = 104370[Pa] is imposed. The airfoil walls are treated as no-slip adiabatic walls. Due to the fine refinement region near the wall (y + ≈ 5), the use of wall function has been avoided. A snapshot of the instantaneous Mach, pressure and temperature fields (LES data) can be visualized in Fig. 2. The Mach number field shows the essential flow features of the simulated case: the flow is subsonic at the inlet of the domain with the local Mach number equal to M in = 0.15, matching the experimental conditions.

Then, the flow is accelerated up to sonic condition at the throat of the nozzle guide vane passage. Passed the throat the flow interacts with compression and expansion waves emitted from the shed vortices at the trailing edge (seen from Fig. 2b). These waves are also impinging on the suction side and then are reflected back towards the normal direction of the blade surface. The flow accelerates further on the suction side just above sonic conditions according to a Prandtl-Mayer expansion. Towards the end of suction side, one can observe a system of weak shocks present, unsteady and moving with its own dynamics. The results from the 2D RANS (used for the initialization) and the 3D LES are then compared with the experimental data. In Fig. 3a andb, the NGV blade loading and the circumferential pressure distribution at x/C ax = 40% downstream from the blade are presented in terms of the isentropic Mach number M is , which is defined as

M is = 2 γ -1 ⎡ ⎣ p 01 p γ -1 γ -1 ⎤ ⎦ ; (1)
where p 01 is the inlet total pressure, p is the local static pressure and γ is the ratio of the specific heats. The blade loading is expressed as a function of the axial normalized coordinate x/C ax while the circumferential pressure distribution is expressed as a function of the vane phase, as presented in the work by Yasa et al. [START_REF] Yasa | Performance of a nozzle guide vane in subsonic and transonic regimes tested in an annular sector[END_REF]. For the pressure side, notably a very good agreement is found between the experiments and both the numerical simulations. The matching is also fairly good on the suction side, with a slight underprediction after the throat region. The pressure distribution along the blade involves different features due to the physics of the flow field. On the pressure side the flow is accelerated up to sonic conditions at almost the end of the blade; the entire pressure sides behaves like a converging portion of a nozzle. On the suction side, the flow accelerates up to chocking conditions until x/C ax = 0.75. Then, the Mach number is increased due to a supersonic expansion on the final part of the blade. Finally, the Mach number drops just below sonic conditions due to the presence of the weak normal shock. The location of the shock is just slightly overpredicted by the 2D RANS. Regarding the pressure distribution downstream of cascade, Fig. 3b, the topology of the domain allows, by construction, to reproduce only one passage, and the subsequent resulting trend is periodic. The experimental data are available for almost three complete passages. It can be noticed that the LES data are in perfect agreement in phase 2; 2D RANS also show an overall good agreement in that region. The mismatch for few other locations is due to the fact that experimental data are not perfectly periodic. Then, the total pressure ratio is defined as

= p 0 p 01 , ( 2
)
where p 0 is the local total pressure. This ratio is plotted in Fig. 4 for two locations downstream of the passage, at x = 7.1%C ax and at x = 48.4%C ax respectively. The same circumferential coordinate y/Pitch as Yasa et al. [START_REF] Yasa | Performance of a nozzle guide vane in subsonic and transonic regimes tested in an annular sector[END_REF] is adopted. The increasing wake region extension downstream the NGV passage can be clearly seen from the pictures. The total pressure losses are also reduced thanks to the enhanced mixing in the wake region. The asymmetric profile at the closest location to the trailing edge is very well captured by the 3D LES. The resulting profile from the 2D RANS appears slightly more symmetric. In the furthest location, the 3D LES exhibits a much closer profile to the experimental data than the 2D RANS, for which the wake width is underpredicted. Overall, a fair agreement is found for both 3D LES and 2D RANS results in terms of wake width at 7.1%C ax downstream the passage. Instead, At 48.4%C ax the 3D LES results agree better with experimental data than the RANS. Finally, regarding the maximum pressure losses, both simulation results overpredict them at 7.1%C ax while the 3D LES results are closer to the experimental ones at 48.4%C ax .

Forced Case

Indirect combustion noise results from non-isentropic temperature fluctuations that are accelerated at the combustor outlet [START_REF] Polifke | Constructive and destructive interference of acoustic and entropy waves in a premixed combustor with a choked exit[END_REF].

In the present LES, temperature disturbances are periodically forced at inlet as planar waves for simulating entropy waves, according to:

T1 = T 1 + T 1 A sin(2πf 1 t), ( 3
)
where T1 is the oscillating static temperature, T 1 , the inlet static temperature of the baseline case, A the amplitude of the disturbance and f 1 the temporal frequency. The amplitude A has been set equal to 4.8% and the the frequency f 1 to 1000 [Hz]. The choice of the frequency is dictated by the high computational cost of the 3D LES. However, this frequency is still in the limit range of real case scenarios; previous studies in literature indicated that combustion noise has its peak intensity typically in the low frequency range around 200 -1000[Hz] [START_REF] Dowling | Combustion noise[END_REF][START_REF] Hassan | Scaling of combustion-generated noise[END_REF][START_REF] Kotake | Combustion noise: effects of the shape and size of burner nozzle[END_REF][START_REF] Ihme | On the generation of direct combustion noise in turbulent non-premixed flames[END_REF]. The compact assumption states that the wavelength of the incoming entropy waves is much larger (by a factor of 10) than the blade axial chord. Still, similar simulations have been conducted by Papadogiannis et al. [START_REF] Papadogiannis | Assessment of the indirect combustion noise generated in a transonic high-pressure turbine stage[END_REF], which found that the limiting frequency for the validity of the compact assumption [START_REF] Cumpsty | The interaction of entropy fluctuations with turbine blade rows; a mechanism of turbojet engine noise[END_REF] is around 2200 [Hz]; which results in a less strict condition for a NGV blade in real case scenarios. Snapshots of the instantaneous entropy wave field in a spanwise section can be visualized in Fig. 5 for both the baseline and the forced case. The entropy wave, denoted by the superscript s , is defined as w s = s /c p = T /T -(γ -1)p /γ p, where the primed quantities are respectively the entropy (s ), temperature (T ) and pressure (p ) fluctuation fields; the overlined ones are the time averaged temperature (T ) and pressure (p) field, and c p and γ are respectively the specific heat at constant pressure and the specific heats ratio. For the baseline case in Fig. 5a, the flow accelerates up to sonic conditions at the throat. Then the "cross"-like pattern of waves, shed from the trailing edge vortices, is present just after the throat. Finally a system of weak shocks is visible at the end of the blade. Figure 5b shows, instead, the features of the flow-field with planar entropy waves injected at the inlet. The planar forcing is revealed at the inlet by the presence of vertical stripes in the entropy fluctuation field. Once the waves approach the blades they get distorted by the large gradients of mean flow quantities across the vane, becoming aligned with the flow direction at the end of the passage. Then, one axial chord downstream from the blade, it is possible to see almostplanar entropy waves propagating downstream. A temporal fast Fourier transform (FFT) has been conducted for a set of two-dimensional snapshots of velocity, temperature and pressure fields, sampled at a frequency equal to four times the vortex shedding frequency, in order to avoid aliasing. This frequency is widely sufficient to capture possible features of the flow-field appearing at f 1 = 1000 [Hz]. For the baseline and the forced case, the amplitude of the acoustic wave field w ± = p /γ p is pictured in Fig. 6 from the FFT at the inlet forcing frequency f 1 . The forced case, Fig. 6b, shows several regions of interest with noticeable amplitudes highlighted by the respective markers. No relevant amplitudes were instead detected for the baseline case of Fig. 6a at f 1 = 1000[Hz], with the exception of a small region on the suction side of the blade, characterized by the presence of compression and expansion waves at the particular frequency f 1 . For the forced case, a glowing yellowwhite area can be seen on the suction side, towards the end of the blade (white rectangle in Fig. 6b). This area is overlapped with the system of weak shocks, visible in the instantaneous Mach, pressure and temperature instantaneous flow-field of Fig. 2. This area is also visible in the baseline case of Fig. 6a, but with a less extended area of high pressure amplitudes. The convergent section of the vane passage also exhibits new notable amplitudes up to the throat section with respect to the baseline case (white circle of Fig. 6b). At the outlet (white ellipse of Fig. 6b), the acceleration of the entropy spots through the NGV passage and the acoustic waves generated on the suction side and propagated downstream, are strong enough to produce a significant pressure fluctuation trace with non-zero modulus [START_REF] Papadogiannis | Assessment of the indirect combustion noise generated in a transonic high-pressure turbine stage[END_REF]. A possible explanation can be that entropy waves, interacting with strong mean pressure and mean velociy gradients, give rise to an entropy-acoustic modal interaction. Considering for instance a weakly non linear Kovasznay decomposition [START_REF] Chu | Non-linear interactions in a viscous heat-conducting compressible gas[END_REF][START_REF] Garnier | Large Eddy Simulation for Compressible Flows[END_REF] of all compressible flow variables, acoustic production due to acoustic-entropy coupling can arise due to scattering: i.e. the acoustic waves are scattered by the convected entropy spots passing through the nozzle. According to Kovasznay decomposition [START_REF] Chu | Non-linear interactions in a viscous heat-conducting compressible gas[END_REF][START_REF] Garnier | Large Eddy Simulation for Compressible Flows[END_REF] the acoustic production scattering term is ∂ t (∇ • s e u p ), where s e is entropy perturbation of the entropy mode and u p is the velocity perturbation of the acoustic mode. Although no acoustic forcing is imposed by the entropy waves injected at the inlet, distinctive higher amplitudes of the acoustic waves in the whole domain reveal that a pressure mode with a distinctive peak at the frequency f 1 appears, confirming the indirect noise generation mechanism [START_REF] Papadogiannis | Assessment of the indirect combustion noise generated in a transonic high-pressure turbine stage[END_REF]. In the forced case, new regions of high amplitudes pressure fluctuations appear, which are not present in the baseline case. A dark line (white arrow in Fig. 6b) the neighbourhood of sonic throat identifies also a region of lower amplitude. In the convergent region, only the reflected waves are present because the waves created downstream of the sonic line cannot travel upstream. The attention will move now to the glowing area with the highest amplitude, in Fig. 6b. Entropy waves seem to strongly interact with the system of weak normal shocks, therefore a sampling line has been traced in this particular zone as sketched in Fig. 7. For this line, the pressure fluctuation field has been sampled in time for a number of 100 spatial points with the resolved frequency range of 50 [KHz]. Then, the sampled field is projected on a space-time grid where the spatial coordinate is represented by the non-dimensional distance x/C ax and the temporal coordinate by a non dimensional time τ * = tU/D, U being the velocity calculated from the Mach number M is = 0.95 and D the trailing edge diameter, equal to 0.0022 [m]. With this method the characteristic dynamics of the system of normal shocks can be identified and the present analysis will be conducted for both the baseline and the forced case. The pressure fluctuation field is expressed in terms of its non dimensional counterpart defined as the acoustic wave w ± . Figure 8 shows the non dimensional acoustic wave projected on the space-time diagram for both the baseline and the forced case. Three regions of interest can be identified in the graph. Region 1 corresponds to the zone where the position of the shock changes slowly during time; this region has fairly similar characteristics for both the baseline and the forced case. The frequency content in Region 1 is in the range of 200-900[Hz]. Region 2, instead, corresponds to the glowing zone for which the amplitude of the pressure fluctuations is the highest in Fig. 6. For the forced case in Fig. 8b, yellow-white stripes at regular intervals of time are visible in Region 2. The frequency of these intervals is evaluated by the means of a FFT at a fixed position in space: it corresponds exactly to the inlet forcing frequency f 1 = 1000[Hz]. This feature does not appear for the baseline case in Fig. 8a and the shock movement is not organized at any particular frequency of interest; its frequency content is instead in the range of 500 -3000[Hz]. Finally, Region 3 corresponds to the initial part of the wake, where compression and expansion waves are emitted by vortical structures and appear with higher frequencies. The structure of Region 3 is similar for 1 , both at plane 1. On the other hand, the acoustic transmission coefficient T r + , denoted by the superscript + , is defined as the ratio between the transmitted acoustic wave w + 2 at plane 2 and the incoming entropy wave w s 1 at plane 1. Finally, the entropy transmission coefficient T r s , denoted by the superscript s , is defined as the ratio between the transmitted entropy wave w s 2 at plane 2 and the incoming entropy wave w s 1 at plane 1. The definition of the entropy and acoustic waves follows the one showed in Section 5. Table 1 summarizes the results for these coefficients.

Conclusions

The indirect combustion noise generation mechanism has been addressed by the means of 3D LES for a simplified nozzle guide vane passage of a high pressure turbine. The solver adopted for the simulations was first validated against experimental data available for offdesign transonic conditions without temperature forcing. Overall the numerical results are in very good agreement with experimental data. The larger discrepancies were only found in terms of maximum total pressure losses which might be due to the simplified 2D-extruded geometry adopted, for which the wake mixing might not be enhanced as in the real case scenario. In the second part of the work, a pulsating entropy wave-train was imposed at the inlet of the domain to simulate temperature non-uniformities coming from the combustion chamber. The resulting flowfield was analysed by the FFT of two-dimensional sections sampled during the simulation. At the forcing frequency imposed at the inlet, certain zones in the fluid respond with high amplitudes of the pressure fluctuations, which are not present in the baseline case. An interesting feature was observed regarding the system of shocks emitted from the trailing edge: their dynamics in a particular region (Region 2) was completely characterized by a frequency content at 1000 [Hz], which is the inlet forcing frequency of the entropy waves. Finally, reflection and transmission coefficients were compared to analytical models available in literature [START_REF] Cumpsty | The interaction of entropy fluctuations with turbine blade rows; a mechanism of turbojet engine noise[END_REF]. While the theoretical approach overpredicts the noise generated, the 3D simulations show that the entropy waves are highly distorted by the passage and lose their strength. This can be addressed from the entropy wave transmission coefficient which is reduced by 55%. The reflected acoustic wave, instead, is far weaker than the one predicted by the compact theory. Finally, the reflected acoustic wave is slightly lower than the analytical predicted one, but its amplitude is still strong. It is therefore the main contributor to the generated indirect combustion noise.

Fig. 1

 1 Fig. 1 Mesh topology a and blocking strategy b for the NGV linear cascade as a model of the 3D annular NGV cascade at mid-span

Fig. 2

 2 Fig. 2 Instantaneous Mach number a, pressure b and temperature c fields at midspan, LES data

Fig. 3

 3 Fig. 3 NGV blade loading a and downstream pressure distribution at x/C ax = 40%, line ( ) in Fig. 1b b; (•) experimental data; ( ) 2D RANS; ( ) 3D LES

Fig. 4

 4 Fig. 4 Total pressure ratio: at x/C ax = 0.071, line ( ) in Fig. 1b a and at x/C ax = 0.484, line ( ) in Fig. 1b b; (•) experimental data; ( ) 2D RANS; ( ) 3D LES

Fig. 5

 5 Fig. 5 Instantaneous entropy wave field for the baseline a and for the forced b case

Fig. 6

 6 Fig. 6 Amplitude of the acoustic at f 1 = 1000 [Hz] for the unforced baseline case (a) and for the forced case (b)

Fig. 7 Fig. 8

 78 Fig. 7 Line for the analysis of shock dynamics

Fig. 9

 9 Fig.9 Azimuthal planes adopted for the evaluation of reflection and transmission coefficients; ( ) upstream plane 1, ( ) downstream plane 2

Table 1

 1 Reflection and transmission coefficients.

			analytical	CFD
	Comparison between numerical simulations analytical model	|Ref -|	0.073	0.011
		|T r + |	0.088	0.068
		|T r s |	1	0 . 4 5

https://sourceforge.net/p/openfoam-extend/wiki/Home/

Acknowledgements

The simulations were performed at the PDC Center for High Performance Computing (PDC-HPC) and at the High Performance Computing Center North (HPC2N), on resources provided by the Swedish National Infrastructure of Computing (SNIC).

Compliance with Ethical Standards

Conflict of interests

The authors declare that they have no conflict of interest.

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


N Rougemaille

F El Gabaly

R Stumpf

A K Schmid

K Thürmer

N C Bartelt

J De La Figuera

Labyrinthine Island Growth during Pd/Ru(0001) Heteroepitaxy

Keywords: numbers: 68.55.Ac, 61.72.Bb, 61.72.Ff, 68.37.Nq

Using low energy electron microscopy we observe that Pd deposited on Ru only attaches to small sections of the atomic step edges surrounding Pd islands. This causes a novel epitaxial growth mode in which islands advance in a snakelike motion, giving rise to labyrinthine patterns.

Based on density functional theory together with scanning tunneling microscopy and low energy electron microscopy we propose that this growth mode is caused by a surface alloy forming around growing islands. This alloy gradually reduces step attachment rates, resulting in an instability that favors adatom attachment at fast advancing step sections.

At sufficiently high temperature, epitaxial growth is expected to occur by the flow of atomic steps: atoms deposited on flat terraces on the surface of a crystal diffuse across the surface until they encounter steps and are incorporated at low energy kink sites, causing the steps to advance. Since the utility of thin films grown in this way often depends on the surface step structure, much work has been devoted to characterizing the morphology of steps. For example, if diffusion along the step edges is fast step edges become smooth, whereas slow step diffusion produces dendritic step edges [1]. In single component systems the basic atomic processes governing the morphology of moving steps are increasingly well understood [START_REF] Pimpinelli | Physics of Crystal Growth[END_REF]. However, in heteroepitaxy, new types of cooperative behavior can arise [START_REF] Hwang | [END_REF][4][5]. In this work we report the discovery of a novel way in which surface alloying and step motion couple, leading to a distinctive labyrinthine island morphology during growth in which only fast moving step edges are good sinks for deposited atoms. By unraveling the cooperative atomic scale processes leading to these striking 100-nm-scale patterns we make an inroad into understanding the intricate connections between surface morphology and alloying.

We use the complementary capabilities of real-time low energy electron microscopy (LEEM) and atomic resolution scanning tunneling microscopy (STM) to study the evolution and microscopic structure of Pd deposited on the Ru(0001) surface. Pd and Ru are somewhat soluble in each other (the maximum solubilities are a few percent at 900 K [START_REF]Binary Alloy Phase Diagrams[END_REF]) and thus at high enough temperature alloying during growth might be expected. We grow submonolayers of Pd by physical vapor deposition onto well cleaned Ru(0001) substrates at 840 K at base pressures below 1×10 -10 torr. Previous diffraction studies have shown that monolayer (ML) films of Pd on Ru are pseudomorphic [START_REF] Kolaczkiewicz | [END_REF]. A sequence of LEEM images [8] of the growing surface during Pd deposition at a substrate temperature of 840 K is shown in Fig. 1. Initially Pd starts to decorate the Ru substrate steps, forming a uniform narrow band at step edges. However this uniform step flow quickly stops and instead growth continues only from a few widely separated locations leading to a snakelike motion of islands, as shown in Fig. 1 The observed growth mode is puzzling at first. All step edges should be good sinks for deposited Pd atoms, yet Pd is incorporated overwhelmingly at a few special regions of the step edge, causing them to advance. Also, growing islands should be easily able to coalesce. What atomic properties cause this growth behavior? One hypothesis is that the narrow islands are an example of a stress domain pattern; i.e., a narrow island width is chosen as a compromise between the energetic cost of creating step edges and the relaxation of surface stress at the island boundaries. For example, Pb overlayer islands grown on Cu(111) can exhibit a similar morphology [9], which has been characterized quantitatively as a stress domain pattern in thermodynamic equilibrium. We will show, however, that an entirely different, kinetic, mechanism is at work here.

To address the reasons for the snakelike growth we imaged the atomic structure of the surface with STM. We grew the Pd islands under the same conditions as in the LEEM experiments and then rapidly quenched the surface to room temperature before imaging it with STM. As seen in Fig. 2(a), the same meandering island shape was observed. Atomic resolution images with chemical contrast [10] show that there is a distinct difference between the atomic structure of the Ru terraces adjacent to steps that have apparently moved fast and those that have not: the terraces adjacent to the stationary step edges are heavily alloyed with Pd, with densities up to about 10%. On the other hand, very little alloying is observed adjacent to active, moving Pd step edges. These data immediately suggest that alloyed Ru regions have slowed the attachment of Pd adatoms on the Ru terraces to the edges of the Pd islands. This suggestion raises two questions: By what mechanism does the alloying impede Pd incorporation at steps, and why is the alloying restricted to the region adjacent to the stationary island edge? To answer these questions we performed density functional theory (DFT) calculations [START_REF]Computational unit cells contain 1×1 and 3×3 surface periodicities on six-and seven-layer Ru slabs[END_REF] within the Perdew-Burke-Ernzerhof [START_REF] Perdew | [END_REF] generalized gradient approximation of the energetics of Pd on Ru(0001) using the VASP code [13].

First, we explain why alloying is rare on terraces far away from Pd step edges. Consider the energetics of the process labeled "A" in Fig. 3(a). The calculated energy difference between dispersed Pd substitutional atoms on the Ru terraces and Pd atoms in the overlayer is E f = 0.24 eV. At the experimental temperature of 840 K this formation energy would give an equilibrium alloy density of approximately exp(-E f / kT) = 0.04. However, the lowest energy barrier we could find for a Pd adatom to exchange with a Ru atom in the top substrate layer [process "C" in Fig. 3(a)] is 2.48 eV. The resulting substitutional Pd and Ru adatom is 1.01 eV higher in energy than the Pd adatom initial state. Consequently, it is likely that the Ru adatom reexchanges with a substitutional Pd, because of the 1.01 eV energy gain and a barrier that is only 2.48 eV -1.01 eV = 1.47 eV. The reexchange is avoided only if the Ru adatom is trapped at a step instead, where it ultimately gains the Ru adatom formation energy of 1.38 eV [process "B" in Fig. 3 We identify two mechanisms by which the existence of a step edge alloy can hinder the attachment of Pd adatoms at Pd step edges. First, the alloy on the terraces might impede the diffusion of Pd toward steps. Supporting this scenario, our calculations show a repulsion between Pd adatoms and Pd atoms substituted into the Ru substrate: only three of the six threefold hollow sites adjacent to an embedded Pd are local minima with a binding energy that is reduced by 0.14 eV over clean Ru(0001). This repulsion impedes diffusion of Pd near Pd islands. The alloy concentration ρ near the step edges observed in STM is, at most, 10%. Our calculations [START_REF]The diffusion coefficient D of Pd adatoms was computed by creating a random arrangement of embedded Pd in a 500×500 array of Ru substrate atoms. A Pd atom placed on this surface hopped to adjacent threefold sites according to the Metropolis Monte Carlo algorithm with the energy landscape from our DFT calculations[END_REF] for this density give a decrease in the diffusion coefficient D by a factor of 2. If the alloying region is only a few nm wide as observed in STM, this decrease in D is by itself not sufficient to effectively block the step edges. A second alloying effect that slows the Pd step advancement is that covering alloyed Pd atoms by the advancing Pd layer costs energy compared to advancing over the pure Ru substrate.

The energy penalty is 0.5 eV per embedded Pd, almost equal to the typical binding energy of Pd to the Pd step of 0.61 eV. Given the low Pd adatom diffusion barrier of 0.19 eV it is reasonable to expect that local equilibrium along steps can be established: Pd adatoms that attach to a Pd island step located over a heavily alloyed substrate region are more likely to detach because of their higher energy; i.e., the Pd adatom density in equilibrium with a step edge section advancing over an alloyed terrace would be higher than over a clean Ru terrace. This would cause concentration gradients on the terraces leading to a Pd flux away from alloyed step edges and toward the edges that are in clean regions. Thus growth of Pd over alloyed regions will be avoided. Notice that the higher Pd adatom concentration produced by this mechanism would reinforce the mechanism discussed above by promoting the alloying of Pd atoms into the Ru substrate.

In any case, the degree of alloying near a stationary step edge, and hence the barrier to Pd attachment, depends on time. We checked experimentally the time dependence of the Pd attachment barrier by starting the growth of narrow Pd islands and then interrupting the Pd dosing for different periods of time. Stopping the deposition causes the growth to cease immediately. If we restart deposition three or more minutes later, none of the existing active regions resume growing. Instead, new islands are nucleated and new active step edges emerge from them, much as during the early state of the experiments. With shorter interruptions the original active regions of the step edges resume growth. This suggests that the alloy concentration necessary to block step edges takes some time to build up.

Can this time dependent alloy concentration profile near step edges account for the observed growth mode? When step edges are moving slowly, the alloy concentration should behave similarly to stationary steps: the alloy concentration next to the step edges increases until it reaches a relatively large (close-to-equilibrium) value of about 5% and then gradually extends onto the adjoining terraces. As shown schematically in Fig. 3(b), this situation leads to an instability in the growth front: the flux of Pd attaching at the step edge will be largest where the alloy concentration is (by random fluctuations, initially) the smallest. This increased attachment and higher step velocity causes the alloy concentration in front of the step to be even lower (because the time spent by the step in any one position is less) and Pd attachment to be even faster. If the speed of the steps becomes sufficiently large, one can envision a steady-state situation where the alloy concentration near the step edge remains at a constant low value below the equilibrium density. We propose that these slow and fast regimes are separated by a distinct critical speed. Below the critical speed, attachment at the Pd step edges becomes increasingly unfavorable with time as the alloy concentration builds up. Above the critical speed, the step edge moves faster than the speed at which the increased alloy concentration spreads onto the terrace, so that concentration remains low and Pd attachment stays easy.

The existence of such a critical velocity has direct experimental support. In Fig. 1 for example, no step propagates at a velocity less than 45 nm/s. Remarkably, the propagation speed is always very close to this minimum velocity, despite strong variations in the environment and size of the moving front. To rationalize this result we note that moving more quickly than the critical speed (given that the local flux remains constant) would result in a narrowing of the growth front (to conserve mass). This shrinking would be resisted by various effects. For example, a narrow growth front would require adatom concentration gradients perpendicular to the growth direction. These gradients would be unstable and would broaden (and thus slow) the growth front. The increased energy associated with large step curvatures would also oppose narrowing. Thus the growth front appears to propagate at the slowest possible sustainable speed [16].

In this scenario, the critical velocity should only depend on the rate of alloying but not on the overall deposition flux. This prediction is tested in LEEM experiments that use different Pd deposition rates. We have measured the growth velocity of deposition rates of 1/650, 1/1400, and 1/6300 ML/s. The measured active step velocities are 50, 44, and 30 nm/s, respectively. So, an order of magnitude change in flux leads only to a relatively small change in velocity, roughly consistent with our proposed model. Given approximately the same number of active regions per unit area, mass conservation then requires that the higher the flux, the wider the active regions. This effect is clearly observed in Fig. 4. Also the active regions tend to get narrower, but move at the same speed, at the late stages of growth when there are more active regions and the flux per active region is less. Notice that the flux dependence of the width rules out that the labyrinthine patterns can be explained in terms of equilibrium surface stress domains. Corroboration that surface alloying is responsible for the novel growth process comes by codepositing Ru with the Pd. According to the DFT calculations described above, one would expect the blocking Pd atoms embedded in the Ru substrate to be displaced by deposited Ru adatoms [START_REF]The process for a Ru adatom exchanging with a substitutional Pd atom has a barrier of only 1.48 eV[END_REF]. Consistent with this we find that even small simultaneous codeposition of Ru [with rates of 1/650 ML/s (Pd) and 1/1500 ML/s (Ru)] causes Pd to be uniformly incorporated into all step edges.

In summary, we have found strong evidence that surface alloying can impede step-flow growth during Pd/Ru epitaxy at around 840 K, leading to distinctive snakelike motion that gives rise to labyrinth patterns. The conclusion that an instability causes the observed growth mode might be applicable to other heteroepitaxial systems with surface alloying, in particular if the deposited species has weaker bonds and lower surface energy than the substrate species. In this general case the higher adatom formation energy of the substrate species should kinetically limit alloying on substrate terraces to regions close to step edges with all its consequences. The kinetic mechanism for this pattern formation is very different from the explanation of similar patterns in terms of surface stress domains, as in, e.g., Pb/Cu(111). It is plausible that labyrinthine growth also could occur in systems where impurity gas adsorption [START_REF] Engelhardt | [END_REF]

3 FIG. 1 .

 31 FIG. 1. Sequence of LEEM images of the Pd growth process on Ru(0001) at 840 K. Field of view is 6 µm. The monoatomic Ru steps (the four parallel lines that cross the images) are separated by approximately 1.5 µm. The deposition rate is 1 ML in 421 s. After first uniformly decorating the step edges in (a), the Pd grows only from several distinct locations. In (b) and (c) the self-avoiding nature of different growth fronts is evident: Colliding growth fronts turn aside or split rather than coalesce.

FIG. 2 .

 2 FIG. 2. STM images after deposition of 0.2 ML of Pd. (a) View of an entire Pd island. (b),(c) Atomic resolution images of particular regions of the islands shown in (a). Close to the step section (b) the Pd island is alloyed with Ru atoms (circled) and the Ru terrace is alloyed with Pd atoms (one is marked by a square). The degree of alloying is much lower at the active step section (c).

  (a)]. Thus the alloying of Pd into the Ru surface by the adatom process is much more likely near steps than far away [14]. Consistent with this picture is the observation by STM of Ru atoms embedded in the Pd overlayer [Fig. 2(b)] with a density that parallels the density of alloyed Pd on the nearby Ru terrace. Furthermore, the density of Pd in the Ru terraces adjacent to the step edges is about 5%, on the order of what one would expect on the basis of the DFT formation energy.

FIG. 3 .

 3 FIG. 3. (a) Various processes involved in the PdRu alloy formation. (b) Schematic of the instability of Pd step attachment caused by variations in alloy concentration in front of the step edges. The dotted lines represent lines of equal concentration of Pd alloyed into the Ru terraces.Regions with locally less alloy (left) become better sinks of Pd, leading to larger growth rates (right) and a further decrease in alloy concentration.

FIG. 4 .

 4 FIG. 4. LEEM images of Pd deposited on Ru(0001) at deposition rates and times of (a) 1/650 ML/s and 42 s and (b) 1/6300 ML/s and 190 s. The deposition temperature is 840 Κ; the field of view is 6 µm. The averaged speed of the advancing active regions is 50 nm/s and 30 nm/s, respectively.

  creates step attachment barriers. This research was partly supported by the Office of Basic Energy Sciences, Division of Materials Sciences, U.S. Department of Energy under Contracts No. DE-AC04-94AL85000 and No. DE-AC02-05CH11231, and by the Spanish Ministry of Science and Technology through Project No. MAT2006-13149-C02-02.


Antoine Waret
email: antoine.waret@grenoble-inp.org

Megumi Kaneko
email: megkaneko@nii.ac.jp

Alexandre Guitton
email: alexandre.guitton@uca.fr

Nancy El

LoRa Throughput Analysis with Imperfect Spreading Factor Orthogonality

Keywords: LoRa, Spreading Factor, Uplink Throughput, Imperfect Orthogonality

LoRa is one of the promising techniques for enabling Low Power Wide Area Networks (LPWANs) for future Internet-of-Things (IoT) devices. Although LoRa allows flexible adaptations of coverage and data rates, it is subject to intrinsic types of interferences: co-SF interferences where end-devices with the same Spreading Factors (SFs) are subject to collisions, and inter-SF interferences where end-devices with different SFs experience collisions. Most current works have considered perfect orthogonality among different SFs. In this work, we provide a theoretical analysis of the achievable LoRa throughput in uplink, where the capture conditions specific to LoRa are included. Results show the accuracy of our analysis despite approximations, and the throughput losses from imperfect SF orthogonality, under different SF allocations. Our analysis will enable the design of specific SF allocation mechanisms, in view of further throughput enhancements.

I. INTRODUCTION

As the amount of mobile data traffic will rapidly increase during the upcoming years (studies forecast 50 billion Internet of Things (IoT) devices by 2020), new spectrum access strategies adapted to high device densities are ever more crucial. LoRa [START_REF]LoRa Modulation Basics -AN1200.22[END_REF] is one of the prominent candidates for Low Power Wide Area Networks (LPWANs), providing wide communication coverage with low power consumption, at the expense of data rate. Operating in license-free ISM bands (i.e., 868MHz in Europe), the LoRa PHY layer uses a chirp spread-spectrum modulation where different Spreading Factors (SFs) tune the chirp modulation rates. Lower SFs such as SF7 allow for higher data rates but reduced transmission range, whereas higher SFs such as SF12 provide longer range at lower data rates. On top of the LoRa PHY layer, the higher layers were defined by the LoRa Alliance and referred as LoRaWAN [START_REF]LoRaWAN -What is it?, A technical overview of LoRa and LoRaWAN[END_REF]. In particular, the MAC protocol is based on a pure ALOHA access with duty cycle limitations. The LoRaWAN network architecture is a star-like topology where end-devices communicate with gateways over several channels.

Most studies on LoRa scalability so far assumed a perfect orthogonality among SFs, thereby creating virtual channels where multiple users with different SFs could simultaneously operate in the same channel and hence boost the achievable system throughput. Thus, a number of works have considered the effect of co-SF interference only, where end-devices Fig. 1. LoRa system setup -Case of SF-distance allocation using the same SF on the same channel are subject to collisions [3][4]. In particular, the outage probability of a LoRa system under co-SF interference was analyzed in [START_REF] Georgiou | Low Power Wide Area Network Analysis: Can LoRa Scale?[END_REF], where a signal could be captured if its Signal-to-Interference-plus-Noise Ratio (SINR) was higher than 6 dB. As the number of devices increased, it was shown that those co-SF interferences were causing a scalability limit. However, recent studies have pointed out the fact that SFs were not perfectly orthogonal among themselves [START_REF] Croce | Impact of Spreading Factor Imperfect Orthogonality in LoRa Communications[END_REF]. Thus, the effect of inter-SF collisions was investigated through computer simulations and/or experiments. Namely, [START_REF] Croce | Impact of Spreading Factor Imperfect Orthogonality in LoRa Communications[END_REF] [START_REF] Zhu | Evaluation of LoRa Receiver Performance under Co-technology Interference[END_REF] showed that inter-SF interferences could considerably decrease LoRa performance, especially for high SFs where frames have a greater time on air.

In this work, we propose a theoretical analysis of the achievable throughput on the uplink of a LoRa network, encompassing the effects of co-and inter-SF interferences. To ensure a successful transmission, a packet must thus satisfy three conditions: 1) its SNR is above the reception threshold, 2) its SINR under co-SF interference is above the co-SF capture threshold, and its SINR under inter-SF interference is above the inter-SF capture threshold. Considering two different types of SF allocations, we theoretically derive the achievable throughput expressions for both perfect and imperfect SF orthogonality. Simulation results show the accuracy of our analytical expressions despite the necessary approximations, as well as the impact of the various types of interferences and SF allocations on the overall system performance.

II. SYSTEM MODEL

We consider one cell of radius R with one gateway located at its center, as depicted in Fig. 1. There are N end-devices uniformly distributed within the cell. We denote by d i the distance from end-device i to the gateway. Since the goal of our analysis is to derive the achievable rate by LoRa, we assume that all end-devices transmit in a single channel of bandwidth BW = 125kHz and that they all have packets to transmit. This corresponds to the pure ALOHA access as in LoRaWAN with saturated traffic 2 . We consider M = 6 SFs, for m = {m min , ..., m max }, with m min = 7 and m max = 12, with symbol times

T m = 2 m BW . The bit-rate R m of SF m is [1] R m = m × CR 2 m BW ,
where CR is the coding rate defined as 4/(4 + n) with n ∈ {1, 2, 3, 4}. Lower SFs allow higher data rate but lower communication range whereas higher SFs provide longer range at the expense of data rate (see Table I). Two types of SF allocation will be investigated. In the first one, the spreading factors are uniformly distributed, i.e., every end-device has a probability p m = 1 M of selecting SF m . We refer to this allocation as SF-random. In the second type of allocation referred to as SF-distance, spreading factors are assigned according to the distance d i . A device located inside the annulus defined by the smaller and larger circle radii l m-1 and l m , respectively, has SF m . The distance threshold

l m for SF m is given by l m = ( P0A(fc) θrx m ) 1 α , where A(f c ) = (f 2 c ×10 -2.8
) -1 is the deterministic loss in the path loss model

L i = A(fc) d α i
as in [START_REF]Propagation data and prediction methods for the planning of indoor radiocommunication systems and radio local area networks in the frequency range 300 MHz to 100 GHz[END_REF], f c the carrier frequency and α the path loss exponent. θ rxm is the receiver sensitivity of SF m (see Table I). All nodes transmit at equal power P 0 . We assigned to l 6 and l 12 the origin of the cell and its radius respectively, i.e., l 6 = 0 and l 12 = R. The ranges for each SF are given in Table I. The probability of selecting SF m for the SF-distance allocation is then given by p m = ∫ lm lm-1 h(r)dr, where h(r) is the pdf of the position of an end-device in the cell at distance r from the gateway. For uniform distribution of devices within a cell of radius R, we get h(r) = 2r R 2 . The instantaneous SNR γ i of end-device i is defined as

γ i = P 0 |h i | 2 L i /σ 2
n , where |h i | 2 is the channel gain between end-device i and the gateway (for Rayleigh fading, h i ∼ CN (0, 1)). σ 2 n = -174 + NF + 10 log(BW ) [dBm] is the AWGN power and NF, the receiver noise figure.

Based on [START_REF] Croce | Impact of Spreading Factor Imperfect Orthogonality in LoRa Communications[END_REF], it is assumed that in the event of a collision between frames of different SFs, one signal is received successfully if its SINR is higher than its "InterSF capture threshold" in Table I. Moreover, if there are several signals with equal SFs transmitting on the same frequency simultaneously, the gateway is able to successfully receive one of them if its SINR is higher than 6 dB, for any SF m [3][8]. Therefore, both types of interferences will be considered.

III. PROPOSED THROUGHPUT ANALYSIS

An end-device's packet transmission in uplink is successfully received at the gateway if the three following conditions are fulfilled:

1) Reception condition: Signal power must be above the SF-specific threshold q SFm ,

P i caprx m = P (γ i ≥ q SFm ), (1)
which is the probability that a received signal from end-device i at a distance d i from the gateway has a SNR γ i above the threshold q SFm (Table I).

2) Co-SF capture condition:

The co-SF SINR of end-device i is defined as

SINR i coSF = γ i ∑ K k=1 γ k + 1 , ( 2
)
where K is the number of end-devices with the same SF as end-device i. A signal is successfully received if its SINR i coSF is above q coSF . Therefore, the second condition is given by

P i capcoSF = P ( γ i ∑ K k=1 γ k + 1 ≥ q coSF ) . ( 3
)
3) Inter-SF capture condition: As shown in [START_REF] Croce | Impact of Spreading Factor Imperfect Orthogonality in LoRa Communications[END_REF][6], SFs are not perfectly orthogonal: a signal at SF m faces interferences from all N -K signals on other SFs. The inter-SF SINR of end-device i is defined as

SINR i intSF = γ i ∑ N p=K+1 γ p + 1 . ( 4
)
A transmission is successful if SINR i intSF is higher than the threshold q iSFm , thus the third condition is given as

P i capintSF m = P ( γ i ∑ N p=K+1 γ p + 1 ≥ q iSFm ) . ( 5
)
Therefore, the uplink throughput τ can be expressed as

τ = mmax ∑ m=mmin R m × P success (SF m ), (6)
where R m is the bit-rate of SF m and P success (SF m ), the probability of a successful transmission. Next, we analyze the throughput under perfect and imperfect SF-orthogonality, for the SF-distance case.

A. Perfect Orthogonality

We assume first that SFs are perfectly orthogonal, i.e., no end-device suffers inter-SF interferences. Hence, the probability of a successful transmission P success (SF m ) is given by

P success (SF m ) = N ∑ j=1 ( N j ) p j m (1-p m ) N -j P (cap rx , cap coSF ), (7)
where j denotes the total number of end-devices at SF m , ( N j ) p j m (1 -p m ) N -j is the probability of having j enddevices among N at SF m and P (cap rx , cap coSF ) is the joint probability for reception condition and co-SF capture.

1) For j = 1: the end-device is not subject to co-SF interferences, thus only the reception condition needs to be satisfied, P (cap rx , cap coSF ) = P i caprx m .

First, we determine P i caprx m for the SF-distance case. Given our assumptions, the SNR γ i is modeled as an exponential random variable with mean γ i . Therefore,

P i caprx m = P (γ i ≥ q SFm |γ i ) × P (γ i )
, where q SFm is the specific threshold of SF m . Defining

γ i = c r α i , where c = P0×A(fc) σ 2 n
is the path-loss constant, we can now rewrite,

P i caprx m = ∫ lm lm-1 exp ( - q SFm r α i c ) × 2r i R 2 dr i . ( 8
)
Although this integral cannot be expressed in closed form, it can be efficiently determined by numerical methods.

2) For j ≥ 2: both the reception and co-SF conditions must be fulfilled. As q SFm ≤1 in linear for all SFs whereas q coSFm = 4 (6 dB) for all SFs as explained in Section II, if co-SF capture is satisfied, so is the reception condition, hence

P (cap rx , cap coSF ) = P i capcoSF .
In case of co-SF interferences, there are j-1 interferers,

P i capcoSF = P ( γi ∑ j-1 k=1 γ k +1 ≥ q coSF )
, which is developed using random instantaneous SNR variables γ k and random average SNR (position) variables γ k as P ( γi ∑ j-1 k=1 γ k +1 ≥ q coSF |γ 1 , ..., γ i , ..., γ j-1 , γ 1 , ..., γ j- 1) ×P (γ 1 , ..., γ i , ..., γ j-1 , γ 1 , ..., γ i-1 , γ i+1 , ..., γ j-1 ).

Marginalizing over γ 1 , ..., γ j-1 and making the change of variable γ i = c r α i , we get by independency of user channels,

P i capcoSF = ∫ lm lm-1 exp(- q coSF r α i c ) × j-1 ∏ p=1 p̸ =i ∫ lm lm-1 h(r p )dr p 1 + q coSF ( ri rp ) α h(r i )dr i . We define I(r i ) = ∫ lm lm-1 h(r)dr 1+qcoSF ( r i rp
) α . Nodes are uniformly distributed within the cell, therefore ∏ j-1 p=1 I(r i ) = [I(r i )] j-1 . As a result, the expression becomes,

P i capcoSF = ∫ lm lm-1 exp ( - q coSF r α i c ) [I(r i )] j-1 h(r i )dr i . (9)
In particular, for α = 4, the primitive function of I(r i ) is

J(r i , r) = ( r R ) 2 - ( r i R ) 2 √ q coSF arctan ( 1 ( ri r ) 2√ q coSF ) . ( 10
)
In case of co-SF interferences, the interferers are the enddevices with the same SF as end-device i, thus with the same distance boundaries. The expression of I(r i ) becomes

I(r i ) = J(r i , l m ) -J(r i , l m-1 ). (11)
Therefore, ( 7) can be written for SF-distance allocation with perfect SF orthogonality as,

P success (SF m ) = ( N 1 ) (1 -p m ) N -1 × P i caprx m + N ∑ j=2 ( N j ) (1 -p m ) N -j × P i capcoSF , ( 12
)
with P i caprx m given in ( 8) and P i capcoSF in (9).

B. Imperfect Orthogonality

In reality, spreading factors are not perfectly orthogonal, so all three capture conditions are to be satisfied to achieve a successful transmission. Thus, P success (SF m ) becomes

P success (SF m ) = N ∑ j=1 ( N j ) p j m (1 -p m ) N -j ×P (cap rx , cap coSF , cap intSF ), (13)
where P (cap rx , cap coSF , cap intSF ) is the joint probability for reception condition, capture co-SF and capture inter-SF. 1) For j = 1: the end-device is only subject to inter-SF interferences and the reception condition. As shown in Table I, the inter-SF condition is dominant over the reception one, especially as the number of end-devices increases. Thus, P (cap rx , cap coSF , cap intSF ) = P i capintSF m . The expression of P i capintSF m is similar to P i capcoSF , only with different thresholds and number of interferers. When dealing with inter-SF interferences, one must consider the end-devices within the cell that are not in the annulus corresponding to SF m, i.e., the end-devices with a different spreading factor. If j is the number of end-devices at SF m , then there are N -j end-devices with other spreading factors. Therefore,

P i capintSF m = ∫ lm lm-1 exp ( - q iSFm r α i c ) [ I(r i )] N -j h(r i )dr i . (14) Here, I(r i ) = ∫ R\Rm h(r)dr 1+qcoSF ( r i rp
) α , where R \ R m denotes the whole cell area excluding the area corresponding to SF m . Using (10), we obtain

I(r i ) = J(r i , R)-J(r i , 0)-[J(r i , l m ) -J(r i , l m-1 )] . (15)
2) For j ≥ 2: all capture conditions are to be considered. As in the perfect orthogonality case, the reception condition derives from the two others, thus

P (cap rx , cap coSF , cap intSF ) = P (cap coSF , cap intSF ). (16)
Because of the difficulty to find an exact expression of (16), as one condition might be dominant over the other depending on j, we approximate it as P (cap coSF , cap intSF ) ≈ min(P i capintSF m , P i capcoSF ). Finally, (13) can be written for SF-distance allocation with imperfect orthogonality as,

P success (SF m ) = ( N 1 ) × P i capintSF m + N ∑ j=2 ( N j ) (1 -p m ) N -j min(P i capintSF m , P i capcoSF ), (17)
where P i capintSF m is given in ( 14) and P i capcoSF in (9). Due to lack of space, the analytic throughput expressions for the SF-random case are not given, as they can be obtained similarly as for SF-distance.

IV. NUMERICAL RESULTS

System throughput was evaluated for perfect/imperfect orthogonality as well as for both types of allocations, to assess the validity of our analytical expressions. The main parameters are f c = 868 MHz, BW = 125 kHz and the end-devices' transmit power P 0 = 14 dBm. The path loss exponent was set to α = 4 as in [START_REF] Georgiou | Low Power Wide Area Network Analysis: Can LoRa Scale?[END_REF] (urban) and R = 1 km. Fig. 2 shows throughput performance obtained by simulations and by our theoretical derivations under SF-distance allocation for both perfect and imperfect orthogonality against varying numbers of end-devices transmitting simultaneously. We observe that our derived throughput expressions approach almost perfectly the simulations results with only a small interstice for small values of N (around 5 end-devices) in the imperfect orthogonality case. Therefore, despite approximations, we obtain accurate throughput expressions for both orthogonality cases. Next, we can see that inter-SF interferences cause an early decrease of performance compared to the perfect orthogonality case. However, as the number of devices increases, co-SF interferences always lead to a scalability limit. These results show the impact of imperfect SF orthogonality over the system throughput, up to 50% loss.

Next, the throughput performance for both types of allocations with imperfect orthogonality are shown in Fig. 3. First, we can see that our analysis provides a good approximation of the achievable throughput under SF-random allocation. Then, we notice that for a small amount of end-devices (less than 25), better throughput efficiency is achieved in SF-distance allocation case (up to 100% gain), since devices are more likely to satisfy the specific threshold q SFm . On the other hand, SF-random allocation performs slightly better for greater values of N . Given the higher density of co-SF end-devices under SF-distance, these results suggest that even a simple SFrandom policy provides a higher throughput as the number of end-devices increases. This is because the SF-random allocation favors the case where a limited number of devices randomly choose a small SF, by decreasing their collision probability. Moreover, small SFs lead to larger throughput than large SFs. Thus, our analysis will be useful to devise new allocation policies under various conditions and environments. V. CONCLUSION We have considered the uplink of a single gateway LPWAN based on LoRa physical layer, for which a theoretical throughput expression was derived. Unlike most previous works, our analytical expression encompasses all three conditions required for successful frame transmission: SNR reception level, SINR level for co-SF capture, and SINR level for inter-SF capture. Results have shown the non-negligible impact of SFs' imperfect orthogonality, as well as the drastic effects of SF allocations on the overall throughput. Our analytic framework hence provides a precious tool for designing tailored SF allocations depending on environments and requirements, by predicting their impact on system performance.

Fig. 2 .

 2 Fig. 2. Throughput performance of SF-distance allocation for both perfect and imperfect orthogonality -R = 1km, α = 4

Fig. 3 .

 3 Fig. 3. SF-random and SF-distance allocation throughput performances with imperfect orthogonality -R = 1km, α = 4

This work was supported by the Grant-in-Aid for Scientific Research (Kakenhi) no. 17K06453 from the Ministry of Education, Science, Sports, and Culture of Japan and by the NII Research Grant.