N J C Spooner
email: n.spooner@sheffield.ac.uk

P Majewski

D Muna

D P Snowden-Ifft

Simulations of the Nuclear Recoil Head-Tail Signature in Gases Relevant to Directional Dark Matter Searches

Keywords: Dark matter, Directional detector, Gas detector, WIMP, neutralino, TPC, Simulation, SRIM

We present the first detailed simulations of the so-called head-tail effect of nuclear recoils in gas, the presence of which is vital to directional WIMP Dark Matter searches. We include comparison simulations of the range and straggling of carbon, sulphur and fluorine recoils in low pressure gas. However, the prime focus is a detailed investigation of carbon and sulphur recoils in 40 Torr negative ion carbon disulfide, a gas proposed for use in large scale directional detectors. The focus is to determine whether the location of the majority of the ionization charge released and observed from a recoil track in carbon disulfide is at the beginning (tail) of the track, at the end (head) or evenly distributed. We used the SRIM simulation program, together with a purpose-written Monte Carlo generator to model production of ionizing pairs, diffusion and basic readout geometries relevant to potential real detector scenarios, such as under development for the DRIFT experiment. The results indicate the likely existence of a head-tail track asymmetry but with a magnitude critically influenced by several competing factors, notably the W-value assumed, the drift distance and diffusion, and the recoil energy.

Introduction

Measurement of the direction of a low energy (¤1 keV/amu) nuclear recoil track and the ionization charge distribution along it, resulting from the elastic scattering of a target nucleus by an incoming WIMP (Weakly Interacting Massive Particle), provides a route towards the unambiguous identification of WIMPs as being responsible for the galactic dark matter [1]. Amongst current radiation detection technologies only Time Projection Chambers filled with low pressure gas, appear capable of such a measurement. This concept has been utilized by the directional Dark Matter search experiments DRIFT [2], NEWAGE [3], MIMAC [4] and DM-TPC [5] using Time Projection Chambers (TPCs) filled with low pressure CS 2 and CF 4 . In these detectors attempts are currently made to reconstruct the orientation of the low energy nuclear recoil tracks, typically of a few millimetres in length. However, due to the character of the electronic and nuclear stopping powers of low energy nuclear recoils in gas, an asymmetric ionization charge distribution along their tracks may also be expected. Thus additional information on the absolute direction of recoils might be available. Such potential information on the track sense is termed the head-tail effect, for instance describing the location of the majority of the ionization charge as being either at the beginning half (tail) or at the end half (head) of the track. It is known that if this information can clearly be extracted in a detector this would, by breaking the forward-back degeneracy in the direction of recoils from WIMPs, have a dramatic impact on the potential directional sensitivity to dark matter. In particular, it could mean that of order 10 fewer detected WIMPs would be needed to identify a directional signal [6].

Realisation of this gain would greatly increase the feasibility of building a large dedicated low pressure TPC detector capable of a definitive identification of dark matter. It is thus vital to understand whether, and to what magnitude, the headtail effect is likely to be present, in both principle and practice. Recently some experimental evidence for head-tail asymmetry has been observed for F recoils in 100-380 Torr CF 4 [5]. Observation of the effect there was possible due to the use of relatively high energy recoils (I200 keV) and of a short drift gap to minimize diffusion. However, a realistic directional detector for dark matter will require sensitivity to much lower energies and with much larger volumes, where diffusion will likely be a significant issue. It is recognized that use of negative ion CS 2 gas can provide a route to large detectors here since in this case the charge carrier in the gas is CS 2 ions rather than electrons. The result is much lower diffusion during transit to the readout, so that drift distances can be longer and hence detector volumes much larger (order m 3 ) [2].

Measurement of the effect in a more realistic scenario for dark matter searches such as this is more difficult to perform. This is due to the limited spatial resolution imposed by the current readout technologies used with larger detectors, namely the Multiwire Proportional Counters (MWPCs) as in DRIFT, but also due to the greater ionization charge diffusion tolerated in order to maximize the size of the detectors. Some preliminary experimental measurements have nevertheless been possible so far with DRIFT itself [7]. These results do indicate a small head-tail effect even at S-recoil energies below 40 keV. However, in view of the importance of the issue for directional dark matter detection and the need to understand better the physical processes, so as to help optimise implementation of possible head-tail discrimination in a large detector with improved readout resolution, it is useful to undertake detailed simulations of the head-tail response and start to make theoretical predictions. This is the objective of the work presented here. In this light it should be noted that although comparison with the limited available head-tail data from DRIFT is of interest (see sec. 3) the purpose here is not to detail a complete head-tail simulation for a specific detector such as DRIFT II but to explore the issue in a more generic fashion.

So far some basic theoretical predictions of the head-tail effect in binary gases have been performed by A. Hitachi, i.e. [START_REF] Hitachi | A Bragg-like curve for the dark matter searches; binary gases[END_REF]. These results were based on Linear Energy Transfer calculations and ion projected range estimations from the SRIM (Stopping and Range of Ions in Matter) software package [START_REF] Ziegler | The Stopping and Range of Ions in Matter[END_REF]. The obtained Bragglike curves indicated that more ionisation charge is produced at the tail than at the head of the track, regardless of the ion type or energy. However, that work did not attempt to account for straggling, or address the issue of diffusion and the influence of track reconstruction geometry. In this work we present first detailed results of simulations of the energy loss and ionization charge distribution along tracks of carbon and sulfur ions in 40 Torr CS 2 (ρ = 1.67 10 4 gcm 3 ), together with, for completeness, some comparison results for the ionization energy loss and straggling for fluorine ions in CF 4 (ρ = 4.73 10 4 gcm 3 ). The head-tail effect is studied as a function of ion energy and diffusion in C and S but also in relation to the effect of straggling.

Simulation procedures and SRIM results

The basis of this work is firstly the SRIM software package and its incorporated TRIM (Transport of Ions in Matter) track generator Monte Carlo program. SRIM is a code commonly, and primarily, used for predictions of parameters relevant to ion implantation, sputtering and transportation in solid materials. However, it is also successfully used to simulate interactions in gases. Example comparison between ion range measurements in gases and results from SRIM calculations are shown in Table 1. SRIM was used here to produce tables as a function of energy of both the electronic (S e ) and nuclear (S n ) energy loss, range, longitudinal and lateral straggling of nuclear recoils (ions) in the gaseous target chosen. In this context straggling is defined, as used within TRIM, as the second moment of the distribution. The TRIM component is then used to generate many individual simulated tracks. This allows initial investigation of the energy loss distribution along individual raw recoil ion tracks and their secondaries, on an event-by-event basis. To complete the process, to give relevance to real detector scenarios, it is then necessary to convert, using the known average energy to create electron-ion pairs (W-value), the energy loss along each recoil track into ionization charges, or number of Negative Ion Pairs (NIPs), distributed along the track. Finally, it is necessary to account for the effects of diffusion of the tracks through a realistic drift volume (taken here to be up to 50 cm) and of the effects of projection onto readout axes. An in-house Monte Carlo was written for these latter stages, with appropriate values for the W to NIPs conversion factor, changing as a function of ion energy.

In this work we calculate parameters for the selected recoil ions using energies up to 500 keV. Fig. 1 and Fig. 2 show SRIM results for the energy loss, lateral and longitudinal straggling, plotted against recoil ion range for ion energies selected at 50 keV intervals. As can be seen in these figures, at higher energies the electronic energy loss is dominant and decreases with decreasing energy, whereas at lower values the nuclear energy loss becomes greater than the electronic and increases with decreasing energy, at least initially. For example, for sulphur ions the nuclear energy loss starts to dominate at 100 keV. This means that the energy loss along the ion track is not continuously decreasing, as might naively be expected.

Due to continuous collisions with gas atoms the direction of moving ions also deviates from their original path. This causes fluctuations in the recoil ions range described by the lateral and longitudinal straggling as a function of ion energy. This parameter, divided by the range, is shown in Fig 2 . Here one can see that as the ions slow down they experience increasing straggling relative to the drift length they have ahead, rising sharply at the very end of the track.

The amount of energy loss due to nuclear interactions as a function of ion energy is presented in Fig 3 . The difference in nuclear energy loss between the different ions is clearly seen here. For carbon ions 50% or more of the total energy is lost via nuclear interactions for energies E H 20 keV, whereas in fluorine and sulphur the equivalent factor occurs at much higher energy, 90 and 200 keV, respectively.

Tables of parameters such as those above, calculated using SRIM, are used in the Monte Carlo generator component of SRIM called TRIM to generate tracks. In addition to on-line histogramming of various quantities, TRIM results are also recorded event by event in two different formats. In the first of these the current energy, position in 3D and type of secondary recoils with their initial energy is recorded for every ion collision. For the latter, information about the energies and positions of all created nuclear recoils is added. TRIM also creates an output file which delivers information on the primary ion electronic and nuclear energy loss along the track with high spatial resolution. This format was chosen here for analysis of the head-tail effect.

On this basis, for each recoil ion with a given energy, 10 4 events were generated, recorded and further analysed. The energies used were for fluorine: from 100 to 400, with a step of 50 keV; sulphur: from 10 to 100, with a step of 10 keV, then up to 250 with a step of 25 keV and up to 400 with a step of 50 keV, and carbon from 10 to 100, with a step of 10 keV, then up to 200 keV with a step of 25 keV and up to 300 with a step of 50 keV. Each track event in the TRIM output file contained on average 100 collision points with information on the current ion energy, the 3D coordinate of the collision point and the electronic energy loss.

