F Montupet-Leblond

L Corso

M Payet

R Delaporte-Mathurin

E Bernard

Y Charles

J Mougenot

S Vartanian

E A Hodille

C Grisolia

Permeation and trapping of hydrogen in Eurofer97

Keywords: Permeation, Trapping, Eurofer97, Hydrogen, Thermal Desorption Spectrometry, MHIMS, surfaces

Diffusion and trapping of hydrogen isotopes in fusion materials need to be fully described in order to evaluate permeation and retention in fusion reactors walls and breeding blankets. Hydrogen gas permeation experiments have been conducted on Eurofer97 with pressures ranging from 10 1 to 10 5 Pa and temperatures between 473 K and 673 K, resulting in solubility

. Trapping parameters have been investigated using thermal desorption spectrometry of deuterium-loaded samples coupled with parametric optimization, leading to detrapping energies E dt,1 = 0.51 eV, E dt,2 = 1.27 eV, E dt,3 = 1.65 eV and densities Nt,1 = 6.01•10 25 m -3 , Nt,2 = 6.44•10 22 m -3 , Nt,3 = 3.88•10 23 m -3 . This parametric optimization is performed using a kinetic surface model: the contribution of this model is compared to the results given by solubility and recombination rate models.

Introduction

Tritium permeation and inventory in the walls of fusion reactors need to be assessed to comply with safety and environmental regulations. The foreseen design for DEMO plasma facing components and breeding blankets includes Eurofer as a structural material [START_REF] Kohyama | Low-activation ferritic and martensitic steels for fusion application[END_REF][START_REF] Federici | An overview of the EU breeding blanket design strategy as an integral part of the DEMO design effort[END_REF], which calls for investigation of the transport and trapping characteristics of this low-activation steel.

To that end, hydrogen permeation experiments were conducted on Eurofer97 samples in order to measure diffusivity, permeability and trapping parameters and evaluate solubility of hydrogen in this material. These results show that hydrogen transport parameters are affected by trapping in the 400 K -700 K range. As permeation measurements only give an integrated vision of trapping (in the occurence of multiple trapping sites, permeation is not able to discriminate between trapping sites), thermal desorption spectrometry (TDS) experiments were performed on deuterium-loaded Eurofer97 in order to give a more detailed understanding of trapping in this material. The experimental results were analyzed using a macroscopic diffusion/reaction model implemented in the MHIMS code [START_REF] Hodille | Macroscopic rate equation modeling of trapping/detrapping of hydrogen isotopes in tungsten materials[END_REF] including a kinetic surface model [START_REF] Hodille | Simulations of atomic deuterium exposure in self-damaged tungsten[END_REF]. The detrapping energies and densities of trapping sites, as well as parameters that render the kinetics of hydrogen on the surface, are obtained from these simulations. These results altogether are shown to be in agreement, giving a complete set of parameters describing hydrogen permeation and trapping in Eurofer97. The improvement brought to the model by the kinetic surface model is shown by comparing the results with a recombination rate model.

Gas permeation

In a gas permeation experiment, the sample (typically disk-shaped) is exposed to a stepwise increase of gas pressure on one side while the other side is under vacuum. This unbalance drives the particles to permeate through the sample, creating a pressure rise on the downstream side. This downstream pressure rise can be analyzed using the timelag method (see section 2.4) to yield permeation parameters such as diffusivity and permeability (see section 2.1). The work presented here aims at describing a new hydrogen permeation setup, which was used to measure the transport parameters of hydrogen in Eurofer97.

Theoretical aspects

Hydrogen transport parameters such as diffusivity D and permeability Φ can be evaluated based on the pressure rise in the case of diffusion-limited regime (see section 2.4). Hydrogen solubility K, also called Sieverts' constant, is then calculated as the ratio of these quantities using the defining relation

Φ(T ) = D(T ) • K(T ) (1)
where T is the sample temperature, given in K.

Hydrogen isotopes permeation is affected by trapping: once in the bulk, dissolved atoms are likely to remain trapped at specific locations, such as dislocations or vacancies, which slows down the migration of hydrogen. Therefore, permeation experiments do not yield D and K, which correspond respectively to interstitial diffusion and interstitial dissolution, but effective diffusivity D eff and effective solubility K eff , defined as follows [START_REF] Esteban | Hydrogen transport and trapping in Eurofer97[END_REF]:

D(T ) = D 0 exp ! - E D k B T " (2)
K(T ) = K 0 exp ! - E K k B T " (3)
D eff (T ) = D(T ) • # 1 + N t N l exp ! E b k B T " $ -1 (4)
K eff (T ) = K(T ) • # 1 + N t N l exp ! E b k B T " $ (5)
D 0 and K 0 are the pre-exponential diffusivity and solubility factors (resp. m 2 .s -1 and mol.m -3 .Pa -1 /2 ); E D (eV) is the energy barrier for the diffusion between interstitial sites, E K (eV) is the solution energy of H in interstitial sites (defectless material) and E b (eV) is the binding energy of hydrogen to traps. Detrapping energy is usually taken as the sum of binding energy E b and diffusive energy E D , as shown in figure 1. N t and N l (m -3 ) are the trap sites and lattice diffusion sites concentrations. Equations ( 4) and ( 5) bear the assumption that a single trapping site affects diffusivity and solubility of hydrogen in the considered range of temperature. If other trapping sites (non-dominant) are present, they do not affect D eff . Another measurement technique called thermal desorption spectrometry (TDS) detailed in part 3 can help us determine the properties of those secondary traps. In particular, if high energy traps are present, they would be filled during the first exposure to H and become transparent for permeation in the following measurement cycles.

Combining equations ( 1) and ( 2)-(3) gives the defining relation for permeability:

Φ(T ) = Φ 0 exp ! - E Φ k B T " (6)
where

Φ 0 = D 0 • K 0 (mol.m -1 .Pa -1 /2 .s -1 ) E Φ = E D + E K (eV)
Equations ( 2)- [START_REF] Esteban | Hydrogen transport and trapping in Eurofer97[END_REF] show in particular that the influence of traps is stronger at lower temperatures. Schematically, D eff and K eff are driven by trapping effects at lower During the permeation analysis, D eff and Φ are measured at several temperatures, after which D eff (T ) is fitted numerically using equation ( 4) with D 0 , E D , N t and E b as free parameters. The interstitial sites density is N l = 5.15 • 10 29 sites.m -3 , assuming that H sits in tetrahedral interstitial positions of a bcc structure in Eurofer97 [START_REF] Hodille | Simulations of atomic deuterium exposure in self-damaged tungsten[END_REF]. Once Φ and D eff are known, K eff is calculated using equation (1).

State of the art

Gas permeation experiments performed on Eurofer97 between 423 K and 723 K yielding diffusivity, solubility and permeability of hydrogen isotopes have been presented in [START_REF] Esteban | Hydrogen transport and trapping in Eurofer97[END_REF], [START_REF] Aiello | Hydrogen isotopes permeability in Eurofer97 martensitic steel[END_REF] and [START_REF] Chen | Deuterium transport and retention properties of representative fusion blanket structural materials[END_REF], respectively with hydrogen and deuterium. Through the evolution of effective diffusivity versus temperature, the results show that trapping has an impact on the transport parameters in the 400 K -700 K range. Effective diffusivity is consequently fitted using equation [START_REF] Hodille | Simulations of atomic deuterium exposure in self-damaged tungsten[END_REF], where E b is seen as the mean interaction energy between hydrogen isotopes and trapping sites. Based on the value of this energy, trapping in Eurofer97 is assumed to be caused by high angle grain boundaries [START_REF] Esteban | Hydrogen transport and trapping in Eurofer97[END_REF] or interfaces between martensitic laths and carbide precipitates or dislocations [START_REF] Aiello | Hydrogen isotopes permeability in Eurofer97 martensitic steel[END_REF][START_REF] Baluc | Status of reduced activation ferritic/martensitic steel development[END_REF].

Experimental setup

The diagram of the Hypertomate setup (HYdrogen PERmeation in TOkamak-relevant MATErials) is given in figure 2. In this setup, two high-vacuum chambers are separated by a thin sample. During a run of the experiment, hydrogen is introduced stepwise in one of the chambers (called the upstream part), creating a flux through the sample. This permeation flux induces a pressure rise on the downstream part. The dynamics of this pressure rise contain information about the sample material, such as diffusivity and permeability.

On both sides, the pressure is recorded using MKS ® Baratron gauges (1000 Torr full gauge on the upstream side, 0.1 Torr full scale downstream), making the measurement independent of the gas.

Hydrogen is introduced stepwise in the upstream part, where the pressure is controlled by a MKS ® Baratron pressure gauge (P 1 , 1000 Torr full scale). Volume has been minimized in order for small permeation fluxes to translate into relevant pressure signals.

Upstream and downstream pressure measurements are sampled on a multiplexer in order to measure the timelag with the method presented in figure 3. As the permeation flux is measured through the pressure rise only, it is not possible to discriminate between hydrogen and other species in the permeation flux. Thorough leak detections, careful preliminary conditioning of the vacuum chambers and background recordings are performed at each step and the background noise is systematically checked to be negligible compared to the permeation flux.

Data analysis

The first step of the analysis is the characterization of the permeation regime, which can be either diffusion-limited or surface-limited [START_REF] Esteban | Diffusion and surface-limited transport parameters of deuterium in Incoloy 800[END_REF]. The distinction is made using the exponential relationship between the steady-state permeating flux J (Pa.m -3 .s -1 ) and the upstream pressure p up (Pa):

-If J ∝ p 1 /2
up , permeation is diffusion-limited. In this case, the timelag method (see 2.4) can be used to calculate D eff , and Φ is linked to p up by

J = Φ e p
1 /2 up [START_REF] Chen | Deuterium transport and retention properties of representative fusion blanket structural materials[END_REF] where e (m) is the thickness of the material.

-If J ∝ p up , permeation is surface-limited and

J = 1 2 σk 1 p up (8)
where σk 1 (m -3 .s -1 ) is the adsorption rate constant on the upstream side of the sample. In that case, the pressure rise is influenced by the surface effects taking place on the upstream side and the bulk effects are negligible; therefore, the bulk parameters D eff and Φ cannot be extracted from the pressure rise.

This characterization step has to be performed at each temperature. After the permeation regime has be proven to be diffusion-limited, Φ is directly calculated using equation [START_REF] Chen | Deuterium transport and retention properties of representative fusion blanket structural materials[END_REF]. Effective diffusivity is evaluated using the timelag method with

