N Strelkov

A Chavent

A Timopheev

R C Sousa

I L Prejbeanu

L D Buda-Prejbeanu

B Dieny

Impact of Joule heating on the stability phase diagrams of perpendicular magnetic tunnel junctions

Keywords: Spintronics, MRAM, STT switching, phase diagrams, Joule heating, LLG linearization

Measured switching voltage-field (V-H) diagrams of perpendicular magnetic tunnel junctions (pMTJ) exhibit unexpected behavior at high voltages associated with significant heating of the storage layer. The boundaries deviate from the critical lines corresponding to the coercive field, which contrasts with the theoretically predicted behavior of a standard macrospin-based model. Combining recent experimental studies of the temperature dependence of spin polarization and perpendicular magnetic anisotropy (PMA) we are proposing a modified model. Our approach takes into account the Joule heating during the writing pulse, which reduces the spin polarization and the anisotropy, thereby reducing the spin torque efficiency and the coercive field during the switching. The numerical macrospin simulations based on this model are in agreement with our experimental measurements and consistent with the results derived from the linearization of Landau-Lifshitz-Gilbert equation.

INTRODUCTION

A large number of studies related to spin transfer torque (STT) switching [1,2] in magnetic tunnel junctions (MTJ) were conducted over the past few years in support of the development of spin transfer torque magnetic random access memory (STT-MRAM). This memory uses as storage elements MgO-based magnetic tunnel junctions with perpendicular (i.e. out-of-plane) magnetization (pMTJs) [3,4]. The common stacks involve MgO/FeCoB interfaces allowing to reach the technological requirements in terms of magnetoresistance variation around 200% [5,6] and thermal stability factor at room temperature around 60 [7,8]. The main properties of the memory MTJ cell can be extracted from their writing voltage-field (V-H) diagram [9]. Experimental observations and theoretical modelling of such diagrams in collinear [6,[10][11][12] and non-collinear geometries [13] give an almost complete description of STT-induced switching processes in MTJ. However, at high voltages (typically above 0.4 V-0.8 V depending on the junction resistance), most experimental diagrams show deviation from the theoretically predicted behavior. A plausible origin of this effect is the heating due to Joule dissipation around the tunnel barrier and the resulting variation of the storage layer magnetic properties. It is well-known that the coercive field of magnetic layers decreases with temperature [14,15]. STT reduction due to the heating was also noticed earlier in experiments on current-induced domain wall motion [16,17] and it was observed both in spin valves [18] and MTJs [19]. Bandiera et al. proposed a concept of thermally assisted STT-MRAM in which the heating due to the STT write current is also used to induce an anisotropy reorientation from out-ofplane to in-plane [20]. Moreover, recent experiments have shown that the MgO/CoFeB interface with thin MgO tunnel barrier have a much lower (more than two orders of magnitude) thermal conductivity than what could be expected from the bulk value [21]. In double barrier MTJs in which the storage layer is sandwiched between two MgO barriers, this further enhances the Joule heating of the storage layer.

In this paper, we present a series of experimental studies and their theoretical interpretation using the macrospin model that takes into account the heating effects described above. These results point out the strong influence of the Joule heating on switching abilities of MTJs at high voltages. Considering the importance of these heating effects, it is mandatory to include these heating effects in MRAM design tools.

EXPERIMENTS

We performed a number of experiments using different pillars based on pMTJ stack with different compositions and sizes. Typical switching diagrams are summarized in FIG. 1 for various kind of samples varying the composition of the stack and the nanofabrication procedure (see TAB I). The wafers in FIG. 1(a,b,d) were grown by dc and rf magnetron sputtering on a thermally oxidized Si substrate. The sample from FIG. 1(c) was grown on pre-patterned conducting non-magnetic pillar without post-deposition etching [START_REF] Nguyen | IEEE International Electron Devices Meeting (IEDM)[END_REF].

The storage layer from FIG. 1(a,b,d) consists of two ferromagnetic parts separated by a thin W film and sandwiched between two MgO layers to increase the interfacial anisotropy: the main MgO barrier which provides the TMR and a thinner MgO layer with lower resistancearea (RA) product. Next to the MgO layer, there is a V o l t a g e ( V ) thin W capping layer is inserted to absorb the B away from the MgO interfaces upon annealing [6,12,[START_REF] Chatterjee | [END_REF]. The magnetization of the bottom reference layer is pinned by a synthetic antiferromagnet (SAF).

The measurement process to derive the V-H diagrams is the same as described in Ref. [11]. A magnetic field is applied out-of-plane, i.e. along the easy axis of magnetization of the storage layer and slowly varied (typically 50 kOe/s) after each write-read sequence. A 100 ns writing voltage pulse of defined amplitude is applied at each magnetic field point. The resistance is measured under low constant (dc) current after applying the writing pulse. The resistance-field (R-H) loop is repeated several times, and the final resistance value is defined by averaging all these loops.

In the FIG. 1, one can see V-H diagrams with various shapes. Three main stability areas are identified: the high-resistance area with antiparallel (red, AP) configuration of storage and reference layer, the low-resistance area with parallel (blue, P) configuration and bistable (green, P/AP) area where both mentioned configurations are possible. The center of the bistable area is shifted with respect to the zero field. This shift corresponds to an uncompensated stray field H stray from the SAF.

The shape of the diagram in FIG. 1(a) is predictable by the theory developed in Ref. [11] -there are straight vertical boundaries corresponding to the coercive field. But the three other diagrams show deviations from the theoretically expected shape. This unpredictable behavior can be due to Joule heating effects inside the storage layer. A number of recent experiments supports this interpretation. For example, it was observed earlier that the STT write current flowing through pMTJs could induce a decrease of both PMA and coercivity due to Joule heating around the tunnel barrier [20]. In FIG. 1(c), it can be clearly seen that the critical lines lean back towards the H stray value at high voltages and completely coincide with H stray when the absolute value of the applied voltage exceeds ∼0.6 V. This could mean that the PMA and spin polarization in storage layer are reduced by Joule heating to the point where the critical lines then correspond to the applied field value which compensates the SAF stray field H stray . Further comparing FIG. 1(a) and FIG. 1(b), a stronger thermal influence can be noted in the smaller junction. This can be explained by considering that smaller samples should be more sensitive to the temperature variation, as the stability factor of MTJ in macrospin approximation ∆ = K eff • Ω/(k B T ) depends on the volume Ω of the storage layer, K eff being the effective perpendicular anisotropy constant. For a given voltage, one might expect that the longer the writing pulse, the higher the temperature reached in the storage layer and therefore stronger the impact of the heating. To check this effect on our samples, the measurement procedure was changed. A constant (dc) voltage was applied to the sample and the applied field varied step by step from -1400 Oe to 1400 Oe and vice versa. At each field step, the current was measured during 500 µs and the resulting resistance calculated as R = V bias /I. In FIG. 2, the normalized experimental V-H diagram measured using this procedure is shown. As can be seen, there are no more vertical straight boundaries corresponding to field equal to coercive field. Since here the current is steadily flowing through the sample, the magnetic tunnel junction cannot cool down between each resistance measurement, so that a significant decrease of the anisotropy due to Joule heating takes place as the bias voltage is gradually increased.

To estimate the temperature variation and its voltage dependence in MTJs, we performed an experiment by the method described in Ref. [24] based on the temperature dependence of RKKY coupling, using the sample similar to that from FIG. 1(d) but with larger diameter D=170 nm. We measured the V-H diagram under dc voltage and extracted the voltage dependence of spin-flop field H SF (white solid line in FIG. 3(a)), at which both SAF layers become magnetized in parallel configuration versus applied voltage. This dependence is parabolic, and its expression can be derived from FIG. 3(b). We then measured the variation of H SF field versus controlled temperature which is found to be almost linear FIG. 3(c). By combining the data from FIG. 3(b) and FIG. 3(c), the temperature elevation ∆T due to current flow versus applied voltage can be derived. A quadratic variation is obtained given by ∆T = 308 • V 2 and shown in FIG. 3(d

) (T in K, V in volt).
On the V-H diagrams in FIG. 1 one can note some asymmetry between the critical lines for positive and negative voltage polarities. For example, the absolute value of switching voltage V sw at H=H stray (black vertical line) corresponding to zero applied field, is different for positive (V + sw ) and negative (V - sw ) voltages. One of the possible reasons for such asymmetry might be due to the different resistances for P-AP and AP-P transitions. If the heating power at fixed voltage reads as V 2 /R, then lower resistance causes more heating and, as a result, a greater decrease of PMA and spin polarization. In FIG. 1(a), the TMR amplitude is 124%, V + sw =0.37 V (AP-P transition with high resistance and lower temperature), V - sw =0.33 V (P-AP transition with low resistance and higher temperature). The switching voltage difference for this sample is 12%. In FIG. 1(b), the TMR amplitude is 139%, V + sw =0.37 V, V - sw =0.3 V and difference is 23%. For the sample with TMR amplitude 66% in FIG. 1(d),V + sw =0.54 V, V - sw =0.52 V and the difference is 4%. So, we can state that the sample with larger TMR amplitude have more pronounced asymmetry in its V- V o l t a g e ( V ) H diagram due to larger Joule heating in P state (low resistance) than in AP state (high resistance).

H S F = 1 . 0 -0 . 0 0 1 9 2 ⋅ ∆T d ) ∆T = 3 0 8 ⋅ V 2 ∆T ( K ) V o l t a g e ( V )
In FIG. 1(c), the sample has a very low TMR amplitude, so the mechanism of asymmetric heating described above is negligible. But one can still see a small difference in switching voltage: V + sw =0.32 V, V - sw =0.34 V. This polarity dependence in the shape of V-H diagram might originate from the heating asymmetry induced by tunneling current [25]. Most of Joule energy is released after the electrons tunneling in the receiving magnetic electrode adjacent to the MgO barrier. Therefore, considering that the thermal conductivity of thin MgO barrier is much lower than that of bulk MgO [21], the storage layer reaches up to 10% higher temperature if the current polarity is such that the storage layer receives the tunneling electrons. This effect becomes more significant in the case of MgO-based capping layer since the storage layer appears to be thermally insulated between two MgO layers with low thermal conductivity.

Another source of asymmetry in the V-H diagram can result from the voltage dependence of the torque amplitude associated with elastic tunneling in MTJ [26]. Indeed, the parallel component of the torque T varies as T ∼ a 1 V + a 2 V 2 . The transport parameter a 2 is known to vanish for MTJs with ferromagnetic electrodes having large spin splitting or being half-metallic. Such MTJs are expected to exhibit high TMR amplitude. In such MTJs, if heating effects are not taken into account, the switching voltages for P-AP and AP-P transitions are therefore close to opposite (V + sw =V - sw ). On the other hand, for low TMR MTJs, the parallel torque exhibits a nonsymmetric quadratic voltage behavior so that switching voltages do not coincide (V + sw = V - sw ). To get deeper insight on the above mentioned features of the experimental V-H diagrams, we developed a model accounting for the heating effects in the storage layer. The standard phenomenological macrospin-based model was modified by including a dependence of the storage layer temperature on voltage and a resulting dependence of magnetization, spin polarization and anisotropy constant on the applied voltage pulse.

MODEL

The first adjustment of the model concerns the change in anisotropy energy while changing the bias voltage. There are several possible mechanisms affecting the uniaxial anisotropy constant. A first one is a direct voltage induced variation of anisotropy due to interfacial charge modulation [27][28][29][30], but this effect is relatively weak compared to the influence of the temperature variation associated with Joule heating. A second mechanism -is a temperature induced variation of the anisotropy due to Joule heating, which is the most significant for our samples. Concerning other thermoelectric effects, as it was shown by Flipse et al [31] and later by Avery and Zink [32] the Peltier coefficient is of order Π ∼10 -3 V and thus the Peltier power at R ∼10 3 Ω and

V ∼1 V is Q Π ∼ Π • V /R ∼10 -6 W, while the Joule power Q J = V 2 /R ∼10 -3
W is larger by three orders of magnitude. Seebeck effect is also negligible as its coefficient S is expressed through Π by Thomson-Onsager relation Π = S • T .

We first estimated the temperature dependence on the bias voltage in MTJ starting from the first law of thermodynamics and describing the MTJ with a simple one dimensional model according to the following differential equation:

V 2 R -Q(T -T 0 ) = C dT dt , (1)
where V 2 /R is a Joule power, R -resistance, Q is the heat transfer coefficient between the tunnel barrier and the conductive bottom and top electrodes assumed to be at a fixed temperature T 0 , T 0 is the room temperature and C is the MTJ heat capacity. Equation (1) gives a solution for the temperature T as a function of pulse duration t p :

T = T 0 + V 2 R Q 1 -e -(Q/C)tp . ( 2
)
We assume that the temperature of the storage layer follows this square dependence on voltage pulse amplitude (2) and all other parameters are kept constant for a given sample. Then, a second step in the model development consists in assuming that the temperature dependence of the uniaxial anisotropy constant is follows the Calen-Calen law [33]:

K(T ) = K 0 M s (T ) M s0 ξ , (3)
where K 0 is the uniaxial anisotropy constant at zero temperature, M s0 is the spontaneous magnetization of the storage layer at zero temperature. From previous studies, the value of the exponent ξ was experimentally estimated to lie in the range between 2 and 3 [20,[34][35][36][37]. The temperature dependence of the storage layer magnetization is supposed to follow a Bloch law:

M s (T ) = M s0 1 -(T /T c ) 3/2 , (4)
where T c is the Curie temperature of the storage layer.

The exponent "3/2" is slightly different for thin deposited films, and its experimentally determined value equals to 1.73 [34]. Taking into account ( 3) and ( 4) and assuming that pulse durations are long enough so that the temperature varies with the voltage as T = T 0 +k V V 2 according to (2), the voltage dependences of the magnetization and anisotropy constant reads:

M s (V ) = M s0 1 - T 0 + k V V 2 T c 1.73 , K(V ) = K 0 1 - T 0 + k V V 2 T c 1.73 ξ , (5)
where k V is a parameter dependent on the writing pulse duration, the thermal capacity of the storage layer, electrical resistance and the thermal resistance of the stack around the storage layer. The parameter k V depends also on the polarity of the applied voltage as mentioned earlier [25]. Actually in the following of our model, the heating asymmetry was not taken into account and k V was kept independent on current polarity. We want to underline that in case of pulsed V-H diagram, the observed back switching anomalies are not only caused by the voltage induced anisotropy reduction but also by the reduction of spin polarization. As the anisotropy amplitude is restored at the end of the voltage pulse, the fact that the storage layer had time to switch or not during the pulse has a great importance. If STT efficiency becomes too weak due to the heating, the storage layer may not be able to switch any more. This is the reason why we can cross the critical lines twice on V-H diagram at a given applied field (e.g. FIG. 1(c) increase voltage at 150 Oe from 0 V to 0.6 V).

The phenomenological macrospin model used here is the same as described in Ref. [11,13] but it takes into account the voltage dependence of several parameters. Landau-Lifshitz-Gilbert (LLG) equation for the magne-tization dynamics is written as:

dm dt = -γ (m × µ 0 H eff (V )) + α m × dm dt -γa (V )V [m × (m × p)] + γa ⊥ (V )V 2 (m × p) , (6)
where m is the unit vector of the magnetization of the storage layer, α is the Gilbert damping factor, γ is the gyromagnetic ratio of free electrons, p is the unit vector along the spin current polarization (magnetization of the reference layer), a (V ) and a ⊥ (V ) are phenomenological transport parameters depending on the applied voltage. The effective field H eff (V ) in our model depends on applied voltage and derives from the Gibbs free energy E(V ) density functional:

E(V ) = -K(V )(u K • m) 2 -µ 0 M s (V )H ext (V ) • m + 1 2 µ 0 M 2 s (V ) N x m 2 x + N y m 2 y + N z m 2 z , (7)
where H ext is the external applied field, u K is the unit vector along the easy axis (here Oz), N x,y,z are diagonal terms of the demagnetizing tensor, M s (V ) and K(V ) are defined from (5). The first term in (7) represents the anisotropy energy E K (V ), so the corresponding effective anisotropy field reads:

H K (V ) = 2K(V ) µ 0 M s (V ) = H K0 1 - T 0 + k V V 2 T c 1.73 ξ , (8)
where H K0 = 2K 0 /(µ 0 M s0 ). The dependence of normalized effective anisotropy field (8) on applied voltage for different values of exponent ξ is presented in FIG. 4. In the inset, the same dependence is plotted versus the temperature reached by the storage layer during the voltage pulse T = T 0 + k V V 2 . This dependence is almost linear as it was shown in the experiment [34,38]. Assuming the ground state of storage layer to be out-of-plane, the field-like torque term ∼ a ⊥ (V ) in ( 6) is neglected here since the Slonczewski term ∼ a (V ) is known to play a dominant role in the switching process. The phenomenological transport parameter a (V ) is proportional to the spin polarization which follows the magnetization variation given by (5) (see solid black line in FIG. 4) as reported in Ref. [19]. Consequently, due to the Joule heating, the STT prefactor has the following dependence on voltage:

a (V ) = a 0 1 - T 0 + k V V 2 T c 1.73 , (9)
where a 0 represents the transport parameter value at zero temperature.

Having set the model as described above, the numerical switching V-H diagrams are computed following the experimental procedure. The magnetization of the storage FIG. 4. Normalized effective anisotropy field versus applied voltage according to (8) for different values of the exponent ξ using the parameters: T0=300 K, kV =1900 K/V 2 , Tc=1200 K. The curve for ξ=2 (solid black line) repeats the normalized magnetization Ms(V )/Ms0. In the inset the same function versus resulting temperature

T = T0 + kV V 2 .
layer is relaxed under the given applied field. After that a writing pulse is applied for 100 ns. The total integration time is 1 µs. The amplitude of applied field is changed from its minimum to maximum value and vice versa to reproduce a hysteresis loop. The initial room temperature value was set at T 0 =300 K. The thermal fluctuations are described by random field added to the total effective field as a white noise field with the following properties:

H TR (t) =0, H TR (t)H TR (t ) = 2αk B T γµ 0 M s Ω δ(t -t ), ( 10
)
where k B is the Boltzmann constant and Ω is the sample volume. Each hysteresis loop at fixed voltage was repeated 10 times resulting an averaged V-H loop similar to the experimental diagram. The numerical V-H diagrams computed for a disk of 100 nm diameter and 1.6 nm thick are presented in FIG. 5 for several values of the exponent ξ. The shape of the diagram evolves progressively with the value of the exponent ξ. Two regimes are identified. At low values of ξ (around 2), the numerical diagram is very similar to the experimental diagram plotted in FIG. 1(d). In contrast, when ξ approaches the upper limit of 3, a strong similarity with the experimental diagram FIG. 1(b,c) is observed. In this second regime, two specific critical lines are identified: a first one, STT driven, at low voltages and a second one, heat driven, at high voltages. The parameter ξ accounts for the strength of the dependence of the material parameters on the voltage. The impact of its value on the shape of the V-H diagrams is quite drastic, confirming the important role of Joule heating on switching abilities of MTJs. To better understand the underlying parameters behind these critical lines and eventually find how to control them, we developed an analytical model based on the stability analysis of the linearized LLG equation. Numerical stability V-H diagrams for different values of exponent ξ at room temperature T0=300 K, kV =1900 K/V 2 , Tc=1200 K, Ms0=10 6 A/m, K0=778319 J/m 3 , α=0.01, a ⊥0 =0, a 0 =8•10 -3 T/V, p=(0, 0, -1), Nx=Ny=0.025, Nz=0.95 using 10 averaged loops for each voltage with resolution 64x64 points. Critical lines (white solid curves) calculated from equations ( 18) and superposed with the numerical V-H diagrams. Vertical dashed lines correspond to the coercive field calculated from the solution of the equations ( 17) at zero applied voltage. For ξ=3.0 the intersection points of two solutions are circled.

ANALYTICAL MODEL

The critical lines of the numerically obtained V-H diagrams presented in FIG. 4 (white lines) can be extracted analytically at zero temperature using the generic model of a nonlinear auto-oscillator proposed by Slavin and Tiberkevich [39]. Following the Holstein-Primakoff transformation the components of the unitary magnetization vector m in (6) are replaced by the canonical variables c and c * as proposed in Ref. [40]:

c = m x -ım y 2(1 + m z ) . (11)
It is convenient to express the modified LLG equation ( 6) in complex variables with four terms: precession, damping, field-like and damping-like:

dc dt = dc dt prec + dc dt damp + dc dt FL + dc dt DL . (12)
Assuming that the polarization vector p=(0,0,p z ) and the applied field H ext =(0,0,H z ) are both perfectly out-ofplane, the full expression of each term is explicitly given below:

dc dt prec = -ıc[ω H + (ω M -ω A )(2|c| 2 -1)], dc dt damp = -αc(1 -|c| 2 )× ω H + (ω M -ω A )(2|c| 2 -1) , dc dt FL =ı αω + ω ⊥ p z c, dc dt DL = -(ω -αω ⊥ )(1 -|c| 2 )p z c, (13)
where

ω H = γ 1 + α 2 H z , ω A = γ 1 + α 2 2K u µ 0 M s , ω = γ 1 + α 2 a V µ 0 , ω ⊥ = γ 1 + α 2 a ⊥ V 2 µ 0 , ω M = γ 1 + α 2 (N z -N x ) M s . (14)
Since the complex variable c has an amplitude p and a phase φ associated by the relation c(t) = p(t)e φ(t) , one can then derive two equations:

dp dt = -2α ω H + ω /α -ω ⊥ p z + (ω M -ω A )(2p -1) (p -1)p, dφ dt = -ω H -(αω + ω ⊥ )p z + (ω M -ω A )(2p -1) . (15)
Looking for the stationary solutions dp/dt=0, two static solutions are found: i) p 0 = 0 meaning magnetization "up" m z = +1 and ii) p 0 = 1 meaning magnetization "down" m z = -1.

