B E Peigney

O Larroche

V Tikhonchuk

Fokker Planck kinetic modeling of suprathermal α particles in a fusion plasma

Keywords: Fokker-Planck equation, fusion reactions, kinetic effects, inertial confinement fusion plasma, suprathermal particles, multi-scale coupling, explicit schemes

We present an ion kinetic model describing the ignition and burn of the deuterium-tritium fuel of inertial fusion targets. The analysis of the underlying physical model enables us to develop efficient numerical methods to simulate the creation, transport and collisional relaxation of fusion reaction products (α-particles) at a kinetic level. A twoenergy-scale approach leads to a self-consistent modeling of the coupling between suprathermal α-particles and the thermal bulk of the imploding plasma. This method provides an accurate numerical treatment of energy deposition and transport processes involving suprathermal particles. The numerical tools presented here are validated against known analytical results. This enables us to investigate the potential role of ion kinetic effects on the physics of ignition and thermonuclear burn in inertial confinement fusion schemes.

Purpose of the study

Inertial confinement fusion (ICF) is a process of energy production obtained from the nuclear fusion reaction between deuterium (D) and tritium (T) ions. It is a promising and abundant energy source for future power plants. The fusion reactions D + T → α + n + 17.56 MeV take place in a hot and dense plasma compressed and heated by intense laser radiation. The thermonuclear burn of the deuterium-tritium (DT) fuel is supported by energetic αparticles, which are created by fusion reactions at the energy 3.52 MeV. Those suprathermal particles subsequently transfer their energy to the fresh fuel through Coulomb collisions.

In the case of Inertial Confinement Fusion [START_REF] Lindl | Inertial Confinement Fusion -The quest for ignition and energy gain using indirect drive[END_REF][START_REF] Atzeni | The physics of inertial fusion[END_REF], a spherical DT shell is compressed to densities of the order of a few hundred g/cc by the ablation pressure. Fusion reactions start in a central zone characterized by a density ρ ∼ 50 g.cm -3 and a high "ignition" temperature T ≈ 7 -10 keV. The surrounding shell is 10 times colder than the hot spot (T ≈ 0.7 keV). The density of the central "hot spot" is such that the mean free path λ α of fast α-particles is roughly equal to the hot spot radius R [START_REF] Fraley | [END_REF]. This allows the self-heating of the hot spot fuel which serves as a spark that subsequently burns the surrounding colder and denser shell.

The design of ICF targets and the interpretation of ICF experiments rely on numerical simulations based on hydrodynamic Lagrangian codes where kinetic effects are only considered as corrections included in the transport coefficients [START_REF] Lindl | Inertial Confinement Fusion -The quest for ignition and energy gain using indirect drive[END_REF][START_REF] Atzeni | The physics of inertial fusion[END_REF]. The fluid description is relevant if the mean free path of plasma particles, namely electrons and ions, is smaller than the characteristic length scale. Although this condition is reasonably fulfilled during the implosion stage, it does not apply to fast particles, in particular to fusion products near the ignition threshold. Thus, an accurate kinetic modeling is required.

The purpose of the present work is to propose an ion-kinetic description of suprathermal fusion products, treated self-consistently with the ion-kinetic modeling of the thermal imploding plasma. The difficulty lies in the coupling of ion populations characterized by two different energy scales:

• Thermal particles D,T , which form the bulk of the imploding plasma and whose kinetic energy is in the keV range.

• Suprathermal α-particles, created at 3.52 MeV by fusion reactions.

Such a strong disparity in energy scales makes it difficult to build viable kinetic models of fusion reactions. Existing ion kinetic codes can describe the implosion of DT targets in sub-ignition conditions [4,5,6], but the energy release from the fusion reactions is not accounted for in a self-consistent manner. Several simplified methods compatible with hydrodynamic codes have been developed. Haldy and Ligou [8] apply the moment method to model ion energy deposition in a hot and dense homogeneous plasma, but only a stationary case has been considered. A variety of methods based on diffusion models applied to charged-particle transport problems have also been developed. Those methods are of considerable interest, since results on energy deposition profiles can be obtained with a low computational effort. Nevertheless, diffusion methods rely on the assumption that the fast particle mean free path is smaller than the characteristic scale length of the energy deposition zone. This hypothesis does not hold for a typical ICF target near ignition. Corman et al [9] derive a multi-group diffusion model from the Fokker-Planck equation to describe fast ion transport in a fusion plasma. However, they introduce heuristically a flux limiter in order to prevent unphysical behavior when particle flux approaches the free-streaming limit. Pomraning [10] develops a more sophisticated flux limiter scheme based on the Chapman-Enskog expansion. However, the flux limited diffusion smoothes artificially energy deposition profiles, especially in situations where ion sources are localized [START_REF] Honrubia | Nuclear fusion by inertial confinement -A comprehensive treatise[END_REF]. This may lead to significant errors in the calculation of ignition thresholds and energy gains. Such diffusion models are employed in all major present-day fluid codes because of their compatibility with the underlying hydrodynamic module.

Several exact methods can be employed to solve the Fokker-Planck equation in a general way, but they are too much time consuming. Monte Carlo algorithms are applied to model charged particle transport in Refs. [START_REF] Lapeyre | Monte Carlo methods fo transport and diffusion equations[END_REF][START_REF] Chaland | [END_REF].

In such an approach, distribution functions are represented by a sum of Dirac measures. Monte Carlo particles are characterized by their numerical weight, their position and their velocity. Those quantities evolve in time according to the Vlasov-Fokker-Planck equation while the tracking of Monte Carlo particles is performed through the spatial mesh. The accuracy of Monte Carlo methods is proportional to N -1/2 , N being the number of Monte Carlo particles, so that N ≫ 1 and variance reduction techniques are usually employed to reduce numerical noise. A significant deficiency of Monte Carlo methods for the investigation of kinetic effects is that the tails of the distribution functions are not described accurately. Moreover, the coupling between suprathermal particles and the thermal bulk is usually treated in a rough manner, by removing the suprathermal particles that are slowed down below a given energy threshold and injecting the removed particles in the thermal bulk. Therefore, the thermalization process is not described with a sufficient precision.

S n methods are also used to solve the Fokker-Planck equation deterministically. They are based on the determination of the angular flux of suprathermal particles at a set of discrete directions, each one associated with a quadrature weight [14,15,16]. Although they are more accurate than diffusion methods and can be extended to highly anisotropic particle distribution functions, the weakly collisional limit is not described accurately and the thermalization process is treated approximately with the same strategy as in Monte Carlo methods. S n methods are usually used to simulate neutron transport and require high computational efforts. For the application of S n methods to suprathermal α-particles transport, we refer to Ref. [START_REF] Honrubia | Nuclear fusion by inertial confinement -A comprehensive treatise[END_REF].

In the present paper we develop a kinetic modeling of suprathermal fusion products in the thermal imploding plasma. We extend the existing code FPion [4,5,6] so as to treat α-particles, for which two scales of energy are considered, namely a suprathermal and a thermal one. Since the developments made to reach this goal have been substantial, they have actually lead to the creation of an entirely new kinetic code called Fuse for FPion Upgrade with two Scales of Energy. This code is able to investigate kinetic effects related to fusion reaction products on the ignition of the hot spot and on the subsequent propagation of the thermonuclear burn wave through the dense fuel. We present here the numerical methods specially designed for the kinetic modeling of α-particles and their validation in several representative tests. Simulations are preformed for a typical ICF DT target, assuming a spherical symmetry in configuration space and axial symmetry in velocity space around the mean velocity. Distribution functions thus depend on one space variable (radius) and two velocity components (radial and azimuthal or perpendicular), depending on the chosen parametrization.

The paper is organized as follows: firstly, we present in Sec. 2 the Vlasov-Fokker-Planck modeling of the fast α-particle transport and collisional relaxation. A specific formalism, based on a two-scale approach with respect to energy is then introduced in Sec. [START_REF] Fraley | [END_REF]. It provides a self-consistent modeling of the coupling between suprathermal and thermal plasma species. Section 4 presents the algorithms devised to solve the two-scale coupling. A finite volume method is applied to the Fokker-Planck equation governing the suprathermal α-particle distribution function. Fast algorithms are then specially designed to solve the discretized model efficiently. Section 5 presents some numerical results regarding the α-particle distribution function evolution and its coupling with the thermal bulk. We show how the methods developed here provide a refined description of the thermalization process. Simulations are carried out in conditions relevant for typical ICF targets. Conclusions are finally presented in Sec. 6.

Physical model for the transport and collisional relaxation of α-particles

Once created by fusion reactions, suprathermal α-particles are transported through an inhomogeneous plasma and slowed down through Coulomb collisions with electrons and thermal ions D and T. Besides, pressure gradients give rise to an electrostatic field E( r, t) that may accelerate or decelerate α-particles. To give an accurate description of the particle transport, as well as the non-local energy and momentum exchange that occur between α-particles and the thermal bulk, a kinetic modeling based on the Vlasov-Fokker-Planck equation is required.

Vlasov-Fokker-Planck equation for the α-particles

The distribution function f α ( r, v, t) of α-particles characterized by a charge Z α e and a mass m α is governed by the Vlasov-Fokker-Planck equation:

∂ f α ∂t + v • ∂ f α ∂ r + Z α e E m α • ∂ f α ∂ v = i ∂ f α ∂t αi + ∂ f α ∂t αe + ∂ f α ∂t fuse . (1)
The first two terms at the right hand side of this equation describe the collisional relaxation of α-particles:

• ∂ f α /∂t| αe stands for the collisions of α-particles with electrons,

• i ∂ f α /∂t| αi describes the collisions of α-particles with thermal ion species. Since thermal species densities are significantly higher than the fast α-particle density (at least at the beginning of the ignition and burn processes), non-linear term corresponding to fast-α/fast-α scattering is neglected. The coupling between the thermalized α particles and the suprathermal ones is naturally included.

We focus now on the collisional part of Eq. ( 1). The Vlasov part of the equation modeling the transport in space and the acceleration due to the electrostatic field is considered separately in Sec. 4. In a fully ionized plasma such as the one considered here, large angle scattering are much less likely than the net large-angle deflection due to a cumulative effect of many small-angle collisions that the projectile experiences along its path [17]. Each of the collision terms in right hand side of Eq. ( 1) can then be expressed as a Fokker-Planck operator in velocity space, which amounts essentially to an advection-diffusion form. More precisely, the slowing down of α-particles on a thermal ion species i can be written as:

∂ f α ∂t αi = 4πΓ αi ∂ ∂ v • m α m i f α ∂S i ∂ v -∇ 2 v T i • ∂ f α ∂ v , (2)
where S i and T i are the so-called Rosenbluth potentials [17] associated to the target ions i. They are defined by a set of Poisson equations in velocity space:

∆ v S i = f i , ∆ v T i = S i . (3)
The coefficient Γ αi = (4πZ 2 α Z 2 i e 4 /m 2 α ) ln Λ αi is proportional to the Coulomb logarithm ln Λ i j (for any species i, j including electrons) related to the Coulomb potential screening and taking quantum effects into account:

Λ i j = λ D / max{λ bar , ρ ⊥ }. The Debye length λ D =         4πn e e 2 /T e + n j=1 4πn j Z 2 j e 2 /T j         -1/2
depends on the temperature T j , which is expressed in energy units. T j is related to the thermal ion distribution function f j by the relation:

T j = m j 3n j (v -V j ) 2 f j ( v) d 3 v,
where n j = f j ( v) d 3 v is the density of ion species j and

V j = n -1 j v f j ( v) d 3 v
is their mean velocity. The characteristic lengths ρ ⊥ and λ bar are the classical and quantum impact parameters:

ρ ⊥ = Z a Z b e 2 /m i j u 2 i j , λ bar = /m i j u i j
where m i j = m i m j /(m i +m j ) is the reduced mass and u i j = √ 3(T i /m i +T j /m j ) 1/2 is an average relative velocity between the particle species i and j. The Coulomb logarithm is thus a particular function of hydrodynamic quantities. It is symmetric with the respect of particle species, Λ i j = Λ ji , which is related to the energy and momentum conservation during the collision.

The effect of electrons on the slowing down of α-particles is modeled by another Fokker-Planck term, in which the electron distribution function is approximated by a Maxwellian characterized by a density n e , a mean velocity u e and a temperature T e :

∂ f α ∂t αe = 1 τ eα ∂ ∂ v • ( v -u e ) f α ( v) + T e m α ∂ f α ∂v α ( v) , (4)
where τ eα is a characteristic e -α collision time defined by:

τ eα = 3 4 √ 2π m α T 3/2 e n e Z 2 α e 4 m 1/2 e ln Λ αe . (5)
Equation ( 4) is obtained by a truncated expansion of the full ion-electron Fokker-Planck operator with respect to the small constant ǫ = (m e /m i ) 1/2 ∼ 0.022 [4,6]. The last term in (1) stands for the creation of α-particles by fusion reactions. The source term is supposed to be isotropic and is given by:

∂ f α ∂t fuse = R DT ( r, t) δ(v -v h ) 4πv 2 , (6)
where v h = 1.3 × 10 9 cm.s -1 is the initial velocity of suprathermal α-particles whose initial energy is 3.52 MeV. R DT is the fusion reaction rate expressed as a function of the distribution functions of D and T, respectively:

R DT ( r, t) = n D n T σv DT = f D ( r, v D , t) f T ( r, v T , t) | v D -v T | σ DT (| v D -v T |) d 3 v D d 3 v T . (7)
The distribution functions f D and f T are solutions of the Vlasov-Fokker-Planck equation written on the deuterium and tritium species, respectively, and they are not necessarily Maxwellian functions. Integrals in Eq. ( 7) are taken over the three-dimensional velocity space.

Dealing with electrons

Since the characteristic time of the considered problem is close to the ion-ion collision time τ ii >> 1/ω pe , ω pe being the electron plasma frequency, and the characteristic length is of the order of the ion collisional mean free path λ i >> λ De , λ De being the electron Debye length, the quasi-neutrality assumption is relevant. We then have:

n e = i Z i n i + Z α n S T α , V e = i Z i n i V i + Z α n S T α V S T α , (8)
where the contribution of suprathermal α-particles is naturally included, n S T α , V S T α being the density and mean velocity of fast α-particles respectively.

Besides, due to a very small ratio of the masses of electrons and ions, the electron equilibration time τ ee is significantly smaller than the mean ion-ion collision time τ ii . According, for example to [START_REF] Braginskii | Transport Processes in a Plasma[END_REF], we have the following ordering of characteristic times: τ ee ∼ ǫτ ii . As a consequence, the electron kinetic equation reduces to a fluid equation. Only an equation for the temperature (or, equivalently, the energy density) is actually needed since the electron density and velocity are known from the quasi-neutrality conditions (8).

In the one-dimensional spherical problem considered here, the electron energy density W e is governed by the following conservation equation :

∂W e ∂t + 1
where κ e is Spitzer's thermal conductivity [19] in the presence of several ion species (see also [7] and Appendix in [20]), the collision time τ e j has been defined in Eq. ( 5) where α is replaced by the considered ion species j. The electron energy density W e and pressure P e are given by an equation of state taking into account Fermi degeneracy [6].

The last term on the right hand side of (9) accounts for the radiation losses of electrons.

Relative importance of electrons and ions on the slowing down of α-particles

3.52 MeV α-particles are created in fusion reactions isotropically, in the system of reference associated with the thermal bulk. Then, they are slowed down through Coulomb collisions with electrons, according to Eq. ( 4), and with thermal ions, according to Eq. ( 2). The relative importance of electrons and ions on the slowing down of α-particles can be estimated by retaining only the dynamical friction terms from the Fokker-Planck equations ( 4) and [START_REF] Atzeni | The physics of inertial fusion[END_REF]. The ratio R i/e between the ion slowing down and the electron one can thus be approximated by:

R i/e = ∂ f α ∂t αi / ∂ f α ∂t αe ∼ T 3/2 e v 3 m 1/2 e m i ∼ T 3/2 e v 3 m 3/2 i ǫ .
The ratio R i/e is thus defined by a characteristic threshold velocity:

v c = ǫ -1/3 (T e /m i ) 1/2 , ( 10
)
so that R i/e ∼ (v c /v) 3 .
The beginning of the slowing-down of α particles is thus governed nearly exclusively by electrons. Then, as v ∼ v c , the effect of ions and electrons on the α relaxation become comparable. Eventually, the final stage of αparticle thermalization is essentially influenced by collisions with thermal ions. Supposing T i ∼ T e , we have the following estimate v c ∼ ǫ -1/3 v th i ∼ 3.6 v th i , v th i being the typical thermal velocity of D and T ions. The effect of thermal ions on the α relaxation dominates when the α velocity is below v c ∼ 3.6 v th i . We shall refer to such α-particles as "moderately suprathermal".

Two-component description of the α distribution function

Physical discussion

From the previous discussion, we know that 3.52 MeV α-particles are firstly slowed down essentially by electrons. The first stage of the α slowing down is thus described by:

∂ f α ∂t coll = 1 τ αe ∂ ∂ v • ( v -u e ) f α ( v) + T e m α ∂ f α ∂ v ( v) . (11)
When v >> u e , the dynamic friction term (first term on the right hand side of ( 11)) dominates so that the α distribution evolves with respect to:

∂ f α ∂t coll ≈ 1 τ αe 1 v 2 ∂ ∂v • v 3 f α (v) . (12)
The stationary solution of (12) behaves as f α ∼ 1/v 3 , where v is suprathermal α-particle velocity. Consequently, as long as fast α-particles remain far from the thermal velocity region, their distribution function varies smoothly over the whole suprathermal velocity region. The associated velocity scale v S T α , defined by:

v S T α ∼ f S T α / ∂ f S T α ∂v , (13)
is in particular greater than the target thermal velocity v th i . Then, when slowed down α-particles get closer to the thermal region but still remain suprathermal, thermal ions tend to dominate the end of the relaxation process, which is then governed by the equation:

∂ f α ∂t coll = i 4πΓ αi ∂ ∂ v • m α m i f α ∂S i ∂ v , (14)
where only the dynamical friction term is retained for the present discussion. We shall deal with the diffusion part separately. Qualitatively, one can consider that the distribution function of the thermal target species i appears highly localized in velocity space, from the suprathermal α-particle point of view. One thus can write:

f i ( v) = n i δ 3 ( v)
(assuming that the mean velocity is zero). Besides, the divergence with respect to velocity that appears on the right hand side of Eq. ( 14) can be expanded as follows:

∂ ∂ v • ∂S i ∂ v f α ≃ ∂S i ∂ v • ∂ f α ∂ v + f α ∆ v S i .
Using the approximation f i ( v) = n i δ 3 ( v), which is valid from the suprathermal α-particle point of view, the first Rosenbluth potential associated to the target ions i can be calculated explicitly: S i (v) ∼ -n i /(4πv). Then, by calculating its derivative, the slowing down of α particles can be modeled by:

∂ f α ∂t coll = i 4πΓ αi m α m i ∂ f α ∂ v • n i 4πv 2 e v + f α f i . (15)
The two terms on the right hand side of Eq.( 15) have a clear physical sense. The first term ∼ ∂ f α /∂ v varies slowly and smoothly far from the thermal velocity region. It can be characterized by a suprathermal velocity scale v S T α , which is greater than the typical thermal ion velocity v th i = (T i /m i ) 1/2 . Actually, the term ∼

n i 4πv 2 ∂ f α ∂ v corresponds to
a conservative convection towards v = 0. The associated convective rate n i 4πv 2 increases as v tends to 0 so that the solution of: (16) tends to a constant f 0 corresponding to the stationary state of (16). The part of the α distribution driven by ( 16) is then stretched and smoothed out as it approaches the thermal velocity region.

∂ f α ∂t coll = i 4πΓ αi m α m i ∂ f α ∂ v • n i 4πv 2 e v
The second term ∼ f α f i appears highly localized in the thermal region of velocity space and behaves qualitatively as a δ-function from the suprathermal α-particle point of view. This term actually leads to the formation of a condensate of width v th i ≪ v S T α . This qualitative analysis shows intuitively how the two-component feature of the α distribution function builds up. It is made of a superposition of two components evolving on two different velocity scales, namely:

• a suprathermal component, fed by fusion reactions and evolving on a large velocity scale, greater than the target thermal velocity.

• A thermal component, corresponding to the thermalized part of the α distribution function, evolving on the same velocity scale as the thermal bulk of the plasma. Note that this component is not fully thermalized since the source term is proportional to i 4πΓ αi f i . There remains a final stage of collisional relaxation between the thermal components of D,T and α ions respectively.

Figure 1 illustrates schematically those processes. From this phenomenological discussion, we can draw a more formal and more rigorous description of the slowing-down which naturally leads to the building of a new multi-scale algorithm solving the initial problem given by Eq. ( 1). The suprathermal component of the α distribution (red) varies on a velocity scale v S T α ≫ v th i . The electron distribution function (green) has a Maxwellian shape with a characteristic width v th e ≫ v th i . The thermal ion component (blue) varies on the thermal ion energy scale ∼ v th i . The contrast between the thermal and suprathermal scales has been reduced artificially for the sake of clarity.