τ l = e 2 6D eff (9)
This well-established method first introduced in [10] uses the asymptote of the steady-state pressure, shown in red in figure 3. The timelag is defined as the time interval between the upstream pressure introduction and the intersection of this asymptote with the time axis. In order for the timelag measurement to be valid, the upstream pressure rise has to be as close as possible to a step function.

τ l Pressure (normalised) time p up p down Intercept

Results

Measurements following the procedure described earlier have been performed on two disk-shaped samples of Eurofer97 with diameter 20 mm and thicknesses 974 and 945 µm, cut from the same rod, with temperatures between 473 and 573 K and loading pressures ranging from 10 2 to 10 5 Pa. The samples were analyzed using X-ray fluorescence and their composition is given in table 1. Prior to the measurements, the samples were mechanically polished to mirror finish and cleaned in ethanol. No annealing was performed in order to measure the properties of coarse Eurofer97.

Before each measurement, the pressure is checked to be around 1•10 -4 Pa on each side on the sample. Valve 2 is then closed and the subsequent background pressure rise is recorded in order to be subtracted from the permeation signal. In most configurations, this background signal is well under 1% of the permeation-induced pressure rise. This method results in error bars under 10% for both D eff and Φ at each temperature, showing good reproducibility for this experimental setup. Once the permeation is ruled to be diffusion-limited, the timelag method and equation ( 7) are used to determine respectively D eff and Φ. The results are shown in figures 5a (effective diffusivity D eff ), 5b (permeability Φ) and 5c (effective solubility K eff ). The plots also display previously obtained values for Eurofer97 taken from [START_REF] Esteban | Hydrogen transport and trapping in Eurofer97[END_REF] and [START_REF] Aiello | Hydrogen isotopes permeability in Eurofer97 martensitic steel[END_REF], showing that Hypertomate results are in good agreement with literature. In particular, [START_REF] Aiello | Hydrogen isotopes permeability in Eurofer97 martensitic steel[END_REF] gives permeation parameters for deuterium in Eurofer97. For the comparison to be relevant, diffusivity and permeability are corrected using a mass ratio factor given by [START_REF] Marchi | Permeability, solubility and diffusivity of hydrogen isotopes in stainless steels at high gas pressures[END_REF]: 10 -3

D H D D = ' m D m H = √ 2 (10)
10 -2 1000/T (1/K) K eff (mol.m -3 .Pa -1 /2
) Hypertomate Interpolation Esteban [START_REF] Esteban | Hydrogen transport and trapping in Eurofer97[END_REF] Aiello [START_REF] Aiello | Hydrogen isotopes permeability in Eurofer97 martensitic steel[END_REF] Chen [START_REF] Chen | Deuterium transport and retention properties of representative fusion blanket structural materials[END_REF] 400 350 300 250 200

T (°C) (c

) Hydrogen effective solubility in Eurofer97

Figure 5: Eurofer97 hydrogen permeation parameters. Fitting experimental results using equations ( 4) and ( 6) yield respectively D eff and Φ while K eff is calculated as their ratio Φ /D eff following (1).

The corresponding quantities describing hydrogen transport in Eurofer97 (D eff (T ), K eff (T ) and Φ(T )), measured using Hypertomate, compare well with the previously published values [START_REF] Esteban | Hydrogen transport and trapping in Eurofer97[END_REF][START_REF] Aiello | Hydrogen isotopes permeability in Eurofer97 martensitic steel[END_REF]. However, the transport parameters given in [START_REF] Chen | Deuterium transport and retention properties of representative fusion blanket structural materials[END_REF] differ from the other ones, especially at low temperatures. This issue was pointed out in [START_REF] Chen | Deuterium transport and retention properties of representative fusion blanket structural materials[END_REF] and calls for further work.

The corresponding Arrhenius parameters, obtained by weighted least-squares fitting, are given in table 2.

Thermal desorption spectrometry

Thermal Desorption Spectrometry (TDS) consists in heating up a previously loaded sample following a given temperature ramp-up. Desorption from the trapping sites ensues, measured by a mass spectrometer. The spectrum typically shows peaks which can be attributed to trapping sites. The work presented here includes a TDS spectrum for Eurofer97 loaded with deuterium by gas exposure, as well as the identification of the trapping and surface parameters through parametric optimization performed with the code MHIMS.

State of the art

Previously performed thermal desorption spectrometry experiments on hydrogen-or deuterium-loaded Eurofer97 show a main desorption peak around 425 K -500 K [START_REF] Chen | Deuterium transport and retention properties of representative fusion blanket structural materials[END_REF][START_REF] Malitckii | Comparative study of hydrogen uptake and diffusion in ods steels[END_REF][START_REF] Yagodzinskyy | Hydrogen effects on tensile properties of Eurofer 97 and ODS-Eurofer steels[END_REF][START_REF] Xu | Transport parameters and permeation behavior of hydrogen isotopes in the first wall materials of future fusion reactors[END_REF] for heating rates varying from 6 to 20 K/s and several loading techniques (plasma implantation, gas loading and electrochemical loading). These spectra also show a similar behavior in the low-temperature range, with a desorption rate close to zero until around 400 K.

Although the loading conditions and heating rates differ in [START_REF] Martynova | Impact of kr and ar seeding on d retention in ferritic-martensitic steels after high-fluence plasma exposure[END_REF], [START_REF] Ogorodnikova | Surface modification and deuterium retention in reduced-activation steels under low-energy deuterium plasma exposure. part I: undamaged steels[END_REF] and [START_REF] Maday | Hydrogen effects on the tensile and fatigue properties of Eurofer97[END_REF], the corresponding Eurofer97 TDS spectra present up to three desorption bands. Single peak Eurofer97 TDS spectra are witnessed in the case of plasma loading ( [START_REF] Xu | Transport parameters and permeation behavior of hydrogen isotopes in the first wall materials of future fusion reactors[END_REF]), electrochemical loading ( [START_REF] Malitckii | Comparative study of hydrogen uptake and diffusion in ods steels[END_REF][START_REF] Yagodzinskyy | Hydrogen effects on tensile properties of Eurofer 97 and ODS-Eurofer steels[END_REF][START_REF] Hollingsworth | Comparative study of deuterium retention in irradiated Eurofer and Fe-Cr from a new ion implantation materials facility[END_REF]) and short (1h) gas loading [START_REF] Chen | Deuterium transport and retention properties of representative fusion blanket structural materials[END_REF]: if these surnumerary traps are located deep enough in the bulk, short gas loading or plasma implantation do not fill them right away, but the waiting time between loading and TDS lets the atoms migrate in the material and fill these two traps.

Sample preparation and measurement

In order to supplement the analysis with a more thorough set of trapping parameters, thermal desorption spectrometry (TDS) is performed on Eurofer97 samples using the device described in [START_REF] Hodille | Hydrogen trapping in carbon film: From laboratories studies to tokamak applications[END_REF].

Similarly to permeation, the sample is not annealed in order to be representative of coarse Eurofer97. Loading is performed by heating up the sample at 400 • C in a deuterium gas atmosphere of 9•10 4 Pa for 6 hours. The conditions of this gentle loading are similar to those of permeation, which allows us to compare the two experiments. At the end of loading, the vial containing the sample and the deuterium atmosphere is quenched in room-temperature water.

The TDS experiment is conducted in a high vacuum chamber equipped with a Bayard-Alpert pressure gauge and a mass spectrometer. After each measurement, the mass spectrometer is calibrated by injecting D 2 in the chamber and recording the corresponding pressure.

Prior to heating up the sample, the residual pressure in the chamber is p bg ∼ 1•10 -6 Pa dominated by residual H 2 O. During the experiment, the loaded sample is heated following a ramp-up of 1 K/s up to 873 K. The maximum temperature is capped in order to avoid a phase change in Eurofer97 [START_REF] Oliveira | Thermodynamic simulation of reduced activation ferritic-martensitic Eurofer97 steel[END_REF]. The sample temperature is controlled and recorded by a thermocouple inserted in the sample, which guarantees a constant heating rate and an accurate measurement of the sample temperature. At the end of the TDS measurement, the sample is let to cool down in the vacuum chamber before performing the measurement a second time. This step is necessary to ensure that deuterium has fully desorbed in the temperature range covered by the experiment and also serves as a background measurement, proving that the sample heating does not interfere with the spectrometer measurements. The resulting spectrum is presented in figure 6. This desorption spectrum shows a prominent peak centered on 510 K. A second peak is appearing at the end of temperature range, interrupted.

A temperature ramp-up performed on an empty sample showed that this high temperature behaviour is not an artifact of the measurement technique. The shape of this TDS can be interpreted as a first trap with relatively low detrapping energy and high density, and a second trap with a higher detrapping energy and a seemingly lower density. The final desorption rate is not zero, meaning that some deuterium remains in the sample at the end of the experiment and that the trap in the 850 K region was not fully emptied.

