L Couëdel

M Mikikian

L Boufendi

A A Samarian

Residual dust charges in discharge afterglow

Keywords: numbers: 52, 27, Lw

An on-ground measurement of dust particle residual charges in the afterglow of a dusty plasma was performed in a rf discharge. An upward thermophoretic force was used to balance the gravitational force. It was found that positively-charged, negatively-charged and neutral dust particles coexisted for more than one minute after the discharge was switched off. The mean residual charge for 200 nm radius particles was measured. The dust particle mean charge is about -5e at pressure of 1.2 mbar and about -3e at pressure of 0.4 mbar.

I. INTRODUCTION

Dusty or complex plasmas are partially ionized gas composed of neutral species, ions, electrons and charged dust particles. In laboratory experiments, these particles can be either injected or grown directly in the plasma. Injected dust particles are usually micron-size particles. Due to their mass, they are confined near the bottom electrode where the electric force counterbalances gravity. Microgravity condition is necessary to study dust clouds of micrometer size particles filling the whole plasma chamber [1]. In laboratory, dense clouds of submicron particles light enough to completely fill the gap between the electrodes can be obtained using reactive gases such as silane [2,3] or using a target sputtered with ions from plasma [4][5][6][7] Dust particle charge is a key parameter in complex plasma. It determines the interaction between a dust particle and electrons, ions, its neighboring dust particles, and electric field. The determination of the dust particle charge is so one of the basic problems in any complex plasma experiments. The knowledge of dust charge will allow us to understand the basic properties of dusty plasma, particle dynamics in dust clouds, and methods to manipulate the particles. The particle charge controls the dust dynamics both in laboratory and technological plasma reactors, and also in space plasmas. Thus one of the main dusty plasma challenges is to understand the dust charging in a wide range of experimental conditions, which simulates industrial and space plasmas. For example, in most industrial processes in the microelectronics industry which uses silane as the reactive gas, dust contamination is a vital problem [START_REF] Bouchoule | Dusty Plasmas: Physics, Chemistry and Technological impacts in Plasma Processing[END_REF]. The dust dynamics and coagulation in space plasma are also phenomena governed by dust particle charges [START_REF] Belov | [END_REF]10]. There are many publications reporting on the investigation of dust charging in discharge plasma, see [START_REF] Vladimirov | Physics and Applications of Complex Plasmas[END_REF][START_REF] Tomme | [END_REF][13][14][15][16][17][18] and references therein. However there are only few papers devoted to dust charging, or decharging to be more specific, in the discharge afterglow. In [19] a diffusion of * Electronic address: lenaic.couedel@univ-orleans.fr fine dust particles (8-50nm) in the afterglow of a dusty plasma has been studied. From the afterglow decay of dust number density the charged fraction of particle have been measured. It was found that some particles carried a very small residual charge and some were neutral. A model has been proposed to explain charges on particles in the late afterglow of a dusty plasma. An observation of decay, or "decharging", of a rf plasma with dust particles after switching off discharge power was reported in [20]. Experiments were performed in the PKE-Nefedov reactor [1] under microgravity conditions on board International Space Station. An existence of negative residual charges has been shown. A simple theoretical model that describes data obtained was proposed. The residual charges were attributed to the presence of an excitation electric field that was used for dust charge measurements. Limitation of the Space Station experiment routine could not prove this hypotheses. Theoretical predictions on an effect of rf plasma parameters on residual charges has neither been verified. So any new experimental evidence of residual dust charges in afterglow plasma will lead to a better understanding in the decharging of complex plasma. Indeed the value and nature of residual charges after plasma extinction is of great importance. The dust particle charge in afterglow plasma could induce problem in future single-electron devices where a residual charge attached on deposited nanocrystals would be the origin of dysfunction. It could make easier industrial plasma processing reactor decontamination thanks to the use of specially designed electric field. The residual charges on dust particles in fusion reactors (such as ITER) also can make the cleaning process much easier. In this paper, we report the first on-ground experiment on the residual charge measurement of dust particles after decay of a dusty plasma. The experiment was performed in the PKE-Nefedov reactor where the dust particles were physically grown in the plasma. Temperature gradient was introduced in the chamber to create an upward thermophoretic force to balance gravity. The residual charges were determined from analysis of dust oscillations, which were excited by applying sinusoidal bias to the bottom electrode. It was found coexistence of positively and negatively charged dust particles as well as non-charged ones for more than one minute after the discharge was switched off. ticles have been measured for two different pressures. It was revealed that dust particles kept the residual charges only when the discharge was abruptly switched off. In the case when the discharge power is decreased slowly until the plasma disappeared, there was no residual charge on dust particles. It was also shown that the presence of the low frequency excitation electric field did not play any role in keeping the charge on dust particles in the afterglow plasma.

The article is organized as following way. Section II is devoted to general description of experimental set-up. In section III we analyze the forces acting on dust particles, paying most attention to the thermophoretic force.

In section IV, we discuss a procedure of residual charge measurement and present experimental data obtained, in section V we analyze discharging of dust particle in afterglow plasma, and section VI is conclusion.

II. EXPERIMENTAL SETUP

The work presented here is performed in the PKE-Nefedov (Plasma Kristall Experiment) chamber designed for microgravity experiments [1]. It is a rf discharge operating in push-pull excitation mode. It consists of 4 cm diameter parallel electrodes separated by 3 cm. The injected power varies in the range 0 -4 W . Dust particles are grown in an argon plasma (0.2 -2 mbar) from a sputtered polymer layer deposited on the electrodes and coming from previously injected dust particles (3.4 µm, melamine formaldehyde). A detailed description of this experiment and previous results are presented in Ref [1,4,5]. For the study concerning residual charges, the top electrode was cooled (Fig. 1). An upward thermophoretic force was applied to dust particles in order to counterbalance gravity [21] when the plasma is off. To study particle charges, a sinusoidal voltage produced by a function generator with amplitude ±30 V and frequency of 1 Hz was applied to the bottom electrode. Induced low frequency sinusoidal electric field E(r, t) generated dust oscillations if they kept a residual electric charge. A thin laser sheet perpendicular to the electrodes illuminates dust particles and the scattered light is recorded at 90 • with standard charge coupled device (CCD) cameras with 25 images per second. Video signals were transferred to a computer via a frame-grabber card with 8 bit gray scale and 560 × 700 pixel resolution. In order to avoid edge effect, a field of view over 8.53 × 5.50 mm 2 restrained to the center of the chamber was used for residual charge measurement. By superimposition of video frames particles trajectories have been obtained. The coordinates of the particles were measured in each third frame. The amplitude of the oscillations was figured out from the measured particle positions. Absolute values for the oscillation amplitude were obtained by scaling the picture pixels to the known size of the field of view.

III. FORCES ACTING ON DUST PARTICLES

Dust particles in laboratory plasmas are subjected to various forces that confine them in the plasma or drag them out (to the walls or pump outlet) [START_REF] Bouchoule | Dusty Plasmas: Physics, Chemistry and Technological impacts in Plasma Processing[END_REF]: the confining electrostatic force, the gravity, the thermophoretic force, the ion drag force and the neutral drag force. After the discharge is switched off, the forces still acting on dust particles are the gravity F g , the thermophoretic force F T , the neutral drag force F dn and, if dust particles keep residual charges, the electric force F E due to the electric field induced by the sinusoidal voltage applied to the lower electrode. In this experiment the thermophoretic force balanced the gravity force, and electrostatic force was the reason for the dust oscillations while the neutral drag force damp the last ones. Expression for the gravity force is:

F g = 4 3 πr 3 d • ρg (1)
where g is the gravity acceleration, r d is the dust particle radius, and ρ its mass density. The neutral drag force was taken as [START_REF] Shukla | Introduction to Dusty Plasma Physics[END_REF]:

F dn = - 8 3 √ 2πr 2 d m n n n v T n (1 + α ac π 8 )(v d -v n ) (2)
where m n is the neutral atom mass, n n the neutral atom density, v T n = 8k B T n /πm n the thermal speed with T n = T the neutral gas temperature, k B the Boltzmann constant, α ac the accomodation coefficient, v d dust particle speed and v n a mean speed of neutral atoms. In our experiment v n can be taken as 0 because there is no gas flow.

The electric force can be written as

F E (r, t) = Q d E(r, t) (3)
where Q d is the dust particle electric charge and E(r, t) is the electric field between the electrodes after the plasma is switched off. In the late afterglow the plasma density is very small and E(r, t) can be approximated by the vacuum field above a charged disk. Taking into account size of the camera field of view, size of the electrode and the fact that camera was directed to the region near center of the chamber, we are interested in the electric field vertical component E(r, t) ≈ E z (z, t) -→ e z ≡ E(z, t) -→ e z only. The calculated value of E(z, t) at the central axis is presented in Fig. 2 The expression for the thermophoretic force must be choosen carefully. Indeed, it depends strongly on Knudsen number Kn = l/r d [START_REF] Zheng | [END_REF] where l is the mean free path of buffer gas species. In our experiment, we work at an operating pressure around 1 mbar. In a previous paper [4], the size of grown dust particles was reported between 200 nm and 800 nm. It gives, using results from Varney [24] for an atom-atom cross section, Knudsen number 250 < Kn < 1000. Consequently we operate in free molecular regime where a dust particule is similar to a very large molecule. Many theories have been developped [START_REF] Zheng | [END_REF][25][26][27][28] and used [21,[START_REF] Zheng | [END_REF]29] for thermophoresis in the free molecular regime. The most commonly used equation is the Waldmann equation [28] which has been verified experimentally [30,31]:

F T = - 32 15 r 2 d k tr v T n ∇T ( 4
)
where ∇T is the temperature gradient in the gas, and k tr the translational part of the thermal conductivity given for a monoatomic gas by [32]:

k tr = 15k B 4m n µ ref T T ref ν (5)
where µ ref is the reference viscosity at reference temperature T ref = 273 K and the exponent ν results from a best fit of experimental viscosity near the reference temperature. For argon, µ ref = 2.117 • 10 -5 P a • s and ν = 0.81 [32].

Another important effect that must be considered in the estimation of thermophoretic force is an influence of the finite volume of gas. If pressure is low enough, gas mean free path can become comparable to a length scale of experimental apparatus and the gas can no longer be treated as a continuous medium. Under such condition, an additional Knudsen number must be added [START_REF] Zheng | [END_REF] Kn L = l/L where L is the length scale of the reactor. In this experiment, the length between electrodes is L = 3 cm giving Kn L ∼ 5 • 10 -3 which means the gas can be considered as a continuous medium. Temperature gradient between the electrodes was calculated by FEMLAB c code (steady state analysis of heat transfer through convection and conduction with heat flux, convective and temperature boundary condi- tions using Lagrange-Quadratic element). The temperatures on the electrodes were measured by thermocouple and used as boundary conditions for the problem. The countour-plot obtained is presented in Fig. 3. It shows that the vertical component of the temperature gradient is constant near reactor center. The value of temperature gradient is about 2K/cm for our experiment. There is also a small horizontal component of the temperature gradient. This is the reason for particle drift in horizontal direction. Such drift allowed us to resolve particle trajectories (see next section) and made particle charge measurement more convenient.

IV. EXPERIMENTAL PROCEDURE AND RESULTS

In this section, dust particle size estimation and residual charge measurement are presented. The residual charge measurements have been performed by the following routine. First the chamber was pumped down to lowest possible pressure (base pressure ∼ 2 • 10 -6 mbar) and the cooling system was turned on. After this, argon was injected up to the operating pressure, the discharge was started and particles were grown, forming familiar structures such as a void ( see for example [5]). Then, the discharge was switched off and bottom electrode was biased by sinusoidal voltage. In afterglow plasma, dynamics of dust particles is determined by temperature gradient and excitation electric field. Fig. 4 presents superimposition of images taken after discharge had been switched off. There are two different types of motion observable. Dust particles drift upwards, downwards and to the side due to existing temperature gradients and they oscillate due to electrostatic force. It is obvious that thermophoretic force acts on any dust particle in chamber, while the electrostatic force acts only if particles have charge in afterglow. Thus presence of oscillating particles (see Fig. 4) clearly indicates that dust particles do have residual charges after the discharge has been switched off. Dust particles oscillated in opposite phases as well as non-oscillating dust grains have been observed indicating that negatively charged, positively charged and non-charged dust particles coexist after plasma extinction.