☛ ❝ ❡❛t✁✂✄ ❛t ✸☎✆✝ ▼❡❱ ❢ ✭✈✮ ☛ ✞ ✐ ❞✂✟✁✄❛t❡ ❢ ❙❚ ✠ ❢ ❚ ✠ ❢ ✡ ☞ ✌❤ ✍ ☛ ✞ ✎ ❞✂✟✁✄❛t❡ ✏ ☞ ✌❤ ✑ ✈ ✒✓

Splitting of the Fokker-Planck operator

From the previous analysis, it seems natural to write the α distribution function as follows:

f α ( v, t) = f S T α ( v, t) + f T α ( v, t), ( 17
)
where: f S T α designates the suprathermal component. It is defined on a large velocity domain, spreading to the MeV range. Its typical velocity variation scale v S T α is greater than the thermal ion velocity v th i ; f T α is the thermal component. It is localized in the region of velocity space corresponding to target thermal ion distribution functions and vanishes in the suprathermal velocity domain. The component f T α is designed to describe accurately the final stage of thermalization of the slowed down α-particles. This final relaxation occurrs on a velocity scale ∼ v th i . Let us emphasize that two components defined in Eq. ( 17) do exist in the whole velocity space, the relevant physical quantity being the full α distribution function f α ( v, t).

The idea is then to deal with each component separately. The original Fokker-Planck operator given in Eq. ( 2) is then transformed into a system of two coupled equations governing the two components f S T α and f T α , respectively:

∂ t f S T α αi = Γ αi n i v 2 ∂ v f S T α -n i Γ αi f S T α δ(v) v 2 , ∂ t f T α αi = 4πΓ αi ∂ v • f T α ∂ v S i + 4πΓ αi f i f S T α (v = 0). ( 18
)
The above equations are written in the system of reference associated with the thermal ions. System [START_REF] Braginskii | Transport Processes in a Plasma[END_REF] describes the coupling between the suprathermal component and the thermal one, the coupling function being ∼ f S T α f i , which is subtracted from the equation on the suprathermal component f S T α and appears as a source term in the equation governing the thermal component f T α . The coupling function can actually be approximated for each of the components of the α distribution function in two different ways, depending on the considered velocity scale:

• From the suprathermal component point of view, we have f S T α f i ∼ n i f S T α δ 3 ( v) since thermal target ions appear highly localized. The first Rosenbluth potential S i associated to thermal ions can then be approximated by its temperature-vanishing form.

• From the point of view of the thermal component, we can consider f S T α f i ∼ f S T α (0) f i since the suprathermal component is almost constant on the thermal velocity scale v th i . The term ∼ f S T α (0) f i appears as a source term for the thermal component. It corresponds to a feeding by the suprathermal component.

In Eq. ( 18), we have disregarded the process corresponding to a feeding of the suprathermal component by the thermal one, which could be the case if we modeled large angle collisions, such as α S T + D → α + D S T . Such collisions would build up a suprathermal component for species D and T . This could be naturally included in the formalism that we describe here, but this is a process of second order since the probability of large angle scattering is ∼ 1/ ln Λ times smaller than the pitch-angle collisions modeled by the Fokker-Planck operator.

Diffusion part of the Fokker-Planck operator

We study now the effect of the second term on the right hand side of Eq. ( 2) corresponding to a diffusion in velocity:

∂ f α ∂t αi = - i 4πΓ αi ∂ ∂ v • ∇ 2 v T i • ∂ f α ∂ v . ( 19
)
T i is the second Rosenbluth potential associated to the thermal target ions. The notation ∇ 2 v ( . ) stands for the Hessian ∂ 2 αβ ( . ). Let us define the field J αi , representing the slowing-down current of α-suprathermal particles:

J αi = - i 4πΓ αi ∇ 2 v T i ∂ f α /∂ v, ( 20
)
Using the Dirac-function approximation for the thermal target distribution functions, we can approximate T i by its temperature-vanishing form,

T i (v) ∼ -n i v/(8π).
The approximation is relevant from the suprathermal component point of view. The Hessian ∇ 2 v T i can then be calculated explicitly:

∇ 2 v T i ∼ - n i 8πv Id - v ⊗ v v 2 . (21)
By taking advantage of a polar representation of the velocity v = v e v , where ( e v , e θ ) is the polar local basis of velocity space, the Hessian (21) simplifies to:

∇ 2 v T i ∼ - n i 8πv e θ ⊗ e θ . (22)
The slowing down current defined in Eq. ( 20) expresses the diffusion in velocity associated to the slowing-down process. It is essentially transverse, that is, perpendicular to the local velocity v. Therefore, one can write:

J αi ∼ - Γ αi 2 n i v 2 ∂ f α ∂θ e θ . (23)
The diffusive slowing-down current is thus highly anisotropic in velocity space and it intensifies as α-particles approaches the thermal bulk region of velocity space. Qualitatively, the collisional relaxation of α-particles on thermal target ions is thus characterized by:

• a pure advection in velocity space at a constant rate, modeled by Eq. ( 14), which tends to accumulate α-particles in the thermal ion velocity region.

• An anisotropic diffusion in velocity space, expressed by Eq. ( 23), which tends to make the distribution isotropic when slowed-down α-particles get closer to the final stage of thermalization.

Algorithms for the transport and collisional relaxation of fast fusion products

In this section, we present the numerical methods developed to solve Eq. ( 1) and Eq. ( 18). Those equations govern the time evolution of suprathermal α-particles. Firstly, we show how to deal with the two-component nature of the α distribution function. We then develop a finite volume approach to discretize the equation on the α suprathermal component. An efficient explicit algorithm is then applied to model the time evolution of the suprathermal component with relatively low computational time. We finally present how to simulate accurately the complete thermalization process of α-particles.

Co-existence of two velocity grids

The two-component nature of the α distribution function naturally leads to the co-existence of two velocity grids, namely:

• A suprathermal grid, designed to represent the evolution of the suprathermal component of the α distribution function f S T α . It covers a large domain in velocity, extending to the range v ≃ v h ≃ 1.3 10 9 cm/s, which is the velocity corresponding to the α particles created by fusion reactions. Moreover, since the suprathermal component varies smoothly, we can use a relatively coarse grid to discretize it. f S T α varies significantly on a velocity scale v S T α ≫ v th i , so that the suprathermal grid resolution is typically of the order of one thermal velocity v th i .

• A thermal grid, on which the thermal component of the α distribution f T α is discretized. This grid is designed to capture the final stage of collisional relaxation of the almost-thermalized component of the α distribution on the other thermal ion species D and T. This process entails a velocity resolution much smaller than the local thermal velocity scale v th i . The thermal grid makes use of a cylindrical parametrization (v r , v ⊥ ) inherited from the code Fpion [7]. 000000000 000000000 000000000 000000000 000000000 111111111 111111111 111111111 111111111 111111111 The thick shell of width ∼ T i corresponds to the source term due to fusion reactions. The thermal component evolves on the small and refined cylindrical grid. Both meshes are centered on the mean local bulk velocity V 0 ∼ V e ∼ V i . Velocity space is characterized by an axial symmetry around the axis v r .

v th i ∼ T i /m i v r v ⊥ v = v h
The two grids that are shown in figure 2 are centered on the local mean bulk velocity V 0 (r), which is close to the mean electron velocity V e (r). By using two grids specially-tailored to capture the variations of each component, it is possible to build an efficient algorithm modeling the two components of the α distribution.

Dimensionless form of the Vlasov-Fokker-Planck equation

For numerical purposes, we write the Vlasov-Fokker Planck equation governing the evolution of the suprathermal component of the α distribution function f S T α in a dimensionless form, based on a specified unit system given in Table 1. It is chosen to manipulate numbers that are close to unity. This prevents computational errors caused by under or overflow floating numbers. As it was shown in Eq. ( 23), the collision term between suprathermal α-particles and ions takes a simple form expressed in polar coordinates. The slowing down currents are co-linear with the local polar basis vectors e v , e θ of velocity space. In the spherical one-dimensional geometry considered here, it thus seems natural to parametrize the suprathermal distribution function as f S T α (r, v, θ, t), with two velocity components v = v cos θ e r + v sin θ e ⊥ . Then, the dimensionless equation governing f S T α reads:

∂ f S T α ∂t + v cos θ ∂ f S T α ∂r + E α A α cos θ ∂ f S T α ∂v = i Γ αi ∂ ∂ v • n i v 2 A α A i f S T α e v + 1 2 ∂ f S T α ∂θ e θ + 1 τ eα ∂ ∂ v • ( v -u e ) f S T α + T e A α ∂ ∂ v f S T α - i=D,T,α 4π Γ αi A α A i f S T α f T i + R DT ( r, t) δ(v -v h ) 4πv 2 , ( 24
)
where the normalized constant Γ αi = (4πZ 2 α Z 2 β /A 2 i ) ln Λ αi and the effective electrostatic field E i applied to ions of species i is defined by the following expression:

E i = -(Z i / n e ) ∂ P e /∂r.
(

) 25
Here, n e and P e are the dimensionless electron density and pressure, respectively, and

τ eα = 3 √ πA α T 3/2 e 2ǫ √ 2Z 2 α n e ln Λ αe
is the dimensionless electron-ion collision time.

Table 1: Units defined from reference values of the particle density n 0 and particle thermal energy T 0 .

Quantity Unit density n 0 (arbitrary reference value) thermal energy T 0 (arbitrary reference value) time

τ 0 = T 3/2 0 m 1/2 p /4πe 4 n 0 length λ 0 = (T 0 /m p ) 1/2 τ 0 = T 2 0 /4πe 4 n 0 velocity v 0 = (T 0 /m p ) 1/2 = λ 0 /τ 0 distribution function f 0 = n 0 /v 3 0 first Rosenbluth pot. S 0 = n 0 /v 0 second Rosenbluth pot. T 0 = n 0 v 0 electric field (E i ) E 0 = m p v 2 0 /λ 0 = m p λ 0 /τ 2 0 heat flux Q 0 = n 0 T 3/2 0 /m 1/2 p
Let us consider the third term on the right hand side of (24). From the point of view of suprathermal α-particles, it can be approximated by:

i 4π Γ αi A α A i f S T α f i ≃ 4π i Γ αi A α A i f S T α n i δ 3 ( v), ( 26
)
supposing that v ≫ v th i , V 0 . The term ( 26) is thus highly peaked with respect to velocity in the thermal component region and leads to the formation of a thermalized condensate that cannot be described on the coarse suprathermal grid. That justifies our approach of subtracting this singular term from (24), so that the variations of f S T α remain everywhere smooth and may be described on the suprathermal grid. The term ( 26) is then re-introduced as a feeding term in the equation governing the thermal component, so that the original Fokker-Planck equation governing the complete α distribution function f α = f S T α + f T α is recovered. To solve the full Vlasov-Fokker-Planck equation (24), we use the same general splitting scheme as in the code FPion, namely we treat the advection, the acceleration and the collisional stages separately. We describe now the method developed to solve the collisional part of (24).

