J.-M André

Ph Jonnard

Effective mass of photons in one dimensional photonic crystal

Keywords: numbers: 78, 67, Bf, 78, 67, Pt

Photons propagating within photonic crystals can be regarded as quasi-particles having an effective mass. Methods of calculation of this mass from the photonic density of mode are presented for the case of one-dimensional photonic crystal, both from infinite and finite structure. The model is implemented to calculate the influence of the number of unit cells and of absorption, for bilayered Bragg structure used as reflectors in the soft-x-ray domain. The physical meaning of the effective photon mass concept is discussed.

Introduction

Photons propagating within photonic crystals (PCs) [START_REF] Joannopoulos | Photonic Crystals, Molding the flow of light[END_REF] can be regarded as quasi-particles having an effective mass similar to the one of electrons in a crystal. The nature of this effective mass can appear somewhat scheming and it is exciting to imagine experiments revealing this effective photon mass. Gorelik has studied the effective photon mass in globular photonic crystals [START_REF] Gorelik | [END_REF] and more recently, reported experimental results establishing the existence of heavy photons in PCs [START_REF] Gorelik | Effective mass of photons and the existence of heavy photons in photonic crystals[END_REF]. The concept of effective photon mass appears also in the microcavities [START_REF] Skolnick | Strong coupling phenomena in quantum microcavity structures[END_REF] such as these ones implemented to obtain the Bose-Einstein condensation of polaritons [START_REF] Levrat | Condensation phase diagram of cavity polaritons in Ga-N based microcavities: Experiment and Theory[END_REF]. The calculation of effective photon mass can be performed from band structure computations by means of codes using expansion in Fourier series of the dielectric function; these codes are available in the literature [START_REF] Sukhoivanov | Photonic crystals[END_REF]. In principle such calculations assume infinite structure. Extension of the coupled-mode theory applied to deep optical grating should also provide the group velocity and the group-velocity dispersion, that is the inverse effective mass of photon, in one-dimensional photonic crystal [START_REF] De Sterke | Coupled-mode theory for light propagation through deep nonlinear gratings[END_REF].

In this communication we focus on calculations of the effective photon mass in onedimensional crystal (1D-PC), especially Bragg reflector, based upon the calculation of the photonic density of mode (DOM). We consider the case of infinite and finite structure. The paper is organized as follows: section 2 is devoted to the formalism of the effective photon mass (epm) in PCs and its relationship with the photonic DOM. Section 3 presents the calculation of the epm in an infinite 1D-PC and in section 4, we consider the case of a finite 1D-PC. Numerical examples are given in section 5 for the case of bilayered Bragg structure used as reflectors in the soft-x-ray domain. Finally in section 6, the physical foundation of the concept of epm is discussed and perspectives are presented.

Effective mass of photons in PCs and density of mode

The analogy between electromagnetism in a PC and quantum mechanics of solid-state physics has been reviewed in details in [START_REF] Joannopoulos | Photonic Crystals, Molding the flow of light[END_REF]. The propagation of a photon in a PC is governed by the following master equation deduced from Maxwell equations:

! ˆ " H(r) = # H(r) (1)
where ! H(r) is the spatial component of the harmonic magnetic field ! H(r,") = H(r) exp #i" t

(

) .

! ˆ " stands for the following linear operator:

! ! ˆ " = # $ 1 %(r) # $ & ' ( ) * + (2)
This operator for a PC is periodic since the dielectric constant ε(r) has the periodicities of the PC lattice. The equation ( 1) forms a traditional eigenvalue problem similar to the Schrödinger equation that the electron satisfies in a crystalline solid. The eigenvalue Ω is analogous to the electron energy and is proportional to the squared frequencies ω of the harmonic modes:

! " = # c $ % & ' ( ) 2 (3)
By virtue of the periodicity of ! ˆ " , the Bloch-Floquet theorem applies to Eq.( 1). The solutions are Bloch waves:

!

H(r) = h K (r) exp(iK.r) ( 4
)
where K is the so-called Bloch wave-vector and the quantities ! h K (r) have the periodicities of the PC lattice. The eigenvalue Ω is a function of the Bloch wave-vector K, so that the mode frequency ω depends on K:

ω = ω(K) (5).
This relationship called dispersion equation, can be explicitly found by inserting Eq.(4) in the master equation.

To the Bloch vector K, one can associate a quasi-impulse P by the relationship:

! P = ! K (6)
This quantity is important for the photon propagation in the PC because it satisfies the law of momentum conservation. Similarly by virtue of Planck's equation, one can find the photon energy ! E = ! "(K) which is also a function of the quasi-impulse:

E = E(P) (7)
This photon energy defined by Eq.( 7) satisfies the law of energy conservation. From Eq.( 7) one can introduce the concept of effective photon mass. If P 0 is an extremum point of E(P), one can expand this function in Taylor series:

! E(P) = E(P 0 ) + d t E(P) dP t " # $ % & ' P 0 (P ( P 0 ) t t ! t=1 ) * (8)
For t = 1, the corresponding term is a vector and for t=2, the corresponding term is a 2-rank tensor.

Since by assumption P 0 is an extremum, the derivative of rank 1 vanishes so that, in terms of quasiimpulse, one has the expansion:

! E(P) = E(P 0 ) + d 2 E dP 2 (P " P 0 ) 2 2 + ... (9)
From Eq.( 9) one introduces the effective mass tensor m* by:

! m* = { d 2 E dP 2 } "1 (10)
where {} -1 denotes the inverse tensor. Provided that the higher terms of the expansion are negligible, the dispersion function E(P) takes the parabolic shape:

!

E(P) = E 0 + 1 m* (P " P 0 ) 2 2 with E 0 = E(P 0 ) (11)
The first term E 0 of the expansion is nothing else but the photon energy of a "free" photon given by Planck's law.

For a 1D-PC, the Bloch vector and the quasi-impulse become mere scalar quantities so that the effective photon mass is also a scalar. One writes:

! m * = 1/ d 2 E dP 2 = ! / d 2 " dK 2 (12)
In the following we shall only consider 1D-PC. In this case, it is convenient to introduce the concept of photonic DOM ρ(ω) defined by

! " # ( ) = dK d# (13)
where K stands for the Bloch wave number. This quantity which is analogous to the inverse of a group velocity, can be calculated by different way as shown in the following sections.

Simple calculus shows that the scalar epm m* can be expressed by means of the DOM from the following relationship:

! m * = ! " 3 " ' (14)
where the prime denotes derivative with respect to ω.

Case of an infinite 1D-PC

An 1D-PC has a periodic structure of period d, along one axis (say z-axis), and consequently the dielectric constant ε(z) of an infinite 1D-PC can be expanded in Fourier series:

! "(z) = " l l= #$ $ % exp i lgz ( ) (15)
where g= 2π/d, d being the period of the PC along z. It has been shown [START_REF] André | Onde électromagnétique sinusoïdale dans un milieu périodiquement stratifié[END_REF] that in this case, the dispersion equation ω = ω(K) can be expressed by means of an infinite determinant as follows:

! sin 2 K d 2 " # $ % & ' = ( 0 ( ) sin 2 ) c J 0 d 2 " # $ % & ' (16)
where Δ(K) is the following infinite determinant:

! "(K) = ... ... 1 J #1 k 2 J 0 k 2 # K # g ( ) 2 J 1 k 2 J 0 k 2 # K 2 1 J #2 k 2 J 0 k 2 # K # g ( ) 2 J #3 k 2 J 0 k 2 # K # g ( ) 2 J #1 k 2 J 0 k 2 # K 2 J #2 k 2 J 0 k 2 # K 2 J 2 k 2 J 0 k 2 # K + g ( ) 2 J 1 k 2 J 0 k 2 # K + g ( ) 2 ... J q k 2 J 0 k 2 # K + q g ( ) 2 1 J #1 k 2 J 0 k 2 # K + g ( ) 2 ... 1 ... ... (17)
and J l are the components of the Fourier expansion of

! J(z) = k z (z) k where ! k z (z)
is the modulus of the normal component of the wave-vector associated with the two counter-propagating plane waves propagating in the whole structures [START_REF] André | Onde électromagnétique sinusoïdale dans un milieu périodiquement stratifié[END_REF][START_REF] Pardo | X-UV synthetic interference mirrors: theoretical approach[END_REF].

Let us now consider a mere Bragg bilayered structure consisting in N bilayers made up with two layers 1 and 2 of thickness d 1 and d 2 and dielectric constant ε 1 and ε 2 , respectively. If the structure is infinite, then N tends towards infinity and in this case, the dispersion equation K(ω) is given by [START_REF] Pardo | X-UV synthetic interference mirrors: theoretical approach[END_REF][START_REF] Yeh | Optical Waves in Layered Media[END_REF]:

! cos(Kd) = 1 1" r(#) 2 cos$ " r(#) 2 cos% ( ) (18)
with d = d 1 + d 2 , r(ω) is the Fresnel coefficient of reflection at the interface 1,2; Φ and Δ are the sum and the difference of the phase in the medium 1 and 2, respectively:

! " = k # 1 d 1 + k # 2 d 2 ! " = k # 1 d 1 $ k # 2 d 2 (19)
where ! k " j are the modulus of the normal component of the wave-vector associated with the two counter-propagating plane waves propagating in the layer j=1, 2. From Eq.( 16) and (17), one can calculate the DOM for an infinite 1D-PC and then by means of Eq.( 14), one deduces the epm.

Let us emphasize that the Eq.( 16) makes it possible to deal with the problem of any 1D-PC while Eq.( 18) is restricted to the case of bilayered structure.

Case of a finite 1D-PC

We consider now the case of a 1D-PC made with N unit cells. The unit cell can be any multilayered stack although the most common stack are bilayered or trilayered. An efficient strategy to calculate the epm in this case is to compute the DOM using the method proposed by Bendickson et al. [START_REF] Bendickson | Analytical expressions of the electromagnetic mode density in finite, one-dimensional, photonic band-gap structures[END_REF] generalized to lossly media in [START_REF] André | Stokes reciprocity equations and density of modes for absorbing stratified media[END_REF]. This strategy is based upon the use of the transmission coefficient t N of the whole structure:

! t N = T N exp(i" N ) (20)
T N being the transmittance and j N being the total phase accumulated as the radiation goes through the N unit cells of the PC. Hence one can introduce the effective wave-number k N and can write

! " N = k N d (21)
Following [START_REF] Pardo | X-UV synthetic interference mirrors: theoretical approach[END_REF][START_REF] Yeh | Optical Waves in Layered Media[END_REF], we define the N-unit cell DOM as:

! " N = dk N d# = 1 d d$ N d# (22) Simple calculations show that ! " N = 1 d z N ' 1+ z N 2 with z N = Arctan(# N ) (23)
the prime standing for the derivate with respect to ω. It has been shown in [START_REF] Bendickson | Analytical expressions of the electromagnetic mode density in finite, one-dimensional, photonic band-gap structures[END_REF][START_REF] André | Stokes reciprocity equations and density of modes for absorbing stratified media[END_REF] how to calculate ρ N , then ϕ N using Eq.( 23), by means of the transfer matrix method.

Let us outline that ρ N differs from r but that these two quantities can be related as shown hereafter. To obtain the photonic DOM ρ defined by Eq.( 13) from ρ N , one makes use of:

! " = " N dK dk N (24)
K being the Bloch wave-number of an infinite 1D-PC introduced in section 3. To get the derivative

! dK dk N
, one can use the following formula [START_REF] Bendickson | Analytical expressions of the electromagnetic mode density in finite, one-dimensional, photonic band-gap structures[END_REF]:

! z N = arctan(k N d) = z tan N K ( ) cot g K ( ) (25)
with z = tan (ϕ), ϕ being the phase accumulated through a unit cell.

Some manipulations give:

! dK dk N = z cos(dk N ) N cot g dK ( ) cos 2 dK N ( ) " tan dK N ( ) sin 2 dK ( ) # $ % & ' ( (26)
so that the photonic DOM reads:

! " = " N cot g(#) cos 2 (# N ) N cot g dK ( ) cos 2 dK N ( ) $ tan dK N ( ) sin 2 dK ( ) % & ' ( ) * (27)
From Eq.( 7) and Eq.( 14), it is possible to calculate the epm m* N of a N-unit cell 1D-PC.

Numerical results

In this section we calculate the epm m* for a periodic stack of N bilayers forming a Bragg reflector for the soft-x-ray domain. We consider a mirror consisting of unit cells made up with a Mo layer of thickness d 1 = 1.41 nm and of B 4 C of thickness d 2 = 1.976 nm designed to reflect at the first Bragg order a radiation of vacuum wavelength 6.77 nm (or photon energy 182 eV) at an incident angle of 30° corresponding to the first diffraction order. The materials are supposed to be absorbing and their complex optical indices are taken from [13]. This kind of device is now standard for x-ray optics [START_REF] Pardo | X-UV synthetic interference mirrors: theoretical approach[END_REF]. We first consider the infinite case that is, N tends towards infinity. In this case the dispersion equation is given by Eq.(18). Figure 1 shows the normalized photonic DOM ! " versus photon energy E between 180 and 210 eV. Note that the normalization consists in choosing ! " equal to unity outside the Bragg domain (180 eV for instance). Figure 2 presents on the same energy domain the epm m* given in kg. The epm is cancelled in the Bragg domain, is negative in the low energy domain before the Bragg gap and is positive in the high energy domain after the gap. This kind of behaviour was also reported by Gorelnik [START_REF] Gorelik | Effective mass of photons and the existence of heavy photons in photonic crystals[END_REF] who used the same theoretical approach. We know consider the case of finite structure. Figures 345show the normalized DOM for N=20, 200 and 2000 versus photon energy E and figures 6-8 present the corresponding variation of the epm m*. First one notes that for large values of N (say 2000 in our example) our model of finite structure is in fair agreement with the infinite model. For small values of N, the infinite model fails and the shape of the curve m* versus E changes considerably with the value of N. At intermediate (realistic) values of N that is typically between 100 and 500, oscillations appear outside the Bragg gap which are due interferences within the multilayer structures. For very low numbers of bilayers (see figure 6, N=20), the epm curve presents "spikes" which correspond to the zeros of the DOM derivatives. The influence of photon absorption is illustrated in Figures 9 and 10 by the case of a reflector with a large number of bilayers (N =2000) without absorption; to deal with this case, calculations are carried out with the imaginary part of the optical index enforced to zero. The transparency of the Bragg structure makes it possible to obtain strong interferences which are not attenuated by the photon attenuation. These interferences give rise to oscillations in the vicinity of the Bragg domain outside the bandgap both in the DOM and the epm curves. These oscillations are so closed that they are not resolved in the figures in spite of the huge number of points (10 5 ) retained in the calculation. Let us mention that for an infinite structure these oscillations disappear and one obtains an epm quite similar to the one presented in figure 2.

Discussion and conclusion

It has been shown that for finite 1D-PC, accurate calculations of the epm required a model which takes into account the limited number of unit cells and the absorption in the media. This paper gives for the first time such a model. The Bendickson et al. model generalized to lossy media [START_REF] Bendickson | Analytical expressions of the electromagnetic mode density in finite, one-dimensional, photonic band-gap structures[END_REF][START_REF] André | Stokes reciprocity equations and density of modes for absorbing stratified media[END_REF] appears to be an efficient method of calculation of m* via the photonic DOM. We have applied it to pure periodic structure but it is important to mention that it works well for quasi-periodic one such as Fabry-Perot resonator where defect states can occur [START_REF] Joannopoulos | Photonic Crystals, Molding the flow of light[END_REF][START_REF] André | Stokes reciprocity equations and density of modes for absorbing stratified media[END_REF]. This case will be considered in a forthcoming paper.

The concept of epm is likely very useful when considering scattering processes involving photons propagating in PCs since it should be introduced in the law of energy conservation. Among the involved scattering processes, there is the Raman scattering in microcavity. Optical-phonon Raman scattering enhancement has been reported in semiconductor planar cavity [START_REF] Fainstein | Raman scattering enhancement by optical confinement in a semiconductor planar microcavity[END_REF] and Bragg-Raman scattering has been observed in the soft-x-ray domain with a La/B 4 C Bragg superlattice [START_REF] André | J X-ray Raman scattering with Bragg diffraction in a La-based superlattice[END_REF]. The main features of the latter process have been accounted for by considering the conservation of momentum and energy but the outgoing Raman photon was consider as a "massless" photon. A calculation done with the introduction of the epm of the outgoing Raman photon in the value of its energy (Eq.( 8)) is under way. We guess that this study can be a test for the validity of the epm concept.

Figure 1 :

 1 Figure 1: Normalized density of mode for the Mo/B 4 C Bragg reflector versus photon energy E with an infinite number of bilayers. The details of the Bragg reflectors are given in the text.

Figure 2 :

 2 Figure 2: Effective photon mass for the Mo/B 4 C Bragg reflector versus photon energy with an infinite number of bilayers. The details of the Bragg reflectors are given in the text.

Figure 3 :

 3 Figure 3: Normalized density of mode for the Mo/B 4 C Bragg reflector versus photon energy with a finite number of bilayers: N=20. The details of the Bragg reflectors are given in the text.

Figure 4 :

 4 Figure 4: Normalized density of mode for the Mo/B 4 C Bragg reflector versus photon energy with a finite number of bilayers: N=200. The details of the Bragg reflectors are given in the text.

Figure 5 :

 5 Figure 5: Normalized density of mode for the Mo/B 4 C Bragg reflector versus photon energy with a finite number of bilayers: N=2000. The details of the Bragg reflectors are given in the text.

Figure 6 :

 6 Figure 6: Effective photon mass for the Mo/B 4 C Bragg reflector versus photon energy with a finite number of bilayers: N=20. The details of the Bragg reflectors are given in the text.

Figure 7 :

 7 Figure 7: Effective photon mass for the Mo/B 4 C Bragg reflector versus photon energy with a finite number of bilayers: N=200. The details of the Bragg reflectors are given in the text.

Figure 8 :

 8 Figure 8: Effective photon mass for the Mo/B 4 C Bragg reflector versus photon energy with a finite number of bilayers: N=2000. The details of the Bragg reflectors are given in the text.

Figure 9 :

 9 Figure 9: Normalized density of mode for the Mo/B 4 C Bragg reflector versus photon energy with a finite but large number of bilayers: N=2000 without absorption. The details of the Bragg reflectors are given in the text.

Figure 10 :

 10 Figure 10: Effective photon mass for the Mo/B 4 C Bragg reflector versus photon energy with a finite but large number of bilayers: N=2000 without absorption. The details of the Bragg reflectors are given in the text.


F Mottez

V Génot

Electron acceleration by an Alfvénic pulse 1 propagating in an auroral plasma cavity 2

With the help of a 2.5D Particle In Cell (PIC) simulation code, we investigate the physics of the acceleration of auroral electrons, through the interaction of an isolated Alfvén wave packet with a plasma density cavity. The cavity is edged by density gradients perpendicular to the magnetic field. We show that a single passing of an isolated wave packet over a (infinite) cavity creates an electron beam. It triggers local current and beamplasma instabilities, and small scale coherent electric structures. The energy flux of downgoing electrons is significantly increased, whereas upgoing electrons are also accelerated, even if no beam is formed. Accelerated electrons remain after the passage of the Alfvénic pulse, allowing the observation of energetic particles without any significant electromagnetic perturbation. The dependence of this process on the electron to ion mass ratio is consistent with its control by inertial effects.

Introduction

Early models of auroral particle acceleration were based on quasi-static parallel electric fields, accelerating electrons and ions in opposite directions, with roughly the same narrow distribution in energy [START_REF] Block | The role of magnetic-field-aligned electric fields in auroral acceleration[END_REF]. Other evidences of accelerated electrons, with broader distributions in energy and direction, have shown the importance of time-varying electric fields [START_REF] Hultqvist | Simultaneous observation of upward moving field-aligned energetic electrons and ions on auroral zone field lines[END_REF], possibly carried by Alfvén waves [START_REF] Kletzing | Electron acceleration by kinetic Alfven waves[END_REF]. Deep plasma cavities above the Earth auroral zone are a privileged place for electron acceleration [START_REF] Hilgers | Probe measurements of low plasma densities with applications to the auroral acceleration region and auroral kilometric radiation sources[END_REF] and the subsequent turbulence, characterized by electrostatic coherent structures such as double layers and solitary waves [START_REF] Bostrom | Characteristics of solitary waves and weak double layers in the magnetospheric plasma[END_REF][START_REF] Eriksson | A statistical survey of auroral solitary waves and weak double layers 2. Measurement accuracy and ambient plasma density[END_REF]. With the Freja spacecraft, is was shown that these regions are pervaded by Alfvén wave packets (termed as Solitary Kinetic Alfvén Waves, SKAW).

Measurements have shown that these waves carry a Poynting flux large enough to accelerate electrons at auroral energies if this flux is to be efficiently dissipated through a relevant mechanism [START_REF] Louarn | Observation of kinetic Alfven waves by the Freja spacecraft[END_REF]Volwerk et al., 1996]. Coincident satellite measurements of fields and particles demonstrate that, as functions of increasing auroral activity, 25-39% of the total electron energy deposited in the ionosphere may be attributed to the action of Alfvén waves [START_REF] Chaston | How important are dispersive Alfven waves for auroral particle acceleration?[END_REF]. [START_REF] Hull | Time development of field-aligned currents, potential drops, and plasma associated with an auroral poleward boundary intensification[END_REF] performed a case study of the development of an acceleration region based on multi-point Cluster observations in the high altitude auroral zone, across the plasma sheet and into the polar cap. They first identify an Alfvén wave dominated system, with its typical broad spectrum of accelerated particles, in the vicinity of a plasma cavity. After Alfvénic acceleration the signature (the so-called inverted V) of particles accelerated by a quasi-static structure, such as a strong double layer is observed. This suggest that the Alfvénic acceleration process could act as a precursor to the quasi-static process. These phenomenons were associated to a poleward boundary intensification of the auroral arcs seen by an ultraviolet camera. [START_REF] Hull | Time development of field-aligned currents, potential drops, and plasma associated with an auroral poleward boundary intensification[END_REF] performed an analysis of the small-scale Alfvénic currents and noticed that, in the inertial dispersive range, they are damped. In their conclusion, the authors suggest that this might be due to Landau damping, near and above the high altitude acceleration region. In the present paper we propose a different interpretation: the current attenuation would be the result of a transfer of energy from the wave to the electrons, with the local heating of electrons, and formation of accelerated electron beams, associated to the generation of a parallel (to the ambient magnetic field) electric field.

There is a precise reason why this Alfvénic current dissipation through electron acceleration would occur at the inertial scale. According to the MHD theory, an Alfvén wave could not accelerate particles along the ambient magnetic field ⃗ B 0 because the parallel electric field (E ∥ ) is null. Within the MHD framework, wavelengths are long compared to the inertial length c/ω pe and to the ion Larmor radius ρ i , i.e. k ⊥ c/ω pe << 1 or k ⊥ ρ i << 1,

where k ⊥ is the perpendicular (to ⃗ B) wave vector, ω pe is the electron plasma frequency, c is the speed of light, and ρ i is the ion Larmor radius (ρ i = m i v ti /eB, m i is the ion mass). Beyond the MHD approximation, there exists two regimes of parameters where

Alfvén waves can carry a parallel electric field. This is when k ⊥ c/ω pe ∼ 1 or k ⊥ ρ i ∼ 1 [START_REF] Goertz | Kinetic Alfvén waves on auroral field lines[END_REF]. In the Earth auroral zone, the plasma beta β = 2µ 0 p/B 2 0 ≪ m e /m i , and one deduces that

ρ i ω pe /c = (βm i /m e ) 1/2 < 1; the condition k ⊥ c/ω pe ∼ 1 is then reached before k ⊥ ρ i ∼ 1.
At such scales, an Alfvén wave carries a parallel electric field, given for a plane wave by the relation:

D R A F T May 25, 2011, 5:48am D R A F T MOTTEZ & G ÉNOT: ALFV ÉN PULSE IN CAVITY X -5 E ∥ E ⊥ ∼ k ∥ k ⊥ ( c ωpe k ⊥ ) 2 1 + ( c ωpe k ⊥ ) 2
(1)

This parallel electric field has a functional form similar to those of the other fields, i.e.

sinusoidal, with the same wavelength.

Therefore, most of the theories of auroral acceleration by waves are based on inertial Alfvén waves (see for instance [START_REF] Kletzing | Electron acceleration by kinetic Alfven waves[END_REF]; [START_REF] Lysak | On the kinetic dispersion relation for shear Alfven waves[END_REF]; Watt et al. [2008]).

Bridging with the early models cited at the beginning of the introduction, we can see now that the time varying electric fields associated to the non stationary acceleration structures [START_REF] Hultqvist | Simultaneous observation of upward moving field-aligned energetic electrons and ions on auroral zone field lines[END_REF] are, in many cases, identified as inertial Alfvén wave electric fields.

Before concluding that such waves can accelerate electrons, one must question the origin of the oblique component k ⊥ of the wave vector. An explanation based on the properties of plasma cavities has been proposed [START_REF] Génot | A study of the propagation of Alfven waves in the auroral density cavities[END_REF]. When an Alfvén wave, initially in pure parallel propagation (k ⊥ = 0), propagates upon a perpendicular density gradient, the wave front is bent (phase mixing) as propagation is faster in the low density region than the denser one. In that case, there is a parallel electric field quite similar to the one of an inertial Alfvén wave: in the above formula, k ⊥ may be replaced by ∂ x ln n (n is the plasma density, and ∂ x is a spatial derivative perpendicularly to the magnetic field). The role of density cavities on the formation of narrow-scale Alfvén waves has been confirmed recently with a three dimensional linear model of the auroral flux tube [START_REF] Lysak | Propagation of kinetic Alfven waves in the ionospheric Alfven resonator in the presence of density cavities[END_REF].

Therefore the parallel electric field could be partly explained by a propagation effect of Alfvén waves on the borders of the plasma cavities, where the perpendicular density gradients are large.

D R A F T May 25, 2011, 5:48am D R A F T X -6 MOTTEZ & G ÉNOT: ALFV ÉN PULSE IN CAVITY
A set of numerical simulations showed that parallel electric fields first develop on the (rather long) Alfvén wavelength and accelerates strongly a minority of electrons [START_REF] Génot | Electron acceleration by Alfven waves in density cavities[END_REF]. Then an unstable electron beam is formed which triggers current and beam instabilities. The non-linear evolution shows the formation of electron holes structures and associated coherent (small scale) electrostatic structures observed in auroral cavities by many spacecraft [Génot et al., 2001[START_REF] Génot | Alfvén wave interaction with inhomogeneous plasmas: acceleration and energy cascade towards smallscales[END_REF]. More recently a comparative study conducted with FAST data [START_REF] Chaston | Ionospheric erosion by Alfvén waves[END_REF] revealed that predictions of the present model were actually observed (electron acceleration on the density gradients, wave focusing in the cavity, ...). On the simulation side, analog results were found independently by Tsiklauri et al. [2005]. However these simulations were initialized with a sinusoidal AW (one wavelength covering the parallel length of the simulation domain), whereas observed SKAW propagate as isolated wave packets (see for instance observations by the Freja spacecraft in [START_REF] Louarn | Observation of kinetic Alfven waves by the Freja spacecraft[END_REF]; Volwerk et al. [1996]). To obtain a quantification of the acceleration process according to the previous scenario for more realistic conditions, it is therefore crucial to wonder whether a localized input (the Alfvén pulse) may still trigger sufficient acceleration to power the aurora, i.e. whether energy transfer from the wave to the electrons may still take place over a reduced distance along field lines. The present paper is addressing this broad question by analysing the propagation of Alfvén wave packets upon a density gradient. Investigating the process at work in such a configuration with the help of a kinetic code which retains non-linearities (for a description of the PIC code see [START_REF] Mottez | A new guiding centre PIC scheme for electromagnetic highly magnetized plasma simulation[END_REF]), we shall set a bridge between the work of [START_REF] Génot | A study of the propagation of Alfven waves in the auroral density cavities[END_REF] which included a localized Alfvénic pulse, but considered only the bi-fluid

D R A F T May 25, 2011, 5:48am D R A F T X -7
approach (via linearized equations), and the self-consistent work of [START_REF] Génot | Alfvén wave interaction with inhomogeneous plasmas: acceleration and energy cascade towards smallscales[END_REF] where Alfvén waves were purely sinusoidal.

The numerical method and the sets of simulation parameters are given in Section 2 whereas the pulse propagation is tested in Section 3. Electron and ion acceleration are evidenced in Section 4.1, followed in Section 4.2 by an analysis of the plasma turbulence generated by accelerated electrons. In Section 4.3, we determine which part of the wave spectrum contributes most efficiently to the acceleration. Section 5 is devoted to the discussion of our results compared to other similar, but contradicting, simulations [Tsiklauri, 2007], and we conclude in Section 6.

Numerical method and parameters

Following [START_REF] Génot | Electron acceleration by Alfven waves in density cavities[END_REF]Génot et al. [ , 2001[START_REF] Génot | Alfvén wave interaction with inhomogeneous plasmas: acceleration and energy cascade towards smallscales[END_REF] the numerical simulations make use of an electromagnetic PIC code that takes into account the motion of the electron guiding center and the full ion motion [START_REF] Mottez | A new guiding centre PIC scheme for electromagnetic highly magnetized plasma simulation[END_REF]].

The 2D simulation domain is defined by the rectangular coordinate system (x, y). The direction x is the direction of the ambient magnetic field B 0 and corresponds to the longest side of the simulation domain. The vectors are tri-dimensional, with a component z perpendicular to the simulation box.

The physical variables are reduced to dimensionless variables. Time (the inverse of) and frequencies are normalized by the electron plasma frequency ω p0 which corresponds to a reference background electron density n 0 . Velocities are normalized to the speed of light c, and the magnetic field is given in terms of the dimensionless electron girofrequency ω ce /ω p0 . The mass unit is the electron mass m e . Therefore, the units are c/ω p0 for distances, ω p0 /c for wave vectors, e for charges, en 0 for the charge density,

cω ce /ω p0 D R A F T May 25, 2011, 5:48am D R A F T
for the electric field, and ce/ω p0 for the magnetic moment µ of the electrons. In the following parts of the paper, all equations, numerical values and figures are expressed in this system of units. Let us note that the reference density n 0 (or plasma frequency ω p0 ) is a free parameter and may therefore be fixed arbitrarily; other plasma parameters are consequently deduced from this choice. For instance, setting n 0 = 10 cm -3 leads to B 0 = 4µT for a magnetization given by ω ce /ω p0 = 4 typical of the auroral region.

We present results from six simulations : run A to run F. Run A is set initially with an Alfvén wave packet which propagates upon an homogeneous plasma. The other simulations are set with Alfvén waves (a wave packet or a sinusoidal wave) which propagate in an infinite plasma cavity extended along the magnetic field direction, and delimited by two strong density gradients perpendicular to the magnetic field.

The simulations are initialized with the sum of 8 (or 1) sinusoidal waves with a maximum at x = 0. The waves have a right-hand circular polarization. In run A and run B, δB y (x) at t=0 is a sum of eight sinus of equal amplitude and δB z (x) the sum of eight cosines of the same amplitude, given by the wave magnetic field δB = 0.032 × B 0 where B 0 = 4.0 is the background magnetic field amplitude (or ω ce /ω p0 = 4). In run D and run E a single sinusoidal wave is initially set with an amplitude δB = 0.0904 × B 0 , such that the runs A, B, D, E have the same initial magnetic energy. (The sum Σ i δB 2 i has the same value in the four simulations, where i is the index of the monochromatic waves that contribute to the wave packet.) Setting the same magnetic energy in these simulations allows for a better comparison of the energy transfer from the waves to the particles, because, at least in the MHD approximation, the wave Poynting flux of AW,

S = S ∥ = δB 2 V A /µ 0 , is proportional to the magnetic energy density δB 2 . D R A F T May 25, 2011, 5:48am D R A F T MOTTEZ & G ÉNOT: ALFV ÉN PULSE IN CAVITY X -9
In run A and run B, wavelengths are λ = λ 0 /m where λ 0 = 409.6 and m varies from 1 to 8. Only the larger wavelength of run B is kept in run D, while the smallest is kept in run E. The phase velocities vary from 0.4529 for the shortest wave to 0.2257 for the longest, i.e. the dispersion is non-negligible. The ideal MHD Alfvén velocity is

V A = B/ √ µ 0 ρ = 0.
2 is smaller; it does not represent accurately the propagation of these waves for which inertial effects are effective. The polarization of the 8 sinusoidal waves is described in details in the appendix of [START_REF] Mottez | A guiding centre direct implicit scheme for magnetized plasma simulations[END_REF].

The size of the whole simulation domain is 4096∆x × 128∆y, where ∆x = ∆y = 0.1 is the size of the grid cells.

Mottez [2003] have shown, in the context of Vlasov-Maxwell formalism, that the gradients of the auroral plasma cavity can be modelized as tangential discontinuities, in spite of a non-scalar pressure tensor that cannot be characterized through the MHD theory [START_REF] Mottez | The pressure tensor in tangential equilibria[END_REF]. Run B contains a plasma cavity whose area of largest depth, n min /n max = 0.2 is a channel of infinite length, and 12∆y broad. The edges of the cavity are smoothed with a gaussian profile. The electron thermal velocity (in and outside the cavity) is v te = 0.1 (ion and electron temperatures are equal). The ion to electron mass ratio is varied from m i /m e = 100 (run C) to m i /m e = 400 (in runs A, B, D, E). There are ∼ 26 × 10 6 particles of each species corresponding to an average of 50 particles per cell. There are 2048 time steps in the simulation, defined by ∆t = 0.2, corresponding to a time lapse t max = 409.6.

Test bench for the pulse propagation

The pulse propagation in a uniform plasma has been tested with the simulation run A. Figure 1 shows the E y (x, y(x, t), t) transverse component of the electric field along the direction x, as a function of time for this simulation. Initially, the line defined by

D R A F T May 25, 2011, 5:48am D R A F T X -10 MOTTEZ & G ÉNOT: ALFV ÉN PULSE IN CAVITY
y(x, 0) = 5 is in the middle of the density gradient. Displacing it with the same velocity as the plasma allows to keep it inside the region of transverse density gradient. Therefore, considering the MHD approximation, we move this line with the velocity

⃗ v = ⃗ E × ⃗ B/B 2 . Practically, y(x, t) is the solution of d t y(x, t) = E z (x, y, t)/B x (x, y, t).
We can see that the superposition of the 8 waves is like a single Alfvénic pulse (and a residual short wavelength sinusoidal wave). It propagates with dispersion, but on the timescale of the simulation, we can still clearly identify a wave packet. (A longer simulation box and more computing time would allow for longer wavelengths, and a less dispersive Alfvén wave packet, as in MHD.) Initializing the wave packet with a sum of 8 waves with an amplitude distribution fitting a Gaussian shape, instead of a sum of 8 waves of equal amplitudes, would make a smoother and better localized wave packet. Nevertheless, wave packets observed on-board satellites, such as in [START_REF] Louarn | Observation of kinetic Alfven waves by the Freja spacecraft[END_REF], do not specifically exhibit a Gaussian shape.

In run A, the initial density is uniform. It is not the case for other simulations. The electron density in run B is shown on the Fig. 2, for three different times t i corresponding to the beginning, (approximately) the middle and the end of the simulation. To help the reader to understand the Fig. 3, we have drawn the line y(x, t i ) on Fig. 2. Figure 2 shows that the perpendicular displacement of the cavity maximizes where the wave packet amplitude also maximizes. This motion is magnified by the aspect ratio of the figure (corresponding to dimensions 409.6 × 12.8), it is actually less important than it looks on the figure. We can see on Fig. 2 that to a good approximation, the iso-density areas follow the two y(x, t) lines superimposed to the gradient. This proves that the motion of the cavity is mainly controlled by the ⃗ E × ⃗ B/B 2 drift of the plasma induced by the (Alfvénic

D R A F T May 25, 2011, 5:48am D R A F T X -11
pulse) electromagnetic field as expected. Is the cavity depth and profile modified by the Alfvénic pulse ? On Fig. 2, we can compare, for instance the cavity at time t = 0 (for any value of x since it is initially uniform in that direction), and at time t = 409.6. At this late time, the Alfvénic pulse is in the middle of the simulation box (around x = 200), and we can see that for x < 200, where the pulse has already passed, the cavity has the same depth and the same transverse shape (i.e. along the y axis) as at t = 0. In other words, the passage of the wave packet does not destroy the plasma cavity. This confirms the results shown in previous works, conducted with only one sinusoidal wave [Génot et al., 2001[START_REF] Génot | Alfvén wave interaction with inhomogeneous plasmas: acceleration and energy cascade towards smallscales[END_REF]. It is at odds with Sydorenko et al. [2008] in which numerical simulations of an Alfvén wave packet of high amplitude propagating in the ionospheric Alfvén resonator are conducted. Contrary to the present work, they initialized a wave with a transverse wave vector. The particle motion is computed under more restrictive simplifications as in our code, but, reciprocally, their simulations include a real mass ratio and longitudinal density gradients, where our simulations are initially homogeneous along the ambient magnetic field direction. As expected in the inertial Alfvén wave regime, acceleration occurs. This acceleration is accompanied, for small enough wavelength (its does not work with the fundamental mode), by the creation of a density depletion in a small area at low altitude, where the accelerating parallel electric field reaches its highest amplitude.

The density depletion is caused here by the accelerated plasma being expelled from the acceleration region (a depression of about 50% may be formed). In our simulation, the acceleration process does not destroy the cavity, but it does not dig it either. This may be a consequence of the homogeneity of the plasma in the longitudinal direction, that is also the direction of plasma acceleration. The significant difference between run A and run B concerns the parallel electric field.

In run A (not shown) it is absent. In run B, no parallel electric field is set initially, but it develops naturally, as a consequence of the phase mixing on the density gradients. A map of the parallel electric field E x (x, y) is shown on Fig. 5 for three different times. From an initial vanishing value, a localized structure develops as the Alfvénic pulse propagates.

The E x electric field is mainly localized in and around the density depletion, and its structure is time dependent. This evolution may also be traced on Fig. 6 that displays E x (x, y(x, t), t). At time zero, it is null everywhere, but a bipolar structure soon emerges with the characteristic size of the wave packet. In the linear regime, the scale of the parallel electric field is given by the characteristic scale of the incoming field, E y , as it can

D R A F T May 25, 2011, 5:48am D R A F T X -13
be seen on Fig. 3. After time t = 200, small scale structures appear. This is a signature of a non-linear evolution of the system. We can notice first a peak of parallel electric field, that propagate at the same speed as the wave packet, soon followed by a bunch of less intense structures that propagate slightly slower, and preceded by small amplitude fast structures. All these structures emerge from the wave packet. The parallel electric field is able to accelerate electrons along the ambient magnetic field (x direction), as can be seen on Fig. 7. This figure shows the electron distribution function f e (x, v x ) integrated over y. The figure is in log-scale, it provides a good way to show the creation of a minority of fast electrons, with velocities v x up to 8v te .

According to [START_REF] Semeter | Persistent quasiperiodic precipitation of suprathermal ambient electrons in decaying auroral arcs[END_REF], the intensity of an auroral arc is quasi-proportional to the flux of electron kinetic energy across a surface, defined as an integral over the velocity space

F (x) = ∫ v ′2 v ′ x f (⃗ v ′ , x, y ′ )d 3 ⃗ v ′ dy ′ . ( 2
)
To evaluate the efficiency of the acceleration process seen in the simulations, we have evaluated this flux. Practically, we compute a sum over the macro-particles,

F (x) = Σ S v 2 v x ( 3
)
where S is a subset of particles. Let x ′ , y ′ be their coordinates, S is restricted to the particles such as x min < x ′ < x max . In order to have sufficient statistics, but a reasonable spatial resolution, we have chosen occurs in the center of the box, the plasma cavity region. This confirms the role of the plasma cavity in the acceleration process.

x
The comparison of the upward and of the downward fluxes shows that the acceleration occurs in the two directions, but mainly in the direction of propagation of the Alfvénic pulse, i.e. downward. If this was a resonant process of acceleration, we would have only acceleration for particles with a velocity close to the Alfvénic pulse velocity (i.e. downward). The observation of an upward flux shows that the process involved here is not purely resonant. Nevertheless, the fact that downward acceleration is more efficient is somehow resonant, as it is due to the fact that the particles with an initially downward velocity see the accelerating electric field carried by the wave during a longer time than the upgoing particles.

Figure 9 is a similar plot for the proton energy fluxes. We can see that the ions are accelerated. The flux of the downward protons increases by a factor 1.6. It is 2.7 times less than for the electrons. Proton distribution functions (not shown) do not exhibit any ion beam; this is consistent with a weaker acceleration efficiency for the protons. On Fig. 9, we can also notice that the difference between upward and downward energy fluxes is smaller than for the electrons. This is due to the fact that the difference of upward and downward ion velocities is smaller than for the electrons. When it is compared to the Alfvén wave packet velocity, we have 2v te /v P ∼ 0.36 for the electrons, and 2v ti /v P ∼ 0.018 for the ions, where we have taken v P = 0.5 for the wave packet. Therefore, the difference of time spent by the upgoing and the downgoing ions in the acceleration region (that moves with the series of simulations, in the context of Solar physics, that presents similarities with the simulations presented here and in previous papers (see [START_REF] Mottez | Comment on "PIC simulations of circularly polarised Alfvén wave phase mixing: a new mechanism for electron acceleration in collisionless plasmas[END_REF] for a detailed comparison). In [Tsiklauri et al., 2005], plots of the ion distribution functions show a broadening that is not seen in our simulations. In Run A for instance, we can see only an increase by 5% of the ion perpendicular thermal speed. Actually, when Tsiklauri et al.

[2005] analyzed the contribution of this broadening to the energy budget, they found no contribution, because the kinetic energies in the directions y and z oscillate in anti-phase.

Therefore, they conclude that there is no ion acceleration in their simulations and that the ion distribution broadening is due to the usual velocity perturbations associated to the AW. Why is there a distribution broadening in their simulations, and not in ours ?

In [START_REF] Mottez | A guiding centre direct implicit scheme for magnetized plasma simulations[END_REF], the perturbation of the velocity fields associated to the AW is given, in the frame of the cold bi-fluid plasma theory. This is the perturbation used for the initial conditions of the simulations of the present study. The ion velocity perturbations (given in code units) are

V py (x) = Sc p1 B y (x) and V pz (x) = c p1 B z (x) (4)
where c p1 = -Sm e /(ω + Sω ci )kωm p where S = ±1 depends on the polarization (right or left-handed). We can see that for a weak ion to electron mass ratio (16 in [Tsiklauri et al., 2005]) this perturbation is strong, while for a larger mass ratio (400 in the present study), the ion velocity perturbation is lower. Therefore, we can conclude that there is no significant ion acceleration in the present simulations, and this result agrees with the analysis conducted by other authors.

Where are the accelerated electrons ? On Fig. 7 there is a peak of high energy electron if the Alfvén wave velocity is comparable to the electron velocity (here v A ∼ 3v te ). This does not happen in the high altitude auroral zone (10 000 km) as seen for instance from the Viking, Cluster or Polar spacecrafts, where v A is of the order of c, but it is possible at lowest altitudes (1000 km) in the regions that have been explored by Freja, for instance.

Coherent structures

We have seen on Fig. 6 that after time 200, the cascade from large scale to small scale parallel electric field structures becomes visible. This effect was already seen in previous

D R A F T May 25, 2011, 5:48am D R A F T X -18 MOTTEZ & G ÉNOT: ALFV ÉN PULSE IN CAVITY
studies, with a sinusoidal Alfvén wave [Génot et al., 2001, b]. When looking at Fig. 7, we can see that the small scale electric field structures are associated to vortices in the phase space. These vortices have a suprathermal velocity comparable to the wave packet velocity. Similarly to the analysis developed in [START_REF] Génot | Alfvén wave interaction with inhomogeneous plasmas: acceleration and energy cascade towards smallscales[END_REF], we infer that they are a consequence of the nonlinear evolution of an electron beam plasma instability. The beam is formed by the fastest electrons which quit the acceleration region (the packet) ahead of it. To assess the triggering of this instability we use the criterion derived by [START_REF] Gary | Electrostatic instabilities in plasmas with two electron components[END_REF]:

v the-beam v drif t-beam ≪ ( n e-beam n e
) 1/3

(5)

where v the-beam is the beam thermal velocity, v drif t-beam is the beam drift velocity, n e-beam is the beam density and n e is the core density. Using the electron distribution function at time 179.2 and x =102.4 (see the second panel of Figure 7), for which small scale structures are already present, we deduce the following parameters: v the-beam = 0.04, v drif t-beam = 0.69, n e-beam = 0.026, and n e = 1.21, which readily satisfy the above criterion.

Similarly to the small amplitude case studied in [START_REF] Génot | Alfvén wave interaction with inhomogeneous plasmas: acceleration and energy cascade towards smallscales[END_REF] we do not observe signatures of a Buneman instability. This instability is a consequence of a large velocity drift of the electron distribution with respect to the ion one and is triggered in the presence of a parallel electric field of sufficiently large amplitude and installed on long spatial scale.

This does not happen with a localized pulse of reasonable amplitude.

Efficiency

D R A F T May 25, 2011, 5:48am D R A F T X -19
The wave packet is built initially as a sum of sinusoidal waves of various wavelengths with the same amplitude. This choice is arbitrary but facilitates the study of each wave mode efficiency in the acceleration process.

Run D was initialized with a single sinusoidal wave with the largest wavelength. Its amplitude was chosen in order to have the same initial wave energy as in runs A and B.

No electron acceleration is observed in run D. On the contrary run E was initialized with a single sinusoidal wave, with the smallest wavelength. Fig 10 shows that it is an efficient electron accelerator. These two simulations indicate that the short Alfvén wavelengths contribute more efficiently to the electron acceleration. The issue is that the Alfvén wave with short wavelengths is more dispersive (the non dispersive waves correspond to the case of the long MHD wavelengths), and do not favour the coherence of the Alfvén wave packet over long distances. Therefore, for efficient acceleration by Alfvén wave packets, there is a balance to find, in terms of the size of the packets, between the dispersive effects, and the efficiency of the electron acceleration.

Influence of the ion to electron mass ratio

Within the MHD framework (ω ≪ ω ci ) Alfvén waves constitute a single mode. It is actually the degeneration of two different modes, with right-hand and left-hand circular polarizations, whose dispersion relations become quite distinct at high frequencies.

Tsiklauri was that the acceleration process is less efficient with a higher ion to electron mass ratio, and they concluded that the acceleration is not controlled by the ion polarization drift amplitude (as with inertial Alfvén waves) but by the ion cyclotron frequency. This is contrary to the interpretation developed in the present paper, and we explore below this contradiction.

The inertial effect is caused by the ion polarization drift, that is proportional to the ion mass m i , and by the electron inertia, that is proportional to the inverse of the electron mass m -1 e . Therefore, the inertial effects are a growing function of the ion to electron (m i /m e = 100). In our simulations with right-hand polarized waves, it appears that when we increase the mass ratio, the growth of the parallel electric field is more significant.

This property is compatible with the predominance of the inertial effect on the Alfvén wave propagation. Figure 12 shows the energy flux of the downward electrons. Here again, we observe a larger electron acceleration for the largest mass ratio (400). The comparison of the acceleration for mass ratios of 100 and 200 is not so simple. Anyway, we do not observe that the acceleration is more efficient with a decreasing mass ratio.

Our results are therefore different from those of Tsiklauri [2007]. Nevertheless, there is no contradiction in the physics, since we do not explore the same branch of wave propagation.

Let us note that both our work and those by Tsiklauri [2007] were conducted with reduced ion to electron mass ratio. We can expect that for the real mass ratio, the efficiency of the acceleration by right-hand polarized waves would be better than in run B. On the contrary, the acceleration provided through a cyclotron wave would be weaker, loosing importance among the auroral acceleration processes. Of course, modifying the mass ratio also changes the ion cyclotron period, as well as the Alfvén velocity to electron thermal velocity ratio, that are important in both kinds of acceleration processes. Therefore, this simple test alone does not assess whether the process underlying acceleration in the A more thorough study of the influence of the mass ratio will be proposed in a forthcoming study.

Discussion and conclusion

In this paper, we have confirmed that an Alfvénic pulse propagating downward along a plasma cavity is able to generate transverse small scales, and thanks to inertial effect, a parallel electric field. This parallel electric field causes a broad spectrum of accelerated electrons, mainly in the downward direction. The electrons are accelerated within the wave packet, but they can propagate faster, and be present at lower altitudes. The acceleration process increases the downward flux of kinetic energy of the electrons, and this should favour the triggering of optical auroras in the ionosphere. Short Alfvénic pulses (small wavelengths) accelerate the electrons more efficiently than long ones. A short study of the efficiency of the acceleration process as a function of the ion to electron mass confirms that the accelerating parallel electric field is compatible with inertial effects (inertial Alfvén wave triggered acceleration).

Watt et al.

[2010] performed a series of numerical simulations to evaluate the efficiency of the shear Alfvén waves in providing the electron energy flux that powers the aurora.

They used a self consistent drift-kinetic code which allows to simulate a significant portion of an auroral magnetic field line, taking into account the variations of the magnetic field strength and plasma density with altitude. By comparison, our simulation box is much smaller than that, and we do not take into account the convergence of the lines of force. box, which corresponds to those of an inertial Alfvén wave packet. Further effects depend on the propagation of this wave packet along the converging magnetic field lines, and the dynamics of the electrons. The authors have computed the electron energy flux, and they conclude that the input of energy given by the wave may be sufficient to excite detectable auroras (in three cases, over the four studied), for a detection threshold set at 1 mW.m -2 .

The amplitude of the wave packet injected at the upper boundary is quite strong, with an associated parallel electric potential drop of several hundreds of Volts. However their

Figure 7 shows that about 80% of the Poynting flux associated with the wave disappear in the very first cells of their simulation box. The initialization procedure may be in cause: it seems that the the Alfvén wave electromagnetic field is included in the initial conditions, but not the associated perturbations of the electron and ion velocities. Therefore, we can expect that with a more consistent initialization of the Alfvén wave packet, the efficiency of the electron acceleration would be even better than claimed by the authors. The correct Alfvén polarization is established quickly, but at the expense of the initial amplitude of the wave electromagnetic field. Therefore, the true amplitude of the Alfvén wave in their simulation does not corresponds to the one set at time 0, and used as a reference.

Let us notice that the initial electron distribution in their simulations already contains an energetic tail, as it is modeled with a kappa distribution function. With the consideration of the plasma cavity interaction with the wave, the high energy tail of the electron distribution appears self-consistently. In both studies it is shown that an Alfvén wave packet can provide a significant increase of the electron energy flux, favoring the excitation of observable polar auroras.

The acceleration of electrons through the interaction of Alfvén waves and plasma cavities is not bound to happen only in the Earth auroral zone although the Earth environment is the only region where, up to now, Alfvénic processes have been observed in situ. Indeed, since the pioneering work of [START_REF] Heyvaerts | Coronal heating by phase-mixed shear Alfven waves[END_REF], many theoretical works have been devoted to Alfvénic acceleration on transverse density gradients in the Sun corona.

Recent works have been specifically devoted to non-MHD effets [Tsiklauri et al., 2005;[START_REF] Mcclements | Inertial Alfvén Wave Acceleration of Solar Flare Electrons[END_REF][START_REF] Bian | Parallel electric field amplification by phase mixing of Alfven waves[END_REF]. There are also observational evidences that Alfvén waves can accelerate electrons in the vicinity of Io and Jupiter [START_REF] Hess | Jovian S burst generation by Alfven waves[END_REF] sometimes combined with other acceleration structures [START_REF] Hess | Effect of electric potential structures on Jovian S-burst morphology[END_REF].

However in the Jovian case, it is not obvious that the plasma cavities are the sources of the inertial effects [START_REF] Mottez | Explanation of dominant oblique radio emission at Jupiter and comparison to the terrestrial case[END_REF] and it is possible that they result from wave filamentation at the border of the Io plasma torus [START_REF] Hess | Power transmission and particle acceleration along the Io flux tube[END_REF].

Up to now, we are not aware of any measurement of the length scale of the auroral cavities along the magnetic field. In our simulations, it has been considered as infinite.

We plan to simulate the effect of Alfvénic pulses propagating over finite size auroral cavities, and evaluate the efficiency of the acceleration mechanism according to the ratio between the characteristic length scales of the cavity and of the Alfvénic pulse.

  Figure3shows E y (x, y(x, t), t) for run B. The waves are the same as in run A, but they

  min = x-10∆x and x max = x+10∆x. In order to focus on the acceleration process, that happens essentially in and around the density depletion, computed fluxes for particles with y ′ in the middle half of the simulation box, (between 0.25L y and 0.75L y ). There is no selection of the particles based on their velocities. All the particles (thermal and suprathermal) are included in the computation of the flux.We have computed four distinct components of the electron kinetic energy flux. They are F eu , computed with the upgoing particles only, F ed with downgoing particles only, F ecu with the upgoing particles only with y ′ in the middle half of the simulation box (in and around the density channel), and F ecd with downgoing particles only in the middle half of the box.

Figure 8

 8 Figure 8 is a plot of the time evolution of these four components, for x = 102.4, in the

D

  waves) is weak compared to the case of the electrons, and the upgoing and the downgoing ion accelerations have approximately the same efficiency.Tsiklauri et al. [2005] made a

  density at the same (moving) location as the wave packet. In order to localize the (x, y)) and parallel (E x (x, y)) electric fields at the same times. The electrons can reach a velocity 0.8 = 8v te , and the wave packet velocity is v P ∼ 0.5, therefore many of the accelerated electrons propagate, at least temporarily, in front of the wave packet. Actually, a comparison with Fig. 4 shows that the wave packet, as seen with E y is very extended. On the contrary, the parallel electric field E x , seen on Fig. 5, has a very localized structure with peaks of high amplitude spread over a distance ∼ 50 only (the smaller peaks of parallel electric field are discussed in section 4.2). At time 409, they are localized around x = 200, and this is precisely where we can see, on Fig. 7, the largest number of accelerated particles. The population of accelerated particles extends on the two sides of this area, with the most energetic particles ahead of it. In that situation, a space probe passing below an acceleration region would measure the energetic downgoing electrons slightly before the parallel electric field of the Alfvén wave at the origin of their acceleration. This is possible only

  [2007] performed simulations of acceleration by ion cyclotron waves, in the continuation of the left-hand polarized Alfvén wave branch. Our simulations are based on the propagation of right-hand polarized Alfvén waves, in the low frequency part of the whistler branch. We made this choice because the right-hand polarized waves are less dispersive, more like the MHD waves usually involved in the literature. In particular with ÉNOT: ALFV ÉN PULSE IN CAVITY the right-hand polarized mode, nothing special happens at frequencies close to the ion cyclotron frequency. The observations of Solitary Kinetic Alfén Waves encourage us to think that, in the auroral zone, Alfvén wave packets are only weakly dispersive. For left-hand polarized waves with a small wavelength -i.e. ion cyclotron waves-the case is different.The dispersive nature of the propagation constitutes an important aspect of the problem, because the ion cyclotron frequency corresponds to a resonance. Tsiklauri[2007] has shown a coupling between the parallel electric field and the perpendicular electric field driven by the Alfvén wave, that, for the left-hand polarized mode is strongly influenced by ion gyrofrequency effects, with important consequences on the acceleration process. This is what he observed in his simulation, based on left-hand polarized waves. As the ion cyclotron frequency is inversely proportional to the ion to electron mass ratio, in order to characterize this resonnant effect, Tsiklauri[2007] performed several simulations in which the ion to electron mass ratio was varied. A result

  mass ratio. Therefore, a test on the dependency of the acceleration on the mass ratio helps to discriminate the control by the gyrofrequency (growing with m e /m i ) from the effects (growing with m i /m e ). Following Tsiklauri[2007], we have conducted a series of numerical simulations with the same parameters as in run B, except for the mass ratio. Fig 11 shows the time evolution of the maximum value of the parallel electric field E x , for simulations run B (m i /m e = 400), run F (m i /m e = 200), run C

  simulations of the present paper is indeed caused by inertial effect. It simply gives a clue, ÉNOT: ALFV ÉN PULSE IN CAVITY in response to Tsiklauri [2007]'s numerical experiments.

  Watt et al. [2010] took into account the ion polarization drift which induces the Alfvén wave inertial effects, and produces, in their formulation, a scalar electric potential drop.Electrons are accelerated by a parallel electric force injected at the top of the simulation

  Comparing the work ofWatt et al. [2010] and ours it appears that our simulations offer an insight into the origin of the parallel electric field and the microscopic aspects of the ÉNOT: ALFV ÉN PULSE IN CAVITY Alfvénic acceleration that cause auroras, whileWatt et al. [2010] focus mainly on their macroscopic effects.

  We are also interested to see if waves of stronger amplitude can contribute to regenerate or to

Figure 1 .

 1 Figure 1. Run A: temporal stack plot of the transverse electric field E y (x, t) along a

Figure 2 .

 2 Figure 2. Run B: map of the density at times (a) t =0, (b) t =179, and (c) t =409.

Figure 3 .Figure 4 .

 34 Figure 3. Run B: transverse electric field E y (x, t) along the curve y = Y (x, t) shown

Figure 5 .

 5 Figure 5. Run B: Map of the parallel electric field E x (x, y) at the same times as in Fig.

Figure 6 .Figure 7 .

 67 Figure 6. Run B: parallel electric field E x (x, t) along the curve y = Y (x, t) shown on

Figure 8 .

 8 Figure 8. Run B: energy fluxes of electrons (see Eq. 3) for x = 102.4. The flux F eu

Figure 9 .

 9 Figure 9. Run B: energy fluxes of ions. Same computations as in Fig. 8, but for the

Figure 10 .

 10 Figure10. Run E: electron distribution function f e (x, v x ), as in Fig.7.

Table 1 .

 1 Simulation parameters. δB max /B 0 is set in order to have the same initial magnetic energy in all simulations.

	X -31

MOTTEZ & G ÉNOT: ALFV ÉN PULSE IN CAVITY

Acknowledgments. The authors acknowledge support from the Institut National des Sciences de l'Univers (INSU/CNRS) of France under the program PNST. The numerical simulations were performed at the computing center (DIO) of the Paris-Meudon observatory.

Journal of Geophysical Research (Space Physics) 106, 12863-12874. Sydorenko, D. and Rankin, R. and Kabin, K., Oct. 2008 Volwerk, M., Louarn, P., Chust, T., Roux, A., de Feraudy, H., Holback, B., Jun. 1996.


Ludovic Dutreve
email: ludovic.dutreve@liris.cnrs.fr

Alexandre Meyer
email: alexandre.meyer@liris.cnrs.fr

Sada Bouakaz
email: saida.bouakaz@liris.cnrs.fr

Computer Animation and Virtual Worlds Easy Acquisition and Real-Time Animation of Facial Wrinkles

Keywords: facial animation, wrinkle acquisition, fine-scale animation, RBF interpolation Facial animation, wrinkle acquisition, fine-scale animation, RBF interpolation

Facial animation details like wrinkles or bulges are very useful for the analysis and the interpretation of facial emotions and expressions. However, outfitting a virtual face with expression details for real-time applications is a difficult task. In this paper, we propose a mono-camera acquisition technique of facial animation details and 2 Page 2 of 28 http://mc.manuscriptcentral.com/cavw -For Peer Review Computer Animation and Virtual Worlds a technique which add a wrinkle map layer (fine-scale animation) to a skinning layer (large-scale animation) for real-time rendering of a virtual 3D face. The acquisition is based on ratio image computed from two pictures of a same face, with and without expression. The real-time dynamic wrinkles technique is based on a small set of reference poses. These two methods offer an easy and low-cost way to capture facial animation details and use it for real-time facial animation.

Introduction

Facial animation details like wrinkles or bulges are very useful for the analysis and the interpretation of facial emotions and expressions [START_REF] Courgeon | Impact of expressive wrinkles on perception of a virtual character's facial expressions of emotions[END_REF]. However, outfitting a virtual face with expression details for real-time applications is a difficult task. The main reason is the difficulty to capture or synthesize wrinkle phenomena which human used to see every day, especially in a real-time context. Since many works have been proposed for large-scale animations and deformations [START_REF] Ersotelos | Building highly realistic facial modeling and animation: a survey[END_REF][START_REF] Gain | A survey of spatial deformation from a user-centered perspective[END_REF], only a few of them deal with real-time animation of small-scale details on face.

We denote small-scale details in an animation context, i.e. wrinkles and bulges appearing while muscles contractions, instead of the micro-structures of the skin independent to facial expressions. An interesting approach was proposed by Oat [START_REF] Oat | Animated wrinkle maps[END_REF], it consists on blending wrinkle maps to render animated details on a human face. This technique provides good results at an interactive time by using common bump-mapping technique for the rendering.

However, it requires manual tuning of blending coefficients, leading to manual efforts for each new animation, and requires wrinkle maps which may be difficult to create.

Capturing these facial animation details is still a challenging issue. On one hand, current capture methods focus on capturing the whole 3D geometry of a face with heavy active system based on lasers, structured lights or gradient-based illumination, with multiple cameras [START_REF] Weise | Fast 3d scanning with automatic motion compensation[END_REF][START_REF] Alexander | The digital emily project: photoreal facial modeling and animation[END_REF][START_REF] Wilson | Temporal upsampling of performance geometry using photometric alignment[END_REF]. Even if recent researches propose passive systems [START_REF] Beeler | High-quality single-shot capture of facial geometry[END_REF][START_REF] Bradley | High resolution passive facial performance capture[END_REF], stereo aspect is still the minimum, with often six or more high resolution cameras. On the other hand, many appli-4 Page 4 of 28 http://mc.manuscriptcentral.com/cavw -For Peer Review Computer Animation and Virtual Worlds cations want to apply captured motions to other faces than the one filmed. Thus, they do not need 3D high resolution meshes of the actor from which desired features are hard to extract and to transfer.

In this paper, we propose two techniques to both capture details and apply them on a real-time facial animation. To avoid complex motion capture system, we propose an easy and low-cost technique to capture facial details. On several frontal photos of different facial expressions, we capture corresponding normal maps with a few manual effort of landmarking on pictures. By considering only frontal single-camera views, we address the problem of normal reconstruction using the shading information in an inverse problem, where classical facial capture approaches use stereo. Our dynamic wrinkle technique consists on adding a wrinkle map layer to a skinning layer. The skinning layer is a common bone-based largescale animation technique on which we add the fine-scale details thanks to a local-area based wrinkle map layer. In order to apply wrinkle map, we use a small set of reference poses.

A reference pose is a pair of a large-scale deformation (a skeleton pose) and its associated small-scale details (a wrinkle map). During the animation, the current skeleton pose is compared with the reference poses and wrinkle maps coefficients are automatically computed with the use a non-linear function. Notice that comparison is done at a bone level resulting in local blends and independence between areas of the face.

5

Since our approach is based on single camera for the normal maps acquisition, we focus on single image shape reconstruction, and specially on inverse shading approaches which have to be be classified in the class of Shape from Shading (SfS) problems. The SfS problem has been widely studied in the computer vision area [START_REF] Durou | Numerical methods for shape-from-shading: A new survey with benchmarks[END_REF]. It is known as difficult because of its ill-posedness [START_REF] Prados | Shape from shading: a well-posed problem ?[END_REF]. Consequently, few approaches have been tried on real photographs [START_REF] Vogel | Making shape from shading work for real-world images[END_REF].

By considering the normal-map reconstruction instead of the 3D shape we relax a bit the difficulty of the problem [START_REF] Wu | Interactive normal reconstruction from a single image[END_REF]. Moreover, normal maps are well suited to real-time rendering and to our wrinkles-synthesis-from-example method. Wu et al. propose an interesting interactive normal reconstruction from a single image [START_REF] Wu | Interactive normal reconstruction from a single image[END_REF]. By aiming at reconstructing only local features such as wrinkles, we can limit the phase of user interaction which is mainly dedicated to give global information. Considering photo of face gives an a priori on the shape form which helps to solve the ill-posedness aspect and allows to reconstruct the large scale aspect of a face [START_REF] Smith | Facial shape-from-shading and recognition using principal geodesic analysis and robust statistics[END_REF]. To our knowledge, state of the art on SfS (understanding mono-camera) does not show convincing results on skin wrinkles acquisition. Thus, they are often procedurally synthesized [START_REF] Bickel | Multi-scale capture of facial geometry and motion[END_REF].

While some techniques were proposed to generate wrinkles on arbitrary surface [START_REF] Bando | A simple method for modeling wrinkles on human skin[END_REF][START_REF] Wang | Example-based wrinkle synthesis for clothing animation[END_REF], few approaches focused on real-time facial applications. Larboulette et al. proposed a technique to simulate dynamic wrinkles [START_REF] Larboulette | Real-time dynamic wrinkles[END_REF]. The user defines wrinkling areas by drawing a perpendicular segment of wrinkles (wrinkling area), following by the choice of a 2D discrete 6 Page 6 of 28 http://mc.manuscriptcentral.com/cavw -For Peer Review Computer Animation and Virtual Worlds control curve (wrinkle template). The control curve conserves its length while the mesh deformation, generating amplitude variations. Wrinkles are obtained by mesh subdivision and displacement along the wrinkle curve. Many methods need a high resolution mesh or a on-the-fly mesh subdivision scheme to generate or animate wrinkles [START_REF] Lo | Capturing facial details by space-time shape-from-shading[END_REF][START_REF] Na | Hierarchical retargetting of fine facial motions[END_REF]. Due to realtime constraints and the resource sharing, these techniques may be difficult to use in an efficient way. GPU computing allows to render efficiently fine details by using bump maps for interactive rendering. Oat presented in [START_REF] Oat | Animated wrinkle maps[END_REF] a GPU technique to easily blend wrinkle maps applied to a mesh. Maps are subdivided in regions, for each region, coefficients allow to blend the wrinkle maps. This technique requires few computational and storage costs, three normal maps are used for a neutral, a stretched and a compressed expressions. Furthermore it is easily implemented and added to an existing animation framework. The main drawback of this method is that it requires manual tuning of the wrinkle maps coefficients for each region. Similarly, De Melo et al. synchronize a normal map with a pseudo-muscular model to avoid manual coefficient tuning [START_REF] De Melo | Expression of emotions using wrinkles, blushing, sweating and tears[END_REF]. Normal maps are obtained by a manual image manipulation step. They directly use luminance information on a single photography and manually fit details to the character's texture. Our details extraction approach is more formal and requires less manual intervention. Our animated details method aims at providing a real-time dynamic wrinkling system for skinned face which generates automatically wrinkle maps coefficients and does not require a mask map.

7

In this section, we detail how we capture the skin details normal maps of several facial expressions. Since we do not need to capture the entire geometry of the face, our approach uses the illumination variation between two pictures of a same face with different expressions. Let I be the image of the neutral pose, and I ′ the image of the expression pose where the details we want to extract are. We seek to compute the normal map which will be applied during the real-time facial animation to produce details.

Normal map generation

The first step of our method consists of computing the ratio image to extract light variation [START_REF] Liu | Expressive expression mapping with ratio images[END_REF]. In order to work on pixel at correct positions, we start to deform I ′ into I ′ d to fit I. This allows us to obtain a mapping between each pixel, and so, to know the difference of illumination at each point of the face surface. The deformed image I ′ d is obtained by using a Radial Basis Functions (RBF) scattered data interpolation [START_REF] Siddiqui | A locally constrained radial basis function for registration and warping of images[END_REF] between two sets of landmarks manually picked on the two images I and I ′ . The ratio image R is computed by

dividing I to I ′ d : R = I ′ d I (1)
Then, we perform a blurring filter to R to reduce acquisition noise. To avoid artefacts due to eventual registration errors, we give the possibility to define regions of interest with 8 a simple painting tool. Furthermore, these masks help to avoid defining a large number of landmarks to obtain a perfect mapping, even in the areas that the user does not need to extract features. Figure 1 shows an example of pictures of different expressions with the associated ratio image.

Let C d the intrinsic color of the skin, C l the color of the light, N = (N x , N y , N z ) t (resp.

N ′ = (N ′ x , N ′ y , N ′ z ) t
) the normal at the surface of the skin of the neutral image (resp. the interpolated expression image 1 ), and L the light-direction vector. If we assume that skin is a diffuse surface, under the Lambertian model,

I = C d ×C l × N • L and I ′ d = C d ×C l × N ′ • L.
We will retrieve N ′ in two step: the normal tilt N z and the normal orientation ( N x , N y ). The Formula 1 gives:

R = I ′ d I = C d × C l × N ′ • L C d × C l × N • L = N ′ • L N • L (2)
By using the flash light of the camera, we can approximate the light direction to L = (0, 0, -1), given R = N ′ z /N z . We make the approximation that regions of interest are locally plane at the neutral pose. By default, we assume that these planes are parallel to the camera plane, resulting in N z = 1, but the user can define different values for N z if considered necessary, to refine the result.

Equation 2 shows that the ratio image intensity is dependant to the orientation of the surface. More R is dark, more N ′ z is small and more the wrinkles are deep. Notice that this use of the ratio image avoids the painful problem of finding the intrinsic color from a photo. 1 Notice that we work into the tangent coordinate system of the neutral image. We should now compute the normal orientation ( N x , N y ) for a given pixel of the wrinkles at each pixel. For this, N ′ z is not enough. We assume that the gradient of the image ratio is a good approximation of the normal orientation, which is mainly similar to the assumption done in [START_REF] Woodham | Photometric method for determining surface orientation from multiple images[END_REF][START_REF] Robles-Kelly | Shape-from-shading using the heat equation[END_REF]. Let G x and G y be the 2D normalized gradient of R for the two axes x and y. We want to find N ′ x and N ′ y such that N ′2

x + N ′2 y + N ′2 z = 1, and N ′ x = αG x and N ′ y = αG y with α a normalization coefficient. We obtain the formulae:

N ′ x = G x × 1 -N ′2 z N ′ y = G y × 1 -N ′2 z (3)
We then convert N ′ into a normal map norm. Figure 1 illustrates a normal map obtained with a ratio image.

3D character adaptation

The last step we perform for details acquisition is the interpolation between the capture space to character texture space. It aims at ranging the computed normal map to the character. We use the same method as for the interpolation between the neutral and the expression images, two sets of landmarks are manually defined on both normal map and character's texture map, and a RBF based interpolation deformed the first to the second (see Figure 2).

Reference Poses

A reference pose is a pair of a large-scale deformation (a skeleton pose) and its associated small-scale details (stored in the form of a normal map we call wrinkle map). To create a reference pose, we deform the facial skeleton to obtain the expression on which we want the wrinkles appear. This expression should be strongly pronounced to provide better results.

Then, we associate to the pose a wrinkle map extracted in the acquisition step. Since influences work with bones or group of bones, it is possible to define a pose where wrinkles appear in various areas of the face. Having details in different areas will not cause that all of these details appear at the same time at an arbitrary frame. For example, details around the mouth would appear independently with forehead details, even if they are described in a same reference pose. 11

Pose Evaluation

While runtime, the first step is to compute the influence of each reference pose on each bone.

This consists in finding how the bone position at an arbitrary frame looks like its position at the reference poses 2 . Computing these influences at the bone level instead of a full pose level allows to determine regions of interest. This offers the possibility to apply different reference poses at a same time. Resulting in the need of less reference poses (actually, only 2 may be sufficient for face: a stretched and a compressed expression).

We define the influence α jf k of the pose Π k for the bone P j at the position P jf for the arbitrary pose Π f by the following equation:

α jf k = min   1, max   0, AB • AC AC     (4)
where AB = P jf -P j0 , AC = P jk -P j0 , • denotes the dot product and ... denotes the Euclidean distance. α jf k is the size of the segment between the orthogonal projection of P jf on the segment [ P jk , P j0 ] normalized by the size of the segment [ P jk , P j0 ].

Bones Masks

Once we know the reference poses influences for each bone, we could use them for the per-pixel lighting. The main idea is to use skinning weights to compute the influence of 2 Notice that we deal with bones positions in the head root bone coordinates system, so we could assume that head rigid transformation will not cause problems while pose evaluation.

12 Page 12 of 28 http://mc.manuscriptcentral.com/cavw -For Peer Review

Computer Animation and Virtual Worlds reference poses for each vertex, and by the interpolation done during the rasterization phase, for each fragment (Fig. 3). Since wrinkles and expressive details are greatly related to the face deformations, we can deduce that these details are related with bones transformations too. So we associate bones influence with reference poses influences. A mesh vertex v i with a skinning weight of w ij with the bone P i is influenced by the reference pose Π k at the arbitrary pose Π f such as:

β if k = n i=1 λ ij × θ(α jf k ) (5)
with

• λ ij = 1 if w ij > 0 and λ ij = 0 else;
• θ(x) the modulation function described in the next section.

Modulation function

Wrinkles do not appear on face linearly according to bones displacements. Indeed, two bones moving from the neutral position to the compressed one will create wrinkles more intensively at this end of the movement than at the beginning. For this purpose, we introduce during the synthesis phase a non-linear function θ such as θ : [0, 1] → [0, 1] which rules this phenomena by modifying the distance used during the normal map blending. To propose a function based on real wrinkles apparition, we run our normal map acquisition algorithm on each frame of a video sequence where a face start from a neutral pose to an expressive 13 one. It provides for each pixel and for a given bones configuration, a vector of normal map displacement. The length of this vector indicates the intensity of the wrinkles. Since the wrinkle synthesis is done on the GPU, it would have been inefficient to store for each pixel a set of discrete values for each bones configurations. We preferred fit an analytic curve:

θ (x) = cos (π (1 + x)) + 1 2 (6)
Figure 4 shows the segment curve we capture and the curve we propose. The cosine function seems to be an interesting base if we apply to it a translation and a normalization to check our constraints. We try to depict a smooth and progressive acceleration of winkles amplitude.

We give priority to highest percentage of deformation, when wrinkles are deepest, whereas the distance between the two curves, which seems high at the begin of the deformation, produces only a weak visual impact. Notice that user can define a new function and/or capture his own data to better fit his needs if necessary.

Details Blending

The final step of our method is to apply wrinkle maps to our mesh by using coefficients θ (β if k ) computed at the previous step. Since normal maps only contain deformations associated with the expression, we should apply a finest blending to avoid loss of static details located on the neutral normal map (i.e. details such as pores, scares and other static fine details). For example, with a simple averaging, a fragment influenced by 100 percents of a wrinkle map will be drawn without using the neutral map, resulting in the fact that details 14 of the neutral map will not appear. Let n f the final normal, n i the ρ normals provided by the wrinkle maps, their coefficients w i and n the normal provided by the default normal map, we first compute weighted normals η i of the wrinkle maps and then compute the final normal:

η i = (0, 0, 1) t + w i × n i -(0, 0, 1) t (7)
w f =           n x + o i=1 ρη i x n y + ρ i=1 η i y n z × ρ i=1 η i z           (8)
The addition of the two first coordinates makes a simple averaging between the direction of the two normals, given the desired direction. The z components are multiplied, this leads to increase details obtained from the two normals. More z is small, and more the surface will be bumpy, multiplication allows to add the neutral surface variation to the wrinkled surface.

Results and Discussion

Our acquisition technique provides good results with a few effort and material. Despite the lighting model approximation, which does not take into account the specular reflection of the skin, we obtain very interesting visual results. Our method provides an alternative to 15 others facial motion captures ystems if the goal is to obtain only animation fine-scale details and not the full geometry of the face. Figure 5 shows an example of use of a captured normal map by our monocular acquisition technique.

Our 3D test models Proog 3 is about 10000 triangles. The faces is rigged with 21 bones.

Animation runs at more than 100 fps on a common personnal computer with an Athlon 64 3800+ CPU, 3Go RAM and a NVidia Geforce GTX 275 GPU. The Figure 6 shows the neutral and the two reference poses used for demonstration. The Figure 7 shows the 6 basis expressions of Proog with and without our dynamic wrinkles and the Figure 5 shows the possibility to apply dynamic details on different area of the face even if details are given on a same wrinkle map.

Large scale skinning-based deformation are not modified. The addition of our method to an existing implementation is easy. No additional rendering pass is required. Our choice to use skinning weights as bones masks offers many advantages. They allow us to relate large scale and small scale deformations, and so, we do not need additional mask textures.

They ensure that vertices influenced by reference poses are vertices which are displaced accordingly with bones too.

Although our technique easily provides wrinkles on facial animation, several steps are important to obtain good results. First, a good rigging is primary since we directly use skinning weights as bones masks, and so, it defines how each vertex will be influenced by the Problems occur when capturing faces with beard where mouth details can not be obtain because of the incoherence in the ratio computation. Our approach is more suited to soft skin subjects.

Conclusion

We 28

9 Page 9

 99 of 28 http://mc.manuscriptcentral.com/cavw -For Peer Review Computer Animation and Virtual Worlds

10 Page

 10 Our goal is to add a wrinkle map layer (fine-scale animation layer) to a bone-based facial animation (large-scale animation). We use the wrinkles data captured at the acquisition step to generate on the fly visual dynamic fine details such as wrinkles or bulges on the face surface. We explain in this section how we use a small set of reference poses in real-time and at arbitrary poses.

3 c

 3 Copyright 2006, Blender Foundation / Netherlands Media Art Institute / www.elephantsdream.org 16 Page 16 of 28 http://mc.manuscriptcentral.com/cavw -For Peer Review Computer Animation and Virtual Worlds different reference poses. Second, the reference poses should consist on an orthogonal set of skeleton poses as much as possible, to avoid an over fitting. Notice that blending reference poses in a same area is possible but it becomes a problem if similar bones displacements lead to different fine-details. Finally, detail maps quality greatly influences the visual results.

Figure 1 : 22 PageFigure 2 : 23 Page

 122223 Figure 1: Ratio image computation using two pictures of a same face with two expressions

Figure 3 : 24 PageFigure 4 : 25 PageFigure 5 :Figure 6 :

 32442556 Figure 3: The two left images show the skinning influences of two right eyebrow bones.

Figure 7 :

 7 Figure 7: The first row shows the 6 basis expressions (joy, sad, disgust, surprise, fear,

Page 14 of 28 http

://mc.manuscriptcentral.com/cavw -

For Peer Review Computer Animation and Virtual Worlds

For Peer Review Computer Animation and Virtual Worlds of

  Our dynamic animation wrinkles runs in real-time, the use of per-pixel lighting allows us to dispense with high-resolution meshes or costly subdivision techniques. the skin in our acquisition step. We will also try to improve our image-based details

	Page 17 of 28
	extraction to a video-based one, allowing to apply winkles in real-time from a captured
	face. Regarding the 3D aspect, last advances in graphics hardware provide local tesselation
	ability, it should be easy to adapt our bump-map method to a local mesh refinement with
	very similar input textures.
	Furthermore, it is based on widely-used techniques such as skinning and bump mapping. Its
	implementation does not present technical difficulties and does not modify usual animation
	and rendering pipeline. In futur works, it may be interesting to consider specular reflection
	17

have presented two techniques. The first aims at capturing facial animation details from two pictures of a same face with and without expression. It requires few manual work and does not need specific material. A simple camera with its flash light was used. The second is a technique which use reference poses to generate non-linearly in real-time wrinkles and fine-details appearing with an arbitrary skinned face animation. In addition to providing interesting visual results, the requirements that we considered necessary and/or important have been met. http://mc.manuscriptcentral.com/cavw -