Naïma Gaudel
email: naima.gaudel@univ-lorraine.fr

Sébastien Kiesgen

De Richter

Effect of vibrations on granular material flows down an inclined plane using DEM simulations

Keywords: Granular material, Vibrations, DEM simulations, Rheology, Inclined plane, Powders conveying

The influence of transverse mechanical vibrations on dense flows of granular material down inclines with a rough bed is analyzed using 3D DEM simulations and compared with experimental results. The vibrations make appear two distinct behaviors: a gravity-driven regime and a vibration-driven regime. In the gravity-driven regime, our results are consistent with previous studies from the literature. The vibrations induce particles velocity fluctuations in the vibrationdriven regime, which generate a granular temperature profile along the depth of the packing related to the velocity profile. The velocity profile is consistent with a creeping flow which makes appear a critical length scale related to a nonlocal rheology. Our results suggest that this length scale emerges from the characteristic damping length of transverse vibration waves which propagate in the bulk of the granular material. These results could be of interest for the optimization of powders conveying.

Introduction

Industries are looking for tests, correlations, or even rheological models to predict the flow of granular matter during transport, in relation with some measurable physical features. It is also of prime importance that correlations and models take into account the impact of vibrations (amplitude, frequency) on the rheological behavior and the structure (rupture, granulation) of granular flows.

Optimizing their transport requires a better understanding of the influence of a large number of parameters (size and morphology of particles, viscosity of the interstitial fluid, volume fraction, intensity of the applied vibrations, geometrical and confinement effects, ...) on the rheological properties of vibrated granular matter. From a fundamental point of view, studying the influence of vibrations on the flowing behavior of granular dispersions represents a fascinating interplay between the jamming transition in dense granular dispersions and rheology.

The flow of granular dispersions in the absence of vibrations has extensively been studied in the literature, both experimentally [START_REF] Pouliquen | Scaling laws in granular flows down rough inclined planes[END_REF][START_REF] Malloggi | Nonlocal effects in sand flows on an inclined plane[END_REF][START_REF] Artoni | Effective wall slip in chutes and channels: experiments and discrete element simulations[END_REF][START_REF] Jop | Initiation of granular surface flows in a narrow channel[END_REF][START_REF] Deboeuf | Flow rule, selfchannelization, and levees in unconfined granular flows[END_REF][START_REF] Midi | On dense granular flows[END_REF][START_REF] Ancey | Dry granular flows down an inclined channel: Experimental investigations on the frictional-collisional regime[END_REF] and numerically [START_REF] Artoni | Scaling laws for the slip velocity in dense granular flows[END_REF][START_REF] Staron | Flow and jamming of a two-dimensional granular bed: Toward a nonlocal rheology?[END_REF][START_REF] Staron | Correlated motion in the bulk of dense granular flows[END_REF][START_REF] Baran | Velocity correlations in dense gravity-driven granular chute flow[END_REF][START_REF] Lois | Numerical tests of constitutive laws for dense granular flows[END_REF][START_REF] Silbert | Granular flow down a rough inclined plane: Transition between thin and thick piles[END_REF][START_REF] Silbert | Granular flow down an inclined plane: Bagnold scaling and rheology[END_REF][START_REF] Hanes | Simulations and physical measurements of glass spheres flowing down a bumpy incline[END_REF]. In these papers, it has been shown that a visco-plastic µ(I) rheology describes very well the flow of such dispersions.

The dynamical friction coefficient µ, and therefore the yield stress, depends on the dimensionless number I which is the ratio of the internal reorganization time of the dispersion to the typical timescale of the flow. When vibrations are applied, the behavior of granular dispersions is drastically altered. In particular, previous studies have shown that the rheology is modified by changing local blockages properties at the grain scale [START_REF] Gaudel | Bulk and local rheology in a dense and vibrated granular suspension[END_REF][START_REF] De Richter | Vibration-induced compaction of granular suspensions[END_REF][START_REF] Hanotin | Viscoelasticity of vibrated granular suspensions[END_REF][START_REF] Capozza | Triggering Frictional Slip by Mechanical Vibrations[END_REF][START_REF] Kiesgen De Richter | Heterogeneous dynamics of a granular pack under vertical tapping[END_REF][START_REF] Richard | Analysis by x-ray microtomography of a granular packing undergoing compaction[END_REF] and therefore the conditions for flow arrest. For instance, the constitutive law for vibrated granular flows down inclines should have a friction coefficient which depends on the shear rate, the dispersion properties, and the vibration intensity through the relative acceleration Γ = A(2πf ) 2 /g; with A and f the amplitude and the frequency of the vibrations, respectively. In particular, vibrations modify the deposits by permitting grains reorganizations close to the flow arrest (see Gaudel et al. [22]).

In this paper, we numerically study how vibrations influence the limit between the jammed and the unjammed regime as well as the apparent rheology µ(I). The Fig. 1 is a hypothetical phase-diagram of what is expected when vibrations are added to a gravitational flow of a granular dispersion tilted from an angle θ. The addition of vibrations make possible flows below the critical angle θ c , defined as the angle below which no flow can occur without vibrations.

The higher the tilt angle, the smaller the value of Γ required to trigger the flow. Contrary to experiments, with DEM simulations we can easily extract the flow profiles in the bulk of the granular material for a wide range of parameters (vibration intensity, inclination angle, particle properties). In this paper, we exhibit (i) how vibrations make appear possible flows below the angle of avalanche, (ii) how vibrations control the deposits, the critical angle and thus the apparent rheology of granular material, and (iii) how vibration induced flows are related to a characteristic length scale which depends on the relative acceleration Γ.

The results are compared with the ones previously obtained experimentally (see Ref. [START_REF] Gaudel | Granular avalanches down inclined and vibrated planes[END_REF]).

Simulation methods

Presentation of the model

In this work, 3D discrete element simulations (DEM) are performed. This model is based on the method developed by Cundall and Strack [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. It consists in a Lagrangian tracking of each particle. The particle-particle and particlewall contacts are modeled using the soft-sphere approach, i.e. the overlapping is allowed to compute the forces. The particle motion is assessed in a deterministic way, i.e. by solving the Newtonian equations of motion. The DEM scheme used in this study is based on the open-source software LIGGGHTS (LAMMPS for Improved General Granular and Granular Heat Transfer Simulations) [START_REF] Goniva | An open source CFD-DEM perspective[END_REF][START_REF] Kloss | Models, algorithms and validation for opensource DEM and CFD-DEM, Progress in Computational Fluid Dynamics[END_REF].

The aim of this paper is to simulate monodisperse glass beads (ρ = 2500 kg/m 3 , d = 6 mm) flows down a rough and vibrated plane (ρ = 2500 kg/m 3 , λ = 18 and 30 mm). In the framework of such flows [START_REF] Artoni | Scaling laws for the slip velocity in dense granular flows[END_REF][START_REF] Baran | Velocity correlations in dense gravity-driven granular chute flow[END_REF][START_REF] Guillard | Lift forces in granular media[END_REF][START_REF] Zhou | An experimental and numerical study of the angle of repose of coarse spheres[END_REF][START_REF] Frankowski | Calibration and validation of DEM rolling and sliding friction coefficients in angle of repose and shear measurements[END_REF], the particle properties chosen for our DEM simulations are the following: Coulomb friction µ = 0.5, Young modulus E = 5 × 10 6 N/m 2 , Poisson ratio ν = 0.45, restitution coefficient e = 0.6, coefficient of friction particle/particle and particle/wall µ pp = µ pw = 0.5, and coefficient of rolling friction particle/particle and particle/wall µ rpp = µ rpw = 0.02. The integration time step is chosen to be (i) a multiple of the period of the vibrations T = 1/f and (ii) much smaller than T . The order of magnitude is δt ∼ 10 -6 s T . Data are saved every X time step, such as X × δt = T /4. Physical quantities presented in this paper are then averaged over the period of the vibrations.

The specific feature of this work is the addition of vibrations in the simulations. For that, the xy-plane is caused to move by oscillations. The oscillations are initiated by choosing a period. This period depends on the chosen frequency f , and an amplitude A according to the y-direction (see Fig. 2). The amplitude range studied in this paper is A = [0.001 to 6.40] d 1 .

A detailed description of the Hertz-Mindlin model [START_REF] Mindlin | Elastic spheres in contact under varying oblique forces[END_REF][START_REF] Di Renzo | Comparison of contact-force models for the simulation of collisions in DEM-based granular flow codes[END_REF] used in this work can be found in the Appendix A.

Geometry and particle properties

A typical snapshot from simulation is presented in Fig. an angle θ with respect to the direction of gravity g. The flow is induced in the x-direction, and can reach a steady state thanks to a rough bottom. This rough bed is created by fixing a thin layer of particles, with the same properties as the bulk ones and a diameter λ. The relative roughness of the studied system is defined as d/λ. These particles are distributed over the base of the system with a volume fraction similar to 45% (see dark particles in Fig. 2). The steady state is defined as proposed in some previous studies [START_REF] Silbert | Granular flow down a rough inclined plane: Transition between thin and thick piles[END_REF][START_REF] Silbert | Granular flow down an inclined plane: Bagnold scaling and rheology[END_REF], i.e. by checking (i) the stationarity of the surface velocity in time and (ii) the uniformity of the velocity profile. Since the simulations employ periodic boundary conditions in the y-direction, they are not influenced by side walls effects that are known to change the profiles significantly [START_REF] Ancey | Dry granular flows down an inclined channel: Experimental investigations on the frictional-collisional regime[END_REF][START_REF] Hanes | Simulations and physical measurements of glass spheres flowing down a bumpy incline[END_REF]. Therefore, the simulations are similar to measurements obtained in the central part of a real chute flow. A summary of the simulations that have been performed can be found in the Appendix B.

Results and discussions

Phase diagram

The phase behavior of granular particles in chute flow is determined by the height of the flow h and the tilt angle θ, as proposed in 2D by Pouliquen and

Renault [START_REF] Pouliquen | Onset of granular flows on an inclined rough surface: dilatancy effects[END_REF]. This study exhibits the existence of a flow line, defined as the boundary between flow and no-flow regions. This flow line is characterized by

Pouliquen [START_REF] Pouliquen | Scaling laws in granular flows down rough inclined planes[END_REF] and later verified by Daerr and Douady [START_REF] Daerr | Two types of avalanche behaviour in granular media[END_REF]. Using empirical argument, Pouliquen suggested the following analytic relation:

tan θ = tan θ 1 + (tan θ 2 -tan θ 1 ) exp - h stop Ld (1)
with θ 1 the angle where h stop (θ) diverges, θ 2 the angle where h stop (θ) vanishes, and L a characteristic length scale. The values of the deposits h stop are determined by keeping the angle constant, and gradually decreasing the height h, i.e. the system size N . The deposit h stop is defined as the height at which the system jams and flow ceases upon reducing h back down: in this case, only a few particles at the free surface rearrange, but this does not lead to a bulk motion [START_REF] Silbert | Granular flow down a rough inclined plane: Transition between thin and thick piles[END_REF][START_REF] Silbert | Granular flow down an inclined plane: Bagnold scaling and rheology[END_REF].

In Fig. 3, the stability diagram is represented in the vibration-free case Γ = 0 and for Γ = 0.21. In the vibration-free case, the simulations data points are well fitted to Eq. 1. The fitting parameters are θ 1 = 20.92°, θ 2 = 33.56°, and L ≈ 4.73. In the following, we define θ 1 in the vibration-free case as θ c the critical angle below which no flow can occur without the addition of an external energy. In Fig. 3(inset), the deposits are plotted as a function of the rescaled angle θ/θ 1 for our DEM simulations and experimental data extracted from Gaudel et al. [START_REF] Gaudel | Granular avalanches down inclined and vibrated planes[END_REF], with θ 1exp = 24.7°, and d/λ = 1/3. The discrepancy concerning the value of θ c may be attributable to the diameter of the particles:

in experiments, beads of hundred microns are humidity-sensitive, increasing the cohesion between them. In our DEM simulations, the cohesion is not taken into account. Except this, results from DEM simulations are in good agreement with experimental data.

The value of h stop decreases with increasing Γ, for a given angle θ (Fig. 3).

The Eq. 1 gives θ 1 = 20.47°, θ 2 = 32°, and L ≈ 2.19. Hence, the minimum angle at which flow occurs decreases when vibration intensity increases. The vibrations affect the system by flowing a particle thickness (h

stop N V -h stop ),
where h stop N V is the thickness of the deposits with no applied vibration and h stop is the thickness of the deposits [START_REF] Gaudel | Granular avalanches down inclined and vibrated planes[END_REF].

The results obtained in DEM simulations are in agreement with those ob-tained in experiments. We thus can define two regimes separated by a critical angle θ c , whose value depends on the properties of the granular system studied.

Above θ c , this is a gravity-driven regime. Below this angle, no flow can occur without adding vibrations: it is a vibration-driven regime. In the following sections, we study in detail depth profiles and propose scaling laws in the two regimes, considering the effect of the vibrations. the boundaries), as predicted in the framework of µ(I) [START_REF] Midi | On dense granular flows[END_REF][START_REF] Artoni | Scaling laws for the slip velocity in dense granular flows[END_REF].