Discretization of the collisional term

The collisional part of ( 24) can be written as:

∂ f st α ∂t coll = 1 v 2 ∂ ∂v v 2 J v + 1 v sin θ ∂ ∂θ sin θ J θ , (27)
where the polar components of the slowing down current J are given by:

J v = f S T α v τ eα + Γ αi A α A i n i v 2 + 1 τ eα T e A α ∂ f S T α ∂v , (28)
and

J θ = 1 v ∂ f S T α ∂θ Γ αi n i 2v + 1 τ eα T e A α , (29)
The slowing-down current J takes the general advection-diffusion form in velocity space:

J v J θ = f u v u θ + K vv K vθ K θv K θθ •                  ∂ f ∂v 1 v ∂ f ∂θ                  (30)
where the components of the tensors u and K are related to the Rosenbluth potentials S and T (associated to the target ion species) as follows:

          u v u θ           =                  ∂S ∂v 1 v ∂S ∂θ                  and           K vv K vθ K θv K θθ           =                     ∂ 2 T ∂v 2 ∂ ∂v 1 v ∂T ∂θ ∂ ∂v 1 v ∂T ∂θ 1 v 2 ∂ 2 T ∂θ 2 + 1 v ∂T ∂v                    
which reduces to:

u v u θ = v/ τ eα + i=D,T Γ αi n i /v 2 0 and K vv K vθ K θv K θθ = T e / τ eα A α 0 0 i=D,T Γ αi n i /(2v) . ( 31
)
Note the simplifications implied by using a polar parametrization of velocity space: the dynamical friction coefficient u is indeed co-linear with the radial velocity basis vector e v and the diffusion tensor is diagonal in the basis e v , e θ . We then integrate (27) with respect to velocity on a given cell δV k j of the polar suprathermal velocity grid, subscripts k and j referring to the θ and v directions respectively (see figure 3).

The cell δV k j is defined by its boundaries

θ k-1 2 , θ k+ 1 2 and v j-1 2 , v j+ 1 2 , for 1 ≤ k ≤ k max and 1 ≤ j ≤ j max . We call f n k j = f S T α (v = v j , θ = θ k , t = t n
) the value of the suprathermal distribution function in the cell δV k j at time t n . Integrating Eq. ( 27) over the cell area δV k j , we obtain the following conservative discretized form:

f n+1 k j -f n k j ∆t = 1 v 2 j v 2 j+1/2 J v k j+1/2 -v 2 j-1/2 J v k j-1/2 2δv 3 j + 3v j δv j 2δv 3 j sin θ k+1/2 J θ k+1/2 j -sin θ k-1/2 J θ k-1/2 j δµ k ( 32
)
where discrete elementary volumes are defined by:

δv 3 j = v 3 j+ 1 2 -v 3 j-1 2 , δv j = v j+ 1 2 -v j-1 2 , δµ k = cos θ k+ 1 2 -cos θ k-1 2 .
The centered radial velocity v j that appears in Eq. ( 32) is defined as v j = (v j+ 1 2 + v j-1 2 )/2. In those notations, the discrete volume of the cell δV k j is given by : Besides, a straightforward centered-difference and explicit discretization of the slowing-down current leads to:

δV k j = δV k j 2πv 2 sin θ dv dθ = 4π 3 δv 3 j δµ k . δV jk dv θ k-1/2 θ k+1/2 v ⊥ v r v j dθ v j-1/2 v j+1/2 u 0
J v k j+1/2 = u v k j+1/2 2 ( f n k j+1 + f n k j ) - K vv k j+1/2 δv j+1/2 ( f n k j+1 -f n k j ) (33) J θ k+1/2 j = K θθ k+1/2 j v j δθ k+1/2 ( f n k+1 j -f n k j ), (34)
where the slowing-down coefficient u and the diffusion coefficients K are explicitly given by (31) as functions of velocity. The time varying coefficients in (31) involving thermal ions and electrons are evaluated at the previous time step t = t n .

A Locally Split Explicit scheme 4.4.1. Need for an explicit approach

The slowing-down and diffusion coefficients given in Eq. (31) are thus very inhomogeneous in velocity space, being highly peaked in magnitude near the thermal component region. Besides, the diffusion term is strongly anisotropic (essentially transverse) outside of the thermal component region. In such a situation, the usual implicit schemes may involve the solution of a very large and ill-conditioned linear system that will only give an approximated solution of the non-stationary problem. In this section, we demonstrate how it is possible to take advantage of the strong inhomogeneity of the slowing down current to build an efficient and simple explicit scheme that describes the non-stationary α distribution function time evolution naturally. This approach stems from ideas that were introduced in [21].

The Von Neumann stability condition for the scheme (32) in the case of constant homogeneous slowing-down coefficient u and diffusion tensor K reads as:

(u δt) 2 ≤ 2Tr(K) δt ≤ δv 2 , ( 35
)
where δv is the velocity mesh size. When the slowing-down coefficient u and diffusion tensor K are inhomogeneous (which is the case for our problem), we can apply (35) locally in each cell δV jk of the suprathermal polar velocity grid. Besides, since the scheme (32) is bi-dimensional and parametrized in polar coordinates, (35) actually leads to two stability conditions, corresponding to the radial direction v and the angular direction θ, respectively. Treating these directions separately, the stability condition for (32) can be written for a given cell δV jk as:

• in the radial v direction:

u v j δt δv j 2 ≤ 2(K vv j )δt δv 2 j ≤ 1 (36)
• in the angular θ direction:

2(K θθ j )δt v 2 j δθ 2 k ≤ 1. ( 37
)
Note that the slowing-down coefficient u as well as the diffusion tensor K given in (31) depend only on v.

The idea is then to use the explicit scheme (32) with the stability conditions (36) and (37) applied locally in each cell of the suprathermal grid. Indeed, the discrete scheme (32) corresponds to the finite volume formulation of a conservation equation where the time evolution of the α distribution function defined at the mesh centers is driven by the difference between the numerical fluxes calculated at the boundaries. The fluxes depend on the value of the distribution function in the neighboring cells. If the fluxes are applied during a time step ∆t which is too large with respect to the absolute values of the fields in the neighboring cells, numerical instabilities occur. The idea is then to apply fluxes during a limited time step ∆t ′ , possibly smaller than the imposed time step ∆t. The time interval ∆t ′ is chosen such that the variation of the fields in the neighboring cells remain below their initial absolute values. Fluxes and fields are updated consistently at the frequency 1 ∆t ′ , until the imposed time step ∆t is reached.

Stability and positivity

These conditions impose the stability of the explicit scheme (32), but not necessarily its positivity. Indeed, we have noticed that applying the explicit scheme (32) with the stability conditions (36) and (37) may lead to negative values of f S T α and thus lead to the development of numerical instabilities. This is especially true in the velocity region where the slowing-down coefficient u is large, which may occur for example in the suprathermal region where α-particles are created.

A possible remedy is to introduce an "adaptative de-centering" in the discretization of the radial slowing-down current. We then go back to Eq. ( 33) and introduce the parameters η j such as:

J v k j+1/2 = 1 2 u v k j+1/2 (1 -η j ) f n k j+1 + (1 + η j ) f n k j - K vv k j+1/2 δv j+1/2 ( f n k j+1 -f n k j ). (38)
The choice η j = 0 leads to the centered scheme (33), while η j = 1 leads to a pure upwind scheme. The decentering defined in (38) may also be seen as a perturbation of the discretized diffusion term. Indeed, Eq. 38 can be written in the following form:

J v k j+1/2 = 1 2 u v k j+1/2 ( f n k j+1 + f n k j ) -K vv f n k j+1 -f n k j δv j+1/2 . ( 39
)
The stability condition (36) applied with the modified coefficient diffusion K vv = K vv k j+1/2 + 1 2 u v k j+1/2 η j δv j+1/2 instead of the original K vv defined in (31) leads to the stability condition:

1 2 |u v k j+1/2 | 2 δt ≤ K vv k j+1/2 + 1 2 u v k j+1/2 η j δv j+1/2 and δt δv 2 j+1/2 2K vv k j+1/2 + u v k j+1/2 η j δv j+1/2 ≤ 1.
Besides the positivity condition written in the case of an initial field f S T α localized in one velocity cell leads to:

K vv k j+1/2 + 1 2 u v k j+1/2 η j δv j+1/2 ≥ 0 and 1 δv j+1/2 2K vv k j+1/2 + 1 2 u v k j+1/2 η j δv j+1/2 ≥ |u v k j+1/2 |.
The minimal value of u v η ensuring positivity is thus:

u v k j+1/2 η j = max 0, |u v k j+1/2 | -2K vv k j+1/2 /δv j+1/2 . ( 40
)
To ensure stability as well as positivity, we calculate the radial flux with respect to (39) with η j given by (40) in each velocity cell. Actually, this amounts to using the scheme (32) with the radial diffusion coefficient K vv replaced by:

K vv = max{K vv , |u v |δv/2} (41)
and apply the conditions (36). Note that in (36), the condition imposed on the slowing-down coefficient |u v |δt ≤ δv is automatically fulfilled as soon as the one imposed on the (modified) diffusion coefficient K vv is satisfied.

Applying the stability condition locally

We describe now the accurate implementation of the algorithm, named Locally Sub-cycled Explicit LSE algorithm that solves the problem of collisional relaxation of α-suprathermal particles. The idea is to apply the explicit scheme (32) with the stability conditions (36) and (37) applied locally in each cell of the suprathermal grid.

Knowing the values of the distribution function f n jk in any cell of the suprathermal velocity at time t = t n , we apply the following strategy:

First step -Local time steps calculation For each cell δV jk of the suprathermal velocity grid, we calculate a local time step ∆t jk such that the stability conditions in the θ and v directions (37) are fulfilled. To find ∆t jk , the global time step, namely ∆t, is halved until (37) and (36) are satisfied. The local time step ∆t jk is then:

∆t jk = min(∆t θ jk , ∆t v jk ), (42)
where:

∆t θ jk = 2 -nsplit θ jk ∆t, (43)
and

∆t v jk = 2 -nsplit v jk ∆t, (44)
nsplit θ jk (resp. nsplit v jk ) is the number of times the global time step has to be halved to fulfill the stability condition in the θ (resp. v) direction.

