Véronique Brouet

David Leboeuf

Ping-Hui Lin

Joseph Mansart

Amina Taleb-Ibrahimi

Patrick Le Fèvre

François Bertran

Anne Forget

Dorothée Colson

ARPES view of orbitally resolved quasiparticle lifetimes in iron pnictides

Keywords: numbers: 79.60.-i, 71.18.-y, 71.30.-h

I. INTRODUCTION

The discovery of iron-based superconductors in 2008 has brought to physicists both a new family of high temperature superconductors and a new type of correlated systems [1]. They have a multiband electronic structure, where the three Fe 3d t 2g orbitals cross the Fermi level and form small hole and electron pockets. This is unusual among correlated systems and it has been argued that the value of Hund's couplings between electrons in different orbitals, rather than simply the Coulomb repulsion, controls the strength of electronic correlations in these systems [2]. These correlations may change significantly for different families or different orbitals, d xy being more correlated in many cases [3,4]. The role of the different orbitals in the electronic properties is a major issue. Angle resolved photoemission (ARPES) should be an ideal tool to discuss these questions, as it has a unique ability to probe the properties of each band and/or orbital individually. There has been a huge effort to map the electronic structure in different families of iron pnictides with ARPES and compare the results with band structure calculations [5][6][7]. However, very little is known about the quasiparticles (QP) lifetimes τ of the different orbitals, which could, in principle, be extracted from the ARPES linewidth δν ∝ 1/τ [8,9]. Despite a few studies [7,10,11], there are no systematic investigation of ARPES linewidths for the different bands of one system. It is often assumed that the lines are too broad -especially in the case of Co doped BaFe 2 As 2 [12,13] -to perform such analysis. In this paper, we argue that meaningful information can be extracted from the ARPES linewidths and that they reveal an intriguing and universal difference in behavior between d xy and d xz /d yz orbitals.

Whether the metallic behavior in iron pnictides can be described within a Fermi liquid (FL) framework is an essential question. The QP lifetimes would then obey the following formula.

/τ = γ ( ω) 2 + (πk B T ) 2 (1)
However, it is noteworthy that the resistivity of BaFe 2 As 2 becomes linear near 8% Co doping [14][15][16], deviating from the canonic T 2 dependence of a FL. This has often been taken as a sign of non FL behavior, bearing analogies with cuprates [17] or organic superconductors [18]. However, the interpretation of this linear resistivity is not straightforward, due to the contribution of holes and electrons with different orbital characters that could have different characteristics. Indeed, Rullier-Albenque et al. argued that the scattering rates remain proportional to T 2 in the entire Co doping range of BaFe 2 As 2 [15]. This is precisely where ARPES could help clarifying the situation by resolving directly the contribution of each type of carrier.