By linearizing the equation dp/dt in (15) around the solution p 0 and keeping only the first order terms, one can derive an equation for the power deviation δp:

d(δp) dt p0 =                        -2α ω H + (ω /α -ω ⊥ )p z +(ω M -ω A ) δp, p 0 = 1 +2α ω H + (ω /α -ω ⊥ )p z -(ω M -ω A ) δp, p 0 = 0 (16)
If the pre-factor in front of δp becomes positive, the deviation δp will diverge and the static equilibrium becomes unstable. From this condition, by setting the pre-factors to zero, one can derive the boundaries of stability under static field and/or voltage. Two critical lines for the applied field without any applied voltage are given by:

ω H + (ω M -ω A ) = 0, ω H -(ω M -ω A ) = 0. ( 17
)
The two critical lines for dc field and voltage are:

ω H + (ω /α -ω ⊥ )p z + (ω M -ω A ) = 0, ω H + (ω /α -ω ⊥ )p z -(ω M -ω A ) = 0. (18)
Substituting M s , K u and a in ( 14) by the corresponding functions of the bias voltage given by expressions ( 5) and ( 9), the solution of equations ( 17) and (18) gives critical the curves H(V ) which coincide well with the boundaries of the numerical diagrams. In FIG. 5, the critical lines are superposed to the numerical V-H diagrams. Vertical dashed lines correspond to the coercive field calculated from the solution of equations ( 17) at zero applied voltage, which correspond to a finite pulse duration. If the voltage pulse is infinite (dc) the solution (17) reproduces the boundaries like in FIG. 2(a) -coercive field decreases with applied bias voltage. The small difference between the numerical and analytical coercive field is due to the fluctuating term (10). It increases the probability of magnetization switching before the field reaches its coercive value and yields this 10% deviation. The critical lines under dc field and dc voltage are rather complex and partially recover the numerical boundaries at low ξ and low voltage pulses. At ξ value close to 3, the critical lines at high voltages are closer to that in the dc diagram from FIG. 2(a). The analytical dc curves reproduce well the two intersection points through which the numerical critical lines pass.