In order to generate the number of electron-ion pairs produced along a recoil track it is essential to know the W-value and its dependence on ion energy. In this way, for each ion of a given energy, an appropriate W-value can be selected at every collision point according to its current energy. Unfortunately, for carbon and sulphur ions data on W-values is currently scarce with values not being well determined. Therefore, to ensure a range of possible W-value outcomes are covered, it was decided for this work, as shown in Fig. 4, to use two sets of values. The first set are measured values taken from [START_REF] Snowden-Ifft | [END_REF] and the second are predicted values derived by A. Hitachi [11]. These are the only data available as known to the authors and come with no clear information on uncertainties. As seen in Fig. 4 there is only partial agreement between the data and simulation, and in both cases values of W at energy below 20 keV are not determined. To allow for this and to provide a range of values at all energies of interest, it was decided to use both sets in this work and, as shown in Fig. 4, to use a fit to the points in both cases of the form:

f `xe f `xe a b
x c where a, b and c are free parameters, including an extrapolation below 20 keV. Values of the fit parameters used are shown in Table 2.

Using the W-value fit of Fig. 4, the ionisation energy loss and ionisation charge spatial distributions, and hence head-tail effects, were studied for two different geometric configurations relative to the ion interactions points: (i) along the track, and (ii) projected onto an axis aligned with the initial direction of the ion. These geometries allow exploration of the influence of readout on head-tail sensitivity. In particular, the first configuration represents the ideal case of perfect 3D track reconstruction. The second configuration is relevant to a simpler 1D type readout. Fig. 5 shows an illustrative plot with projected lines for readout geometries for an example 100 keV sulphur ion track. On this plot the small circles along the track represent interaction points of energy loss with the size being proportional to the magnitude of the energy loss. It is important to consider that different range straggling will occur for ions of the same energy. Hence, for comparison purposes in this plot and elsewhere, the position of each interaction point has been normalised to the projected range of the ion. This is discussed further in Sec. 2.2.

Ion Ranges

Based on the processes described above a study was performed first to compare ion ranges, before diffusion, as a function of energy and NIPs. Fig. 6 shows results for sulphur and carbon in CS 2 vs. energy in comparison with data from [START_REF] Snowden-Ifft | [END_REF]. Here, the projected range was used, calculated using the coordinates of the first and last point of the ion track. It can be seen that there is a very good agreement for sulphur. For carbon ions however, agreement is less compelling, the results from this work yielding less than half the values in [START_REF] Snowden-Ifft | [END_REF]. This discrepancy, which is the subject of further experimental study, although significant does not effect the general conclusions of this work, since sulphur and fluorine have the dominant response for dark matter scattering. For the conversion of the energy loss into ionization charge, the number of electrons attached, with no loss, to CS 2 molecules creating CS 2 was summed for each carbon and sulphur ion energy. Fig. 7 shows, as a function of NIPs, the resulting range of nuclear recoils together with the range of alpha particles and electrons calculated with TRIM and ESTAR [START_REF]ESTAR: stopping power and range tables for electrons[END_REF], respectively. The NIP values were calculated here for the latter species using a W-value of 19 eV [START_REF][END_REF]. It can be seen from Fig. 7, comparing electron recoils with nuclear recoils yielding the same number of NIPs, that the difference in range is typically one order of magnitude greater for electrons than for nuclear recoils. This illustrates the expected great electron (gamma) rejection power of a TPC operating with 40 Torr CS 2 , as has been experimentally confirmed [15]. However, it also shows that the range of nuclear recoils and alpha particles of similar low energy are close to each other, the difference being not greater than 1-2 millimeters for tracks a few millimetres long. Note that for dark matter searches with a gaseous directional detector the internal background from alpha particles is in fact easily rejected by side veto regions and because their energy is generally far above the region of interest (I MeV). Ensuring high radiopurity of materials is sufficient to deal with potential events that may result from alpha particles that lose most energy within detector materials and then leak out into the sensitive volume with energy H 100 keV.

To check further the consistency of the SRIM procedure used here, and developed for the head-tail analysis below, some comparisons was made with existing measurements. In particular, the ranges of 103 keV 206 Pb recoils in several pure and binary gases at STP were calculated and compared with data from [16]. The results are shown in Table 1. As can be seen in most cases there is less than a 10% difference. Similar differences between measurements and simulation were observed for long range 5-6 MeV alpha particles from Rn and Po decays in 40 Torr CS 2 [17].

Head-Tail

Along with the event-by-event information on ion energy loss, TRIM also delivers histograms accumulated on-line and shows distributions of average values of several quantities as a function of target depth. With the default initial TRIM settings in place the target depth is defined by projection onto the axis of the initial direction of ion motion. One of the histograms shows distributions of the average ionization energy loss by the primary ion recoils and secondary recoils separately as a function of target depth. An example of such distributions is shown in Fig. 8 for 1000 sulphur ions with energy of 100 keV. Here one can see the ionization energy loss by the recoil ion, by the secondary recoils produced, and by their sum. The latter is compared in Fig. 8 with results derived from analysis of the TRIM output file using the event by event analysis where the position of those track charge depositions have been scaled to be appropriate for the target depth direction. Reassuringly, both show the same form of ionization energy loss vs. target depth, indicating here that more energy is lost at the beginning of the track (tail) than at the end (head).

However, whilst such averaged SRIM plots are applicable for use in, for instance, ion damage calculations, caution is needed in interpreting or using them directly for studying head-tail discrimination in a dark matter detector. This is because the distribution is an average provided as a function of target depth regardless of the ion range that, due to straggling, varies from ion to ion. For the experiments of interest here, recoils are detected on an event-by-event basis with ionization charge distribution created along the track and then projected onto an axis. For this situation, as pursued here, use of the track generation output file alone is more appropriate.

Based on this, Figs 9-11 show the distributions of total, nuclear and electronic energy loss for sulphur, carbon and fluorine ions along the track and projected onto the initial direction of the ion motion and normalised to its projected range R. Calculations were performed for ion energies up to 400 keV. As expected from the stopping power shape for both fluorine and carbon ions with high energies, where the electronic energy loss dominates, the total energy loss distribution has mainly a clear negative slope. This indicates greater energy loss at the tail of the track than at the head. However, due to ion straggling at the very end of the track, the energy loss there accumulates in a small volume that, when projected, creates a clearly visible sharp rise at the end of the energy loss distribution. For lower ion energies the slope at the beginning of distributions gets smaller and for energies below 200 keV and 50 keV for fluorine and carbon respectively it becomes flat. Only the sharp end remains, indicating now that more energy is lost at the head of the track than at the tail. For all studied energies of sulphur ions the energy loss distribution remains largely flat but, like the others, ends with a sharp rise. Clearly this distribution for sulphur differs from the one given in Fig. 8. It shows an opposite effect, because now the influence of straggling is properly included and this can dominate over the stopping power characteristic at the head of the track. Next, converting from energy loss to ionisation (NIPs) charges, the distributions of ionization charge were studied for sulphur and carbon ions in CS 2 . This was performed for the geometric scenario when ionization is projected onto the axis of the initial ion motion. The charge diffusion in gas was accounted for by redistributing the NIPs in space according to the thermal diffusion formula:

σ p 2k B T e p L E
where k B is the Boltzmann constant, T the ambient temperature, e the electron charge, E the drift electric field and L the drift distance. Note that the thermal behavior of diffusion in carbon disulphide negative ion gas has been well determined and measured [START_REF] Ohnuki | [END_REF]. After numerical simplifications the diffusion formula used was of form: σ mm¥ 0.72 q L m¥ E keV cm 1 ¥ for three different drift lengths L: 0.0, 0.3 and 0.5 m and electric field E = 0.58 kVcm 1 as used in the DRIFT IIb detector [15].

The head-tail effect information was quantified on an event-by-event basis using the same technique as formulated for the DRIFT II data analysis described in [7], through the quantity Ratio

î 50%R 0 N IP s î 100%R 50%R N IP s
vs. NIPs. The charge being simply integrated along the two halves of the projected range R. Fig. 12 shows results of the head-tail effect calculated in this way separately for sulphur and carbon tracks. Results are given as a function of NIPs for the three different W-value sets described earlier, those in Fig. 4 plus the electron values. To illustrate the impact of diffusion, for instance as might be expected in an experimental run of DRIFT, results are given for tracks drifted for a random distance up to a maximum of 0, 30 and 50 cm. Points with Ratio values above 1 are for ion tracks with greater tail than head and those below for ions with the reverse. Errors on the plots are dominated by systematic errors, including the assumed W-value, discussed below. Random errors on the points arising from the simulations are less than 1 percent (within the size of the symbols). However, note that in a practical detector the degree to which a value different from 1 can be observed as a head-tail asymmetry will be determined by several specific additional factors, specific to the detector and so outside the remit of this work. In particular, the readout position resolution, energy resolution and background. Nevertheless, we estimate that deviations from 1 of greater than 5% should reasonably be distinguishable.