Solid packing fraction. Depth profiles of the solid packing fraction φ are computed for θ = 28°and different values of system sizes N (Fig. 4(c)). The profiles exhibit no significant variation with depth throughout the bulk. The mean value of φ in the bulk is the same for a given θ, whatever the system size, as expected for thick systems [START_REF] Staron | Correlated motion in the bulk of dense granular flows[END_REF][START_REF] Silbert | Granular flow down an inclined plane: Bagnold scaling and rheology[END_REF]. Moreover, the profiles are exactly the same even when vibrations are applied for a given system size N . However, the mean value of φ decreases with the slope (Fig. 4(d)) due to a growing dilatant region near the free surface (i.e. a lower packing fraction), as already observed in previous studies [START_REF] Staron | Correlated motion in the bulk of dense granular flows[END_REF][START_REF] Silbert | Granular flow down an inclined plane: Bagnold scaling and rheology[END_REF][START_REF] Ertaş | Gravity-driven dense granular flows[END_REF]. Our data are well fitted to the relation proposed by Ertas et al. [START_REF] Ertaş | Gravity-driven dense granular flows[END_REF] in the 3D case:

φ(θ) = φ max -c 3D (θ -θ c ) (2)
with the imposed value φ max = 0.595, and the fitted value c 3D = 0.009.

Velocity profiles. Bagnold [START_REF] Bagnold | Experiments on a Gravity-Free Dispersion of Large Solid Spheres in a Newtonian Fluid under Shear[END_REF][START_REF] Bagnold | The Flow of Cohesionless Grains in Fluids[END_REF] proposed a simple 2D model describing the collisional process between particle layers during granular flow. The relations predict a shear stress proportional to the square of the shear rate:

σ ∝ γ2 , with γ = ∂vx(z) ∂z (3)
In a local and Newtonian case, we typically use the relationship σ = η γ, with η a viscosity. If Bagnold's argument is applied to the case of bulk granular flow, with no-slip condition at the base, the velocity depth-profile is such as:

v x (z) ∝ h 3/2 -(h -z) 3/2 (4)
for a system size N of total height h, where z is measured from the bottom of the stack. The velocity profiles are plotted in Fig. 5(a), for various θ, system 145 sizes N , and Γ values. In the following, the flow velocity is denoted v x = v. As expected, all the velocity profiles are well fitted to a Bagnold profile (Eq. 4).

In the µ(I) framework, the main results are retrieved: (i) for a given angle θ and a given Γ value, the velocity increases with increasing the system size N and (ii) for a given system size N , the velocity increases with increasing the tilt 150 angle.

Nevertheless, for a given θ and a given system size N , increasing the Γ value does not change the profile nor the values of the velocity. This observation suggests that the local law µ(I) is not impacted by the vibrations in the steadystate.

Rheology. We present here results obtained for Γ = 0.21 which is the maximum value of the relative acceleration used in this regime in accordance with the experiments. 3 Given this, we have plotted the evolution of the shear stress with z in Fig. 5(b), only for the free-vibration case. The local values of the shear rate are obtained by deriving the raw data of the velocity profile, as defined in Eq. 3.

From these results, the flow curves are plotted for various systems (Fig. 5(c)), allowing us to check the assumption of a Bagnold rheology (Eq. 3).

However, there is a transition away from this Bagnold rheology as the stack becomes thinner. Indeed, the data in the case θ = 28°and N = 12500, i.e. thin pile (h → h stop ), are well fitted to a linear profile. This leads to an increase of the strain rate near the surface. This behavior has already been observed and may be explained by the emergence of nonlocal effects, resulting from competing length scales in these thin piles [START_REF] Malloggi | Nonlocal effects in sand flows on an inclined plane[END_REF][START_REF] Deboeuf | Flow rule, selfchannelization, and levees in unconfined granular flows[END_REF][START_REF] Silbert | Granular flow down an inclined plane: Bagnold scaling and rheology[END_REF][START_REF] Mills | Model for a stationary dense granular flow along an inclined wall[END_REF].

Deposits. To go further, the evolution of the deposits with respect to the Γ value is studied. To properly compare results obtained by DEM simulations and those by experiments, the evolution of the deposits is plotted as a function of Γ/Γ c (Fig. 6), where Γ c is defined as the critical relative acceleration for which h stop vanishes, and for the same value of θ/θ c . These results are in agreement with the observations in experiments [START_REF] Gaudel | Granular avalanches down inclined and vibrated planes[END_REF]. Indeed, adding vibrations in the gravity-driven regime only modifies the values of the deposits, not the dynamics in the steady-state regime.

Velocity scaling

Pouliquen have shown in experiments that granular flows down inclines obey a simple scaling relation [START_REF] Pouliquen | Scaling laws in granular flows down rough inclined planes[END_REF], which was checked in many experimental and nu-merical studies:

Fr = v √ gh ∝ β h h stop (θ) (5)
where β = 0.136 is a fitted coefficient, and the dimensionless Froude number is the ratio of the mean velocity and a velocity of reference √ gh, with g the gravity and h the height of the grains pile.

Our recent study about granular avalanches down rough and vibrated inclines put in evidence that vibrations affect the apparent scaling law, leading to the following relation [START_REF] Gaudel | Granular avalanches down inclined and vibrated planes[END_REF]:

Fr = v √ gh = β h h ef f stop (θ) -1 (6)
where βexp 0.1, and h ef f stop (θ) = h stop N V (θ)×d/λ is an effective stopping height defined as the product of the height of the deposits for a given tilt angle when no vibration are applied, and the relative roughness. In Fig. 7(inset), the ratio of the mean velocity v over the surface velocity v s is plotted as a function of h/h ef f stop , for DEM simulations and experimental results. This ratio is about v/v s 0.6, in agreement with a Bagnold rheology, and falls when h → h stop due to nonlocal effects.

The evolution of v s / √ gh is plotted as a function of h/h ef f stop (Fig. 7) for DEM simulations and experimental results extracted from Gaudel et al. [START_REF] Gaudel | Granular avalanches down inclined and vibrated planes[END_REF].

All the data points collapse on the same straight line which confirms that the influence of the bead size d, the tilt angle θ, and the roughness condition λ, are included in the single function h ef f stop . The data are well fitted to v s / √ gh = 0.18 × h/h ef f stop -1 . Taking into account the value of the ratio v/v s , it leads to the scaling law of Eq. 6, with βsimu = 0.108. This value is in good agreement with the experimental result βexp obtained in the case of vibrated inclines, and also with β determined by Pouliquen.

The scaling law (Eq. 6), derived from the scaling relation proposed by Pouliquen (Eq. 5), appears to be robust to predict flow properties of granular flows, even under vibrations. However, when the gravity is not sufficient to trigger the flow, i.e. in the vibration-driven regime, the apparent scaling law is changed.

Vibration-driven regime

In this section, we are interested in the vibration-driven regime, i.e. θ < θ c

where the flow is induced by vibrations. For the first time, we explore the behavior of granular flows in this regime only for θ = 20°. In the range of 0.02 to 0.21 of Γ values, granular flows occur but velocity profiles are linear.

