Corentin Pigot
email: corentin.pigot@st.com

Fabien Gilibert

Marina Reyboz

Marc Bocquet

Paola Zuliani

Jean-Michel Portal

Phase-Change Memory: A Continuous Multilevel Compact Model of Subthreshold Conduction and Threshold Switching

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

The rapidly growing market of the internet-of-things requires the use of embedded non-volatile memory (e-NVM) presenting ultra-small area, ultra-fast access time and ultra-low consumption. The mainstream solution, the NOR flash technology, needs high potentials hardly compatible with high-k metal gate of advanced CMOS processes (like FinFET or FD-SOI) and thus requires costly advanced process modifications. Other back-end resistive memories are then investigated, and among them, phase-change memory (PCM) is one of the most mature technologies and has reached a pre-production level. 1,2) The memory principle relies on the phase transition of an active element between two phases (amorphous or crystalline) presenting resistance levels separated by two or three orders of magnitude. 3) The state of the cell is determined by reading the resistance in low field area. The PCM has unique multilevel capabilities, because the resistance can vary continuously between the full RESET and the full SET state. [START_REF] Papandreou | IEEE Int. Conf. Electron. Circuits, Syst[END_REF] For some states including the RESET state, the I-V characteristics of the cell exhibits a threshold switching, above which the amorphous phase becomes suddenly conductive. The nature of this threshold switching has been a long-term discussion and relies classically on two main hypotheses. First of all, the mechanism has been reported to be mainly electronic, [START_REF] Adler | [END_REF][6][7][8][9] but recent studies brought evidences in favor of a thermal activation in some nanoscale cell. 10,11) Although it has nothing to do with phase transition, this behavior is central in the PCM's functioning and crucial for circuit designers to determine sense margin. Indeed, the read operation has to be done under the threshold, otherwise phase transition could happen. However, fast reading means setting the reading voltage as high as possible in order to maximize the current difference between two states and to speed up the read process. Moreover, the lower the SET programing current, and the lower the SET resistance achieved, so the better the programming window, [START_REF] Kluge | IMW[END_REF] which is interesting in terms of multilevel programing. To fully exploit this unique feature, designers need a trustworthy model to verify by simulation the validity of their design. Mainly, designers want to validate that no switching occurs when supplying the array during the read phase and that peripherals are well designed and provide proper biasing during programming phase with the continuous simulation of the state of the cell including the threshold switching process.

In this context, a compact model requires to be fast, robust and accurate. The threshold switching is a difficult part of the PCM modeling due to its intrinsic non-linearity and abrupt transition regarding transient simulation. Thus, it may generate convergence problems in the electrical simulators used in the CAD tools. Lots of compact models of PCM have appeared through the years using various modeling strategies. SPICE macro-models have been developped, [START_REF] El-Hassan | [END_REF][14][15] other more physical models based on a crystalline fraction have been implemented in verilog-A, [16][17][18] but most of them dedicate themselves to the phase transition and attach too few importance to the DC behavior. Among those, some use a negative resistance area, 19) some use a Fermi-like smoothing function, 20,21) others use switches. 22) In this work, it has been modeled for the first time using exclusively self-heating mechanism of the cell. This original approach has been validated through I-V measurements for a large set of intermediate states. The simplicity and the continuity for all regimes (below and above the threshold voltage) of the approach is highly interesting in terms of simulation time and convergence ease required in compact modeling. This paper expands the abstract presented on the 2017 International Conference on Solid State Devices and Materials, 23) justifying deeper the validity of the proposed compact model, and exhibiting new simulation results. First, the measurement setup, followed by the modeling method are presented. The correlation between experimental and modeling results is then detailed, and the good convergence is validated with additional simulations. Finally, comments on the coherency of such modeling approach is discussed and the compliance with a new cell-metric is shown.

Experimental Setup and Modeling Method

Experimentation

Measurements have been performed on a test structure manufactured on a 90nm CMOS node with embedded PCM option. This test structure is composed of a PCM stack serially connected to a MOS transistor, the latter being used to limit the current flowing through the cell.

A TEM cross-section along with a 3-D equivalent schematic of the memory cell is shown on Electrode (TE) and a heater with a wall structure shape. 3) The size of the amorphous dome that can be seen on Fig. 1 reflects the state of the cell, so the goal of the measurements is to let this thickness ua vary in order to highlight threshold switching for all the states where it happens. Tuning the WL voltage from 1V to 2V, a resistance continuum between 125kΩ and 1.3MΩ can be achieved. The current-voltage characteristics is then obtained by reading the Bit Line current while applying a 1V/ms ramp (0 to 2V) on the top electrode. During this read phase, the WL voltage is set to 1.2V to limit the current and thus the PCM stress. In order to avoid any drift effect, 24) and to ensure similar measurement conditions whatever the resistance level, a fixed delay has been introduced between every SET pulse and read ramping.

Compact Model

It is widely known that the amorphous part of the subthreshold transport is a hopping conduction of Poole-Frenkel type. [25][26][START_REF] Shih | ICSSICT[END_REF][START_REF] Ielmini | [END_REF] In this work however, for compact model purpose, a limited density of traps is assumed and only a simplified form 29) is considered, given by,

𝐼 𝑃𝐹 = 𝐴 * 𝐹 * exp (- 𝛷-𝛽√𝐹 𝑘𝑇 ) with 𝐹 = 𝑉 𝑢 𝑎 and 𝛷 = 𝐸 𝑎 0 - 𝑎𝑇 2 𝑏+𝑇 ( 1
)
where k is the Boltzmann constant, β a constant of the material linked to its permittivity and 𝐴 is a fitting parameter. T is a global temperature inside the active area and F the electric field across the amorphous phase. It is calculated through this simplified equation under the assumption of a negligible voltage drop inside the crystalline GST, allowing the access to the amorphous thickness ua, straightly linked to the state of the memory (Fig. 1). V is the PCM's voltage, and Φ is the activation energy of a single coulombic potential well. It follows the Varshni's empirical law for its temperature dependence, 11,30) with 𝐸 𝑎 0 the barrier height at 0K, a and b material-related fitting parameters.

The threshold switching is modeled as a thermal runaway in the Poole-Frenkel current triggered by the self-heating of the cell. Any elevation of the temperature in the material being due to the Joule Effect, the temperature is calculated, under the assumption of a short time constant, as, 𝑇 = 𝑇 𝑎𝑚𝑏 + 𝑅 𝑡ℎ * 𝑃 𝐽 where 𝑃 𝐽 = 𝑉 𝑃𝐶𝑀 * 𝐼 𝑃𝐶𝑀 (2) with Tamb the ambient temperature, Rth an effective thermal resistance, taking amongst other the geometry of the cell into account. 𝑃 𝐽 is the electrical power dissipated inside the PCM.

As it depends on the current flowing through the cell, the calculation of the temperature implies a positive feedback responsible of the switching inside the amorphous phase. Once it has switched, a series resistance of 6kΩ corresponding to the heater resistance limits the current.

Extending the field approximation as long as some amorphous phase exists in the active areaneglecting the voltage drop outside the areathe same Poole-Frenkel current is applied to all the intermediate states as well. ua parameter carries the state as it varies from 0nm to the maximum thickness ua,max extracted from the full RESET state. The crystalline resistance is said to be semiconducting-type, so it can be expressed as, 31) 𝑅 𝑐𝑟𝑦 = 𝑅 𝐶 0 * e -𝐸 𝑎 𝑐 (

1 𝑘𝑇 𝑎𝑚𝑏 - 1 𝑘𝑇 ) (3)
where 𝐸 𝑎 𝑐 is an activation energy and 𝑅 𝐶 0 = 𝑅 𝑐𝑟𝑦 when 𝑇 = 𝑇 𝑎𝑚𝑏 ; they are both treated as fitting parameters.

Results and Discussion

Subthreshold conduction and threshold switching modeling

The comparison of the I-V characteristics between model and measurements for a full range of resistance values is presented on Fig. 4. It shows a very good agreement between data and simulations for two decades of current. The measured resistance is extracted at a constant voltage of 0.36V during the slow ramping procedure. The current at high applied voltage is fitted by the modeling of the serially connected MOS transistor. (b) (a)

800K

The model card parameters are summarized in Table I. Rth is in good agreement with the commonly accepted value for a high thermal efficiency nanoscale PCM cell. 20) A value of relative dielectric constant εr = 10 11) implies, accordingly to Poole-Frenkel's theory, 25) β = 24μeV.V -0.5 .m 0.5 . a and b parameters have been chosen to fit the self-heating inside the GST but they were kept close to the one found in the literature. 11) Similarly, the couple of parameters (RC0, Eac) has been chosen to fit the self-heating of the material in crystalline phase for high current density so it is not surprising that it is found higher than previous values. 32) Fig. 5. Resistance of the cell as a function of the amorphous thickness

The only model parameter that varies from one state to another is the amorphous thickness ua. Reversing Eq. ( 1), amorphous thicknesses have been calculated as a function of the measured resistance for each state. The resistance level as a function of the amorphous thickness is given Fig. 5 and one can verify that there is an excellent correlation between the simulated and the measured-based calculated ua. It means that ua can indeed be used as state parameter for further model-measurement correlation. This allows the computation of the threshold field (cf. Eq. ( 1)), which is plotted Fig. 6, along with the threshold power (see Eq.

(

)), as a function of the resistance of the cell. The threshold is defined as the value of voltage and current where the current in one voltage step of 10mV exceeds a given value of 1μA.

Based on this definition, states that are less resistive than 0.45MΩthat have an amorphous dome smaller than 28.8nmdo not present the threshold switching.

Robustness and coherency of the model

The method used for the measurements and simulations presented Fig. 4 was to apply voltage steps on the top electrode and read the current, it is not possible to see a snapback this way.