Second step -Sorting the cells Then, the cells of the suprathermal velocity grid are sorted with respect to their local time step ∆t jk calculated above. This can for instance be done with an efficient algorithm (e. g., 'Heapsort' [START_REF] Press | Numerical recipes in C[END_REF]), which takes on the order of N ln N operations for each time step where N is the number of cells of the suprathermal velocity grid. This sorting stage then allows cells to be visited by the algorithm only when they actually need to be updated, and is thus an essential step for an computationally efficient algorithm, as shown in ref. [21].

Third step -Sub-cycling Each cell has to be advanced in both directions v and θ over a time ∆t with respect to its local time-step ∆t jk , this procedure ensuring stability. We thus have to perform a sub-cycling for each cell. The effective computation proceeds through a loop over the smallest local time-step. Inside the loop, the fields (evaluated at the center of the cell) and the flux (evaluated at the borders) are updated consistently with the local time step of the considered cell. More precisely, we perform the following iterations:

f p+1 k j -f p k j ∆t jk = 3v j δv j 2δv 3 j sin θ k+1/2 J θp k+1/2 j -sin θ k-1/2 J θp k-1/2 j δµ k + 1 v 2 j v 2 j+ 1 2 J vp k j+ 1 2 -v 2 j-1/2 J vp k j-1/2 2δv 3 j , (45)
where the superscript p refers to the sub-cycled iterations. The sub-cycling starts with f p=0 k j = f n k j and ends after p max jk iterations where ∆t = p max jk ∆t jk . During the process, the flux J θ k+1/2 j (resp. J v k j+1/2 ) defined in (34) (resp.(39) and (40)), are updated with a frequency corresponding to 1/∆t θ jk (resp. 1/∆t v jk ). For more details on the sub-cycling method, we refer to [21]. This strategy guarantees stability and positivity everywhere on the suprathermal velocity grid.

By applying the local sub-cycling described above, we are able to treat the collisional part of the Vlasov-Fokker-Planck equation governing the suprathermal component of the α distribution function using a tractable explicit approach that does not lead to prohibitive computational time.

To illustrate the efficiency of the LSE algorithm, we present in figure 4 the map of nsplit θ jk and nsplit v jk defined in (43) and (44) on the suprathermal velocity grid. We consider two locations corresponding to the hot spot and the dense shell of a typical imploding capsule taken 1 ns before stagnation. We note that the sub-cycling is more expensive in the dense shell region than in the hot spot. Indeed, the high density and low temperature of the shell imply smaller time step.

Furthermore, considering the maps of nsplit θ jk represented at the bottom of figure 4, we note that to advance the fields in θ, we mainly have to sub-cycle the most central cells, where the local time step imposed by the stability condition is the smallest since the local cell size v j δθ is small close to the center. For the outermost velocity cells, no sub-cycling is actually needed.

v ⊥ v r v 3/2 θ k-1/2 θ k+1/2 +f (1, k) -i 4πΓ αi f ST α f i

Coupling with the thermal component

We now discuss the implementation of the coupling strategy between the suprathermal and the thermal components, as described by system (18) in Sec. 3.2.

From the suprathermal point of view

From the point of view of suprathermal α-particles, the coupling with the thermal component is made by the third term in Eq. ( 26) on the right-hand side of (24). It induces a time variation of the suprathermal distribution given by the following equation:

∂ f S T α ∂t S T →T = - i 4π Γ αi A α A i f S T α f i ≃ - i 4π Γ αi A α A i f S T α n i δ 3 ( v). ( 46
)
The time evolution of the suprathemal distribution function in central velocity meshes is then governed by:

∂ f st α ∂t coll = 1 v 2 ∂ ∂v v 2 J v + 1 v sin θ ∂ ∂θ sin θ J θ - i Γ αi A α A i f S T α n i δ(v) v 2 , (47)
where the slowing-down currents J v and J θ are given by Eq.(28) and Eq.( 29) respectively. As slowed down α-particles approach the thermal velocity region, the transverse diffusion current J θ intensifies so that the distribution function is almost isotropic in the central velocity meshes. Eq.(47) simplifies to:

∂ f st α ∂t coll = 1 v 2 ∂ ∂v v 2 J v - i Γ αi A α A i f S T α n i δ(v) v 2 , ( 48
)
where the slowing-down current J v can be approximated by:

J v ≃ Γ αi A α A i n i v 2 f S T α .
We then integrate Eq.(48) over a central mesh ( j = 1, 1 ≤ k ≤ k max ) of the suprathermal velocity. The suprathermal component in the central meshes corresponding to j = 1 are then calculated as follows (see Fig. 5):

f n+1 k1 -f n k1 ∆t v 3 3/2 3 = i n i Γ αi ( f n k3/2 -f n k1 ). (49)
In such a way, the distribution function remains stable in the most central part of the suprathermal velocity grid.

From the thermal point of view

To recover the full Fokker-Planck equation on the physical α distribution function f α = f T α + f S T α , we define an α thermal component f T α , which evolves on the thermal velocity grid defined above. This is also the grid on which the thermal ion D, T distribution functions evolve. This grid is actually inherited from the code FPion, so that we use the same cylindrical parametrization as explained in [7] for the α thermal component: f T α (r, v r , v ⊥ ), v r and v ⊥ being the radial and tangential components of the velocity, respectively.

The term [START_REF] Peigney | [END_REF] subtracted from the suprathermal component equation reappears as a source term in the Vlasov-Fokker-Planck equation governing the thermal component of the α distribution function f T α , so that the relaxed suprathermal component feeds the thermal one and no α particle is lost in the process:

∂ f T α ∂t + v r ∂ f T α ∂r + v ⊥ r v ⊥ ∂ f T α ∂v r -v r ∂ f T α ∂v ⊥ + E α A α ∂ f T α ∂v r = i 4π Γ αi ∂ ∂ v • A α A i f T α ∂S i ∂ v -∇ 2 T i ∂ f T α ∂ v + 1 τ eα ∂ ∂ v • ( v -u e ) f T α + T e A α ∂ ∂ v f T α + i 4π Γ αi A α A i f S T α f i . ( 50
)
The source term coming from the slowing down of the suprathermal component appears in the last term on the right-hand side of (50). From the point of view of the thermal component, the suprathermal component f S T α appears relatively constant over the whole thermal velocity grid since it varies significantly on the coarse suprathermal velocity grid whose mesh size is of the order of the thermal velocity. That is why we use the following estimate:

i 4π Γ αi A α A i f S T α f i ∼ f S T α (V 0 ) i 4π Γ αi A α A i f i , (51)
V 0 being the mean ion velocity. This procedures guarantees an exact mass conservation: the number of particles that are removed form the suprathermal component are injected into the thermal component. Note that the source term feeding the α thermal component depends on the thermal distribution functions of all thermal ion species. To solve (50), we use algorithms inherited from the code FPion. Their numerical implementation are for example discussed in [7].

Transport and acceleration of the suprathermal component

We discuss in this section the algorithm developed to solve the Vlasov part of Eq. ( 24)), namely:

∂ f S T α ∂t + v • ∇ r f S T α + E α A α • ∂ ∂ v f S T α = 0 (52)
We deal with the advection and acceleration separately.

Advection

In this stage, we solve the pure advection equation on the suprathermal component f S T α for a given velocity v:

∂ f S T α ∂t + v • ∇ r f S T α = 0, ( 53
)
whose exact solution is given by:

f S T α ( r, v, t + ∆t) = f S T α ( r -v∆t, v, t). ( 54
)
Thus, solving (53) amounts to interpolating (54) on the whole phase space. We thus start with a given point (r, v, θ) of the phase space, v, θ being chosen on the polar suprathermal velocity grid. We have to compute the transformation of the suprathermal phase space coordinates r, v, θ during one time step ∆t. Since the suprathermal velocity grid is centered on the mean bulk velocity V 0 , we firstly project the polar velocity coordinates on the cylindrical basis:

v r = V 0 + v cos θ, v ⊥ = v sin θ. ( 55
)
Then, we apply the following transformations on r, v r , v ⊥ over one time step ∆t:

r(t -∆t) = r(t) 2 -2r(t)v r (t)∆t + v 2 ∆t 2 1/2 , v r (t -∆t) = r(t)v r (t) -v 2 ∆t r(t -∆t) , v ⊥ (t -∆t) = r(t)v ⊥ (t) r(t -∆t) , ( 56
)
which gives us the advected point in phase space. For the interpolation in space, we have to find the two consecutive nodes r i 0 and r i 0 +1 of the spatial mesh such that r i 0 ≤ r(t -∆t) ≤ r i 0 +1 . Then, for each spatial nod r i 0 (respectively r i 0 +1 ), we have to carry out an interpolation of (56) on the polar suprathermal velocity grid centered on the local mean bulk velocity V 0 (r i 0 ) (respectively V 0 (r i 0 +1 )). We thus calculate:

v(t -∆t) = (v r (t -∆t) -V 0 (r i )) 2 + v 2 ⊥ (t -∆t) 1/2 , θ(t -∆t) = cos -1 v r (t -∆t) v(t -∆t) , ( 57
)
for i = i 0 and i = i 0 + 1. We then interpolate (57) on the nodes of the suprathermal velocity grid centered on V 0 (r i ), using a simple linear interpolation method. This gives us the advected points:

f i 0 = f S T α (r i 0 , v(t -∆t), θ(t -∆t), t -∆t), f i 0 +1 = f S T α (r i 0 +1 , v(t -∆t), θ(t -∆t), t -∆t). ( 58
)
The final stage is a cubic interpolation with respect to space:

f S T α (r(t -∆t), v(t -∆t), θ(t -∆t), t -∆t) = f i 0 + pδr f ′ i 0 + p 2 [3δ f -δr(2 f ′ i 0 + f ′ i 0 +1 )] + p 3 [δx( f ′ i 0 + f ′ i 0 +1 ) -2δ f ] with δr = r i 0 +1 -r i 0 , p = r(t -∆t)/δr, δ f = f i 0 +1 -f i 0 .
In this equation, the spatial gradients f ′ i 0 and f ′ i 0 +1 are evaluated by finite differences. The slopes are limited to prevent unphysical over/undershoots in the interpolation process.

Acceleration

The electric field effect on the α suprathermal component is modeled by:

∂ f S T α ∂t + E α A α ∂ f S T α ∂ v = 0 ( 59
)
where the effective electrostatic field E α is defined by Eq. ( 25). Here again, we use a method of characteristics to solve (59) since an acceleration can be seen as an advection in velocity. The situation gets simpler here, since we only have to carry out an interpolation in velocity on the suprathermal velocity grid. The process is repeated independently in each spatial cell.

Chain of algorithms to solve the suprathermal Vlasov-Fokker-Planck problem