D0 ED Φ0 EΦ K0 EK Nt E b (m 2 .s -1 ) (eV) (mol.m -1 .Pa -1/2 .s -1 ) (eV) (mol.m -3 .Pa -1/2 ) (eV) (m -3 ) (eV

Modeling

To estimate the kinetic trapping and detrapping parameters, a macroscopic rate equation model coupling hydrogen diffusion between interstitial sites and trapping/detrapping based on the McNabb and Forster equations [START_REF] Mcnabb | A new analysis of the diffusion of hydrogen in iron and ferritic steels[END_REF] was used. Several codes are available to describe retention and desorption from plasma facing materials and to simulate TDS experiments (MHIMS [START_REF] Hodille | Macroscopic rate equation modeling of trapping/detrapping of hydrogen isotopes in tungsten materials[END_REF], TESSIM [START_REF] Schmid | Comparison of hydrogen retention in W and W/Ta alloys[END_REF], TMAP [START_REF] Longhurst | TMAP7 user manual[END_REF], FESTIM [START_REF] Delaporte-Mathurin | Finite element analysis of hydrogen retention in ITER plasma-facing components using FESTIM[END_REF], HIDT [START_REF] Sun | Modeling and simulation for surface helium effect on hydrogen isotopes diffusion and trapping/detrapping behavior in plasma facing materials[END_REF]). These codes share a similar model for diffusion and trapping of hydrogen in a given metal. The difference lies in the boundary condition used to describe the interaction of hydrogen with the surface. The usual approach is to use a recombination coefficient [START_REF] Waelbroeck | Influence of bulk and surface phenomena on the hydrogen release out of clean metallic surfaces[END_REF]. In particular, TMAP relies on the ratedep boundary condition [START_REF] Longhurst | TMAP7 user manual[END_REF], which corresponds to the recombination rate model presented in subsection 4.2. However, the use of this recombination coefficient is based on the assumption that bulk to surface processes are fast. By using a code that implements more fundamental aspects of the surface processes such as MHIMS [START_REF] Hodille | Macroscopic rate equation modeling of trapping/detrapping of hydrogen isotopes in tungsten materials[END_REF] or TESSIM [START_REF] Schmid | Comparison of hydrogen retention in W and W/Ta alloys[END_REF], the simulations go beyond this assumption.

The simulation is done in four steps:

1. Gas loading: constant temperature of 673 K and deuterium loading pressure of 9•10 4 Pa. 2. Resting phase: temperature is lowered to 300 K.

Consequently, deuterium pressure decreases to 4•10 4 Pa. This step represents the quenching performed during the experimental preparation of the sample. 3. Pumping down: temperature remains at 300 K and deuterium pressure drops to 0 Pa, corresponding to the insertion of the sample in the vacuum enclosure. 4. TDS phase: deuterium pressure is 0 Pa on both sides, corresponding to the vacuum atmosphere of the device. Temperature is increased by 1 K/s, starting at 300 K.

The output of the simulation steps is the desorption rate as a function of temperature. This numeric result is compared to the experimental spectrum, from which a new set of parameters is inferred.

The evolution of mobile (C m ) and trapped (C t,i ) particles in the bulk is described by equations ( 11)-( 12):

∂C t,i ∂t = -ν 0 C t,i exp - E dt,i k B T +ν 0 exp - E D k B T N t,i N l C m ! 1 - C t,i N t,i " (11)
∂C m ∂t = D(T ) ∂ 2 C m ∂x 2 - ( ∂C t,i ∂t (12)
where ν 0 (s -1 ) is the pre-exponential frequency for trapping, E dt,i (eV) is the detrapping energy for trap i and is N t,i (m -3 ) the associated density.

As the sample is exposed on both sides during implantation, and as the desorption is equivalent on both sides of the sample, half of the sample is simulated and a Neumann boundary condition following ∂Cm ∂t = 0 renders the symmetry.

To simulate gas exposure, several boundary conditions are available:

-Detailed surface model decomposing the surface processes as proposed in [START_REF] Pick | A model for atomic hydrogen metal interactions -application to recycling, recombination and permeation[END_REF] and detailed in subsection 3.5;

-Sieverts' law C m0 = K √ p load where C m0 is the mobile concentration at the first layer of the bulk in mol.m -3 , giving a Dirichlet boundary condition that corresponds to equilibrium in the surface model;

-Recombination model, imposing a recombination flux J r = K r C 2 m0
and, as a counterpart, a dissociation flux J d = K d P up . The recombination rate K r can be linked to the surface rate constants as shown in [START_REF] Hodille | Simulations of atomic deuterium exposure in self-damaged tungsten[END_REF]:

K r (T ) = ν 2 bs (T )ν des (T ) ν 2 sb (T ) (13)
where ν bs is the resurfacing rate, ν d is the desorption rate and ν sb is the absorption rate, which are detailed in 3.5. These two rates follow Arrhenius laws:

K r (T ) = K r0 exp % -E Kr k B T & K d (T ) = K d0 exp % -E K d k B T with K r0 in m 4 .s -1 , K d0 in mol 2 .m -2
.Pa -1 .s -1 and recombination/dissociation energies in eV. Assuming that the surface processes are in steady-state, recombination and dissociation rates are linked by solubility:

K(T ) = ) K d (T ) K r (T ) (14)

Dynamic surface model

This model takes into account all the fluxes taking place on the surface from H 2 dissociation to H absorption in the bulk [START_REF] Hodille | Simulations of atomic deuterium exposure in self-damaged tungsten[END_REF][START_REF] Pick | A model for atomic hydrogen metal interactions -application to recycling, recombination and permeation[END_REF]. The equations initially given in [START_REF] Hodille | Simulations of atomic deuterium exposure in self-damaged tungsten[END_REF] are completed with the gas loading flux first presented in [START_REF] Pick | A model for atomic hydrogen metal interactions -application to recycling, recombination and permeation[END_REF], giving equations ( 15) -( 16):

∂C s ∂t = φ gas -φ des -φ s→b + φ b→s (15) λ ! ∂C m ∂t " x=0 = φ s→b -φ b→s -φ diff ( 16
)
where C s is the surface concentration and λ (m) is the distance between two adjacent interstitial sites, and

-φ gas = 2Γ • s • p load (1 -θ) 2
is the part of the incident flux of gas molecules that dissociates and sticks on the surface. Γ = 4.67

• 10 24 / √ M T (Pa -1 .m -2 .s -1
) is a coefficient depending on molecular mass and temperature taken from kinetic gas theory, s = s 0 exp -2E diss /kT is the initial sticking probability of the impinging atoms with s 0 (dimensionless) the exponential pre-factor for this probability and E diss (eV) half the dissociation energy barrier for the D 2 molecules on the surface. p load (Pa) is the loading pressure and θ = Cs /ns is the surface coverage, with n s (m -2 ) the number of sites on the surface.

φ des = 2ν d (T )C 2 s is the desorption flux of atoms leaving the surface as molecules. The desorption rate is ν d (T ) = ν 0,d λ 2 d e -2E des kT where ν 0,d (s -1 ) is the desorption frequency , λ d = 1

√ ns (m) is the jumping distance between two surface adsorption sites and E des (eV) is half the energy barrier corresponding to the desorption of a D 2 molecule.

φ s→b = ν sb (T )C s models the absorption of a HI atom from the surface to the bulk, assuming that the mobile concentration is low i.e. C m ≪ N l . The absorption rate is ν sb (T ) = ν 0,sb e -E A/kT where ν 0,sb (s -1 ) represents the absorption frequency and E A (eV) is the absorption energy.

φ b→s = ν bs C m0 (1 -θ) represents the release of a HI atom from the bulk to the surface. The resurfacing rate is ν bs (T ) = ν 0,bs λ bs e -E R/kT where ν 0,bs represents the frequency associated to the process (s -1 ), E R (eV) is the resurfacing energy and λ bs = ns /N l (m) the jumping distance between the first lattice site encountered by the HI and the initial adsorption site.

-

φ diff = -D(T ) % ∂Cm /∂x & x=0
models the diffusion of the absorbed atoms from the first lattice site in the bulk to deeper lattice sites.

Parametric optimization

In the optimization, all three boundary condition options have been explored in order to show the influence of surface processes on the desorption of deuterium from Eurofer97. In the fit of the experimental TDS signal, the following parameters are left free, with energies taken as constants:

-Detrapping energies E dt,i

-Trapping densities N t,i -Resurfacing energy E R -Absorption energy E A -Desorption energy E des
The last three parameters only concern the kinetic surface model. For this optimization to remain consistent with the permeation results, the free parameters E R , E A and E des are constrained by equation [START_REF] Hollingsworth | Comparative study of deuterium retention in irradiated Eurofer and Fe-Cr from a new ion implantation materials facility[END_REF], which is yielded by steady-state in ( 15)-( 16) as shown in [START_REF] Pick | A model for atomic hydrogen metal interactions -application to recycling, recombination and permeation[END_REF]:

K 0 = 1 λ abs ) Γs 0 ν 0 λ 2 des (17) E K = E A -E des -E R (18)
The rest of the parameters are fixed and given in table 3. In particular, s 0 is calculated using equation ( 17) in order to fit the solubility measured by Hypertomate and given in table 2.

Detailed surface model

In this part, the experimental TDS is fitted using the kinetic surface model described in 3.5. Results are shown in figure 7.

In this optimization, diffusivity and solubility are taken from Hypertomate. The first guess on detrapping energies and traps densities is based on the previously obtained trapping parameters, namely E b and N t .

The first round of optimization is performed using two traps. This configuration gives an acceptable fit regarding the peaks width and centering along the temperature axis, showing a good prediction of trapping energies. However, this fit does not predict the desorption taking place around 750 K, as shown by the discrepancy between the green plot and the experimental points on figure 7. In order to bridge this gap, a third trapping site needs to be added to improve the fit in the 750 K region. For the sake of comparison, the result obtained with only one trap is also shown.

The parameters values obtained with the surface model fit are given in table 4.

Discussion

Permeation and TDS coupling

Permeation experiments have resulted in values for diffusivity D(T ), solubility K(T ) and permeability Φ(T ) that have been shown to be in good agreement with permeation literature. These experiments also yielded trapping parameters corresponding to the dominant trapping site in Eurofer97, namely E b and N t :

E b = 0.42 eV N t = 6.08 • 10 25 m -3
As shown in figure 1, detrapping, binding and diffusion energies are linked by

E dt = E D + E b
Using the values determined by permeation (E D = 0.16 eV and E b = 0.42 eV), we obtain E dt, permeation = 0.58 eV This value is close to the detrapping energy corresponding to the main trapping site (E dt,1 = 0.51 eV) visible on the TDS of Eurofer97. The trapping densities found with both experiments are also very close, showing that the two experiments are in agreement. Furthermore, we can see that N t,1 ≫ N t,2 and N t,1 ≫ N t,3 , which confirms the assumption made in equations ( 4) and ( 5) that trapping is dominated by one site in Eurofer97.

Recombination rate and solubility models

This part presents simulations performed with the recombination rate and solubility surface models. Traps energies and densities keep the values found in the previous part. In the solubility model, the boundary condition follows Sieverts' law with K taken from Hypertomate. In the recombination model, recombination rate K r is evaluated using equation [START_REF] Yagodzinskyy | Hydrogen effects on tensile properties of Eurofer 97 and ODS-Eurofer steels[END_REF], meaning that the kinetic surface processes are assumed to be in steady-state.The line for the sake of comparison. The solubility simulation shows a satisfying fit in the 600 K -800 K region, but fails to replicate the initial desorption taking place in the 300 K -500 K range. This higher initial desorption also lowers the main peak, as particles get desorbed at lower temperatures. Recombination and solubility solutions coincide in the 300 K -900 K temperature range presented in figure 8, which was to be expected, as both models correspond to the steady-state of the surface model. In particular, this shows that steady-state of the kinetic surface model cannot be assumed to simulate properly this TDS experiment. Based on the activation energies used in the kinetic surface model, the limiting step for desorption is the bulk to surface transition, which energy barrier is 0.65 eV. This energy barrier inhibits the bulk to surface processes under 400 K: in this range, desorption is limited by the bulk to surface flux, i.e. the assumption that all processes are in steady-state does not apply. Thus, the solubility and steady-state recombination rate models cannot be used to describe the initial desorption profile seen in the experimental results.