On the contrary, the current-driven simulations for different sizes of amorphous dome shown As the threshold switching is highly dependent on the temperature calculation, the subthreshold conduction in full amorphous state (ua = 48nm) as a function of the temperature has been plotted Fig. 8. The threshold power has been extracted based on the same criteria as before, and it is plotted against the ambient temperature in the inset of the Fig. 8. The subthreshold conduction has a strong temperature dependence, but the threshold seems here again to happen at fixed power. These simulations are fully coherent with previous experiments about threshold switching 11) and validate that the switching is indeed triggered by a thermal runaway in the model. I-V characteristics simulated for ambient temperature ranging from 0°C to 85°C. The inset plot the threshold power as a function of the ambient temperature, showing a constant trend. The main roadblock preventing a good multilevel programing of the phase-change memory is for now on the resistance drift, due to relaxation inside the amorphous phase. 2) The resistance of the cell tends to increase with time, preventing to correctly read the state of the cell after a while. To get rid of the inconvenience, Sebastian et al. purposed a new metric M that is less drift-dependent than the low-field resistance for the cell reading. 33) M is defined as the voltage needed to reach a reference current under the application of a linear ramping voltage. M as a function of ua is plotted Fig. 9

Conclusion

This work presents a compact modeling of the threshold switching in phase-change memory based solely on self-heating in the Poole-Frenkel's conduction. This new approach presents the advantage of modeling the current characteristic in a fully continuous way, even the non-linearity of the threshold switching, which eases the convergence and speed-up the simulation time. It has been shown that the model presents a good correlation with measurements

Fig. 1 .

 1 Fig. 1.TEM cross-section (a) and 2D equivalent schematic (b) of the test structure.

Fig. 1 .

 1 Fig. 1. The 50nm-thick phase-change material (GST225) layer has been inserted between Top

Fig. 2

 2 Fig. 2 represents the chronogram of the measurement protocol for each programmed state. A reset pulse of 2V is applied during 200ns before any programming pulse. The word line (WL) bias is then tuned in order to modulate the bit line (BL) current during an 800ns pulse of 2V on the top electrode, resulting in a wide range of intermediate states (see Fig. 3).

Fig. 3 .Fig. 2 .

 32 Fig. 3. Cell resistance as a function of the programming gate voltage, displaying the continuously distributed states achieved

Fig. 4 :

 4 Fig. 4: Cell current versus applied voltage (a) in logarithmic scale and (b) in linear scale for several intermediate states with model (line) and measurements (symbol).

Fig. 7 ,

 7 Fig. 7, exhibit the snapback behavior without any convergence trouble, which illustrates the robustness of the model. The 19.2nm and 48nm simulated curves correspond respectively to the minimum and maximum size of amorphous dome measured. Those values are coherent with the height of the deposited GST layer of about 50nm. The snapback appears for amorphous domes larger than 28.8nm, this limit corresponding to the minimum size of dome where the threshold is observed.

Fig. 6 .

 6 Fig. 6. Threshold field (up) and power (down) versus the resistance of the cell.

Fig. 7 .

 7 Fig. 7. Current-driven simulation for cell states corresponding to the range measured at 273K. The snapback is observable for states that have an amorphous dome larger than 28.8nm.

Fig. 8 .

 8 Fig. 8. I-V characteristics simulated for ambient temperature ranging from 0°C to 85°C. The inset plot the threshold power as a function of the ambient temperature, showing a constant trend.

  for a detection current of 1µA. The simulation presents an excellent correlation with the measurement and both exhibits a linear relationship between M and ua. This proportionality is in agreement with Sebastian's et al., even though the amorphous thickness is not extracted the same way as in article because of the different expression of the subthreshold current used. It confirms the relevance of the choice of ua as a state parameter of the model and validates that the model is suitable for multilevel programing.

Fig. 9 .

 9 Fig. 9. New metric M as a function of the amorphous thickness of the state

Table I .

 I Model card parameters

	Parameter	Value
	Rth	2.0K.μW -1
	A	1.45.10 -4 Ω -1
	β	24µeV.V -0.5 .m 0.5
	𝐸 𝑎 0	0.3eV
	Rheater	6kΩ
	𝑅 𝐶 0	10kΩ
	𝐸 𝑎 𝑐	0.1eV
	a	1.2meV.K -1
	b


Luc Le Magoarou

Stéphane Paquelet

PARAMETRIC CHANNEL ESTIMATION FOR MASSIVE MIMO

Keywords: Cramér-Rao bound, Channel estimation, MIMO

Channel state information is crucial to achieving the capacity of multiantenna (MIMO) wireless communication systems. It requires estimating the channel matrix. This estimation task is studied, considering a sparse physical channel model, as well as a general measurement model taking into account hybrid architectures. The contribution is twofold. First, the Cramér-Rao bound in this context is derived. Second, interpretation of the Fisher Information Matrix structure allows to assess the role of system parameters, as well as to propose asymptotically optimal and computationally efficient estimation algorithms.

INTRODUCTION

Multiple-Input Multiple-Output (MIMO) wireless communication systems allow for a dramatic increase in channel capacity, by adding the spatial dimension to the classical time and frequency ones [START_REF] Telatar | Capacity of multi-antenna gaussian channels[END_REF][START_REF] Tse | Fundamentals of wireless communication[END_REF]. This is done by sampling space with several antenna elements, forming antenna arrays both at the transmitter (with nt antennas) and receiver (with nr antennas). Capacity gains over single antenna systems are at most proportional to min(nr,nt).

Millimeter wavelengths have recently appeared as a viable solution for the fifth generation (5G) wireless communication systems [START_REF] Theodore S Rappaport | Millimeter wave mobile communications for 5g cellular: It will work![END_REF][START_REF] Swindlehurst | Millimeter-wave massive mimo: the next wireless revolution?[END_REF]. Indeed, smaller wavelengths allow to densify half-wavelength separated antennas, resulting in higher angular resolution and capacity for a given array size. This observation has given rise to the massive MIMO field, i.e. the study of systems with up to hundreds or even thousands of antennas.

Massive MIMO systems are very promising in terms of capacity. However, they pose several challenges to the research community [START_REF] Rusek | Scaling up mimo: Opportunities and challenges with very large arrays[END_REF][START_REF] Erik G Larsson | Massive mimo for next generation wireless systems[END_REF], in particular for channel estimation. Indeed, maximal capacity gains are obtained in the case of perfect knowledge of the channel state by both the transmitter and the receiver. The estimation task amounts to determine a complex gain between each transmit/receive antenna pair, the narrowband (single carrier) MIMO channel as a whole being usually represented as a complex matrix H ∈ C nr ×nt of such complex gains. Without a parametric model, the number of real parameters to estimate is thus 2nrnt, which is very large for massive MIMO systems. Contributions and organization. In this work, massive MIMO channel estimation is studied, and its performance limits are sought, as well as their dependency on key system parameters. In order to answer this question, the framework of parametric estimation [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF] is used. A physical channel model is first presented, with the general considered observation model, and the objective is precisely stated. The Cramér-Rao bound for is then derived, which bounds the variance of any unbiased estimator. Then, the interpretation of the bound allows to precisely assess the role of system design on estimation performance, as well as to propose new computationally efficient channel estimation algorithms showing asymptotic performance equivalent to classical ones based on sparse recovery.

PROBLEM FORMULATION

Notations. Matrices and vectors are denoted by bold upper-case and lower-case letters: A and a (except 3D "spatial" vectors that are denoted -→ a ); the ith column of a matrix A by: ai; its entry at the ith line and jth column by: aij or Aij. A matrix transpose, conjugate and transconjugate is denoted by: A T , A * and A H respectively. The image, rank and trace of a linear transformation represented by A are denoted: im(A), rank(A) and Tr(A) respectively. For matrices A and B, A ≥ B means that A-B is positive semidefinite. The linear span of a set of vectors A is denoted: span(A). The Kronecker product, standard vectorization and diagonalization operators are denoted by vec(•), diag(•), and ⊗ respectively. The identity matrix, the m×n matrix of zeros and ones are denoted by Id, 0m×n and 1m×n respectively. CN (µ,Σ) denotes the standard complex gaussian distribution with mean µ and covariance Σ. E(.) denotes expectation and cov(.) the covariance of its argument.

Parametric physical channel model

Consider a narrowband block fading channel between a transmitter and a receiver with respectively nt and nr antennas. It is represented by the matrix H ∈ C nr ×nt , in which hij corresponds to the channel between the jth transmit and ith receive antennas.

Classically, for MIMO systems with few antennas, i.e. when the quantity nrnt is small (up to a few dozens), estimators such as the Least Squares (LS) or the Linear Minimum Mean Squared Error (LMMSE) are used [START_REF] Biguesh | Training-based mimo channel estimation: a study of estimator tradeoffs and optimal training signals[END_REF].

However, for massive MIMO systems, the quantity 2nr nt is large (typically several hundreds), and resorting to classical estimators may become computationally intractable. In that case, a parametric model may be used. Establishing it consists in defining a set of np parameters θ (θ1,...,θn p ) T that describe the channel as H ≈ f (θ) for a given function f , where the approximation is inherent to the model structure and neglected in the sequel (considering H = f (θ)). Channel estimation then amounts to estimate the parameters θ instead of the channel matrix H directly. The parametrization is particularly useful if np ≪ 2nrnt, without harming accuracy of the channel description. Inspired by the physics of wave propagation under the plane waves assumption, it has been proposed to express the channel matrix as a sum of rank-1 matrices, each corresponding to a single physical path between transmitter and receiver [START_REF] Akbar | Deconstructing multiantenna fading channels[END_REF]. Adopting this kind of modeling and generalizing it to take into account any three-dimensional antenna array geometry, channel matrices take the form

H = P p=1 cper( --→ ur,p).et( --→ ut,p) H , (1)
where P is the total number of considered paths (no more than a few dozens), cp ρpe jφp is the complex gain of the pth path, --→ ut,p is the unit vector corresponding to its Direction of Departure (DoD) and --→ ur,p the unit vector corresponding to its Direction of Arrival (DoA). Any unit vector -→ u is described in spherical coordinates by an azimuth angle η and an elevation angle ψ. The complex response and steering vectors er( -→ u ) ∈ C nr and et( -→ u ) ∈ C nt are defined as (ex( -→ u ))i = 1