It worth to mention that in order to observe dust oscillations the discharge must be switched off abruptly. It was shown that if the power was decreased slowly until the plasma disappears there are no residual charges (no oscillations of dust particles in the sinusoidal electric field were observed). Another interesting fact is that residual charge on dust particles has long relaxation time and does not depend on time when the excitation electric field was applied to. Dust oscillations were observed for more than one minute after plasma extinction and in both cases when the function generator was switched on during the disharge or few seconds after the discharge is turned off.

As it can be seen in Fig. 4, there are dust particles falling after the discharge is switched off. These particles are too big to be sustained by the thermophoretic force. Other particles are horizontally drift at constant height, this means for these particles gravity force is balanced by thermophoretic force. These particles have been used to measure residual charges. It is clear from Fig. 4, that use of a large field of view camera gives us nice overall picture of decaying dusty plasma but it is not suitable for residual charge measurement because edge effects can not be neglected. Thus the camera with the small field of view was used for the charge measurement (see Fig. 4). The superimposition of images from this camera is presented in Fig. 5. These images give us clear track of the dust oscillations so dust grain trajectories can be reconstructed (Fig. 6).

A. Size and mass of levitating dust grains in discharge afterglow

Particle size (mass) and residual charge measurement are strongly related in this experiment. Charge, size and mass of dust particles have so to be determined. Considering that dust particles levitating in reactor at a constant height after plasma extinction are the ones for which the gravity is exactly balanced by the thermophoretic force, the dust particle radius can be found using Eq.1 and Eq.4:

r d = - 8 5πρg k tr v th ∇T (6)
Dust grains are supposed to be spherical and mainly made of carbon (sputtering of a carbonaceous polymer material) [4] so that the mass m d can be deduced from Eq.6:

m d = 4 3 πr 3 d • ρ ( 7
)
where ρ is the mass density of graphite. In our experimental condition, the radius of levitating dust grains is estimated to r d ≃ 190 nm and their mass is m d ≃ 6.5 • 10 -17 kg.

∇T -190 K • m -1 -177 K • m -1 r d 194 nm 180 nm m d 6.9 • 10 -17 kg 5.4 • 10 -17 kg γ 1.36 • 10 -13 kg • s -1 0.39 • 10 -13 kg • s -1 Q dresmax +2e +2e Q dresmin -12e -4e Q dresmean -5e -3e
TABLE I: Measured maximum, minimum and mean dust particle residual charges for two operating pressures.

B. Dust charge measurement

From the measurement of oscillation amplitude, the residual charge on a dust particle can be obtained. As the gravity is compensated by the thermophoretic force, the equation of motion for one dust particle, neglecting its interactions with other dust particles, can be reduced to:

m d z = F E (z, t) + F nd ( ż) (8)
Taking E(t) = E 0 (z mean ) cos(ωt) (the amplitude of the electric field E 0 is the one at the mean dust levitation height z mean ) and using Eq.3, Eq.2 and Eq.8, the dust particle oscillation amplitude b can be obtained [START_REF] Landau | Physique Théorique[END_REF]:

b(ω, Q d , E 0 ) = Q d E 0 (z mean ) m d ω ω 2 + 4γ 2 /m 2 d (9)
where γ = (4/3) √ 2πr 2 d m n n n v T n (1 + α ac (π/8)) is the damping coefficient and ω = 2πf where f is the frequency imposed by the function generator. Eq.9 can easily be inversed to derive the dust residual charge:

Q dres = m d b(ω, Q d , E 0 (z mean ))ω ω 2 + 4γ 2 /m 2 d E 0 (z mean ) (10)
The sign of the dust particle charge is deduced from the phase of the dust particle oscillation with respect to the excitation electric field. Oscillation amplitudes up to 1.1 mm have been measured (depending on the operating pressure) and charges from -12e to +2e are deduced where e is the elementary charge. It has been found that at high pressure dust particles keep a higher mean residual charge (Tab.I) but uncertainties are of ∼ 2e for each measurement.

V. DISCUSSION

The charging (discharging) process of dust particle in a plasma is governed by the contributions of all currents entering (or leaving) the dust surface, involving plasma electron and ion currents, photoemission and thermionic emission currents, etc:

dQ d dt = I a + I l (11)
where I a and I l are currents absorbed and emitted by the particle (with appropriate sign). In most cases for discharge plasmas we can ignore the emission current and kinetics of the particle charge can be expressed as:

dQ d dt = J i -J e = -πer 2 d [n e v Te e -ϕ -n i v Ti (1+ Te ϕ)] (12)
where J e and J i are the fluxes of electrons and ions onto the particle, v Ti(e) = 8k B T i(e) /πm i(e) the thermal velocity of ions (electrons), Te = T e /T i is the electron to ion temperature ratio, n i(e) is the density of ions (electrons) and ϕ = -eQ d /4πǫ 0 k B r d T e is the dimensionless surface potential of a dust particle where ǫ 0 is vacuum dielectric permitivity. According to Eq.12, charge on dust particle depends on the electron-ion masses, temperatures and density ratios m e /m i , n e /n i , T e /T i . Thus to analyze the decharging of dust particle in afterglow plasma one need to consider the kinetics of plasma decay.

The plasma diffusion loss and electron temperature relaxation determine kinetics of the discharge plasma decay [START_REF] Raizer | Gas Disharge Physics[END_REF]. In presence of dust particles, plasma loss is due to diffusion onto the walls completed by surface recombination on dust particles. Equations for plasma density and electron temperature relaxation are [START_REF] Raizer | Gas Disharge Physics[END_REF][START_REF] Tsytovich | [END_REF]:

dñ dt = - ñ τ L ( 13
)
d Te dt = - Te -1 τ T (14)
where ñ = n i,e /n 0 with n 0 the initial plasma density, τ L the time scale of the plasma loss, and τ T is the time scale for electron temperature relaxation. The expressions for the time scales are [START_REF] Raizer | Gas Disharge Physics[END_REF][START_REF] Tsytovich | [END_REF]:

τ -1 L = τ -1 D + τ -1 A (15) τ -1 D ≃ l in v Ti 3Λ 2 (1 + Te ) ≡ 1 2 (1 + Te ) 1 τ ∞ D (16) τ -1 A ≃ πr 2 d n d v Ti (1 + ϕ Te ) ≡ 1 + ϕ Te 1 + ϕ 1 τ ∞ A (17) τ -1 T = π 2 m e m i v Ti l en Te ≡ √ T e τ ∞ T (18)
where τ D is the ambipolar diffusion time scale onto the walls, τ A the particle absorption time scale, l i(e)n the mean free path of ion (electron)-neutral collision, Λ is the characteristic diffusion length (Λ ∼ 1 cm in this experiment). The ∞ exponent stands for the limit at very long time. For charging time scale lower than plasma decay or temperature relaxation time scales the charge on dust particle is in equilibrium, i.e. ion and electron fluxes balance each other, ϕ ≃ ϕ eq and using Eq.12, ϕ eq is given by:

n e n i Te e -ϕeq = m e m i (1 + Te ϕ eq ) (19)
In this case, expressions for charge fluctuation and charge fluctuation time scale τ Q are: where λ i0 = ǫ 0 k B T i /n 0 e 2 is the initial ion Debye length.

dQ d dt ≃ - Q d -Q deq τ Q (20) τ -1 Q ≃ v Ti r d 4λ 2 i0 (1 + ϕ eq )ñ ≡ ñ τ 0 Q (21)
It should be noted that the time scale for dust charge fluctuations strongly depends on plasma density and can vary from microsecond for initial stages of plasma decay up to seconds in case of almost extinguished plasma. Taking into account Eq.13 and Eq.21, the time dependence of τ Q can be expressed as [20]:

τ -1 Q = 1 τ 0 Q exp(-t/τ L ) (22)
To understand the dusty plasma decharging dynamics we have to compare different time scales. In Tab.II, time scales for this experiment are presented. It can be seen that the initial charge fluctuation time scale is the shortest. The temperature relaxation time scale is shorter or becomes comparable (for 0.4mbar) to the plasma density decay time scale, and plasma losses mainly determined by diffusion. The latest means that for our experimental conditions dust particles did not affect plasma decaying at initial stage.

Q d n / i n e n 0 T e / T i n / i n e T e / T i n 0 Q d t p t c t d
The fig. 7 presents the qualitative dependence of the main plasma and dust particle parameters during the afterglow. Four stages of the dust plasma decay can be labeled. As we can see the first stage of the plasma decay (t < τ T ) is characterized by the electron temperature T e drop down to the room temperature, while the plasma density (especially in case τ 0 T < τ 0 L ) is slightly decreased. As the charging time scale is almost independent on Te (Eq.21), the charge is still determined by its equilibrium value (Eq.12). During the temperature relaxations stage the particle charge should decrease to the value [20]:

Q rT = 1 Te0 ϕ eq (1)
ϕ eq ( Te0 )

Q 0 ≃ Q 0 62 ≃ -15e (23)
where Q 0 is the initial dust charge in the plasma and Q rT the value of dust residual charge at the end of first decay stage. The dust charge in the plasma Q 0 was estimated as Q 0 = -950e solving numerically the equation 19, with given parameters T i = 300 K and T e ≃ 3 eV , for argon plasma with n i ∼ = n e . At the next stage of decay, electron temperature is stabilized while the plasma density is still decreasing (see fig 7). So τ Q continues increasing according to Eq.12 and Eq.21. When τ Q becomes comparable to τ L , the particle charge can not be considered as in equilibrium and to determine particle charge we should use the Eq.12. The time scale when the particle charge starts sufficiently deviating from the equilibrium can be estimated as (Eq.15 and Eq.22):

t d ∼ -τ ∞ L ln 8 3 • λ i0 Λ 2 • l in r d ∼ 6τ ∞ L (24)
However, according to Eq.12, as long as plasma is neutral (n e = n i ) the charge on dust particle does not change. The plasma will keep quasineutrality until decaying rates for the electrons and ions are the same. It will be true in the case of ambipolar diffusion. When the nature of diffusion changes, electrons and ions start to diffuse independently that will lead to changing n e /n i ratio and consequently the dust charge change.