We conclude this section by summarizing the sequence of algorithms that have been developed to solve the whole problem of creation, transport and collisional relaxation of α suprathermal particles, consistently with a ion-kinetic treatment of the plasma thermal bulk. In particular, we show how the algorithms related to the suprathermal components are linked with those dealing with electrons and thermal ion distribution functions. This constitutes the main loop of our kinetic code Fuse. For a global time step ∆t, we apply the following splitting sequence:

Step 1 -Electron conductivity We solve the conduction part of ( 9), which takes the form of a pure diffusion (or heat) equation during the time ∆t/2.

Step 2 -Acceleration We accelerate ion thermal distribution functions for species D, T, α over the time ∆t/2, and at the same time we solve the convective part of (9), which enables us to improve the energy conservation between ions and electrons (see [6]). Then, we accelerate the suprathermal α component.

Step 3 -Advection We carry out the advection of thermal components for every ion species D, T, α as well as the suprathermal α component over the time ∆t/2.

Step 4 -Feeding the suprathermal component The suprathermal α component is fed by the fusion reaction according to (6) applied over the whole time step ∆t.

Step 5 -Suprathermal collisional relaxation We next solve the collisional part of ( 24) applying the Locally Split Explicit (LSE) algorithm over the time step ∆t.

Step 6 -Feeding the thermal component We apply the feeding term (51) of the α thermal component by the suprathermal one over the time step ∆t.

Step 7 -Thermal collisional relaxation We perform the collisional relaxation of every ion thermal distribution functions (for ion species D, T, α) on thermal ions and on electrons, applying the same algorithms as in Fpion. Note that the collisional relaxation of ion distribution functions on themselves is non-linear and is solved using Crank-Nicholson iterations with an ADI scheme (see Appendix of [20]).

Step 8 -Advection Step 3 is repeated for another ∆t/2.

Step 9 -Acceleration Step 2 is repeated for another ∆t/2.

Step 10 -Electron conduction Step 1 is repeated for another ∆t/2.

After each modification of the ion distribution functions (thermal or suprathermal), the ion moments as well as the slowing-down and diffusion coefficients are updated consistently.

Validation of the code by test problems

In this section, we apply the algorithms developed to model the collisional relaxation and thermalization of αparticles in simplified configurations where analytical results are known.

Isotropic time-dependent test problem

In this first test problem, we consider the collisional relaxation of fast α-particles in an homogeneous and steady plasma made of one mean ion species Z i = 1, A i = 2.5 and electrons. The reference density is n i = n e = 10 22 particles/cm 3 , and the temperature is 1 keV. We keep those conditions constant during the test problem calculation. Suprathermal α particles are then injected isotropically at the energy 3.52 MeV at a steady rate S 0 (particles.cm -3 .s -1 ), so that the suprathermal component remains isotropic during the slowing down process. Following our two-scale approach, the α distribution function

f α (v, t) = f S T α (v, t) + f T α (v r , v ⊥ , t)
is the solution of:

∂ t f S T α = Γ αi n i v 2 ∂ v f S T α + 1 τ αe v 2 ∂ v v 3 f S T α -4πn i Γ αi f S T α δ(v) 4πv 2 + S 0 δ(v -v h ) 4πv 2 , ∂ t f T α = ∂ t f T α αi + ∂ t f T α αe + 4πΓ αi f i f S T α (0). (60)
∂ t f T α αi (resp. ∂ t f T α αe ) corresponds to the collisional terms of the thermal ions (resp. electrons) on the α-thermal particles.

In those conditions, we have the characteristic velocity scales, expressed in cm/s:

v th i ∼ 3.0 × 10 7 << v c ∼ 1.1 × 10 8 << v h ∼ 1.3 × 10 9 < v th e ∼ 4.2 × 10 9 (61)
For v > v c (v c given in Eq. ( 10)), the slowing down of α-particles is mainly due to the Coulomb collisions with electrons. The suprathermal component f S T α (v, t) then tends to the stationary solution of:

∂ t f S T α = 1 τ αe v 2 ∂ v v 3 f S T α + S 0 δ(v -v 0 ) 4πv 2 . ( 62
)
The stationary solution is given by:

f 1 (v) = S 0 τ αe v 3 H(v 0 -v), v > v c , (63)
where v h is the velocity corresponding to the injected α-particles at 3.52 MeV, which corresponds to v h ∼ 1.3 × 10 9 cm/s, and H is the Heaviside distribution. We plot f S T α (v, t) calculated by Fuse at different times as well the stationary analytical solution given by (63) (see Fig (6)). The numerical solution agrees with (63) as long as

v > v c . When v < v c , 0 3e+08 6e+08 9e+08 1.2e+09 1.5e+09
Velocity (cm/s) ions tend to dominate the slowing down of the α-particles and the suprathermal component solution of (60) tends to a stationary state that is almost constant close to thermal ions. This is due to the removal of the term ∝ f S T α n i δ 3 ( v) in the collision term governing the slowing down of f S T α . The suprathermal component actually feeds the thermal one, the feeding process being driven by the source term ∝ f S T α (v = 0) f i . The thermal component subsequently evolves towards a Maxwellian characterized by the total density n α of α particles injected in the system, and the reference temperature T 0 (which is kept constant during the test problem calculation):

M α (v) = n α m α 2πT 0 3/2 exp - m α v 2 2T 0 . (64)
The total density is given by:

n α = τ s 0 S 0 dt, (65)
τ s being the time when the source is shut down. The convergence to the Gaussian (64) is represented on Fig (6). Note that this convergence is calculated on the refined thermal grid. The α thermal component is fed by a source term ∝ f i , of width ∼ √ T 0 /m i , and relaxes on the thermal grid towards the Gaussian (64) of width ∼ √ T 0 /m α .

Anisotropic time-dependent test problem

We next consider the following anisotropic test problem. We consider an initial condition for the α suprathermal component highly localized in velocity space. Namely, we take:

f S T α (v, θ, t = 0) = n α δ(v -v 0 ) 4πv 2 δ(cos θ -cos θ 0 ), (66)
with v h = 1.3 × 10 9 cm/s and θ 0 = π/4. We then let the suprathermal α distribution slow dow on electrons and on thermal ions. As previously, the thermal plasma is homogeneous and made of one ion species Z i = 1, A i = 2.5 and electrons. The temperature of the thermal plasma is kept constant during the calculation: we take T 0 = 5 keV. In those conditions, the characteristic velocity scales are (in cm/s):

v th,i ∼ 6.9 × 10 7 << v c ∼ 2.4 × 10 8 << v h ∼ 1.3 × 10 9 < v th,e ∼ 9.4 × 10 9 (67)
The evolution of the α distribution function is represented in Fig. 7. As long as v > v c , the momentum and energy losses by the fast ions to the background plasma electrons are the dominant process. The distribution function remains highly localized in velocity space around a velocity v b (t) that declines due to the slowing down on electrons. The velocity of the bulk v b (t) can be calculated analytically [START_REF] Rax | Physique des Plasmas[END_REF]:

v b (t) = [(v 3 0 + v 3 c ) exp - 3t τ αe -v 3 c ] 1/3 (68)
The comparison between the code and the exact solution is represented on Fig. 8 and reveals a pretty good agreement, as long as v > v c . Then, as v ≤ v c , the energy diffusion process as well as the perpendicular diffusion due to the thermal ions become significant. The α distribution function is scattered in the θ direction, due to the diffusion on the thermal ions, that intensifies as v → 0. Consequently, as v → 0, the α suprathermal distribution tends to become isotropic while feeding the thermal component. Finally, the thermal component then converges towards the Gaussian, as in the first test problem. To model properly what happens in the vicinity of the thermalization, for v ∼ v th i , we solve the full Coulomb operator applied to the α thermal component f T α that evolves on the thermal refined grid. This guarantees a proper modeling of the thermalization of the α distribution function, as it slows down, scatters and diffuses in energy in joining up with the background thermal ions.

Energy conservation

We finally consider a full collision relaxation process, starting from an isotropic α suprathermal component that slows down through collisions on the electrons and the thermal ions. In this test problem, the electron (res. ion) temperatures evolve consistently with the slowing down of the suprathermal particles. More precisely, as v > v c , suprathermal particles slow down essentially on electrons. The electron temperature thus increases. Then, due to the collisional relaxation of thermal ions with electrons, the thermal ion temperature increases. When the suprathermal particles reach the thermal velocity region, the α thermal component builds up and a collisional relaxation between electrons and thermal ions (including the α thermal component) brings the system to a stationary state. The aim of this test problem is to illustrate that the way we solve the coupling between the suprathermal component and the thermal background ensures the conservation of mass and energy. We check that the total mass remains constant (with a numerical error less than 1% due to the finite size of the velocity mesh). We plot the time evolution of the temperatures (electrons, thermal background ions and α-thermal component) on Fig. (9). We show how the system evolves naturally to a stationary state calculated by the algorithm described above. The total energy variation of the system between the initial state and the final stationary state is less than 1%. Our original algorithm based on a 2-scale approach to model the collisional relaxation between suprathermal particles and the thermal background is thus validated in simplified test problems where exact results are known. Besides, the mass and energy conservation principles are fulfilled at a discrete level. We can consider that our code Fuse is reliable. We then apply it on real target configurations.

Application on the ignition and thermonuclear burn of typical ICF capsules

We apply the numerical scheme presented in Sec. 4 to model a typical spherical implosion of a cryogenic DT capsule. Our code allows us to study ion-kinetic effects during the ignition stage and the beginning of the thermonuclear burn stage.

Initial conditions

We consider the same fluid reference simulation as in [6] corresponding to an ICF target with parameters typical of ignition capsules designed for the LMJ and NIF laser [START_REF] Saillard | [END_REF] and [24]. Namely, we consider a 0.3 mg cryogenic DT layer deposited on the inner surface of a CH shell of a 1 mm (inner) radius. The kinetic calculation is started at t = 17 ns after the beginning of the implosion, when the main converging shock reaches the center of the target. The boundary condition is taken from the hydrodynamic simulation. The densities, temperatures and velocities are recorded on the fuel/pusher interface in the fluid simulation.

The kinetic simulation considers three ion species, namely D, T and α. Initially, only thermal species D and T are present. They give birth to suprathermal α particles in the fusion reactions. The relaxation of the suprathermal α component then leads to the creation of an α thermal component interacting with the other thermal ion distribution functions (D and T, respectively). Note that the thermal bulk is described in more details than in [6] where a single mean ion species with a mass number of 2.5 was considered.