√ nx e -j 2π λ --→ a x,i . -→ u for x ∈ {r,t}. The set { --→ ax,1,..., ---→ ax,n x } gathers the positions of the antennas with respect to the centroid of the considered array (transmit if x = t, receive if x = r). In order to lighten notations, the matrix Ax 2π λ ( --→ ax,1,... , ---→ ax,n x ) ∈ R 3×nx is introduced. It simplifies the steering/response vector expression to ex( -→ u ) = 1

√ nx e -jA T x -→ u , where the exponential function is applied component-wise. In order to further lighten notations, the pth atomic channel is defined as Hp cper( --→ ur,p).et( --→ ut,p) H , and its vectorized version hp vec(Hp) ∈ C nr nt . Therefore, defining the vectorized channel h vec(H), yields h = P p=1 hp. Note that the channel description used here is very general, as it handles any three-dimensional antenna array geometry, not only Uniform Linear Arrays (ULA) or Uniform Planar Arrays (UPA) as is sometimes proposed.

In short, the physical channel model can be seen as a parametric model with θ = {θ (p) (ρp,φp,ηr,p,ψr,p,ηt,p,ψt,p), p = 1,...,P }. There are thus 6P real parameters in this model (the complex gain, DoD and DoA of every path are described with two parameters each). Of course, the model is most useful for estimation in the case where 6P ≪ 2nrnt, since the number of parameters is thus greatly reduced.

Note that most classical massive MIMO channel estimation methods assume a similar physical model, but discretize a priori the DoDs and DoAs, so that the problem fits the framework of sparse recovery [START_REF] Mallat | Matching pursuits with timefrequency dictionaries[END_REF][START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF][START_REF] Bajwa | Compressed channel sensing: A new approach to estimating sparse multipath channels[END_REF]. The approach used here is different, in the sense that no discretization is assumed for the analysis.

Observation model

In order to carry out channel estimation, ns known pilot symbols are sent through the channel by each transmit antenna. The corresponding training matrix is denoted X ∈ C nt×ns . The signal at the receive antennas is thus expressed as HX + N, where N is a noise matrix with vec(N) ∼ CN (0,σ 2 Id). Due to the high cost and power consumption of millimeter wave Radio Frequency (RF) chains, it has been proposed to have less RF chains than antennas in both the transmitter and receiver [START_REF] Ayach | Spatially sparse precoding in millimeter wave mimo systems[END_REF][START_REF] Alkhateeb | Channel estimation and hybrid precoding for millimeter wave cellular systems[END_REF][START_REF] Heath | An overview of signal processing techniques for millimeter wave mimo systems[END_REF][START_REF] Akbar | Millimeter-Wave MIMO Transceivers: Theory, Design and Implementation[END_REF]. Such systems are often referred to as hybrid architectures. Mathematically speaking, this translates into specific constraints on the training matrix X (which has to "sense" the channel through analog precoders vi ∈ C nt , i = 1,...,nRF, nRF being the number of RF chains on the transmit side), as well as observing the signal at the receiver through analog combiners. Let us denote wj ∈ C nr , j = 1,...,nc the used analog combiners, the observed data is thus expressed in all generality as

Y = W H HX+W H N, (2)
where W (w1,... ,wn c ) and the training matrix is constrained to be of the form X = VZ, where Z ∈ C n RF ×ns is the digital training matrix.

Objective: bounding the variance of unbiased estimators

In order to assess the fundamental performance limits of channel estimation, the considered performance measure is the relative Mean Squared Error (rMSE). Denoting indifferently H(θ) f (θ) or H the true channel (h(θ) or h in vectorized form) and H( θ) f ( θ) or Ĥ its estimate (h( θ) or ĥ in vectorized form) in order to lighten notations, rMSE is expressed

rMSE = E H-Ĥ 2 F . H -2 F = Tr cov ĥ Variance + E( Ĥ)-H 2 F Bias . H -2 F , (3)
where the bias/variance decomposition can be done independently of the considered model [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF]. The goal here is to lower-bound the variance term, considering the physical model introduced in the previous subsection. The bias term is not studied in details here, but its role is evoked in section 3.3.

CRAM ÉR-RAO LOWER BOUND

In this section, the variance term of eq. ( 3) is bounded using the Cramér-Rao Bound (CRB) [START_REF] Calyampudi | Information and the accuracy attainable in the estimation of statistical parameters[END_REF]18], which is valid for any unbiased estimator θ of the true parameter θ. The complex CRB [START_REF] Van Den Bos | A cramér-rao lower bound for complex parameters[END_REF] states,

cov g( θ) ≥ ∂g(θ) ∂θ I(θ) -1 ∂g(θ) ∂θ H , with I(θ) E ∂logL ∂θ ∂logL ∂θ H the Fisher Information Matrix (FIM),
where L denotes the model likelihood, and g is any complex differentiable vector function. In particular, regarding the variance term of eq. ( 3),

Tr cov h( θ) ≥ Tr ∂h(θ) ∂θ I(θ) -1 ∂h(θ) ∂θ H , (4)
with ∂h(θ) ∂θ = ∂h(θ) ∂θ 1 ,..., ∂h(θ) ∂θn p

. A model independent expression for the FIM is provided in section 3.1, and particularized in section 3.2 to the model of section: 2.1. Finally, the bound is derived from eq. ( 4) in section 3.3.

General derivation

First, notice that vectorizing eq. ( 2), the observation matrix Y follows a complex gaussian distribution,

vec(Y) ∼ CN (X T ⊗W H )h(θ) µ(θ) ,σ 2 (Idn s ⊗W H W) Σ .
In that particular case, the Slepian-Bangs formula [START_REF] Slepian | Estimation of signal parameters in the presence of noise[END_REF][START_REF] Bangs | Array Processing With Generalized Beamformers[END_REF] yields:

I(θ) = 2Re ∂µ(θ) ∂θ H Σ -1 ∂µ(θ) ∂θ = 2α 2 σ 2 Re ∂h(θ) ∂θ H P ∂h(θ) ∂θ , (5)
with

P σ 2 α 2 (X * ⊗W)Σ -1 (X T ⊗W H ) where α 2 1
ns Tr(X H X) is the average transmit power per time step. Note that the expression can be simplified to P = 1 α 2 (X * X T ) ⊗ (W(W H W) -1 W H ) using elementary properties of the Kronecker product. The matrix W(W H W) -1 W H is a projection matrix onto the range of W. In order to ease further interpretation, assume that X H X = α 2 Idn s . This assumption means that the transmit power is constant during training time ( xi 2 2 = α 2 , ∀i) and that pilots sent at different time instants are mutually orthogonal (x H i xj = 0, ∀i = j). This way, 1 α 2 X * X T is a projection matrix onto the range of X * , and P can itself be interpreted as a projection, being the Kronecker product of two projection matrices [22, p.112] (it is an orthogonal projection since P H = P).

Fisher information matrix for a sparse channel model

Consider now the parametric channel model of section 2.1, where h = P p=1 hp, with hp = cpet( --→ ut,p) * ⊗er( --→ ur,p). Intra-path couplings. The derivatives of h with respect to parameters of the pth path θ (p) can be determined using matrix differentiation rules [START_REF] Brandt Petersen | The matrix cookbook[END_REF]: ∂φp , ∂h(θ) ∂ηr,p , ∂h(θ) ∂ψr,p , ∂h(θ) ∂ηt,p , ∂h(θ) ∂ψt,p , the part of the FIM corresponding to couplings between the parameters θ (p) (intra-path couplings) is expressed as

•
I (p,p) 2α 2 σ 2 Re ∂h H ∂θ (p) P ∂h ∂θ (p) . (6)
Let us now particularize this expression. First of all, in order to ease interpretations carried out in section 4, consider the case of optimal observation conditions (when the range of P contains the range of ∂h(θ) ∂θ ). This allows indeed to interpret separately the role of the observation matrices and the antenna arrays geometries. Second, consider for example the entry corresponding to the coupling between the departure azimuth angle ηt,p and the arrival azimuth angle ηr,p of the pth path. It is expressed under the optimal observation assumption as since Ar1n r = 0 and At1n t = 0 by construction (because the antennas positions are taken with respect to the array centroid). This means that the parameters ηr,p and ηt,p are statistically uncoupled, i.e. orthogonal parameters [START_REF] David | Parameter orthogonality and approximate conditional inference[END_REF]. Computing all couplings for θ (p) yields

I (p,p) = 2ρ 2 p α 2 σ 2      1 ρ 2 p 0 01×2 01×2 0 1 01×2 01×2 02×1 02×1 Br 02×2 02×1 02×1 02×2 Bt      , (7)
where

Bx = 1 nx   A T x --→ vη x,p 2 2 --→ vη x,p T AxA T x ---→ v ψx,p ---→ v ψx,p T AxA T x --→ vη x,p A T x ---→ v ψx,p 2 2   , (8)
with x ∈ {r, t}. These expressions are thoroughly interpreted in section 4.

Global FIM. Taking into account couplings between all paths, The global FIM is easily deduced from the previous calculations and block structured, 2) ... I (1,P ) I (2,1) I (2,2) . . . . . . I (P,1) I (P,P )   , where I (p,q) ∈ R 6×6 contains the couplings between parameters of the pth and qth paths and is expressed I (p,q) 2α 2 σ 2 Re ∂h ∂θ (p) H P ∂h ∂θ (q) . The off-diagonal blocks I (p,q) of I(θ), corresponding to couplings between parameters of distinct paths, or inter-path couplings, can be expressed explicitly (as in eq. ( 7) for intra-path couplings). However, the obtained expressions are less prone to interesting interpretations, and inter-paths couplings have been observed to be negligible in most cases. They are thus not displayed in the present paper, for brevity reasons. Note that a similar FIM computation was recently carried out in the particular case of linear arrays [START_REF] Garcia | Optimal robust precoders for tracking the aod and aoa of a mm-wave path[END_REF]. However, the form of the FIM (in particular parameter orthogonality) was not exploited in [START_REF] Garcia | Optimal robust precoders for tracking the aod and aoa of a mm-wave path[END_REF], as is done here in sections 4 and 5.