The nature of the plasma diffusion changed when the particle volume charge can not be ignored or when the plasma screening length becomes comparable to the chamber size. In first case the ion diffusion will be influenced by the negatively charged dust particles, while the electrons will be free to go. In second case, large density differences appear over distances less than the screening length and electrons and ions diffuse independantly. Lets us estimate the characteristic times for both cases. The influence of the overall particle charge is determined by the value of Havnes parameter P e = -N Q d /en e . Based on the model discussed, the qualitative evolution of the Havnes parameter P e in dusty plasma afterglow can be plotted as shown in Fig. 8. The initial value of P e is small (∼ 0.06 with an estimated dust density N ∼ 2 • 10 5 cm -3 and n 0 ≃ n e0 ∼ 5 • 10 9 cm -3 ) and there is no influence of dust. At the first stage of decaying (temperature relaxation stage) P e decreases due to dramatic decrease of dust charge while the plasma density decrease of a factor 1.1. At τ T , P e reaches its minimum value about 1 × 10 -3 . After this P e starts increasing. During this stage dust particle charge changes slowly while plasma number density decays fast (see Fig. 7). The time at which P e becomes ∼ 1 can be estimated as (Eqs. [13][14][15][16][17][18]:

t p ∼ τ ∞ L ln T e0 T n -en 0 Q d N ∼ 8τ ∞ L ( 25
)
The screening length becomes comparable to the chamber size, i.e. λ i (ñ c ) ∼ Λ when the density drops down to ñc = λ 2 i0 /Λ 2 . This occurs at [20]:

t c ∼ τ ∞ L ln ñ-1 c ∼ 10τ ∞ L (26)
At time min[t p , t c ], electrons start running away faster than ions and the ratio n i /n e grows. For our experimental conditions t p < t c , thus the neutrality violation due to presence of dust particles happens before Debye length exceeds the chamber size. So the third stage of dusty plasma decay starts at t p . During this stage the charge on dust particles is changed due to changing of n e /n i ratio. At this stage t d < t p , thus the kinetic Eq.12 should be used for estimations of the charge variation. The upper limit of the charge change can be estimated ignoring the electron current and considering the time interval between t p and t c ,

dQ dt < J i    tc tp (27)
< πer 2 d n i (t p )v Ti 1 - e 4πǫ 0 k B r d T i Q d    tc tp (28)
Solving Eq.28, the charge should evolve following:

Q d = Q dT - 1 α exp -Kα∆t + 1 α (29)
where α = e/4πǫ 0 k B r d T i ∼ 0.28/e, K = πer 2 d n i (t p )v Ti ∼ 190e, and thus (Kα) -1 ∼ 20 ms ∼ ∆t = (t c -t p ). Therefore, the charge during the third stage decreases to -3e. At the forth stage of plasma decay, t > t c , the plasma density decreases such that any further changes of dust charge become negligible and charge remains constant for a while. Thus the final residual charge for our condition is expected to be about Q dres ∼ -3e which is well correlated with the charges measured in the experiment (Tab.I). The observation of neutral and positively charged particles in our experiments can be explained by the particle charge distribution. The charge distribution in dust particle ensemble is due to charge fluctuation on every individual dust particle. It has been shown [36][37][38][39][40] that rms of stochastic charge fluctuations varies as δ • √ < Q d >, where δ is a parameter depending on plasma conditions and close to 0.5. Thus for mean charge around -3e the charge distribution should lay mainly between -1e and -5e. However in our case we can expect broader distribution taking into account two reasons. First for small particles with a charge of about few electrons the image force (polarizability) correction starts playing role, see [37]. And for the electron temperature less then 0.1eV the correction factor can be up to 3. That is significantly broadening distribution so dust particles can experience fluctuations to neutral and positive charges, see fig. 9. Another reason for the broadened distribution could be the earlier freezing of the charge distribution so it remains the same while the mean charge decreases. Indeed as t > t d the charge is not in equilibrium with surrounding plasma so the existing charge distribution should remain unchanged or changed slightly during the later stages. Unfortunately, in this experiment the ratio between positively and negatively charged particles has not been measured due to experimental limitations concerned with the fact that dust grains do not stay in the laser sheet for a sufficiently long time due to their out of plane thermal motion. Thus to elaborate the true reason of the presence of residual charges of different polarities one need additional measurement of charge distributions.

The considered decharging model predicts a higher residual charge for lower pressure. For low pressure the temperature relaxation time scale exceeds plasma density decay time (see characteristics times for 0.4 mbar in Table II). This leads to quick increase of charging time and makes charge freezing at early stage of decay keeping charge at relatively high value. In fact bigger charges were measured for the higher pressure. A possible expla-nation for this is an electron re-heating in late afterglow from metastable atoms [41]. The presence of a significant amount of argon metastables in PKE discharge have been established by plasma modeling [42] and by experimental measurements [START_REF] Mikikian | Proc. 30 th EPS Conf. on Contr. Fusion and Plasma Phys[END_REF]. Due to re-heating the electrons can have temperature few times more than the room temperature which will lead to the increase of residual charge on dust particles. As re-heating is more effective for highpressures thus bigger charges are expected.

VI. CONCLUSION

Residual dust particle charges have been measured in late afterglow of a dusty plasma. Dust particles have been grown directly in the plasma by sputtering of a polymer material previously deposited on electrodes. Dust particles with a radius of few hundreds of nanometers are levitating after the power of the discharge is switched off. The gravity was balanced by an upward thermophoretic force. Residual charges were determined from analysis of dust oscillations, which were excited by applying sinusoidal bias to the bottom electrode. Positive, negative and non-charged dust particles have been detected. Mean residual charge for 200 nm radius particles was measured. The particle charge is about -5e at pressure of 1.2 mbar and about -3e at pressure of 0.4 mbar. A model for the dusty plasma decay was exploited to explain the experimental data. According to this model the dust plasma decay occurs in four stages: temperature relaxation stage, density decay stage, dust charge volume stage, and frozen stage (ice age IV). The main decrease of the dust charge happens during the first stage due to cooling of the electron gas. The final residual charge established during the third stage when the density of ions exceeds the density of electrons and the plasma density is still high enough to affect the charge. Measured values of the dust residual charges are in a good agreement with values predicted by the model. However the residual charge dependence on discharge conditions and detection of positively charged particles show that a more detailed model taking into account various phenomena (electron re-heating, electron release [START_REF] Berndt | [END_REF], afterglow chemistry) in decaying plasma have to be developed for better understanding of dusty plasma afterglow.

  FIG. 2: (Color online) Vertical component of the electric field in the camera field of view

FIG. 3 :

 3 FIG. 3: Temperature profile and gradient

FIG. 4 :

 4 FIG.4: Superimposition of video frames taken with a large field of view CCD camera. Arrows 1 and 2 represent respectively the vertical and the horizontal components of the temperature gradient. Edge effects as well as falling dust particles can be seen.

FIG. 5 :

 5 FIG. 5: Superimposition of video frames 10 seconds after plasma extinction. Dust particle oscillations can clearly be seen. The temperature gradient has a slight horizontal component. Therefore, oscillations are in the 2D laser plane.

FIG. 6 :

 6 FIG.6: (Color online) Oscillation of 6 dust grains 10 seconds after plasma extinction. Non oscillating dust grain and opposite phase oscillations are observed.

FIG. 7 :

 7 FIG. 7: (Color online) Qualitative time evolution of dust charge, plasma density and electron temperature during the afterglow Four stages of the dust plasma decay can be identified: I -temperature relaxation stage up to tT II -plasma density decay stage up tp, III-dust charge volume stage tc, IV -frozen stage.

FIG. 8 :

 8 FIG.8: Qualitative time evolution of the Havnes parameter after the power is switched off.

FIG. 9 :

 9 FIG. 9: Charge distribution with and without image force correction. Dotted line represents distribution with ∆Q = 0.5 √ Q. Solid line represents distribution with ∆Q = 0.87 √ Q [37].

  The residual charges for 200 nm radius par-

			COOLING SYSTEM
			ELECTRODE
	laser sheet	GLASS	plasma & particles	GLASS
			ELECTRODE	field of view
	FIG. 1: (Color online) Schematic of the experimental appa-
	ratus

TABLE II :

 II Values of the different time scales for two operating pressures.

		P = 0.4 mbar P = 1.2 mbar
	τ 0 D	30 µs	90 µs
	τ ∞ D	1.5 ms	5 ms
	τ 0 A	0.4 ms	0.4 ms
	τ ∞ A	30 ms	30 ms
	τ 0 L	30 µs	90 µs
	τ ∞ L	1.4 ms	4.3 ms
	τ 0 T	100 µs	34 µs
	τ ∞ T	1 ms	340 µs
	τ 0 Q	4 µs	4 µs
	tc	20 ms	60 ms
	τQ(tc)	1.4 s	1.4 s

Acknowledgments

The authors would like to thank S. Dozias and B. Dumax for electronic support, J. Mathias for optical support and Y. Tessier for experimental support. The PKE-Nefedov chamber has been made available by the Max-Planck-Institute for Extraterrestrial Physics, Germany, under the funding of DLR/BMBF under grants No.50WM9852. This work was supported by CNES under contract 793/2000/CNES/8344. This paper is partially supported by the French-Australian integrated research program (FAST) of the French foreign affairs ministry and the International Science Linkages established under the Australian Government's innovation statement, Backing Australia's Ability'.


J S Brest

V Keryvin

P Longère

Y Yokoyama

Insight into plasticity mechanisms in metallic glasses by means of a Brazilian test and numerical simulation

Keywords: Metallic glasses, Diametral compression, Digital Image Correlation, Finite Element Modelling, Plasticity, Shear bands, Constitutive models

Multiaxial and heterogeneous mechanical experiments, by means of the diametral compression test (Brazilian test), were carried out on Zr-based Bulk Metallic Glass. Attention was notably focused on obtaining the displacement field in the area of interest (large strains) employing a Digital Image Correlation device, in addition to the measurement of the usual global load-displacement curve. In order to reproduce the behaviour of the BMG at stake, several constitutive equations were considered with growing complexity: von Mises, Drucker-Prager (pressure dependence), free volume based model as well as a viscoplastic Coulomb-Mohr type model. The two latter were implemented as user-material in a Finite Element computation code. The results of these investigations including experiments and computational simulations are discussed.

Introduction

Bulk metallic glasses (BMG) are a quite recent class of materials that exhibit exceptional mechanical properties including strength, resilience, hardness or fracture toughness as well as large yield strains [1,2,3]. The mechanisms of deformation, damage and failure of these materials must be fully identified for designing compositions of metallic glasses even more resistant and ductile. The numerical prediction of the mechanical resistance for these new materials requires an advanced modelling of their behaviour. A better understanding of the plastic deformation (in terms of localised thin shear bands), the consideration of temperature, viscosity and pressure effects are necessary parameters.

Obtained via a rapid solidification from the liquid state, BMG belong to the class of amorphous alloys characterized by the absence of structural long range ordering. Their macroscopic mechanisms of plasticity mainly proceed from the propagation and interaction of very thin shear bands. From the microscopic viewpoint, BMG plasticity does not result from the accumulation and motion of dislocations, as it is the case for crystalline alloys, but may be depicted, depending on the authors, by the theories of shear transformation zones (STZ) [4] or free volume (FV) [5]. Starting from the aforementioned models (STZ and FV), several authors have tried to describe the macroscopic response of BMG to various loading cases, see e.g. Schuh et al. [1] and Gao [6] (FV model), respectively. Following another way, Anand and Su (AS model) [7] developed a Coulomb-Mohr type constitutive theory to reproduce the elastic-viscoplastic behaviour of BMG.

The present work aims at investigating in a first step the mechanical behaviour on Zrbased BMG by means of a multiaxial and heterogeneous test. BMG exhibiting an apparent brittle behaviour in tension as well as in compression at room temperature, the information, namely plasticity, provided by the usual tension or compression tests on these materials are de facto far from being suitable. To get some insight on plasticity mechanisms, the indentation test is commonly carried out [8,9,10,11], so are bending [12] or constrained tests [13]. There exist other tests, though being less known for this type of material, which may also be of particular interest. Among them is diametral compression tests [14], also called Brazilian test, which is extensively used on brittle materials like concrete,. It allows one to generate a tensile loading in the sample section transversally to the compression axis; it also generates heterogeneous mechanical fields in the sample, some parts of the sample remaining almost free of loading. The present approach differs from previous indentation studies since it allows to record extensive plastic fields during the test in addition to the mean response (the load-displacement curve) and it creates multiaxial stress states different from those underneath an indenter. For the indentation test, they are highly compressive and localised underneath the indenter tip, while for the Brazilian test, they are compressive close to the loading conditions, tensile close to the free boundaries and close to pure shear at the centre.

In a second step, computational simulations employing the engineering Finite Element computation code Abaqus are conducted, involving configurations with different BMG constitutive models. The models which do not belong to the standard model library of the computation code were implemented as user materials. The numerical results are commented and compared to experimental ones.

The material and experimental procedures are displayed in Sect. 2. with the corresponding results. The constitutive models and the numerical procedures are presented in Sect. 3. The given results are then compared to experiments in Sect. 4.

Materials and experiments

Materials and samples

The Zr 55 Cu 30 Al 10 Ni 5 (at. %) BMG is studied in this paper. By a very careful tilt-casting process and with highly purified (50 ppm of oxygen) iodine refining processed Zr, specimens with less than 300 ppm of oxygen and no micrometric crystalline defects are prepared (see [15] for details). The material is supplied in the form of 60 mm long and 8 mm diameter rods. The samples used are cut and polished to obtain 6 or 15 mm long cylinders.

Experiments were carried out at room temperature and atmospheric pressure employing a hydraulic tension-compression machine (INSTRON 8803) equipped with a 500 kN capacity load cell with a 8 N resolution. A 0.8 µm resolution mechanical extensometer was used to measure the global displacement between the compression anvils, and two cameras were installed to monitor the displacement (with a 0.25 mm resolution) of a cloud of dots sprayed on a side of the sample (see Figure 1), and further derived the related strain field using a Digital Image Correlation (DIC) software [16] (Aramis). The tests were carried out at a displacement rate of 0.05 mm.min -1 . For questions of reproducibility, the same test was made on four different samples times from the same batch. The load-displacement curves were the same within less than 5%. An equivalent shear strain used in this paper is defined as :

ǫ eq = 2 3 e ∼ : e ∼ (1)
Where e ∼ is the deviatoric part of the strain tensor ǫ ∼ and ":" denotes the doubley contracted product.

Experimental results

The load-displacement (P-u) curve is presented in Fig. 2. The carried out test is a sequence of two successive loadings: the first one on a 15 mm sample, and the second one on the same sample cut to 6 mm long, after unloading of the former, to be able to reach fracture without damaging the compression anvils (100 kN limit load). Therefore P stands for a load per unit length. On Fig. 2 one can clearly distinguish two distinctive steps during the deformation process, a stiff part followed by a more compliant one. The numerical results presented in section 4 provide an explanation to this phenomenon. A higher magnification of the last part of the plot displays the presence of regularly spaced serrations, already reported elsewhere for different tests [13,17]. They are assumed to be the consequences of the sequential formation and emergence of shear bands inside and at the surface of the sample.

DIC allows one to capture the deformation pattern of the sample (see Fig. 3). During the first part of the loading (below ∼ 11 kN/mm) it concentrates in a column underneath the load with a maximum equivalent strain of ∼ 0.2 close to the anvil (indentation-like feature) and at the centre of the sample (Brazilian-like feature). At the end of the second loading, before fracture, most of the deformation fields are localised in a cross with a maximum intensity reached in its centre. The equivalent strain in this cross is more than 0.3, while in the centre the deformation is higher than 0.6. These strains are average values since the shear bands thickness is 10 to 100 nm [18] while the DIC's resolution is 0.25 mm. Much larger strains are expected inside the shear bands.

Constitutive modelling and numerical simulations

Principles of the different approaches

Metallic glasses are often compared to polycrystalline metallic alloys for their chemical composition and atomic bonds, with polymers for their behaviour at high temperature [19] and with ceramics for their brittleness [20]. Some models exist aiming at describing each of these materials behaviour, yet for amorphous metals their number is much smaller because there remain many unknowns in the understanding of the mechanisms governing their deformation. One can quote two specific theories proposed in the late seventies and which are references in metallic glasses behaviour: Spaepen's theory of free volume [5] and Argon's theory of shear transformation zones [4]. Depending on the underlying mechanisms assumed to be at the origin of the BMG deformation and on the final application, different state variables have to be accounted for in the constitutive modelling. For this reason, our modelling approach consisted in starting from very simple models and progressively increase their complexity. The constitutive models used in the following are accordingly: von Mises (VM) plasticity model involving the deviatoric stress; Drucker-Prager (DP) model involving the pressure in addition to the latter (some indentation studies found this model to be quite suitable [21,10,11]); the FV theory described by Gao [6] and AS model [7]. There exist several other models that have not been considered here [22,23,24].

Constitutive equations

The main equations describing each model are presented here. Let consider beforehand the decomposition of the total strain rate ǫ ∼ into a reversible part, namely ǫ ∼ e , and an inelastic part, namely ǫ ∼ p . In the small strains assumption :

ǫ ∼ = ǫ ∼ e + ǫ ∼ p (2)
In the context of linear isotropic elasticity, the following relationship holds:

ǫ ∼ e = 1 + ν E σ ∼ - ν E tr σ ∼ δ ∼ (3)
where σ ∼ represents the second order stress tensor. Note that σ ∼ = s ∼ -p δ ∼ , with s ∼ being the deviatoric part of σ ∼ , p the pressure and δ ∼ the second order identity tensor. E and ν are the Young's modulus and Poisson's ratio, respectively. The absence of strain hardening [1] lead us to use DP and VM as elastic perfectly plastic models (associated flow). This effect may be described using von Mises yield criterion which is written as:

f M (σ eq ) = σ eq -σ 0 ≤ 0 ( 4
)
where σ eq = 3 2 s ∼ : s ∼ represent the equivalent shear stress, σ 0 being the tensile yield strength. The asymmetry observed on BMG behaviour between yield strengths in tension and compression before failure may be described using the pressure dependent DP model widely used in soils mechanics:

f DP (σ eq , p) = σ eq -p tan ϕ -σ 0 ≤ 0 ( 5
)
where ϕ is the internal friction angle.

The next two models are specifically dedicated to amorphous materials with the purpose of representing inhomogeneous deformation under the form of localised shear bands. These works are inspired from two theories born in the late seventies. They argued that the specific behaviour of amorphous metals founds its origins at the atomic scale. Spaepen [5] assumed that the deformation is intimately related to the evolution of the free volume appearing during the cooling process leading to the formation of the amorphous metal. The free volume is frequently defined as the difference in volume between a crystalline alloy and its amorphous equivalent. The theory proposed by Argon [4] states that the deformation occurs by local shearing of clusters of atoms. According to his theory, tens of atoms slip on each other leading to local intense shearing. This phenomenon was reproduced using computational simulations aiming at describing interactions at the atomistic scale [25]. The two models presented hereafter have no yield criteria, the plastic strain is always present, but negligible during the so-called elastic part of the deformation.

For the FV model [6], the inelastic strain rate is written as:

ε ∼ p = exp - 1 v f sinh σ eq σ re f s ∼ σ eq ( 6
)
where v f represents the free volume whose evolution law is given by:

vf = 1 α exp - 1 v f 3(1 -ν) E σ re f χ v f cosh σ eq σ re f -1 - 1 n D ( 7
)
where σ re f is the reference stress, α is a geometrical factor, χ is a constant and n D is the number of atomic jumps needed to annihilate a free volume.

The last model considered in the present work is proposed by Anand and Su [7]. It aims at reproducing the elastic-viscoplastic response of pressure-sensitive and plastically-dilatant isotropic materials. The slip systems (m (α) and g (α) are the normal to the α th slip plane and its slip direction, respectively) are assumed to be related to the principal directions of stress, giving six of them in an isotropic material. The slip directions also depend on the angle of internal friction. The inelastic strain rate takes the form:

ε ∼ p = 6 α=1 γ(α) G ∼ (α) + β M ∼ (α) G ∼ (α) = s y m g (α) ⊗ m (α) M ∼ (α) = m (α) ⊗ m (α) (8)
Where α represents one of the slip systems considered (⊗ is the dyadic tensor product). β is the dilatancy parameter that controls the plastic flow direction:

β = g 0 1 - η η c v p (9)
The initial value of β is g 0 thus the inelastic strain is a combination of plastic slip and dilatation; and tends to zero when the plastic volumetric strain η (Eq. 10) reaches its limit η c v (the maximum density reduction as a result of plastic deformation) leading to plastic slip without dilatation. g 0 is a parameter that allows to dissociate the tension (g 0t ) and compression (g 0c ) behaviour. The plastic volumetric strain rate ( η) reads (γ is the total shearing rate):

η = β γ ( 10
)
Let γ(α) be the shearing rate on the slip system α. It is expressed via a viscoplastic Norton-like law with a Mohr-Coulomb type resistance coefficient:

γ(α) = γ0 τ (α) c + µσ (α) 1/m > 0 γ = 6 α=1 γ(α) (11)
Where τ (α) = σ ∼ : G ∼ (α) and σ (α) = -σ ∼ : M ∼ (α) represent the resolved shear stress and the normal stress on the slip system α, respectively. The quantity γ0 is a reference shearing rate. m is the strain-rate sensitivity, taken small for a nearly rate independent response. The constant µ is the internal friction coefficient. The evolution law of the cohesion c is taken in the form:

c = c c v + b 1 - η η c v q c(0) = c c v + b ( 12
)
Where c(0) is the initial cohesion, composed of the final cohesion c c v with an overshoot b. The parameters q and p are positive material constants of order unity.

Numerical procedures

The engineering Finite Element computation code Abaqus was used to conduct the computational simulations. VM and DP models belong to the material library of the code. FV and AS models (see Sec. 3) were implemented as user-material in Abaqus. The calculation is first made in the corotational framework via equations presented in Sec. 3.2, then the software takes into account the finite strains. Taking advantage of the symmetries of the boundary value problem, the sample is modelled by a quarter of circle composed with 4-nodes finite elements with reduced integration and assuming plane strain conditions (CPE4R). The compression anvil is modelled using an analytical rigid part and is submitted to a displacement at constant velocity with a frictionless contact.

The different coefficients used for the simulations are given in Tab. 1-4. The elastic constants have been determined using ultra-sound techniques [26]. For DP model, the yield strengths were chosen corresponding to the material's behaviour [27]. For VM model, the mean value of the two was used. Those two are elastic perfectly plastic models since there is no strain hardening in the material. For FV and AS models, the coefficients are taken from literature [6,7] and correspond to a Zr-based BMG.

Identification of material parameters

With the aim of obtaining a better match between experimental and numerical results, an automated material parameters identification procedure based on the hybrid method (Levenberg-Marquardt, gradient and Newton-Raphson algorithms) [28,29] is used. In doing so, only a fraction of the parameters (see Tab. 1-4) has been allowed to vary. Namely for DP model: the compressive and tensile yield stresses, σ 0c and σ 0t respectively; for AS model: the cohesion and the internal friction, c and µ respectively.

Computational results and discussion

The computational simulation offers the possibility to observe the inner stress fields, unavailable experimentally. The origin of the inflection point (Fig. 2 at ∼ 7 kN.mm -1 ) can then be linked to the coalescence of the high stress fields up and down the sample (Fig. 4). The fields of plastic strain developed at the contact join, forming a constant vertical field. Besides, one can notice the von Mises stress decreasing from the centre to the sides which are almost free of loading. That is why the shear bands do not propagate in this area (as we can observe numerically on Fig. 5).

The equivalent strain map obtained with FV (see Fig. 5) or AS models localise the strain in fine numerical shear bands [30]. Note that the simulations were conducted with a very fine meshing (120.000 elements) using a variable free volume for FV model and a variable cohesion for AS model. The shear bands directions obtained numerically correspond to those observed experimentally, their number and intensity obviously depend on the mesh size. Corresponding to each serration, one can notice (not shown here) on the deformation pattern the creation of a new band with high deformation. The magnitude of the serrations is closely related to the mesh size. The finer the latter, the smoother the former. This the drawback of models with softening like AS and FV: local fields depend on the mesh size and therefore can not be compared one-to-one to experiments. Let us point out finally that the computational simulations were made in plane strain conditions while the experimental data gave a plane stress condition map. Nevertheless, it appears that the most intense numerical shear band is one branch of the pattern observed experimentally (Fig. 3). The resolution obtained with the DIC device does not allow us to observe thin shear bands like in the numerical results, but a semi-local map of the deformation process.

The superposition of the experimental and numerical overall load-displacement curves are plotted in Fig. 6. From a qualitative point of view, the curves obtained with the VM and DP models follow an evolution similar to the experiment. Before the inflection point VM and DP models have a slope very close to the experimental one, while the amorphous dedicated models overestimate it. After this point, the slopes of the remaining curve are the same at the beginning, except for FV model. The latter exhibits a plateau which does not exist experimentally. For higher deformations, the experimental load grows faster than the numerical ones.

Optimizing the plasticity coefficients for DP and AS models (they have proven to be much better than the other models studied) has given interesting results (Fig. 7). The first part of the curve (label 1) is well described by DP, and the slope is a little overestimated by AS. For the second part (label 2), the slopes obtained are the same, the reloading is well described by the two models. On the third part (label 3) one can notice that the slopes of the two models slightly underestimate the experiments. The results obtained with those models give a possible set of plasticity coefficients presented in Tab. 5. The simulations not taking into account the stiffness of the compression anvils, the coefficients obtained are probably underestimated, and so represent a first approach before the use of a more accurate numerical procedure. As for DP model the parameters differ from those obtained by indentation techniques [21].

In terms of local information, and this is a key point of this paper, it is highlighted in Fig. 5 that two models that predict load-displacement curves very similar (see Fig. 2 give very different deformation fields.

Concluding remarks

In this study, we have carried out diametral compression or Brazilian tests on a Zr-based BMG. This test allows us to get an insight on the plastic behaviour of the material, unlike conventional tests (tension, compression) because of the brittleness of the material. The present approach differs from previous constrained tests studies such as indentation since it allows to record in situ extensive plastic fields during the test, via Digital Image Correlation techniques (local information), in addition to the mean response (the load-displacement curve, global information). Besides, for the indentation test, the stress fields are highly compressive and localised underneath the indenter tip, while for the Brazilian test, they are compressive close to the loading conditions, tensile close to the free boundaries and close to pure at the centre. Intense plasticity was found to develop in the sample, first like in an indentation test, then localised in a cross-like pattern. Values of the equivalent strain, prior to failure, was as high as more than 0.8, in the centre of this cross.

Through the use of constitutive models with growing complexity, we have observed that all models may reproduce the global information (load-displacement curve) with identified material parameters. However, they predict different contrasted levels of the local information (equivalent strain).

The search for a relevant model able to describe the plasticity of metallic glasses under various mechanical loadings will therefore require to take into account this local information, as well as other constrained tests like indentation or scratching [31] and uniaxial tests (yield strengths). E (MPa) ν (-) σ 0 (MPa) 85000 0.37 1700

Figure 1 :

 1 Figure 1: Schematics of the experimental set-up and its instrumentation.

1 )

 1 per unit length, P (kN.mm -

Figure 2 :

 2 Figure 2: Diametral compression results -Plot of load per unit length vs displacement -Magnification of the serrations at the end of the curve.

Figure 3 :

 3 Figure 3: Map of equivalent strain (logarithmic measure log(1+ǫ eq )) obtained by DIC ; at the end of the first loading (∼ 11 kN.mm -1 , top) and at the end of the second loading up to fracture (bottom, strain reinitialised to zero).

Figure 4 :

 4 Figure 4: Mises stress fields obtained numerically below the inflection point of the curve in Fig. 2 with VM model.

Figure 5 : 1 )

 51 Figure 5: Map of equivalent strain obtained for a 2 mm displacement: DP model (left) and FV model (right). In this latter case, the whole sample was meshed. The material parameters used are presented in Tab. 1-4. Note that a similar colormap was used for comparing.

Figure 6 : 1 )Figure 7 :

 617 Figure 6: Comparison between experimental and numerical results using the material parameters presented in Tab. 1-4.