In our kinetic simulation, the position of each spatial meshes is updated after each time step with respect to the imposed boundary condition and to the fixed number of spatial meshes i max . This updating is performed before each advection phase. This means that the position of a given spatial cell r i 0 , with 1 ≤ i 0 ≤ i max is time dependent, decreasing with the size of the imploding system. To represent in a satisfactory manner both the dense region where the fluid simulation grid is the finest and the central zone where it is rather coarse, we employ 78 cells with a geometrically varying mesh size (with the ratio 0.97) so that the mesh size δr is decreasing from 20 µm near the center to less than one micron near the outer boundary. The thermal velocity space (v r , v ⊥ ) is discretized into 129 × 64 cells, whereas the suprathermal velocity grid (v, θ) makes use of 100 × 60 cells. The reference time-step value is 0.05 ps.

Comparison with Fpion and FCI1

To validate the thermal part of our code Fuse, we compare the density, velocity and temperature profiles with the hydrodynamic code FCI1 as well as with the kinetic code FPion at two different times of the implosion:

• at t = 17.1 ns, that is to say 100 ps after the beginning of the implosion. We find a pretty good agreement between the Fuse kinetic calculation and the FCI1 fluid simulation (Fig. 10). The kinetic modeling reveals a significant anisotropy on the ion temperatures (and pressures), as the one observed with FPion [6]. The anisotropy then tends to disappear during the implosion.

• At t = 17.65 ns, in the vicinity of the target stagnation, Fuse and FCI1 are still in good agreement. However, we note that the compression zone near the inner interface of the dense fuel lies closer to the target center in the kinetic calculation (see the negative velocity gradient region about r = 70 µm on the right part of Fig. 11). This result has already been obtained with FPion and discussed in [6]. This is related to a higher ion heat flux, which tends to increase the rate of ablation of the cold fuel by the hot spot.

As long as t ≤ 17.65 ns, the α-particles number is small, so that the above comparisons between the codes Fuse and Fpion (which does not take α-particles into account) are relevant and tend to validate the methods programmed in Fuse regarding the thermal background (thermal ions and electrons).

Transport of α particles

We analyze the transport of suprathermal α particles throughout the capsule. Figure 12 shows the spatial density profiles during the implosion for the suprathermal and thermal components of α-particles. At early times, suprathermal α-particles are produced in the hot central region of the capsule and deposit their energy in the surrounding cold shell. The region corresponding to the suprathermal α energy deposition is indicated by a sharp decreasing of the suprathermal density profile. This occurs at a distance which corresponds to the collisional mean free path of suprathermal α particles. Meanwhile, the slowing down of suprathermal α particles feeds the thermal component, that process corresponding to the bump observed in the thermal α density profiles (Figure 12-right).

During the implosion process, the α collisional mean free path decreases, so that the α suprathermal particles are trapped in a smaller radius. In the mean time, the production of suprathermal α-particles intensifies due to the increasing ion temperature. As a result, the suprathermal α density increases. f S T α (r i 0 (t), v, θ, t) is presented in figure 13.

The suprathermal distribution function is rather anisotropic. It is highly peaked toward positive velocities v r > 0. This can be explained by the inhomogeneous fusion reaction source term, which strongly depends on the ion local temperature. Since T i is more peaked towards the center of the capsule, as it can be seen in the temperature profiles in figure 10, an observer located outside of the highly emissive central region sees the suprathermal α-particles passing from the center to the outside. That leads to a local distribution shape shown in the top panel of figure 13. The spatial gradient of the fusion reaction source term (6) thus accounts for the anisotropy of the suprathermal α distribution function.

Let us consider the cell i 0 with the radius such that r i 0 (t) = λ α (ρ(t)), where λ α is the collisional mean free path of a suprathermal α particle and ρ the mean density of the capsule. As α-particles deposit their energy in the considered spatial cell i 0 , which corresponds to the sequence shown in figure 13, the suprathermal α distribution function slows down significantly towards the thermal velocity region. During this slowing down process, the distribution function tends to spread over a wider domain in the polar angle θ. This is a consequence of the diffusion part of the Fokker-Planck equation, which leads to a mainly transverse slowing-down current that intensifies close to the thermal velocity region.

To check that the collisional module of the code behaves correctly in a real target configuration, we artificially do not calculate the effect of the advection and acceleration on the α-suprathermal component, so that the time evolution is driven by the collisions on electrons and thermal ions only. The corresponding time evolution is represented in Fig. 14. This numerical test is closed the third test problem presented in Sec.4.9.3, but is carried out in thermodynamic conditions corresponding to real ICF target configuration. The suprathermal particles are initially distributed anisotropically in velocity space with respect to Fig. 14 (top-left). For v ≥ v c ∼ 3 -4v th i , fast ions mostly slow down by collisional drag on the background electrons with very little pitch-angle scattering. The fast ions stay mostly in their original pitch-angle direction. For v ≤ v c , the suprathermal particles slow-down predominantly on the thermal background ions and scatter in pitch-angle. The suprathermal distribution function tends to be isotropic as it approaches the thermal velocity region. The suprathermal grid resolution is fine enough to represent the variations of the suprathermal component, that tends to be constant as it gets closer to the thermal velocity region.

Feeding the thermal component

When the slowed down suprathermal α-particles reach the thermal velocity region, a fraction of α-particles is removed from the suprathermal component, to feed the thermal component according to Eq. (50). The sequences represented in figures 13-14 illustrates this coupling from the suprathermal component point of view. The distribution function remains stable, while the particles are accumulating in the vicinity of the thermal region. Without the removal of the term (26) on the right hand side of Eq. ( 24), the suprathermal distribution function would have become unstable as v → V 0 . The evolution of the thermal component of the α-particle distribution function represented in figure 15. It shows how the thermal component builds up. 1. Times refer to beginning of the kinetic calculation.

Ignition and burning wave propagation

We finally give the density, velocity and temperature profiles calculated by Fuse and compare the results with the fluid code at the time t = 17.85 ns (Fig. 16) . After that time, corresponding to the arrival of the flame near the outermost cells, the kinetic simulation may not be relevant since the boundary condition (which comes from the hydrodynamic calculation) may not be consistent with the pressure calculated by the kinetic code. In the kinetic calculation, the heating of the hot spot appears to be faster than in the fluid code. This is consistent with the differences observed during the implosion phase, where the dense zone corresponding to the ablated cold fuel was imploding faster in the kinetic calculation. Besides, the kinetic ion temperature profile displays a preheating wave ahead of the main temperature front. This is specially visible on the ion temperature profiles of (Fig. 16). This structure is related to the Bragg peak of the D,T ions located in the dense cold fuel cold. Suprathermal α-particles are created mainly in the central hot spot and deposit their energy and momentum near the inner interface of the cold fuel, where the thermal ion heating occurs. This interpretation will be examined more closely with future kinetic calculations of different target designs (that may be less efficient than the one considered here). By applying the efficient algorithm (based on a 2-scale approach) exposed and validated in Sec.4 on real target configurations (that could not be solved analytically), the code Fuse is able to simulate the fuel of real ICF targets at a kinetic level over a time corresponding to 1 ns after the start of the implosion. One thus models the ignition and the beginning of the burning wave propagation. Besides, by making use of a parallelization method of the collisional part of the code (which is possible since we can calculate the effect of collisions in each spatial cell independently from the others), it takes less than 1 day of computation time, which is roughly twice as long as the usual simulations performed by Fpion (corresponding to the implosion phase without α-particles).

Summary and perspectives

We have developed a numerical strategy to model fast α-particles produced by fusion reactions at a ion kinetic level. A two-scale approach has been specially-tailored to represent the two-component nature of the α distribution function and simulate the thermalization process accurately.

Efficient algorithms have been designed to simulate the time evolution of the fast α component, driven by the transport in the inhomogeneous thermal plasma as well as the Coulomb collisional relaxation on electrons and ions. The energy and momentum exchange between fast fusion products and the thermal plasma are thus calculated at the kinetic level. The methods have been tested in thermodynamic conditions corresponding to typical DT targets close to ignition. It has been shown that a locally split explicit scheme can be used to describe the fast α population evolution in non-prohibitive computational time. Besides, the algorithms presented here are easily parallelizable to take advantage of present-day multi-core architectures.

The ion-kinetic code Fuse, built as an extension of the former code FPion, is thus able to model a full DT target implosion, including the ignition and burn processes, at a ion-kinetic level. Investigating in more details the role of kinetic effects of fusion products in the ignition and burn of DT targets is the purpose of ongoing work and will be published elsewhere [START_REF] Peigney | [END_REF]. We may have in view to study implosions in the vicinity of the ignition threshold, where kinetic effects should be enhanced and may modify the energy gain.

Finally, the algorithms developed here may be naturally extended to add the effect of Boltzmann-type large angle scattering, that would feed a suprathermal component for the D,T ions. Neutron momentum and energy deposition may be modeled in a similar way.

Figure 1 :

 1 Figure 1: Schematic representation of the collisional relaxation of suprathermal α-particles on thermal target ions i. The suprathermal component of the α distribution (red) varies on a velocity scale v S Tα ≫ v th i . The electron distribution function (green) has a Maxwellian shape with a characteristic width v th e ≫ v th i . The thermal ion component (blue) varies on the thermal ion energy scale ∼ v th i . The contrast between the thermal and suprathermal scales has been reduced artificially for the sake of clarity.

Figure 2 :

 2 Figure2: Schematic representation of the two velocity grids used to model the α suprathermal and thermal components respectively. The suprathermal component evolves on the coarse polar grid, covering a wide domain extending to the MeV region. The thick shell of width ∼ T i corresponds to the source term due to fusion reactions. The thermal component evolves on the small and refined cylindrical grid. Both meshes are centered on the mean local bulk velocity V 0 ∼ V e ∼ V i . Velocity space is characterized by an axial symmetry around the axis v r .

Figure 3 :

 3 Figure 3: The suprathermal velocity grid.

Figure 4 :

 4 Figure 4: Map of nsplit v (top) and nsplit θ (bottom) represented in the suprathermal velocity grid in 2 locations. On the left, we consider a point in the hot spot where: n e ∼ 10 21 cm -3 and T i ∼ T e ∼ 0.5 keV. On the right, we focus on a point taken in the dense shell where: n e ∼ 10 24 cm -3 and T i ∼ T e ∼ 0.01 keV. Those conditions correspond to a typical implosion 1 ns before stagnation. Illustrations are given for a global time step ∆t = 0.1 ps

Figure 5 :

 5 Figure 5: Central mesh of the suprathermal velocity grid

Figure 6 :

 6 Figure 6: Time evolution of the α distribution function corresponding to the isotropic test problem. Distribution functions are expressed in cgs units, namely in cm -4 .s -1 . The α source term is shut down after τ s = 1 ns.

Figure 7 :

 7 Figure 7: α suprathermal distribution solution of the anisotropic test problem at different times. Final stages of collisional relaxation. The values of the distribution function are expressed in cgs units.