However, if equation ( 13) is no longer taken into account, K r0 and E Kr become free parameters that can be optimized in order to replicate the kinetic surface model with a recombination rate model. The spectra resulting from this optimization, performed using the optimization routine first introduced in [START_REF] Delaporte-Mathurin | Parametric optimisation based on TPD experiments for rapid and efficient identification of hydrogen transport materials properties[END_REF], are presented in figure 9. This plot shows that the recombination rate can be used to replicate the experimental desorption flux in the 300 K -500 K range. The optimized recombination energy is now positive (meaning that the recombination rate will increase with the temperature) and switches from -0.54 eV to 0.73 eV. Thus, the kinetic surface model and the optimized recombination rate model both require the limiting step to have an activation energy of about 0.7 eV. However, the optimized recombination rate is just a fitting parameter with little physical meaning, while the kinetic surface model corresponds to a more careful description of the elementary processes taking place at the surface. Thus, the recombination rate needs to be handled with caution as discussed in [START_REF] Schmid | On the use of recombination rate coefficients in hydrogen transport calculations[END_REF].

Conclusion

Using the Hypertomate setup, the gas permeation technique has been used on Eurofer97 to determine its hydrogen diffusivity, solubility and permeability. The results were shown to be in good agreement with literature and showed that effective diffusivity and solubility are affected by traps in the investigated temperature range. In order to get a more detailed understanding of this trapping behaviour, TDS experiments on deuterium-loaded samples were performed to further investigate trapping in Eurofer97. The resulting spectrum presents several traps, similarly to [START_REF] Martynova | Impact of kr and ar seeding on d retention in ferritic-martensitic steels after high-fluence plasma exposure[END_REF]. In this work, the multiple-peak TDS spectrum was fitted using the MHIMS code including a kinetic surface model describing the elementary processes of hydrogen on Eurofer97 surface: dissociation, adsorption, desorption, absorption and (re-)surfacing. This parametric fit showed that three traps are necessary to model the trapping behaviour of deuterium in Eurofer97, which energies and densities were evaluated. The results of both experiments were shown to be in good agreement. In order to assess the improvement brought to the simulations by the kinetic surface model, fits with a recombination boundary condition were performed. As the initial value of K r evaluated from the surface parameters failed to describe the initial desorption, a fit was performed with K r left as a free parameter, showing that the recombination parameter can be tuned to fit the TDS, although the corresponding values differ widely from the initial ones and no longer describe the surface processes. A set of parameters describing adsorption and absorption of deuterium on the surface was also given, based on the TDS fit.

Figure 1 :

 1 Figure 1: Energy diagram presenting diffusion energy E D , trapping energy Et, binding energy E b and detrapping energy E dt . This diagram illustrates the assumption that E D = Et and explicits the relation between trapping, detrapping and diffusion : E dt = E D +E b

Figure 2 :

 2 Figure 2: Diagram presenting the Hypertomate hydrogen permeation experiment

Figure 3 :

 3 Figure 3: Timelag method illustration

  Particles flux (D/m"/s)

Figure 6 :

 6 Figure 6: Thermal desorption spectrometry spectrum of deuterium in Eurofer97. The sample was loaded under deuterium atmosphere of 9•10 4 Pa at 673 K for 6 hours. TDS was then performed with a 1 K/s temperature ramp-up.

Figure 7 :

 7 Figure 7: Experimental and simulated TDS spectra. Optimization is performed using Hypertomate solubility and the kinetic surface model. Fits obtained with one, two and three traps are displayed, showing that the trapping behaviour is best represented with three traps.

Figure 9 :

 9 Figure 9: Experimental and simulated TDS spectra comparing the kinetic surface model with the optimized Kr model

Table 1 :

 1 Eurofer97 detailed composition

	Element Z Concentration (mass ratio)
	Fe	26	82.88 %
	Cr	24	8.84 %
	Cu	29	2.46 %
	W	74	1.11 %
	Na	11	0.59 %
	Mn	25	0.53 %
	V	23	0.19 %
	Se	34	547 ppm
	Co	27	213 ppm
	K	19	201 ppm
	Zn	30	193 ppm

Diffusion-limited behaviour is validated by plotting log(J) versus log(p up ) and checking that the slope is close to 0.5. An example of this validation, performed at 300°C, is given in figure

4

. up

Figure 4: Permeation regime characterisation

Table 2 :

 2 Arrhenius permeation parameters for hydrogen isotopes in Eurofer97

	)

Table 3 :

 3 Simulation parameters used in the MHIMS code to replicate Eurofer97 TDS

	Name	Symbol Value	Unit	Reference
	Diffusivity pre-exponential factor Diffusion energy	D0 ED	2.52 • 10 -7 0.16	m 2 .s -1 eV	this work this work
	Solubility pre-exponential factor Solubility energy	K0 EK	1.76 • 10 -1 0.27	mol.m -.Pa -1/2 this work eV this work
	Sticking probability pre-exponential factor Dissociation energy	s0 E diss	1.61 • 10 -4 0	dimensionless [27] eV [28]
	bcc iron lattice constant	a	2.856	Å	[29]
	Eurofer97 density Number of adsorption sites on the surface Number of interstitial sites in the bulk Distance between two adjacent interstitial sites	ρEu ns N l λ ν 0,d	m -3 m -2 6 • ρEu = 5.15 • 10 29 m -3 2 /a 3 = 8.59 • 10 28 ρ 2/3 Eu = 1.95 • 10 19 a /2 √ 2 = 101 • 10 -12 m	[4] [4] [4]
	Pre-exponential desorption, absorption, resurfacing and detrapping frequencies	ν 0,sb ν 0,bs	1 • 10 13	s -1	[4]
		ν0

Table 4 :

 4 Detrapping energies, traps densities and surface energies E A , E R and E des resulting from the TDS fit with MHIMS and a kinetic surface model values taken for E A , E R , E des are those obtained in the kinetic surface model fit. The resulting values for K r and K d are given in table 5 and a visualization of the simulated TDS spectra corresponding to each of the three models is given in figure8.

	Name	Symbol	Value	Unit
		E dt,1	0.51
	Detrapping energies E dt,2	1.27	eV
		E dt,3	1.65
	Traps densities	Nt,1 6.01 • 10 25 Nt,2 6.44 • 10 22 m -3 Nt,3 3.88 • 10 23 EA 1.32
	Surface energies	ER	0.65	eV
		2E des	0.80

Table 5 :

 5 Recombination and dissociation rates used in the recombination rate simulations. Surface model equivalent is based on[START_REF] Yagodzinskyy | Hydrogen effects on tensile properties of Eurofer 97 and ODS-Eurofer steels[END_REF] and surface energies obtained from the fit, optimized ones come from the parametric optimization described in 4.2.

Acknowledgments

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement no. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.


Rémi Delaporte-Mathurin
email: remi.delaporte-mathurin@cea.fr

Mykola Ialovega

Etienne A Hodille

Jonathan Mougenot

Yann Charles

Elodie Bernard

Céline Martin

Christian Grisolia

Influence of exposure conditions on helium transport and bubble growth in tungsten

Helium diffusion, clustering and bubble nucleation and growth is modelled using the finite element method. The existing model from Faney

In fusion devices, extreme fluxes of helium (He) and hydrogen (H) are expected. These fluxes will be mostly located on the tungsten (W) divertor which will also exhaust the plasma He ashes.

Due to a strong W-He repulsion, interstitial He in W tend to react and form clusters [START_REF] Becquart | A density functional theory assessment of the clustering behaviour of He and H in tungsten[END_REF][START_REF] Becquart | Migration energy of He in W revisited by ab initio calculations[END_REF] (see Fig. 1). Eventually, when mobile clusters reach a critical size, trap-mutation (also called self-trapping) will occur. Frenkel pairs (selfinterstitial W atom and vacancy) will be produced and the mobile clusters will sit in the created vacancies [START_REF] Boisse | Modeling of the self trapping of helium and the trap mutation in tungsten using DFT and empirical potentials based on DFT[END_REF] . At this point, the clusters are immobile and will continue to grow by absorbing small helium clusters and act as nuclei for bubble formations. When the cluster is over-pressurised, additional Frenkel pairs will be created. He bubbles can then form in W and their morphology depends on the exposure conditions [START_REF] Ialovega | Hydrogen trapping in tungsten: Impact of helium irradiation and thermal cycling[END_REF][START_REF] Taylor | Investigating helium bubble nucleation and growth through simultaneous in-situ cryogenic, ion implantation, and environmental transmission electron microscopy[END_REF][START_REF] Qin | Helium bubble nucleation at grain boundaries and its influence on intergranular fracture[END_REF][START_REF] Lemahieu | He irradiation on tungsten exposed to ELM-like thermal shocks[END_REF] . Such bubbles have been observed using Molecular Dynamics (MD) [START_REF] Boisse | Modeling of the self trapping of helium and the trap mutation in tungsten using DFT and empirical potentials based on DFT[END_REF][START_REF] Hammond | Helium flux effects on bubble growth and surface morphology in plasma-facing tungsten from large-scale molecular dynamics simulations[END_REF][START_REF] Hamid | Molecular dynamics simulations of helium clustering and bubble growth under tungsten surfaces[END_REF][START_REF] Bergstrom | A molecular dynamics study of subsurface hydrogen-helium bubbles in tungsten[END_REF][START_REF] Maroudas | Helium segregation on surfaces of plasma-exposed tungsten[END_REF][START_REF] Sefta | Surface response of tungsten to helium and hydrogen plasma flux as a function of temperature and incident kinetic energy[END_REF] and Object Kinetic Monte Carlo [START_REF] Valles | Temperature dependence of underdense nanostructure formation in tungsten under helium irradiation[END_REF][START_REF] De Backer | Modeling of helium bubble nucleation and growth in austenitic stainless steels using an Object Kinetic Monte Carlo method[END_REF] . He bubbles can alter the mechanical properties of W [START_REF] Qin | Helium bubble nucleation at grain boundaries and its influence on intergranular fracture[END_REF][START_REF] Das | Hardening and strain localisation in helium-ion-implanted tungsten[END_REF][START_REF] Nguyen | Modeling the effects of helium-vacancy clusters on the stress-strain response of a grain boundary in iron by a mechanistic finite element approach informed by molecular dynamics data[END_REF] and reduce the thermal properties of components [START_REF] Qu | Degradation of thermal conductivity of the damaged layer of tungsten irradiated by helium-plasma[END_REF][START_REF] Cui | Thermal conductivity reduction of tungsten plasma facing material due to helium plasma irradiation in PISCES using the improved 3-omega method[END_REF] . Eventually, when over-pressurised bubbles are close to the surface, bursting can occur as shown by Sefta et al. in MD simulations [START_REF] Sefta | Surface response of tungsten to helium and hydrogen plasma flux as a function of temperature and incident kinetic energy[END_REF][START_REF] Sefta | Helium bubble bursting in tungsten[END_REF] . Bursting greatly modifies surface morphologies by increasing the roughness and producing craters [START_REF] Ialovega | Hydrogen trapping in tungsten: Impact of helium irradiation and thermal cycling[END_REF] and W-fuzz [START_REF] Lemahieu | He irradiation on tungsten exposed to ELM-like thermal shocks[END_REF][START_REF] De Temmerman | Effect of high-flux H/He plasma exposure on tungsten damage due to transient heat loads[END_REF][START_REF] Mccarthy | Enhanced fuzzy tungsten growth in the presence of tungsten deposition[END_REF][START_REF] Khan | Helium irradiation effects on the surface modification and recrystallization of tungsten[END_REF][START_REF] Bernard | Temperature impact on the micro structure of tungsten exposed to He irradiation in LHD[END_REF] . Moreover, He exposure can alter H retention in W [START_REF] Ialovega | Hydrogen trapping in tungsten: Impact of helium irradiation and thermal cycling[END_REF][START_REF] Markelj | Hydrogen isotope accumulation in the helium implantation zone in tungsten[END_REF][START_REF] Ogorodnikova | Deuterium retention in tungsten exposed to low-energy pure and helium-seeded deuterium plasmas[END_REF][START_REF] Baldwin | Effect of He on D retention in W exposed to low-energy, high-fluence (D, He, Ar) mixture plasmas[END_REF][START_REF] Miyamoto | Microscopic damage of tungsten exposed to deuterium-helium mixture plasma in PISCES and its impacts on retention property[END_REF][START_REF] Ueda | Simultaneous irradiation effects of hydrogen and helium ions on tungsten[END_REF] .