In the case where ξ equals 2 or 3 and the magnetization follows the Bloch law (4), one can find the solution for intersection points of the critical lines defined by equations (18) (see yellow circles in FIG. 5 for ξ=3.0). Considering that the magnetization depends on voltage due to Joule heating and introducing the function f B (V ) in ( 5) defined by f B (V ) = M s (V )/M s0 , neglecting the field-like torque term ω ⊥ =0 and assuming the polarization is "up" p z =1, one gets the following equations from ( 14) and ( 18):

H z + a 0 V µ 0 α f B (V )+ M s0 (N z -N x )f B (V ) -H K0 f ξ-1 B (V ) = 0, H z + a 0 V µ 0 α f B (V )- M s0 (N z -N x )f B (V ) -H K0 f ξ-1 B (V ) = 0, (19)
where H K0 = 2K 0 /(µ 0 M s0 ). Now it is clear that to satisfy these equations (19) at the same value of applied field H z , the last term in brackets should be equal to zero. If the exponent ξ=2, then the solution is trivial: f B (V )=0. The intersection points correspond to heating the storage layer up to the Curie temperature:

V c0 = ± (T c -T 0 )/k V .
If the exponent ξ=3, then besides the trivial solution, there is another one which reads:

V c1 = ± T c k V (1 -η) 2/3 - T 0 T c , (20)
where η = M s0 (N z -N x )/H K0 . As one can see, the ratio V c0 /V c1 depends on the factor η and Curie temperature T c . The thermally induced slope that passes through these intersection points also depends on these parameters but does not dependent on the heating coefficient k V :