I(θ) =   I (1,1) I (1,

Bound on the variance

The variance of channel estimators remains to be bounded, using eq. ( 4). From eq. ( 5), the FIM can be expressed more conveniently only with real matrices as

I(θ) = 2α 2 σ 2 DT P D, with D Re{ ∂h(θ) ∂θ } Im{ ∂h(θ) ∂θ } , P Re{P} -Im{P} Im{P} Re{P} ,
where P is also a projection matrix. Finally, injecting eq. ( 5) into eq. ( 4) assuming the FIM is invertible, gives for the relative variance (this is actually an optimal SNR, only attained with perfect precoding and combining). Optimal bound. The first inequality in eq. ( 9) becomes an equality if an efficient estimator is used [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF]. Moreover, the second inequality is an equality if the condition im ∂h(θ) ∂θ ⊂ im (P) is fulfilled (this corresponds to optimal observations, further discussed in section 4). Remarkably, under optimal observations, the lower bound on the relative variance is directly proportional to the considered number of paths P and inversely proportional to the SNR, and does not depend on the specific model structure, since the influence of the derivative matrix D cancels out in the derivation. Sparse recovery CRB. It is interesting to notice that the bound obtained here is similar to the CRB for sparse recovery [START_REF] Ben | The cramér-rao bound for estimating a sparse parameter vector[END_REF] (corresponding to an intrinsically discrete model), that is proportional to the sparsity of the estimated vector, analogous here to the number of paths.

Tr cov h( θ) . h -2 2 ≥ σ 2 2α 2 Tr D( DT P D) -1 DT . h -2 2 ≥ σ 2 2α 2 Tr D( DT D) -1 DT . h -2 2 = σ 2 2α 2 h 2 2 np = 3P SNR , (9)

INTERPRETATIONS

The main results of sections 3.2 and 3.3 are interpreted in this section, ultimately guiding the design of efficient estimation algorithms. Parameterization choice. The particular expression of the FIM allows to assess precisely the chosen parameterization. First of all, I(θ) has to be invertible and well-conditioned, for the model to be theoretically and practically identifiable [START_REF] Thomas | Identification in parametric models[END_REF][START_REF] Kravaris | Advances and selected recent developments in state and parameter estimation[END_REF], respectively. As a counterexample, imagine two paths indexed by p and q share the same DoD and DoA, then proportional columns appear in ∂h(θ) ∂θ , which implies non-invertibility of the FIM. However, it is possible to summarize the effect of these two paths with a single virtual path of complex gain cp +cq without any accuracy loss in channel description, yielding an invertible FIM. Similarly, two paths with very close DoD and DoA yield an ill-conditioned FIM (since the corresponding steering vectors are close to colinear), but can be merged into a single virtual path with a limited accuracy loss, improving the conditioning. Interestingly, in most channel models, paths are assumed to be grouped into clusters, in which all DoDs and DoAs are close to a principal direction [START_REF] Adel | A statistical model for indoor multipath propagation[END_REF][START_REF] Jensen | Modeling the indoor mimo wireless channel[END_REF][START_REF] Michael | A review of antennas and propagation for mimo wireless communications[END_REF]. Considering the MSE, merging close paths indeed decreases the variance term (lowering the total number of parameters), without increasing significantly the bias term (because their effects on the channel matrix are very correlated). These considerations suggest dissociating the number of paths considered in the model P from the number of physical paths, denoted P φ , taking P < P φ by merging paths. This is one motivation behind the famous virtual channel representation [START_REF] Akbar | Deconstructing multiantenna fading channels[END_REF], where the resolution at which paths are merged is fixed and given by the number of antennas. The theoretical framework of this paper suggests to set P (and thus the merging resolution) so as to minimize the MSE. A theoretical study of the bias term of the MSE (which should decrease when P increases) could thus allow to calibrate models, choosing an optimal number of paths P * for estimation. Such a quest for P * is carried out empirically in section 5.

Optimal observations. The matrices X and W (pilot symbols and analog combiners) determine the quality of channel observation. Indeed, it was shown in section 3.3 that the lowest CRB is obtained when im ∂h(θ) ∂θ ⊂ im (P), with

P = 1 α 2 (X * X T ) ⊗ (W(W H W) -1 W H ) .
In case of sparse channel model, using the expressions for ∂h (θ) ∂θ derived above, this is equivalent to two distinct conditions for the training matrix: ux,p) with x ∈ {r, t} and ξ ∈ {η, ψ}. These conditions are fairly intuitive: to estimate accurately parameters corresponding to a given DoD (respectively DoA), the sent pilot sequence (respectively analog combiners) should span the corresponding steering vector and its derivatives (to "sense" small changes). To accurately estimate all the channel parameters, it should be met for each atomic channel. Array geometry. Under optimal observation conditions, performance limits on DoD/DoA estimation are given by eq. ( 8). The lower the diagonal entries B -1

x , the better the bound. This implies the bound is better if the diagonal entries of Bx are large and the offdiagonal entries are small (in absolute value). Since the unit vectors --→ vη x,p and ---→ v ψx,p are by definition orthogonal, having AxA T x = β 2 Id with maximal β 2 is optimal, and yields uniform performance limits for any DoD/DoA. Moreover, in this situation, β 2 is proportional to

1 nx nx i=1 --→ ax,i 2 2
, the mean squared norm of antenna positions with respect to the array centroid. Having a larger antenna array is thus beneficial (as expected), because the furthest antennas are from the array centroid, the larger β 2 is. Orthogonality of DoA and DoD. Section 3.2 shows that the matrix corresponding to intra-path couplings (eq. ( 7)) is block diagonal, meaning that for a given path, parameters corresponding to gain, Algorithm 1 Sequential direction estimation (DoA first)

X H et( -→ v 1 ) X H et( -→ v 1 ) 2 |...| X H et( -→ vn) X H et( -→ vn) 2
6: Find the index ĵ of the maximal entry of er(

-→ u î ) H YKt, set - → ut ← -→ v ĵ (O(n) complexity)
phase, DoD and DoA are mutually orthogonal. Maximum Likelihood (ML) estimators of orthogonal parameters are asymptotically independent [START_REF] David | Parameter orthogonality and approximate conditional inference[END_REF] (when the number of observations, or equivalently the SNR goes to infinity). Classically, channel estimation in massive MIMO systems is done using greedy sparse recovery algorithms [START_REF] Mallat | Matching pursuits with timefrequency dictionaries[END_REF][START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF][START_REF] Bajwa | Compressed channel sensing: A new approach to estimating sparse multipath channels[END_REF]. Such algorithms can be cast into ML estimation with discretized directions, in which the DoD and DoA (coefficient support) are estimated jointly first (which is costly), and then the gain and phase are deduced (coefficient value), iteratively for each path. Orthogonality between the DoD and DoA parameters is thus not exploited by classical channel estimation methods. We propose here to exploit it via a sequential decoupled DoD/DoA estimation, that can be inserted in any sparse recovery algorithm in place of the support estimation step, without loss of optimality in the ML sense. In the proposed method, one direction (DoD or DoA) is estimated first using an ML criterion considering the other direction as a nuisance parameter, and the other one is deduced using the joint ML criterion. Such a strategy is presented in algorithm 1. It can be verified that lines 3 and 6 of the algorithm actually correspond to ML estimation of the DoA and joint ML estimation, respectively. The overall complexity of the sequential directions estimation is thus O(m+n), compared to O(mn) for the joint estimation with the same test directions. Note that a similar approach, in which DoAs for all paths are estimated at once first, was recently proposed [START_REF] Noureddine | A two-step compressed sensing based channel estimation solution for millimeter wave mimo systems[END_REF] (without theoretical justification).

PRELIMINARY EXPERIMENT

Let us compare the proposed sequential direction estimation to the classical joint estimation. This experiment must be seen as an example illustrating the potential of the approach, and not as an extensive experimental validation. Experimental settings. Consider synthetic channels generated using the NYUSIM channel simulator [START_REF] Mathew | 3-d millimeterwave statistical channel model for 5g wireless system design[END_REF] (setting f = 28 GHz, the distance between transmitter and receiver to d = 30 m) to obtain the DoDs, DoAs, gains and phases of each path. The channel matrix is then obtained from eq. ( 1), considering square Uniform Planar Arrays (UPAs) with half-wavelength separated antennas, with nt = 64 and nr = 16. Optimal observations are considered, taking both W and X as the identity. Moreover, the noise variance σ 2 is set so as to get an SNR of 10 dB. Finally, the two aforementioned direction estimation strategies are inserted in the Matching Pursuit (MP) algorithm [START_REF] Mallat | Matching pursuits with timefrequency dictionaries[END_REF], discretizing the directions taking m = n = 2, 500, and varying the total number P of estimated paths. Results. Table 1 shows the obtained relative MSE and estimation times (Python implementation on a laptop with an Intel(R) Core(TM) i7-3740QM CPU @ 2.70 GHz). First of all, for P = 5, 10, 20, the estimation error decreases and the estimation time increases with P , exhibiting a trade-off between accuracy and time. However, increasing P beyond a certain point seems useless, since the error re-increases, as shown by the MSE for P = 40, echoing the trade-off evoked in section 3.3, and indicating that P * is certainly between 20 and 40 for both methods in this setting. Finally, for any value of P , while the relative errors of the sequential and joint estimation methods are very similar, the estimation time is much lower (between ten and twenty times) for sequential estimation. This observation validates experimentally the theoretical claims made in the previous section.

CONCLUSIONS AND PERSPECTIVES

In this paper, the performance limits of massive MIMO channel estimation were studied. To this end, training based estimation with a physical channel model and an hybrid architecture was considered. The Fisher Information Matrix and the Cramér-Rao bound were derived, yielding several results. The CRB ended up being proportional to the number of parameters in the model and independent from the precise model structure. The FIM allowed to draw several conclusions regarding the observation matrices and the arrays geometries. Moreover, it suggested computationally efficient algorithm which are asymptotically as accurate as classical ones.