In order to understand the dynamics due to vibrations in this regime, velocity profiles presented in this section are obtained for Γ values in the range 2.13 to 6.40. For this range of vibrations, no experiment have been performed due to experimental constraints. The following results are thus only obtained by DEM simulations.

Depth profiles

Stress profiles and solid packing fraction. Depth profiles as functions of z/d are computed for the normal stress ratio σ zz /σ xx , the shear stress ratio σ zx /σ zz , and for the solid packing fraction φ (Fig. 8(a) and (b), respectively),

for various system sizes N and Γ values, and for θ = 20°. One can make the same remarks as in the gravity-driven regime (cf Section 3.2.1): the normal stress ratio equals one and the shear stress ratio equals tan(θ), except at the bottom and the surface of the flow (i.e. the boundaries). The depth profiles exhibit no significant variation with depth throughout the bulk, and the mean value of φ in the bulk is the same for a given θ, whatever the system size N . It dismisses the assumption of a plug flow that would be created by the vibrations, otherwise we shall observe an abrupt change in the volume fraction with height, denoting the interface between the flowing and the plug regions [START_REF] Brewster | Plug flow and the breakdown of Bagnold scaling in cohesive granular flows[END_REF][START_REF] Campbell | Chute flows of granular material: some computer simulations[END_REF]. The bulk value of φ keeps constant with the Γ value, and seems to have reached its maximal value φ = 0.59 (Fig. 8(d)), in agreement with the expected value of φ max (Eq. 2).

Rheology. The depth profiles of the shear stress are plotted for θ = 20°and

various system sizes N and Γ values (Fig. 9(a)). The shear stress is linear with the depth and only depends on the system size, as already observed in the regime θ < θ c (cf Fig. 5(b)). However, the flow curves presented in Fig. 9(b) differ from the case gravity-driven regime. Actually, the shear stress seems to be proportional to the square root of the shear rate σ ∝ γ1/2 .

Velocity profiles. This difference comes from a velocity profile which is not of the Bagnold form anymore (Fig. 9(c)).

In general, for a given Γ value, the velocity decreases with increasing the system size N .

And for a given value of N , the velocity increases with increasing Γ. These behavior can be explained by an increase of the particle mobility, all the more important that the energy of vibrations provided to each particle is high.

The velocity profiles exhibit an exponential form, with a fast increase of the velocity near the bottom of the plane before reaching a quasi constant value.

This kind of behavior has already been observed, in flow of upper-convected

Maxwell fluid over a plate [START_REF] Hayat | Series solution for the upper-convected Maxwell fluid over a porous stretching plate[END_REF][START_REF] Hayat | MHD stagnation-point flow of an upperconvected Maxwell fluid over a stretching surface[END_REF], or also in creeping flows down planes [START_REF] Jop | Initiation of granular surface flows in a narrow channel[END_REF][START_REF] Komatsu | Creep Motion in a Granular Pile Exhibiting Steady Surface Flow[END_REF]. Most particularly, this behavior has been put in evidence in flows of a 2D granular bed close to the jamming transition (i.e. for I ≤ 2 × 10 -2 ) by Staron et al. [START_REF] Staron | Flow and jamming of a two-dimensional granular bed: Toward a nonlocal rheology?[END_REF], and the velocity profile is given by:

v(z) √ gd = v s / √ gd 1 -exp -H/d ξ/d 1 -exp - z/d ξ/d (7)
with v s = v(z = H) depending on the vibration intensity and the system size, and ξ a characteristic length.

The evolution of this dimensionless characteristic length ξ/d is plotted as a function of the acceleration of the vibrations Γ for various system sizes N (Fig. 10, blue markers). The value of ξ/d can be considered as a constant for a given Γ value, whatever the system size. However, this value decreases with increasing the intensity of vibrations, and seems to be an inverse function of Γ such as ξ/d ∝ 1/Γ. An interpretation of this characteristic length is proposed in the next section 3.3.2. In this figure, the case N = 37500 is not presented because the data do not collapse on the others. A possible explanation is that, for thin systems, confinement effects are expected to become significant. This would need to be established in further work.

Velocity scaling

Staron et al. [START_REF] Staron | Flow and jamming of a two-dimensional granular bed: Toward a nonlocal rheology?[END_REF] propose that, in the case of an exponential profile for the velocity, the shear rate is actually a function of the granular temperature T .

The granular temperature is defined as the square of the mean particles velocity fluctuations:

T (z) = δv 2 = δv 2 x (z) + δv 2 y (z) + δv 2 z (z) (8)
The rescaled shear rate is plotted as a function of the rescaled granular temperature in Fig. 11, for θ = 20°and various values of N and Γ. The collapse of all data on the same master curve over several orders of magnitude is suitable for a power-law, and the best fit is given by: γ(z)

d g = α × T (z) gd β (9)
with α = 1.33 and β = 2/3. This scaling and the values of α and β are in agreement with previous observations of grains dynamics in granular bed flows [START_REF] Staron | Flow and jamming of a two-dimensional granular bed: Toward a nonlocal rheology?[END_REF], in simulated plan shear flows [START_REF] Xu | Stabilization of nonlinear velocity profiles 23 in athermal systems undergoing planar shear flow[END_REF], and in annular shear configurations [START_REF] Losert | Particle dynamics 475 in sheared granular matter[END_REF].

This result strongly suggests a nonlocal behavior.

This scaling also implies that an exponential decay is a good approximation for the decrease of the temperature along the z-direction. Since the derivative of the Eq. 7 has to be equal to the expression of the shear rate in Eq. 9, the depth temperature profile is given by:

T (z) gd ∝ T 0 gd exp - z βξ (10)
with T 0 a constant depending on the value of Γ. The depth profiles of the 265 granular temperature are plotted in Fig 12 . Except close to the boundaries, the data are well fitted to the above relation (Eq. 10). This confirms the scaling law proposed in Eq. 9.

However, the origin of the higher temperature at the bottom and its exponential decrease is still not clear. In this paper, we make the assumption that the exponential profile would come from the damping of the waves created by the vibrations. In our work, transverse vibrations are imposed to the granular system at the bottom of the plane. We thus propose to draw an analogy with the flow of a viscoelastic fluid close to a plane oscillating transversely. In this case, the velocity parallel to the plane, i.e. v y , verifies: [START_REF] Baran | Velocity correlations in dense gravity-driven granular chute flow[END_REF] with f (z) a function of the depth z, ω a typical frequency, and ϕ the phase offset.