dV dH ξ=3 = V c1 -V c0 H c1 = - 1 η µ 0 α a 0 1 - 1 -T 0 /T c (1 -η) 2/3 -T 0 /T c . (21)
Therefore, this slope should not depend on the voltage pulse duration as it does not depend on k V . On the other hand, if the exponent ξ=2, the thermally induced slope can be calculated as an inverted derivative of the function H z (V ) extracted from (19) taken at the point V =V c0 :

dV dH ξ=2 = dH dV V =Vc0 -1 = - 1 3
a 0 µ 0 α 1 - T 0 T c + M s0 (N z -N x ) 1 -η η k V T c 1 -T 0 /T c -1 . (22)
The slope slightly increases with decreasing writing pulse duration. In FIG. 6, the critical lines extracted from the experimental diagrams for different values of the writing pulse duration are presented. The switching voltage is almost constant at pulse durations more than 100 ns and increases sharply for shorter pulses as it was observed earlier [13]. The thermally induced boundaries remain almost constant at pulses longer than 30 ns and moves upward for shorter pulses with slight increase of the slope. The macrospin model used in this study gives an overestimate value for the coercive field which can be found from equations (17) at zero applied voltage:

H c = ± M s0 (N z -N x )× 1 - T 0 T c x 1 η 1 - T 0 T c x ξ-2 -1 . (23)
For the FeCoB-based storage layer with a diameter of 100 nm, theoretical H c in macrospin approximation should be of the order of ∼1 T (10 4 Oe). But in the experiment, its value is usually ten times smaller ∼0.1 T (10 3 Oe) due to the non-uniform switching mode (not macrospin). Similarly, the STT switching is not macrospin for the used MTJ diameters. Therefore, to fit the experimental diagrams with this model, one has to use an effective value for the uniaxial anisotropy constant to account for the fact that the barrier height for switching is different from that of macrospin uniform switching. The result of the fitting is presented in FIG. 7. The first sample in FIG. 7 As the TMR amplitude is larger than for the sample of FIG. 7(a), the asymmetry of STT boundaries is more pronounced; the STT boundary corresponding to the P-AP transition is shifted towards the zero voltage, because more heat is produced in the low resistance state (P) in the storage layer which reduces the switching voltage. To take into account the heating asymmetry in our model,