This paper is obviously only a first step toward a deep theoretical understanding of massive MIMO channel estimation. Apart from more extensive experimental evaluations and optimized algorithms, a theoretical study of the bias term of the MSE would be needed to calibrate models, and the interpretations of section 4 could be leveraged to guide system design.

  Regarding the complex gain cp = ρpe jφp , the model yields the expressions ∂h(θ) ∂ρp = 1 ρp hp and ∂h(θ) ∂φp = jhp. • Regarding the DoA, ∂h(θ) ∂ηr,p = Idn t ⊗diag(-jA T r --→ vη r,p ) hp and ∂h(θ) ∂ψr,p = Idn t ⊗diag(-jA T r ---→ v ψr,p ) hp, where --→ vη r,p and ---→ v ψr,p are the unit vectors in the azimuth and elevation directions at --→ ur,p, respectively. • Regarding the DoD, ∂h(θ) ∂ηt,p = diag(jA T t --→ vη t,p )⊗Idn r hp and ∂h(θ) ∂ψt,p = diag(jA T t --→ v ψt,p )⊗Idn r hp, where --→ vη t,p and --→ v ψt,p are the unit vectors in the azimuth and elevation directions at --→ ut,p, respectively. Denoting ∂h ∂θ (p) ∂h(θ) ∂ρp , ∂h(θ)

2α 2 σ 2

 2 Re ∂h(θ) vη r,p = 0,

2 h 2 2 σ 2

 222 where the second inequality comes from the fact that P being an orthogonal projection matrix, P ≤ Id ⇒ DT P D ≤ DT D ⇒ ( DT P D) -1 ≥ ( DT D) -1 ⇒ D( DT P D) -1 DT ≥ D( DT D) -1 DT (using elementary properties of the ordering of semidefinite positive matrices, in particular [26, Theorem 4.3]). The first equality comes from the fact that Tr D( DT D) -1 DT = Tr(Idn p ) = np. Finally, the second equality is justified by np = 6P considering the sparse channel model, and by taking SNRα