An effort has been made to assess He transport in W using atomistically informed macroscopic models called cluster dynamics models implemented using finite differences. Some examples of such implementations are the work of Faney et al. [START_REF] Faney | Spatially dependent cluster dynamics modeling of microstructure evolution in low energy helium irradiated tungsten[END_REF][START_REF] Faney | Spatially dependent cluster dynamics model of He plasma surface interaction in tungsten for fusion relevant conditions[END_REF] and the Xolotl code [START_REF] Blondel | Continuum-scale modeling of helium bubble bursting under plasmaexposed tungsten surfaces[END_REF] . The simulated results are promising but these codes require substantial computational resources considering the thousands of equations that need to be solved. The current study therefore proposes an approach to further simplify these models so that they can be easily implemented in finite element codes and later coupled to H transport modelling codes such as FESTIM [START_REF] Delaporte-Mathurin | Finite element analysis of hydrogen retention in ITER plasma facing components using FESTIM[END_REF] . This could serve as a base for H-He coupled simulations to better assess H transport in plasma facing materials and couple with additional physics like heat transfer.

The simplified model presented in this work is applied on a simple case and compared with existing results of the literature to ensure the model is not over-simplified. Then, the influence of exposure conditions is investigated by running a parametric study varying temperature and implanted particle flux. The results of this parametric study are analysed using a regression method previously employed [START_REF] Delaporte-Mathurin | Parametric study of hydrogenic inventory in the ITER divertor based on machine learning[END_REF] . Experiments are then conducted to quantitatively assess the He bubble density and size in He irradiated W. The current model is finally compared to these experimental results.

Methodology

This section describes the He transport model and the grouped approach employed to simplify it.

Helium clustering model.

This model describes the evolution of the concentrations of pure interstitial He clusters (He x ) and mixed He-vacancies clusters (He x V y ) that are formed by trap-mutation events.

The spatio-temporal evolution of each species of size i is defined by:

In Equation 1, the first term of the right hand-side is the diffusion term where

D = D 0 • exp(-E diff /(k B • T))
is the thermally activated diffusion coefficient expressed in m 2 s -1 with E diff the diffusion activation energy in eV, k B the Boltzmann constant in eV K -1 and T the temperature in K. If a species i is assumed to be immobile, its diffusion coefficient D i is zero. Ŵ i is the external production rate of species i.

The term R i is the coupling term due to reactions between species. A simple reaction between two species can be described as:

The forward rate constant k + A,B is the clustering rate and is calculated using the theory of diffusion-limited reactions [START_REF] Goldstein | Diffusion limited reactions[END_REF] :

where r A and r B are the capture radii and D A and D B are the diffusion coefficients of species A and B respectively. The backward rate constant k - A,B is the dissociation rate and is obtained using chemical equilibrium principles [START_REF] Goldstein | Diffusion limited reactions[END_REF] :

where ρ is the atomic density in m -3 ( ρ = 6.3 × 10 28 m -3 for W), k B is the Boltzmann constant in eV K -1 , T is the temperature in K and E b is the binding energy for the reaction AB → A + B in eV. The reaction term R i is the coupling term between concentrations and is expressed as:

In Equation 5, c i is the concentration of clusters of size i in m -3 . The first term corresponds to the reactions producing clusters of size i. The second one corresponds to the ones reacting with clusters of size i. The third term accounts for bigger clusters dissociating. Finally, the last term corresponds to clusters of size i dissociating.

Grouped approach.

Extending this clustering model to clusters containing millions of helium extremely increases the computational cost. A grouped approach proposed by Faney et al. [START_REF] Faney | Spatially dependent cluster dynamics modeling of microstructure evolution in low energy helium irradiated tungsten[END_REF] for reducing the number of equations will therefore be employed. This technique consists in grouping the big clusters that have a similar behaviour in a single equation while explicitly accounting for smaller clusters.

The clustering equations can be written as follows:

(1) where N is some threshold required for the grouping technique.

∂c i ∂t = ∇ • (D i ∇c i ) + Ŵ i + R i (2) A + B⇋ k + A,B k - A,B AB (3) k + A,B = 4π(r A + r B )(D A + D B ) (4) k - A,B = ρk + A,B e -E b k B T (5) R i = m k + m,i-m c m c i-m -c i m k + i,m c m + k - i+1 c i+1 -k - i c i (6a) ∂c 1 ∂t = ∇ • (D 1 ∇c 1 ) + Ŵ + N i=2 k - i c i -2k + 1,1 c 2 1 - N i=2 k + 1,i c 1 c i - ∞ i=N+1 k + 1,i c 1 c i
In order to simplify this model, the following quantities are defined:

Clusters with more than N He ( c b ) will be referred as "bubbles" in the following. Equation 6 therefore reads:

The mean radius of pure He clusters [START_REF] Faney | Spatially dependent cluster dynamics model of He plasma surface interaction in tungsten for fusion relevant conditions[END_REF] is given by: with r He 1 = 0.3 nm.

Several assumptions are made:

• The average radius is assumed to be a function of i b :

(6b) ∂c 2 ∂t = ∇ • (D 2 ∇c 2 ) -k + 1,2 c 1 c 2 + k + 1,1 c 2 1 -k - 2 c 2 + k - 3 c 3 (6c)
. . .

∂c i ∂t = -k + 1,i c 1 c i + k + 1,i-1 c 1 c i-1 -k - i c i (6d) ∂c i+1 ∂t = -k + 1,i+1 c 1 c i+1 + k + 1,i c 1 c i . . . (7) c b = ∞ i=N+1 c i : total concentration of clusters containing more than N He (8) �i b � = 1 c b ∞ i=N+1 ic i : average He content in c b (9) �r b � = 1 c b ∞ i=N+1 r i c i : average radius in c b (10) �k + b � = 1 c b ∞ i=N+1 k + 1,i c i = 4πD 1 (r 1 + �r b �) : average clustering rate in c b (11a) ∂c 1 ∂t = ∇ • (D 1 ∇c 1 ) + Ŵ + N i=2 k - i c i -2k + 1,1 c 2 1 - N i=2 k + 1,i c 1 c i -�k + b �c 1 c b (11b) ∂c 2 ∂t = ∇ • (D 2 ∇c 2 ) -k + 1,2 c 1 c 2 + k + 1,1 c 2 1 -k - 2 c 2 + k - 3 c 3 (11c)
. . .

∂c N ∂t = -k + 1,N c 1 c N + k + 1,N-1 c 1 c N-1 -k - N c N (11d) ∂c b ∂t = k + 1,N c 1 c N (11e) ∂(�i b �c b ) ∂t = (N + 1)k + 1,N c 1 c N + �k + b �c 1 c b (12) r He x = r He 1 + 3 4π a 3 0 10 x 1/3 - 3 4π
a 3 0 10 1/3 ( 13
)
�r b � = r(He �i b � V �i b �/4 ) = r He 0 V 1 + 3 4π a 3 0 2 �i b � 4 1/3 - 3 4π
a 3 0 2 1/3
with a 0 = 0.318 nm the lattice parameter and r He 0 V 1 = a 0 √ 3/4 . The average radius r b is assumed to be only dependent on i b . The number of vacancies in bubbles is assumed to be i b /4 . This assumption is motivated by MD computations showing that trap mutation events occur for every four additional helium in large vacancy-helium clusters. Moreover, theoretical models for He bubbles growth in metal suggest a similar trend [START_REF] Hammond | Theoretical model of helium bubble growth and density in plasma-facing metals[END_REF] .

• Dissociation of large clusters is neglected (i.e. k - i = 0 for i > N ). Indeed, the activation energy for trap mutation events is lower than that of He or vacancy emission [START_REF] Boisse | Modeling of the self trapping of helium and the trap mutation in tungsten using DFT and empirical potentials based on DFT[END_REF] . Dissociation of large clusters by vacancy or He emission is therefore negligible.

The current implementation further simplifies Faney's model 29 :

• Interactions with self-interstitial atoms or pre-existing vacancies are not taken into account. In this work, the only dissociations are He emissions from small mobile clusters and trap-mutation for large clusters. It was showed that this assumption did not have an impact on the results (see Fig. 2). • The only clusters explicitly computed are He x≤6 (i.e. N = 6 ) whereas Faney's work explicitly accounted for clusters up to V 50 He 250 and solved a bigger system of equations. The influence of this threshold N above which clusters are integrated in the quantity c b is discussed in section "Influence of N". • Clusters containing more than six He atoms are assumed to be immobile (i.e. D i = 0 for i > 6 ) due to trap mutation events. This assumption is motivated by DFT and MD results suggesting that the self-trapping energy is below the binding energy of one He atom in a pure He cluster for clusters containing more than five He atoms [START_REF] Boisse | Modeling of the self trapping of helium and the trap mutation in tungsten using DFT and empirical potentials based on DFT[END_REF] . For smaller clusters ( He 1 , He 2 , ..., He 6 ) the diffusion coefficient and the dissociation by He emission energy vary with the number of He atoms in the cluster (see Table 1). This He transport model was implemented in Python and solved using the finite element solving platform FEniCS [START_REF] Alnaes | The FEniCS project version 1.5[END_REF] . All plots in this work were generated with Matplotlib [START_REF] Hunter | Matplotlib: A 2D graphics environment[END_REF] .