SUMMARY

To conclude, we measured the stability V-H diagrams in a variety of samples of different composition and diameters in which significant Joule heating take place during STT writing. We explained the high voltage anomaly of the V-H diagrams as due to heating effects in the storage layer. Two mechanisms of the asymmetry of the critical lines for positive and negative voltages were suggested: heating asymmetry induced by the asymmetric inelastic relaxation of tunneling electrons and heating asymmetry due to the difference of the resistance in P and AP states. The second mechanism makes the largest contribution to the heating asymmetry versus bias voltage when the TMR is large. We also developed a modified macrospin-based model by including a dependence of the material parameters on the applied voltage, such as uniaxial anisotropy constant (K u ), spontaneous magnetization (M s ) and STT prefactor a . The model was used to numerically simulate the V-H diagrams. Good agreement with the experimental results was obtained. Furthermore, we derived explicit analytical expressions for the critical lines that exactly follow the contours of the numerical diagram. The obtained expressions can be applied to fit also the experimental diagrams, but some effective values for K u and M s should be used, as the switching of the real samples occurs in a non-uniform regime (for example, by domain wall propagation) which reduces the coercive field.

The presented model can be used in pMTJ-based MRAM circuits design to take into account such unusual behavior of the MTJ properties at high bias voltages. Moreover, the predictions presented in this study can be extrapolated to similar multilayer stacks used for other applications.

FIG. 1 .

 1 FIG. 1. Experimental stability V-H diagrams for the different samples based on pMTJ stacks. The color scale corresponds to the resistance of pMTJ. Field shift with respect to the H=0 Oe corresponds to an uncompensated stray field Hstray from the SAF (a) 12 R-H loops, Hstray=192 Oe, TMR=124% (b) 12 R-H loops, Hstray=352 Oe, TMR=139% (c) 12 R-H loops, Hstray=72 Oe, TMR=9% (d) 10 R-H loops, Hstray=379 Oe, TMR=66%. See sample details in TAB. I.

FIG. 2 .

 2 FIG. 2. Normalized experimental stability V-H diagram measured under constantly applied (dc) voltage. Sample is the same as in FIG. 1(b) (see details in TAB. I).

2 N

 2 o r m a l i z e d s p i n -f l o p f i e l d V o l t a g e s q u a r e d ( V 2 ) c ) N o r m a l i z e d s p i n -f l o p f i e l d ∆T ( K )

FIG. 3 .

 3 FIG. 3. (a) Experimental stability V-H diagram measured under dc voltage for the sample similar to that from FIG. 1(d) with D=170 nm and extracted spin-flop field-voltage dependence (white solid line). (b) Linear fit of normalized spin-flop field-voltage squared dependence. (c) Linear fit of normalized spin-flop field-temperature variation dependence. (d) Temperature elevation due to current flow versus voltage derived by combining the fits from FIG. 3(b) and FIG. 3(c).

  p e r a t u r e ( K )

  FIG. 5.Numerical stability V-H diagrams for different values of exponent ξ at room temperature T0=300 K, kV =1900 K/V 2 , Tc=1200 K, Ms0=10 6 A/m, K0=778319 J/m 3 , α=0.01, a ⊥0 =0, a 0 =8•10 -3 T/V, p=(0, 0, -1), Nx=Ny=0.025, Nz=0.95 using 10 averaged loops for each voltage with resolution 64x64 points. Critical lines (white solid curves) calculated from equations (18) and superposed with the numerical V-H diagrams. Vertical dashed lines correspond to the coercive field calculated from the solution of the equations (17) at zero applied voltage. For ξ=3.0 the intersection points of two solutions are circled.

FIG. 6 .

 6 FIG.6. Critical lines extracted from the experimental diagrams of the sample from FIG.1(c) for different lengths of applied voltage pulse. Thermally induced boundaries remain almost constant at pulses longer than 30 ns and move upwards for shorter pulses with slight increase of the slope. See sample details in TAB. I.

  (a) has a low TMR value and the critical lines are almost symmetric. FIG. 7(b) shows a diagram of the sample from the same wafer as the sample in FIG. 1(b) measured under finite voltage pulse.

FIG. 7 .

 7 FIG. 7. Experimental diagrams fit with the analytical solution. (a) Sample from FIG. 1(c) with low TMR, ξ=2.16, Ms0=900 kA/m, K0=570 kJ/m 3 , kV =1200 K/V 2 , Tc=700 K, a 0 =1.4 mT/V, Nx=Ny=0.01, Nz=0.98. Other parameters the same as in FIG. 5. (b) Sample from the same wafer as the sample in FIG. 1(b) with high TMR=142%, ξ=2.3, Ms0=900 kA/m, K0=580 kJ/m 3 , kV =340 K/V 2 , Tc=700 K, a 0 =2 mT/V. Dashed line corresponds to the critical curve at kV =500 K/V 2 . Other parameters the same as in FIG. 5.

TABLE I .

 I Stack composition of the used samples

	# Figure	Sample description
	S1 FIG. 1(a) SAF/MgO/FeCoB(1.3nm)/W(0.3nm)/
		FeCoB(0.5nm)/MgO/Capping,
		D=110 nm, RA=10 Ω • µm 2
	S2 FIG. 1(b), SAF/MgO/FeCoB(1.3nm)/W(0.3nm)/
	S2b FIG. 2,	FeCoB(0.5nm)/MgO/Capping,
	FIG. 7(b)	D=80 nm, RA=10 Ω • µm 2
	S3 FIG. 1(c),	SAF/MgO/FeCoB(1.5nm)/
	FIG. 6	W(2nm)/Capping,
	FIG. 7(a)	D=300 nm, RA=9 Ω • µm 2
	S4 FIG. 1(d) SAF/MgO/FeCoB(0.9nm)/W(0.1nm)/
	S4b FIG. 3	Co(0.5nm)/W(0.1nm)/FeCoB(0.8nm)/
		Cappng, D=70 nm, RA=9 Ω • µm 2

Acknowledgements: This work was supported by ERC Advanced Grant MAGICAL n o 669204.


Matthieu Lequesne
email: matthieu.lequesne@inria.fr

Jean-Pierre Tillich
email: jean-pierre.tillich@inria.fr

Attack on the EDON-K Key Encapsulation Mechanism

The key encapsulation mechanism EDON-K was proposed in response to the call for post-quantum cryptography standardization issued by the National Institute of Standards and Technologies (NIST). This scheme is inspired by the McEliece scheme but uses another family of codes defined over F 2 128 instead of F2 and is not based on the Hamming metric. It allows significantly shorter public keys than the McEliece scheme.

In this paper, we give a polynomial time algorithm that recovers the encapsulated secret. This attack makes the scheme insecure for the intended use. We obtain this result by observing that recovering the error in the McEliece scheme corresponding to EDON-K can be viewed as a decoding problem for the rankmetric. We show that the code used in EDON-K is in fact a super-code of a Low Rank Parity Check (LRPC) code of very small rank (1 or 2). A suitable parity-check matrix for the supercode of such low rank can be easily derived from for the public key. We then use this parity-check matrix in a decoding algorithm that was devised for LRPC codes to recover the error. Finally we explain how we decapsulate the secret once we have found the error.