Table 1 :

 1 Materials parameters used for VM model.

	18

Table 2 :

 2 Materials parameters used for DP model.

	E (MPa) ν (-) σ re f (MPa) n D (-) α (-) χ (-)
	85000	0.37	100	3	0.15	1

Table 3 :

 3 Materials parameters for FV model.

	20

Table 4 :

 4 Materials parameters used for AS model.

	21

Table 5 :

 5 Materials parameters obtained after identification.

Acknowledgements

Jean-Sébastien Brest acknowledges the Brittany region (ARED CODYVER) and DGA for his PhD grant ; the European Brittany University (EPT COMDYNVER, avec le soutien du Fonds Européen de Développement Régional) and the Brittany region (DCICOLL) for their financial support as well as Pr P. Pilvin for fruitful discussion.

List of Figures


Mohamed Jellouli
email: mohamed.jellouli@fsm.rnu.tn

Michel Mehrenberger

Optimal decay rates for the stabilization of a string network

, and obtain optimal decay rates when the lengths are all equal. We also define a classical space semi-discretization and compare the results with the exact method introduced in [3,4].

Introduction

We recall that the dissipation condition ∂ x u(t, 0) = α∂ t u(t, 0) at the origin of a vibrating elastic string fixed at its end point, stabilizes this string. More precisely, if E(t) denotes the energy of the solution of the wave equation ∂ 2 t u(t, x) -∂ 2 x u(t, x) = 0 on (0, ℓ) subject to the initial condition u(0, .) = a and ∂ t u(0, .) = b, then E(t) ≤ Ce -γαt where γ α = 1 ℓ log 1+α 1-α . Moreover, we remark that if α = 1, then E(t) = 0, ∀t ≥ 2ℓ. Thus, the value 1 is the best choice of α that makes the system is an equilibrium state. The situation is completely different in the case when we consider a network of strings. The problem of stabilization of nondegenerate network of strings was studied by K. Ammari, M. Jellouli and M. Khenissi in [START_REF] Ammari | Stabilization of Generic Tree of String[END_REF], [START_REF] Ammari | Stabilization of Star-Shaped Networks of Strings[END_REF] and [START_REF] Ammari | Remark on Stabilization of Tree-Shaped Networks of Strings[END_REF]. When the strings are coupled at a common end in a star-shaped configuration, it is proven in [START_REF] Ammari | Stabilization of Star-Shaped Networks of Strings[END_REF] that the solutions are not exponential stable in the energy space, in the nondegenerate case. In the particular case of two strings, which is equivalent to the pointwise stabilisation of one string, such a problem has been studied in [START_REF] Ammari | Asymptotic behaviour of the solutions and optimal location of the actuator for the pointwise stabilization of a string[END_REF]; exponential stabilization is obtained if and only if the lengths satisfy ℓ1 ℓ1+ℓ2 = p q , with p and q odd numbers, and the best decay rate, when fixing the total length ℓ 1 + ℓ 2 is obtained when the lengths are equal (ℓ 1 = ℓ 2 = ℓ), with best decay rate γ = ln (3) ℓ . Such problem has been also considered in [START_REF] Nicaise | A remark on the stabilization of the 1-d wave equation[END_REF], with different boundary conditions; in the degenerate case, the energy limit was identified and it was proven that the solution decays exponentially to that limit. We consider here the case of a degenerate network of vibrating elastic strings when the pointwise feedback acts in the root of the tree (tree shaped network). Note that in the nondegenerate case, it is proven in [START_REF] Ammari | Stabilization of Generic Tree of String[END_REF] that the solutions are not exponential stable in the energy space. We calculate the limit energy E ∞ := lim t→+∞ E(t) and we show that the decrease from E(t) to E ∞ is exponential, giving the best decay rate, when the lengths are equal. Finally we give numerical results in Section 3, which confirm the theoretical results.