v y (z, t) = f (z) cos(ωt + ϕ) = f (z)e iωt (
The velocity v y is plotted with respect to z at different times for a typical simulation (Fig. 13(a)). For one period of vibration defined as T = 1/f , the mean values of the velocity in time and in space are equal to zero. The parallel velocity propagates along the z-direction with damping oscillations. Indeed, the vibrations propagate in the bulk and are dissipated at each layer of particles.

The function describing the oscillations seems to be of the form of Eq. 11 but there is a change in the period of the oscillations: it decreases with z. An example is represented in Fig. 13(a), and the observation is exactly the same for all our DEM simulations. This may come from a coupled effect with the gravity but further investigations are necessary to obtain this information. In all cases, it is well known that the dissipative mechanism of the system is included in the damping coefficient extracted from the exponentially decaying envelope:

v y (z) √ gd ∝ γ exp -κ z d ( 12
)
with γ a constant depending on the parameters of the simulation, and κ the damping coefficient. A typical example is plotted in Fig. 13 important that Γ is high, and thus unjammed the system.

Conclusions

In order to test the influence of vibrations on free surface granular flows, we have performed particle-based DEM numerical simulations of stationary chute flows, for a large range of vibration intensities and inclination angles. As expected in the phase-diagram (Fig. 1), adding vibrations allows flow under the angle of avalanche and modifies the transition between the jammed and the unjammed regime. We identify two regimes, in agreement with experimental results.

A gravity-driven regime appears for angles larger than the angle of avalanche, where a collisional rheology is recovered at no variation and various vibration intensities. Vibrations modify the depth of the deposits and thus the apparent rheology. All data can be rationalized by introducing an effective deposit length

h ef f stop .
Below the angle of avalanche, stationary flows can be obtained by applying controlled vibrations. In this vibration-driven regime, the Bagnold rheology fails and exponential velocity profiles are obtained. We show that the rheology in this regime is related to an exponential profile of the granular temperature along the depth of the granular packing. A characteristic length scale appears which is directly related to the damping factor of transverse vibrations waves. These waves propagate in the bulk of the packing and make the rheology nonlocal, with γ(z) ∝ T (z) 2 /3 . This scaling is in agreement with results previously obtained in the literature when I → 0, where velocity fluctuations control the rheology.

In a future work, it could be of interest to perform more simulations in order to obtain the relationship between the minimum relative acceleration necessary to trigger the flow and the inclination angle.

F c = F n + F t F n = k n δ n n + γ n v n F t = k t δ t t + γ t v t , if |F t | ≤ µ |F n | = µ |F n | , else
where v n and v t are the normal and the tangential components of the relative velocity of the particles and µ is the friction coefficient. The normal and tangential stiffness k n and k t , and the damping coefficients γ n and γ t are defined in Tab. 1.

Concerning the calculation of the rolling coefficients, a rolling torque resistance model for spherical particles is applied. This model applies a constant directional torque on two particles in contact, depending on their rotational velocity [START_REF] Zhou | An experimental and numerical study of the angle of repose of coarse spheres[END_REF][START_REF] Zhou | Rolling friction in the dynamic simulation of sandpile formation[END_REF][START_REF] Ai | Assessment of rolling resistance models in discrete element simulations[END_REF]:

M µ = µ r d 4 F n ωr |ωr|
where M µ is the frictional torque, µ r is the rolling friction coefficient and ω r is the relative rotating velocity.

Appendix B Summary table presenting the performed simulations

Here is presented a summary of the simulations that have been performed in this paper, with their properties (Tab. 2).

Notes

1 In order to compare DEM simulations and experimental results presented in Ref. [START_REF] Gaudel | Granular avalanches down inclined and vibrated planes[END_REF], a range of vibration amplitudes is selected for DEM simulations, ensuring the same kinetic energy E k by grains in both cases:

E k DEM = E k EXP ⇔ mg DEM (Af ) 2 DEM = mg EXP (Af ) 2

EXP

where mg is the mass of one grain defined as mg = 4/3πR 3 ρ, with R the radius and ρ the density of the particle under consideration. For instance, we used glass beads in experiments Inset: rescaled stability diagram hstop/d -θ/θ 1 , with θ 1 obtained through Eq. 1, for DEM simulations and experimental data extracted from Gaudel et al [START_REF] Gaudel | Granular avalanches down inclined and vibrated planes[END_REF]. θ 1 is the critical angle below which no flow occurs for a fixed value of Γ. All data correspond to a relative roughness d/λ = 1/3.

σ zz /σ xx 0,8 0,9 1 1,1 1,2 z/
d/λ = 1/3, θ = 22 • d/λ = 1/3, θ = 24 • d/λ = 1/3, θ = 28 • d/λ = 1/5, θ = 22 • Experiments
Figure 7: Froude number vs/ √ gh as a function of h/h ef f stop for DEM simulations and experimental data extracted from Gaudel et al [START_REF] Gaudel | Granular avalanches down inclined and vibrated planes[END_REF]. The dashed line corresponds to the fit given by Eq. 6. Inset: Velocity ratio v/vs as a function of h/h ef f stop , for the same data.

  2. The studied system is made of N monodisperse spherical particles (defined as the system size in this paper) of diameter d and mass m flowing down an inclined chute geometry. The simulations volume consists of a rectangular box with a chute base area 2 equals 20d × 100d, a rough bed in the xy-plane, an open top, and periodic conditions in both x and y directions. This simulation cell is tilted by

3. 2 .

 2 Gravity-driven regime 3.2.1. Depth profiles Stress ratio. Depth profiles as functions of z/d are computed for the normal stress ratio σ zz /σ xx (Fig. 4(a)) and the shear stress ratio σ zx = σ zz (Fig. 4(b)) for various inclinations angles and N = 112500. The normal stress ratio is close to one at intermediate values of the incline angle. The normal stress ratio is no longer one for high values of θ because kinetic contribution becomes more important than the frictional one. However, when θ is near the avalanche angle, there are large fluctuations in the contact network because flows are very slow (creep flow). Indeed, large fluctuations in the normal stress appear. The shear stress ratio equals tan(θ), except at the bottom and the surface of the flow (i.e.

  (b). The damping coefficient is extracted for each simulation, i.e. for θ = 20°, all system sizes and Γ values studied. This damping coefficient can be defined as the inverse of a penetration depth of the parallel velocity κ -1 = δ/d, meaning the penetration depth of the vibrations in the bulk. Taking into account that the granular temperature is actually created by the vibrations, it should exist a relation of proportionality between T (z) and v y (z). This directly leads to the fact that δ/d = βξ/d. Hence, the evolution of the extracted penetration depth δ/βd is plotted as a function of the intensity of the vibrations Γ for θ = 20°and various values of N (Fig.10, red markers). The data collapse on the same master curve and δ/d seems to be an inverse function of Γ in the range of values studied here. The waves propagate more efficiently than the intensity of vibrations is high. Moreover, the values of δ/βd extracted from the exponential envelopes are the same than the values of ξ/d extracted from the velocity profiles. These results validate our assumption, and point at the velocity fluctuations due to the vibrations at the bottom as the agent allowing for nonlocality. The vibrations create fluctuations that enable rearrangements at the grain scale, all the more

Figure 1 :

 1 Figure 1: Hypothetical phase-diagram representing the influence of the vibrations on the jamming transition in gravitational granular flows. θc represents the critical angle below which no flow is possible without a minimum acceleration of vibrations Γa.

Figure 2 :Figure 3 :

 23 Figure 2: Typical snapshots from simulation, filled and empty box, (3D and bottom views) with a chute base area equals 20d × 100d, that is inclined at θ to the horizontal. The rough bed is fixed and denoted by dark particles of diameter λ. The flowing particles of diameter d are in light color. Here, d/λ = 1/3, and H/d 25 corresponding to N = 62500 particles.

Figure 4 := 24 •= 24 •Figure 5 :Figure 6 :

 4242456 Figure 4: (a) Profiles of average normal stress ratio for N = 112500 and different values of θ, the dashed line represents σzz = σxx. (b) Profiles of average shear stress ratio for N = 112500 and different values of θ. (c) Depth profiles of the packing fraction for θ = 28°and different values of system sizes N , solid lines are data extracted for Γ = 0 and markers for Γ = 0.21. (d)Tilt dependence of the packing fraction, the solid line is fit to Eq. 2[START_REF] Ertaş | Gravity-driven dense granular flows[END_REF]. All data correspond to a relative roughness d/λ = 1/3.

26 σFigure 8 :

 268 Figure 8: (a) Profiles of average normal stress ratio for θ = 20°and different values of N and Γ, the dashed line represents σzz = σxx. (b) Profiles of average shear stress ratio for θ = 20°a nd different values of N and Γ, the dashed line represents tan(20 • ). (c) Depth profiles of the packing fraction for θ = 20°and different values of N and Γ. (d) Γ dependence of the packing fraction, for θ = 20°and different values of N . All data correspond to a relative roughness d/λ = 1/3. For clarity, legend is the same for all figures.

Figure 9 :

 9 Figure 9: For θ = 20°and for various Γ values and system sizes N , (a) Profiles of the shear stress σxz = σ as functions of the normalized depth z/d. (b) Flow curves σ -γ. The dashed line is a guide for the eyes: σ ∝ γ1/2 . (c) Velocity profiles. Markers are data extracted from DEM simulations and solid lines are fits to exponential form (Eq. 7). All data correspond to a relative roughness d/λ = 1/3. For clarity, legend is the same for all figures.

Figure 10 :Figure 11 :

 1011 Figure 10: Evolution of the characteristic length ξ/d (blue markers) and the penetration depth δ/βd (red markers) with respect to the value of Γ, for θ = 20°and various system sizes N . All data correspond to a relative roughness d/λ = 1/3. The dashed line is a guide for the eyes: ξ/d, δ/d ∝ 1/Γ.

Figure 13 :

 13 Figure 13: Evolution of the transverse velocity vy/ √ gd along z/d for θ = 20°, N = 62500, and Γ = 2.13. (a) At different times t during one simulation, t = 24.5 ms corresponding to one period of vibration. Values of mean velocities in time and space are equal to zero. (b) Superposition of the waves. The solid lines represent fits to exponential envelopes (Eq. 12), here vy/ √ gd = ±0.115 exp (-0.112z/d).

Table 1 :

 1 Expressions of the coefficients of the contact force model used in the DEM simulations. Here, E is the Young modulus, ν the Poisson ratio, e the restitution coefficient, and β = ln(e) √ ln 2 (e)+π 2 the damping ratio.

	Coefficient	Particle/particle	Particle/wall
	k * n = k n / k * t = k t / √ √ δ n δ n γ n /δ 1/4 n γ t /δ 1/4 t	expression E 3(1-ν 2 ) √ d 3(1-ν) 2-ν k * n -5 2 β mk * n -5 3 β mk * t	expression 2E 3(1-ν 2 ) d/2 3(1-ν) 2-ν k * n √ 5β √ 2mk * n -10 -3 β √ 2mk * t

Table 2 :

 2 Summary of the performed simulations with the properties that have been used: the relative roughness d/λ, the inclination angle θ, the system size N and the relative acceleration Γ.

		Gravity-driven regime θ > θ c	Vibration-driven regime θ < θ c
	d/λ θ (°)	N	Γ	d/λ θ (°)	N	Γ
		21	112500-125000 137500-150000	0-0.021 0.21			37500	2.13 4.26
	1/3	22 24	37500-50000 75000-112500 12500-25000-37500 0-0.021-0.043 0-0.021-0.043 0.13-0.21 62500-75000-112500 0.13-0.21	1/3	20	62500	6.40 2.13 4.26 6.40
		28	5000-12500-37500 62500-87500-112500	0-0.021-0.043 0.13-0.21			112500	2.13 4.26
	1/5	22	87500	0-0.13				6.40

Acknowledgments

This study is conducted in the framework of the "PowderReg" project, funded by the European programme Interreg VA GR within the priority axis 4 "Strengthen the competitiveness and the attractiveness of the Grande Région / Großregion". The authors thank F. Qi (University of Luxembourg) for fruitful discussions and for his wise advice.

Appendix A Description of the model

In this paper, the contact model used in the DEM simulations is a Hertzian-Mindlin model, with friction, rolling, and dissipation, but with no cohesion [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF][START_REF] Kloss | Models, algorithms and validation for opensource DEM and CFD-DEM, Progress in Computational Fluid Dynamics[END_REF][START_REF] Mindlin | Elastic spheres in contact under varying oblique forces[END_REF][START_REF] Di Renzo | Comparison of contact-force models for the simulation of collisions in DEM-based granular flow codes[END_REF]. For two grains in contact of diameter d and mass m separated by a distance r, the normal direction n is defined as the line connecting them centers, and the tangential one t is in the plane perpendicular to n. The contact force F c is then zero when r > d. When r < d, one can define the normal overlap δ n = d -r, and the tangential one δ t corresponding to the relative displacement of the grains from the onset of the contact. The force is then decomposed in its normal and tangential parts, F n and F t , respectively:

(ρ = 2500 kg/m 3 ) with a diameter d = 100 µm, submitted to a frequency f = 30 Hz, and a range of amplitude of vibrations A EXP = [0.5 to 4] d, the order of magnitude of the amplitude of vibrations that has to be used in DEM simulations is then A DEM = [0.001 to 0.01] d. 2 The study of the propagation of the vibrations along the z-direction requires a large dimension size in the y-direction, i.e. along the vibrations direction. 3 We have checked that these results remain fair for Γ < 0.21.


Vianney Motte

Dominique Gosset

Thierry Sauvage

Hélène Lecoq

Nathalie Moncoffre

Helium apparent diffusion coefficient and trapping mechanisms in implanted B 4 C boron carbide

Keywords: Boron carbide, helium diffusion, nuclear reaction analysis, transmission electron microscopy 1 Introduction

When boron carbide is irradiated in nuclear reactors, large helium quantities are produced due to the (n,α) neutron absorption reactions. In a recent paper, we have identified the trapping sites for helium, grain boundaries and damaged zones. In this paper, we propose the determination of an apparent diffusion coefficient for helium. 3 He implantations then annealing were performed in B4C samples of different grain sizes. The 3 He( 2 H,α) 1 H nuclear reaction was used to profile helium before and after annealing. Helium profiles with two superimposed components were observed. The narrow component is attributed to helium trapped in the implanted, damaged zone, either in clusters too small to be seen by TEM or as helium-defect complexes. The large component evolution is bounded by the distribution of the grain boundaries surrounding the implanted zone, confirming the high trapping efficiency of the grain boundaries. From an Arrhenius plot, an activation energy for helium diffusion in grains of 2.0 ± 0.2 eV was deduced in a large grain size material.

2 rhombohedra with a most probable C-B-C composition; this leads to a mean composition close to B4C .

Under neutron irradiation at high temperature, amorphisation of boron carbide has never been observed. However, recent studies [START_REF] Gosset | Evidence of amorphisation of B4C boron carbide under slow, heavy ion irradiation[END_REF][START_REF] Gosset | Amorphisation of boron carbide under slow heavy ion irradiation[END_REF][START_REF] Victor | Structural modifications induced by ion irradiation and temperature in boron carbide B4C[END_REF] have shown that amorphisation can occur under ion implantation at and below room temperature when the damage reaches about 8 dpa (displacements per atom). In addition, neutron irradiation produces large quantities of helium due to (n,α) absorption reactions on boron-10. Above 500°C, most observations show helium accumulates as highly pressurised bubbles that lead to swelling and microcracking [START_REF] Jostsons | Defect structure of neutron irradiated boron carbide[END_REF][START_REF] Jostsons | Microstructure of boron carbide after fast neutron irradiation[END_REF][START_REF] Bykov | Effect of high temperature annealings on swelling and helium release in boron carbide irradiated with neutrons[END_REF][START_REF] Shcherbak | Radiation damage in neutron irradiated boron carbide[END_REF][START_REF] Hollenberg | In-reactor measurement of neutron absorber performance[END_REF][START_REF] Hollenberg | Effect of irradiation temperature on the growth of helium bubbles in boron carbide[END_REF][START_REF] Stoto | Swelling and microcracking of boron carbide subjected to fast neutron irradiations[END_REF].

In order to avoid the difficulties associated with neutron irradiation, we have studied the behaviour of helium clusters that form after ion implantation and annealing in implanted boron carbide [START_REF] Motte | Helium clusters nucleation and growth in implanted B4C boron carbide[END_REF].

Increasing post-implantation annealing temperature from 500°C to 1500°C, different clusters form in the implanted zone: nanometric spherical bubbles, lenticular agglomerates of nanometric bubbles, platelets along the (111) or (001) rhombohedral planes, tri-dimension facetted bubbles. Dual-beam implantations, leading to a helium upon atom displacement ratio of about 1/200, similar to fast neutron irradiation [START_REF] Gosset | Evaluation of damage in neutron irradiated boron carbide[END_REF], have been performed. They show a higher cluster density in the implanted zone as compared to the single beam implantation. These observations are in agreement with those performed on material irradiated in a fast neutron reactor [START_REF] You | Helium bubbles and trace of lithium in B4C control rod pellets used in JOYO experimental fast reactor[END_REF]. We have also shown helium clusters only form in the grain boundaries and in the damaged zones (Figure 1). In fact, it appears that helium bubbles nucleate only in damaged or faulted zones. (10 16 cm -², peak concentration 6.10 20 /cm 3 ) at room temperature and annealed at 700°C for 1 h [START_REF] Motte | Helium clusters nucleation and growth in implanted B4C boron carbide[END_REF].

However, helium behaviour in boron carbide has been studied mostly in fast neutron irradiated samples. In these conditions, the effect of the different parameters that could influence the 3 behaviour of the gas can be hardly identified. Therefore, the helium diffusion mechanisms and parameters are not well known. Clayton et al. [START_REF] Clayton | The release of helium from slightly irradiated boron carbide and boron carbide-silicon carbide plates[END_REF] have irradiated high density B4C samples in the BNLR (Brookhaven National Laboratory Research) reactor at low temperature, possibly in a thermal neutron flux. The 10 B burnup is estimated about 0.05 at.% (i.e. about 5x10 19 cm -3 ). They estimated a diffusion coefficient from helium release measurements during post-irradiation annealing. The classical expression D = D0 exp(-Ea/kT) was used, with Ea(activation energy) equal to 1.26 eV and D0 (pre-exponential factor incorporating the diffusion mechanism yield and frequency) equal to 4.4x10 -9 cm².s -1 . Similar analyses were performed by Hollenberg from 92% density boron carbide samples irradiated at low fluence (burnup 5x10 20 cm -3 , i.e. ~0.4 at.%) and high temperatures (630 -950°C) in the EBR-II fast neutron reactor in sealed capsules [START_REF] Hollenberg | Helium diffusion in irradiated boron carbide[END_REF]. It is worth noting that the fast neutron scattering induces a high atom displacement rate. Slightly different diffusion coefficients were derived, with Ea = 1.46 eV and D0 = 1.40x10 -8 cm² s- 1 . But those values correspond to helium diffusion in materials that have been irradiated at different temperatures,, burnups and for different neutron spectra. Moreover, integral methods are used: annealing is performed on as-irradiated samples. The values are then probably biased because of interactions with irradiation damage, helium clusters or grain boundaries. Kovyrshin [START_REF] Kovyrshin | Liberation of helium in the heating of irradiated boron carbide[END_REF] has irradiated boron carbide powders of different grain sizes in a VVR-M reactor (water-cooled reactor, thermal neutron fluences from 1.2x10 19 to 9x10 19 cm -²) at low temperature. In those conditions, the structural damage arises only from the slowing down of the He and Li resulting from the neutron absorption reaction. After thermo-desorption experiments, the author observed different release steps which analysis leads to different activation energies, 0.3±0.1 eV (attributed to grain surface desorption), 1.2±0.2 eV and 2.3±0.3 eV. A nearly complete helium release was obtained after heating at 1050°C. Re-analysing those data with a different physicomathematical model, Svetukhin [START_REF] Svetukhin | Determining the Energy Parametersof the Helium Accumulation Centersin Irradiated Boron Carbide from Thermodesorption Spectra[END_REF] obtained quite different and very low diffusion energies (~0.7 eV). Nevertheless, it seems that the migration energy of helium is much lower than the B or C interstitial ones, about 3 to 4.5 eV [START_REF] You | First-principles investigation of neutron-irradiation induced point defects in B4C, a neutron absorber for sodium-cooled fast nuclear reactors[END_REF]. It can be then inferred that the migration of free helium will occur at much lower temperature than the healing of the structural damage, which arises around 1200 -1500°C [START_REF] Gosset | Evolution of thermo-physical properties and annealing of fast neutron irradiated boron carbide[END_REF].

On the other hand, Schneider et al. have simulated helium behaviour by DFT calculations [START_REF] Schneider | Stability and kinetics of helium interstitials in boron carbide from first principles[END_REF]. They first determined the most likely interstitial sites for helium. Two different sites were found, respectively located around (6 sites "2") and above (2 sites "3") the central chain. The migration of helium was then studied with the NEB (Nudged Elastic Band) method. Three migration paths were identified between the two most likely interstitial sites, respectively between two sites "2" in a given cell (type "A"), between site "2" and site "3" of adjacent cells (type "B") and between two sites "3" of adjacent cells (type "C").From those migration paths, it was deduced that two kinds of migration exist: planar diffusion in the (111) plane of the rhombohedral structure.

-A 3D migration: combination of types B and C and an energy barrier of 2.2 eV.

Very similar results were recently obtained by [START_REF] You | First-principles investigation of neutron-irradiation induced point defects in B4C, a neutron absorber for sodium-cooled fast nuclear reactors[END_REF]. Those values are thus very close to the ones determined by Kovyrshin. However, the issue with the present DFT calculations is that they consider one single helium atom in a B4C cell without neither point defect nor electron charge effect. In irradiated boron carbide, the material has undergone neutron damage leading to a high density of point defects (about 300 defects per created helium in a fast neutron flux [START_REF] Gosset | Evaluation of damage in neutron irradiated boron carbide[END_REF]) and helium atoms that have been created can agglomerate. Hence, it is difficult to compare those results.

In that context, we have launched a parametric study to determine the parameters that influence the behaviour of implanted helium. In this paper, an apparent diffusion coefficient is determined and the diffusion mechanisms of helium are discussed.

Experiments

Experimental procedure

Pure B4C powders of different grain sizes were sintered either by hot pressing or by SPS (Spark Plasma Sintering) leading to small grain size (0.3-1 µm), medium grain size (1-6 µm) and large grain size (20-60 µm) materials. The large grain size material is intended to determine the apparent helium diffusion coefficient. The small and medium size materials allow studying the influence of the microstructure. The two sintering processes lead to cylindrical pellets with nearly fully dense materials (over 98% of the theoretical density: Figure 2). The B4C pellets were then cut and polished with diamond pastes until the material is mirror-polished. As boron carbide is a very hard material, the polishing steps induce a high density of mechanical defects up to depths of about 500 nm that could interfere with the implanted helium. To heal those defects or reduce their density, subsequent annealing at 1600°C (i.e. close the admitted ductile-brittle temperature threshold) for 3 h was performed. (20-60 µm). Densities are over 98%, most of the porosities correspond to grain pulling out during polishing.

Helium ions were then implanted at different fluences in the discs. This allows controlling some key parameters such as the He concentration, the implanted range and depth and the temperature. Moreover, irradiation defects, whose consequence on diffusion could be important, are produced during helium implantation at relatively low concentration as compared to fast neutron irradiation (see here under). Subsequent annealing was then performed. The heating and (initial) cooling rates were 100°C/min. Different samples were used for each implantation and annealing condition. This is required since low temperature mechanisms could influence the high temperature ones: see discussion (e.g., clusters nucleation, helium-defects complexes dissociation…).

The study was carried out using two main techniques:

• NRA (Nuclear Reaction Analysis), which is a well-known ion beam analysis technique widely used to profile helium [START_REF] Sauvage | Profile measurements of helium implanted in UO2sintered pellets by using the 3 He(d, α) 1 H nuclear reaction analysis technique[END_REF][START_REF] Trocellier | IBA studies of helium mobility in nuclear materials revisited[END_REF]. For this purpose, 3 He implantations were performed on the DIADDHEM (CNRS-Orléans, France) facility [START_REF] Chamssedine | DIADDHEM set-up: New IBA facility for studying the helium behavior in nuclear glasses[END_REF].

• TEM (Transmission Electron Microscopy), to observe helium clusters. For these experiments, 4 He implanted B4C samples were prepared at the JANNUS (CEA-Saclay, France) facility [START_REF] Beck | irradiation and radiation effect characterization at the JANNUS-Saclay triple beam facility[END_REF].

A synthesis of the experimental conditions as a function of the probing technique is presented in Table 1. 2.2 Helium profiling by using the 3 He( 2 H,α)H nuclear reaction 3 He in-depth profiling was obtained thanks to the 3 He( 2 H, 4 He) 1 H nuclear reaction. Hence, the 3 He isotope was implanted to simulate the presence of helium in B4C. The 3 He in-depth distribution was deduced from the 4 He energy spectrum. The apparent helium diffusion coefficients in B4C were then determined from the 3 He depth profile evolution as a function of annealing conditions, temperature and time.

The basic principle of this analysis is the following: deuterons impinge the sample and 4 He particles emitted from the nuclear reaction are detected by a PIPS annular detector placed at 178° from the beam direction and at a distance of 104 mm from the sample surface (detector surface 100 mm², detection solid angle 6.27 mSr). For a single detection configuration, the 4 He signal is convoluted with signals inherent to the presence of parasite nuclear reactions such as (d,p) reactions on boron isotopes. To suppress all interferences in the region of interest of the 4 He signal, we used the coincidence technique. Whenever a particle is detected by the 178° detector, the event is recorded if and only if an event is detected by a second detector placed at 0° with a 1200 mm² active surface and at a distance of 60 mm from the rear surface of the sample (detection solid angle 0.27 Sr). For the 3 He( 2 H, 4 He) 1 H nuclear reaction, the 4 He signal is systematically recorded because simultaneously to the 4 He particle emission, the kinematic of the reaction induces a proton is emitted in the opposite direction of 4 He. The proton energy is high enough (16 MeV) to be transmitted through the 500 µm thick B4C sample and then detected by the 0° detector. The logic coincidence gate is adjusted to 0.2 µs. The proton signal emitted from the 10/11 B(d,p) 11/12 B reactions is rejected because the output of the logic gate is meant to be 0. This coincidence configuration allows reducing the background noise by about 10 4 (Figure 3). Counts (Q=1000 µC) coincidence for the 2x10 15 3 He.cm -2 as-implanted sample

Energy (keV)

The 3 He incident energy was 1200 keV and the angle from the sample normal was 60°, in order to broaden the reference helium profile over the depth experimental resolution, this increasing the depth resolution. The projected range estimated from the SRIM software [START_REF] Ziegler | SRIM -The Stopping and Range of Ions in Matter[END_REF] is 1.37 µm with a range straggling of 0.087 µm (peak FWHM = 0.20 µm). The range is a compromise chosen first to allow a significant 3 He profile broadening before having to take into account surface effects and second to facilitate the FIB-sample preparation for the TEM observations. The implantations of all samples of a given material were simultaneously performed in order to ensure that the fluence is the same for a given analysis. The 2 H incident energy for the helium analysis was 800 keV which corresponds to a smooth evolution of the 3 He( 2 H, 4 He) 1 H cross section from the surface to about 3 µm depth as shown inFigure 4. It allows a satisfactory compromise between cross section and depth resolution.

Figure 4: Blue: cross section curve of the 3 He( 2 H, 4 He) 1 H reaction for 800 keV deuterons depending on the depth in B4C. Red: depth resolution for the DIADDHEM CEMHTI facility [START_REF] Chamssedine | DIADDHEM set-up: New IBA facility for studying the helium behavior in nuclear glasses[END_REF].Green: as-implanted 3He profile (arbitrary units, from SRIM).

The depth resolution was calculated using the RESNRA software. At the projected range (1.37 µm), the depth resolution is close to 0.1 µm, to be compared to the full width at half maximum of about 0.20 µm deduced from SRIM calculation. The calculation and fitting process to determine the helium depth resolution and profile is mentioned in a previous paper [START_REF] Martin | Accurate automated nonresonant NRA depth profiling: Application to the low 3He concentration detection in UO2 and SiC[END_REF].

We chose to implant helium at a low fluence (about 2x10 15 at.cm -2 , i.e. about 0.07 He at.% at the maximum of the distribution). This is a good compromise to have a correct detection yield in acceptable experiment duration and to avoid the formation of helium clusters which were observed above this value by TEM [START_REF] Motte | Helium clusters nucleation and growth in implanted B4C boron carbide[END_REF] as shown in Figure 5. Two samples were implanted at room temperature at different fluences then annealed at 800°C for 1h: such conditions allow the formation of helium bubbles both in grain boundaries and in the implanted zone. TEM specimens have then been prepared by FIB. At a fluence of 2x10 15 at.cm -2 (maximum concentration ~10 20 at.cm -3 ), no bubbles are seen in the implanted zone. However, bubbles form in the grain boundaries, as far as 1 µm from the implanted zone, associated with high strain contrasts. At 4x10 15 at.cm -2 (maximum concentration ~2x10 20 at.cm -3 ), small bubbles form in the implanted zone. However, even if bubbles are not observed in the implanted zone at 2x10 15 cm -², it should be noted that too small clusters cannot be observed by TEM. On the other hand, SRIM simulations [START_REF] Ziegler | SRIM -The Stopping and Range of Ions in Matter[END_REF] show that each helium atom creates about 50 atom displacements before stopping, mainly in their immediate vicinity (the damage peak is nearly superimposed to the implantation peak: Figure 6): this could perturb the subsequent helium diffusion. The He analysis by NRA was performed at the DIADDHEM facility [START_REF] Chamssedine | DIADDHEM set-up: New IBA facility for studying the helium behavior in nuclear glasses[END_REF] with a 2x2 mm 2 deuterons beam size and a 60 nA current. The beam charge was modulated to obtain a 4 He signal integral higher than 3000 events, necessary to reduce the helium concentration error bars.

An example of 3 He profile is shown in Figure 6 and compared with SRIM calculations. The analysed depth differs from the calculated one by about 50 nm, which can be considered as the depth resolution of the method. More importantly, the shapes and the widths of the depth profiles are identical.

Results

Results are divided into three parts according to the sample grain sizes (large, medium, and small grain size).

Large grain size (20-60 µm)

In order to determine an apparent helium diffusion coefficient, helium trapping at grain boundaries had to be minimised. For that purpose, B4C samples with large grain sizes were 3 He implanted: in this material, the grain boundary density is about one hundred times lower than in the medium grain size one. Annealing was performed from 600°C to 800°C (precision ± 10°C) for 1 to 10 h. The helium profiles after annealing are compared to the as-implanted one in Figure 7. Figure 7 shows that after annealing the profiles have broadened, which characterises a diffusion process. After heat treatments, the profiles have a complex shape that can be decomposed into two

Gaussian curves (e.g. Figure 8) centred at nearly the same depth, giving two distinct contributions:

-A narrow one located at the same depth as the implanted zone, we call Ft in the following. Its width is close to the one in the as-implanted sample, we then assume it corresponds to helium trapped in the implanted zone. In all Gaussian adjustments, the same width as in the non-annealed sample was enforced for this component.

-A broad one we call Fm in the following, we assume corresponding to mobile helium that actually diffused in the material. In order to analyse helium diffusion, we then considered the sole mobile fraction. In such a large grain material, the helium profiles remain much narrower than the grain size. We then assumed that the grain boundary trapping is negligible so that the helium volume diffusion is not perturbed. In that case, an apparent diffusion coefficient can be derived, first assuming it does not depend on the crystal structure. To calculate a diffusion coefficient, the width of the mobile fraction for each profile was determined from the Gaussian fittings. A consequence of the Fick's law for Gaussian distributions allows writing:

2 • • (Equation 1)
Where σT is the standard deviation of the Gaussian corresponding to the mobile fraction after a heat treatment at a temperature T and a duration t. σ0 is the standard deviation of the as-implanted profile before annealing (the two σ then include the experimental resolution). DT is the apparent diffusion coefficient at a temperature T and is defined by:

• exp • (Equation 2)
With D0 the pre-exponential factor; Ea the activation energy; k the Boltzmann constant (8.617x10 -5 eV.K -1 ).

DT was calculated for each temperature and duration of annealing (Table 2) and the values were reported into an Arrhenius diagram (Figure 9). In Figure 9, the dots appear relatively aligned; we then deduce that only one diffusion mechanism is active in this temperature range. We then obtain:

12 -D0 = 9x10 -4 ± 3x10 -4 cm².s -1 .

-Ea = 2.0 ± 0.2 eV.

Medium grain size (1-6 µm)

B4C samples were 3 He implanted at the same fluence, 2x10 15 at.cm -2 (Table 1). The samples were then annealed from 600 to 1000°C for 1 h and analysed by NRA. Helium profiles are shown in Figure 10. We observe that the profile areas are nearly equal whatever the thermal treatment, meaning that most helium remained near the implanted zone. As mentioned above, this means that a longerrange diffusion is negligible. and annealed up to 1000°C for 1 h (one different sample for each annealing and NRA analysis).

As observed in the large grain size material, the helium profiles broaden from temperatures above 700°C, showing that helium begins to significantly diffuse. However, for temperatures above 800°C

and annealing duration of 1 h, the helium profiles do no evolve anymore. Figure 11a gives a better look for high temperature annealing and Figure 11b shows helium profiles in implanted B4C samples annealed at 800°C for more than 1 h. Note in Figure 11 In Figure 11a, the profiles corresponding to samples annealed at 900°C and 1000°C for 1 h are similar whereas from our calculated diffusion coefficient, the '1000°C-1h' profile should be about twice as wide as the '900°C-1h' one. Figure 11b shows that the samples annealed at 800°C show similar profiles whatever the annealing time between 1 and 10 h. These results clearly show that helium is trapped. Indeed, we have shown in [START_REF] Motte | Helium clusters nucleation and growth in implanted B4C boron carbide[END_REF] that grain boundaries are very efficient traps for helium and there was no helium diffusion along and across the grain boundaries. This trapping effect could explain why helium profiles remain unchanged for longer time annealing at high temperatures. To check this point, B4C samples with small grain sizes were studied.

Small grain size (0.3-1 µm)

B4C samples with small grain sizes (0.3-1 µm) were implanted with 500 keV 4 He ions at a higher fluence of about 10 16 at.cm -2 (Table 1) in order to allow the formation of clusters. Figure 12 shows TEM images of an implanted sample after a heat treatment at 800°C / 1 h. Bubbles appear in the grain boundaries (Figure 12a), in agreement with results obtained in [START_REF] Motte | Helium clusters nucleation and growth in implanted B4C boron carbide[END_REF].

Grain boundaries were distinguished between two types: with bubbles (delimited in green) and without bubbles (in red). Figure 12b presents a triple point where bubbles are visible in the green grain boundaries and not in the red ones. In Figure 12a, we observe that the grain boundaries that contain helium bubbles all surround the implanted zone. This shows that, in those annealing conditions, there is no helium bubble formation beyond the grain boundaries surrounding the implanted zone.

Discussion

Diffusion at grain boundaries

From the results presented in sections 3.2 and 3.3, it can be concluded that there is no diffusion in and across the grain boundaries, at least up to 1000°C. As already emphasized in [START_REF] Motte | Helium clusters nucleation and growth in implanted B4C boron carbide[END_REF], this is an important result, which cannot be inferred from experiments on specimens irradiated in reactor. In this latter case, helium is homogeneously produced in the whole material. The trapping effect at the grain boundaries manifests only through the existence of denuded zones in their vicinity, about 30 -50 nm wide, without helium bubbles. However, evidence of these zones does not mean there is no diffusion through or along grain boundaries. To go further, we analysed helium implanted in the HD20 small grain B4C after annealing at 800°C for 1 h (in such annealing conditions, we have shown that the helium profile is stabilized in the medium grain size material: Figure 11). In Figure 13, the obtained He profile is compared with the profile obtained in the large grain size material (ERM) in the same annealing conditions. Superimposed is the distribution of the grain boundaries circling or crossing the implanted zone for the HD20 material (green curve), as deduced from its log-normal grain size distribution.

The HD20 helium profile looks clearly narrower than the ERM one. This means that helium cannot diffuse freely. Here again, the profile looks like constituted by two contributions, the narrow one is located at the depth of the implanted zone. In order to interpret this complex profile; we have considered the distribution of the sole grain boundaries circling or crossing the implanted zone (with a 100 nm width). The metallographic observations of the hot-pressed HD20 material lead to a lognormal grain size distribution (logarithmic variance and median respectively σ = -0.2; µ = 0.6). From this distribution, we have calculated the density of grain boundaries circling the implanted zone as a function of the distance to the implanted zone or intersecting it. This leads to a bell-shaped curve (figure 13) which closely superimposes to the HD20 helium distribution, excepted in the implanted zone. The distribution of the grain boundaries circling the implanted zone reproduces the broad component of the helium profile; the density of grain boundaries intersecting the implanted zone reproduces only partially the narrow component of the helium profile. It follows that at the helium fluence here chosen, helium trapping (e.g. in clusters) in the implanted zone or into the grains is actually low (see section 4.2): there is almost no helium outside the grain boundaries. Finally, due to a much higher density of grain boundaries in the HD20 material, the trapped fraction in the implanted zone is much higher than in the ERM material. Thus, we can assume that, as soon as annealing has been processed for sufficient time and temperature conditions (see also Figure 11), the low-concentration implanted helium profiles in boron carbide with small or medium grain sizes fit the distribution of the grain boundaries circling the implanted zone. There however remains a small fraction trapped in the implanted zone. This confirms helium does not diffuse through or in the grain boundaries in the temperature range we consider here.

In our previous papers [START_REF] Motte | Helium clusters nucleation and growth in implanted B4C boron carbide[END_REF][START_REF] Motte | Helium behaviour in implanted boron carbide[END_REF] helium behaviour at high temperatures was followed by TEM. After annealing at 1200°C, small helium bubbles were observed near triple grain boundaries, pointing outside the implanted zone, as shown in Figure 14. This means that in this temperature range, helium can either diffuse in the grain boundaries or evaporate from the bubbles into the material and then re-condense in the vicinity of the grain boundaries. In addition, at 1500°C, the grain boundaries can be opened, this meaning growth and interconnection of the intergranular bubbles.

Helium distribution

The 3 He implantations were performed at helium concentrations close but lower than the bubble formation threshold determined by TEM [START_REF] Motte | Helium clusters nucleation and growth in implanted B4C boron carbide[END_REF], in order to favour the mobile fraction. However, the helium distributions present two superimposed Gaussian-like components, both centred close to the implantation maximum. For the low and medium grain size materials, the so-called trapped component mainly corresponds to helium precipitated into the grain boundaries, as calculated in section 4.1, leading to the He bubbles free zone observed by TEM (denuded zones), about 30 to 50 nm large. Nevertheless, for the large grain size material, the grain boundary density is about 100 times lower than in the mean grain size material; this trapped component should then be considered as negligible. As a result, we have to assume that in the implanted zone, a small fraction of helium is trapped even if bubbles are still not visible. Apart from few bubbles, the trapped species cannot here be identified. Either they form helium clusters too small to be observed by TEM, or they form low mobile helium-defect complexes since about 40 atoms are displaced by each implanted helium atom (SRIM estimation), with the damage profile nearly superimposed to the helium profile (Figure 6).

Diffusion coefficient and activation energy

As mentioned in section 3.2, the helium concentration had to be high enough so that we could determine the profiles with a good resolution in a reasonable experimental time. A 2x10 15 cm -² fluence was chosen, leading to a maximum implanted concentration of about 10 20 cm -3 . This value appears close to the cluster formation threshold. In Figure 5, no helium clusters are visible. As mentioned above, this however cannot preclude the formation of helium-defect complexes. On the other hand, in such a rhombohedral structure, an anisotropic diffusion cannot be excluded (see herebelow), that the sintered, non-textured microstructure of the samples would hide.

In a preliminary work [START_REF] Motte | Helium behaviour in implanted boron carbide[END_REF], an apparent diffusion coefficient had already been estimated and activation energy of 0.52 eV was derived, strongly different from the present result. This can be explained by the depth resolution and the yield of this seminal experiment, too low to allow evidencing the two components we observed here.

The apparent diffusion coefficients estimated in the present study are quite different (about 2 decades higher) from the ones published by Clayton [START_REF] Clayton | The release of helium from slightly irradiated boron carbide and boron carbide-silicon carbide plates[END_REF] and Hollenberg [START_REF] Hollenberg | Helium diffusion in irradiated boron carbide[END_REF] (Figure 15). As mentioned in the introduction, such a difference can tentatively be explained. Clayton and Hollenberg performed helium desorption experiments on materials irradiated in reactor. As a result, the diffusion coefficient deduced from those experiments possibly recovers different mechanisms such as helium-defect complexes migration, trapping and detrapping from the helium bubbles, interactions with grain boundaries, surface desorption, porosity degassing, and so on. Moreover, helium diffusion occurs in a material with a modified composition (boron burning, lithium formation) and with different defect configurations (defect creation by neutron scattering in a fast neutron spectrum is about six times higher than by the sole He slowing-down [START_REF] Gosset | Evaluation of damage in neutron irradiated boron carbide[END_REF]). The D values determined in the present study are certainly less biased by such effects. This suggests the lower the defects density the faster the helium diffusion. This could be checked by performing infusion then desorption experiments.

Figure 15. Helium diffusion coefficients in boron carbide. Clayton [START_REF] Clayton | The release of helium from slightly irradiated boron carbide and boron carbide-silicon carbide plates[END_REF] and Hollenberg [START_REF] Hollenberg | Helium diffusion in irradiated boron carbide[END_REF]: helium produced by neutron irradiation, analyses performed by desorption. Motte (this study): implanted helium, NRA profiling.

The ab -initio calculations performed by A. Schneider and co-workers lead to two activation energies:

1.2 eV for a 2D diffusion along the (111) planes and 2.2 eV for a 3D diffusion [START_REF] Schneider | Stability and kinetics of helium interstitials in boron carbide from first principles[END_REF]. Those values are very close to those obtained by Kovyrshin [START_REF] Kovyrshin | Liberation of helium in the heating of irradiated boron carbide[END_REF] who irradiated boron carbide powders and who identified two release steps by thermo-desorption. The activation energy we determined here is close to the higher one; that is for a 3D diffusion. EBSD (Electron Back-Scattered Diffraction) observations (not presented in this paper, see [START_REF] Motte | Behavior of helium implanted in B4C boron carbide[END_REF]) have evidenced that the materials of the present study are nearly isotropic, with no noticeable texture. Figure 16 illustrates the difference between the actual diffusion length along the (111) direction and the apparent one towards the surface or the bulk. Integrating this effect along all possible (111) orientations in the material would lower the apparent diffusion lengths and distort the helium profile. This does lead to lower the apparent activation energy, but still far from the 1.2 eV value calculated by Schneider et al. A 2D diffusion, associated to low activation energy, should then be visible, for example with large wings on both sides of the diffusion profiles. Part of them could be hidden in the background noise.

However, as mentioned above, the He content measurement by NRA corresponds to the implanted fluence, so this masked long-range diffusion should be very weak. Mechanisms able to inhibit this 2D diffusion should be then identified. As already mentioned, helium implantation creates a high density of defects, very close to the implanted helium zone (cf. Figure 6; note that this defects to helium ratio would be similar in a thermal neutron flux but about 6 times higher in a fast neutron flux [START_REF] Gosset | Evaluation of damage in neutron irradiated boron carbide[END_REF]).

They could prevent migration along the (111) planes or lead to the formation of helium -defect complexes with a higher diffusion energy similar to the 3D value. At last, the micrographs of the samples show twin boundaries: in that case, the free diffusion along the (111) planes cannot take place anymore. Either He is trapped at the twin boundaries or the diffusion paths make zig-zags.

Those defects are always present in B4C crystals, which density depends on the elaboration process and the stoichiometry [START_REF] Xie | Atomic-Level Understanding of "Asymmetric Twins" in Boron Carbide[END_REF]: this effect should be analysed. On the other hand, in-progress studies show that most defects induced by ion implantation are not stable at 500°C [START_REF] Gosset | Amorphisation of boron carbide under slow heavy ion irradiation[END_REF][START_REF] Victor | Structural modifications induced by ion irradiation and temperature in boron carbide B4C[END_REF]. The role of the implantation defects could then be evaluated by performing the 3 He implantation at 500°C (note that the profile broadening would be negligible).

Diffusion lengths

From the apparent diffusion coefficient values, it is possible to estimate the diffusion lengths for different annealing temperatures and durations (Figure 17). The diffusion lengths are very low at 500°C. At higher temperatures, the diffusion lengths are much larger. For example, after a 800°C/10h annealing, helium can diffuse up to 5 µm, which is larger than the medium grain size of this study, thus large enough for helium to encounter grain boundaries and create the bubble walls as observed in [START_REF] Motte | Helium clusters nucleation and growth in implanted B4C boron carbide[END_REF][START_REF] Stoto | Behaviour of implanted helium in boron carbide in the temperature range 750 to 1720°C[END_REF]. Extrapolating those values to the irradiation performed in fast neutron reactors (that is for about one year in the temperature range 500 -1000°C) lead to sub-millimetric diffusion lengths.

Those distances are comparable to the separation lengths of the microcracks forming under irradiation. Free helium should then escape through the microcrack network as soon as it forms: this release acceleration is clearly observed in [START_REF] Hollenberg | Helium diffusion in irradiated boron carbide[END_REF]. However, the retention rates below this cracking threshold (burnup about 5x10 21 neutron captures per cm 3 ) are high [START_REF] Hollenberg | In-reactor measurement of neutron absorber performance[END_REF], meaning the trapping mechanisms (helium bubbles in the grain boundaries and in the grains) are very efficient.

Conclusion

Results presented in this paper contribute to a better understanding of helium diffusion in B4C and of the influence of the initial microstructure. They are complementary of results dealing with helium cluster nucleation and growth [START_REF] Motte | Helium clusters nucleation and growth in implanted B4C boron carbide[END_REF]. Moreover, all those results are obtained in materials with a relatively low defect concentration: this allows a comparison with neutron irradiated materials then an estimation of the induced specific effects, such as helium clusters, composition variation (boron burning and lithium formation) or ballistic damage induced by fast neutron scattering.

NRA analyses of helium-implanted then annealed B4C samples were performed. It was shown that the helium profiles obtained in B4C with a large grain size broadened with increasing temperature and annealing duration, meaning that diffusion occurred. However, in the case of medium grain size samples, above 800°C, or for durations longer than 1h at those same temperatures, the profile widths remain constant. At last, comparing the grain size distribution with the implanted helium profiles of boron carbide with small grain size, allowed us concluding that in those conditions helium is trapped in the grain boundaries. There is no helium diffusion along and across the grain boundaries for temperatures at least up to 1000°C, in agreement with the assumptions we made in [START_REF] Motte | Helium clusters nucleation and growth in implanted B4C boron carbide[END_REF].

To minimize grain boundary trapping, helium was implanted in large grain size B4C. From the evolution of the profile widths between 600 and 800°C, an apparent diffusion coefficient was evaluated, with activation energy of 1.96 ± 0.2 eV. Unlike ab-initio calculations, only one diffusion mechanism is observed, corresponding to the ab-initio isotropic diffusion: the low energy, anisotropic diffusion along the (111) planes is not observed, possibly hidden by the formation of helium-defect complexes or interactions with the twin boundaries.

Further studies, such as desorption or infusion experiments could be performed to improve the understanding of the diffusion mechanisms and obtain a more accurate value of the helium apparent diffusion coefficient. Such experiments should be performed on monocrystals, which are unfortunately not available. Ab-initio calculations are also in progress to study the behaviour of several helium atoms in a damaged B4C matrix and the influence of charged point defects [START_REF] Gillet | The influence of irradiation induced vacancies on the mobility of helium in boron carbide[END_REF].

Figure 1 :

 1 Figure 1: TEM micrograph of helium implanted boron carbide (10 16 cm -², peak concentration

Figure 2 :

 2 Figure 2: FSE ESM micrographs of boron carbide sintered with different powders. a) small grain size (0.3-1 µm); b) medium grain size(1-6 µm); c) large grain size(20-60 µm). Densities are over 98%, most