=

  diag(-jA T x --→ v ξx,p )ex( --→

Acknowledgments. The authors wish to thank Matthieu Crussière for the fruitful discussions that greatly helped improving this work.

This work has been performed in the framework of the Horizon 2020 project ONE5G (ICT-760809) receiving funds from the European Union. The authors would like to acknowledge the contributions of their colleagues in the project, although the views expressed in this contribution are those of the authors and do not necessarily represent the project.


Quentin Bramas
email: bramas@unistra.fr

The Stability and the Security of the Tangle

In this paper we study the stability and the security of the distributed data structure at the base of the IOTA protocol, called the Tangle. The contribution of this paper is twofold. First, we present a simple model to analyze the Tangle and give the first discrete time formal analyzes of the average number of unconfirmed transactions and the average confirmation time of a transaction.

Then, we define the notion of assiduous honest majority that captures the fact that the honest nodes have more hashing power than the adversarial nodes and that all this hashing power is constantly used to create transactions. This notion is important because we prove that it is a necessary assumption to protect the Tangle against double-spending attacks, and this is true for any tip selection algorithm (which is a fundamental building blocks of the protocol) that verifies some reasonable assumptions. In particular, the same is true with the Markov Chain Monte Carlo selection tip algorithm currently used in the IOTA protocol.

Our work shows that either all the honest nodes must constantly use all their hashing power to validate the main chain (similarly to the bitcoin protocol) or some kind of authority must be provided to avoid this kind of attack (like in the current version of the IOTA where a coordinator is used).

The work presented here constitute a theoretical analysis and cannot be used to attack the current IOTA implementation. The goal of this paper is to present a formalization of the protocol and, as a starting point, to prove that some assumptions are necessary in order to defend the system again double-spending attacks. We hope that it will be used to improve the current protocol with a more formal approach.

Introduction

Since the day Satoshi Nakamoto presented the Bitcoin protocol in 2008 [START_REF] Nakamoto | Bitcoin: A peer-to-peer electronic cash system[END_REF], the interest in Blockchain technology has grown continuously. More generally, the interest concerns Distributed Ledger Technology, which refers to a distributed data storage protocol. Usually it involves a number of nodes (or processes, or agents) in a network that are known to each other or not. Those nodes may not trust each-other so the protocol should ensure that they reach a consensus on the order of the operations they perform, in addition to other mechanism like data replication for instance.

The consensus problem has been studied for a long time [START_REF] Pease | Reaching agreement in the presence of faults[END_REF][START_REF] Fischer | Impossibility of distributed consensus with one faulty process[END_REF] providing a number of fundamental results. But the solvability of the problem was given in term of proportion of faulty agents over honest agents and in a trustless network, where anyone can participate, an adversary can simulate an arbitrary number of nodes in the network. To avoid that, proof systems like Proof of Work (PoW) or Proof of Stake (PoS) are used to link the importance of an entity with some external properties (processing power in PoW) or internal properties (the number of owned tokens 1 in PoS) instead of simply the number of nodes it controls. The solvability of the consensus is now possible only if the importance of the adversary (given in terms of hashing power or in stake) is smaller than the honest one (the proportion is reduced to 1/3 if the network is asynchronous).

In Bitcoin and in the other blockchain technologies, transactions are stored in a chain of blocks, and the PoW or PoS is used to elect one node that is responsible for writing data in the next block. The "random" selection and the incentive a node may have to execute honestly the protocol make the whole system secure, as it was shown by several formal analysis [START_REF] Nakamoto | Bitcoin: A peer-to-peer electronic cash system[END_REF][START_REF] Garay | The bitcoin backbone protocol: Analysis and applications[END_REF]. Usually, there are properties that hold with high probability i.e., with a probability that tends to one quickly as the time increases. For instance, the order between two transactions do not change with probability that tends to 1 exponentially fast over the time in the Bitcoin protocol, if the nodes executing honestly (or rationally) the protocol have more than a third of the hashing total power.

In this paper we study another distributed ledger protocol called the Tangle, presented by Serguei Popov [START_REF] Popov | The tangle. white paper[END_REF] that is used in the IOTA cryptocurrency to store transactions. The Tangle is a Directed Acyclic Graph (DAG) where a vertex, representing a transaction, has two parents, representing the transactions it confirms.

According to the protocol a PoW has to be done when adding a transaction to the Tangle. This PoW prevents an adversary from spamming the network. However, it is not clear in the definition of the Tangle how this PoW impacts its security.

When a new transaction is appended to the Tangle, it references two previous unconfirmed transactions, called tips. The algorithm selecting the two tips is called a Tip Selection Algorithm (TSA). It is a fundamental parts of the protocol as it is used by the participants to decide, among two conflicting transactions, which one is the valid. It is the most important part in order for the participants to reach a consensus. The TSA currently used in the IOTA implementation uses the PoW contained in each transaction to select the two tips.

Related Work Very few academic papers exist on this protocol, and there is no previous work that formally analyzes its security. The white paper behind the Tangle [START_REF] Popov | The tangle. white paper[END_REF] presents a quick analysis of the average number of transactions in the continuous time setting. This analysis is done after assuming that the evolution converge toward a stationary distribution. The same paper presents a TSA using Monte Carlo Markov Chain (MCMC) random walks in the DAG from old transactions toward new ones, to select two unconfirmed transactions. The random walk is weighted to favor transactions that are confirmed by more transactions. There is no analysis on how the assigned weight, based on the PoW of each transaction, affects the security of the protocol. This MCMC TSA is currently used by the IOTA cryptocurrency.

It is shown in [START_REF] Popov | Equilibria in the tangle[END_REF] that using the default TSA correspond to a Nash equilibrium. Participants are incite to use the MCMC TSA, because using another TSA (e.g. a lazy one that confirms only already confirmed transactions) may increase the chances of seeing their transactions unconfirmed.

Finally, the tangle has also been analyzed by simulation [START_REF] Kuśmierz | The first glance at the simulation of the tangle: discrete model[END_REF] using a discrete time model, where transactions are issued every round following a Poisson distribution of parameter λ. Like in the continuous time model, the average number of unconfirmed transactions (called tips) seems to grow linearly with the value of λ, but a little bit slower (≈ 1.26λ compared to 2λ in the continuous time setting) Contributions The contribution of our paper is twofold. First, we analyze formally the number of tips in the discrete time setting, depending on the value of λ by seeing it as a Markov chain where at each round, their is a given probability to obtain a given number of tips. Unlike previous work, we here prove the convergence of the system toward a stationary distribution. This allows use to prove the previous results found by simulations [START_REF] Kuśmierz | The first glance at the simulation of the tangle: discrete model[END_REF] that the average number of tips is stationary and converge towards a fixed value.

Second, we prove that if the TSA depends only on the PoW, then the weight of the honest transactions should exceed the hashing power of the adversary to prevent a double-spending attack. This means that honest nodes should constantly use their hashing power and issue new transactions, otherwise an adversary can attack the protocol even with a small fraction of the total hashing power. Our results is interesting because it is true for any tip selection algorithm i.e., the protocol cannot be more secure by simply using a more complex TSA.

The remaining of the paper is organized as follow. Section 2 presents our model and the Tangle. In Section 3 we analyze the average confirmation time and the average number of unconfirmed transactions. In Section 4 we prove our main theorem by presenting a simple double-spending attack.

Model

The Network

We consider a set N of processes, called nodes, that are fully connected. Each node can send a message to all the other nodes (the network topology is a complete graph). We will discuss later the case where the topology is sparse.

We assume nodes are activated synchronously. The time is discrete and at each time instant, called round, a node reads the messages sent by the other nodes in the previous round, executes the protocol and, if needed, broadcast a message to all the other nodes. When a node broadcast a message, all the other nodes receive it in the next round. The size of messages are in O(B), where B is an upper bound on the size of the messages when all the nodes follow the protocol honestly i.e., the size of the messages is sufficiently high not to impact the protocol when it is executed honestly. Those assumptions are deliberately strong as they make an attack more difficult to perform.

The DAG

In this paper, we consider a particular kind of distributed ledger called the Tangle, which is a Direct Acyclic Graph (DAG). Each node u stores at a given round r a local DAG G u r (or simply G r or G if the node or the round are clear from the context), where each vertex, called site, represents a transaction. Each site has two parents (possibly the same) in the DAG. We say a site directly confirm its two parents. All sites that are confirmed by the parents of a site are also said to be confirmed (or indirectly confirmed) by it i.e., there is a path from a site to all the sites it confirms in the DAG (see Figure 1). A site that is not yet confirmed is called a tip. There is a unique site called genesis that does not have parents and is confirmed by all the other sites. For simplicity we model a DAG simply by a set of sites G = (s i ) i∈I Two sites may be conflicting. This definition is application-dependent so we assume that there exists a function areConf licting(a, b) that answer whether two sites are conflicting or not.

If the Tangle is used to store the balance of a given currency (like the IOTA cryptocurrency), then a site represents a transaction moving funds from a sender address to a receiver address and two sites are conflicting if they try to move the same funds to two different receivers i.e., if both executing transactions results in a negative balance for the sender. The details of this example are outside the scope of this paper, but we may use this terminology in the remaining of the In each site, the first number is its score and the second is its cumulative weight. The two tips (with dashed border) are not confirmed yet and have cumulative weight of 1.

paper. In this case, signing a transaction means performing the PoW to create the site that will be included in the Tangle, and broadcasting a transaction means sending it to the other nodes so that they can include it to their local Tangle. At each round, each node may sign one or more transactions. For each transaction, the node selects two parents. The signed transaction becomes a site in the DAG. Then, the node broadcast the site to all the other node.

DAG extension

Definition 1. Let G be a DAG and A a set of sites. If each site of A has its parents in A or in the tips of G, then we say that A is an extension of G and G ∪ A denotes the DAG composed by the union of sites from G and A. We also say that A extends G.

One can observe that if A extends G, then the tips of G forms a cut of G ∪ A. Definition 2. Let A be a set of sites extending a DAG G. We say A completely extends G (or A is a complete extension of G) if all the tips of G ∪ A are in A. In other word, the sites of A confirm all the tips of G.

The local DAG of a node may contain conflicting sites. If so, the DAG is said to be forked (or conflicting). A conflict-free sub-DAG is a sub-DAG that contains no conflicting sites.

Weight and Hashing Power When a transaction is signed and become a site in the DAG, a small proof of work (PoW) is done. The difficulty of this PoW is called the weight of the site. Initially, this PoW has been added to the protocol to prevent a node from spamming a huge number of transactions. In order to issue a site of weight w a processing power (or hashing power ) proportional to w needs to be consumed.

With the PoW, spamming requires a large amount of processing power, which increases its cost and reduces its utility. It was shown [START_REF] Popov | The tangle. white paper[END_REF] that site should have bounded weight and for simplicity, one can assume that the weight of each site is 1.

Then, this notion is also used to compute the cumulative weight of a site, which is the amount of work that has been done to deploy this site and all sites that confirm it. Similarly, the score of a site is the sum of all weight of sites confirmed by it i.e., the amount of work that has been done to generate the sub-DAG confirmed by it, see Figure 1 for an illustration.

Tip Selection Algorithm

When signing a transaction s, a node u has to select two parents i.e., two previous sites in its own version of the DAG.According to the protocol, this is done by executing an algorithm called the tip selection algorithm (TSA). The protocol says that the choice of the parents must be done among the sites that have not been confirmed yet i.e., among tips. Also, the two selected parents must not confirm, either directly or indirectly, conflicting sites.

We denote by T the TSA, which can depend on the implementation of the protocol. For simplicity, we assume all the nodes use the same algorithm T . As pointed by previous work [START_REF] Popov | The tangle. white paper[END_REF], the TSA is a fundamental factor of the security and the stability of the Tangle.

For our analysis, we assume T depends only on the topology of the DAG, on the weight of each site in the DAG and on a random source. It is said to be statefull if it also depends on previous output of the TSA by this node, otherwise we say it is stateless.

The output of T depends on the current version of the DAG and on a random source (that is assumed distinct for two different nodes). The random source is used to prevent different nodes that has the same view from selecting the same parents when adding a site to the DAG at the same time. However, this is not deterministic and it is possible that two distinct nodes issue two sites with the same parents.

Local Main DAG The local DAG of a node u may contain conflicting sites. For consistency, a node u can keep track of a conflict-free sub-DAG main u (G) that it considers to be its local main DAG. If there are two conflicting sites a and ā in the DAG G, the local main DAG contains at most one of them.

The main DAG of a node is used as a reference for its own view of the world, for instance to calculate the balance associated with each address. Of course this view may change over the time. When new transactions are issued, a node can change its main DAG, updating its view accordingly (exactly like in the bitcoin protocol, when a fork is resolved due to new blocks being mined). When changing its main DAG, a local node may discard a sub-DAG in favor of another sub-DAG. In this case, several sites may be discarded. This is something we want to avoid or at least ensure that the probability for a site to be discarded tends quickly to zero with time.

The tip selection algorithm decides what are the tips to confirm when adding a new site. Implicitly, this means that the TSA decides what sub-DAG is the main DAG. In more detail, the main DAG of a node at round r is the sub-DAG confirmed by the two sites output by the TSA. Thus, a node can run the TSA just to know what is its main DAG and even if no site has to be included to the DAG.

One can observe that, to reach consensus, the TSA should ensure that the main DAG of all the nodes contain a common prefix of increasing size that represents the transactions everyone agree on.

Adversary Model

Among the nodes, some are honest i.e., they follow the protocol, and some are byzantine and behave arbitrarily. For simplicity, we can assume that only one node is byzantine and we call this node the adversary. The adversary is connected to the network and receive all the transactions like any other honest node. He can behave according to the protocol but he can also create (and sign) transactions without broadcasting them, called hidden transaction (or hidden sites). To make the results stronger, we can assume that the adversary cannot sign a message using another node's identity. Here, even if we use the term "signing", sites may not have identity attached to them so it is actually not relevant to decide whether or not adversary can steal the identity of honest nodes.

When an honest node issues a new site, the two sites output by T must be two tips, at least in the local DAG of the node. Thus, one parent p 1 cannot confirm indirectly the other p 2 , because in this case the node is aware of p 2 having a child and is not a tip. Also, a node cannot select the same site as parent for two different site, thus the number of honest children cannot exceed the number of nodes in the network. This implies the following property.

Property 1. In a DAG constructed by n honest nodes using a TSA, a site cannot have one parent that confirms the other one. Moreover, the number of children of each site is bounded by n.

The first property should be preserved by an adversary as it is easy for the honest nodes to check and discard a site that does not verify it. However the adversary can issue multiple sites that directly confirm the same site and the honest nodes have no way to know which sites are honest.

Assiduous Honest Majority Assumption

The cumulative weight and the score can be used by a node to select its main DAG. However, even if it is true that a heavy sub-DAG is harder to generate than a light one, there is no relation yet in the protocol between the weight of sites and the hashing power capacity of honest nodes.

We define the assiduous honest majority assumption as the fact that the hashing power of honest nodes is constantly used to generate sites and that it is strictly greater than the hashing power of the adversary. In fact, without this assumption, it is not relevant to look at the hashing power of the honest nodes if they do not constantly use it to generates new sites.

Thus, under this assumption, the cumulative weight of the honest DAG grows according to the hashing power of the honest nodes, and the probability that an adversary generates more sites than the honest nodes in a given period of time tends to 0 as the duration of the period tends to infinity. Conversely, without this assumption, an adversary may be able to generates more sites than the honest nodes, even with less available hashing power.

Average Number of Tips and Confirmation Time

In this section we study the average number of tips depending on the rate of arrival of new sites. In this section, like in previous analysis [START_REF] Popov | The tangle. white paper[END_REF], we assume that tip selection algorithm is the simple uniform random tip selection that select two tips uniformly at random.

We denote by N (t) the number of tips at time t and λ(t) the number of sites issued at time t. Like previously, we assume λ(t) follows a Poisson distribution of parameter λ. Each new site confirms two tips and we denote by C(t) the number of sites confirmed at time t, among those confirmed sites, C tips (t) represents the number of tips i.e., the number of sites that are confirmed for the first time at time t. Due to the latency, if h > 1, a site can be selected again by the TSA. We have:

N (t) = N (t -1) + λ(t) -C tips (t)
We say we are in state N ≤ 1 if there are N tips at time t. Then, the number of tips at each round is a Markov chain (N (t)) t≥0 with an infinite number of states [1, ∞). To find the probability of transition between two states (given in Lemma 2) we first calculate the probability of transition when the number of new site is known. Proof. If k new sites are issued, then there are up to 2k sites that are confirmed. This can be seen as a "balls into bins" problem [START_REF] Mitzenmacher | Probability and computing: Randomized algorithms and probabilistic analysis[END_REF] with 2k balls thrown into N bins, and the goal is to see how many bins are not empty i.e. how many unique sites are confirmed. First, there are N 2k possible outcome for this experience so the probability of a particular configuration is 1 N 2k . The number of ways we can obtains exactly C = N -N + k non empty bins, or confirmed transaction (so that there are exactly N tips afterward) is the number of ways we can partition a set of 2k elements into C parts times the number of ways we can select C bins to receive those C parts.

P N k →N = N ! N 2k (N -k)! 2k N -N + k
The first number is called the Stirling number of the second kind and is denoted by

2k N -N +k . The second number is N ! (N -k)! .
Then, the probability of transition is a direct consequence of the previous lemma Lemma 2. The probability of transition from N to N is

P N →N = N k=|N -N | P(Λ = k)P N k →N = N k=|N -N | N !λ k e -λ N 2k (N -k)!k! 2k N -N + k
Proof. We just have to observe that the probability of transition from N to N is null if the number of new sites is smaller than N -N (because each new site can decrease the number of tips by at most one), smaller than N -N (because each site can increase the number of tips by at most one), or greater than N (because each new site is a tip). Lemma 3. The Markov chain (N (t)) t≥0 has a positive stationary distribution π.

Proof. First, it is clear that (N (t)) t≥0 is aperiodic and irreducible because for any state N > 0, resp. N > 1, there is a non-null probability to move to state N + 1, resp. to state N -1. Since it is irreducible, we only have to find one state that is positive recurrent (i.e., that the expectation of the hitting time is finite) to prove that there is a unique positive stationary state.

For that we can observe that the probability to transition from state N to N > N tends to 0 when N tends to infinity. Indeed, for a fixed k, we even have that the probability to decrease the number of tips by k tends to 1:

P N k →N -k = N ! N 2k (N -2k)! = (1 - 1 N )(1 - 2 N ) . . . (1 - 2k + 1 N ) (1)
lim

N →∞ P N k →N -k = 1 (2)
So that for any ε > 0 there exists a k ε such that P(Λ ≥ k ε ) < ε/2 and from (2) an N ε such that ∀k < k ε , P

N k →N -k
< ε 2kε so we obtain:

A Nε = N >N ≥Nε P N →N = P(N (i + 1) > N (i + 1) | N (i) ≥ N ) (3)
< P(Λ ≥ k ε ) + k<kε P N k →N -k (4)
< ε

So that the probability A N to "move away" from states [1, N ] tends to 0 when N tends to infinity.

In fact, it is sufficient to observe that there is a number N 1/2 such that the probability to "move away" from states [1, N 1/2 ] is strictly smaller than 1/2. Indeed, this is a sufficient condition to have a state in [1, N 1/2 ] that is positive recurrent (one can see this by looking at a simple random walk in one dimension with a mirror at 0 and a probability p < 1/2 to move away from 0 by one and (1 -p) to move closer to 0 by 1). From the irreducibility of (N (t)) t≥0 , all the states are positive recurrent and the Markov chain admit a unique stationary distribution π.

The stationary distribution π verifies the formula π N = i≥1 π i P i→N , which we can use to approximate it with arbitrary precision.

When the stationary distribution is known, the average number of tips can be calculated N avg = i>0 iπ i , and with it the average confirmation time Conf of a tip is simply given by the fact that, at each round, a proportion λ/N avg of tips are confirmed in average. So Conf = N avg /λ rounds are expected before a given tip is confirmed. The value of Conf depending on λ is shown in Figure 3.

With this, we show that Conf converges toward a constant when λ tends to infinity. In fact, for a large λ, the average confirmation time is approximately 1.26, equivalently, the average number of tips N avg is 1.26λ. For smaller values of λ, intuitively the time before first confirmation diverges to infinity and N avg converges to 1.

A Necessary Condition for the Security of the Tangle

A simple attack in any distributed ledger technology is the double spending attack. The adversary signs and broadcast a transaction to transfer some funds to a seller to buy a good or a service, and when the seller gives the good (he consider that the transaction is finalized), the adversary broadcast a transaction that conflicts the first one and broadcast other new transactions in order to discard the first transaction. When the first transaction is discarded, the seller is not paid anymore and the funds can be reused by the adversary.

The original motivation behind our attack of the Tangle is as follow: after the initial transaction to the seller, the adversary generates the same sites as the honest nodes, forming a sub-DAG with the same topology as the honest sub-DAG (but including the conflicting transaction). Having no way to tell the difference between the honest sub-DAG and the adversarial sub-DAG, the latter will be selected by the honest nodes at some point. This approach may not work with latency in the network, because the sub-DAG of the adversary is always shorter than the honest sub-DAG, which is potentially detected by the honest nodes. To counter this, the adversary can generate a sub-DAG that has not exactly the same topology, but that has the best possible topology for the tip selection algorithm. The adversary can then use all its available hashing power to generate this conflicting sub-DAG that will at some point be selected by the honest nodes. For this attack we use the fact that a TSA selects two tips that are likely to be selected by the same algorithm thereafter. For simplicity we captured this with a stronger property: the existence of a maximal deterministic TSA.

Definition 3 (maximal deterministic tip selection algorithm).

A given TSA T has a maximal deterministic TSA T det if T det is a deterministic TSA and for any DAG G, there exists N G ∈ N such that for all n ∈ N the following property holds:

Let A det be the extension of G obtained with N G + n executions of T det . Let A be an arbitrary extension of G generated with T of size at most n, conflicting with A det , and let G = G∪A∪A det . We have:

P(T (G ) ∈ A det ) ≥ 1/2
Intuitively this means that executing the maximal deterministic TSA generates an extension this is more likely to be selected by the honest nodes, provided that it contains more sites than the other extensions. When the assiduous honest majority assumption is not verified, the adversary can use this maximal deterministic TSA at his advantage.

Theorem 1. Without the assiduous honest majority assumption, and if the TSA has a maximal deterministic tip selection, the adversary can discard one of its transaction that has an arbitrary cumulative weight.

Proof. Without the assiduous honest majority assumption, we assume that the adversary can generate strictly more sites than the honest nodes. Let W be an arbitrary weight. One can see W has the necessary cumulative weigh a given site should have in order to be considered final. Let G 0 be the common local main DAG of all node at a given round r 0 . At this round our adversary can generate two conflicting sites confirming the same pair of parents. One site a is sent to the honest nodes and the other ā is kept hidden.

The adversary can use T det the maximal deterministic TSA of T to generate continuously (using all its hashing power) sites extending G ∪ {ā}. While doing so, the honest nodes extend G ∪ {a} using the standard TSA T . After r W rounds, it can broadcast all the generated sites to the honest nodes. The adversary can choose r W so that (i) the probability that it has generated N G more sites than the honest nodes is sufficiently high, and (ii) transaction a has the target cumulative weight W .

After receiving the adversarial extension, by definition 3, honest nodes will extend the adversarial sub-DAG with probability greater than 1/2. In expectation, half of the honest nodes start to consider the adversarial sub-DAG as their main DAG, thus the honest nodes will naturally converge until they all chose the adversarial sub-DAG as their main DAG, which discard the transaction a.

If the bandwidth of each channel is limited, then the adversary can start broadcasting the sites of its conflicting sub-DAG at round r W , at a rate two times greater than the honest nodes. This avoids congestion, and at round r W + r W /2 all the adversarial sub-DAG is successfully received by the honest nodes. Due to this additional latency, the number of sites in the honest sub-DAG might still be greater than the number of sites in the adversarial sub-DAG, so the adversary continues to generate and to broadcast sites extending its conflicting sub-DAG and at round at most 2r W , the adversarial extension of G received by the honest nodes has a higher number of sites than the honest extension.

So the same property is true while avoiding the congestion problem. Now that we have our main theorem, we show that any TSA defined in previous work (especially in the Tangle white paper [START_REF] Popov | The tangle. white paper[END_REF]) has a corresponding maximal deterministic TSA. To do so we The rectangle site conflicts with all site in A det so that when executing the TSA on G ∪ A ∪ A det , tips either from A or from A det are selected. The strategy to construct A det can be either to increase the number of children of G tips or to increase their weight ; both ways are presented here. can see that to increase the probability for the adversarial sub-DAG to be selected, the extension of a DAG G obtained by the maximal deterministic TSA should either increase the weight or the number of direct children of the tips of G as shown in Figure 4. We now prove that the three TSA presented in the Tangle white paper [START_REF] Popov | The tangle. white paper[END_REF], (i) the random tip selection, (ii) the MCMC algorithm and (iii) the Logarithmic MCMC algorithm, all have a maximal deterministic TSA, which implies that the assiduous honest majority assumption is necessary when using them (recall that we do not study the sufficiency of this assumption).

The Uniform Random Tip Selection Algorithm

The uniform random tip selection algorithm is the simplest to implement and the easiest to attack. Since it chooses the two tips uniformly at random, an attacker just have to generates more tips than the honest nodes in order to increase the probability to have one of its tips selected. Lemma 4. The Random TSA has a maximal deterministic TSA.

Proof. For a given DAG G the maximal deterministic T det always chooses as parents one of the l tips of G. So that, after n + l newly added sites A det , the tips of G ∪ A det are exactly A det and no other extension of G of size n can produce more than n + l tips so that the probability that the random TSA select a tip from A det is at least 1/2. Corollary 1. Without the assiduous honest majority assumption, the Tangle with the Random TSA is susceptible to double-spending attack.

The MCMC Algorithm

The MCMC algorithm is more complex than the random TSA. It starts by initially putting a fixed number of walker on the local DAG. Each walker performs a random walk towards the tips of the DAG with a probabilistic transition function that depends on the cumulative weight of the site it is located to and its children. In more details, a walker at a site v has a probability p v,u to move to a child u with

p v,u = exp(-α(w(v) -w(u))) c∈Cv exp(-α(w(v) -w(c)) (6)
Where the set C v is the children of v, and α > 0 is a parameter of the algorithm.

The question to answer in order to find the maximal deterministic TSA of MCMC algorithm is: what is the best way to extend a site v to maximize the probability that the MCMC walker chooses our sites instead of another site. The following Lemma shows that the number of children is an important factor. This number depends on the value α.

Lemma 5. Suppose a MCMC walker is at a site v. There exists a constant C α such that if v has C α children of weight n, then, when extending v with an arbitrary set of sites H of size n, the probability that the walker move to H is at most 1/2.

Proof. When extending v with n sites, one can choose the number h of direct children, and then how the other sites extends those children. There are several ways to extends those children which changes their weights w 1 , w 2 , . . . , w h . The probability p H for a MCMC walker to move to H is calculated in the following way:

S H = h 1 exp(-α(W -w i )) S H = C α exp(-α(W -n)) S = S H + S H p H = S H /S
The greater the weight the greater the probability p H . Adding more children, might reduce their weights (since H contains only n sites). For a given number of children h, there are several way to extends those children, but we can arrange them so that each weight is at least n -h + 1 by forming a chain of length n -h and by connecting the children to the chain with a perfect binary tree. The height l i of a children i gives it more weight. So that we have w i = n -h l i . A property of a perfect binary tree is that h 1 2 l i = 1. We will show there is a constant C α such that for any h and any l 1 , . . . , l h , with h 1 2 l i = 1, we have

S H ≥ S H C α exp(-α(W -n)) ≥ h i=1 exp(-α(W -w i )) C α ≥ h i=1 exp(-α(h -l i )) (7)
Surprisingly, one can observe that our inequality does not depend on n, so that the same is true when we arrange the sites when extending a site v in order to increase the probability for the walker to select our sites. Let f h : (l 1 , . . . , l h ) → e -αh h 1 exp(αl i ). So the goal is to find an upper bound for the function f h that depends only on α.

The function f h is convex (as a sum of convex functions), so the maximum is on the boundary of the domain, which is either (l 1 , . . . , l h ) = (1, 2, . . . , h) or (l 1 , . . . , l h ) = ( log(h) , . . . , log(h) , log(h) , . . . , log(h) )

For simplicity, let assume that h is a power of two so that the second case is just ∀i, l i = log(h).

In the first case we have

f h (1, . . . , h) = exp(-αh) exp(α(h + 1)) -exp(-α) exp(-α) -1 = exp(α) -exp(-α(h + 1)) exp(-α) -1
which tends to 0 when h tends infinity, so it admits a maximum C 1 α

In the second case, we have

f h (1, . . . , h) = exp(-αh)h exp(α log(h))
which again tends to 0 when h tends infinity, so it admits a maximum C 2 α By choosing C α = max(C 1 α , C 2 α ) we have the inequality [START_REF] Kuśmierz | The first glance at the simulation of the tangle: discrete model[END_REF] for any value of h.

Lemma 6. The MCMC tip selection has a maximal deterministic TSA.

Proof. Let G be a conflict-free DAG with tips G tips .

Let T be the number of tips times the number C α defined in Lemma 5. The T first executions of T det select a site from G tips (as both parents) until all site from G tips as exactly C α children.

The next executions of T det selects two arbitrary tips (different if possible). After T executions, only one tip remains and the newly added sites form a chain.

Let N G = 2T . N G is a constant that depends only on α and on G. After N G + n added sites, each site in G tips has a C α children with weight at least n. Thus, by Lemma 5, a MCMC walker located at a site v ∈ G tips moves to our extension with probability at least 1/2. Since this is true for all sites in G tips and G tips is a cut. All MCMC walker will end up in A det with probability at least 1/2. One can argue that this is not optimal and we could have improved the construction of the extension to reduce the number of sites, but we are mostly interested here in the existence of such a construction. Indeed, in practice, the probability for a walker to move to our extension would be higher as the honest sub-DAG A is not arbitrary, but generated with the TSA. Our analyze shows that even in the worst configuration, the adversary can still generate an extension with a good probability of being selected.

Corollary 2. Without the assiduous honest majority assumption, the Tangle with the MCMC TSA is susceptible to double-spending attack.

The Logarithmic MCMC Algorithm

In the Tangle white paper, it is suggested that the MCMC probabilistic transition function can be defined with the function h → h -α = exp(-α ln(h)). In more details, a walker at a site v has a probability p v,u to move to a child u with

p v,u = (w(v) -w(u)) -α c∈Cv (w(v) -w(c)) -α (8)
Where the set C v is the children of v, and α > 0 is a parameter of the algorithm. The IOTA implementation currently uses this function with α = 3. With this transition function, the number of children is more important than their weight.

Lemma 7. The logarithmic MCMC tip selection has a maximal deterministic TSA.

Proof. Let G be a conflict-free DAG with tips G tips .

Let T be the number of tips.

T det always selects two sites from G tips in a round-robin manner. After kT executions (k ∈ N), each site from G tips has 2k children.

Let n be the number of sites generated with T det and A an arbitrary extension of G. Let v ∈ G tips and C v be the number of children of v that are in A and C det = 2n/T be the number of children of v generated by T det .

With w(v) the weight of v, we have that w(v) ≤ 2n and C det ≤ w(v) -n ≤ n Let p be the probability that a walker located at v chooses a site generated by T det . We have

p ≥ C det (w(v) -1) -α C v (w(v) -n) -α ≥ C det (2n) -α C v (C det ) -α = C det C v T α = 2n T 1-α C v
With T a constant and C v bounded, we have that p tends to infinity as n tends to infinity. This is true for each site of G tips , so after a given number of generated site N G , the probability that a LMCMC walker located at any site of G tips moves to a site generated by by T det is greater than 1/2.

Corollary 3. Without the assiduous honest majority assumption, the Tangle with the Logarithmic MCMC TSA is susceptible to double-spending attack.

Discussion

Sparse Network

The case of sparse networks intuitively gives more power to the adversary as it is more difficult for the honest nodes to coordinate. Let assume that the communication graph is arbitrary. It can be a geometric graph (eg. a grid), a small-world graph with a near constant diameter, or anything else. In order for the protocol to work properly, we assume that for each link between two nodes, there is enough bandwidth to send all the sites generated by the honest nodes and that the usage of each link is a constant fraction of its available capacity. For simplicity we can assume that no conflict occurs from the honest nodes so that, without adversary, the local DAG of each node is its main DAG. Due to multi-hop communications, the local DAGs of the honest nodes may differ but only with respect to the sites generated during the last D rounds, where D is the diameter of the communication graph. Take an arbitrary node u in this graph. We connect our adversary to node u, so that every sites received by u at round r is received by our adversary at round r + 1.

In this arbitrary network, the attack is exactly the same, except for the number of rounds r W that the adversary waits before revealing the conflicting sub-DAG. Indeed, r W should be greater to take into account the propagation time and ensure that that the first adversarial site a has a cumulative weight greater than W in all the honest nodes (typically we should wait D more rounds compared to the previous case). As in the previous case, the adversary can broadcast its conflicting sub-DAG while ensuring not to induce congestion, for instance, at a rate two times greater than the honest. The topology of the network does not change the fact that after at most 2r W + D rounds, all the honest nodes have a greater probability to choose the adversarial sub-DAG for their next transactions.

Trusted Nodes

The use of trusted nodes is what makes IOTA currently safe against this kind of attacks. Indeed, a coordinator node regularly add a site to the DAG, confirming an entire conflict-free sub-DAG. Trusted sites act like milestones and any site confirmed by a trusted site is considered irreversible. However if the trusted node is compromised or is offline for too long, the other nodes are left alone.

The current implementation of IOTA uses a trusted node called the coordinator and plans to either remove it, or replace it by a set of distributed trusted nodes.

One can observe that the crypto-currency Byteball [START_REF] Churyumov | Byteball: a decentralized system for transfer of value[END_REF] uses a special kind of trusted nodes called witnesses. They are also used to resolve conflicts in the DAG.

An important question could be: is the use of trusted node necessary to secure a distributed ledger based on a DAG?

Avoiding Forks

In the current protocol, conflicting sites cannot be confirmed by the same site. Part of the community already mentioned that this can cause some problem if the latency is high (i.e., if diameter of the communication graph is large). Indeed, by sending two conflicting sites to two different honest nodes, half of the network can start confirming one transaction, and the other half the other transaction. When a node finally receives the two conflicting transactions ( assuming that the convergence of is instantaneous) a single site will be assumed correct and become part of all the main DAG of all the honest node. However all the sites confirming the other conflicting site are discarded, resulting in wasted hashing power and increase in confirmation time for those site as a reattachment must occur. The cost for the adversary is small and constant and the wasted hashing power depends on the maximum latency between any two nodes.

On way to avoid this is to include conflicting sites in the DAG (like in the Byteball protocol [START_REF] Churyumov | Byteball: a decentralized system for transfer of value[END_REF] for instance), by issuing a site confirming the two conflicting sites and containing the information of what site is considered valid and what site is considered invalid. This special site called decider site would be the only site allowed to confirm directly or indirectly two conflicting sites. This has the advantage that all the site confirming the invalid one remains valid and do not need to be reattached. However, the same thing can happen if the adversary send two conflicting decider sites to two end of the network. But again, a decider site could be used to resolve any kind of conflict, including this one. Indeed, this may seems like a circular problem that potentially never ends, but every time a decider is issued, a conflict is resolved and the same conflict could have happened even without the decider site. So having decider site should no change the stability of the Tangle and only help avoiding reattaching sites.

Conclusion

We presented a model to analyze the Tangle and we used it to study the average confirmation time and the average number of unconfirmed transaction over the time.

Then, we defined the notion of assiduous honest majority that captures the fact that the honest nodes have more hashing power than the adversarial nodes and that all this hashing power is constantly used to create transactions. We proved that for any tip selection algorithm that has a maximal deterministic tip selection (which is the case for all currently known TSA), the assiduous honest majority assumption is necessary to prevent a double-spending attack on the Tangle.

Our analyze shows that honest nodes cannot stay at rest, and should be continuously signing transactions (even empty ones) to increase the weight of their local main sub-DAG. If not, their available hashing power cannot be used to measure the security of the protocol, like we see for the Bitcoin protocol. Indeed, having a huge number of honest nodes with a very large amount of hashing power cannot prevent an adversary from attacking the Tangle if the honest nodes are not using this hashing power. This conclusion may seem intuitive, but the fact that it is true for all tip selection algorithms (that have a deterministic maximal TSA) is something new that have not been proved before.

Figure 1 :

 1 Figure 1: An example of a Tangle where each site has a weight of 1. In each site, the first number is its score and the second is its cumulative weight. The two tips (with dashed border) are not confirmed yet and have cumulative weight of 1.

Lemma 1 .

 1 If the number of tips is N and k new sites are issued, then the probability P N k →N of having N tips in the next round is:

  where a b denotes the Stirling number of the second kind S(a, b).

Figure 2 :

 2 Figure 2: Stationary distribution of the number of tips, for different values of λ. For each value of λ, one can see that the number of tips is really well centered around the average

Figure 3 :

 3 Figure 3: Expected number of round before the first confirmation, depending on the arrival rate of transaction. We see that it tends to 1.26 with λ. Recall that Conf = N avg /λ where N avg refers to the average number of tips in the stationary state.

Figure 4 :

 4 Figure4: A and A det are two possible extensions of G. The rectangle site conflicts with all site in A det so that when executing the TSA on G ∪ A ∪ A det , tips either from A or from A det are selected. The strategy to construct A det can be either to increase the number of children of G tips or to increase their weight ; both ways are presented here.