I. INTRODUCTION

The syndrome decoding problem is a fundamental problem in complexity theory, since the original paper of Berlekamp, McEliece and van Tilborg [START_REF] Berlekamp | On the inherent intractability of certain coding problems[END_REF] proving its NP-completeness for the Hamming distance. The same year, McEliece proposed a public-key cryptosystem based on this problem [START_REF] Mceliece | A Public-Key System Based on Algebraic Coding Theory[END_REF] and instantiated it with binary Goppa codes. This scheme was for a long time considered inferior to RSA due to its large key size. However, this situation has changed drastically when it became apparent in [START_REF] Shor | Algorithms for quantum computation: Discrete logarithms and factoring[END_REF] that RSA and actually all the other public-key cryptosystems used in practice could be attacked in polynomial time by a quantum computer. There are now small prototypes of such computers that lead to think that they will become a reality in the future and in 2016, the National Institute of Standards and Technology (NIST) announced a call for standardization of cryptosystems that would be safe against an adversary equiped with a quantum computer. Four families of cryptosystems are often mentioned as potential candidates: cryptosystems based on error correcting codes, lattices, hash functions and multivariate quadratic equations [START_REF]Post-Quantum Cryptography[END_REF]. All of these are based on mathematical problems that are expected to remain hard even in the presence of a quantum computer.

The key encapsulation mechanism EDON-K [START_REF] Gligoroski | first round submission to the NIST post-quantum cryptography call[END_REF] was proposed by Gligoroski and Gjøsteen in response to the call issued by the NIST. This scheme is inspired by the McEliece scheme but uses another family of codes defined over F 2 128 instead of F 2 . This choice leads to very short keys for a codebased scheme. The metric used for the decoding is not properly defined and the security relies on an ad-hoc problem named finite field vector subset ratio problem supposedly hard on average.

In this paper, we show that the metric used for EDON-K is in fact equivalent to the well-known rank metric. This metric was first introduced in 1951 as "arithmetic distance" between matrices over a field F q [START_REF] Hua | A theorem on matrices over a sfield and its applications[END_REF]. The notion of rank distance and rank codes over matrices was defined in 1978 by Delsarte [START_REF] Delsarte | Bilinear forms over a finite field, with applications to coding theory[END_REF]. He introduced a code family, named maximum rank distance (MRD) codes, that attains the analogue of the MDS (maximum distance separable) bound for the rank metric. Gabidulin suggests in [START_REF] Mukhamedovich | Theory of codes with maximum rank distance[END_REF] to consider a subfamily of such codes that are linear over an extension field F q m . It provides a vectorial representation of these codes and allows to represent them in a much more compact way. This is the main reason why the rank metric based McEliece schemes achieve significantly smaller key sizes. Moreover this vectorial representation allows to view the known families of MRD codes as rank metric analogues of Reed-Solomon codes and to obtain an efficient decoding algorithm for them [START_REF] Mukhamedovich | Theory of codes with maximum rank distance[END_REF]. There are also rank metric analogues for other families of codes. For instance, the Low Rank Parity-Check (LRPC) codes introduced in [GMRZ13] can be considered as analogues of Low Density Parity-Check (LDPC) codes. Just like their binary cousins, they enjoy an efficient decoding algorithm that is based on a low rank parity-check matrix of such a code.

Here, we prove that the code used in EDON-K is a actually a super-code of an LRPC code of rank 2. What is more, this LRPC code is itself a subspace of codimension 1 of another LRPC code of rank 1. It turns out that parity-check matrices of rank 2 for the first super-code and rank 1 for the second one can easily be derived from the public key. In both cases, this allows us to decode the ciphertext without the secret key. This gives a way to recover the encapsulated secret and breaks completely the EDON-K system.

The paper is organized as follows. First, we recall some basic definitions and properties of rank metric and LRPC codes in Section II. In Section III we present the scheme of EDON-K. Then we explain the general idea of our attack in section IV. In Section V, we detail how we reconstruct a parity-check matrix of the code and in Section VI how we decode the ciphertext. In Section VII, we explain how we derive the encapsulated secret from the error. Finally in Section VIII we discuss the cost of this attack and its consequences.

II. RANK METRIC CODES

A. Notation

In the following document, q denotes a power of a prime number. In the case of EDON-K, we will have q = 2. F q denotes the finite field with q elements and, for any positive integer m, F q m denotes the finite field with q m elements. We will sometimes view F q m as an m-dimensional vector space over F q .

We use bold lowercase and capital letters to denote vectors and matrices respectively.

We denote x 1 , . . . , x k K the K-vector space generated by the elements {x 1 , . . . x k }.

B. Definitions

Definition 1 (Rank metric over F n q m ). Let x = (x 1 , . . . , x n ) ∈ F n q m and (β 1 , . . . , β m ) be a basis of F q m viewed as an mdimensional vector space over F q . Each coordinate x j ∈ F q m is associated to a vector of F m q in this basis:

x j = m i=1 m i,j β i .
The m × n matrix associated to x is given by M(x) := (m i,j ) 1≤i≤m,1≤j≤n .

The rank weight wt(x) of x is defined as :

wt(x) := Rank M(x).
The associated distance d(x, y) between elements x and y of F n q m is defined by d(x, y) := wt(xy). Definition 2 (Support of a word). Let x = (x 1 , . . . , x n ) ∈ F n q m . The support of x, denoted Supp(x), is the F q -subspace of F q m generated by the coordinates of x:

Supp(x) := x 1 , . . . , x n Fq .
We have dim(Supp(x)) = wt(x).

Definition 3 (F q m -linear code). An F q m -linear code C of dimension k and length n is a subspace of dimension k of F n q m . C can be represented in two equivalent ways: by a generator matrix G ∈ F k×n q m such that C = {xG | x ∈ F k q m } and by a parity-check matrix H ∈ F (n-k)×n q m such that C = {x ∈ F n q m | Hx = 0 n-k }.
The decoding problem in the rank metric can be described as follows.

Problem 1 (Decoding problem for the rank metric). Let C be an F q m -linear code of dimension k and length n. Given y = c + e where c ∈ C and e ∈ F n q m is of rank weight ≤ r find c and e.

C. LRPC codes

Definition 4 (LRPC code). A Low Rank Parity Check (LRPC) code of rank d, length n and dimension k over F q m is a code that admits a parity-check matrix H = (h i,j ) ∈ F (n-k)×n q m such that the vector space of F q m generated by its coefficients h i,j has dimension at most d.

LRPC codes can be viewed as analogues of LDPC codes for the rank metric. In particular, they enjoy an efficient decoding algorithm based on their low rank parity-check matrix. Roughly speaking, Algorithm 1 of [START_REF] Gaborit | Low rank parity check codes and their application to cryptography[END_REF] decodes up to d errors when rd ≤ n -k in polynomial time (see [GMRZ13, Theorem 1]). It uses in a crucial way the notion of the linear span of a product of subspaces of F q m Definition 5. Let U and V be two F q subspaces of F q m . We denote by U • V the linear span of the product of U and V :