Examining these plots we can note several features. Firstly, it can be seen that in general there is evidence for head-tail asymmetry in most cases and that the magnitude of this (that is the deviation from Ratio = 1.0) increases with increasing energy (NIPs), and increases in such a way as to give ion tracks with greater tail than head. Whilst this work is not aimed at simulating directionality in a specific detector design, it is interesting to note that the recent measurements by DRIFT II using directed neutrons also suggest a small positive effect with greater head than tail [7]. Secondly, as expected, the deviations from a ratio of 1.0 are generally greater when less diffusion is assumed (lower plots in Fig. 12). This reflects the greater smearing of charge along the tracks when diffusion is greater. Thirdly, there is a large dependence on the assumed W-value, particularly for sulfur. For this higher mass nucleus use of the Hitachi electronic W values results in a very strong head-tail asymmetry, with greater tail than head. However, this highlights again the need for improved measurements of W, which provides here the dominant uncertainty.

Examining the results more closely we see that the asymmetry is found to be most significant for the carbon recoils (right column of Fig. 12). For carbon although the head-tail asymmetry yields a ratio greater than 1.0 for most energies, the value does reach unity in some scenarios around 1000 NIPs and in the case of zero diffusion for the Hitachi W-values dips below 1.0, suggesting a reversal of the head-tail asymmetry. This phenomena of reduced asymmetry at lower energy and possible reversal of asymmetry appears more pronounced for the case of sulphur. Here, for the W values assumed and extrapolated from the Snowden-Ifft data, the ratio appears to reach 1.0 typically near 5000 NIPs but then rises again below 3000 NIPs. The Hitachi W values result in a steeper decrease in asymmetry with decreasing energy, sufficient to cause a reversal of asymmetry below 3000 NIPs. On the other hand for the electronic values the asymmetry ratio remains positive and large at all energies.

Focussing further now on the sulphur recoils, sulphur being the species of most interest for dark matter searches here (left column), although the asymmetry is greater than 10% for the electronic case (greater than 20% for zero drift distance) the deviation from ratio of 1.0 is less when the Hitachi and Snowden-Ifft W-values are used and less than that for carbon recoils. At high NIP values, for instance above about 6000 NIPs for the Hitachi points and 10,000 for the Snowden-Ifft points, the ratio does deviate by more than 10%. However, for the NIPs range of 1000-5000, the main region of interest for dark matter searches with sulphur recoils, for drift distances of 30 and 50 cm, it can be seen that although the plots do show deviations from a ratio of 1.0, these deviations are less than 5%. Somewhat greater deviations are seen in this region when no drift is assumed (bottom left plot).

Whilst this level of deviation is smaller than the level estimated as observable in practice, the possibility that it could be observed, with sufficient statistics in a real detector, can not be ruled out. Nevertheless, particularly if the Snowden-Ifft W-values are assumed, the results suggest that detection of a head-tail asymmetry for low energy sulphur recoils (below 200 keV), would be challenging. However, this latter conclusion appears somewhat at odds with the recent measurements of headtail asymmetry by DRIFT [7] where an observable effect is found in the relevant energy range for sulphur recoils. The quantity of test neutron events gathered in that experiment may mean that the effective sensitivity to the ratio is better than the 5% assumed in this work. Alternatively, the result may be indicating that there is a significant component of electronic stopping power, sufficient to boost the ratio above 1.0. These issues can only be resolved by further study.

Conclusions

In this work preliminary results of the head-tail effect of carbon and sulphur ions in low pressure CS 2 gas, such as used in the DRIFT directional dark matter detectors, have been presented. Results were obtained using the SRIM/TRIM package with a purpose-written Monte-Carlo using W-values available from both existing data and simulations. The shape of the ionisation charge distribution along an ion track, without effect from ion range straggling, is found to be similar to that previously presented by Hitachi. It shows that for ion energies down to 10 keV more ionization charge is created at the tail than at the head, irrespective of the ion type. However, when full account is taken of straggling at the end of tracks, diffusion, particularly in relation to the geometry of the readout, and the uncertainty in the W-values, the situation becomes complex.

In summary, we note some main conclusions. Firstly, for the scenarios explored here, some form of head-tail asymmetry is present at most energies explored. This is encouraging news for directional dark matter detectors. The asymmetry occurs usually with greater tail than head, as appears to have been seen in the limited experimental data available. However, for the assumed Hitachi W-values and lower assumed diffusion, there is evidence of a reversal of asymmetry at low energy, particularly for sulfur recoils. For the assumed W-values from Snowden-Ifft the asymmetry is less clear for sulfur, particularly in the region 2000-5000 NIPs. This likely reflects uncertainties in the W-values, for instance a very high asymmetry is predicted for sulfur using a pure electronic W-value.

Finally, we note that the nature of the sense of the asymmetry (whether head or tail is dominant), can not be predicted by simple simulations of the average ionisation loss, but rather is critically influenced by details of the straggling, readout used, diffusion and W-values at low energy. More experimental effort, underway by the DRIFT collaboration and others, is needed to untangle these issues further.

Acknowledgements

Authors would like to thank prof. A. Hitachi for valuable comments and sharing his calculation results of the W-values for Sulfur and Carbon ions in CS 2 and prof. J.F. 15 , from Snowden-Ifft ( ) [START_REF] Snowden-Ifft | [END_REF] and electronic (j). Points with ratio above 1 indicate that more ionization charge is created at the beginning (tail) while those below that level, at the end (head) of the track. See the text for an explanation of the ion energy sampling used.