7 Figure 3 .

 73 Figure 3. Signals recorded by the 178° annular PIPS detector with (blue line) and without (black line)

Figure 5 .

 5 Figure 5. TEM observations of large grain size B4C implanted with 3 He (implantation depth 1.37 µm, surface on the left) then annealed at 800°C 1h. a: 2x1015 at.cm -² (by chance there were grain boundaries in the FIB sample, here decorated with He bubbles). b: 4x1015 at.cm -²: small bubbles are visible in the implanted band.

Figure 6 : 3

 63 Figure 6: 3 He depth profile in boron carbide ( 1200 keV, 1.6x10 15 at.cm -2 , room temperature) ) In blue, the experimental profile with step widths corresponding to the NRA depth resolution. In red, the He profile from SRIM calculations. In dashed green, the damage profile.

Figure 7 :

 7 Figure 7: Helium distribution profiles in B4C (large grains) implanted with 1200 keV 3 He at RT and 2x10 15 at.cm -2 and annealed from 600°C to 800°C for 1 to 10h (width of histogram bars: depth resolution; a different sample for each profile).

Figure 8 :

 8 Figure 8: Helium distribution profiles in B4C (large grains) implanted with 1200 keV 3 He at RT and 2x10 15 at.cm -2 and annealed at 800°C for 1h. The profile after annealing is decomposed into two Gaussian curves: in green the "mobile" fraction; in red the "trapped" fraction (horizontal error bars: depth resolution).