Figure 8 :

 8 Figure 8: Time evolution of the velocity corresponding to the maximum of the α suprathermal distribution function solution corresponding to the anisotropic test problem.

Figure 9 :

 9 Figure 9: Time evolution of electron and thermal ion temperatures corresponding to a pure collisional relaxation test problem.

Figure 10 :

 10 Figure10: Profiles of the density, velocity and of the electron and total ion temperatures in a DT ignition target at the time t = 17.1 ns, which corresponds to 100 ps after the beginning of the kinetic calculation and roughly 1 ns before the target stagnation.

Figure 11 : 5 . 4 .Figure 12 :

 115412 Figure11: Profiles of the density, velocity and of the electron and total ion temperatures in a DT ignition target at the time t = 17.65 ns, which corresponds to 650 ps after the beginning of the kinetic calculation. This time is also just before the target stagnation.

Figure 13 :

 13 Figure 13: α suprathermal distribution observed in a given mesh of the imploding hot spot at different times. The simulation takes into account the creation, the transport and the collisional relaxation of α particles. The values of the distribution function are expressed in cgs units. Times refer to beginning of the kinetic calculation.

Figure 14 :

 14 Figure 14: α suprathermal distribution observed in a given mesh of the imploding hot spot at different times, when only the collisional relaxation is considered, starting from a given anisotropic initial state. Times refer to beginning of the kinetic calculation.

Figure 15 :

 15 Figure 15: Thermal component of the α distribution function observed in a given cell of the imploding hot spot at different times. This component is fed by the relaxation of the suprathermal component. The values of the distribution function are expressed in the units presented in table 1. Times refer to beginning of the kinetic calculation.

Figure 16 :

 16 Figure 16: Profiles of the density (n), velocity (u) (top), of the electron (T e ) and total (T i ) ion temperatures (bottom) in the fuel at t = 17.85 ns in a DT ignition target implosion; left : fluid simulation results, right : kinetic simulation results.

  Time evolution of the α suprathermal component for t ≤ τ s . The chosen time interval between 2 consecutive curves is 0.05 ns. The exact solution is represented in dashed lines. Time evolution of the α thermal component for t ≤ 4 ns. The suprathermal source term is shut down after 1 ns, the α thermal component relaxes towards the Gaussian.

	0.0001			0.02	Analytical stationary solution
	1e-07 1e-06 1e-05 Suprathermal α component	Numerical transient solution Numerical stationary state Analytical solution	Thermal α component (cgs)	0.005 0.01 0.015	Numerical solution
	1e-08
	0 (a) -1e+08 (b)	-5e+07	0 Velocity (cm/s)	5e+07	1e+08

Acknowledgments. The authors are grateful to Professors Xavier Blanc, Josselin Garnier, Rémi Sentis and Gerald Samba for fruitful discussions on the subject.


Jean-Luc Menet
email: jean-luc.menet@univ-valenciennes.fr

Ion-Cosmin Gruescu

COMPARATIVE LIFE CYCLE ASSESSMENT OF TWO INSULATED WALLS

Keywords: Life Cycle Assessment, LCA, wall, environmental impacts, climate change indicator 1

The present study presents some results obtained by applying the LCA methodology to evaluate the environmental footprint of a wall. The approach is motivated by the increasing number on markets of new materials such the mono-wall blocks, which are offering interesting insulation properties.

The proposed study focuses on a comparison between two walls with the same thermal resistance. The environmental footprint of a wall made by using mono-wall blocks is then compared to a classical wall made with breeze blocks and glass wool.

The environmental impacts are evaluated by using the Eco-invent 2.0 data-base, with midpoint indicators associated to the CML method.

The main result of the study is very interesting: it appears that the mono-wall block is not better for the environment than a classical wall correctly insulated.

INTRODUCTION

Building industry must overcome nowadays major difficulties, mainly because of the very important energy consumptions and of the big quantities of raw materials employed to ensure hygro-thermal and acoustic properties or different other comfort facilities and economical performances to the main's building elements such the walls. The passive buildings seem to be the right solution permitting a decrease or at least a stabilization of climate change indicators. Because of this fact, national legislations are progressively introduced; the lecturer can refer for more details to the so called RT 2012 French thermal regulation [START_REF]relatif aux caractéristiques thermiques et à la performance énergétique des constructions[END_REF].

Note also that the requirements and the legislation relative to the building's envelope performance vary from a country to another and they are all very restrictive. The construction systems must consequently be adapted to each specific decision or design context.

The Life Cycle Assessment methodology, which is a standardized approach [START_REF]Environmental management -Life cycle assessment -Principles and framework[END_REF][START_REF]Environmental management -Life cycle assessment -Requirements and guidelines[END_REF] allows the quantification of the environmental footprint for goods, services and processes. It is then possible to identify some main points allowing the diminution of the environmental impacts since their earlier design stage. This methodology is regularly used in the field of building industry, for example in [START_REF] Frenette | Using life cycle assessment to derive an environmental index for light-frame wood wall assemblies[END_REF]. We have recently used the LCA methodology for performing some case studies concerning walls; it has been shown that the "common sense" generally leads the user to do very bad choices in the materials selection, from an environmental point of view [START_REF] Gruescu | Environmental footprint of a wall assembly by Life Cycle Assessment[END_REF][START_REF] Menet | A comparative life cycle assessment of exterior walls constructed using natural insulation materials[END_REF][START_REF] Menet | Environmental footprint of walls using natural materials coming from agricultural wastes[END_REF][START_REF] Menet | Environmental footprint of building elements using Life Cycle Analysis methodology, 21è[END_REF]. Among different solutions which are proposed to the builders, one consists to use adapted building elements, such the specific bricks called mono-wall blocks, which have the particularity to ensure both the function of a brick (or a breeze block), and the one of insulation. Besides, this permits a time gain because these elements are easy to implement. In summary, this element is a brick made of terra-cotta and has an alveolus structure material which ensures the insulation (figure 1). In the present work, two walls are considered: a classical wall made of breeze blocks and a wall built with mono-wall blocks. This paper presents the results of a LCA study made on the walls following the standardized methodology reported on figure 2.

GOAL AND SCOPE DEFINITION

The goal of the present study is the comparison of environmental impacts of the two considered walls.

The study is made by considering 1m 2 of wall insuring a thermal insulation resistance of 2.6 K/W. The time life of the wall is supposed to be 100 years (usual replacement time considered in the building area in France), but the environmental footprint is calculated only for 10 years. The joins are taken into account. The different tools used for the building are not taken into account. In the two cases, the blocks are supposed to be produced at a distance of 300 km from the distributor, which is supposed to be 10 km far from the user. Only the wood pallets are considered for the transportation of the blocks. All the steps of the LCA methodology are considered. The packaging of the glass wool is not considered.

INVENTORY ANALYSIS

Only the main data are reported in the following. More details about the inventory analysis can be found in ref [START_REF] Ghillebaert | Manson Analyse de cycle de vie Parpaing & Brique monomur[END_REF]. The parameters used for the two studied walls are presented in the Table 1.

Wall made of mono-wall blocks

Every block weights 16.5 kg, and has the following dimensions: 30cm x 30 cm x 21.2 cm. No vertical join is necessary for the building of mono-wall blocks. The horizontal join is supposed to be 1mm thick; table 1 shows the main characteristics of the mortar used forit.

Wall made of breeze blocks

Every block weights 17 kg and has the following dimensions: 20cm x 20 cm x 50 cm. The vertical and horizontal joins are about 1.5 cm thick. In Table 1 are presented the main characteristics of the mortar used for the join. To achieve the same thermal behavior, a glass wool panel must be added to the present wall (thickness: 80mm).

IMPACT ASSESSMENT

In order to evaluate the environmental impacts from the inventory analysis, the Eco-invent 2.0 data-base is used [10]. The environmental impacts were estimated using midpoint indicators associated with the CML method [11]. Since each indicator has a different scale and unit of measurement, the relative impacts are normalized, according to the specifications of ISO procedures [START_REF]Environmental management -Life cycle assessment -Principles and framework[END_REF][START_REF]Environmental management -Life cycle assessment -Requirements and guidelines[END_REF]. Each indicator is then expressed in points (see Figures 3 and4). A point represents the potential impact for a given indicator divided by the value of the same impact for a mean European people during one day. The chosen indicators and the corresponding units are reported in table 2.

CONCLUSION (INTERPRETATION PHASE)

The comparison of the two walls leads to two main conclusions: -In both cases, the production phase leads to the more important environmental impacts than the classical wall, except for aquatic ecotoxicity, which is mainly due to the mortar used for the joins; -Considering the whole life cycle, there is no significant difference in the environmental impacts whatever the indicator. Consequently, because of the errors margins, it is in fact impossible to affirm that the wall made with mono-wall blocks is better for the environment than a classical wall correctly insulated. This main conclusion seems to show equally that the systematic choice of monowall blocks for new buildings is due to the prices and the easiness of implementation rather than the environmental characteristics of this material.

This fact is the more important, since it is easily possible to replace glass wool used in the conventional configuration by natural insulation materials, which will lead to an environmental footprint lower to the one of the wall made with mono-wall blocks.

Figure 1 :

 1 Figure 1: A typical mono-wall brick Figure 2: Different phases of the LCA methodology [2, 3]

Figure 3 :

 3 Figure 3: Environmental impacts for each life phase (a) : breeze blocks and glass wool (b) mono-wall blocks

Table 1 :

 1 Main characteristic of the mortar and of the considered wall area

	Properties	Mono-wall block	Breeze block
	Density (Kg/m 3 )	1300	1300
	Height of the wall area (m)	1.065	1.05
	Lenght of the wall area (m)	0.939	0.952
	Volume of join (m 3 )	1.41 . 10 -3	1.41 . 10 -2
	Mass of joins (kg)	1,831	24.9

Table 2 -

 2 Chosen indicators for the calculation of the environmental impacts

	Letter	Indicator	Unit
	A	non-renewable energy consumption	MJ eq.
	B	resources depletion	kg Sb eq.
	C	100 year Global Warning Potential	kg CO 2 eq.
	D	Acidification	kg SO 2 eq.
	E	Eutrophication	kg PO4 2-eq.
	F	photochemical pollution	kg C 2 H 4 eq.
	G	aquatic toxicity	kg 1.4-DB eq.
	H	human ecotoxicity	kg 1.4-DB eq.

ACKNOWLEDGEMENTS

This study has been partly made using the free Bilan Produit ® software [12] jointly developed by the French Agency for the Environment and the Energy (ADEME) and the Cergy-Pontoise University.

The authors thank the students C. Ghillebaert, A. Doulcet and S. Manson which have worked on the subject during their university projects [9].