Results

In this section, the current implementation is first compared with the one from Faney [START_REF] Faney | Spatially dependent cluster dynamics model of He plasma surface interaction in tungsten for fusion relevant conditions[END_REF] to ensure the additional assumptions do not produce different results. A standard half-slab case is then described and a parametric study is performed by varying the exposure conditions. Finally, the model is compared against experimental data. Table 1. Pure He clusters properties in W. Diffusion properties are taken from Faney et al. [START_REF] Faney | Spatially dependent cluster dynamics model of He plasma surface interaction in tungsten for fusion relevant conditions[END_REF] and binding energies are taken from Becquart et al. [START_REF] Becquart | Microstructural evolution of irradiated tungsten: Ab initio parameterisation of an OKMC model[END_REF] . Tendril case. The tendril application case described in [START_REF] Faney | Spatially dependent cluster dynamics model of He plasma surface interaction in tungsten for fusion relevant conditions[END_REF] was simulated in 1D with the current implementation and the results were compared. The domain size is 30 nm and the volumetric source term is described as follows:

Cluster D 0 (m 2 s -1 ) E diff (eV) E b (eV) He 1 2.95 × 10 -8 0.
where ϕ imp = 5 × 10 25 m -2 s -1 is the implanted He flux and f(x) is a Gaussian distribution with a mean value µ = R p = 1.5 nm and a standard deviation σ = 1 nm which corresponds to a 100 eV He implantation based on SRIM computations [START_REF] Ziegler | The stopping and range of ions in matter[END_REF] . Mobile He clusters concentrations were set to zero at the tendril's surfaces ( x = 0 nm and x = 30 nm).

Concentration profiles computed by the current implementation showed good agreement with the ones obtained by Faney et al. [START_REF] Faney | Spatially dependent cluster dynamics model of He plasma surface interaction in tungsten for fusion relevant conditions[END_REF] (see Fig. 2). The discrepancies are likely due to a difference in the set of dissociation energies that have been used. Indeed, at low temperature, where dissociation is not activated, the discrepancies were very small whereas at high temperature, differences increased because dissociation became more dominant. When c b is small compared to c He 1 , the equilibrium of He 1 is independent of these dissociation energies and the profiles for He 1 are identical.

Moreover, increasing the temperature tended to inhibit bubble formation in the tendril. This was explained by a greater increase in the dissociation rate and in losses at surfaces than the increase in the clustering rate. This observation is in agreement with MD results simulating He implantation in tendrils [START_REF] Wei | A better nanochannel tungsten film in releasing helium atoms[END_REF][START_REF] Wei | Understanding the release of helium atoms from nanochannel tungsten: A molecular dynamics simulation[END_REF] . The current implementation and the additional assumptions that were made are therefore valid.

Half-slab case. He transport was simulated in a 1D semi-infinite W slab. This case is the standard case describing the main quantities of interest of the parametric study performed in section "Influence of N".

The domain size is 0.1 mm which is much greater than the penetration depth of He in the simulations. 100 eV He were implanted in the first 1.5 nm as in section "Tendril case". The implanted flux was 1 × 10 22 m -2 s -1 and the temperature was 1000 K.

At low fluences, He diffused really quickly into the bulk (see Fig. 3) and the bubbles concentration c b was found to be zero. As the fluence increased, bubbles started to appear and acted as strong sinks for mobile He. This lead to a great decrease in the mobile He concentration profile.

It is worth noticing the maximum of c b was not located at the maximum of c He 1 which is the implantation depth R p . This was explained by the diffusion of small mobile clusters as shown by analytical models [START_REF] Krasheninnikov | On helium cluster dynamics in tungsten plasma facing components of fusion devices[END_REF] . As He clusters, small mobile clusters diffuse deeper into the bulk until trap-mutation occurs and bubbles nucleons (clusters with more than 6 He) are created. From that point, bubbles are formed relatively far from the surface. Because He is implanted in the first nanometres, c He 1 is maximum at R p = 1.5 nm and interactions with bubbles are stronger in this region. This tends to draw the maximum location of c b towards the surface.

The He content in bubbles i b and the radius r b were computed. After 10 s of implantation, bubbles located in the near surface contained up to 3 × 10 7 He. The maximum of r b was found to be very close to the surface at approximately 2 nm (see Fig. 4). This is explained by the high concentration of mobile He in this near surface region. Moreover, a bursting zone can be defined by the region where r b is greater than the depth of the bubble. In this region, bubble of this size would have likely burst.

From this average He content in bubbles and from Equations 13 and 12 expressing the clusters radii, the average radius r can be computed as:

The average content of He in all clusters i is computed in a similar way: These two quantities are comparable to the ones obtained by Faney et al. [START_REF] Faney | Spatially dependent cluster dynamics model of He plasma surface interaction in tungsten for fusion relevant conditions[END_REF] . After 100 s of exposure, the average radius 50 nm below the surface was above 10 nm (see Fig. 5). Moreover, the location of the maximum of these quantities move towards the exposed surface.

(14) Ŵ(x) = ϕ imp f (x) (15) �r� = ∞ i=1 c i r i ∞ i=1 c i = N i=1 c i r i + c b �r b � N i=1 c i + c b
The average radius r cannot be easily compared to experimental observations for it includes contributions from very small mobile He x clusters which are not visible experimentally (only bubbles with a radius greater than 1-3 nm are observable depending on the observation technique).

Influence of exposure parameters on He bubble growth. The impact of He flux and temperature T was studied on the case described in section "Half-slab case" in order to identify trends. Behaviour laws are identified and can be used to obtain information on He transport without needing to run any simulation.

Parametric study. A parametric study was performed by varying the implanted flux ϕ imp between 1 × 10 17 m -2 s -1 and 5 × 10 21 m -2 s -1 and the sample temperature T between 100 K and 1200 K.

Several quantities were computed. First the bubbles inventory is defined as:

The total helium inventory is calculated by: Figure 5. Average helium content i and average radius r in all clusters (mobile and bubbles) in W exposed to 100 eV He at 10 22 m -2 s -1 and 1000 K.