Figure 10 :

 10 Figure 10: Helium distribution profiles in B4C implanted with 1200 keV 3 He at RT and 2x10 15 at.cm -2

  the two "as implanted" profiles correspond to two different samples. The slight difference we observe corresponds to the experimental reproducibility.

Figure 11 :

 11 Figure 11: Helium distribution profiles in B4C implanted with 1200 keV 3 He at RT and 2x10 15 at.cm -2 . a: annealed at 900 and 1000°C for 1 h. b: annealed at 800°C for 1, 3 and 10 h. One different sample for each profile.

Figure 12 :

 12 Figure 12: TEM micrographs of B4C (small grain size) implanted with 500 keV He ions at RT with a fluence of 8x10 15 at.cm -2 and annealed at 800°C / 1 h. Orange: implanted zone. Green: grain boundaries that contain bubbles. Red: grain boundaries that do not contain any bubble. a: general view. b: zoom of the centre of a.

Figure 13 :

 13 Figure 13: Helium profiles (step-wise curves) obtained with large grain size (ERM, blue curve) and small grain size (HD20, orange curve) materials after annealing at 800°C for 1h (fluence 2.10 15 /cm²).

Figure 14 .

 14 Figure 14. Helium bubbles in grain boundaries pointing outside the implanted zone (mean grain size material). After annealing at 1200°C for 1 h, bubbles can be seen in the bottom left grain boundary, while the implanted zone is towards the upper right corner (arrow).