U • V := uv : u ∈ U, v ∈ V Fq . III. THE EDON-K KEM EDON-K [GG17]
is a key encapsulation mechanism proposed by Gligoroski and Gjøsteen for the NIST post-quantum cryptography call. Here we describe the key generation, encapsulation and decapsulation, omitting some details that are not relevant for the attack. We refer to [START_REF] Gligoroski | first round submission to the NIST post-quantum cryptography call[END_REF] for the full description.

A. Parameters and notations

The parameters for EDON-K are given in Table I. In this paper we often refer to the parameters of edonk128ref, the reference version proposed for 128 security-bits. The scheme makes use of a hash function H (•) corresponding to standard SHA2 functions (SHA-256 or SHA-384 depending on the parameters). We will denote

H i (•) := H(. . . H(•)) i times
Given a binary matrix P = (p i,j ) and two non-zero elements a = b of F 2 m , P a,b = (p i,j ) denotes the matrix of the same size with coefficients in F 2 m where pi,j = a if p i,j = 0 and pi,j = b if p i,j = 1.

In particular, if P is orthogonal:

P c,d = (P a,b ) -1 (1)
where c := a a 2 +b 2 and d := b a 2 +b 2 . For two vectors (or matrices) x and y, we will denote x||y their concatenation.

B. Key generation

Given the security level and the appropriate parameters.

• a, b $ ← F 2 m non-zero elememts such that a = b.

• P $ ← F N ×N 2 an orthogonal matrix. • H $ ← F R×N 2 such that H = [H T ||H B ]
where H B is an R × R orthogonal matrix and H T has columns of even Hamming weight.

• c := a a 2 +b 2 , d := b a 2 +b 2 . • g $ ← F ν 2 m . • V g := Support(g). • G $ ← V K×N g such that GH = 0 K×R .
(2)

• G pub := GP c,d .

• Return (PublicKey := G pub , SecretKey := (a, b, P, H)).

C. Encapsulation

Given the PublicKey and the public parameters.

• m

$ ← F K 2 m . • ẽ ∈ F L 2 m generated as follows: -(ẽ 0 , ẽ1 ) $ ← F 2 m ; -for 1 ≤ i ≤ L 2 -1, (ẽ 2i , ẽ2i+1 ) = H (ẽ 2i-2 ||ẽ 2i-1 ). • V e := Support(ẽ). • e $ ← V N e . • c := mG pub + e. • (s 0 , s 1 ) := H (ẽ L-2 ||ẽ L-1 ). • SharedSecret := H (s 0 ||s 1 ||H (c)). • h := H (s 1 ||s o ||H (c)). • Ciphertext := (c, h). • Return (Ciphertext, SharedSecret).

D. Decapsulation

Given Ciphertext, SecretKey and the public parameters.

• Recover e by decoding the c using the private matrix H := HP a,b . • Deduce V e the vector space spaned by the coefficients of the vector e.

• For all (λ, ν) ∈ V e × V e , for 1 ≤ i ≤ L 2 -1: -(s 0 , s 1 ) := H i (λ||µ||H (c)); -if H (s 1 ||s 0 ||c) = h: Return SharedSecret := H (s 0 ||s 1 ||H (c)).

IV. OUTLINE OF THE ATTACK ON EDON-K Our attack is based on three observations

• The ciphertext is a vector c such that c = mG pub + e.

(3)

This error e is of low rank, since its rank is at most • This code C pub generated by G pub is a subcode of an LRPC code, namely the code C with parity-check matrix H := HP a,b . This code is indeed an LRPC code of rank 2 since all the entries of H belong to a, b F2 . We have

C pub ⊂ C (4)
since

G pub H = GP c,d (HP a,b ) = GP c,d P a,b H = GH (from (1)) = 0 K×R (from (2)).
This equation also appears as Corollary 1 of [GG17, p.19]. We have given its proof here for the convenience of the reader. Let K = N -R be the dimension of C . • If we recover a parity-check matrix of rank 2 for C we will be able to recover mG pub and e from c. Indeed, mG pub ∈ C and we can decode C using a variation of Algorithm 1 of [START_REF] Gaborit | Low rank parity check codes and their application to cryptography[END_REF] and the knowledge of the parity-check matrix, provided wt(e) ≤ L < (N -K )/2 = R/2 is verified, which is the case for the parameters of EDON-K. Hence we will proceed in three steps: 1) constructing and solving a linear system of equations to find a parity-check matrix for the code C (detailed in Section V); 2) decoding the ciphertext using a slight variation of Algorithm 1 of [START_REF] Gaborit | Low rank parity check codes and their application to cryptography[END_REF] (see Section VI); 3) recovering the secret from the error vector (explained in Section VII).

V. RECONSTRUCTING THE PARITY-CHECK MATRIX

A. Compressed public key

In order to reduce the public key size, the designers of EDON-K chose to represent the public key in a compressed form. They took advantage of the fact that all the coefficients of G pub live in the vector space V g,c,d := cg 1 , . . . , cg ν , dg 1 , . . . , dg ν F2 of dimension 2ν. Hence, the compressed public key consists in two parts: first the basis gc,d := (cg 1 , . . . , cg ν , dg 1 , . . . , dg ν ) ∈ F 2ν 2 m of the vectorspace V g,c,d , then the entries of the matrix G pub such that each entry is represented by its coefficients in the basis gc,d . For example, if an entry x of G pub is equal to c ν i=1 γ i gi + d ν i=1 δ i gi with γ i , δ i ∈ F 2 , x will be represented by (γ 1 , . . . , γ ν , δ 1 , . . . , δ ν ) ∈ F 2ν 2 . There is another subtlety in the compression that we will not mention here.

B. Finding a basis

The attacker does not have access to the value of a and b but can deduce the value of ab -1 = cd -1 = (cg 1 )(dg 1 ) -1 from gc,d as mentioned in paragraph 7.2.2 of the documentation of EDON-K [START_REF] Gligoroski | first round submission to the NIST post-quantum cryptography call[END_REF].

Let us bring in α := ab -1 .

We notice that H" := b -1 H is also a parity-check matrix of the LRPC code C . This matrix has all its coefficients in 1, α F2 . We use this information to reconstruct such a parity-check matrix of the code C by solving a linear system, similarly to what is done in [GRS16, Section IV B]. This system is derived from the following facts:

(i) G pub H = 0 K×R ; (ii) the entries of H belong to 1, α F2 . In other words, the possible rows x = (x 1 , . . . , x N ) of H are solutions of the following system

G pub x = 0 K x i ∈ 1, α F2 for all i ∈ {1, . . . , N }. ( 5
)
This system is linear over F 2 and the solution set is an F 2 -linear subspace which can be computed in O(KmN 2 ) operations. A basis of this subspace can then be used as rows for H . The details on how to transform this system into a linear system over F 2 can be found in the long version of the article [START_REF] Lequesne | Attack on the Edon-K key encapsulation mechanism[END_REF].