(16) �i� = ∞ i=1 c i i ∞ i=1 c i = N i=1 c i i + c b �i b � N i=1 c i + c b (17
The spatial mean helium content in bubbles can be computed as:

The approximation made in Equations 18 and 19 is valid as long as �i b �c b dx ≫ N i=1 ic i dx (i.e. the He inventory is dominated by that of the bubbles). This is the case in these simulations because N = 6 (the influence of this parameter is discussed in section "Influence of N").

More than 160 simulations were performed simulating 1 h of exposure. For each simulation, the quantities of interest described above were computed. A Gaussian regression process [START_REF] Bowman | C-bowman/inference-tools: 0[END_REF] was used to interpolate the data based on Bayesian inference as done in [START_REF] Delaporte-Mathurin | Parametric study of hydrogenic inventory in the ITER divertor based on machine learning[END_REF] (see Fig. 6). The temporal evolution of these quantities was also assessed (see Fig. 7).

After 1 h of exposure, the bubbles inventory I bubbles shows a weak dependence on temperature at high tem- perature and a weak dependence on the implanted flux at low temperature (see Fig. 6a). I bubbles varies from 4 × 10 12 bubbles m -2 at high temperature and low flux to 2 × 10 19 bubbles m -2 at low temperature and high flux.

The He inventory I varies from 8 × 10 16 m -2 at high temperature and low flux to10 25 m -2 at low temperature and high flux (see Fig. 6b). For temperatures above 600 K, the temperature dependence is rather weak compared to the flux dependence.

For temperatures above 300 K, and after 1 h of exposure, the sample temperature does not impact the value of � īb � (see Fig. 6c). The mean He content increases with the implanted flux as expected and varies between 10 3 He at low flux and 5 × 10 8 He at high flux.

The position of the maximum of i b tended to increase with temperature and decrease with implanted flux (see Fig. 6d). After 1 h of exposure, it was found to be really close to the surface down to 0.1 nm at low temperatures and high fluxes. The validity of the model in this region of the parameter space is questionable considering that the bubble radius is greater that the thickness of the ligament between the edge of the bubble and the surface. Such a bubble would therefore have burst before reaching this size. For each simulation point, the temporal evolution of the quantities described above has been computed. To better identify the time series on the ϕ imp , T plane, lines have been coloured according to the parameter c He 1 ,ideal which is a function of both the implanted flux and the temperature (see Equation 20) expressed in m -3 . where ϕ imp is the implanted flux, D is the diffusion coefficient of mobile He 1 in W (see Table 1), R p = 1.5 nm is the implantation depth and T is the temperature in K.

All these quantities showed a similar behaviour in time even though the kinetics were found to be different (see Fig. 7). For instance, for each (T, ϕ imp ) couple, I bubbles first increased as a power law of time before reaching a maximum (see Fig. 7a). The total He inventory I increased with time and for each simulation point but the growth rate decreased at long exposure times (see Fig. 7b). This phenomenon is explained in details in section "Influence of N". Similarly, � īb � could be written as a power law of time described in Eq 21 (see Fig. 7c). The depth of the maximum of i b tended to decrease with time as it was observed in section "Half-slab case" (see Fig. 7a).

Identification of regimes for inventory evolution. For every (T, ϕ imp ) couple, I bubbles increased rapidly at low flu- ences until reaching a maximum at high fluences (see Fig. 7a). On the other hand, the mean He content � īb � was constant at low fluences and increased as a power law of time at high fluences (see Fig. 7c). The evolution of � īb � can be described as:

where N = 6 in this model, a and b depend on (T, ϕ imp ) . The choice of N = 6 in this model is detailed in section "Influence of N". The total He inventory I being the product of these two quantities, two different growth rates were observed (see Figs. 7b and8).

This phenomenon can be attributed to two different regimes. The first regime is the nucleation regime where new bubbles nucleons are created (i.e. c b and I bubbles increase). In the nucleation regime, the bubble concentra- tion c b and the capture radius r b are too low for the He content in bubbles i b to increase significantly (i.e. � īb � is constant). The second regime is the bubble growth regime. In this regime, c b is high enough for interactions between bubbles and mobile He to occur. Implanted interstitial He atoms ( c He 1 ) therefore interact preferably with bubbles rather than clustering with other interstitial He atoms. This means that no additional bubbles nucleons are created (i.e. c b reaches a maximum). Because interactions between bubbles and mobile He are strong, the term �k + b �c 1 c b in Equation 11 becomes significant and the He content increases (i.e. � īb � increases). This is illustrated by the thickness of the interaction arrows in Fig. 1.

Influence of N.

In order to assess the impact of the parameter N in Equation 11, the evolution of the He inventory I, the mean He content in immobile clusters (different from � īb � ) and the bubbles inventory I bubbles was computed with several values of N.

The flux of 100 eV He in this test case was 10 20 m -2 s -1 and the temperature was 1000 K. It was shown that varying N had no impact on these quantities whatsoever (see Fig. 9). This highlights the very quick transition from nucleation regime to growth regime in this model.

The number of equations that need to be solved can therefore be minimised by setting the parameter N to its minimum ( N = 6 ) without losing accuracy in the results. This minimum value corresponds to the number of mobile clusters which have to be explicitly simulated in order to account for all the diffusion mechanisms.

Comparison with experiments.

He implantation experiments were performed on W in the linear plasma device PSI-2 [START_REF] Kreter | Linear plasma device PSI-2 for plasma-material interaction studies[END_REF] . W was irradiated with 75 eV He at 2.3 × 10 22 m -2 s -1 and 1053 K for 13 s. A thin lamella for cross-sectional observations was prepared using the FIB (Focused Ion Beam) technique with a Dual Beam FIB (FEI Helios 600 NanoLab). Prior to FIB cutting, the surface of the sample was coated with a SiO layer for better contrast and then with a protective platinum layer to avoid damaging the surface during the lamella preparation. Cross-sectional observations of the He-implanted W were performed using transmission electron microscopy in a TEM FEI Titan 80-300 apparatus. A typical TEM image of the lamella is presented in Fig. 10. Comparison of under-and over-focused TEM images allowed identification of the bubbles. Bubbles were observed up to 100 nm with larger bubbles closer to the surface and smaller bubbles deeper in the bulk. Open bubbles and holes at the surface were also observed suggesting bursting events occurred. This is in accordance with what was observed in the simulations (see Fig. 4).

A procedure was developed to automate the bubble detection on TEM images using the ImageJ software [START_REF]Fiji: an open-source platform for biological-image analysis |Nature Methods[END_REF] . The area of bubbles were computed as well as their diameter assuming a spherical shape for the bubbles. Bubble density and size as a function of depth was therefore computed using 12 pairs of under-and over-focused TEM images. The bubble density was found to range from 7 × 10 19 m -3 to 2 × 10 20 m -3 and the bubble radius ranged between 1 nm and 10 nm (see Fig. 11). Although the resolution of the TEM is below 1 nm, the number of bubbles with radius below 2 nm is underestimated due to the limited contrast. This experiment was simulated using the same exposure conditions. The simulated bubbles density c b was found to be in accordance with the one measured experimentally. Some discrepancies were found at the near surface.

The bubble radius ( r b ) is however overestimated by an order of magnitude compared to experimental meas- urements. This could imply that the current model linking the He content to the bubble radius is overestimated and that a more accurate one is needed. Finally, it would be worth investigating this further to determine if some saturation in the bubble size occurs and the impact of initial defects.

Conclusion

A simplified model for He clustering in W based on the existing model by Faney et al. [START_REF] Faney | Spatially dependent cluster dynamics modeling of microstructure evolution in low energy helium irradiated tungsten[END_REF][START_REF] Faney | Spatially dependent cluster dynamics model of He plasma surface interaction in tungsten for fusion relevant conditions[END_REF] was presented. The set of coupled equations of this model was solved using the finite element method. Medium sized clusters were not explicitly accounted for and all immobile clusters were described by the grouped quantity c b . Interactions with self interstitial W atoms produced by trap mutation events were neglected. Similarly, dissociations by vacancy emission have been assumed to be negligible compared to dissociations by He emission.

This model was compared to existing numerical results. It was found that the results produced by the current simplified model were very similar. This allows to identify non-dominant processes that can be safely neglected. Moreover, reducing the number of equations (only eight in this model compared to thousands in Faney's work) greatly reduces the computational cost of the solver. As an example, one simulation point in section "Influence of N" took only a few minutes to run on a non-parallelised solver.

This novel simplified implementation was then applied to a semi-infinite W slab exposed to 100 eV He. The implanted flux was 10 22 m -2 s -1 and the exposure time was 10 s. The maximum bubble size was located in the first few nanometres below the surface which was in very good agreement with experimental observations. The maximum average radius was 15 nm located around 50 nm below the surface.

A parametric study was performed to investigate the influence of both He implantation flux and temperature. Several quantities of interest were computed such as bubbles inventory, total He inventory or average He content in bubbles. Two distinct growth regimes were found namely the nucleation regime in which new bubbles are created and the growth regime in which existing bubbles increase in size.

Comparison with experimental measurements was performed. The bubble density predicted by the model was found to be in good accordance with experimental observations. The predicted bubble radius was however overestimated. It suggests that the model linking He content to cluster radius may need to be adapted for large bubbles.

This model however has a few caveats, some of which have already been discussed by Faney et al. [START_REF] Faney | Spatially dependent cluster dynamics model of He plasma surface interaction in tungsten for fusion relevant conditions[END_REF] . First, the He-to-vacancy ratio of 4:1 in this work differs from MD studies [START_REF] Hammond | Theoretical model of helium bubble growth and density in plasma-facing metals[END_REF][START_REF] Hammond | Large-scale atomistic simulations of low-energy helium implantation into tungsten single crystals[END_REF] . The dilute limit approximation may also not be valid as soon a the volume fraction occupied by He becomes comparable to unity. The volume fraction of He is defined by V = 4/3πr 3 i c i . The diffusion limited reaction model may therefore not be valid for very large bubbles and/or very high concentrations. Then, as observable on Fig. 7c, the He content in bubbles i b keeps increasing with time. One could argue that some physical processes would tend to limit or even avoid this growth at some point. This is particularly true for bubbles very close to the surface that would tend to burst and release their He content to the atmosphere. Implementing a bursting model similar to the one developed by Blondel et al. [START_REF] Blondel | Continuum-scale modeling of helium bubble bursting under plasmaexposed tungsten surfaces[END_REF] could help alleviate this issue. He diffusion may also be limited by mechanical stress induced by pressurised bubbles and He presence in the material, suggesting the current model gives only an upper estimation of the He clusters size in the sample. Similarly, this model does not account for bubbles coalescence which would tend to decrease bubbles concentrations while increasing their sizes. Finally, the approximation made by �r b � ≈ f (�i b �) is valid when the bubble size distribution does not deviate much from i b . Faney et al. suggested that this is the case for large bubbles [START_REF] Faney | Spatially dependent cluster dynamics model of He plasma surface interaction in tungsten for fusion relevant conditions[END_REF] .

In future work, this model will be used to estimate bubble production induced by indirect He production (tritium decay and W transmutation). The use of more accurate correlations for He bubbles growth [START_REF] Hammond | Theoretical model of helium bubble growth and density in plasma-facing metals[END_REF][START_REF] Hammond | Large-scale atomistic simulations of low-energy helium implantation into tungsten single crystals[END_REF] would tend to reduce the overall He content. Additional bubble models such as bursting or coalescence will also be investigated.

Figure 1 .

 1 Figure 1. Representation of He clustering in solids. Dissociation is omitted for simplification purposes. The grey arrows thicknesses represent the magnitude of the reaction rate between mobile He 1 and other clusters at the same distance.

Figure 2 .

 2 Figure 2.He clusters concentration profiles in the tendril at 500 K, 1000 K and 1500 K. Comparison between the current implementation (solid) and Faney's results 30 (dashed). Discrepancies at high temperature are likely due to the use of a different set of dissociation energies. At 1000 K and 1500 K the dashed and solid profiles of He 1 overlap.

Figure 3 .

 3 Figure 3. Concentration profiles of He 1 (left) and bubbles (right) in W exposed to 100 eV He at 10 22 m -2 s -1 and 1000 K.

Figure 4 .

 4 Figure 4. Profile of mean bubble radius r b as a function of depth x in W exposed to 100 eV He at 10 22 m -2 s -1 and 1000 K.

Figure 6 .

 6 Figure 6. Evolution of quantities as a function of the implanted flux and temperature after 1 h of 100 eV He exposure. Grey crosses correspond to simulations points.

( 20 ) 7 .

 207 c He 1 ,ideal = ϕ imp R p D(T) (21) � īb � = N + 1 + a t b Temporal evolution of quantities in W exposed to 100 eV He at 10 22 m -2 s -1 and 1000 K for temperatures varying from 120 K to 1200 K and implanted fluxes varying from 10 17 m -2 s -1 to 10 21 m -2 s -1 . Each line corresponds to a simulation point (grey crosses on Fig. 6a and points on d). The lines are coloured according to the parameter c He 1 ,ideal = ϕ imp R p /D(T) with R p = 1.5 nm and D the diffusion coefficient of He 1 in W.

Figure 8 .

 8 Figure 8. Temporal evolution of � īb � , I bubbles and I in W exposed to 100 eV He at 1.59 × 10 18 m -2 s -1 and 1020 K. The dashed grey vertical line represents the transition from nucleation regime to bubble growth regime.

Figure 9 .

 9 Figure 9. Comparison of several quantities of interest for several values of N in W exposed to 100 eV He at 10 20 m -2 s -1 and 1000 K.

Figure 10 .Figure 11 .

 1011 Figure 10. TEM images of W after exposure to 75 eV He at 2.3 × 10 22 m -2 s -1 and 1053 K for 13 s showing bubbles that have burst, large size bubbles at the near surface and small size bubbles in the bulk.

Scientific Reports | (2021) 11:14681 | https://doi.org/10.1038/s41598-021-93542-9

Acknowledgements

The project leading to this publication has received funding from Excellence Initiative of Aix-Marseille University-A*Midex, a French "Investissements d' Avenir" programme as well as from the French National Research Agency (Grant No. ANR-18-CE05-0012). Work performed under EUROfusion WP PFC. This work has been carried out within the framework of the EUROfusion consortium and has received funding from the Euratom

Author contributions

Competing interests

The authors declare no competing interests.


In all tables, data are weighted by socio-demographic variables and population size. The results for 1999 are not presented in order not to make them too cumbersome.

Religious indifference and atheism: Developments over time and national differences in Europe

Pierre Bréchon Sciences Po Grenoble/Pacte CNRS/University of Grenoble-Alps (France) Email Address: pierre.brechon@iepg.fr Sociology has focused much on the analysis of religions and has left in shadow the study of secularity and irreligion.

It is important to distinguish religious indifference from militant atheism (Quack, Shuh, 2018;Bréchon, Zwilling, 2020).

Indifference: I am not interested in religion, I have no knowledge of any major religion, I have no convictions about them, I am neither for nor against, this religious world is foreign to me, I am unconcerned.

The atheist on the contrary has convictions. He is anti-religious. He thinks that the Gods do not exist, that they are an invention of human people.

What is the evolution of secularity in Europe? Is there an increase of it as the theories of secularization say? We will see that the answer to this question must be nuanced for several reasons: -All the indicators of the secularity-religion dimension do not evolve in the same way, -The evolution of secularity is different in different parts of Europe (the West is much more secular than the East), -It is also different according to the dominant religion that has shaped the culture of each country.

This paper uses data from the European Values Survey from 1990 to 2017, on 22 countries for which data are available for the whole period. It is the best tool for my research aims as there are many indicators about religiosity and secularity and the main dimensions of values for a long period of time.

Institutional secularity can be measured by the not belonging to a religion Very strong differences in the level of non-religious affiliation.

People without affiliation are much more rare in orthodox countries but much more numerous in countries with a Protestant or multi-faith tradition, while the level is average in Catholic countries due to very strong differences (for example between Poland and Italy on the one hand, and France on the other).

In countries with a Protestant tradition, there is a big difference between the Scandinavian countries, where membership of Lutheranism is as much an expression of adherence to the nation as an indicator of religiosity, and the Anglican Great Britain, where many people have no religious affiliation.

Subjective identity: religious, "not religious" people or "convinced atheist" Atheism was rare in 1990, but is multiplied by almost 3 in 30 years. Opposition to religions is thus on the rise overall (except in Orthodox countries).

Countries where atheism is the most developed (in 2017): France: 23%, Sweden: 20%, Spain: 15%, Czech Republic and Slovenia: 14%, Great Britain: 13%, Germany: 12% (but much more in the ex-GDR).These are the countries where anti-religious values are the most developed.

Religious indifference was already strong in 1990. Globally, we do not observe an overall rise in secularity, but rather a shift from indifference to atheism. Secularism is in fact progressing essentially in Protestant and multi-faith environments.

If religion is of little or no importance in one's life, this reflects either an indifference or a rejection of the religious dimension. There has been little change over the past 30 years in the direction of a slight decline in secularity. Above all, it is necessary to underline the differences according to the dominant religion of the country. Secularity is much stronger in Europe with a Protestant or multi-confessional culture. The regression of secularity is very strong in Orthodox countries.

Educating children in faith, a rare concern

In the battery of 11 qualities that can be privileged in the education of children, religious faith is the one that is the least cited, and is falling (everywhere except in orthodox countries), a sign of a low valuation (compared to other dimensions).

Tableau 4. Citing religious faith as a quality to be encouraged among children (%)

... while hopes for a future after death are growing

The index of belief in an extra-mundane future is constructed with 3 questions: believing in an afterlife, heaven and hell. Non-belief has regressed by 8 points in Europe, very strongly in Orthodox countries, quite strongly in Catholic and Protestant countries, stable in multi-faith environments. And hopes for an extra mundane future are also on the rise among young people, who are nevertheless the least religious (Table 7). The disbelief in the future is stronger among the elderly almost everywhere. Young people do not follow the classic paths of religiosity, but at least some of them develop a form of psychoreligious belief in a future after death. Beliefs in God and the future are linked, although 19% of believers in the future only believe in God to a limited extent.

How can we explain the slow growth of secularity?

When a binary logistic regression on the level of secularity/religiosity is computed, many variables seem to contribute to explain it. Age is not the most explainable variable. By far, it is the fact of not having been religiously socialized in one's youth that is the most explanatory of indifference and atheism. There is therefore a strong element of reproduction in secularizing/religious attitudes, with however a slow loss of religiosity from one generation to the next.

Then comes the fact of being a man or a woman, the gap between genders does not seem to decrease, which raises the question: do women have more need for religion? Or do their social roles lead them to be more religious?

Third explanation: the dominant religion of the country, which is also related to the level of development (Inglehart, 2018 and2021). Some religions, more rooted in a culture of duty, of tradition, embodying the cult of the nation, are more resilient. Societies marked by

Protestantism, more open to modernity, have been secularizing longer and more strongly.

The most secularized countries are also generally the richest countries and this often corresponds to the regions most marked by Protestant culture. So economic and educational modernization does seem to be a factor in secularization (Inglehart, 2018 and2021).

Modernization is also characterized by a change in the value system, with the growth of individualization values (valorization of individual autonomy). This evolution of values is more pronounced among young people who are also more secularized. And the more individualized one is, the more secularized one is.

Table 1

 1

	. No religious affiliation (EVS)*
	Tradition	Pays	1981	1990	2008 2017-2018 Average 2017
		Iceland	1	2	8
		Demark	6	8	12
	Protestant	Finlande Norvay	-4	11 11	24 20		53
		Sweden	7	19	37
		Great Britain	in	42	45
	Multi-confessional	Germany The Netherlands Estonia	9 39 -	18 54 87	20 53 56		42
		Poland	-	4	4	8
		Lithuania	-	37	15
		Italy	7	16	16
		Portugal	-	28	18
		Austria	-	15	15		35
	Catholic	Slovakia	-	29	26
		Slovenia	-	27	29
		Spain	9	14	25
		Hungary	-	42	47
		France	26	38	51
		Czech Republic	-	60	71
	Orthodox	Bulgaria Romania	--	64 6	26 2	4	9
	European Mean	-	33	29		38

Table 2 .

 2 Feeling non-religious or convinced atheist (in %)

			1990	2008	2017-2018
	Country of tradition	Not	Convinced	Not	Convinced	Not	Convinced
		religious	atheist	religious	atheist	religious	atheist
	Protestant	33	5	46	8	49	13
	Multi-confessional	69	4	37	7	37	12
	Catholic	31	4	26	9	27	12
	Orthodox	33	3	28	3	18	2
	European Mean	37	4	31	8	32	11
	Source : EVS.

Table 3 .

 3 Religion, little or not at all important in my life(in %)

	Country of tradition	1990	2008	2017
	Protestant	66	68	65
	Multi-confessional	77	65	65
	Catholic	50	54	48
	Orthodox	39	34	23
	European Mean	57	55	53
	Source : EVS

Non-belief in God: a small progression...

  Small global progression of non-belief but still with big differences according to the religious areas: strong rise in non-belief in Protestant and multi-faith areas, overall stability in Catholic countries, fall in Orthodox countries.

	Country of tradition	1990	1999	2008	2017
	Protestant	12	10	8	9
	Multi-confessionnal	17	11	10	10
	Catholic	26	27	19	17
	Orthodox	33	43	44	39
	European Mean	20	23	19	15
	Source : EVS

Table 5 .

 5 Index of non-belief in God (in %)*

	Country of tradition	1990	2008	2017-2018
	Protestant	30	36	41
	Multi-confessional	26	31	38
	Catholic	27	22	26
	Orthodox	18	11	7
	European Mean	27	23	30
	* Scores from 6 to 10 on a distribution from 6 to 26.

Table 6 .

 6 Index of non-belief to an extra-mundane future* by to the main religion

	Country of tradition	1990	2008	2017-2018
	Protestant	58	50	51
	Multi-confessional	50	52	52
	Catholic	51	43	41
	Orthodox	50	46	30
	European Mean	53	45	45

*Neither to an afterlife, nor to heaven, nor to hell. In %.

Table 7 .

 7 Index of non-belief in an extra-worldly future* in 2017 by age and main religion of

	country (%)
	Country of tradition	18-34 years 35-54 years 55 years and +	Mean
	Protestant	50	50	52
	Multi-confessional	44	48	59
	Catholic	38	42	42
	Orthodox	24	28	36
	European Mean	41	45	47
	*Neither to an afterlife, nor to heaven, nor to hell.

Secularity/religiosity index

From the numerous indicators of secularity/religion, we can construct an index with 10 highly correlated indicators (member of a religious organization, attending services at least once a month, feeling religious, praying every week, believing in God, believing in a personal God or life force, finding God important in one's life, believing in life after death, in reincarnation, finding important to encourage religious faith among children) (alpha=0.82). We retain here the followers of the secular pole (0-3 out of 10), which declined a bit from 1990 to 2008 and seems stable since then. Secularism is more frequent among young people, except in the Orthodox milieu where religiosity is more or less equivalent whatever the generation.

Consequential dimension of secularity/religiosity

With secularization, one could believe that religiosity no longer impacts the value system. But this is not the case. Religiosity and secularity correspond to two distinct value systems (Table 10). The secularists value very strongly the autonomy of the individual, the liberalism of morals (deciding on one's sexuality, life and death), new family forms, equality between men and women, distrust of duties and authority, less national pride and a stronger relativization of work. On the other hand, whether one has little or much religiosity, the feeling of happiness is the same. Religious belief and the hope of eternal life do not seem to generate a stronger "positivity".

Secularists are more active in political protest, expressing themselves more through petitions and demonstrations, and they more often share a left-wing orientation. They are more tolerant of deviance and somewhat more accepting of outsiders.

When we look at a table, as similar as possible to the previous one, built with 1990 data, we observe many permanencies in the relations over the last thirty years. Some relationships have certainly weakened, for example in the sense of happiness or the demand for authority. But we also find reinforcements, particularly for the culture of individualization, which is increasing very strongly in irreligious circles, whereas religious people are as resistant as before to this desire for individual autonomy.

Conclusions

The slow rise of secularism is more due to an atheism of conviction than of a simple indifference.

This secularizing climate of contemporary societies in Western Europe does not produce a clear rapprochement of values between secular and religious people. The ideologicalreligious divide persists.

There are many nuances between secular and pure religious people. Individualized recompositions ("I take some and I leave some in the universe of my religion") are very numerous, as we have seen particularly with the development of hopes for a future beyond the world among people who do not adopt the grand narrative of religious institutions. There is therefore a lot of deregulation of religious universes, less and less controlled by the great institutions of belief.

Europe is religiously very composite and more divided than before because of the contrasting evolutions between countries resistant to secularity and very secular countries.

To consider all Europe as highly secularized and to speak of a European exception (Davie, 2002) in a world that would be increasingly religious is dubious. It is especially the west of Europe that is very secularized. The European situation is not very different from the US: too often we forget that secularity has progressed strongly in recent years in this country, especially in certain states (church attendance going from 63% in Mississippi to 23% in Vermont and New Hampshire).

Thank you for your attention.