Fig. 1 .

 1 Fig.1. Calculated with SRIM ionization energy loss dE/dx as a function of ion range in the medium. Energy losses in electronic (solid line) and nuclear (dashed line) channels are shown for Sulfur(!) and Carbon(") ions in 40 Torr CS 2 and Fluorine(#) in 100 Torr CF 4 . Ion energy is marked up to 500 keV with an interval of 50 keV.

12

 12

Figure 1 :Fig. 2 .

 12 Figure 1: Ionization energy loss dE/dx as a function of ion range in the medium, calculated with SRIM . Energy losses in electronic (solid line) and nuclear (dashed line) channels are shown for Sulfur ( ) and Carbon ( ) ions in 40 Torr CS 2 and Fluorine ( ) in 100 Torr CF 4 . Ion energy is marked up to 500 keV with an interval of 50 keV.

13

 13

Figure 2 :Fig. 3 . 14 Figure 3 :Fig. 4 .

 231434 Figure 2: Lateral (solid line) and longitudinal (dashed line) straggling/range as a function of range for Fluorine ions in CF 4 ( ), Sulfur ( ) and Carbon ( ) ions in CS 2 . Symbols mark the energy up to 500 keV with an interval of 50 keV.

Figure 4 :

 4 Figure 4: Carbon and Sulfur ion in CS 2 gas as a function of energy. Points marked with ( and ) show results from [10] while (j and Z) values calculated by A. Hitachi. Dashed lines are from a fit used with TRIM output file analysis.

Figure 5 :Fig. 6 .

 56 Figure 5: A single track of 100 keV Sulfur ion obtained from the TRIM output files. In this work the ionization energy spatial distribution and head-tail effect were studied along the ion track line (solid) and along the axis of the initial motion (dotted). The projected range is denoted by the dashed line. The size of the circles along the track is proportional to the magnitude of the energy loss in the collision. The markers TAIL and HEAD show the position of the beginning and the end of the track, respectively.

17

 17

Figure 6 :Fig. 7 .

 67 Figure 6: Sulfur and Carbon ion range in CS 2 as a function of the energy, calculated with SRIM ( and ) and taken from [10] (Z and j).

  18

Figure 7 :Fig. 8 .

 78 Figure7: Sulfur (Z), Carbon (j), alpha ( ) and electron ( ) range in CS 2 as a function of ionization charge (NIPs) created. NIPs for Sulfur and Carbon were calculated in this case using the W-values from Snowden-Ifft shown in Fig.4while for alpha particles and electrons a constant value of W = 19 eV was used. To calculate the range of electrons the ESTAR program was used for energies from 10 to 125 keV with steps 2.5, 5 and 10 up to 20, 60 and 100 keV, respectively. Results for alpha particle ranges are for energies from 10 to 150 with a step of 10 keV.

Figure 8 :Fig. 9 .

 89 Figure 8: Average ionization energy for 100 keV Sulfur ion ( ) and created in the cascade secondary Carbon and Sulfur nuclear recoils ( ) as a function of target depth, obtained from TRIM histogramming. The total energy loss ( ) is compared with scaled results from TRIM output file analysis (Z).

Figure 9 :Fig. 10 .

 910 Figure 9: Energy loss distributions for 400 keV Sulfur ions given along the track (top) and along the direction of the initial ion motion (bottom). Dashed, dotted and dash-dotted lines indicate the total, electronic and nuclear energy loss respectively. The vertical dashed lines indicate the ion energy at those points.

Figure 10 :Fig. 11 .

 1011 Fig.10. As in Fig.9for Carbon ion. Figure10: As in Fig.9for Carbon ions.

Figure 11 :Fig. 12 .Figure 12 :

 111212 Fig. 11. As in Fig.9for Fluorine ion. Figure11: As in Fig.9for Fluorine ions.

Table 1 :

 1 Ziegler for his help in SRIM calculations. P. Majewski would like to thank Ranges of 103 keV 206 Pb recoils in pure and binary gases at STP measured and calculated with SRIM. Values from measurement are taken from[16]. 1σ of the simulated range distributions is 10% of the tabulated mean values for all gases except Xe for which 1σ is 15%. Measurement uncertainty is 2µm.

	Gas Measurement [µm] SRIM [µm]
	Ar	79	73
	Xe	44	36
	CH 4	84	95
	C 2 H 4	58	61
	Air	83	80
	N	80	74

Table 2 :

 2 Values of the parameters from the fit to experimental and theoretical data of the W-values shown in Fig.4


Jonathan El

Laurent Gardes

Stéphane Girard

Armelle Guillou

Estimation of extreme quantiles from heavy and light tailed distributions

Keywords: AMS Subject Classifications: 62G05, 62G20, 62G30 Weibull tail-distributions, Pareto-type distributions, extreme quantile, maximum domain of attraction, asymptotic normality

, a new family of distributions is introduced, depending on two parameters τ and θ, which encompasses Pareto-type distributions as well as Weibull tail-distributions. Estimators for θ and extreme quantiles are also proposed, but they both depend on the unknown parameter τ , making them useless in practical situations. In this paper, we propose an estimator of τ which is independent of θ. Plugging our estimator of τ in the two previous ones allows us to estimate extreme quantiles from Pareto-type and Weibull tail-distributions in an unified way. The asymptotic distributions of our three new estimators are established and their efficiency is illustrated on a small simulation study and on a real data set.

Introduction

Let X 1 , . . . , X n be a sequence of independent and identically distributed random variables with a cumulative distribution function F and let X 1,n ≤ • • • ≤ X n,n denote the order statistics associated to this sample. The Gnedenko theorem [START_REF] Gnedenko | Sur la distribution limite du terme maximum d'une série aléatoire[END_REF] insures that for a large class of cumulative distribution functions, the maximum X n,n (after proper renormalization) converges in distribution to an extreme-value distribution with shape parameter γ. Depending on its sign, three possible maximum domains of attraction for F are possible: Fréchet (γ > 0), Gumbel (γ = 0) and Weibull (γ < 0). Since distributions in the Weibull maximum domain of attraction have a finite right tail, in most applications this maximum domain of attraction is not relevant. In this paper, we focus on the Fréchet and Gumbel maximum domains of attraction.

Distributions in the Fréchet maximum domain of attraction can be characterized through their survival function F = 1 -F as F (x) = x -1/γ L(x) where γ > 0 and L is a slowly varying function at infinity i.e. L(λx)/L(x) → 1 as x → ∞ for all λ ≥ 1. F is said to be regularly varying at infinity with index -1/γ. This property is denoted by F ∈ R -1/γ , (see [START_REF] Bingham | Regular Variation[END_REF] for more details on regular variations theory) and F is called a Pareto-type distribution. To make inference on the distribution tail, most approaches consist of using the k n upper order statistics X n-kn+1,n ≤ • • • ≤ X n,n since the tail information is only contained in the extreme upper part of the sample. Here, (k n ) is an intermediate sequence of integers i.e. such that lim n→∞ k n = ∞ and lim n→∞ k n /n = 0.

A large part of the extreme-value literature is devoted to the estimation of the tail-index γ > 0, the most known estimator being the Hill estimator [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF]. It can be equivalently defined in terms of log-excesses log(X n-i+1,n /X n-kn+1,n ) or in terms of log-spacings log(X n-i+1,n /X n-i,n ):

H n (k n ) = 1 k n -1 kn-1 i=1 log X n-i+1,n X n-kn+1,n = 1 k n -1 kn-1 i=1 i log X n-i+1,n X n-i,n . (1)
At the opposite, there is no simple representation for distributions in the Gumbel maximum domain of attraction. However, an interesting subset of the Gumbel maximum domain of attraction is the Weibull tail-distributions family. It encompasses for instance Weibull, Gaussian and Gamma distributions. Let us recall that a cumulative distribution function F has a Weibull tail if F (x) = exp (-H(x)) , where H ← (t) := inf{x, H(x) ≥ t} ∈ R β . The function H ← is the so-called generalized inverse of H. The tail of such distributions is driven by the shape parameter β > 0 called the Weibull tail-coefficient. Many papers are dedicated to the estimation of the Weibull tail-coefficient. A first family of approaches [START_REF] Beirlant | Semiparametric lower bounds for tail index estimation[END_REF][START_REF] Beirlant | The mean residual life function at great age: Applications to tail estimation[END_REF][START_REF] Broniatowski | On the estimation of the Weibull tail coefficient[END_REF][START_REF] Dierckx | A new estimation method for Weibull-type tails based on the mean excess function[END_REF] is based on the log-excesses while a second one relies on the log-spacings [START_REF] Beirlant | Practical analysis of extreme values[END_REF][START_REF] Diebolt | Bias-reduced estimators of the Weibull tail-coefficient[END_REF][START_REF] Gardes | Comparison of Weibull tail-coefficient estimators[END_REF][START_REF] Gardes | Estimation of the Weibull tail-coefficient with linear combination of upper order statistics[END_REF][START_REF] Girard | A Hill type estimate of the Weibull tail-coefficient[END_REF][START_REF] Goegebeur | Generalized kernel estimators for the Weibulltail coefficient[END_REF][START_REF] Goegebeur | Goodness-of-fit testing for Weibull-type behavior[END_REF]. All these estimators are thus similar to the Hill statistic [START_REF] Artin | The gamma function[END_REF].

In order to understand the similarity between most estimators of the Weibull tail-coefficient and the Hill estimator, a new family of distributions has been proposed in [START_REF] Gardes | Weibull tail-distributions revisited: a new look at some tail estimators[END_REF]. These distributions depend on two parameters τ ∈ [0, 1] and θ > 0. More specifically, letting K x (y) = y 1 u x-1 du where x ∈ R, the considered family of survival distribution functions is given by:

(A 1 (τ, θ)) F (x) = exp(-K ← τ (log H(x))) for x ≥ x * > 0 with τ ∈ [0, 1].
Here, H is an increasing function such that H ← ∈ R θ where θ > 0. The parameter τ allows us to represent a large panel of distribution tails ranging from Weibull-type tails (in this case τ = 0 and θ coincides with the Weibull tail-coefficient β) to distributions belonging to the Fréchet maximum domain of attraction (in this case τ = 1 and θ corresponds to the tail index γ).

In [START_REF] Gardes | Weibull tail-distributions revisited: a new look at some tail estimators[END_REF], an estimator of θ based on the Hill statistic is also introduced:

θ n,τ (k n ) = H n (k n ) µ τ (log(n/k n )) , (2)
with, for all t > 0,

µ τ (t) = ∞ 0 (K τ (x + t) -K τ (t)) e -x dx,
and an estimator of the extreme quantile x pn = F ← (p n ) with p n → 0 as n → ∞ is derived:

x pn, b θn,τ = X n-kn+1,n exp θ n,τ (k n ) [K τ (log(1/p n )) -K τ (log(n/k n ))] . (3)
We refer to [START_REF] Ditlevsen | Distribution Arbitrariness in Structural Reliability, Structural Safety and Reliability[END_REF][START_REF] Smith | Estimating the upper tail of flood frequency distributions[END_REF][START_REF] Embrechts | Modelling extremal events[END_REF] for applications of extreme quantiles respectively in reliability, hydrology and finance. The asymptotic distributions of θ n,τ and x pn, b θn,τ are established in [START_REF] Gardes | Weibull tail-distributions revisited: a new look at some tail estimators[END_REF]. Let us highlight that estimators (2) and (3) are only of theoretical interest since, in practical situations τ is unknown, and therefore, they cannot be used.

This paper builds on the work of [START_REF] Gardes | Weibull tail-distributions revisited: a new look at some tail estimators[END_REF]. We propose an estimator τ n of τ , independent of θ. Replacing τ by τ n in (2) and (3) yields two new estimators which can be computed in practical situations. As a result, we are able to estimate extreme quantiles from Pareto-type and Weibull tail-distributions in an unified way. The asymptotic normality of our three new estimators is also established.

The paper is organized as follows. The estimators are defined in Section 2, their asymptotic properties are established in Section 3. The behavior of the extreme quantile estimator is illustrated on simulated data in Section 4 and on a real data set in Section 5. Proofs of the main results are presented in Section 6 while proofs of auxiliary results are postponed to the Appendix.

Definition of the estimators

Let us first describe the construction of an estimator of τ . Let (k n ) and

(k n ) with k n > k n be two intermediate sequences of integers such that θ n,τ (k n ) P -→ θ and θ n,τ (k n ) P -→ θ. It straightforwardly follows that θ n,τ (k n ) θ n,τ (k n ) = H n (k n ) H n (k n ) µ τ (log(n/k n )) µ τ (log(n/k n )) P -→ 1.
Introducing for all t > t > 0 the function defined by ψ(.; t; t ) : R → (-∞, exp(t -t )) and ψ(x; t, t ) = µ x (t)/µ x (t ), it follows that

H n (k n ) H n (k n ) P ∼ ψ(τ ; log(n/k n ), log(n/k n )). (4)
Moreover, it can be shown (see Lemma 3 in Section 6) that ψ(.;

log(n/k n ), log(n/k n )) is a bijection from R to (-∞, k n /k n ).
As a consequence, the following estimator of τ is considered:

τ n =    ψ -1 Hn(kn) Hn(k n ) ; log(n/k n ), log(n/k n ) if Hn(kn) Hn(k n ) < k n kn u if Hn(kn) Hn(k n ) ≥ k n kn , (5)
where u is the realization of a standard uniform distribution. In practice, only the first situation has to be considered, since Lemma 5 in Section 6 shows that, for n large enough,

H n (k n )/H n (k n ) is almost surely smaller than k n /k n . It is thus possible to plug τ n in (2) to obtain a new estimator of θ: θ n,b τn (k n ) = H n (k n ) µ b τn (log(n/k n )) . (6)
Similarly, replacing τ and θ n,τ by their estimates in (3) yields a new estimator of extreme quantiles:

x pn, b θ n, b τn = X n-kn+1,n exp θ n,b τn (k n ) [K b τn (log(1/p n )) -K b τn (log(n/k n ))] . (7)
The asymptotic behavior of the three new estimators τ n , θ n,b τn and

x pn, b θ n, b
τn is established in the next section.

Asymptotic properties

As a first result, we establish the consistency of τ n under the following assumption:

H ← (t) = t θ (t) = c t θ exp x 1 ε(u) u du (8)
with c a positive constant and ε(s) → 0 as s → ∞. Let us highlight that (8) amounts to supposing that H ← ∈ R θ and that the slowly varying function at infinity is normalised, see [START_REF] Bingham | Regular Variation[END_REF], page 15.

Proposition 1 Suppose that (A 1 (τ, θ)) holds with (.) a normalised slowly varying function. If

(k n ) and (k n ) are two intermediate sequences of integers such that k n /k n → 0, then τ n P -→ τ.
Let us note that the consistency of τ n is established for all θ > 0 and τ ∈ [0, 1] in an unified way.

In this sense, the asymptotic behavior of this estimator is more a consequence of the log-spacings property than a tail behavior (which can be exponential for τ = 0 as well as polynomial for τ = 1). Next, to establish the asymptotic normality of the three estimators ( 5), ( 6) and ( 7), a second-order condition on is necessary:

(A 2 (ρ)) There exist ρ < 0 and b(x) → 0 such that uniformly locally on λ ≥ λ 0 > 0 log (λx) (x) ∼ b(x)K ρ (λ), when x → ∞,
with |b| asymptotically decreasing.

It can be shown that necessarily |b| ∈ R ρ . The second order parameter ρ < 0 tunes the rate of convergence of (λx)/ (x) to 1. The closer is ρ to 0, the slower is the convergence. Condition (A 2 (ρ)) is the cornerstone in all the proofs of asymptotic normality for extreme value estimators. It is used in [START_REF] Beirlant | Tail index estimation and an exponential regression model[END_REF][START_REF] Häusler | On asymptotic normality of Hill's estimator for the exponent of regular variation[END_REF][START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF] to prove the asymptotic normality of several estimators of the extreme value index. Let log 2 (.) = log(log(.)). The first theorem establishes the asymptotic normality of τ n :

Theorem 1 Suppose that (A 1 (τ, θ)) and (A 2 (ρ)) hold. If (k n ) and (k n ) are two intermediate sequences of integers such that k n /k n → 0 and k n b(exp K τ (log n/k n )) → 0, (9)
k n (log 2 (n/k n ) -log 2 (n/k n )) → ∞, (10) log
(n/k n ) (log 2 (n/k n ) -log 2 (n/k n )) → ∞, (11)
then k n (log 2 (n/k n ) -log 2 (n/k n )) ( τ n -τ ) d -→ N (0, 1).
Condition ( 9) is standard in extreme-value theory. It imposes that the bias induced by the slowly varying function is asymptotically negligible. Condition [START_REF] Diebolt | Quasi-Conjugate Bayes Estimates for GPD Parameters and Application to Heavy Tails Modeling[END_REF] forces the speed of convergence of τ n to tend to infinity. Finally, [START_REF] Diebolt | Bias-reduced estimators of the Weibull tail-coefficient[END_REF] is of the same nature as [START_REF] Diebolt | Quasi-Conjugate Bayes Estimates for GPD Parameters and Application to Heavy Tails Modeling[END_REF]: It imposes some minimal spacing between the two sequences (k n ) and (k n ). Besides, if τ = 0, conditions ( 9) and ( 10) imply that xb(x) → 0 as x → ∞. As a consequence, if τ = 0 and ρ > -1, it is not possible to choose sequences (k n ) and (k n ) satisfying the above assumptions. In such a case, only the consistency of τ n can be guaranteed.

In our next result, it is established that θ n,b τn inherits from the asymptotic normality of τ n .

Theorem 2 Suppose the assumptions of Theorem 1 hold. If, moreover,

(log 2 (n/k n ) -log 2 (n/k n )) / log 2 (n/k n ) → 0, ( 12
)
k n (log 2 (n/k n ) -log 2 (n/k n )) / log 2 (n/k n ) → ∞, (13)
then √ k n (log 2 (n/k n ) -log 2 (n/k n )) log 2 (n/k n ) θ n,b τn (k n ) -θ d -→ N (0, θ 2 ).
It appears that the estimation of τ has a cost in terms of rates of convergence. Condition [START_REF] Diebolt | Bias-reduced extreme quantiles estimators of Weibull tail-distributions[END_REF] implies that θ n,b τn converges slower than θ n,τ , see Lemma 7 in Section 6. As previously, ( 13) forces the speed of convergence of θ n,b τn (k n ) to tend to infinity. Note that this condition implies [START_REF] Diebolt | Quasi-Conjugate Bayes Estimates for GPD Parameters and Application to Heavy Tails Modeling[END_REF] in Theorem 1. Similarly as for Theorem 1, it can be shown that, if τ = 0, ( 9) and ( 13) imply

x log(x) b(x) → 0 as x → ∞. Thus, again no sequences (k n ) and (k n ) exist in case τ = 0 and ρ > -1. Now, if τ ∈ (0, 1] or if τ = 0 and ρ < -1, a possible choice for the two intermediate

sequences is log(k n ) = aK τ (log(n)) and log(k n ) = a K τ (log(n))
with the following restrictions on (a, a ) ∈ R:

   0 < a < a < 2ρ/(2ρ -1) if τ = 1 0 < a < a < -2ρ if 0 < τ < 1 2 < a < a < -2ρ if τ = 0. ( 14
)
Finally, in the case where τ = 0 and ρ = -1, the existence of sequences (k n ) and (k n ) depends on the underlying distribution.

The last result is dedicated to the asymptotic distribution of the extreme quantile estimator.

Theorem 3 Suppose the assumptions of Theorem 2 hold. If, moreover,

k n (log 2 (n/k n ) -log 2 (n/k n )) /log 2 (1/p n ) → ∞, (15)
(log(n/k n )) 1-τ [K τ (log(1/p n )) -K τ (log(n/k n ))] → ∞, (16)
log 2 (n/k n )[K τ (log(1/p n )) -K τ (log(n/k n ))] log(1/pn) log(n/kn) log(u)u τ -1 du → 0, (17)
then √ k n (log 2 (n/k n ) -log 2 (n/k n )) log(1/pn) log(n/kn) log(u)u τ -1 du x pn, b θ n, b τn x pn -1 d -→ N (0, θ 2 ).
A sufficient condition to verify ( 16) and ( 17)

is log 2 (1/p n )/ log 2 (n/k n ) → ∞.
This imposes an upper bound on the order p n of the extreme quantile. At the opposite, condition (15) provides a lower bound on p n . A possible choice for the order p n of the extreme quantile is given by log

2 (1/p n ) = [log 2 (n)] α for all α > 1. The finite sample performances of x pn, b θ n, b
τn are illustrated in the next section.

A small simulation study

In this section, our extreme quantile estimator is compared to the Moment estimator of Dekkers et al. [START_REF] Dekkers | A moment estimator for the index of an extreme-value distribution[END_REF] and the Peaks Over Threshold (POT), see for instance [START_REF] Davison | Models for exceedances over high thresholds[END_REF], on simulated data. The POT approach relies on the approximation of the excesses distribution, over a high threshold, by a Generalized Pareto Distribution (GPD). Among the numerous methods available for estimating the GPD parameters, we focus on the moments method which yields the best results in our simulations. Let us emphasize that the Moment and POT estimators are designed to work in any domain of attraction. Let us recall that the extreme quantile estimator [START_REF] Broniatowski | On the estimation of the Weibull tail coefficient[END_REF] requires the numerical inversion of the function ψ. This computation is achieved thanks to a dichotomy procedure since Lemma 3(i) ensures that ψ is increasing.

The comparison is achieved on twelve different distributions:

• five Pareto-type distributions (τ = 1): the absolute value of the standard Cauchy distribution (θ = 1 and ρ = -2), standard Pareto distribution (θ = 2 and ρ = -∞), the absolute value of the Student distribution with two degrees of freedom (θ = 1/2 and ρ = -2) and two Burr distributions (θ = 1/2 and ρ ∈ {-1, -1/2}) with quantile function F ← (x) = (x ρ -1) -θρ ,

x ∈ (0, 1),

• five Weibull tail-distributions (τ = 0): the absolute value of the standard Gaussian distribution (θ = 1/2 and ρ = -1), a Weibull distribution with 2 as shape parameter (θ = 2 and ρ = -∞), a Gamma distribution with 2 as shape parameter (θ = 1 and ρ = -1) and two distributions with quantile function

F ← (x) = (-log x) θ (1+ (ρ+θ)(-log x) ρ ), x ∈ (0, 1) with θ = 1/2 and ρ ∈ {-1/2, -1/4},
• two log-Weibull tail distributions (τ = 1/2), see [START_REF] Gardes | Weibull tail-distributions revisited: a new look at some tail estimators[END_REF], paragraph 2.2: the standard lognormal distribution (θ = √ 2/2 and ρ = 0) and the distribution with quantile function

F ← (x) = exp{ √ 2[(-log x) 1/2 -1]}, x ∈ (0, 1) for which θ = √ 2/2 and ρ = -∞. These distributions represent various situations (τ ∈ {0, 1/2, 1}, θ ∈ {1/2, √ 2/2, 1, 2} and ρ ∈ {-∞, -2, -1, -1/2, -1/4, 0}
) in which the asymptotic normality of our estimator is not always established. In the following, we take p n = 10 -3 and simulate N = 100 samples (X n,j ) j=1,...,N of size n = 500. For each sample X n,j , the estimator x pn, b θ n, b τn is computed for k n = 3, . . . , 500 and k n = ck n with c = 0.1. This value has been chosen on the basis of intensive Monte-Carlo simulations. With such a choice, our estimator and the Moment and POT estimators depend on only one intermediate sequence of integers (k n ). For all the considered distributions, the meansquared errors associated to the estimators are computed as functions of k n and are reported on Figures 123. It appears that our estimator outperforms the Moment and POT estimators for almost all the values of k n . Also for k n ≥ 50, the mean-squared errors associated to our estimator are almost constant as a function of k n , whatever the distribution is.

A real data set

The performance of our estimator [START_REF] Broniatowski | On the estimation of the Weibull tail coefficient[END_REF] is illustrated through the analysis of extreme events on the Nidd river data set, which is common in extreme values studies, see for instance [START_REF] Davison | Models for exceedances over high thresholds[END_REF]. It consists of 154 exceedances of the levels 65m 3 s -1 by the river Nidd (Yorkshire, England) during the period 1934-1969 (35 years). There is no general agreement on a maximum domain of attraction for this data set. In [START_REF] Diebolt | Quasi-Conjugate Bayes Estimates for GPD Parameters and Application to Heavy Tails Modeling[END_REF], a Fréchet maximum domain of attraction is assumed and heavy taileddistributions are considered as a possible model for such data. However, according to [START_REF] Hosking | Parameter and quantile estimation for the generalized Pareto distribution[END_REF], the Nidd data may reasonably be assumed to come from a distribution in the Gumbel maximum domain of attraction. This result was in accordance with [START_REF] Diebolt | Bias-reduced extreme quantiles estimators of Weibull tail-distributions[END_REF] where it has been shown that a Weibull tail-distribution could be considered for such a data. The estimation of τ is therefore of great interest. The estimated values of τ and θ are depicted on Figure 4 as functions of k n (recall that k n = 0.1k n ). It appears that the estimators become stable for k n ≥ 80 with τ n 1 and θ n,b τn (k n ) 0.3. These results indicate that the data may be assumed to come from a distribution in the Fréchet maximum domain of attraction. They are in accordance with the ones obtained by Bayesian methods [START_REF] Diebolt | Quasi-Conjugate Bayes Estimates for GPD Parameters and Application to Heavy Tails Modeling[END_REF], Figure 7 where the tail index is also estimated at 0.3.

The standard quantity of interest in environmental studies is the N -year return level, defined as the level which is exceeded on average once in N years. Here, we focus on the estimation of the 50and 100-year return levels. In Figure 4, our extreme quantile estimator is compared to the Moment and POT estimators, plotting the associated N -year return level as a function of k n for N = 50 and N = 100. It appears that our estimator and the Moment estimator yield similar curves. The POT sample path is different, but for all methods, choosing k n 60, we obtain an estimation of the 50-year return level which belongs approximately to the interval [340m 3 s -1 , 375m 3 s -1 ], and an estimation of the 100-year return in [400m 3 s -1 , 470m 3 s -1 ]. Again, these results are in accordance with the the credibility intervals obtained in [START_REF] Diebolt | Quasi-Conjugate Bayes Estimates for GPD Parameters and Application to Heavy Tails Modeling[END_REF], Table 1.

Proofs

We first give some preliminary lemmas. Their proofs are postponed to the appendix.

Preliminary lemmas

In the following, C is a compact subset such that [0, 1] ⊂ C ⊂ (-∞, 2). The first lemma is a standard result on the behavior at infinity of the Laplace transform.

Lemma 1 Let x ∈ C and h x ∈ C ∞ (R + ). Set i = min {j ∈ N/h (j) x (0) = 0}. If sup x∈C y≥0 h (i+1) x (y) < ∞, then lim t→∞ sup x∈C t i+1 h x (t) -h (i)
x (0) = 0, with h x (t) = +∞ 0 exp(-tu)h x (u)du is the Laplace transform of h x .

Lemma 1 is the key tool for establishing the uniform asymptotic expansions of µ x and ∂µ x /∂x given in Lemma 2:

Lemma 2 For all x ∈ C and t > 0, we have

(i) lim t→∞ sup x∈C µ x (t) t x-1 -1 = 0, (ii) lim t→∞ sup x∈C ∂ ∂x µ x (t) -log(t)µ x (t) t x-2 -1 = 0.
As a consequence of the above expansions, some important properties of ψ can be derived in the two next lemmas:

Lemma 3 For all t > t > 0 and x ∈ R, we have The following lemma ensures that the estimator τ n is well-defined. In [START_REF] Beirlant | Practical analysis of extreme values[END_REF], the situation where

(i) x → ψ(x; t,
H n (k n )/H n (k n ) ≥ k n /k n is eventually impossible. Lemma 5 Let (k n ) and (k n ) be two intermediate sequences of integers such that k n /k n → 0. If θ n,τ (k n )/ θ n,τ (k n ) P -→ 1 then P H n (k n ) H n (k n ) ≥ k n k n → 0 as n → ∞.
We now prove that θ n,τ (k n ) is a consistent estimator for θ when τ is known and under general assumptions.

Lemma 6 Suppose that (A 1 (τ, θ)) holds with a normalised slowly varying function (.).

If (k n ) is an intermediate sequence of integers, then θ n,τ (k n ) P -→ θ.
The asymptotic normality of θ n,τ and x pn, b θn,τ has already been established in the particular case where τ is known by [START_REF] Gardes | Weibull tail-distributions revisited: a new look at some tail estimators[END_REF]. Let us quote two results from this work.

Lemma 7 (Theorem 1, [START_REF] Gardes | Weibull tail-distributions revisited: a new look at some tail estimators[END_REF]). Suppose that (A 1 (τ, θ)) and

(A 2 (ρ)) hold. If (k n ) is an intermedi- ate sequence such that √ k n b(exp K τ (log(n/k n ))) → 0, then k n θ n,τ (k n ) θ -1 d -→ N (0, 1). ( 18
)
Lemma 8 (Theorem 2, [START_REF] Gardes | Weibull tail-distributions revisited: a new look at some tail estimators[END_REF]). Under the assumptions of Lemma 7 and if, moreover,

(log(n/k n )) 1-τ (K τ (log(1/p n )) -K τ (log(n/k n ))) → ∞, then, √ k n K τ (log(1/p n )) -K τ (log(n/k n )) log x pn, b θn,τ x pn d -→ N (0, θ 2 ).
The next lemma establishes that θ can be replaced by θ n,τ (k n ) in [START_REF] Gardes | Weibull tail-distributions revisited: a new look at some tail estimators[END_REF] without changing the asymptotic distribution.

Lemma 9 Suppose that (A 1 (τ, θ)) and (A 2 (ρ)) hold. Let (k n ) and (k n ) be two intermediate sequences of integers such that k n b(exp K τ (log n/k n )) → 0 and k n /k n → 0. We have

k n θ n,τ (k n ) θ n,τ (k n ) -1 d -→ N (0, 1).
The last lemma quantifies the effect of estimating τ in θ n,τ (k n ).

Lemma 10 Suppose that (A 1 (τ, θ)) and (A 2 (ρ)) hold. Let (k n ) and (k n ) be two intermediate sequences of integers such that

k n /k n → 0, k n (log 2 (n/k n ) -log 2 (n/k n )) / log 2 (n/k n ) → ∞, k n b(exp K τ (log n/k n )) → 0 and log(n/k n ) (log 2 (n/k n ) -log 2 (n/k n )) → ∞.
We have

√ k n (log 2 (n/k n ) -log 2 (n/k n )) log 2 (n/k n ) θ n,b τn (k n ) θ n,τ (k n ) -1 d -→ N (0, 1).

Proofs of the main results

Proof of Proposition 1 -For all ε > 0, let us write

P (| τ n -τ | > ε) = P ( τ n > τ + ε) + P ( τ n < τ -ε)
and consider the two terms separately. Clearly, one has

P( τ n > τ + ε) = P { τ n > τ + ε} H n (k n ) H n (k n ) ≥ k n k n P H n (k n ) H n (k n ) ≥ k n k n + P { τ n > τ + ε} ∩ H n (k n ) H n (k n ) < k n k n ≤ P H n (k n ) H n (k n ) ≥ k n k n + P H n (k n ) H n (k n ) ≥ ψ(τ + ε; log(n/k n ), log(n/k n )) . ( 19
)
Focusing on the first probability of ( 19), note that, θ n,τ (k n ) and θ n,τ (k n ) are consistent estimators for θ in probability by Lemma 6. Thus θ n,τ (k n )/ θ n,τ (k n ) P -→ 1 and Lemma 5 entail

P H n (k n ) H n (k n ) ≥ k n k n → 0.
The second probability in [START_REF] Girard | A Hill type estimate of the Weibull tail-coefficient[END_REF] tends to 0 by (4) and Lemma 3(i). Similarly, we can prove that P ( τ n < τ -ε) → 0 which achieves the proof.

Proof of Theorem 1 -Letting v n = √ k n (log 2 (n/k n ) -log 2 (n/k n ))
, our goal is to prove that P (v n ( τ n -τ ) ≤ s) → Φ(s), for all s ∈ R where Φ is the cumulative distribution function of the standard Gaussian distribution. To this aim, let us first remark that Lemma 9 yields θ n,τ (k n )/ θ n,τ (k n ) P -→ 1. Thus, introducing E n (s) = {v n ( τ n -τ ) ≤ s}, we have

P(E n (s)) = P E n (s) ∩ H n (k n ) H n (k n ) < k n k n + P E n (s) H n (k n ) H n (k n ) ≥ k n k n P H n (k n ) H n (k n ) ≥ k n k n =: T (1) n (s) + o(1),
by Lemma 5. From the definition of τ n given in [START_REF] Beirlant | Practical analysis of extreme values[END_REF] and recalling that, from Lemma 3(i), ψ(.; log(n/k n ), log(n/k n )) is an increasing function, we obtain

T (1) n (s) = P H n (k n ) H n (k n ) ≤ ψ(τ + s/v n ; log(n/k n ), log(n/k n )) ∩ H n (k n ) H n (k n ) < k n k n = P H n (k n ) H n (k n ) ≤ min ψ(τ + s/v n ; log(n/k n ), log(n/k n )), k n k n .
Lemma 3(ii) thus yields

T (1) n (s) = P H n (k n ) H n (k n ) ≤ ψ(τ + s/v n ; log(n/k n ), log(n/k n )) = P H n (k n ) H n (k n ) ≤ µ τ +s/vn (log(n/k n )) µ τ +s/vn (log(n/k n )) ,
and, from (2) and the fact that

ζ n := k n θ n,τ (k n ) θ n,τ (k n ) -1 d -→ N (0, 1),
(see Lemma 9), we have

T (1) n (s) = P ζ n ≤ k n µ τ (log(n/k n )) µ τ (log(n/k n )) µ τ +s/vn (log(n/k n )) µ τ +s/vn (log(n/k n )) -1 = P ζ n ≤ k n µ τ (log(n/k n )) µ τ (log(n/k n )) [ψ(τ + s/v n ; log(n/k n ), log(n/k n )) -ψ(τ ; log(n/k n ), log(n/k n ))] .
A first order Taylor expansion leads to

T (1) n (s) = P ζ n ≤ s (log 2 (n/k n ) -log 2 (n/k n )) µ τ (log(n/k n )) µ τ (log(n/k n )) ∂ ∂x ψ(τ 0 ; log(n/k n ), log(n/k n )) ,
where τ 0 = τ + sη 0 /v n with η 0 ∈ (0, 1). Since τ 0 → τ , for n large enough, τ 0 < 2 and Lemma 4 entails

T (1) n (s) = P ζ n ≤ s µ τ0 (log(n/k n )) µ τ0 (log(n/k n )) µ τ (log(n/k n )) µ τ (log(n/k n )) 1 + O 1 log(n/k n ) (log 2 (n/k n ) -log 2 (n/k n )) = P ζ n ≤ s µ τ0 (log(n/k n )) µ τ0 (log(n/k n )) µ τ (log(n/k n )) µ τ (log(n/k n )) (1 + o(1)) .
Besides, Lemma 2(i) entails that

s µ τ0 (log(n/k n )) µ τ0 (log(n/k n )) µ τ (log(n/k n )) µ τ (log(n/k n )) = s log(n/k n ) log(n/k n ) τ0-τ (1 + o(1)) = s exp sη 0 √ k n (1 + o(1)) -→ n→∞ s,
and thus

T (1) n (s) = P(ζ n ≤ s(1 + o(1))) ≤ Φ(s) + sup x∈R |Φ(x) -P(ζ n ≤ x)| → Φ(s)
by [START_REF] Embrechts | Modelling extremal events[END_REF], page 552. This achieves the proof of Theorem 1.

Proof of Theorem 2 -On the first hand, Lemma 7 shows that

θ n,τ (k n ) -θ = ζ n θ √ k n where ζ n d -→ N (0, 1),
and on the other hand, Lemma 10 entails

θ n,b τn (k n ) θ n,τ (k n ) -1 = δ n log 2 (n/k n ) √ k n (log 2 (n/k n ) -log 2 (n/k n )) where δ n d -→ N (0, 1). ( 20
)
Collecting these results, it follows that

θ n,b τn (k n ) -θ = θ n,τ (k n ) θ n,b τn (k n ) θ n,τ (k n ) -1 + θ n,τ (k n ) -θ = θ n,τ (k n ) δ n log 2 (n/k n ) √ k n (log 2 (n/k n ) -log 2 (n/k n )) + ζ n θ √ k n .
Consequently, one immediately has

√ k n (log 2 (n/k n ) -log 2 (n/k n )) log 2 (n/k n ) θ n,b τn (k n ) -θ = θ n,τ (k n )δ n + ζ n θ log 2 (n/k n ) -log 2 (n/k n ) log 2 (n/k n )
and Theorem 2 follows. x pn =:

T (2) n + T (3) n .
By (3), we have

T (2) n = θ n,b τn (k n )[K b τn (log(1/p n )) -K b τn (log(n/k n ))] -θ n,τ (k n )[K τ (log(1/p n )) -K τ (log(n/k n ))] = [K τ (log(1/p n )) -K τ (log(n/k n ))] θ n,b τn (k n ) -θ n,τ (k n ) + θ n,b τn (k n ) [(K b τn (log(1/p n )) -K b τn (log(n/k n ))) -(K τ (log(1/p n )) -K τ (log(n/k n )))] =: T (2,1) n + T (2,2) n .
Focusing on the first term, (20) leads to

T (2,1) n = log 2 (n/k n )[K τ (log(1/p n )) -K τ (log(n/k n ))] √ k n (log 2 (n/k n ) -log 2 (n/k n )) θ n,τ (k n )δ n which implies that √ k n (log 2 (n/k n ) -log 2 (n/k n )) log(1/pn) log(n/kn) u τ -1 log(u)du T (2,1) n = δ n log 2 (n/k n )[K τ (log(1/p n )) -K τ (log(n/k n ))] log(1/pn) log(n/kn) u τ -1 log(u)du θ n,τ (k n ) = o P (1).
Let us now consider T

. Theorem 1 states that

τ n = τ + ξ n √ k n (log 2 (n/k n ) -log 2 (n/k n )) =: τ + ξ n σ n where ξ n d -→ N (0, 1).

Replacing in T

(2,2) n

, we obtain

T (2,2) n = θ n,b τn (k n )[(K τ +σnξn (log(1/p n )) -K τ +σnξn (log(n/k n ))) -(K τ (log(1/p n )) -K τ (log(n/k n )))].
By definition

K τ (log(1/p n )) -K τ (log(n/k n )) = log(1/pn) log(n/kn) u τ -1 du,
and therefore it immediately follows that

T (2,2) n = θ n,b τn (k n ) log(1/pn) log(n/kn) u τ +σnξn-1 -u τ -1 du = θ n,b τn (k n ) log(1/pn) log(n/kn) u τ -1 u σnξn -1 du.
Letting ϕ(x) := exp(x) -1 -x, we deduce that

T (2,2) n = θ n,b τn (k n )σ n ξ n log(1/pn) log(n/kn) u τ -1 log(u)du + θ n,b τn (k n ) log(1/pn) log(n/kn) u τ -1 ϕ(σ n ξ n log(u))du.
≤ cσ n ξ 2 n log(1/pn) log(n/kn) u τ -1 [log(u)] 2 du 2 log(1/pn) log(n/kn) u τ -1 log(u)du ≤ cξ 2 n log 2 (1/p n ) 2 √ k n (log 2 (n/k n ) -log 2 (n/k n )) = o P (1),
and, replacing in

T (2,2) n , we obtain √ k n (log 2 (n/k n ) -log 2 (n/k n )) log(1/pn) log(n/kn) log(u)u τ -1 du T (2,2) n = ξ n θ n,b τn (k n ) + o P (1) d -→ N (0, θ 2 ).
Finally, Lemma 8 states that

T (3) n = n [K τ (log(1/p n )) -K τ (log(n/k n ))] √ k n where n d -→ N (0, θ 2 ),
and thus, under our assumptions,

√ k n (log 2 (n/k n ) -log 2 (n/k n ))
log(1/pn) log(n/kn) u τ -1 log(u)du

T (3) n = o P (1).
Combining the above results and using the delta method, Theorem 3 follows.

Appendix: Proof of auxiliary results

Proof of Lemma 1 -A (i + 1) th order Taylor expansion leads to

h x (u) = u i i! h (i) x (0) + u i+1 (i + 1)! h (i+1)
x (ηu), with η ∈ (0, 1) and consequently,

h x (t) = ∞ 0 exp(-tu) u i i! h (i) x (0)du + ∞ 0 exp(-tu) u i+1 (i + 1)! h (i+1)
x (ηu)du =: T (1) x (t) + T (2) x (t).

It follows that

T (1) x (t) = h (i) x (0) i! ∞ 0 u i exp(-tu)du = h (i) x (0) t i+1 ,
and the change of variable v = tu yields

T (2) x (t) = 1 (i + 1)! 1 t i+2 ∞ 0 exp(-v)v i+1 h (i+1) x η v t dv.
Therefore, we have, uniformly in x ∈ C:

|T (2) x (t)| ≤ sup z∈C y≥0 h (i+1) z (y) 1 (i + 1)! 1 t i+2 ∞ 0 exp(-v)v i+1 dv = 1 t i+2 sup z∈C y≥0 h (i+1) z (y) = O 1 t i+2
and thus

h x (t) = h (i)
x (0)

t i+1 + O 1 t i+2 = 1 t i+1 h (i) x (0) + O 1 t ,
which achieves the proof.

Proof of Lemma 2 -(i) By definition, for all x ∈ R and t > 0, we have

µ x (t) = ∞ 0 (K x (u + t) -K x (t)) exp(-u)du = ∞ 0 u+t t y x-1 dy exp(-u)du.
Using Fubini's theorem, it follows

µ x (t) = ∞ t y x-1 ∞ y-t exp(-u)du dy = exp(t) ∞ t y x-1 exp(-y)dy, (21)
and the change of variable u = y/t -1 yields

µ x (t) = t x ∞ 0 exp(-tu)(u + 1) x-1 du = t x h x (t)
with h x (t) = (t + 1) x-1 . Applying Lemma 1 with i = 0 concludes the proof of (i).

(ii) From ( 21) we obtain, for x ∈ R and t > 0,

µ x (t) = exp(t) ∞ t y x-1 exp(-y)dy = exp(t)Γ(x, t), (22)
where Γ(x, t) is the upper incomplete gamma function. We thus have (see for instance [START_REF] Artin | The gamma function[END_REF]) The result is proved.

∂ ∂x µ x (t) = exp(t) ∞ t exp(-y) log(y)y x-1 dy, (23)
Proof of Lemma 4 -Recall that C is a compact subset such that [0, 1] ⊂ C ⊂ (-∞, 2
Proof of Lemma 5 -Let us remark that

P H n (k n ) H n (k n ) ≥ k n k n = P θ n,τ (k n ) θ n,τ (k n ) ≥ k n k n µ τ (log(n/k n )) µ τ (log(n/k n )) = P θ n,τ (k n ) θ n,τ (k n ) ≥ k n k n log(n/k n ) log(n/k n ) 1-τ (1 + o(1)) by Lemma 2(i). Moreover k n /k n → ∞ implies (log(n/k n )/ log(n/k n )) 1-τ ≥ 1 eventually while θ n,τ (k n )/ θ n,τ (k n ) P -→ 1
θ n,τ (k n ) = 1 µ τ (log(n/k n )) θ 1 k n -1 kn-1 i=1 [K τ (F i + E n-kn+1,n ) -K τ (E n-kn+1,n )] + 1 k n -1 kn-1 i=1 log (exp K τ (F i + E n-kn+1,n )) (exp K τ (E n-kn+1,n )) . (25)
Now, Lemma 2(ii) and Lemma 5 in Gardes et al. (2011) yield

1 µ τ (E n-kn+1,n ) 1 k n -1 kn-1 i=1 [K τ (F i + E n-kn+1,n ) -K τ (E n-kn+1,n )] P -→ 1. (26)
Moreover, since the slowly varying function (.) is assumed to be normalised, it follows that log

(exp(K τ (F i + E n-kn+1,n ))) (exp(K τ (E n-kn+1,n ))) = exp(Kτ (Fi+E n-kn +1,n )) exp(Kτ (E n-kn +1,n )) ε(u) du u , with ε(s) → 0 as s → ∞. Thus for all i = 1, . . . , k n -1, we have exp(Kτ (Fi+E n-kn +1,n )) exp(Kτ (E n-kn+1,n )) ε(u) du u ≤ sup u≥exp(Kτ (E n-kn+1,n )) |ε(u)| exp(Kτ (Fi+E n-kn +1,n )) exp(Kτ (E n-kn +1,n )) du u = o P (1) [K τ (F i + E n-kn+1,n ) -K τ (E n-kn+1,n )] .
Thus, (26) entails

1 µ τ (log(n/k n )) 1 k n -1 kn-1 i=1 log (exp K τ (F i + E n-kn+1,n )) (exp K τ (E n-kn+1,n )) = o P (1). (27)
Collecting ( 25)-( 27), Lemma 6 follows.

from Lemma 2(i). Moreover, under the event A n , using Lemma 2(i) and Theorem 1, we have

µ τ0 (log(n/k n )) µ b τn (log(n/k n )) = (log(n/k n )) τ0-b τn (1 + o P (1)) = (log(n/k n )) η(τ -b τn) (1 + o P (1)) = exp(η(τ -τ n ) log 2 (n/k n ))(1 + o P (1)) = exp - ηξ n log 2 (n/k n ) √ k n (log 2 (n/k n ) -log 2 (n/k n )) (1 + o P (1)) P -→ 1.
We thus obtain

w n θ n,b τn (k n ) θ n,τ (k n ) -1 = ξ n (1 + o P (1)),
and, replacing in T

n (s), it finally follows that

T (4) n (s) = P({ξ n ≤ s(1 + o P (1))} ∩ A n ) = P(ξ n ≤ s(1 + o P (1))) -P({ξ n ≤ s(1 + o P (1))} |A c n )P(A c n ).
Using the same arguments as in the proof of Theorem 1, it is easily seen that P(ξ n ≤ s(1+o P (1))) → Φ(s) and thus T

Proof of Theorem 3 -

 3 Let us consider the following expansion:

Now, there existsc 2 σ 2 n ξ 2 n

 2 c > 0 such that x < log(c) implies |ϕ(x)| < c 2 x 2 . As a consequence, since σ n log 2 (1/p n ) → 0 and σ n log 2 (n/k n ) → 0, for n large enough, we have log(1/pn) log(n/kn) u τ -1 ϕ(σ n ξ n log(u))du ≤ log(1/pn) log(n/kn) u τ -1 |ϕ(σ n ξ n log(u))|du ≤ log(1/pn) log(n/kn) u τ -1 [log(u)] 2 du.Thus, log(1/pn) log(n/kn) u τ -1 ϕ(σ n ξ n log(u))du σ n log(1/pn) log(n/kn) u τ -1 log(u)du

  by assumption. The conclusion follows. Proof of Lemma 6 -According to Lemma 3 in Gardes et al. (2011), we have

( 4 )

 4 n (s) → Φ(s) since P(A c n ) → 0.Combining the above results the proof of Lemma 10 is achieved.

Figure 1 :

 1 Figure 1: Mean-squared errors as a function of k n associated to x pn, b θ n, b τn (dotted line) and to the Moment estimator of Dekkers et al. (bold line) and the POT estimator (solid line). They are computed on 100 samples of size 500. Upper left: absolute value of Cauchy distribution, upper right: absolute value of Student distribution, bottom left: Burr distribution with ρ = -1, bottom right: Burr distribution with ρ = -1/2.

Figure 2 :Figure 3 :

 23 Figure 2: Mean-squared errors as a function of k n associated to x pn, b θ n, b τn (dotted line) and to the Moment estimator of Dekkers et al. (bold line) and the POT estimator (solid line). They are computed on 100 samples of size 500. Upper left: Pareto distribution, upper right: Weibull distribution, bottom left: absolute value of Gaussian distribution, bottom right: Gamma distribution.

Figure 4 :

 4 Figure 4: Results obtained on the Nidd river data set. Top: estimation of τ (left) and θ (right) as functions of k n . Bottom: N -year return levels as functions of k n obtained with x pn, b θ n, b τn (dotted line) and the Moment estimator of Dekkers et al. (bold line) and the POT estimator (solid line) (left: N = 50, right: N = 100).

Acknowledgments

The authors are grateful to the referees for their helpful comments.

and the change of variable u = y/t -1 yields

exp(-tu) exp(-t)t x-1 (u + 1) x-1 log(t(u + 1))tdu = t x ∞ 0 exp(-tu)(u + 1) x-1 log(u + 1)du + log(t)t x ∞ 0 exp(-tu)(u + 1) x-1 du =: t x ∞ 0 exp(-tu)g x (u)du + log(t)µ x (t), with g x (u) := (u + 1) x-1 log(u + 1). Applying Lemma 1 with i = 1, the conclusion follows.

Proof of Lemma 3 -(i) First of all, remark that ψ(x; t, t ) > 0 for all x ∈ R and t > t > 0.

Besides, routine calculations show that

=: ψ(x; t, t

where, by ( 21) and [START_REF] Häusler | On asymptotic normality of Hill's estimator for the exponent of regular variation[END_REF],

Let us remark that Q x is an increasing function on (0, ∞) since

As a consequence t > t implies ∂/∂x(ψ(x; t, t )) = ψ(x; t, t )(Q x (t) -Q x (t )) > 0 and concludes the first part of the proof.

(ii) From ( 22), we have

Besides, considering the following inequalities,

and

and thus M x (t, t ) → 0 as x → ∞. The conclusion follows.

Proof of Lemma 9 -First of all, let us note that the two conditions

). As a consequence, we have

and therefore

It concludes the proof.

Proof of Lemma 10 -A first order Taylor expansion yields

, with τ 0 = τ n + η(τ -τ n ) and η ∈ (0, 1). Let

and recall that Φ is the cumulative distribution function of the standard Gaussian distribution. Our aim is to prove that P (F n (s)) → Φ(s), for all s ∈ R. Keeping in mind that C is a compact subset such that [0, 1] ⊂ C ⊂ (-∞, 2) and A n = {τ 0 ∈ C} ∩ { τ n ∈ C}, we have

where A c n is the complementary event of A n . Let us first consider T

n (s). Theorem 1 states that

Consequently, under the assumptions of Theorem 1,

n (s) → 0. Let us focus on T