To perform our study, we take two examples : LiFeAs and BaFe 2 As 2 doped with 8% Co (called Co8 in the following). We choose them because they are both metals down to relatively low temperatures and both become superconducting below 18K for LiFeAs [19] and 22K for Co8 [20]. On the other hand, LiFeAs appears slightly more correlated in DMFT calculations [3,21]. This is mainly due to the elongation of the Fe tetrahedra along z, which tends to reduce the hopping between neighboring Fe and enhance correlations, especially for the in-plane dxy orbital [3]. Moreover, the two compounds have slightly different structures, with different stackings of the FeAs slabs [1], which introduces differences in their electronic structure that are useful to discuss the intrinsic origin of the behaviors we observe. In particular, there is no disorder in LiFeAs introduced by Co doping and the surface is non-polar, suppressing some of the problems discussed in ref. [12,13]. The general electronic structure has already been studied in some details by ARPES, both for Co8 [5,[22][23][24] and LiFeAs [10,25,26]. We add here more focused information on the separation of hole and electron bands of different orbital character and their linewidths.

In Fig. 1, we summarize the main results of this paper. In panel (a), the evolution of the lifetime in LiFeAs is shown as a function of binding energy ω for the d xz /d yz (red) and d xy (blue) orbitals on hole (open symbols) and electron (closed symbols) pockets. In panel (b), the same information is given for electron pockets in Co8. We do not give results for the hole pockets in Co8, because they significantly overlap and are more difficult to resolve from each other [22,23]. In panel (c), we further show the temperature dependence of the lifetime at the Fermi level in Co8. There are two striking effects. First, the ω dependence of the linewidths in LiFeAs are very similar for the hole and electron bands of the same orbital character. This evidences that the orbital character rather than the hole/electron character or the location on Fermi Surface (FS) determines the lifetime behavior. In this paper, we detail the different steps necessary to obtain the data in Fig. 1. We first locate each band on the FS (section III). This shows that their properties (dispersion and width) can indeed be very clearly resolved from each other. We extract effective masses by comparison to band calculations. They are enhanced by a factor ∼ 2 in Co8 for all bands and of ∼ 2 for d xz /d yz and 3.5 for d xy in LiFeAs (see Table I). This confirms predictions from DMFT [3,21] that the orbitals start to differentiate in LiFeAs, d xy being more renormalized. One potential problem in determining lifetimes from ARPES is the three dimensionality (3D) of the electronic structure, which can induce extrinsic broadening [27,28]. In section IV, we show that 3D effects are present in both compounds, but they are quite different, due to the different stacking of the FeAs slabs, so that we can rule out such an effect as the origin of the difference in d xy and d xz /d yz linewidths. Finally, we return to the discussion of the lifetimes in section V. We first discuss the absolute values of the linewidths, including the role of finite resolution and impurities. We then detail how the linewidths, determined by fitting Momentum Distribution Curves (MDC) in Å-1 , are converted to lifetimes in eV by using the slope of the band dispersion. We conclude that the difference in behavior of d xy and d xz /d yz is robust against all these problems.

II. EXPERIMENTAL DETAILS

Single crystals were grown using a FeAs self-flux method and were investigated in details by transport measurements [15,29]. ARPES experiments were carried out at the CASSIOPEE beamline at the SOLEIL synchrotron, with a Scienta R4000 analyser, an angular resolution of 0.3 • (0.015 Å-1 at 34eV) and an energy resolution better than 15 meV. For these measurements, the linearly polarized light was set in the sample plane along Table I. First and second columns : renormalization and shift needed to fit the experimental dispersion with the theoretical ones (see Fig. 2 and3). Third column : energy position with respect to EF of the top or bottom of the band (for hole bands, we use the position of the fitted calculation). Fourth and fifth columns : number of carriers n per Fe in each band, after integration over kz, from theory and experiment. In this configuration, polarization selection rules [8,30] select orbitals even with respect to the plane containing this axis and the normal to the sample surface and odd with respect to the plane perpendicular to the previous one. Note that the band parity may be different from its main orbital character, because there are 2 inequivalent Fe in the unit cell [30,31]. The photon energy selects a particular k z value, which we estimate through the free electron final state approximation, already heavily used in iron pnictides [24,32,33].

k z = 2m/ 2 * (hν -W + V 0 ) -k 2 // (2)
For the inner potential V 0 , we use 11eV for LiFeAs and 14eV for Co8 (see section IV). Note that because Co8 has a body centered structure [1], two adjacent BZ are shifted in k z , which is not the case in LiFeAs. Band structure calculations were perfomed within the local density approximation (LDA), using the Wien2K package [START_REF] Blaha | WIEN2K: An Augmented Plane Wave + Local Orbitals Program for Calculating Crystal Properties[END_REF], with the experimental structures. The Co doping of 8% was treated in the virtual crystal approximation [23].

III. BAND DISPERSION

In Fig. 2, we present the ARPES measurements of the d xz /d yz and d xy parts of the electron pockets in Co8. A sketch of the distribution of orbital character on the electron pockets expected in LDA calculations is recalled in Fig. 2(b). We have explained in [30], which bands can be observed in this experimental configuration and why. Along (a1), we observe the shallow d xz /d yz band, while along (a2), only the deeper d xy band is detected. Their very different dispersion makes it easy to ensure that only one band is observed in each case. This difference in dispersion is expected from band structure calculations. In Fig. 2(c), we show the experimental dispersions extracted by MDC analysis (symbols), together with the calculated band dispersions renormalized by ∼ 2 and shifted up (see exact values in Table I). The colors correspond to the main orbital character. The dispersion is calculated exactly at the experimental k z , obtained using Eq. 2. The asymmetry of the d xy dispersion is due to the change of k z along the cut (through k // ) and is well reproduced in the calculation. To match the data, it is however further necessary to shift up the bands by 0.1eV before renormalizing them. This corresponds to a "shrinking" of the electron pockets that we described in details before [23]. A concomitant downward shift of the hole bands allows charge conservation.

In Fig. 3, we present the electronic structure of LiFeAs, measured in similar experimental conditions. Again, we can distinguish very well the d xz /d yz and d xy parts of the electron pockets in panels (a1) and (a2). Although the two bands are very clear in this geometry, it was not the case in previous studies, where the d xy band was either completely missing [26] or observed simultaneously with d xz /d yz [25]. Three different bands form the hole pockets at the zone center, as detailed in panels (a3) to (a5). The most outer hole band has d xy character and forms a large squarish pocket [blue contour in Fig. 3(a)].

The d xz and d yz bands form smaller hole pockets around the zone center, one having roughly a circular contour and one barely crossing the Fermi level. This is in good agreement with previous ARPES measurements [25,26].

A big advantage of LiFeAs is that all these hole bands are well separated and each of them can be unambiguously attributed to one orbital. Therefore, we have a unique opportunity here to study in details the contribution of d xz /d yz and d xy , both for hole and electron pockets. In Fig. 3(b), we extract the dispersions of these different bands at k z =0 (we choose the appropriate photon energy in accordance to section IV). We compare them with the calculated bands divided by a factor 2, this value being chosen to get a global view of the situation. Obviously, the hole d xz /d yz bands are shifted down compared to this calculation. On the contrary, the electron bands and the hole d xy band are not significantly shifted from calculation. This contrasts with the global shrinking discussed before in Co8. Here, holes are mainly transferred from d xz /d yz to d xy (see table I and section IV for more details). This transfer is induced by correlations and was correctly predicted by DMFT [3] and already observed by ARPES [25,26]. This puts d xy closer from half-filling, which is one of the reasons why it is predicted to get more correlated [4].

In Fig. 3(c), we fit each band individually by adjusting the shift and renormalization with parameters reported in Table 1. As was already clear in Fig. 3(b), d xy disperses more slowly than the calculation divided by 2, both for hole and electron pockets, meaning some additional renormalization is needed. We find values of 3.3 and 4 for d xy , compared to ∼2 for d xz /d yz . This shows that the renormalization value is fixed by orbital character rather than hole or electron character. In LiFeAs, d xy starts to differentiate from other orbitals by exhibiting larger correlation effects. In Fig. 3(d), we zoom on the dispersion of each band to show in more details the quality of the fit at low energies.

IV. THREE DIMENSIONAL DISPERSION

In a 3D system, the ARPES spectra are often broadened by the photoemission process [27,28] and their linewidths may not reflect the QP lifetimes anymore. The short electron escape depth λ in the photon energy range we use (20-100eV) leads to an integration over some finite k z range (δk z ∝1/λ). This can induce severe broadening of the ARPES spectra if the k z dispersion is large. In this section, we pay particular attention to this point to understand how it could affect our linewidth analysis. a function of photon energy or equivalently of k z using Eq. 2. No variation of k F is detected for d xy , while it is clear and periodic for d xz /d yz (note the different photon energy ranges). This is in good agreement with calculations (solid lines) for which d xy is strictly 2D, while d xz /d yz exhibits clear 3D dispersion. The experimental 3D dispersion is not as strong for d xz /d yz as predicted in the calculation. In fact, the corresponding spectra [Fig. 4 (c)] are good lorentzians only near k z =0 (where we performed most of our analysis) and broaden asymmetrically towards k z =1 (see black spectra for k z =6 and 7). Assuming this distortion of lineshape is due to some averaging over k z , we extracted the variation of k F with k z using the outer edge of the spectra, as shown by the red points. The two dispersions of the electron bands are then summarized in Fig. 4(d). The dispersion extracted in this way is 30% smaller than the theoretical one. From these dependences, we can extract the number of electrons in each band. Assuming the electron pockets can be divided into a squarish inner sheet of d xz /d yz symmetry and an outer sheet of d xy symmetry, we obtain 0.1 electron in each band after integrating over k z . These numbers are compared in Table I with the theoretical ones. The number of holes, for which we did not measure significant 3D dispersion at E F , are also added to this table. We find number of carriers reduced by ∼10-20% compared to calculations and compatible with charge neutrality. We reported the 3D electronic structure of Co8 in [30] and just recall the results in Fig. 4 These differences in k z variations turn out to be useful to estimate the importance of 3D effects in the measured linewidth. In LiFeAs, one could suppose at first that d xy is narrower because it is more 2D. However, d xy is nearly as 3D as d xz /d yz in Co8, and yet, it is also narrower than d xz /d yz . Moreover, the k z dispersion is nearly twice smaller for d xz /d yz in Co8 compared to LiFeAs, but the linewidths are broader (see Fig. 1). This suggests that the contribution of 3D effects in linewidths does not dominate. Similarly, the d xz /d yz hole band in LiFeAs has a much smaller k z dispersion than the electron band (a factor 7 smaller in theory), but a similar linewidth. We conclude that the difference in linewidth behavior between d xz /d yz and d xy cannot be attributed to k z dispersion.

V. LINEWIDTH ANALYSIS A. Linewidth at the Fermi level

We now return to the discussion of Fig. 1. We consider first the different values at the Fermi level for d xz /d yz and d xy . The meaning of these differences may better appear in Fig. 5, where we have gathered values of MDC linewidths measured in different samples, different cleaves and/or different Co contents. This reveals quite a large scattering in absolute values, usually by a factor 2 and even much larger for samples cleaved at high temperatures (blue triangles). Some trend remain however clear, such as the smaller values for d xy compared to d xz /d yz or the increase of both linewidths as a function of Co doping. We do not observe a larger increase of d xy linewidth with Co content in BaFe 2 As 2 , contrary to what was recently reported in Co-doped LiFeAs [7]. On the other hand, d xy was found to be more sensitive to temperature cycles (see Fig. 1), which we attribute to the effect of impurities formed during such cycles. This is consistent with the idea of a larger or different coupling of d xy to impurities. Interestingly, the trend of the variations with Co content scales very well with expectations deduced from residual resistivity [15]. Expecting ρ ∝ 1/nτ , we plot on the right axis the residual resistivity multiplied by the number of electrons estimated in [23] (black line, we neglected the hole contribution, as suggested in ref. [15]). The similar tendency is a very compelling indication that the linewidths indeed reflect intrinsic physics of the compounds, even if there may be an additional constant term of extrinsic origin.

The linewidths in LiFeAs are somewhat smaller than the Co ones, but not dramatically. For the hole d xy band, our value is similar to that of the previous report by Kordyuk et al. [10]. In fact, our data do not support the common idea that ARPES spectra are much broader in Co-doped BaFe 2 As 2 than LiFeAs. To explain the large scattering in the absolute values of the linewidths, we infer that they contain a rather large contribution from impurities mainly created at the surface during the cleave. This also explains the much smaller values of 1/τ near E F measured by optics, namely 7 meV for Co8 [36] and 3 meV for LiFeAs [37]. However, this impurity contribution should be independent of binding energy and temperature [9], so that it can be considered as a mere offset. Indeed, we usually observe similar energy dependences of linewidths for samples exhibiting different absolute values. There is also a constant contribution from experimental resolution that we have not attempted to subtract.

Transport [15] and quantum oscillation experiments [38,39] have often suggested longer mean free paths for electrons than holes. This is not corroborated by these measurements, as we observe similar linewidths for hole and electron pockets for a given orbital. Nevertheless, as the orbital content is different for hole and electron pockets, taking into account different behaviors for d xz /d yz and d xy may reconcile these observations.

B. Dependence with binding energy and temperature

The linewidths are extracted experimentally by fitting the MDCs to lorentzians. This is much easier than fitting directly the Energy Distribution Curves (EDC), whose background is more difficult to model [8,40]. These MDC linewidths, in units of Å-1 , depend both on band properties (such as v F ) and many body effects. To isolate the contribution of many-body effects, we now compute 1/τ by multiplying the MDC linewidth by the slope of the dispersion. As long as the experimental dispersion is proportional to the theoretical one (v bare ), one can define a constant renormalization value Z = v F /v bare that simply relates the lifetime to the imaginary part of the self energy Σ .

/τ [eV ] = slope( ω)[eV. Å] * δν [1/ Å] = Z * Σ (ω, T ) (3)
As we discussed in section III, the experimental dispersions are usually well described by renormalized LDA calculation, so that we can apply this procedure. It is [15] and multiplied by the number of electrons from ref. [23].

straightforward when the dispersion is linear in the considered energy window (the two quantities are simply proportional), but would work as well for a parabolic [41] or arbitrary dispersion. In Fig. 6, we plot the MDC linewidths as a function of ω and the slopes used for normalizing them. The larger deviation from linearity is found for the d xz /d yz electron band, which affects significantly the ω dependence of 1/τ . However, we obtain a flat ω dependence for 1/τ (Fig. 1b) that echoes its flat temperature dependence (Fig. 1c) and gives confidence in this procedure. In LiFeAs, this normalization yields similar dependences for the hole and electron lifetimes on d xz /d yz , despite their very different dispersion, which again sounds correct. For d xy , the ω dependences on the two pockets are indeed similar, except for an offset that we attribute to a better quality of the sample in which holes were measured.

To go one step further, we fit the data in Fig. 1 to the ω 2 or T 2 law expected for a Fermi liquid (Eq. 1), wherever it is possible. For d xy , such fits describe well the data below ∼ 40 meV. This defines a coherent energy scale E coh below which the Fermi liquid is obeyed, as it is generally expected [42]. We obtain γ=10 for LiFeAs and 16 for Co8. The difference in γ is partly due to the more abrupt deviation of the ω 2 behavior in Co8 and may not reflect essential physics (the broadening between E F and 80meV is for example nearly identical in the two compounds). The broadening of the d xy band with increasing temperature in Co8 can be described by the same exponent γ=16. For d xz /d yz , there is no clear FL regime, the data are dominated by a constant term. We could either fit a much smaller γ (such as 2.5 in LiFeAs below 50 meV) or define a much smaller E coh (∼10 meV in Co8).

Recently, the lifetimes have been studied by optical spectroscopy both in BaFe 2 As 2 doped with 10% Co [36] and LiFeAs [37]. In both cases, the extracted lifetimes were found to follow the Fermi liquid behavior below (cd) Slopes of the dispersion used to compute 1/τ (Eq. 3) in LiFeAs (c) and Co8 (d). For electron bands, it is the slope of the renormalized LDA dispersion. For hole bands, it is straight lines, as they give a better or equivalent fit than the LDA (see Fig. 3d).

150K (i.e. πk B T=40 meV) or ω ∼50 meV. This E coh is in excellent agreement with that deduced from ARPES on d xy . For 10% Co, a value γ=4 was obtained from both energy and temperature dependence (a specific factor p should be added to Eq. 1 for the temperature dependence for optics [36,43], which was estimated to be 1.5). For LiFeAs, γ ∼ 7 would be obtained from temperature dependence, if we assume the same p value. These values are smaller than those found here for d xy , but would be compatible if averaged with the nearly flat dependence of d xz /d yz . This reinforces the idea that ARPES measures true QP lifetimes, but emphasizes that it is essential to consider the contribution of the two orbitals to interpret the results.

For a simple metal, the renormalized bandwidth ZW would set a low energy scale. One expects E coh ZW and γ ∼1/ZW [42]. In iron pnictides, the situation is of course more complicated because of the interaction between the different bands, especially the coherence depends on the Hund's couplings J between them [44,45] We are then in a situation where ZW is smaller for d xz /d yz , despite the fact that Z can be larger for these orbitals. One could wonder whether it is possible at all to establish a Fermi liquid regime in such conditions. Indeed, the broad linewidth of d xz /d yz could be interpreted as a totally incoherent behavior. Strictly speaking, QP are not defined if ω<1/τ , which is always the case here for d xz /d yz (of course, this criterion is difficult to apply strictly for ARPES, due to the unknown impurity contribution). A very small E coh ∼10 meV on d xz /d yz would explain the lack of ω dependence, and the large γ∼1/ZW, associated to it, could yield a high saturated value, larger than that of d xy near E F . This explanation has the advantage to explain directly the similar behavior of Co8 and LiFeAs, where the relative ZW of d xz /d yz and d xy are almost identical. However, it remains to be understood why the bands could be decoupled, so that the individual bandwidth becomes the relevant parameter.

VI. CONCLUSION

We have demonstrated that ARPES is an efficient tool to resolve the properties of d xy and d xz /d yz in iron pnictides. Especially, the two orbitals can be clearly separated on the electron pockets, which was rarely used, but allows accurate measurements and easy comparison between systems. Using both hole and electron pockets, we could confirm the increase of renormalization for d xy in LiFeAs compared to d xz /d yz . We further unveil an unexpected difference of behavior between the d xy and d xz /d yz lifetimes, already present in Co8, where the renormalization are similar on the two orbitals. The d xy lifetime can be described by a Fermi liquid behavior below E coh ∼40 meV with γ ∼10, while that of d xz /d yz is nearly flat as a function of temperature or binding energy, implying (i) a smaller γ or (ii) a smaller E coh .

(i) A smaller γ could underline smaller correlations in d xz /d yz compared to d xy . In this case, based on our renormalization values, one would expect an evolution from nearly identical linewidth behaviors in Co8 to different behaviors in LiFeAs. As it is not the case, this explanation is not straightforward. Moreover, if γ was smaller for d xz /d yz , there should be a much stronger impurity contribution to explain its larger linewidth at E F . The possibility of orbital-dependent scattering to impurities is certainly an important point to consider when interpreting such data [13,46]. However, in Co-doped BaFe 2 As 2 , we observe similar broadenings of d xz /d yz and d xy as a function of Co content, which rather suggests similar couplings. In fact, it is the d xy orbital that appears more sensitive to impurities, through temperature cycles in our measurements or Co doping in LiFeAs in ref [7]. Again, this does not support this scenario unambiguously.

(ii) A smaller E coh on d xz /d yz would explain at the same time the larger value at E F and the absence of well defined ω dependence. We propose that the reduced effective bandwidth for d xz /d yz (∼50 meV) compared to d xy (∼150 meV) could play a key role in defining E coh . This would directly explain the similar behavior in Co8 and LiFeAs. We hope our data will stimulate further studies of how the Fermi liquid regime can be established in a multiband systems with small and very different ef-fective bandwidths.

Qualitatively, averaging the d xy and d xz /d yz behavior revealed here reproduces quite well observations from optics or transport. However, our results suggest that it is the orbital content, rather than the hole or electron character, that is the main factor of differentiation.

Figure 2 .

 2 Figure 2. (a) Fermi Surface measured in Co8 at T=25 K, with polarization A along the red arrow and photon energy 34 eV. (a1)-(a2) Energy-momentum images of the two different cuts indicated in (a). (b) Sketch of the orbital characters on the two electron pockets expected by calculation. (c) Experimental band dispersion (symbols) compared with calculated bands shifted and renormalized as indicated in Table1.

Figure 3 .

 3 Figure 3. (a) Fermi Surface measured in LiFeAs at T=25 K, with polarization A along the red arrow and photon energy 34eV. Sketches of the different pockets with orbital characters are indicated as guides to the eye. (a1)-(a5) : energy-momentum images of the different cuts indicated on the FS map. (b) Experimental dispersions of all bands at kz=0 (points). Calculated dispersions renormalized by a factor 2 (lines). (c) Same experimental dispersions with calculation shifted and renormalized as indicated in Table 1. (d) Zoom on the dispersion of the different bands near EF . Color lines are the LDA dispersion of (c) and black line on dyz hole a linear fit.

Fig. 4 Figure 4 .

 44 Figure 4. (a) Spectral weight in LiFeAs integrated in a 5meV window around EF for the dxy electron band as a function of photon energy (cut (a2) in Fig. 3). Lines are theoretical variation for kF . (b) Same for the dxz/dyz electron band (cut (a1) in Fig. 3). Dotted line is the position of dxy. (c) Spectra corresponding to (b). (d) Experimental (points) and calculated (lines) variations of kF with photon energy in LiFeAs. (e) Same as (d) for Co8. The dotted lines indicate kF calculated after a shift of 0.1eV.

  (e). The different stacking of the FeAs slabs in Co8 yields quite different 3D effects, with a sizable energy dependence of d xy , which is nearly as large as that of d xz /d yz in theory (solid lines, the maxima are shifted in k z , because the two bands are measured in two adjacent BZ, where k z is shifted by 1). Our results are quite well reproduced by calculations, both for d xy and d xz /d yz , if we compute k F at -0.1eV (dotted line) rather than E F (solid line), to take into account the FS shrinking. The variation for d xz /d yz is again slightly reduced in experiment. The number of electrons is divided by 2 compared to the theoretical expectation. This very large shrinking in Co8, almost absent in LiFeAs, is an intriguing feature that should be further explored. It was assigned by Ortenzi et al. to interband interactions mediated by spin fluctuations [35] and could then suggest a different role of such fluctuations in the two compounds.

  Figure 6. (a-b) MDC half width at half maximum δν as a function of binding energy in LiFeAs (a) and Co8 (b).(cd) Slopes of the dispersion used to compute 1/τ (Eq. 3) in LiFeAs (c) and Co8 (d). For electron bands, it is the slope of the renormalized LDA dispersion. For hole bands, it is straight lines, as they give a better or equivalent fit than the LDA (see Fig.3d).

  Figure 1. Lifetime /τ (eV) extracted from ARPES linewidths (see section V) in (a) LiFeAs for hole and electron pockets at 25K as a function of binding energy ω (b) Co8 for electron pockets at 25K as a function of binding energy ω (c) Co8 for electron pockets at EF as a function of temperature. As we could never obtain a behavior reversible in temperature for dxy, we indicate width measured when warming up (large symbols) and cooling down (small symbols). The lines are fits to simple Fermi liquid behaviors, as discussed in section V-B. They are calculated from Eq. 1, with the values of γ indicated on the figure and adequate offsets.

	(a)	0.05	γ = 10	LiFeAs	(b)		Ba(Fe 0.92 Co 0.08 ) 2 As 2	(c)	Ba(Fe 0.92 Co 0.08 ) 2 As 2
	h /τ (eV)	0.04 0.03 0.02	γ = 2.5		0.08 0.06 0.04 h /τ (eV)	γ = 16	h /τ (eV)	0.08 0.06 0.04	γ = 16
		0.01 0.00	Holes dxz/dyz Holes dxy			0.02 0.00	Elec dxz/dyz Elec dxy		0.02 0.00	Elec dxz/dyz Elec dxy
			-0.06 -0.04 -0.02	0.00	-0.08 -0.06 -0.04 -0.02 0.00		300	200	100	0
			hω (eV)			hω (eV)		Temperature (K)
	arXiv:1602.07986v1 [cond-mat.str-el] 25 Feb 2016

  Second, there is a contrast in the behavior of d xz /d yz and d xy . There is a clear dependence as a function of ω for d xy below 50meV that is missing or very reduced for d xz /d yz . The fact that the temperature dependence is similarly larger on d xy compared to d xz /d yz in Co8 [Fig. 1(c)], following expectations from Eq. 1, reinforces the idea that the difference is intrinsic. Moreover, the absolute values at the Fermi level E F are generally larger for d xz /d yz than d xy . As the same tendencies are detected for LiFeAs and Co8, this appears as a rather universal feature, which likely has important consequences for the physics of iron pnictides.

  . Nevertheless, it seems important to keep in mind that ZW can be very small in iron pnictides, and, in fact, significantly smaller for d xz /d yz than d xy . Table I recalls that the bottom of the electron bands is located around 50 meV for d xz /d yz and 150 meV for d xy . The top of the hole bands are even closer from E F , and closer for d xz /d yz than d xy .

We thank A. van Roekeghem, S. Biermann and F. Rullier-Albenque for many useful discussions. This work benefitted from financial support from the ANR "Pnictides".


Keywords: elliptic curve, formalising mathematics, primality proving courbe elliptique, formalisation des mathématiques, preuves de primalité

Elliptic curves are fascinating mathematical objects. In this paper, we present the way they have been represented inside the COQ system, and how we have proved that the classical composition law on the points is internal and gives them a group structure. We then describe how having elliptic curves inside a prover makes it possible to derive a checker for proving the primality of natural numbers.

Introduction

Elliptic curves are a very classical topic of interest in number theory and algebraic geometry, and the basis for a very large body of theory, with lots of far-fetched generalisations. Their most famous application is, without any doubt, the proof of Fermat's last theorem. They are also commonly used in algorithmic number theory since the mid 80's for factoring integers [START_REF] Lenstra | Factoring integers with elliptic curves[END_REF] or proving primality [START_REF] Goldwasser | Almost all primes can be quickly certified[END_REF], but have also made their way more recently in cryptography, see e.g. [START_REF] Cohen | Handbook of Elliptic and Hyperelliptic Curve Cryptography[END_REF].

It is then natural to try formalising them inside a proof system. In this work, the main property we are interested in is their group structure. The difficult part of the formalisation is in proving associativity. In our formal proof, we follow an elementary algebraic approach that is particularly suited to proof systems. Our main source of inspiration has been the proof proposed by Stefan Friedl [START_REF] Friedl | An elementary proof of the group law for elliptic curves[END_REF]. The proof is very technical and consists in massaging equalities over an abstract field under various conditions. The only sour point is that some of these equalities (exactly 3) are so huge that they cannot be reasonably handled by a human being. In his paper, Stefan Friedl advocates the use of a computer algebra like COCOA [START_REF] Capani | Cocoa: Computations in commutative algebra[END_REF] to check these equalities. Reflecting these computations inside the proof system is the main difficulty of the formal proof. In this paper, we explain how this has been done inside the proof system COQ [START_REF]The Coq Proof Assistant Reference Manual[END_REF].

Once the group structure has been established, we can use elliptic curves in a very effective way to prove the primality of some natural numbers. For this, we use prime certificates: various primality proving algorithms provide the user with a compact output, a certificate, with the property that it can be easily checked, whereas generating it is usually cpu-intensive.

We have already played with the idea of prime certificates and in particular of Pocklington certificates in [START_REF] Grégoire | A Purely Functional Library for Modular Arithmetic and its Application to Certifying Large Prime Numbers[END_REF][START_REF] Grégoire | A Computational Approach to Pocklington Certificates in Type Theory[END_REF]. This allows us to prove the primality of some large numbers like the millennium prime (a prime number with exactly 2000 digits discovered by John Cosgrave) or the 27 th Mersenne prime 2 44497 -1. The only drawback of the certificates we have been using so far is that they only apply to a restricted class of prime numbers. For example, a Pocklington certificate can be generated for N if N -1 can be factored; in particular, we are not able to generate a Pocklington certificate in the case where (N -1)/2 turns out not to be prime (this can be shown by compositeness tests) but it is impossible to factor it, which is a typical case as soon as N grows large. In contrast, under some deep (typically Crámer's conjecture) conjectures of analytic number theory, elliptic curve certificates exist for any prime and can be found in polynomial time. At the moment, the largest prime for which an elliptic certificate has been generated has 20562 decimal digits. Ten days are needed to verify its certificate [START_REF] Morain | [END_REF].

The paper is organised as follows. In Section 2, we recall what elliptic curves are and explain how their group structure has been formally established. In Section 3, we focus on prime certificate and show how it is possible to formally derive a checker for elliptic curve certificates.

Elliptic curves and group structure

In algebraic geometry, an elliptic curve is defined as follows.

Definition. Let K be a commutative field. An elliptic curve E(K) is a pair (C, P ) where C is a (complete) algebraic curve of genus 1 defined over C and P a point of C(K). This definition is not very tractable in practice, and the following characterization, which follows easily from the Riemann-Roch Theorem [START_REF] Joseph | The Arithmetic of Elliptic Curves[END_REF], is rather used as a definition.

Proposition. Every elliptic curve over K admits an affine plane model of the form

y 2 + a 1 xy + a 3 y = x 3 + a 2 x 2 + a 4 x + a 6 , (1)
where the distinguished point of the definition is the projective point (0 : 1 : 0), i.e. the direction y → ∞. If the characteristic of K is greater than 5, the equation further simplifies to

y 2 = x 3 + Ax + B. (2)
Conversely, a nonsingular curve of the form (1) -or (2) if the characteristic is greater than 5defines an elliptic curve. For (2) the non-singularity assumption is equivalent to 4A 3 +27B 2 = 0. Most of the interest for elliptic curves stems from the fact that they come with a natural group structure. A composition law (P, Q) → P + Q over E(K) can be constructed as follows:

• The neuter element is the point at infinity;

• Let (P Q) be the line through P and Q (the tangent to the curve at P if P = Q).

• Let R be the third point of intersection of (P Q) and K.

• The point P + Q is the symmetrical of R with respect to the line y = 0.

An alternative way of seeing the composition law is that three aligned points add up to zero, which is the point at infinity.

Theorem 1 (E(K), +) is a group.

The only difficult part is associativity. Classically, there are two roads for proving this. The "low road" uses either computer algebra (the one that we shall carefully take later on) or elementary but subtle geometric considerations. The "high road" instead defines a group by taking the quotient of the free abelian group over E(K) by a subgroup of relations (generated by the relations "three aligned points add up to zero"), and proceeds to prove that each class contains one and only one representative consisting of a single point (so that the underlying set is in bijection with E(K)) and that the operation is the same as +. All these are consequences of the (non-trivial) Riemann-Roch Theorem.

In our formalisation, we start from [START_REF] Boyer | Single-threaded objects in ACL2[END_REF]. We take an arbitrary field K of characteristic other than 2 (this is the only assumption that is needed to prove the group structure) and two elements A and B such that 4A 3 + 27B 2 = 0. The type elt that represents the elements of the curve is defined as follows Inductive elt: Set := inf_elt: elt | curve_elt (x: K) (y: K) (H: y^2 = x^3 + A * x + B): elt.

An element of type elt is either 0 (inf_elt) or an element of the curve (curve_elt) that contains an x, a y and a proof H that the point (x, y) is on the curve.

For the operations on the curve, we need to define the opposite and the addition. Given a point p of the curve, we define -p as -If p = 0 then -p = 0 -If p = (x, y) then -p = (x, -y) Given two points p 1 and p 2 on the curve, we define p 1 + p 2 as

-If p 1 = 0 then p 1 + p 2 = p 2 -If p 2 = 0 then p 1 + p 2 = p 1 -If p 1 = -p 2 then p 1 + p 2 = 0 -If p 1 = p 2 = (x, y) and y = 0 then let l = (3x 2 + A)/2y and x 1 = l 2 -2x in p 1 + p 2 = (x 1 , -y -l(x 1 -x)).
-If p 1 = (x 1 , y 1 ) and p 2 = (x 2 , y 2 ) and x 1 = x 2 then let l = (y 2 -y 1 )/(x 2 -x 1 ) and

x 3 = l 2 -x 1 -x 2 in p 1 + p 2 = (x 3 , -y 1 -l(x 3 -x 1 )).
We also use p 1 -p 2 as a shortcut for p 1 + (-p 2 ). Note that, as points of the curve contain a proof that they belong to the curve, the definition of these operations must also incorporate a proof that these operations are internal. For example, the opposite is represented by a function opp whose type is elt → elt and its definition is

Definition opp p := match p with inf_elt => inf_elt | curve_elt x y H => curve_elt x (-y) opp_lem end.
where opp_lem represents the proof that if we know that y 2 = x 3 + Ax + B (this is H) then we have (-y)

2 = x 3 + Ax + B.
We get that 0 is the neuter element, -p is the opposite of p, and the addition is commutative directly from the definitions. Proving associativity is rather technical. A detailed description of the formal proof is given in [START_REF] Théry | Proving the group law for elliptic curves formally[END_REF]. Here, we only sketch the proof and discuss the most interesting issues.

As the definition of addition contains five cases, the proof heavily relies on case analysis. The addition p 1 + p 2 is trivial if p 1 = 0 or p 2 = 0 or p 1 = -p 2 . When the addition is non-trivial we

x 1 -x 2 = 0 ∧ x 4 -x 3 = 0 ∧ x 2 -x 3 = 0 ∧ x 5 -x 1 = 0 ∧ y 2 1 = x 3 1 + Ax 1 + B ∧ y 2 2 = x 3 2 + Ax 2 + B ∧ y 2 3 = x 3 3 + Ax 3 + B ∧ x 4 = (y 1 -y 2 ) 2 /(x 1 -x 2 ) 2 -x 1 -x 2 ∧ y 4 = -(y 1 -y 2 )/(x 1 -x 2 )(x 4 -x 1 ) -y 1 ∧ x 6 = (y 4 -y 3 ) 2 /(x 4 -x 3 ) 2 -x 4 -x 3 ∧ y 6 = -(y 4 -y 3 )/(x 4 -x 3 )(x 6 -x 3 ) -y 3 ∧ x 5 = (y 2 -y 3 ) 2 /(x 2 -x 3 ) 2 -x 2 -x 3 ∧ y 5 = -(y 2 -y 3 )/(x 2 -x 3 )(x 5 -x 2 ) -y 2 ∧ x 7 = (y 5 -y 1 ) 2 /(x 5 -x 1 ) 2 -x 5 -x 1 ∧ y 7 = -(y 5 -y 1 )/(x 5 -x 1 )(x 7 -x 1 ) -y 1 ⇒ x 6 -x 7 = 0 where (x 4 , y 4 ) = p 1 ⊕ g p 2 , (x 5 , y 5 ) = p 2 ⊕ g p 3 , (x 6 , y 6 ) = (p 1 ⊕ g p 2 ) ⊕ g p 3 and (x 7 , y 7 ) = p 1 ⊕ g (p 2 ⊕ g p 3 )
Figure 1: Condition for the x component of the first case. use the notation p 1 ⊕ p 2 . In turn, this addition can correspond to the tangent case, p 1 ⊕ t p 2 , where p 1 = p 2 or the general case, p 1 ⊕ g p 2 , where p 1 = p 2 . For associativity, it is easy to discharge the trivial cases and what remains to be proved is then that p 1 ⊕ (p 2 ⊕ p 3 ) = (p 1 ⊕ p 2 ) ⊕ p 3 . Each of these additions can be tangent or general, this can be further reduced to four main cases:

1. p 1 ⊕ g (p 2 ⊕ g p 3 ) = (p 1 ⊕ g p 2 ) ⊕ g p 3 . 2. p 1 ⊕ g (p 2 ⊕ t p 2 ) = (p 1 ⊕ g p 2 ) ⊕ g p 2 . 3. p 1 ⊕ g (p 1 ⊕ g (p 1 ⊕ t p 1 )) = (p 1 ⊕ t p 1 ) ⊕ t (p 1 ⊕ t p 1 ) 4. p 1 ⊕ g (p 2 ⊕ g (p 1 ⊕ g p 2 )) = (p 1 ⊕ g p 2 ) ⊕ t (p 1 ⊕ g p 2 )
The first three cases are discharged using explicit computations while the last case is derived using auxiliary properties such as the cancellation law, i.e. if p + p 1 = p + p 2 then p 1 = p 2 . To understand what it means to perform cases 1÷3 by explicit computations, let us consider the first case. If we take the x component of equation 1 and unfold the definition of addition, we get the conditional equation to be satisfied given in Figure 1. We can now transform this equation into a polynomial one and then use Buchberger algorithm [START_REF] Buchberger | Introduction to Gröbner Bases[END_REF] to check that the equation is satisfied. This is where the computer algebra system comes into play. Integrating Buchberger algorithm in INRIA a safe way inside a prover is possible and has already been done [START_REF] Créci | Gb: une procédure de décision pour le système Coq[END_REF][START_REF] Harrison | Complex quantifier elimination in HOL[END_REF], but cases 1 ÷ 3 can be decided by a much simpler strategy. To illustrate this strategy, let us consider a more elementary example, i.e. the proof that the addition is internal in the tangent case. It amounts to proving that

y 2 = x 3 + Ax + B ∧ l = (3x 2 + A)/2y ∧ x 1 = l 2 -2x ∧ y 1 = -y -l(x 1 -x) ∧ ⇒ y 2 1 = x 3 1 + Ax 1 + B
The first step is to get rid of division by normalising the equation into an expression of the form N/D = 0 where N and D are two polynomial expressions. We are then left with proving that N = 0. In our example, a naive normalisation gives 2 10 y 8 -2 10 y 6 x 3 -2 10 Ay 6 x -2 10 By 6 = 0

The second step is to replace all the occurrences of y 2 by x 3 + Ax + B.

2 10 (x 3 +Ax+B) 4 -2 10 (x 3 +Ax+B) 3 x 3 -2 10 A(x 3 +Ax+B) 3 x-2 10 B(x 3 +Ax+B) 3 = 0
The final step is to use distributivity, associativity, commutativity and collect equal monomials so that everything cancels out.

Note that even a dedicated computer algebra system like Maple may fail, by using the general Gröbner machinery, to prove such a large identity. The remedy, like here, is to perform "by hand" the reductions modulo the ideal < y 2 1 -x 3 1 -Ax 1 -B, y 2 2 -x 3 2 -Ax 2 -B >, by using the fact that the basis given for this ideal is already a Gröbner basis for the lex ordering

y 1 > y 2 > x 1 > x 2 .
To sum up, three ingredients are needed to automate the proof of the three lemmas corresponding to cases 1 ÷ 3: a procedure to normalise polynomial expressions (last step), a procedure to rewrite polynomial expressions (second step) and a procedure to normalise rational expressions (first step).

The procedure to normalise polynomial expressions was already present in COQ and is described in [START_REF] Grégoire | Proving ring equalities done right in Coq[END_REF]. Its main characteristic is to use an internal representation of polynomials in Horner form. This representation is unique up to variable ordering. Given a polynomial P and a variable x, we write P as P 1 + x i Q 1 where x does not occur in P 1 and is not a common factor in Q 1 and we proceed on P 1 and Q 1 recursively. The normal form is further simplified writing 0 for m0 and P for 0 + P and P + 0. For example 2 10 y 8 -2 10 y 6 x 3 -2 10 Ay 6 x -2 10 By 6 is represented as

y 6 ((B(-2 10 ) + x(A(-2 10 ) + x 2 (-2 10 ))) + y 2 2 10 )
for the ordering y > x > A > B. Once the procedures to add and multiply polynomials in Horner form have been defined, the normal form of a polynomial P is obtained by a simple structural traversal of P . As described in [START_REF] Grégoire | Proving ring equalities done right in Coq[END_REF], this leads to a very effective way of proving ring equalities by just checking that the normal form of the left side of the equality is structurally equal to the normal form of the right side of the equality.

Rewriting has been implemented in a naive way on Horner representation. It uses a simple procedure that, given a monomial m, splits a polynomial P in a pair (P 1 , Q 1 ) in such a way that P = P 1 + mQ 1 . Now rewriting once with the equation m = R is performed by building the polynomial P 1 + RQ 1 . For example, if we want to rewrite the previous polynomial with the equation y 2 x 2 = z + t, we first split the polynomial

y 6 ((B(-2 10 ) + x(A(-2 10 ) + x 2 (-2 10 ))) + y 2 2 10 ) into ( y 6 ((B(-2 10 ) + x(A(-2 10 ))) + y 2 2 10 ) , y 4 x(-2 10 ) )
We then build

(y 6 ((B(-2 10 ) + x(A(-2 10 ))) + y 2 2 10 )) + (z + t) (y 4 x(-2 10 ))
and normalise it. For rewriting with respect to a list of equations, we just iterate the single rewriting operation for all the equations in a fair way. Note that in our specific case we are always rewriting with equations of the type y 2 i = x 3 i + Ax i + B, so we take a special care in choosing a variable ordering for the Horner representation such that the y i are privileged. By doing so, the splitting procedure has only to visit a small part of the polynomial.

Normalising rational equalities is not as simple as for polynomial expressions. When reducing the initial expression to a common denominator, the expressions for the numerator and the denominator can grow very quickly. Our experiments with the proofs of cases 1 ÷ 3 have shown that getting a normalised numerator of a reasonable size was mandatory in order to be able to complete the following rewriting phase. So, when normalising, some reductions are needed. For example, when normalising P 1 /Q 1 + P 2 /Q 2 we can do much better than building the naive fraction

(P 1 Q 2 + P 2 Q 1 )/Q 1 Q 2 .
Unfortunately, gcd algorithms for multivariate polynomials are far more complex to implement than the simple Euclidean algorithm for univariate polynomials. So, we have just implemented a syntactic heuristic. Our heuristic factorises the products that are in common between Q 1 and Q 2 by a simple structural traversal of the lists of products that compose Q 1 and Q 2 . For example, normalising 1/x + 1/xy + 1/y = 0 gives y + 1 + x = 0.

The three procedures (polynomial normaliser, rewriter and rational normaliser) are put together in a single tactic field. This tactic works on arbitrary field structures. Given F a rational expression, calling the tactic field[H 1 H 2 ] attempts to solve the goal F = 0 using the rewriting rules H 1 and H 2 . The three procedures have been defined inside the COQ logic using the two-level approach [START_REF] Boutin | Using Reflection to Build Efficient and Certified Decision Procedures[END_REF]. Their correctness can then be stated inside the proof system and formally proved. This ensures once and for all applications that the field tactic only performs valid simplifications. A key aspect of this tactic is its efficiency. Proving the group law requires non-trivial computation. Having a relatively efficient procedure is mandatory to be able to complete the proof. On an Intel Xeon (2.66 GHz) with 2Gb of RAM, compiling the definition of elliptic curves up to associativity takes one minute and twenty seconds in COQ.

Prime certificate

The group structure on E(Z/pZ) is used in the same way as the group structure on (Z/pZ) * is used in Pocklington's p -1 test, which we recall: Proposition 1 (Pocklington) Let N be an integer. Assume that there exists a coprime to N and s such that

• a s = 1 mod N ;
• for all prime divisor p of s, one has gcd (a s/p -1, N) = 1.

Then, for any prime q|N, one has q = 1 mod s, so that if s ≥ √ N, N is prime.

The main drawback of Pocklington's test is the fact that if one is not able to factor N -1 to some extent (a little thought shows that the first condition implies, if N is actually prime, s|N -1, since the conditions stated imply that s is the order of a mod N), then there is no way one can apply the theorem.

The following theorem is thus much more versatile, the main point being that, under reasonable conjectures, we are always able to find, in polynomial time, a curve over Z/NZ with cardinality which can be factored.

We make use of points in projective coordinates, which is required to give a meaningful statement. This allows to give formulas without division for the group law, which thus remain meaningful in any ring -since Z/NZ can be assumed to be a field only after N has been proved prime... Proposition 2 (Goldwasser-Kilian) Let N be an integer. Assume that there exist an elliptic curve y 2 = x 3 + Ax + B with A, B ∈ Z and gcd (4A 3 + 27B 2 , N) = 1, a point P = (x P : y P : 1) such that y 2 P = x 3 P + Ax P + B mod N, and an integer s such that • s.P = (0 : 1 : 0) mod N ;

• for all prime p|s, (s/p).P = (x p : y p : z p ) mod N with gcd(z p , N) = 1.

Then, for all prime q|N, we have the number of points of E(Z/qZ) is such that #E(Z/qZ) = 0 mod s.

In view of the Hasse-Weil bound [START_REF] Joseph | The Arithmetic of Elliptic Curves[END_REF], which states that for all prime q and elliptic curve E over Z/qZ, one has #E(Z/qZ) ≤ ( √ q + 1) 2 , we see that N is prime as soon as s > √ q + 1. In the sequel, we shall use the slightly weaker but much easier bound #E(Z/qZ) ≤ 2q + 1, which shows that N is prime as soon as s > √ 2q + 1.

An elliptic curve certificate is thus defined recursively by:

C N := (N, A, B, x P , y P , s, (q 1 , C q 1 ), . . . , (q k , C q k )),

where the q i are the prime factors of s, and checking the certificate amounts to checking the hypotheses of the theorem, hence modular arithmetic.

Computing a certificate is a much harder problem. The most efficient way is to use the theory of complex multiplication, after Atkin and Morain. For a more extensive bibliography, the interested reader can refer to the nice survey of primality proving by Morain [START_REF] Morain ; Adleman | La primalité en temps polynomial[END_REF].

The idea behind the formal proof of Proposition 2 is to relate computations with projective coordinates done in Z/NZ with the same computations with standard coordinates done in Z/qZ when q is a prime divisor of N. We use a non-standard definition of projective coordinates where we still single out the zero. A point is either 0 or a triplet (x : y : z). We also define the projection t q as t q (0) = 0 and t q ((x : y : z)) = (x mod q/z mod q, y mod q/z mod q). Note that the projection in the second case is possible only if z mod q has effectively an inverse in Z/qZ for prime divisor q of N, so only if gcd (n, N ) = 1 . In the following, we write the property of n having an inverse as inv N (n).

We have to be careful when computing in Z/NZ since we want the computation in Z/NZ to project faithfully to the same computation in Z/qZ. Suppose that we test to zero a variable n. Obviously if we have n = 0, we also have n mod q = 0, but the converse may not be true: n = 0 does not imply that n mod q = 0. So every time we have a test to zero on a variable n, we have to check the extra condition that inv N (n) to ensure that the computation on Z/qZ would take the same branch. For this reason, all the functions that we have defined for Z/NZ take an extra argument that represents the side conditions so far and return a pair composed of the result plus the new side conditions. We encode all the side conditions in a single number using the property that inv N (mn) implies inv N (m) and inv N (n).

A result (p 1 , sc 1 ) of a computation can be correctly interpreted only if we have inv N (sc 1 ) when p 1 = 0 or inv N (z 1 sc 1 ) when p 1 = (x 1 : y 1 : z 1 ). In the following, we use the notation ⌈p 1 , sc 1 ⌉ to indicate that this property holds. To give a concrete example, consider the function double that doubles a point p 1 under the initial side condition sc 1 . Its definition is

double(p 1 , sc 1 ) = if p 1 = 0 then (0, sc 1 ) else let (x 1 : y 1 : z 1 ) = p 1 in if y 1 = 0 then (0, z 1 sc 1 ) else let m = 3x 2 1 + Az 2 1 and l = 2y 1 z 1 in let l 2 = l 2 and x 2 = m 2 z 1 -2x 1 l 2 in ((x 2 l : l 2 (x 1 m -y 1 l) -x 2 m : z 1 l 2 l), sc 1 )
and its correctness is stated as

∀p 1 sc 1 p 2 sc 2 , let (p 2 , sc 2 ) = double(p 1 , sc 1 ) in ⌈p 2 , sc 2 ⌉ ⇒ t q (p 2 ) = 2.t q (p 1 )
Another key property is that the final side condition includes the initial one

∀p 1 sc 1 p 2 sc 2 , let(p 2 , sc 2 ) = double(p 1 , sc 1 ) in ⌈p 2 , sc 2 ⌉ ⇒ ⌈p 1 , sc 1 ⌉
Knowing that these two properties must hold, we can now explain how the side conditions have been generated in the body of the function double. When p 1 is zero, no condition has to be added. When p 1 is (x 1 : 0 : z 1 ), the result is zero so we include the final side condition that z 1 must be invertible. Finally, when p 1 is (x 1 : y 1 : z 1 ) with y 1 = 0, a defensive final side conditions would be y 1 z 1 sc 1 to insure that y 1 and z 1 are invertible but this is not necessary since if we look at the z component, z 1 l 2 l, of the resulting point, that contains already the product y 1 z 1 .

The function add that adds two points is defined in a similar way. With the two functions double and add, we can define the function scal that adds n times the point p 1 under the side condition sc 1 :

scal (n, p 1 , sc 1 ) = if n = 1 then (p 1 , sc 1 ) else if n is even then let (p 2 , sc 2 ) = double(p 1 , sc 1 ) in scal (n/2, p 2 , sc 2 ) else let (p 2 , sc 2 ) = double(p 1 , sc 1 ) in let (p 3 , sc 3 ) = scal ((n -1)/2, p 2 , sc 2 ) in add(p 1 , p 3 , sc 3 )
and its correctness is stated as

∀n p 1 sc 1 , p 2 sc 2 , let (p 2 , sc 2 ) = scal(n, p 1 sc 1 ) in ⌈p 2 , sc 2 ⌉ ⇒ t p (p 2 ) = n.t p (p 1 )
Given a list l = [(q 1 , α 1 ), . . . , (q n , α n )] such that s = q α 1 1 . . . q αn n , checking that s.a = 0 and (s/q i ).a = 0 for 1 ≤ i ≤ n is done by factorising computations as much as possible. We first compute p 1 = (s/(q 1 . . . q n )).a. The function scalL that, given a point p 1 and a list l such that l = [q 1 , . . . , q n ] and a side condition sc 1 , verifies that (q 1 . . . q n ).p 1 = 0 and (q 1 . . . q n /q i ).p 1 = 0 for every q i is defined as follows

scalL(l, p 1 , sc 1 ) = if length(l ) = 0 then (p 1 , sc 1 ) else let (p 2 , sc 2 ) = scal (hd (l), p 1 , sc 1 ) in let (p 3 , sc 3 ) = scal _list(tl (l), p 1 , sc 2 ) in if p 3 = 0 then (0, 0) else let (x 3 , y 3 , z 3 ) = p 3 in scalL(tl (l), p 2 , z 3 sc 3 )
where the functions hd and tl return the head element and the tail of a list respectively and the function scal_list applied to l and p 1 , where l = [q i , . . . , q n ], computes (q i . . . q n ).p 1 . In every recursive call of the function scalL, p 3 represents one of the (q 1 . . . q n /q i ).p 1 . If p 3 is 0, we put 0 in the side condition in order to invalidate the final result. If p 3 is not 0, t p (p 3 ) must also be non-zero so we add its z component into the side condition.

The actual representation of COQ certificates is a list of elementary certificates rather than a tree. This list must be well-founded in the sense that the primality of a number that is needed by an elementary certificate at position i is justified by another certificate at position j > i. An elementary certificate is either an immediate certificate, a Pocklington certificate or an elliptic certificate. An immediate certificate is for numbers whose primality is given by a predefined theorem. In our library there are predefined theorems for the first 5000 prime numbers. Our elliptic certificates slightly differ from the standard ones to accommodate the certificates produced by the tool PRIMO [START_REF] Martin | PRIMO -primality proving[END_REF]. An elliptic certificate is composed of seven elements: {N, A, B, x 1 , y 1 , n 0 , [(q 1 , α 1 ), . . . , (q n , α n )]}:

-N is the number to be proved prime; -A and B are the parameters of the elliptic curve y 2 = x 3 + Ax + B; -x 1 and y 1 are the coordinates of the initial point, in projective coordinates the initial point is p 0 = (x 1 : y 1 : 1);

n 0 is an initial scalar, so the point whose order is checked is n 0 .p 0 .

-[(q 1 , α 1 ), . . . , (q n , α n )] is a complete factorisation of the order s of n 0 .p 0 , i.e. s = q α 1 1 . . . q αn n where all the q i are prime.

An example of such a certificate is { 329719147332060395689499, -

] } It allows to assess the primality of 329719147332060395689499 knowing that 40165264598163841 is prime with the curve y 2 = x 3 -94080x + 9834496 and the point 8209062.(0, 3136) whose order is 40165264598163841.

The function testCertif that checks an elliptic certificate c is defined as follows:

testCertif (c) = let (N, A, B, x 1 , y 1 , n, l) = c in let p 0 = (x 1 , y 1 , 1) in let sc 0 = 4A 3 + 27B 2 in INRIA let (p 1 , sc 1 ) = scal (n, p 0 , sc 0 ) in let s = mull (l) in let (n ′ , l ′ ) = split(l) in let (p 2 , sc 2 ) = scal (n ′ , p 1 , sc 1 ) in let (p 3 , sc 3 ) = scalL(l ′ , p 2 , sc 2 ) in
if N > 1 and odd (N) and 4N < (s -1) 2 and y 2 1 mod N = (x 3 1 + Ax 1 + B) mod N and p 3 = 0 and gcd (sc3, N) = 1 then true else false where, if l = [(q 1 , α 1 ), . . . , (q n , α n )], the function mull applied to l returns q α 1 1 . . . q αn n and the function split applied to l returns (q α 1 -1 1 . . . q αn-1 n , [q 1 , . . . , q n ]). The correctness of the function testCertif is stated as

∀c, (testCertif (c) = true ∧ (∀q α, (q, α) ∈ getL(c) ⇒ prime(q))) ⇒ prime(getN (c))
where getN and getL are the accessor functions for the first component and the last component of the certificate respectively.

Proving the correctness of the function testCertif is straightforward. It is standard program verification. Getting the side conditions right was the only tricky part. Note that it would be very difficult to detect a missing side condition just by testing. The program only returns a boolean value and it is mainly applied to certificates for which a positive answer is expected. Missing a side condition never invalidates valid certificates. Compared to standard checkers, we had to compromise with the bound on the order s of an element over Z/qZ. We use the naive bound 4N < (s -1) 2 that is easy to establish. Standard checkers use the sharper Hasse-Weil bound. Getting this bound formally is a challenge on its own.

Some benchmarks

Figure 2 shows some benchmarks on prime numbers from 25 to 250 decimal digits. These benchmarks have been done on a processor Intel Xeon (2.66 GHz) with 2Gb of RAM. For each number length, there are nine examples. We generate the i th example of length n by starting with the digit i and repetitively incrementing the previous digit by 1 modulo 10 to get the next digit. Once we have n digits, we just take the first prime number above this number. For example, the nine prime numbers with 25 digits are: To generate certificates, we automatically convert the ones produced by the program PRIMO.

For each example, we give the number of Pocklington certificates and elliptic certificates that compose the overall certificate and the time in seconds that COQ needs to verify the certificate. Out of the 1538 elliptic certificates only two could not be proved inside COQ. They are indicated by a star in Figure 2. This comes from the fact that our condition 4N < (s -1) 2 is stronger than the usual bound given by Hasse Theorem. The two certificates that fail were for N = 896191029759386718510467381341, s = 1501356318335977 and N = 13524493647665641984723, s = 95147824089 respectively. Since in each case the value of N is relatively small, we could replace the elliptic certificate with an alternative Pocklington certificate. Note that this should happen only for quite small values of N, so that the corresponding parts of the certificate can always be replaced by Pocklington certificates. Another strategy, if one has access to the code looking for the certificates, is to modify it very slightly to make it generate only certificates that can be checked using the modified Hasse bound (which, somehow, corresponds to making the search for certificates "less aggressive" and content itself with certificates of "lesser quality").

All the examples and the code are available at http://coqprime.gforge.inria.fr/ They compile with the current version 8.1 of the COQ system. The time to check a certificate grows quickly with the number of digits, for example we need one hour to check a number with 250 decimal digits. Note that the theory predicts that the worst-case complexity for checking a certificate for N with an underlying quadratic multiprecision arithmetic is O((log N) 4 ), which more or less corresponds to what we observe. Note also that it would be easy to gain a small constant factor in the current implementation by optimizing our function scal, eg. by using sliding window methods and/or addition-substraction chains.

In order to be able to check larger numbers, we would need two main improvements in the COQ system. First, all our computations are done inside COQ with user-defined data-structures, i.e. we use the modular arithmetic defined in [START_REF] Grégoire | A Computational Approach to Pocklington Certificates in Type Theory[END_REF] based on a datatype composed of 256 constructors that simulates an 8-bit arithmetic. We would need a direct access to the machine 32-bit arithmetic instead. In [START_REF] Grégoire | A Computational Approach to Pocklington Certificates in Type Theory[END_REF], some tests for Pocklington certificates with a native 32-bits arithmetic exhibit a speed-up of 80. Such speed-up would make it possible, in our case, to verify a 250-digit number in less than a minute. We have also experienced with extracting our code to OCAML [START_REF] Leroy | Objective Caml[END_REF] and linking it with the arbitrary-precision integer library BIGNUM [START_REF] Ménissier-Morain | The CAML Numbers Reference Manual[END_REF]. We got in this case an even bigger speed-up since checking one of our 250-digit numbers only takes 0.2 second and we were capable to check in 36 seconds the millennium prime, a prime with 2000 decimal digits whose certificate generated by PRIMO contains 245 elliptic certificates and 8 Pocklington ones.

Second, we are in a purely functional setting and a large amount of time is thus spent in allocating memory for new numbers. Computations should use constant memory to be really efficient. For this, we would need COQ to accommodate destructive datastructures as proposed for ACL2 in [START_REF] Boyer | Single-threaded objects in ACL2[END_REF] with a monadic approach.

Acknowledgments

Many people made this work possible. Joe Hurd was the first to draw our attention to formalising elliptic curves inside a prover [START_REF] Hurd | Formalized elliptic curve cryptography[END_REF]. After unsuccessful attempts, John Harrison revitalised our interest showing us that he could automatically solve the equation of the x component of the generic case in less than 3 minutes inside HOLLIGHT [START_REF] Harrison | HOL light: A tutorial introduction[END_REF] with his integrated version of Buchberger algorithm. Bruno Barras and Benjamin Grégoire helped us with their expertise in integrating the field tactic into COQ. Finally, Marcel Martin modified the PRIMO system adding an option to disable N + 1 certificates that we are not yet capable to handle.

Conclusion

We believe that what is presented in this paper is a very nice illustration of a key aspect of formal verification. Defining prime numbers from scratch is very simple in a proof assistant. It is done with five definitions in COQ: one for the natural numbers, one for the addition, one for the multiplication, one for the divisibility and finally one for the primality. All these definitions are straightforward, so it is easy to get convinced that the notion of primality has been correctly captured in the proof system. Then, the fact that there exists a formal proof for each of the numbers of Figure 2 ensures that the primality of these numbers follows logically from these five definitions. For the actual proofs, we are free to use any elaborate devices like in our case elliptic certificates. Of course, this is true only if the proof system is sound. It is then crucial to keep the trusted base of the prover as small as possible. This is what we have done here, reflecting the symbolic manipulations and internally checking the prime certificates generated by PRIMO. Another aspect that this example illustrates is how important is for a proof system not only to reason right but also to compute right. Our initial attempts in proving cases 1 ÷ 3 have clearly indicated that we were on the fringe of what was actually possible to prove inside COQ: the parser had problems to handle the huge terms we were generating and using the standard rewriting tactic of COQ was generating titanic proof terms. Reflecting symbolic manipulation was the key point that made everything possible. All we have developed was generic enough to be integrated into the COQ system. Nevertheless, a way to look at what we have done is to consider that we have just been reflecting a specific rewriting strategy. It would be interesting to investigate how we could reflect a generic rewriting strategy instead and get back our symbolic manipulation by simple instantiation.

Finally, the initial motivation that has lead to this work was to see how far we could go in proving primality inside a proof system. With elliptic certificates, we have reached the current state of the art. We can now prove the primality of any number with less than 300 decimal digits. We believe that, with the future improvement of numerical evaluation inside COQ, this limit will quickly increase. Anyway, the infrastructure is there. In particular, this is, to our knowledge, the first time that associativity has been formally proved inside a proof system. We hope that this initial effort will motivate deeper mathematical formalisations with elliptic curves inside proof systems.

RR n°6155


Fabien Teytaud

Olivier Teytaud

QR mutations improve many evolution strategies -a lot on highly multimodal problems

. Quasirandom mutations have many advantages, in particular their application is stable, efficient and easy to use. In this article, we extend this principle by applying quasi-random mutations on several well known continuous evolutionary algorithms (SA, CMSA, CMA) and do it on several old and new test functions, and with several criteria. The results point out a clear improvement compared to the baseline, in all cases, and in particular for moderate computational budget.

Introduction: derandomized evolutionary algorithms

Several works have been devoted to the derandomization of evolutionary algorithms, typically using quasi-random (QR) methods. QR numbers are designed in order to be more uniformly distributed than classic pseudo-random numbers. We here focus on the continuous case. Some works apply QR sequences for the initialization of the individuals [START_REF] Kimura | Genetic algorithms using low-discrepancy sequences[END_REF] or for the mutations [START_REF] Auger | Xse: quasi-random mutations for evolution strategies[END_REF][START_REF] Teytaud | Dcma: yet another derandomization in covariance-matrixadaptation[END_REF]. [START_REF] Teytaud | When does quasi-random work?[END_REF] pointed out that QR mutations are moderately better than pseudo-random ones in several settings. Inspired by [START_REF] L'ecuyer | Recent Advances in Randomized Quasi-Monte Carlo Methods[END_REF], [START_REF] Teytaud | When does quasi-random work?[END_REF] concluded that in order to have good QR mutations it is important to have good QR sequences. This is particularly true for optimization of strongly multimodal functions [START_REF] Chaslot | On the huge benefit of quasirandom mutations for multimodal optimization with application to grid-based tuning of neurocontrollers[END_REF][START_REF] Teytaud | When does quasi-random work?[END_REF]. Another application of QR sequences is for restarts algorithms [START_REF] Teytaud | A rigorous runtime analysis for quasirandom restarts and decreasing stepsize[END_REF] ; in this article, authors show that QR restarts outperform classic restart algorithms for multi-modal optimization. In this article, we point out the generality of the use of QR mutations by extending results from [START_REF] Teytaud | Dcma: yet another derandomization in covariance-matrixadaptation[END_REF] to more algorithms, including CMSA and SA. We also discuss various criteria and include new highly multimodal objective functions. Recently, [START_REF] Berthier | Experiments on the cec 2015 expensive optimization testbed[END_REF] experiments several algorithms on the CEC2015 testbed and notes that all algorithms which perform the best are the ones with QR mutations.

Algorithms

Generating quasi-random sequences in the square [0, 1] N is classical, many toolboxes can be found on the web. Converting from the uniform distribution in [0, 1] N to the standard Gaussian distribution can be done by composition, coordinatewise, with the inverse cumulative distribution function of the standard Gaussian; in Octave/Matlab notation, quasiRandomGaussian = norminv(quasiRandomInT he01Interval). More details around this can be found in [START_REF] Teytaud | Dcma: yet another derandomization in covariance-matrixadaptation[END_REF].

Our modification is generic, in the sense that it can be applied to many evolution strategies. QRCMA refers to CMA with quasi-random mutations. QRSA refers to SA with quasi-random mutations. QRCMSA refers to CMSA with quasirandom mutations. SA is defined in [START_REF] Rechenberg | Evolutionstrategie: Optimierung Technischer Systeme nach Prinzipien des Biologischen Evolution[END_REF][START_REF] Schwefel | Adaptive Mechanismen in der biologischen Evolution und ihr Einfluss auf die Evolutionsgeschwindigkeit[END_REF]. CMA and CMSA are respectively defined in [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF] and [START_REF] Beyer | Covariance matrix adaptation revisited -the CMSA evolution strategy[END_REF].

We also use the (1 + 1)-ES with one-fifth rule, multiplying the step-size by 1.5 in case of success and dividing by 1.5 1 4 otherwise.

Objective functions

In this study we experiment 4 unimodel functions (sphere, cigar, Ellipsoid and Schwefel) and 8 multimodal functions (Griewank, Rastrigin, Schaffer, Hump, Nsine, Randcuts and Randsines). The Schwefelmult function has been rescaled (f rescaled (x) = f (100x)) so that its range matches the usual range of other functions. These functions are basically the same as in the CMAES source code, except that the SchwefelMult function is rescaled by a factor 100 so that its relevant range is similar to the one of other functions. The hump function [START_REF] Singh | Comparison of multi-modal optimization algorithms based on evolutionary algorithms[END_REF] has been slightly modified (the optima are distributed as a Gaussian random variable with standard deviation 5 and one of the optimum has been slightly improved so that there is a unique optimum. The Nsine function [START_REF] Deb | An investigation of niche and species formation in genetic function optimization[END_REF] has been slightly modified so that the optimum is unique. The randcut function is defined

as f (x) = 30 i=1 3b i (a i .(x i -x * i ) -a i .(x i -x * i ) ) 2 ) + N ||x|| 2 10000
and the randsines as

f (x) = 30 i=1 3b i sin(a i .(x i -x * i )) 2 + N ||x -x * || 2 10000

Experiments

In all our experiments, we measure (i) the simple regret, i.e. the fitness at the proposed individual minus the fitness at the optimum (ii) for various thresholds , the frequency at which the simple regret was greater than . All experiments are reproduced 99 times, moreover, in order to improve the statistical significance, for all functions which have a random optimum, the 99 randomly drawn optima are the same for all algorithms (statistical pairing). The starting point is 0, and the optimum is randomly drawn, for each coordinate, as U ([-2, 2]) + σN, where U (I) is the uniform random variable on interval I, N is the standard Gaussian random variable, and σ = 3 is the initial step-size used in all algorithms. Because of space limitations, results are in http://www.lri.fr/ ~teytaud/ qrr.pdf.

Additional results and downloads

The experiments were reproduced with random full rotations (built with Gram-Schmidt transformation) without changing the results (figures available at http: //www.lri.fr/ ~teytaud/withrot.pdf, except for some functions, in the case of the results with highest precision and lowest dimension -for the Cigar function in dimension 5 and 10, CMA and CMA(r) perform better without restarts for the frequency of solving with precision 1 (all other results are preserved).

Conclusion

Previous results have shown that CMA is improved by QR mutations on BBOB. We here extend these results to more algorithms (CMSA, SA, CMA, and CMA with restarts), showing that QR works in most cases; additional functions, showing a bigger difference on highly multimodal functions.

CMA(R) (CMA with restarts), CMA, SA and CMSA are tested with random mutations and with quasi-random mutations. QR typically brings a significant improvement -there are cases in which the use of quasi-random numbers is the main difference between tested algorithms, i.e. all ES with QR outperform all ES without QR. Some multimodal functions from the global optimization literature (such that the hump and the NSine functions), challenging and not that usual in testbeds, and new multimodal functions have been included in tests, as well as randomly constructed multimodal functions -with these highly multimodal functions (hump, Schaffer, SchwefelMult, RandomCuts, RandomSines) QR sometimes switches from 80% failure to a negligible proportion of failure.

Furtherwork. Restarts should also benefit from quasi-random numbers; this is left as further work. SchwefelMult is an interesting stable counterexample -CMA outperforms QRCMA in a stable manner -this is not explained.