E ∞ and best decay rate

Let N ≥ 3. We consider the initial data ((

a j ) 1≤j≤N , (b j ) 1≤j≤N ) ∈ H := N j=1 H 2 (0, ℓ j ) × N j=1 H 1 (0, ℓ j ),
satisfying the compatibility conditions

a ′ 1 (0) = αb 1 (0), a ′ 1 (ℓ 1 ) = N j=2 a ′ j (0) and a j (ℓ j ) = 0, a 1 (ℓ 1 ) = a j (0), j = 2, . . . , N, (1)
and the system of partial differential equations (S) :

                         ∂ 2 t u j (t, x) -∂ 2 x u j (t, x) = 0, t > 0, x ∈ (0, ℓ j ), j = 1, . . . , N, u j (t, 0) = u 1 (t, ℓ 1 ) and u j (t, ℓ j ) = 0, t ≥ 0, j = 2, . . . , N, ∂ x u 1 (t, 0) = α∂ t u 1 (t, 0), t ≥ 0 (α > 0), ∂ x u 1 (t, ℓ 1 ) = N j=2 ∂ x u j (t, 0), t ≥ 0, u j (0, x) = a j (x), ∂ t u j (0, x) = b j (x), x ∈ [0, ℓ j ], j = 1, . . . , N.
We define the energy of the solution of (S) by

E(t) = 1 2 N j=1 ∂ t u j (t) 2 L 2 (0,ℓj ) + 1 2 N j=1 ∂ x u j (t) 2 L 2 (0,ℓj ) (2)
and we prove the following theorems in the case where ℓ j = ℓ, j = 1, . . . , N (we denote by the norm of L 2 (0, ℓ)): 1), then the energy limit is given by