Figure 16 :

 16 Figure 16: 2D diffusion along the (111) planes. DL: diffusion length along the plane for given annealing conditions. ADL: apparent diffusion length from the helium layer.

Figure 17 :

 17 Figure 17: Diffusion lengths of helium in boron carbide estimated from the apparent diffusion coefficient, depending on the temperature and duration of annealing.

Table 1 .

 1 Summary of the experiments.

	Initial powder		ERM	HP	HD20
	Grain size µm		20 -50	1 -6	0.3 -1.0
	Analysis	NRA	TEM	NRA	NRA	TEM
	Objective	Apparent	Fluence	He profile vs microstructure	He bubbles
		diffusion	threshold for	and annealing conditions	distribution
		coefficient	bubbles in
			implanted zone
	Ion	3 He	3 He	3 He	3 He	4 He
	Energy (keV)	1200	1200	1200	1200	500
	Ion beam incidence (°)	60	60	60	60	15
	Range (µm)	1.37	1.37	1.37	1.37	1.21
	Straggling (µm)	0.073	0.073	0.073	0.073	0.054*
	Fluence /10 15 .cm²	2	2 -4	2	2	8
	Max. He c° /10 20 .cm 3	1	1 -2	1	1	5.5
	at.%	0.07	0.07 -0.14	0.07	0.07	0.40
	Annealing °C	600 -800	800	600 -1000	800	800
	h	1 -3	1	1 -10	1	1

* the different stragglings come from different implantation angles. This results in different fluence over maximum concentration ratio.

Table 2 .

 2 Full widths at half-maximum of the helium profiles and apparent diffusion coefficients as a function of annealing temperature and time.

	Annealing	Annealing	FWHM	Apparent
	temperature	time	(µm)	diffusion
	(°C)	(h)		coefficient
				(10 -15 cm²/s)
	As-implanted	-	0.221	-
	600	10	0.419	3.2
	700	1	0.675	102
	700	3	0.884	61
	750	1	1.093	287
	800	1	1.146	317

Figure 9: Arrhenius diagram of the apparent diffusion coefficient of helium in boron carbide (large grains)

Acknowledgments

The experiments at the CEMHTI laboratory have been performed in the frame of the CEA-CNRS NEEDS collaboration contract. We are very grateful to our colleagues of IRCER (Limoges University, France) who prepared the small grain size material by hot-pressing. The FIB equipment was partly financed thanks to the French GENESIS ANR project.

7 References