VI. DECODING STEP

The previous step recovers an R × N matrix H (3) whose entries all belong to 1, α F2 . The matrices H (3) and H share the property that their rows form a basis of solutions of (5). Therefore, there exists an R × R binary invertible matrix Q such that

H (3) = QH . (6)
We use H (3) to decode and recover e from the ciphertext c. The vectors are linked by the equation

c = mG pub + e. (7)
We use here a slight variation of Algorithm 1 of [START_REF] Gaborit | Low rank parity check codes and their application to cryptography[END_REF] to decode. Algorithm 1 would consist in performing the following steps: 1) Compute s := H (3) c and then V := Support(s). Here we typically have

V = Support(e) • 1, α Fq when H (3) is a random matrix. 2) Compute V := V ∩ α -1 V .

This step typically recovers

Support(e) when V = Support(e) • 1, α Fq . 3) Once we have Support(e) we recover e = (e 1 , . . . , e N ) by solving the linear equation H (3) e = s with the additional constraints e i ∈ Support(e) for i ∈ {1, . . . , N }. In our case, due to the special structure of H which contains only a's and b's V is not equal to Support(e) • 1, α Fq . This is due to the following result.

Proposition 1. We have for every e ∈ F N 2 m :

Support(H (3) e ) ⊂ (1 + α)Support(e) + N i=1 e i F2 .
Proof. From (6), we deduce that Support(H (3) e ) = Support(QH e ) = Support(H e ).

Let s := H e . Denote the i-entry of s by s i and the entry of H in row i and column j by h ij . We have:

s i = N j=1 h ij e j = j s.t. h ij =1 e j + j s.t. h ij =α αe j = N j=1 e j + (1 + α) j s.t. h ij =α e j .
This implies the proposition. This proposition directly gives a subspace of dimension L + 1 that contains Support(e) since we deduce from it that

Support(e) ⊂ (1 + α) -1 Support(H e). (8)
A slight modication of Algorithm 1 of [START_REF] Gaborit | Low rank parity check codes and their application to cryptography[END_REF] yields therefore e:

1) compute the syndrome s := H (3) c and then V := (1 + α) -1 Support(s); 2) The space V contains Support(e), so we can recover e = (e 1 , . . . , e N ) by solving the linear equation H (3) e = s with the additional constraints e i ∈ V for i ∈ {1, . . . , N }. Note that we can also skip step 2 and directly look for s 0 and s 1 in the space V of dimension L + 1 instead of decoding exactly the value of e. In fact, this is what is specified in the decapsulation of EDON-K.

VII. RECOVERING THE SHARED SECRET

Once we have recovered the error vector e ∈ F N q m , we need to recover s 0 and s 1 to obtain the value of SharedSecret. We know that the elements of e were picked randomly in V e = Support(ẽ).We proceed just like in the decapsulation algorithm.

We generate Support(e) which is equal to V e with high probability. More exactly, the probability that Support(e) is of dimension < L is L-1 L N . For the parameters of edonk128ref this probability is 2 -37 . In such a case, the attack might fail, but the decapsulation would fail too.

Then, among the 2 L elements of V e , we need to identify a couple of consecutive elements of ẽ to deduce the secret. For all pairs of candidates (λ, µ) ∈ V e × V e , for 1 ≤ i ≤ L 2 -1 we compute (s 0 , s 1 ) := H i (λ||µ||H (c)). If H (s 1 ||s 0 ||c) = h then we have (s 0 , s 1 ) = (s 0 , s 1 ). Finally we recover the secret SharedSecret = H (s 0 ||s 1 ||c). In total this operation requires O(L2 2L ) operations, just like the decapsulation. This is the reason why the value of L needs to remain small, otherwise the decapsulation is not possible.

VIII. CONCLUDING REMARKS

A. Cost of the attack

Let us analyze the cost of the three steps of the attack mentioned in Section IV.

Step 1 and 2 are polynomial in terms of the parameters of the code. Step 1 only uses linear algebra operations and has a complexity at most O(KmN 2 ). The complexity of step 2 is given by Theorem 1 of [START_REF] Gaborit | Low rank parity check codes and their application to cryptography[END_REF] (using n = N, k = N -R, r = L and d = 2), hence is equal to L 2 (16m + N 2 ). The complexity of step 3 is O(L2 2L ). This is not polynomial in L but L is a very small parameter (4 ≤ L ≤ 8 in the proposal). Moreover this third step is the same as the decapsulation algorithm, so L needs to stay small, otherwise the decapsulation would become too costly or even impossible. So L can be considered as a constant ≤ 10 to allow a reasonable decapsulation. Hence the most costly operation appears to be step 1.

B. Without compression of the public key

Our attack takes advantage of the compressed form of the public key that allows a direct access to the value α = ab -1 . One could think that this is the origin of the attack, and decide to express the public key in its uncompressed form to fix the attack. As a consequence, the public key would be of size K × N × m bits instead of K × N × ν bits in the compressed form. In practice the public key for edonk128ref would be 16 times longer (around 288 kbits). This inflation of the key size could be avoided by sending out a random basis of V g,c,d . However, this is not enough. There is an even more direct way to proceed, without the value of α. Instead of looking for a matrix H (3) with entries liyng in 1, α F2 , we can use the following result.

Proposition 2. There exists a full rank (R -1) × N binary matrix H (4) that satisfies G pub H (4) = 0 K×(R-1) .

Proof. Let T be a binary full-rank matrix (R -1) × R matrix that has rows of even Hamming weight. For instance we can choose

T =       1 1 0 • • • 0 0 1 1 0 . . . . . . . . . . . . . . . . . . 0 • • • 0 1 1      
.

We observe now that TH has all its entries in {0, a + b}. This follows directly from the fact that if we sum an even number of elements in {a, b} we either get 0 (if the number of a's is even, and therefore also the number of b's) or a+b (if the number of a's is odd). From this, it follows immediately that H (4) := 1 a + b TH satisties the property. First, it is clear that this is a binary matrix and we also have

G pub H (4) = 1 a + b
G pub H T = 0 K×(R-1) .

Obtaining such a matrix H (4) is straightforward. We just have to use the algorithm given in Section V to recover a basis of dimension R -1 of binary vectors x satisfying G pub x = 0 K .

We then use this matrix H (4) to compute the syndrome s = H (4) c . Since H (4) c = H (4) e we directly obtain with very high probability that Support(e) = Support(H (4) c ).

This reveals the support of the error and from there we can go directly to the last step of the attack to reconstruct the shared secret.

C. Security of the scheme

Considering the attack that we described, there is a way to recover the secret of the edonk128ref scheme from a public key without the private key in polynomial time. In practice, the attack implemented with Sage on a personal computer recovers the secret in less than a minute, so the scheme is far from achieving the 128-bits security claimed in [START_REF] Gligoroski | first round submission to the NIST post-quantum cryptography call[END_REF]. Hence this scheme is insecure for the intended use. Moreover, the cost of this attack is polynomial in terms of the parameters, so there is no proper way to increase the parameters to achieve the intended security level while keeping a reasonably small key size.