Theorem 2.1 If ((a j ) 1≤j≤N , (b j ) 1≤j≤N ) ∈ H verify (
E ∞ = 1 2(N -1) N j=2 N k=j+1 (a ′ k -a ′ j ) 2 + (b k -b j ) 2 . (3)
2 Now, we denote by

α 0 = 2 √ N -1 N , ∆ = 4(α 2 N 2 -4N +4) (α+1) 2
and λ = -

2(N -2) α+1 + √ ∆ 2N
. 1 Theorem 2.2 There is a constant C = C N (α) > 0 such that for any initial data ((a j ) 1≤j≤N , (b j ) 1≤j≤N ) ∈ H satisfying (1) the following inequality holds

0 ≤ E(t) -E ∞ ≤ C a ′ 1 2 + b 1 2 + N j=2 a ′ j 2 + N j=2 b j 2 e -γt |∆| , if α = α 0 , (4)
0 ≤ E(t) -E ∞ ≤ C a ′ 1 2 + b 1 2 + N j=2 a ′ j 2 + N j=2 b j 2 t 2 e -γ0t , if α = α 0 , (5)
where

γ = 1 ℓ log 1 |λ| > 0 and γ 0 = 1 ℓ log N +2 √ N -1 N -2
. The best decay rate is achieved in the sense that for all α > 0, there exists an initial data such that ( 4) and ( 5) become an equivalence:

E(t) -E ∞ ≥ Ce -γt , or, for α = α 0 , E(t) -E ∞ ≥ Ct 2 e -γ0t
, with C > 0, and ∀α = α 0 , we have γ < γ 0 . Remark 1 The case N = 2 (with ℓ 1 , ℓ 2 not necessarily of same size) can be recasted into the case N = 1

(with ℓ = ℓ 1 + ℓ 2 ).
The proofs of these theorems are based on the operator of type τ developed in [START_REF] Ammari | Remark on Stabilization of Tree-Shaped Networks of Strings[END_REF] and [START_REF] Ammari | Méthode numérique pour la décroissance de l'énergie d'un réseau de cordes[END_REF] and will be detailed in [START_REF] Jellouli | Optimal decay rates for the stabilization of a string network[END_REF]. In fact we recall that

E(t) = E(0) -α t 0 |∂ t u 1 (s, 0)| 2 ds (6)
and that for all t ≥ 0,

∂ t u 1 (t, 0) = P • P -1 F (t) + 1 α+1 L + β,ℓ,1 a ′ 1 + L - β,ℓ,1 b 1 , β = α-1 α+1
. In a first step, we evaluate the caracteristic elements of P -1 and we prove that P -1 ≡ (n, δ), where n : n(j) = 2(j -1)ℓ , j ≥ 1 and in the case α = α 0

δ : δ(j) = a(β)λ j-2 + b(β)µ j-2 , j ≥ 2 with µ = - 2(N -2) α+1 - √ ∆ 2N , a(β) = 1 N (1-λ)(β-λ) λ-µ and b(β) = -1 N (1-µ)(β-µ) λ-µ
. In the case α = α 0 ,

δ(j) = δ(2) + 1 λ 0 δ(3) -δ(2) (j -2) λ j-2 0 , ∀j ≥ 2.
In the second stage we calculate 2nℓ 0

|∂ t u 1 (s, 0)| 2
ds, what will allow to find E ∞ :

E ∞ = 1 2 N j=2 a ′ j 2 + b j 2 - 1 4(N -1)   N j=2 (a ′ j -b j ) 2 + N j=2 (a ′ j + b j ) 2   = 1 2(N -1) N j=2 N k=j+1 (a ′ k -a ′ j ) 2 + (b k -b j ) 2 .
In the end, we estimate the difference E(t = 2nℓ) -E ∞ and we show the inequality (4) and ( 5).

1. √ ∆ = i √ -∆ if ∆ < 0. Note that, ∆ < 0 (resp. ∆ = 0) if and only if α ∈]0, α 0 [ (resp. α = α 0 )

Semi-discretization and numerical results

Let N j , j = 1, . . . , d be positive integers. We define h = (h 1 , . . . , h d ), with h j = ℓj Nj , j = 1, . . . , d. We suppose that α > 0 and consider the following classical finite difference space semi-discretization of (S):

                             u ′′ j,k (t) - u j,k-1 (t) -2u j,k (t) + u j,k+1 (t) h 2 j = 0, k = 1, . . . , N j -1, j = 1, . . . , d, u 1,1 (t) -u 1,0 (t) h 1 = αu ′ 1,0 (t), u 1,N1 (t) -u 1,N1-1 (t) h 1 = d j=2 u j,1 (t) -u j,0 (t) h j , u 1,N1 (t) = u j,0 (t) and u j,Nj (t) = 0, j = 2, . . . , d, u j,k (0) = u 0 j,k , u ′ j,k (0) = u 1 j,k , j = 1, . . . , d, k = 0, . . . , N j , t ∈ (0, ∞), (7)
where the initial conditions satisfy for s = 0, 1, the compatibility conditions

u s 1,N1 = u s k,0 = d j=2 u s j,1 hj + u s 1,N 1 -1 h1 d j=1 1 hj and u s k,N k = 0, k = 2, . . . , d.
The solution of ( 7) can be written in the form V (t) = exp(tA)φ, where

V = (U, U ′ ) t , U = (u 1,0 , u 1,1 , . . . , u 1,N1-1 , u 2,1 , . . . , u 2,N2-1 , . . . , u d,1 , . . . , u d,N d -1 ) ∈ R M
and φ ∈ R 2M . The semi-discrete energy is defined by

E h,φ (t) = 1 2 d j=1 h j   Nj -1 k=1 u ′ j,k (t) 2 + Nj -1 k=0 u j,k+1 (t) -u j,k (t) h j 2   , (8)
and is decreasing as it satisfies E ′ h,φ (t) = -α u ′ 1,0 (t) 2 . The following proposition identifies the limit of the discrete energy, which may not be zero, as in the continuous case. Proposition 3.1 For φ ∈ R 2M , we can write φ = φ 1 +φ 2 where φ 1 belongs to the space Λ 1 of eigenvectors relative to the eigenvalues λ of A, with ℜ(λ) = 0 and φ 2 belongs to the space Λ 2 of generalized eigenvectors relative to the eigenvalues λ of A, with ℜ(λ) < 0. We then have

E h,φ,∞ := lim t→∞ E h,φ (t) = E h,φ1 (0), (9)
and this value does not depend on α > 0. Sketch of the proof We first proof that for φ 1 ∈ Λ 1 ⊂ C 2M , we have φ 1,M +1 = 0, that no eigenvalue of A can have positive real part and that generalized eigenvectors with eigenvalue of zero real part are necessarily eigenvectors, using the fact that the energy is decreasing and nonnegative. We then use the decomposition and estimates ≤ 2d E h,φ1 (t) E h,φ2 (t), where w 1 j,k (resp. w 2 j,k ) is the solution corresponding to the initial condition φ 1 (resp. φ 2 ). We know that E h,φ1 (t) = E h,φ1 (0), thus we get (9) since E h,φ2 (t), being a finite sum and product of terms tending to zero, as t → +∞, also tends to zero. This limit does not depend on α, since the only equation where α is present is

|E h,φ (t) -E h,φ1 (t) -E h,φ2 (t)| = d j=1 Nj -1 k=1 h j (w 1 j,k ) ′ (w 2 j,k ) ′ + d j=1 Nj -1 k=0 1 h j (w 1 j,k+1 -w 1 j,k )(w 2 j,k+1 -w 1 j,k ) ≤ d j=1 h j Nj -1 k=1 (w 1 j,k ) ′ 2 Nj -1 k=1 (w 2 j,k ) ′ 2 + 1 h j Nj -1 k=0 w 1 j,k+1 -w 1 j,k 2 Nj -1 k=0 w 2 j,k+1 -w 2 j,k 2
αλφ λ,M +1 = 0 = 1 h 1 (φ λ,M +2 -φ λ,M +1 ) , for each λ; ℜ(λ) = 0.
Numerical results. Finally we give some exact and semi-discrete results for the energy. We take d = 3, initial condition

u 0 j (x) = 0, j = 1, 3, u 0 2 (x) = sin 2 ( πx ℓ 2
) and u 1 j (x) = 0, j = 1, 2, 3.

ℓ j = 1 and N j = N, j = 1, 2, 3. We see on Figure 1 (left) the time evolution of E(t) -E ∞ for α = 1 (exact, N = 100, 400) and for α = α opt , 0.9 (exact). On Figure 1 (right), we plot the time evolution of E(t) -E ∞ (exact, N = 100, 400) taking as before u 0 2 (x) = sin 2 (πx) (sin) but also u 0 2 (x) = sin 2 (8πx) (sin8), in the case α = α opt . The results are conform with Theorem 2.2: taking α = α opt leads to the best growth rate with asymptotic behaviour of the form Ct 2 exp(-γt). We also note that the semi-discrete scheme is able to reproduce the behaviour of the exact scheme until the discretization error dominates, which is higher by taking a higher mode (Figure 1, right). The semi discrete energy seems to converge to a discretization dependent limit, which is coherent with Proposition 3.1.

Figure 1 .

 1 Figure 1. E(t)-E∞ vs time (exact and semi-discrete case) for α ∈ {0.9, 1, αopt} (left) and different initial conditions (right).


Robert Cori
email: robert.cori@labri.fr

Michel Marcus
email: michel.marcus@labri.fr

Gilles Schaeffer
email: schaeffe@lix.polytechnique.fr

Antonio Machì

Odd permutations are nicer than even ones ✩

We give simple combinatorial proofs of some formulas for the number of factorizations of permutations in S n as a product of two n-cycles, or of an n-cycle and an (n -1)-cycle.

Introduction

The parameter number of cycles plays a central role in the algebraic theory of the symmetric group, however there are very few results giving a relationship between the number of cycles of two permutations and that of their product.

The first results on the subject go back to O. Ore, E. Bertram, R. Stanley (see [START_REF] Ore | Some remarks on commutators[END_REF], [START_REF] Bertram | Even permutations as a product of two conjugate cycles[END_REF], [START_REF] Stanley | Factorization of permutations into n-cycles[END_REF]), who proved some existence theorems. These results allowed to obtain informations on the structure of the commutator subgroup of the alternating group.

After these pionneering works, enumerative results were investigated probably for the first time by D. Walkup (in [START_REF] Walkup | How many ways can a permutation be factored into two n-cycles[END_REF]). In order to obtain enumeration formulas, different techniques were used. Some of them are purely combinatorial, while others use character theory. The key point of this second approach is a theorem, often attributed to Frobenius, expressing linearization coefficients for the product of conjugacy classes in the group algebra, as a sum of products of characters of irreducible representations.

The use of character theory was very fruitful to many authors and gave rise to many papers by R. Stanley, D. Jackson, D. Zagier G. Jones, A. Goupil, G. Schaeffer, D. Poulalhon, P. Biane ([15], [START_REF] Jackson | Counting cycles in permutations by group characters with an application to a topological problem[END_REF], [START_REF] Jones | Trees, permutations and characters[END_REF], [START_REF] Zagier | On the distribution of the number of cycles of elements in symmetric groups[END_REF], [START_REF] Goupil | Factoring n-cycles and counting maps of given genus[END_REF], [START_REF] Poulalhon | Factorizations of large cycles in the symmetric group[END_REF], [START_REF] Biane | Nombre de factorisations d'un grand cycle[END_REF]).

The main results using combinatorial methods were obtained by G. Boccara in a paper [START_REF] Boccara | Nombre de représentations d'une permutation comme produit de deux cycles de longueurs données[END_REF] containing many results, but using also some integrals of polynomials with a lot of computations, that can be considered as ugly by pure combinatorists. The most surprizing result in Boccara's paper is that for any odd permutation in S n , the number of different ways to write it as a product of an n-cycle and an (n -1)-cycle is 2(n -2)!. This formula was puzzling, and M. P. Schützenberger asked for a bijective proof for it. Antonio Machì found such a proof (in [START_REF] Machì | On the representation of odd permutations[END_REF]), giving a bijection, in which cycles of length 3 play a central role, while generally combinatorialists use products of transpositions (i.e. cycles of length 2) to represent permutations. His approach was applied by some of his students (see [START_REF] Cangelmi | Factorizations of an n-cycle into two n-cycles[END_REF] and [START_REF] Cangelmi | A representation of even permutations[END_REF]), obtaining new enumerative results.

G. Boccara's method was recently used by R. Stanley (in [START_REF] Stanley | Two enumerative results on cycles of permutations[END_REF]) in order to solve a conjecture of M. Bóna (mentioned in [START_REF] Bóna | The average number of block interchanges needed to sort a permutation and a recent result of Stanley[END_REF]), and to give another proof of a result obtained by D. Zagier (see [START_REF] Zagier | On the distribution of the number of cycles of elements in symmetric groups[END_REF], application 3 of Theorem 1). In his paper R. Stanley asks for purely combinatorial proofs of these results. Our paper gives some answers to his questions. Another approach to the same problem was developed independently by V. Feray and E. Vassilieva (see [START_REF] Feray | Linear coefficients of Kerov's polynomials: bijective proof and refinement of Stanley's result[END_REF]).

Our approach is purely combinatorial, it consists in rederiving combinatorialy Boccara's main result, and using this central theorem to derive many other results. The paper is organized as follows: after some definitions and notation in the first part, the second one gives a bijective proof of Boccara's main result which is different from Machì's one. In the third section we show how to deduce results for even permutations from those for odd ones, thus illustrating the title of our article, we also give a new proof of Zagier's result. The forth section is devoted to a new, purely combinatorial proof of Bóna's conjecture. We conclude in the fifth section with a new proof for a formula of Boccara (see [START_REF] Boccara | Nombre de représentations d'une permutation comme produit de deux cycles de longueurs données[END_REF] Corollary 4.8) enumerating factorizations of even permutations.

Notation and definitions

A permutation α of S n will be considered here as a bijection from {1, 2, . . . , n} on to itself, so that α(i) will denote the image of i by the permutation α. The convention we use is that the product αβ of the two permutations α and β is the permutation mapping i on α(β(i)). Permutations will be represented by their cycles like:

α = (i 1 , α(i 1 ), . . .)(i 2 , α(i 2 ), . . .) • • • (i m , α(i m ), . . .).
The number of cycles of a permutation α will be denoted z(α). Hence α is an n-cycle of S n if and only if z(α) = 1; moreover if β is an (n -1)-cycle of S n then z(β) = 2. Recall that a permutation α of S n is odd if it is the product of an odd number of transpositions, or equivalently if n + z(α) is an odd number, it is even otherwise. Definition 1. A factorization of a permutation σ of S n into two large cycles is a pair of permutations α, β such that: σ = αβ, the permutation α is an n-cycle, and the permutation β is an n-cycle when σ is even and an (n -1)-cycle when σ is odd.

Example 1. There are 4 factorizations of the permutation σ = (1, 2, 3, 4) into two large cycles, they are:

α β α β (1, 3, 2, 4) (1, 3, 2) (4) (1, 2, 4, 3) (2, 4, 3) (1) (1, 4, 2, 3) (1, 4, 3) (2) (1, 3, 4, 2) (1, 4, 2) (3)
We use the following classical notation and definitions on partitions of an integer: A partition of the integer n is a finite sequence of postive integers λ = (λ 1 , λ 2 , . . . , λ m ) such that:

λ 1 ≥ λ 2 ≥ . . . ≥ λ m > 0, and m i=1 λ i = n.
We will write λ ⊢ n and define the length of λ as ℓ(λ) = m. A partition λ ⊢ n will be said to be odd if n + ℓ(λ) is an odd number it to be even otherwise.

The cyclic type of a permutation is the sequence of the lengths of its cycles written in weakly decreasing order, defining hence a partition of n, whose length is the number of cycles of the permutation. Note that a permutation has the same parity as the partition given by its cyclic type. Definition 2. An ELC-factorization of type λ (even factorization into large cycles) is a pair of n-cycles (α, β) such that αβ = σ is of cyclic type λ.

An OLC-factorization of type λ (odd factorization into large cycles) is a pair (α, β) such that, α is an n-cycle, β an (n -1)-cycle and αβ = σ is of cyclic type λ.

The sets of ELC-factorizations and of OLC-factorizations of type λ are denoted E(λ) and O(λ) respectively. We will denote by E n the set of all even factorizations into two large cycles of S n , and O n that of all odd factorizations into two large cycles of S n . Clearly the number

|E n | of elements in E n is equal to (n -1)! 2 and |O n | = n!(n -2)!.

Odd factorizations into two large cycles

In this section we give a new combinatorial proof of the following central result: Theorem 1. For any odd permutation σ in S n the number of OLC-factorizations (α, β) such that αβ = σ is equal to 2(n -2)!.

The proof we give here uses a technique introduced by A. Lehman in order to enumerate the number of maps with one vertex and one face embedded on an orientable surface [START_REF] Lehman | Enumeration of one-vertex one-face rooted maps[END_REF] 1 . The key point is to represent factorizations (α, β) of a permutation σ in 2 large cycles by using a sequence of integers u = b 0 , b 1 , b 2 , . . . , b n-1 , which represents the permutation β, and a directed graph G σ,u which partially represents the permutation α.

A permutation α of S n is classically represented by a (directed) graph, which we denote G α , it has {1, 2, . . . , n} as vertex set and contains n arcs, each one has head i and tail α(i) for 1 ≤ i ≤ n. The circuits of G α correspond to the cycles of α.

When

σ = αβ where β = (b 1 , b 2 , . . . , b n-1 )(b n ) is an (n -1)-cycle, then (since α = σβ -1
), the graph G α may be obtained from σ, by joining b i to σ(b i-1 ) for i = 2, . . . , n-1, then adding an arc form b 1 to σ(b n-1 ) and an arc from b n to σ(b n ). This remark leads to an algorithm allowing to compute factorizations of σ into two large cycles. Lehman's contribution consists in adding an element b 0 , allowing to obtain enumeration formulas. ) is the tail of an arc. Notice that no vertex can be the head of two arcs, however σ(b 0 ) may be the tail of either one or two arcs.

Example 2. The 6 Lehman sequences of the permutation (1, 2, 3) are given below with their graphs: ) are the tail of no arc. The other vertices are the head of exactly one arc, and the tail of exactly one. This graph is the union of two paths, where one may be reduced to a single vertex. We have thus one of the two situations depicted in Figure 2: to G ′ σ,u has two circuits, and the permutation α has two cycles. This implies that σ is the product of two permutations each one having two cycles, contradicting the fact that σ is odd. Hence only the second case has to be considered and α(b

3 1 2 3 1 2 3 1 2 3 3, 2, 1 3, 3, 2 2, 1, 3 1 2 3 1, 1, 3 1, 3, 2 2, 2, 1 1 2 3 1
1 ) = σ(b n-1 ), α(b n ) = σ(b n ), proving that α has only one cycle. ✷ Proposition 2.
The number of Lehman sequences of a permutation σ of S n is equal to n!.

Proof : One may build all the Lehman sequences of σ using the following algorithm, which proceeds as an iteration of n steps:

• Choose any b 0 among {1, 2, . . . , n}, and b 1 different from σ(b 0 ) in the same set,

• For i = 2, . . . n -1 choose b i in {1, 2, . . . , n} different from b 1 , b 2 , . . . , b i-1 and from b ′ i
, the vertex which is the tail of the longest path with head σ(b i-1 ) in the graph G σ,b0b1...bi-1 .

Notice that at each step i the number of possible choices for b i is ni. Indeed it is clear that their are n possible choices for b 0 and n -1 choices for b 1 ; moreover since b ′ i is a vertex of out-degree 0 of G σ,b0b1...bi-1 , it cannot be equal to any of the b j for 1 ≤ j < i, as these vertices are all of out-degree 1. Hence at each step the number of possible choices decreases by one, then the number of sequences built by the algorithm is n!.

Since the algorithm closely follows the definition of a Lehman sequence it is clear that it builds exactly all the Lehman sequences of σ. ✷ So that the graph obtained has shape given in Figure 3.

For each choice of b 1 the number of vertices of the path P from σ(b n-1 ) to b n takes any value from 1 to (n -1), and all these values are obtained. Now b 0 has to be chosen such that the arc (b 1 , σ(b 0 )) does not create a circuit in the graph; hence σ(b 0 ) must be among the vertices of this path P . This gives a total number of choices equal to sum of the lengths of these paths for all the possible choices of b 1 , that is 1

+ 2 + • • • + (n -1). Hence the number of Lehman sequences u such that Λ(u) = β is n(n-1)

2

.

✷

The proof of Theorem 1 then follows from Propositions 1, 2 and 3.

From odd factorizations to even ones

In this section we show how to build O n+1 from E n . We express for a given odd partition µ ⊢ n + 1, the set O(µ) as a union of sets obtained from the E(λ) for some partitions λ ⊢ n closely related to µ. We will then use this result to determine the number of factorisations of even permutations with k cycles, giving a new combinatorial proof of a result of D. Zagier.

Building factorizations of even permutations

We will use transpositions, i. e. permutations of S n with n -2 fixed points and one cycle of length 2; a transposition will be denoted (i, j), where i, j are the two elements which are not fixed.

Remark 1. Let α be a permutation in S n and let i be an integer such that 1 ≤ i ≤ n, then α(i, n + 1) is a permutation of S n+1 obtained from α inserting n + 1 after i in its cyclic representation and (i, n + 1)α is a permutation of S n+1 obtained from α inserting n + 1 before i in its cyclic representation.

Let (α, β) be an element of E n and let i, j be two integers such that 1 ≤ i ≤ n and 1 ≤ j ≤ n + 1. We define the factorization Φ i,j (α, β) = (α ′ , β ′ ) in O n+1 in the following way:

• Insert n + 1 in the cyclic representation of α immediately before i giving a permutation α 1

• Let β 1 ∈ S n+1 be equal to the permutation β to which is added n + 1 as a fixed point.

• Exchange j and n + 1 in α 1 and in β 1 obtaining α ′ and β ′ .

More precisely we have, using the fact that exchanging n + 1 and j is obtained by conjugation by (j, n + 1):

α ′ = (j, n + 1)(i, n + 1)α(j, n + 1), β ′ = (j, n + 1)β(j, n + 1) (1)
Definition 5. Let λ be a partition of n and µ a partition of n + 1; λ is covered by µ, if they have the same length and if λ i = µ i for all 1 ≤ i ≤ ℓ(µ), except for only one value of i. This will be denoted by λ < • µ Remark 2. When λ < • µ then for the unique k such that λ k = µ k we have

µ k = λ k + 1.
Moreover, for a given µ the number of different λ such that λ < • µ is equal to the number of different values taken by the µ i that are larger than 1.

For instance µ = (4, 3, 3, 1) covers 2 partitions: λ = (3, 3, 3, 1) and λ = (4, 3, 2, 1).

Theorem 2. Let (α, β) be an even factorization of type λ ⊢ n. For any 1 ≤ i ≤ n and

1 ≤ j ≤ n + 1, Φ i,j (α, β) is a factorization which type µ is such that λ < • µ. Conversely let (α ′ , β ′ ) ∈ O(µ) with µ ⊢ n + 1, then there exist two integers i, j and a factorization (α, β) ∈ E n such that 1 ≤ i ≤ n, 1 ≤ j ≤ n + 1 and (α ′ , β ′ ) = Φ i,j ((α, β)).
Moreover i, j, α, β are uniquely determined.

Proof : Considering the first part of the theorem, let σ = αβ and Φ i,j ((α, β)) = (α ′ , β ′ ). It is easy to check that α ′ is an (n + 1)-cycle, since (i, n + 1)α is an (n + 1)-cycle, and conjugation does not modify the cyclic type. For the same reason β ′ is an n-cycle, which fixed point is j. The type of Φ i,j (α, β) is the cyclic type of the permutation (j, n + 1)(i, n + 1)αβ(j, n + 1) but, since conjugation does not modify the cyclic type of a permutation, this is also the cyclic type of σ ′ = (i, n + 1)αβ. But from Remark 1 this last permutation is obtained from αβ inserting n + 1 before i in its cycle, showing that its cyclic type is a partition µ covering λ.

We now prove the converse. Let (α ′ , β ′ ) ∈ O(µ) and let j be the point fixed by β ′ . Consider the formula giving α ′ in (1):

α ′ = (j, n + 1)(i, n + 1)α(j, n + 1)
We get α ′ (j) = n + 1 if i = j, and α ′ (j) = i if j = i; this allows to determine uniquely i from α ′-1 (j). We then obtain (α, β) by equation ( 1). Since n + 1 is a fixed point in both α and β, these two permutations may be considered as elements of S n . ✷

Permutations with a given number of cycles

We are now able to give a combinatorial proof of the following result due to D. Zagier (see [START_REF] Zagier | On the distribution of the number of cycles of elements in symmetric groups[END_REF], application 3 of Theorem 1.) Corollary 1. Let n, k be two integers such that 1 ≤ k ≤ n and nk is even. The probability that the product of two n-cycles of S n (taken randomly) has k cycles is equal to the probability that an odd permutation of S n+1 has k cycles.

Proof : We compute the number e n,k of factorizations (α, β) ∈ E n such that αβ has k cycles, for that we consider first the number of elements (α ′ , β ′ ) in O n+1 such that α ′ β ′ has the same number k of cycles. By Theorem 1, this number is equal to 2(n -1)! times the number of elements of S n+1 with k cycles. That is:

2(n -1)!s n+1,k
where s n,p is the unsigned Stirling number of the first kind.

By Theorem 2, for any factorization of (α, β) ∈ E n we obtain by Φ i,j , n(n+1) different factorizations (α ′ , β ′ ) ∈ O n+1 . For all of them, α ′ β ′ has k cycles, and any element of O n+1 with k cycles is obtained in this way. Hence:

e n,k = 2(n -1)!s n+1,k n(n + 1)
Since the number of elements of E n is (n -1)! 2 , the probability for the product of two n-cycles of S n (taken randomly) has k cycles is:

e n,k (n -1)! 2 = 2s n+1,k (n + 1)!
Now the number off odd permutations in S n+1 is (n+1)! 2 and s n+1,k of them have k cycles, proving the result. ✷

Permutations in which two given elements are in the same cycle

In this section we consider the following question conjectured by M. Bóna and solved by R. Stanley [START_REF] Stanley | Two enumerative results on cycles of permutations[END_REF]; Proposition 5 below was obtained in collaboration with Valentin Feray and Amarpreet Rattan. Definition 6. A factorization (α, β) into two large cycles of S n for n ≥ 2 is said to be connecting if 1 and 2 are in the same cycle of αβ.

Let p n be the number of connecting factorizations of E n and q n be the number of those in O n . It is easy to compute q n : Proposition 4. The number of connecting factorizations in O n is given by:

q n = (n -2)! n! 2 + (-1) n (n -2)! (2)
Proof : We first begin by computing the numbers a n and b n of even and odd permutations of S n such that 1 and 2 are in the same cycle. Any even permutation of S n is obtained either from an even one in S n-1 adding n as a fixed point or from an odd one in S n-1 adding n inside a cycle. Since this last operation can be done in n -1 ways, and since a similar statement holds for odd permutations of S n we have for n ≥ 3:

a n = a n-1 + (n -1)b n-1 , b n = b n-1 + (n -1)a n-1
the difference of these two numbers is equal to:

a n -b n = (2 -n)(a n-1 -b n-1 )
Multiplication by the transposition (1, 2) is a bijection between permutations in S n where 1 and 2 are in the same cycle and those where they are not, hence giving;

a n + b n = n! 2
Each odd permutation with 1 and 2 in the same cycle gives 2(n -2)! connecting factorizations in O n hence:

q n = 2(n -2)!b n = 2(n -2)! n! 4 + (-1) n (n -2)! 2
✷ We again find results for even factorizations using our result for odd ones.

Proposition 5. The number of odd connecting factorizations in S n+1 and that of the even ones in S n are related by:

q n+1 = (n + 2)(n -1)p n + 2(n -1)! 2 (3)
Proof : We divide the set of connecting odd factorizations (α ′ , β ′ ) in S n+1 into three subsets X n+1 , Y n+1 and Z n+1 according to the value of the fixed point j of β ′ . Those in X n+1 , Y n+1 have j = 1 and j = 2 respectively, and those in Z n+1 are such that j > 2.

We use Theorem 2 to compute the number z n+1 of elements in Z n+1 . Each factorization (α ′ , β ′ ) in this set is equal to a Φ i,j (α, β) such that (α, β) is a connecting even factorization in E n , and j = 3 ≤ j ≤ n + 1. This gives n -1 possible values for j and since there are n possible values for i we have:

z n+1 = n(n -1)p n .
The computation of the numbers x n+1 , y n+1 (which are equal by symmetry) of elements in X n+1 and Y n+1 is more intricate. Notice that if (α ′ , β ′ ) = Φ i,j (α, β) is such that β ′ (1) = 1 then j = 1 and (α, β) is not necessarily connecting.

Nevertheless we can obtain an inductive formula for x n+1 . We have for (α ′ , β ′ ) ∈ X n+1 :

α ′ β ′ = (1, n + 1)(i, n + 1)αβ(1, n + 1)
for a given i. In order that (α ′ , β ′ ) be connecting we must have one of the two following conditions:

• i = 2
• i = 2 and (α, β) is such that i and 2 are in the same cycle of αβ.

If i = 2, any factorization in E n will give a connecting factorization, the number of such is (n -1)! 2 .

If i = 2, we have to choose a factorization such that i and 2 are in the same cycle of α ′ β ′ , and the number of such factorizations is p n by symmetry. Since there are n -1 possible choices for i this gives x n+1 = (n -1)! 2 + (n -1)p n .

Hence using x n+1 = y n+1 we get:

q n+1 = x n+1 + y n+1 + z n+1 = 2(n -1)! 2 + 2(n -1)p n + n(n -1)p n .
✷ Corollary 2. Taking two n-cycles at random the probability that their product contains 1 and 2 in the same cycle is 1 2 if n is odd and 1 2 -

2 (n+2)(n-1)
if n is even.

Proof : From the two Propositions above we get:

(n + 2)(n -1)p n = (n -1)! 2 2 n 2 + n -2(-1) n -4
Since n 2 + n -2(-1) n -4 = (n + 2)(n -1) -2((-1) n + 1) we obtain:

p n = (n -1)! 2 2 1 -2 (-1) n + 1 (n + 2)(n + 1)
The probability we are seeking for is pn (n-1)! 2 completing the proof. ✷

A formula for the number of elements of E(λ)

We have proved in Section 2 that the number of factorizations of an odd permutation of S n into two large cycles depends only on n. For an even permutation σ this number depends on the conjugacy class of σ which is determined by a partition λ = (λ 1 , λ 2 . . . , λ m ) of n. A formula was given for these numbers by G. Boccara (in [START_REF] Boccara | Nombre de représentations d'une permutation comme produit de deux cycles de longueurs données[END_REF] Corollary 4.8), we give a new proof of this formula here based on Theorem 2.

We give first some notation for such a partition and introduce the function P (λ) which is central in this enumeration.

• Let I be the set of indexes of the elements of λ namely: I = {1, 2, . . . , m}.

• For any subset X of I, denote s X (λ) the sum of the parts of λ whose indexes are in X :

s X (λ) = i∈X λ i
• If µ has k parts equal to 1 then we procced by induction. Above we have proved k = 0. Let µ be a a partition with k parts equal to 1. Let µ ′ be obtained from µ deleting one part equal to 1. By the induction hypothesis we have

λ ′ <• µ ′ c λ ′ u(λ ′ , µ ′ )P (λ ′ ) = 2(n -1)!c µ ′
But there are as many partitions λ ′ such that λ ′ < • µ ′ as those λ such that λ < • µ each of those λ ′ is obtained from a λ deleting one part equal to 1. For these we have

c λ ′ = c λ k n . Moreover c µ ′ = cµ k n+1 and u(λ ′ , µ ′ ) = u(λ, µ). Giving: λ ′ <• µ ′ c λ k n u(λ, µ)P (λ ′ ) = c µ k n + 1 2(n -1)!
Then by Lemma 1, we have (n + 1)P (λ ′ ) = P (λ) and the result. ✷

Proof of Theorem 3

Let λ be an even partition of n and let µ be such that λ < • µ, where µ k = λ k + 1 then we denote by u(λ, µ) the sum of the parts of size λ k in λ, that is u(λ, µ) = i λ k λ k where i p denotes the number of λ j equal to p. Then we will use the following Lemma giving an elementary proof of a formula allowing to compute the e λ = |E(λ)| due to V. Feray and E. Vassilieva [START_REF] Feray | Linear coefficients of Kerov's polynomials: bijective proof and refinement of Stanley's result[END_REF]. This lemma may be considered as a Corollary of Theorem 2. Lemma 3. Let µ be an odd partition of n + 1 then:

2c µ (n -1)! = (n + 1) λ<• µ u(λ, µ)e λ (10)
Where c µ is the number of permutations with cyclic type µ and the range of the sum in the right hand side consists of all the partitions λ covered by µ.

Proof : The left hand side of the above equality is the number of factorizations (α ′ , β ′ ) in O n+1 such that α ′ β ′ has cyclic type µ. By Theorem 2, each of such factorization is the image by some Φ i,j of a factorization (α, β) in E n , such that the cyclic type of αβ is a partition λ covered by µ.

Conversly let λ be a partition covered by µ such that λ k +1 = µ k and let (α, β) ∈ E(λ). Then Φ i,j (α, β) ∈ O(µ) if and only if i belongs to a cycle of length λ k in the permutation αβ. Since there are u(λ, µ) elements in {1, 2, . . . , n} belonging to a cycle of length λ k of αβ the number of possible values for i is u(λ, µ) and the number of possible values for j is n + 1. Hence the number of factorization in O(µ) obtained as a Φ i,j (α, β) for (α, β) ∈ E(µ) is (n + 1)u(λ, µ)a λ proving the result. ✷ Consider equations [START_REF] Jones | Trees, permutations and characters[END_REF] for all µ ⊢ (n + 1) of length m. For m = 1 there is only one equation allowing to obtain for n odd:

e (n) = 2(n -1)! 2 n + 1 .
For a given m > 1 each partition µ = (µ 1 , µ 2 , . . . µ m ) of length m gives an equation:

(n + 1) λ<• µ u(λ, µ)e λ = 2c µ (n -1)!.

In this equation all the partitions λ appearing in the sum, except one of them, contain the same number of highest parts equal to µ 1 . Using induction on the value of the highest part λ 1 of λ, this proves that these equations enable to obtain uniquely the e λ for the even partitions λ of length m. The following result may be useful for the computation of e λ when some of the λ i are equal to 1.

Remark 3. Let λ be an even partition of length m with k > 0 parts equal to 1, let λ ′ be obtained from λ deleting a part of size 1, then the following relation allows to compute e λ inductively: e λ = e λ ′ n(n -1) k

Proof : Notice that from any factorization (α ′ , β ′ ) in E(λ ′ ) one builds n -1 factorizations in E(λ) inserting n before an i in the cycle of β ′ and after that element i in the cycle of α ′ . The factorizations (α, β) obtained in this way are such that n is a fixed point of αβ. Moreover any factorization (α, β) of cyclic type λ and such that n is a fixed point of αβ can be obtained in that way. This proves that the number of these factorizations is equal to (n -1)e λ ′ . Now given k integers 1 ≤ i 1 , i 2 , . . . , i k ≤ n, the number of factorizations (α, β) in E(λ) such that i 1 , i 2 , . . . , i k are fixed points of αβ is independent of i 1 , i 2 , . . . , i k giving :

Definition 3 .

 3 Let σ be a permutation of S n , and let u = b 0 , b 1 , . . . , b k be a sequence of positive integers not greater than n and such that b 1 , b 2 , . . . , b k are all distinct. The graph G σ,u associated to σ, u has for vertices the integers {1, 2, . . . , n} and an arc for each j, 1 ≤ j ≤ k with head b j and tail σ(b j-1 ). Definition 4. A sequence of integers u = b 0 , b 1 , . . . , b n-1 is a Lehman sequence of the permutation σ in S n , if no two b i for i ≥ 1 are equal and if the graph G σ,u has no circuit. Let u be a Lehman sequence of the permutation σ in S n , then the graph G σ,u has n vertices and n -1 arcs, each vertex b 1 , b 2 , . . . b n-1 is the head of an arc and each vertex σ(b 1 ), σ(b 2 ), . . . , σ(b n-2

2 Figure 1 : 3 )Proposition 1 .

 2131 Figure 1: The graphs of the 6 Lehman sequences of σ = (1, 2, 3)

Figure 2 :

 2 Figure 2: The two situations in the proof of Proposition 1.

Proposition 3 . 2 .Figure 3 :

 323 Figure 3: Deletion of two arcs in Gα

  Now let us compare with what we obtained for P (λ) in Proposition 6:λ<• µ c λ u(λ, µ)P (λ) = 2 n!c µThis clearly shows that the c λ P (λ) n(n+1) satisfy the same set of equations as the e λ . Since this set of equations has a unique solution this showsc λ P (λ) n(n + 1) = e λ ✷

This is a manuscript he gave to his students which we obtained thanks to Timothy Walsh

Acknowledgments

The results in Section 4 were obtained in collaboration with Valentin Feray and Amarpreet Rattan.

The authors wish to thank the referee for very hepfull comments concerning the presentation of this paper.

✩ The authors acknowledge the support of ERC under the agreement "ERC StG 208471 -ExploreMap"

• The number P (λ) is the sum given below, where |X| is the number of elements of X: P (λ) = X⊆I (-1) s X (λ)+|X| s X (λ)!(ns X (λ))! (4)

Let c λ denotes the number of permutations of cyclic type λ and e λ the number of elements in E(λ).

Theorem 3. For any even partition λ:

5.1. Some properties of the function P (λ) In order to prove this Theorem we need to prove some properties of the function P (λ).

Lemma 1. If λ has length 1 we have:

Let λ ⊢ n be an even partition with at least one part equal to 1, and let λ ′ ⊢ n -1 be the partition obtained from it deleting one of these parts. Then:

Proof : For λ = (n) the sum giving P (λ) reduces to two elements with X = ∅ and X = I, for each of these we have: s ∅ (λ) = 0 and s I (λ) = n giving (6).

To prove [START_REF] Feray | Linear coefficients of Kerov's polynomials: bijective proof and refinement of Stanley's result[END_REF] we notice that the sum giving P (λ) have twice the number of terms than that giving P (λ ′ ). Each X giving the term u ′

) with the same sign as u ′ X . These two terms sum up to:

Hence ending the proof. ✷ Lemma 2. Let µ ⊢ (n + 1) be an odd partition of length m, where all the m parts are strictly greater than 1. Let for 1 ≤ i ≤ m, λ (i) be the sequence obtained from µ deleting one element in part i. Then for any non empty X strictly included in {1, 2, . . . , m} we have

where s i X = j∈X λ (i) j .

Proof : Let p = s X (µ), since λ (i) i = µ i -1 and λ (i) j = µ j for i = j we have:

But since p = j∈X µ j we have j / ∈X µ j = n + 1p and the result. ✷ Proposition 6. Let µ be an odd partition of n + 1 then:

Proof : We distinguish two cases:

• If µ has no part of length 1. Then we use Lemma 2 and add equations ( 8) multiplied by (-1) |X| for all subsets X of I = {1, 2, . . . , m}. For X = I we have for all i: s i X = n and |X| + s i X = n + m, which is even. For X = ∅ we have s i X = 0 and |X| + s i X = 0, hence we get (since

Reordering any sequence λ (i) we obtain a partition λ such that λ < • µ. Moreover for a given λ such that λ < • µ, the number of the sequences λ (j) which give λ after reordering it, is equal to n i , the number of occurrences of µ i in µ (where i is the index such that λ i = µ i -1).

Hence reversing the order of the two sums in the equation above we obtain: