L Garrigues
email: laurent.garrigues@laplace.univ-tlse.fr.

J Pérez-Luna

Juslan Lo
email: juslan.lo@laplace.univ-tlse.fr.

J M Hagelaar
email: gerjan.hagelaar@laplace.univ-tlse.fr.

J P Boeuf

S Mazouffre
email: stephane.mazouffre@cnrs-orleans.fr.

G J M Hagelaar

Determination of the Axial Electron Mobility Profile in the PPS®X000-LM Thruster

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Electron transport across the magnetic field in

I. Introduction

n a Hall Effect Thruster (HET), the thrust is provided by the acceleration of xenon ions that originates in the efficient ionization of a xenon neutral flow in a E×B field discharge. In HETs, the ion beam velocity can reach few 10's of km.s -1 in comparisons with the velocity of the exhaust beam obtained in chemical thrusters limited to few km.s -1 . The HETs offer considerable advantages for geostationary satellites for station-keeping and orbit control 1 and for exploration probes 2 because the very high exhaust beam velocity leads to less propellant consuming. A new generation of 5kW-class Hall thrusters is now developed in order to achieve orbital manoeuvres for the next generation of large geosynchronous satellites. [START_REF] Mazouffre | Expanding Frontiers : Towards High Power Hall Effect Thrusters for Interplanetary Journeys[END_REF] Nevertheless, the physic involving in such kind of thruster is far from being completely mastered. Especially, there is not yet a clear theory of the so-called "anomalous electron mobility". I Indeed, in a HET configuration, electrons travel from the cathode to the anode through a E×B field region. The cross field electron mobility µ e,⊥ is proportional to the momentum-transfer frequency of electron-particle collisions.

Due to the intense ionization in the thruster channel, the neutral density is strongly depleted in the exhaust region. The momentum transfer frequency only due to collisions between electrons and neutrals is consequently low, and can not explain the measured electron current. Different theories have been proposed to explain this anomalous electron diffusion. It has often been attributed to electron-wall interactions or plasma turbulence, or a combination between the two theories. [START_REF] Hagelaar | Role of Anomalous Electron Transport in a Stationary Plasma Thruster Simulation[END_REF][START_REF] Meezan | Anomalous Coaxial Electron Mobility in a Coaxial Hall Discharge Plasma[END_REF][START_REF] Janes | Anomalous Electron Diffusion and Ion Acceleration in a Low-Density Plasma[END_REF][START_REF] Barral | Wall Material Effects in Stationary Plasma Thrusters. II. Near-Wall and In-Wall Conductivity[END_REF][START_REF] Koo | Modeling of Anomalous Electron Mobility in Hall Thrusters[END_REF][START_REF] Scharfe | Comparison of Hybrid Hall Thruster Model to Experimental Measurements[END_REF][START_REF] Parra | A Two-Dimensional Hybrid Model of the Hall Thruster Discharge[END_REF] Fully kinetic Particle-In-Cell (PIC) simulations in the axial and azimuthal directions have demonstrated the existence of a turbulent azimuthal electric field that leads to a cross-magnetic field transport of the electrons [START_REF] Adam | Study of Stationary Plasma Thrusters using Two-Dimensional Fully Kinetic Simulations[END_REF][START_REF] Adam | Physics, Simulation and Diagnostics in Hall Thrusters[END_REF] . Recently, a collective light scattering diagnostic has evidenced azimuthal fluctuations of the electron density in the same wavelength range as in PIC simulations. [START_REF] Tsikita | Dispersion Relations of Electron Density Fluctuations in a Hall Thruster plasma, Observed by Collective Light Scattering[END_REF] However, a scaling law of the anomalous transport has not yet been proposed. Consequently, most of hybrid models have adopted elementary formulas with tuned parameters to account for the anomalous electron transport. [START_REF] Hagelaar | Role of Anomalous Electron Transport in a Stationary Plasma Thruster Simulation[END_REF][START_REF] Meezan | Anomalous Coaxial Electron Mobility in a Coaxial Hall Discharge Plasma[END_REF][START_REF] Janes | Anomalous Electron Diffusion and Ion Acceleration in a Low-Density Plasma[END_REF][START_REF] Boniface | Anomalous Cross Field Electron Transport in a Hall Effect Thruster[END_REF][START_REF] Bareilles | Critical Assessment of a Two-Dimensional Hybrid Hall Thruster Model : Comparisons with Experiments[END_REF] Over the past 10-years, the time-averaged LIF technique has been extensively used and many measurements of the IVDF have been performed along the thruster channel centerline both inside and outside the thruster over a broad range of electrical power. [START_REF] Hargus | Laser-Induced Fluorescence Measurements of Velocity Within a Hall Discharge[END_REF][START_REF] Dorval | Determination of the Ionization and Acceleration Zones in a Stationary Plasma Thruster by Optical Spectroscopy Study : Experiments and Model[END_REF][START_REF] Hargus | Near Exit Plane Velocity Field of a 200-Watt Hall Thruster[END_REF][START_REF] Hargus | Ion Velocity Measurements Within the Acceleration Channel of a Low-Power Hall Thruster[END_REF][START_REF] Gawron | Influence of Magnetic Field and Discharge Voltage on the Acceleration Layer Features in a Hall Effect Thruster[END_REF][START_REF] Mazouffre | Xe + Ion Transport in the Crossed-Field Discharge of a 5-kW-Class Hall Effect Thruster[END_REF] Recent improvements of the signal-to-noise ratio have been achieved. [START_REF] Mazouffre | Ion Diagnostics of a Discharge in Crossed Electric and Magnetic Fields For Electric Propulsion[END_REF] In parallel, a technique based on the calculations of the moments of the Boltzmann equation has been recently proposed to extract the on-axis distribution of the electric field from well-defined LIF measurements. 23 The purpose The 31st International Electric Propulsion Conference, University of Michigan, USA September 20 -24, 2009 of this work is to combine modelling and experimental approaches to extract useful informations to determine the axial electron mobility profile.

Results show that we are able to propose an analytical constant in time fit of the electron mobility that can reproduces the experimental features such that the electric field and ion distribution functions for a PPS®X000 for discharge voltages varying between 500 and 700 V. The paper is organized as the following, we come back to the previous method used to determine the axial electron mobility in section II. We briefly describe the PPS®X000 thruster in section III. The description of the hybrid model is presented in section III, while the main results are summarized in section IV. We finally end this paper with a conclusion in section V.

II. Direct measurement of the cross-magnetic field electron mobility

The electron mobility perpendicular to the magnetic field µ e,⊥ under the drift-diffusion approximation used in hybrid model is defined as : where e is the elementary electron charge, and n e , T e , and j e,⊥ , are the electron density, temperature, and cross-field current density, respectively. ∇ ⊥ symbolizes the cross-field gradient, and E ⊥ is the electric field perpendicular to the magnetic field.

The elementary theory links the cross magnetic field electron mobility to the total electron momentum transfer collision frequency ν m and to the cyclotron angular frequency ω :

= + = ⊥ ω ν ω ν µ ; ) ( 1 1 2 , (2),
where m e is the electron mass, and B is the magnetic field magnitude. Because ω >> ν m , eq. ( 2) can be rearranged as:

B e Ω = ⊥ 1 , µ (3),
where Ω is called the Hall parameter :

m ν ω = Ω (4).
We can also rearrange eq. ( 1), the electron current being calculated from the current conservation current :

(5).

⊥ ⊥ - = , , / i i d e V en A I j
In eq. ( 5), I d is the discharge current, A is the channel cross section area, n i and V i,⊥ are respectively the density and the mean velocity perpendicular to the magnetic field of ions.

Combining eqs. ( 1), ( 3) and ( 5), assuming a quasineutral plasma and neglecting the electron diffusion term, we can link the Hall parameter Ω (or electron cross field mobility ) to I d , B, E ⊥ n e and V i,⊥ :

⊥ ⊥ - = Ω E en V en A I B e i e d
, / 1 (6).

The 31st International Electric Propulsion Conference, University of Michigan, USA September 20 -24, 2009

Meezan et al. [START_REF] Meezan | Anomalous Coaxial Electron Mobility in a Coaxial Hall Discharge Plasma[END_REF] have combined LIF diagnostic to determine the ion velocity profile and electrostatic probes to measure the electron density and the electric field profiles in order to obtain the Hall parameter variation. Note that the measurements have been achieved along the thruster centerline ; the perpendicular to the magnetic field direction in eq. ( 5) reduces to the axial direction in the region of strong magnetic field magnitude where the magnetic lens is convex.

We report in Figure 1 the measurements of Meezan et al. [START_REF] Meezan | Anomalous Coaxial Electron Mobility in a Coaxial Hall Discharge Plasma[END_REF] for Hall thruster with a channel of 90 mm external diameter and 80 mm length. The xenon mass flow is 2.3 mg.s -1 and the discharge voltage V d varies from 100 to 200 V. The "experimental" inverse Hall parameter profile has been derived from the plasma measured data, while the "classical" inverse Hall parameter profile assumes that only electron neural collisions participate to the momentum transfer [the Hall parameter being derived from eq. ( 4)]. In the estimation of the electron-neutral collision frequency ν e-n , the neutral density profile is deduced from LIF measurement of the neutral atom velocity (combined with the LIF ion velocity data) and the electron-neutral collision momentum rate is calculated assuming a Maxwellian distribution function for an electron temperature measured with probes.

The measurements presented in Figure 1 demonstrate without ambiguity the role of anomalous transport mechanism in the region of strong magnetic field and low neutral density. The authors correlate their measurements with the Bohm formulation of the anomalous transport [START_REF] Bohm | The Characteristics of Electrical Discharges in Magnetic Fields[END_REF] (τ being mean time between collisions). The large error bars observable in Figure 1 can be certainly a consequence of different factors : the number of data required, the use of intrusive probes perturbating the thruster operation, the difficulty to interpret the measurements in the region of strong magnetic field, etc. The accuracy of the measurements is insufficient to extract a quantitative axial profile of the Hall parameter and consequently of the cross field electron mobility that can be useful to simulate the thruster working. HETs used nowadays on board telecommunication satellites are in the range of 1-2 kW power. In the same range of power, the Smart-1 lunar mission of the European Space Agency has demonstrated that HETs can be employed for long duration mission. [START_REF] Koppel | The Smart-1 Electric Propulsion Subsystem in Flight Experience[END_REF] During the operation, the PPS®1350 manufactured by Snecma has covered 100 millions of km during 16 months with 82 kg of xenon for a thrust of 70 mN. The necessary increase of the mass on board telecommunications satellites (e.g. Alphabus platform [START_REF] Roux | Alphabus, the European Platform for Large Communication Satellites[END_REF] ), enabled by the increase of the payload carried by the launchers, requires high power HETs. Several 5 to 10 kW-class HETs have been manufactured and life-tested worldwide to be employed for North/South station keeping for such heavy-mass satellites (namely the NASA-173M in United States [START_REF] Hofer | High-Specific Impulse Hall Thrusters, Part 1 : Influence of Current Density and Magnetic Field[END_REF] , the SPT-140 in Russia [START_REF] Lazurenko | Dual-Mode Operation of Stationary Plasma Thrusters[END_REF] , and the PPS®X000 in Western Europe [START_REF] Duchemin | The Smart-1 Electric Propulsion Subsystem in Flight Experience[END_REF] ). New types of mission like orbit transfer and manoeuvres of the satellites for high power HETs are expected too. We present in Figure 2 a photo and a view of the plume of the PPS®X000. To test the PPS®X000, the pumping system in the PIVOINE facility (4 m long and 2.2 in diameter) in Orleans has been upgraded to reach a pumping speed of 150 000 l/s in xenon. [START_REF] Mazouffre | Expanding Frontiers : Towards High Power Hall Effect Thrusters for Interplanetary Journeys[END_REF] The PPS®X000 has been intensively tested in order to examine the possible dualoperation of the thruster (high specific impulse -low thrust for station keeping and low specific impulse -large thrust for orbit transfer). [START_REF] Duchemin | Stretching the Operational Envelope of the PPS®X000 Plasma Thruster[END_REF] During the study, the range of power was varying between 1 to 7 kW and the magnetic field was kept constant ; the mass flow rates varies from 5 to 15 mg.s -1 and the discharge voltage ranges from 200 to 1000 V. At 5 kW, the maximums of specific impulse or thrust were 2800 s or 260 mN, for a thruster efficiency around 0.5. The main conclusion of this study was that the PPS®X000 was able to operate at high thrust level, as desired. [START_REF] Duchemin | The Smart-1 Electric Propulsion Subsystem in Flight Experience[END_REF] A laboratory version of the PPS®X000, namely the PPS®X000-LM, has been manufactured in order to characterize the discharge properties through LIF measurements. [START_REF] Mazouffre | Xe + Ion Transport in the Crossed-Field Discharge of a 5-kW-Class Hall Effect Thruster[END_REF][START_REF] Mazouffre | Ion Diagnostics of a Discharge in Crossed Electric and Magnetic Fields For Electric Propulsion[END_REF] A hole has been bored through the center of the anode plane so that a laser beam can pass in the direction parallel to the thruster channel. The external ceramic has a slit from 15 mm to the exit plane in order to collect the fluorescence signal inside the thruster channel. In the rest of the paper the exit plane is positioned at x = 0.

IV. Description of the hybrid model

The hybrid model is a transient two-dimensional code. The computational domain takes into account the channel and the near-outside region ; it extends on 10 cm in the axial direction x and 12 cm in the radial direction r. The cathode is positioned in the exit plane at r = 11 cm. The number of grid cells is 50 and 60 in the x and r directions, respectively. We begin this section by shortly describing in the sub-section A. the heavy species transport, in the sub-section B. we briefly come back to the treatment of the electron transport. The ad hoc hypothesis up to now used in the description of the cross magnetic field mobility of the electrons is detailed in sub-section C.

A. Ion and neutral transport

Xe singly-charged ions are generated by the ionization of the neutral atoms. The ionization rate is tabulated as a function of the electron mean energy assuming a Maxwellian electron distribution. We use the macro-particle technique as in Particle-In-Cell model, each macro-ion represents a certain number of individual ions (30 000 macro-ions are used) For simplicity reasons, we use the term ion for macro-ion in the rest of the section. The ion trajectories are integrated, the Newton Lorentz force reduces to the electric field component, ions being unaffected by the magnetic field due to their large mass and high velocity. Ions are followed until their leave the computational domain or impact on walls. Ions impinging walls are neutralized and a neutral atom is re-emitted isotropically assuming a half-Maxwellian distribution function with temperature of 500 K. A certain number of diagnostics such as the velocity distribution function has been also incorporated to compare the calculated distributions with the measured ones.

We use for xenon atoms the same technique as the ions using macro-atoms that account for a large number of individual atoms (typically 20 000 macro-atoms are simulated). The xenon neutral atoms are injected through a hole

The 31st International Electric Propulsion Conference, University of Michigan, USA September 20 -24, 2009 in the anode plane assuming a Maxwellian distribution of the flux in the axial direction for a given temperature equal to 500 K. We also account for the background neutral atoms due to the backpressure tank by injecting an additional neutral flux at the thermal velocity (temperature of 300 K) from the boundaries of the computational domain in the near field region. Neutrals colliding with the walls are isotropically reflected back in the domain.

Neutrals that exit the computational domain are not longer followed.

B. Electron transport

The electron fluid transport is described with the continuity, momentum and energy equations. The equation for the electron conservation is given by :

S t n e e = ⋅ +∇ ∂ ∂ Γ (7),
where Γ e is the electron flux, and S the ionisation source term.

The electron momentum equation is written -in the form of drift-diffusion approximation -as :

( ) (
)

e e e e T n V n ∇ - ∇ =µ Γ (8).
Due to the magnetic field, the mobility µ is not a scalar but a tensor [START_REF] Lieberman | Principles of Plasma Discharges and Materials Processing[END_REF] , where the component parallel to the magnetic field is given by the following relation : The component perpendicular to the magnetic field is given by the eq. ( 2). Combining eqs. ( 7) and ( 8), and assuming quasineutrality ( In eq. ( 10), the plasma density n as well as the ion flux Γ i are deduced from the ion transport. The numerical method used to solve eq. ( 10) is presented by Hagelaar. [START_REF] Hagelaar | Modelling Electron Transport in Magnetized Low Temperature Discharge Plasmas[END_REF] The Boundary conditions are fixed potential at the cathode V c =0 and at the anode V a =V d -where V d is the discharge voltage. The electric potential is fixed to 0 on the open frontier of the computational domain. We impose on the dielectric walls, on the front plane an equality between the electron flux Γ e and the ion flux Γ i .

In order to determine the electron temperature T e (or the electron mean energy ε=3/2βT e -β is the Boltzmann factor) to calculate the ionization rate, we solve an energy equation : (11),

nW Nn V e n e t n e e - - ∇ = ∇ ⋅ ∇ - ⋅ ∇ + ∂ ∂ κ ε ε µ ε ε . ) ( 9 10 ) ( 3
5 ) ( Γ Γ
N being the neutral gas density. The last two terms in the energy equation represent energy loss by collisions with gas particles and with the walls, respectively, where κ and W are effective energy loss coefficients [START_REF] Hagelaar | Two Dimensional Model of a Stationary Plasma Thruster[END_REF] .

C. Ad hoc formulation of the cross field electron mobility

Up to now, in the hybrid model, the total electron momentum transfer frequency ν m is an effective collision frequency including the electron-atom collisions (frequency ν e-n ), electron-ion collisions (frequency ν e-i ), electronwall collisions (frequency ν e-w ), and Bohm-like transport (frequency ν Bohm ). The total electron momentum transfer frequency in eq. ( 2) is written as the sum of each contribution, therefore, ν m = ν e-n + ν e-i + ν e-w + ν Bohm.

4-6, 14-15 The momentum transfer frequency collision between electrons and atoms ν e-n is inferred from the neutral density profile and the electron-atom momentum transfer rate that is assumed constant (k m = 2.5×10 -13 m 3 s -1 ). electron-ion collisions on the momentum transfer is deduced from the Coulomb cross sections assuming a Maxwellian electron distribution function. [START_REF] Cherrington | Gaseous Electronics and Gas Lasers[END_REF] We account for wall effects inside the thruster channel with a collision frequency equal to ν e-w = αν ref (with ν ref = 10 7 s -1 ). Notice that some authors use refine theories to account for secondary electron emission under high energetic electron impacts on the ceramic walls. [START_REF] Barral | Wall Material Effects in Stationary Plasma Thrusters. II. Near-Wall and In-Wall Conductivity[END_REF][START_REF] Scharfe | Comparison of Hybrid Hall Thruster Model to Experimental Measurements[END_REF][START_REF] Parra | A Two-Dimensional Hybrid Model of the Hall Thruster Discharge[END_REF] We assume that the anomalous transport is governed outside the channel by Bohm-like transport with an equivalent collision frequency ν Bohm = KΩ/16. The two parameters α and K are adjusted in order to match experimental integrated measurements (performance, ionization efficiency, current and its time variation, etc.). [START_REF] Bareilles | Critical Assessment of a Two-Dimensional Hybrid Hall Thruster Model : Comparisons with Experiments[END_REF] We can remark that electron-wall collisions does not depend of the magnetic field, consequently, the mobility varies as 1/B 2 like the "classical" collision mobility, while Bohm-like transport, leads to a mobility varying as 1/B.

Scharfe et al. [START_REF] Scharfe | Comparison of Hybrid Hall Thruster Model to Experimental Measurements[END_REF] have compared the plasma modeling results when the Bohm or the experimental (discussed in section II -see Fig. 5) mobility is assumed, the authors conclude that the experimental mobility leads to plasma properties in better agreement with experimental results obtained for the Stanford Hall thruster. Koo and Boyd 8 have also compared the ad hoc Bohm-type assumed anomalous transport with a deduced mobility profile from measurements (combined with some modeling assumptions due to the lack of the electron current density measurements). Koo and Boyd find a better agreement of the overall behavior of the UM/AFRL P5 Hall thruster of 5 kW is observed when the experimental mobility profile is used. Following the same idea, the goal of the paper is to compare the ad hoc hypothesis of the anomalous transport with the reconstructed mobility profile obtained from LIF measurements for the PPS®X000-LM.

V. Results and discussion

In sub-section A, we test the model ad hoc hypothesis on the electron mobility comparing the experimental and calculated ion velocity profiles for the PPS®X000-LM. In sub-section B, we examine the influence of the discharge voltage on the ion velocity distribution functions. In all this study the coil currents (17 A) and the xenon mass flow (6 mg.s -1 ) remain the same.

A. Mobility profile

In Figure 3a, the experimental as well as the simulated time-averaged on-axis profiles of the axial electric field are shown together with the computed ionization source term for a discharge voltage of 500 V. In the calculations indicated as hybrid # 1, we represent the anomalous electron transport mechanisms taking into account wall effects inside the channel with a collision frequency equal to ν e-w = αν ref (with ν ref = 10 7 s -1 ), and outside the channel Bohmlike transport with an equivalent collision frequency ν Bohm = KΩ/16. The adjustable coefficients are α = 1.45 and K = 0.2, as in the work of Boniface et al. [START_REF] Boniface | Anomalous Cross Field Electron Transport in a Hall Effect Thruster[END_REF] The electron mobility profile calculated with eq. ( 1) is plotted in Fig. 3b. Clearly there is no agreement with experiments (see Fig. 3b). The use of empirical laws with 1/B 2 (or 1/B -not shown here) mobility inside the channel and 1/B outside the channel as we have previously done with tunable coefficients do not allow us to reproduce the measured profile of axial ion velocity.

The reason is clearly associated to the misestimating of the electric field profile. The method employed to extract electric field profile from measurements is detailed in the paper of Pérez-Luna et al. 23 In the measurements, the maximum of electric field reaches 35 kV.m -1 , while the calculations give a maximum of electric field in order of only 20 kV.m -1 . We also see that the length of the acceleration layer is typically few millimeters in the experiments, while it extends on more than one centimeter in the calculations. Consequently, the experimental ion velocity profile differs from the calculated one.

In the calculations labeled as hybrid # 2, we have adjusted the anomalous electron mobility profile in order to match the experimental profile of the axial ion velocity (see Fig. 3b). The shape of calculated axial electric profile presented in Fig. 3a is now in agreement with the profile deduced from LIF measurements. Integrating the measured axial electric field profile leads to a potential drop much smaller than the discharge voltage (sheath falls). This explains the shift between the simulated and measured most probable velocity (see Fig. 4). The fitted electron mobility profile (depending on the axial position) shown in Fig. 3b is in qualitative agreement with the measured profile in the P5 thruster [START_REF] Koo | Modeling of Anomalous Electron Mobility in Hall Thrusters[END_REF] . Empirical laws with 1/B 2 or 1/B mobility (see Fig. 3b -hybrid # 1) differ by far from the analytical profile.

B. Ion velocity distribution functions

We present in the Figure 4 the experimental and calculated IVDFs at the end of the acceleration layer at x = 4 cm (the reference axial position x = 0 is at the exit plane). Part of the broadening of the IVDFs is associated to the overlap between the ionization and acceleration layers that is visible in the Figure 3a. We have already noticed that transient evolution of the plasma also influences the shape of the IVDFs. The choice of tuned coefficients in the elementary laws (1/B 2 or 1/B) used so far in hybrid models has strongly influenced the dynamic behavior of the discharge, especially the magnitude of transit-time oscillations. The transit-time oscillations (100-500 kHz) are associated to the time needed by the ions to cross the acceleration layer. The visible consequence was noticed in the shape of IVDFs at the end of the acceleration layer where a large velocity dispersion was obtained [START_REF] Hagelaar | Role of Anomalous Electron Transport in a Stationary Plasma Thruster Simulation[END_REF][START_REF] Bareilles | Critical Assessment of a Two-Dimensional Hybrid Hall Thruster Model : Comparisons with Experiments[END_REF] .

In Figure 4a, we show the calculated and the measured IVDFs of the PPS®X000-LM thruster for discharge voltages in the range of 500-700 V. The empirical mobility profile is the same as in Fig. 3b. The calculated IVDFs now match the measured IVDFs (especially the broadening of the distribution). This suggests that the highfrequency dynamic behavior is correctly represented in the model. The theoretical maximum axial velocity an ion can reach for a given discharge voltage V d is v max = (2eV d /m Xe+ ) 1/2 where m Xe+ is the ion mass. Both model and experiments confirm the existence of very fast ions with a velocity larger than v max in the final part of the acceleration layer. This is clearly attributed to a "wave riding" phenomenon that makes possible for the ions to acquire a kinetic energy larger than the given potential energy e×V d . [START_REF] Hagelaar | Role of Anomalous Electron Transport in a Stationary Plasma Thruster Simulation[END_REF][START_REF] Bareilles | Critical Assessment of a Two-Dimensional Hybrid Hall Thruster Model : Comparisons with Experiments[END_REF][START_REF] Gawron | Influence of Magnetic Field and Discharge Voltage on the Acceleration Layer Features in a Hall Effect Thruster[END_REF][START_REF] Mazouffre | Ion Diagnostics of a Discharge in Crossed Electric and Magnetic Fields For Electric Propulsion[END_REF] Conserving a given profile of the cross field electron mobility leads also to reasonable describe the thruster working.

The 31st International Electric Propulsion

VI. Conclusions and future work

In this paper, we have demonstrated that one single diagnostic such that laser spectroscopy technique is able to provide indirect informations on the axial electron mobility profile. Results show that previous elementary laws (1/B 2 or 1/B) with adjusted coefficients in order to match the correct "macroscopic" properties of the thruster (current, power, performance) used so far in hybrid models fail to reproduce experimental observations such that as the measured ion velocity profile. Fitting the calculated ion velocity profiles with the LIF measurements obtained at 500 V leads to a empirical steady spatial profile of the axial electron mobility that is able to reproduce the thruster properties and the measured ion velocity distribution functions, in the wide range of discharge voltages for the PPS®X000-LM thruster. The hybrid model is still unfortunately non-self consistent.

The adjusted electron mobility profile in this study is in qualitative agreement with the calculated electron mobility from PIC simulations [START_REF] Adam | Physics, Simulation and Diagnostics in Hall Thrusters[END_REF] . We now plan to use the PIC model in order to identify and understand the physical The 31st International Electric Propulsion Conference, University of Michigan, USA September 20 -24, 2009

The 31st International Electric Propulsion Conference, University of Michigan, USA September 20 -24, 2009 10 parameters that govern the anomalous electron transport. We also want to deduce from these calculations scaling laws of the electron mobility as a function of the discharge parameters in order to propose a less empirical description of the electron mobility in the hybrid model.

  Hall effect thrusters is still an open question. Models have so far assumed 1/B 2 or 1/B scaling laws for the "anomalous" electron mobility, adjusted to reproduce the integrated performance parameters of the thruster. We show that models based on such mobility laws predict very different ion velocity distribution functions (IVDF) than measured by induced fluorescence (LIF). A fixed spatial mobility profile, obtained by analysis of improved LIF measurements, leads to much better model predictions of thruster performance and IVDF than 1/B 2 or 1/B mobility laws, for discharge voltages in the 500-700 V range for a 5 kW Hall thruster. Nomenclature A = channel cross section area B = magnetic field magnitude e = elementary electron charge E ⊥ = component of the electric field perpendicular to the magnetic field I d = discharge current j e,⊥ = cross magnetic field electron current density K = Bohm fitting parameter k m = electron-atom momentum transfer rate m a = xenon anode mass flow rate m e , m Xe+ = electron, ion mass N, n, n e , n i = neutral, plasma, electron and ion density r = radial direction S = ionization source term T e = electron temperature V a , V c , V = anode, cathode and electric potential V d = discharge voltage V i,⊥ = ion mean velocity perpendicular to the magnetic field v max = theoretical maximum axial ion velocity W = electron-wall effective energy loss coefficient x = axial direction Xe + = singly charged ion α = adjustable coefficient for electron-wall collision frequency β = Boltzmann factor ε = electron mean energy Γ i , Γ e = ion and electron vector flux κ = inelastic effective energy loss coefficient µ = electron mobility tensor µ e,// , µ e,⊥ = electron mobility parallel and perpendicular to the magnetic field ν e-i , ν e-n = electron-ion, electron-neutral collision frequency ν e-w = electron-wall collision frequency ν m = total electron momentum transfer collision frequency τ = mean time between collisions ω = cyclotron angular frequency Ω = Hall parameter ∇ = gradient ∇ ⊥ = cross magnetic field gradient

Figure 1 :

 1 Figure 1 : Axial profile of the inverse Hall parameter.[START_REF] Meezan | Anomalous Coaxial Electron Mobility in a Coaxial Hall Discharge Plasma[END_REF]

Figure 2 :

 2 Figure 2 : (a) photo of the PPS®X000, (b) plume of the PPS®X000.

  are not interested in the description of cathode and anode sheaths, the electric potential V is calculated solving the following elliptic equation:

1 )Figure 3 :

 13 Figure 3 : Time-averaged profiles along the PPS®X000-LM thruster channel axis of (a) the electric field deduced from LIF measurements and the computed electric field profile and the calculated ionization source term, (b) the electron mobility perpendicular to the magnetic field and most probable ion velocity. Conditions are : V d = 500 V, m a = 6 mg.s -1 . The coil currents is fixed to 17 A. The channel exit is at x = 0.

Figure 4 :

 4 Figure 4 : IVDFs obtained by LIF measurements 22 and calculated with the PPS®X000-LM at x = 4 cm for (a) V d = 500 V, (b) V d = 600 V, and (c) V d = 700 V, and for m a = 6 mg.s -1 . The coil currents is 17 A. The theoretical maximum axial velocity v max is also shown.

Acknowledgments

This work is performed in the frame of the joint-research program CNRS/CNES/Sncema/Universités n° 3161 "Propulsion par Plasma dans l'Espace" and by the TELIOPEH ("Transport ELectronique et IOnique dans les Propulseurs à Effet Hall") project, funded by ANR (National Research Agency), under contract ANR-06-BLAN-0171. J. Pérez-Luna has benefited from a CNES/Snecma PhD fellowship.


L V Simonchik

Th Callegari

J Lo

J Sokoloff

B I Stepanov

DETERMINATION OF ELECTRON DENSITY IN PLASMA USED FOR THE WAVEGUIDE BAND-PASS FILTER CONTROL

Introduction

Most of electromagnetic band gap (EBG) devices are realized by artificial metallic structures such as using wires to produce effective negative permittivity. In some number of papers it was shown that the glow discharge plasmas have high potentials as effective control elements in the EBG devices /1, 2/. Band-pass waveguide microwave filter can be considered as 1D EBGstructure /3/. Microwave filter constitutes waveguide section with diaphragms (metallic rods), which form the resonators with direct (or half wavelength) couplings. The frequency characteristic tuning of such waveguide microwave filter can be carried out by changing of resonators properties, for example, using plasma columns instead some metallic rods. In this case it turns out that the filter frequency characteristic is connected with plasma parameters and the possibility to determine the electron density appears. This possibility we demonstrate in this work.

In general, the permittivity of a plasma ε pe for microwave at frequency ω propagated in plasma is given by

) / 1 ( 1 ) / 1 ( 1 2 0 2 2 2 ω ν ω ε ω ν ω ω ε i m n e i e e pe pe + - = + - = (1)
where ν is the electron collision frequency, e is the electron charge, ε 0 is the permittivity in vacuum and m e is the electron mass. From this relation it is seen that ε pe can be modified from unity to negative values according to n e .

In the given work we have used the band-pass waveguide filter of the 3centimeter wavelength ranges as a 1D periodical structure with gas discharge lamps GSh-5 for its control /4/.

Results and discussions

The 3-centimeter wavelength filter is formed by periodically distributed inhomogeneities in waveguide (with cross-section 23×10 mm 2 and length of 200 mm). In our case, the five open-end holes (diameter 4 mm) were drilling through the opposite wide waveguide walls with period l along waveguide. Different inhomogeneities (metallic or dielectric rods) can be installed in holes (Fig. 1).

The period of the location of rods (length resonators) is about l = Λ 0 /2, where Λ 0 is a wavelength in the waveguide, which is expressed by relation

2 / 1 2 0 ) ) 2 / ( 1 ( - - = Λ a m m λ λ , (2)
where λ m is a wavelength corresponding to the middle of the filter bandwidth; a is the wide waveguide wall dimension (a = 23 mm). Thus, a character of resonance is half-wave and in case of five installed metallic rods (Fig. 1,b) we have four resonators with direct coupling. For the central frequency f 0 = 10.2 GHz (λ m = 29 mm) the resonator lengths is l = 19 mm.

The amplitude-frequency characteristics (transmission spectrums) of the 3centimeter microwave filter was measured using the vector analyzer Anritsu 37369C in the frequency range of 7-15 GHz (Fig. 2). The transmission characteristic of the filter with five metallic rods have one pass band at the frequency of about f 0 = 10.2 GHz. The transmission level for this band of about 10-15 dB lower than for the empty waveguide section. It can be explained by rough-and-ready manufacturing of this filter. Outside the bandwidth the transmission signal is less than -60 dB. For filter with three metallic rods installed in odd holes only (two resonators with l = 38 mm) the amplitude-frequency characteristic has three pass bands with the central frequencies f 01 = 7650 MHz, f 02 = 10240 MHz and f 03 = 13520 MHz (fig. 2).

In the experiments, the lamps GSh-5 was used as a control elements of microwave filters /4/. The total length of the lamp glass tube is 250 mm, and its internal diameter is about 3 mm. The GSh-5 lamp is filled by neon at a pressure of 70 Torr, that gives the neon atoms concentration of about N = 2.5×10 18 cm -3 . If the diffusion current component is negligibly small compared with the drift and if the ion current is negligible compared with the electron one, the electron density n e in the positive column is connected with values of current density j and the drift velocity ν dr (E/N) by following relation: where E/N is reduced electric field in the positive column, e is electron charge. Current density was defined as the ratio of discharge current to the cross-section of the current area flow (internal cross-section area of GSh-5). Since the filamentary oxide cathode is used in this lamp the cathode drop voltage is small (a few tens of Volts) and electric field strength in positive column was defined as a ratio of the voltage between electrodes (about 200 V) and the interelectrode distance (200 mm). The dependence of the drift velocity of electrons from the reduced electric field in the plasma was calculated using BOLSIG + /5/. Electron densities at different currents defined by the formula (3) are presented in Fig. 3. Let's consider the case when the filter is formed by combination of three metal rods and two lamps GSh-5 (Fig. 4,a). If the lamps are turn off, the waveguide section has two pass band in analyzer swath of 7-12 GHz -more narrow at 7.6 GHz and wider at 10.4 GHz (Fig. 4,b). With increasing current of lamps from 0 to 120 mA, a decrease of the transmission level of the bandwidth at 7.6 GHz from -10 to -50 dB was observed (Fig. 4,b) and it shifts toward high frequency.

( / ) e dr j en v E N = , (3)
The results of calculation of the filter transmission in this case for different values of electron density are presented in fig. 4,c. Calculations were carried out by the code HFSS Ansoft where n e and ν were used as variable input parameters. In this calculation it was supposed that ν = 2.8×10 10 s -1 . As it can be seen, a total suppression and shift of the transmission peak at 7.6 GHz is observed even at the electron density of 5×10 12 cm -3 . This is in accordance with electron density at discharge current of ~ 30 mA (n e = 6.6×10 12 cm -3 ) (fig. 3). It should be noted that position of the shifted peak practically does not depend on the collision frequency. Let's use this fact in order to determine electron density in lamps.

First of all we should calibrate a pass band shift in accordance with change of plasma permittivity. For calibration we used the two alumina rods (4 mm in diameter) with known permittivity ε = 9 instead discharge lamp. Corresponding amplitude-frequency characteristics are shown in fig. 2. The increase of changed filter elements permittivity from 1 to 9 (Δε =8) leads to shifts of both the left and right pass bands toward the low frequency at the saving of the central pass band position. Unfortunately, the shift of the left pass band is so large that it is outside of analyzer swath. However we can use the shift of the right pass band (Δf ~ 2200 MHz) for calibration. It corresponds to Δε/Δf ~ 0.0036 MHz -1 .

As it can be seen in fig 4, b, at discharge current of 30 mA the left pass band frequency shift Δf pe is about 200 MHz and at 85 mA -900 MHz. Then the plasma permittivity change Δε pe can be defined as the following:

Δε pe = (Δε/Δf)×Δf pe ×(V d /V pe ), ( 5) where V d /V pe is a ratio of the volumes of dielectric and plasma column in waveguide. Thus we obtain Δε pe (30mA) = -0.03 and Δε pe (85mA) = -4.5. It gives according to (1) and without consideration of the collision frequency the electron density of 7.2×10 11 cm -3 and 3.8×10 12 cm -3 at 30 and 85 mA, correspondingly. If we take into account the collision frequency (2.8×10 10 s -1 ) according to the HFSS calculation, the electron density magnitudes will be n e (30 mA) = 1.1×10 13 cm -3 and n e (85 mA) = 5.7×10 13 cm -3 .

Thus, the electron density determined using resonance properties of direct coupling waveguide microwave filter of 3-cantimiter wavelength range (EBGstructure) is in good agreement with estimation obtained by using reduced electric field in positive column of discharge lamp. Developed method will be applied in future at the usage of the atmospheric pressure glow discharges as controlling elements of waveguide filters.

Financial support of BRFBR grants (F11F-002, F11M-044) is acknowledged.

Fig. 1 .

 1 Fig. 1. Photos of the waveguide section with (a) hole, (b) 5 and (c) 3 metallic rods.

Fig. 2 .

 2 Fig. 2. Transmission of waveguide section without metallic rods (1), with 5 m. rods (2), with 3 m. rods (3), with 3 metallic and 2 dielectric rods (4).

Fig. 3 .

 3 Fig. 3. Electron density in lamp against discharge current.

Fig. 4 .

 4 Fig. 4. Photo of waveguide section with 2 lamps and 3 rods (а) and corresponding experimental (b) and calculated (c) transmittance spectra.

-0 mA;

-6.5 mA;

-29 mA;

-85mA. 1' -0 cm -3 ; 2' -1*10 12 cm -3 ; 3' -5*10 12 cm -3 ; 4' -1*10 13 cm -3 ν = 2.8×10 10 s -1 ;

5' -1*10 13 cm -3 ν = 0 s -1 .


Marwan Hamze

Mehdi Benallegue

Rafael Cisneros

Abdelaziz Benallegue

A Balance-Force Controller for a Legged Robot with Multiple Flexible-Contacts

In this paper, we propose a controller for multicontact legged robots that takes into account the contact compliance. This controller is able to achieve both balance stabilization and force control in the same loop in task space thanks to the use of an explicit contact flexibility model and simplified centroidal dynamics that allow exploiting the redundancy of the kinematic and force feedback. The control problem is formulated as an LQR based on a linearized model of the reduced non-linear model. The performance of the controller with regard to robustness to modeling error and external perturbations has been tested in simulation, and compared to a stabilizer based on the zero-moment-point feedback within a control that generates the desired contact forces with the environment.

I. INTRODUCTION

Legged robots are considered under-actuated because they lack direct actuation to generate translations and rotations of their body in their environment. A robot's unactuated dynamics and its balance in particular depend only on external forces, mostly reduced to contact forces and gravity [START_REF] Wieber | On the stability of walking systems[END_REF]. In particular, the bipedal walking of humanoid robots requires fine control of the contact forces to ensure stability despite the small contact surface. In fact, it is common to exploit additional contact positions allowing the robot to achieve complex transportation motions with increased stability [START_REF] Escande | Contact planning for acyclic motion with tasks constraints[END_REF], [START_REF] Saab | Dynamic Whole-Body Motion Generation Under Rigid Contacts and Other Unilateral Constraints[END_REF], [START_REF] Kudruss | Optimal control for whole-body motion generation using center-of-mass dynamics for predefined multi-contact configurations[END_REF], [START_REF] Vaillant | Multi-contact vertical ladder climbing with an HRP-2 humanoid[END_REF]. However, this usually makes dynamics more complex when the contacts are non-coplanar and of different nature (unilateral, bilateral, point contact, etc.) [START_REF] Perrin | Effective Generation of Dynamically Balanced Locomotion with Multiple Non-coplanar Contacts[END_REF], making it difficult to ensure that all reference contact forces are respected, even when the robot is equipped with force sensors. This is particularly true when dealing with uncertain environments in terms of geometry and compliance.

Indeed, the inaccuracy in force tracking comes from the fact that the contacts that the robot makes with the environment aren't fully rigid. That is not only because the environment is never perfectly stiff, but also because the robot itself has a certain level of compliance. HRP-2 which has flexible rubber bushes, placed at the ankles to absorb foot impacts [START_REF] Kanehira | Design and experiments of advanced leg module (HRP-2L) for humanoid robot (HRP-2) development[END_REF], is a good example of a compliance by design. Likewise, it is common to see robots with a high level of structural compliance, such as COMAN [START_REF] Li | Walking pattern generation for a humanoid robot with compliant joints[END_REF] or Sarcos Primus [START_REF] Stephens | Push Recovery by stepping for humanoid robots with force controlled joints[END_REF]. The presence of such compliance has two consequences. The first one is that forces modify the shape of the environment, and displace the application point of the contact forces, altering the robot's dynamics. The second consequence is that, contrarily to the rigid contact model, the robot is not able to instantly generate the desired force since it depends on its state, which needs then to be modified accordingly. This second consequence is the reason why the dynamics of the Zero Moment Point (ZMP) are often modeled with a first order system [START_REF] Kajita | Biped walking stabilization based on linear inverted pendulum tracking[END_REF] as an approximation, but this model is also inaccurate since the simplest physicsbased model has at least a second order dynamics. Some works explicitly take into account the flexibility model [START_REF] Fahmi | STANCE: Locomotion Adaptation over Soft Terrain[END_REF], [START_REF] Azad | Balance control strategy for legged robots with compliant contacts[END_REF], [START_REF] Vasilopoulos | Monopod hopping on compliant terrains[END_REF], [START_REF] Chang | Learning to jump in granular media: Unifying optimal control synthesis with Gaussian process-based regression[END_REF], but most of them consider only point contacts. Furthermore, most of them lack any force feedback in the control loop, or they just consider the compliance as concentrated on one contact support [START_REF] Benallegue | Humanoid flexibility deformation can be efficiently estimated using only inertial measurement units and contact information[END_REF], [START_REF] Mifsud | Stabilization of a compliant humanoid robot using only Inertial Measurement Units with a viscoelastic reaction mass pendulum model[END_REF].

Classically, to overcome these issues, a second loop of force control is often included, usually as admittance control [START_REF] Stephens | Push Recovery by stepping for humanoid robots with force controlled joints[END_REF], [START_REF] Caron | Biped Stabilization by Linear Feedback of the Variable-Height Inverted Pendulum Model[END_REF], [START_REF] Caron | Stair Climbing Stabilization of the HRP-4 Humanoid Robot using Wholebody Admittance Control[END_REF]. These controllers are usually set up in series with the balance control. This means that the kinematicfeedback-based stabilization sends force references to the force-feedback-based admittance control loop. This is based on the assumption that the dynamics of the force control are fast enough to converge within the time required to maintain balance. However, not only both dynamics remain always coupled, but the admittance control becomes significantly slower with more compliant and uncertain environments. For this reason, there is a need for a controller able to take both dynamics into account in a single loop.

Finally, some parameters in the robot's model or in the environment, such as the stiffness of the contacts, can be very hard to measure correctly. Consequently, we usually resort to simplified models such as a viscoelastic contact pattern, which does not correspond to the physical reality which is nonlinear and even time-varying. Some methods intend to estimate these parameters online [START_REF] Fahmi | STANCE: Locomotion Adaptation over Soft Terrain[END_REF], [START_REF] Alves | A comparative study of the viscoelastic constitutive models for frictionless contact interfaces in solids[END_REF], but this estimation requires to gather dynamic data and may take a lot of time to converge. Therefore we need a controller that is robust to this kind of modeling errors.

We propose in this paper a controller with explicit multicontact support allowing to perform both balance and force control in a single loop. This paper will first introduce the dynamics of the robot and the selected reduced model. Next, it will proceed to describe the synthesis of the proposed controller based on the corresponding linearized model. Finally, it will present some simulations and results.

II. DYNAMICS OF THE ROBOT

In this section we present the Lagrangian dynamics of the robot in contact with a compliant surface. Let us consider a multi-body humanoid robot with n + 6 degrees of freedom (dof), having a configuration described as Ψ = (p B , R, q), where p B ∈ R 3 and R ∈ SO(3) represent the position and orientation matrix of the non-actuated floating base, and q ∈ R n , the joint angles vector. We define α ∈ R n+6 as the robot's configuration velocity vector given by

α = ṗ⊤ B ω ⊤ q⊤ ⊤ , (1)
where ṗB and ω ∈ R 3 are the linear and angular velocities of the base. The angular velocities are such that the derivative with respect to time of the rotation matrix is given by Ṙ = S(ω)R, with S(ω) being the skew-symmetric matrix operator allowing to perform cross-product. The robot's acceleration vector α is the time derivative of α.

The Lagrangian dynamics of the robot are written as

H (Ψ) α + C (Ψ, α) α + G (Ψ) -F = τ, (2)
where n+6) is a matrix accounting for the Coriolis and centrifugal effects, G (Ψ) ∈ R n+6 is the vector of gravitational effects, τ ∈ R n+6 is the input torque vector corresponding to both the under-actuated and actuated dof, and F ∈ R n+6 is the vector of external forces acting through the contacts with the environment, which is calculated using

H (Ψ) ∈ R (n+6)×(n+6) is the robot's inertia matrix, C (Ψ, α) ∈ R (n+6)×(
F = J ⊤ c,1 . . . J ⊤ c,nc        f c,l,1 t c,l,1 . . . f c,l,nc t c,l,nc        , (3)
where J c,i ∈ R 6×(n+6) is the Jacobian matrix of the contact i, and f c,l,i , t c,l,i are respectively the force and torque applied at contact i (i = 1, 2, ..., n c ; n c is the number of contacts), expressed in the local frame of the contact body.

The classic use of this equation is to consider that the environment is rigid. In this case, these forces depend mostly on the torque τ in a way that any force that is feasible (within unilaterality, friction, and torque limit constraints) can be generated instantly. However, a perfectly rigid contact cannot exist, thus, the previous assumption is wrong. The contact forces at a given instant depend on many factors and joint torque is virtually not one of them. In particular, if the forces depend only on the local deformation caused by the interaction with the contact body, we can write the forces as a function of the robot state F (Ψ, α). In this case, we can control the forces indirectly by modifying this deformation through the state of the robot. This requires to explicitly consider the coupled dynamics between the kinematics and the contact forces to ensure the convergence to a given reference.

III. DYNAMIC MODEL FOR THE CONTROLLER

To simply parameterize the contact wrenches and formulate how they affect the position of the Center of Mass (CoM) of the robot, the dynamics are reduced to consider the robot as a single rigid body with viscoelastic contacts. Hereinafter is the description of this model.

A. Viscoelastic Contact Model

A common model for a compliant contact is the viscoelastic approximation allowing to emulate a linear passive interaction. Contacts are defined as the surfaces of the robot's hands and legs applying forces on the environment, each having their own contact frame with a yaw axis that is perpendicular to the environment. Thus, we write

f c,l,i = K f,p,i (p c,i -p c,i,r ) + K f,d,i ṗc,i , (4)
t c,l,i = K t,p,i Θ R c,i,r R ⊤ c,i + K t,d,i ω c,i , (5)
where p c,i ∈ R 3 , ṗc,i ∈ R 3 and ω c,i ∈ R 3 are respectively the position, linear and angular velocities of contact i written in the contact frame and R c,i ∈ SO(3) is the orientation matrix of the contact i. The variables p c,i,r ∈ R 3 and R c,i,r ∈ SO(3) are respectively the rest position and orientation of the contact i, i.e. when forces and torques are zero. The matrices K f,p,i and K f,d,i are the linear stiffness and damping at contact i, while K t,p,i and K t,d,i are the angular stiffness and damping at contact i. The difference in orientation represented by R c,i R ⊤ c,i,r is considered so small that the approximation sinθ ≈ θ is used, thus we define the function Θ as the axis-sine of angle representation. The values of stiffness and damping can be set independently from a contact to another.

B. State and Control Vectors Definition

We define here a reduced state vector focusing on the contacts with the environment. The state variables are the positions, orientations, and velocities of these contacts, in addition to those of the CoM. The contact positions and orientations are going to be defined with respect to the CoM frame, which is the frame centered at the CoM and having the same orientation as the floating base. The state vector is defined by

x ≜ x com x 1 • • • x n ⊤ , (6)
where

x com = ( p ⊤ com ṗ⊤ com Ω ⊤ ω ⊤ ) ⊤ and x i = ( p ⊤ c,i,com Ω ⊤ c,i,b ṗ⊤ c,i,com ω ⊤ c,i,b ) ⊤ .
The vectors p com ∈ R 3 and ṗcom ∈ R 3 are respectively the linear position and velocity of the CoM, Ω and ω are respectively the orientation and the angular velocity of the floating base, all written in the world frame. Similarly,

p c,i,com ∈ R 3 , ṗc,i,com ∈ R 3 , ω c,i,com ∈ R 3
and Ω c,i,com are respectively the linear position, linear velocity, angular velocity and orientation of contact i, all written in the CoM frame. It should be noted that the orientations (Ω and Ω c,i,com ) can be written by using any representation of the orientation, such as the quaternion, the axis-angle, etc. Each Ω has an associated rotation matrix R ∈ SO(3).

As for the control vector, and as we are aiming to control the contacts, the control vector is defined as

u ≜ p⊤ c,1,com ω⊤ c,1,com • • • p⊤ c,nc,com ω⊤ c,nc,com ⊤ , (7)
where pc,i,com ∈ R 3 and ωc,i,com ∈ R 3 are respectively the linear and angular accelerations of contact i, written in the CoM frame. Note that these are immediate time-derivatives of state variables. The dynamics of this state is described hereinafter.

C. Reduced Dynamics

The reduced model is a rigid body with massless legs, which gives a constant inertia tensor I ∈ R 3×3 in the base frame, expressed as RIR ⊤ in the world frame (R here is the orientation matrix of the robot's base). Thus, the angular momentum of the robot is RIR ⊤ ω.

Using Euler's second law, the relation between angular momentum and the total external torque, is expressed by

nc i=1 (S (Rp c,i,com ) f c,i + t c,i ) = S (ω) RIR ⊤ ω + RIR ⊤ ω, (8)
where f c,i and t c,i are respectively the force and torque at contact i written in the world frame, which can be obtained from equations ( 4) and (5) using

f c,i = R c,i f c,l,i , (9)
t c,i = R c,i t c,l,i . (10)
Using Newton's second law, and equation ( 8), the linear and angular accelerations of the CoM and the floating base of the robot can be expressed as

pcom = 1 m nc i=1 f c,i + g, (11
)

ω = RI -1 R ⊤ nc i=1 (S (Rp c,i,com ) f c,i + t c,i ) -RI -1 R ⊤ S (ω) RIR ⊤ ω, (12)
where g = 0 0 -g 0 ⊤ and g 0 is the gravity constant.

Using equations ( 7), ( 9), ( 10), ( 11) and ( 12), we can finally write the non-linear model of the robot as

ẋ = f (x, u) . (13)

IV. PROPOSED BALANCE-FORCE CONTROLLER

A. Linearized Dynamics

The dynamics given by equation ( 13) can be used to derive different types of control. In this paper we propose to balance the robot around a fixed, stationary desired equilibrium state x * . But even in that case, the non-linear dynamics are complex to tackle directly; therefore, to simplify the control synthesis, the local dynamics are approximated by linearization around this desired state. To simplify the notation, we represent the state error x ∆ between the actual state and the desired one by using an operator noted ∆. For the positions and velocities in the state vector, ∆ represents the Euclidean difference. As for the orientations, it is the axis-sine of angle representation of the relative orientation between the orientation matrices. The linearization is achieved as detailed in the Appendix.

It is important to note that this linearization is very different from the one commonly performed with the inverted pendulum. Indeed, no assumptions are made either on the kinematics of the CoM, on the position nor on the orientation of the contacts. So, this linearization is not less precise in the case of multiple non-coplanar contacts with different stiffness and damping (including point and edge contacts).

Having linearized the reduced model ( 13), we define the matrices A and B (given in the Appendix) such that

ẋ∆ = Ax ∆ + Bu ∆ . ( 14
)
This model is well suited for a state feedback control law. However, it would be purely based on kinematics and cannot track reference forces except through the viscoelastic model. To add more reliability to the knowledge of the contacts, the forces measured by force sensors are going to be appended to the positions as explained in the next section.

B. Balance-Force Control Law

One main contribution of the proposed controller is to include force tracking in the same control loop as CoM and kinematic tracking. However, it is not possible for the controller to use all the available actuation for kinematics and force control at the same time; thus there is a need to establish a trade-off between the conflicting kinematics and force control, which is common to have. Therefore, we need to build a new state vector where the contact forces and torques are combined with the positions and orientations of the corresponding contact bodies. Hence, the vector z is defined as

z ≜ 0 1×12 z 1 • • • z n ⊤ , (15)
where

z i = f ⊤ s,i t ⊤ s,i 0 1×3 0 1×3 ⊤ . The vectors f s,i
and t s,i ∈ R 3 are respectively the force and torque at contact i, scaled to the positions and orientations by dividing by the stiffness of the contact, and written in the base frame to make their appendage to the positions and orientations possible. The actual forces and torques are measured with force sensors at the contacts. Now, we can combine both vectors x and z in what we define as the new state error y ∆

y ∆ ≜ (I -W ) x ∆ + W z ∆ . ( 16
)
where I is an identity matrix, and W is a weight matrix, used to set the trade-off between force tracking and position tracking. One way of defining it is as

W ≜ diag (0 1×12 , w 1 , • • • , w nc ) , (17)
where diag(.) is an operator that gives a square matrix, having on its diagonal the values given in between the parenthesis and zeros elsewhere, and

w i ≜ w f i w ti 0 1×3 0 1×3 , (18)
where w f i ∈ R 3 and w ti ∈ R 3 are vectors having values between 0 and 1 that respectively multiply f s,i and t s,i in the vector z.

Considering that the forces and torques in z can be written as functions of the variables in x using ( 4) and ( 5), we can write, after linearization, the vector z ∆ as a matrix M (given in the Appendix) multiplying the vector

x ∆ z ∆ = M x ∆ . (19)
Using equation ( 19) in ( 16), we can write y ∆ as:

y ∆ = N x ∆ , (20)
with N = I -W + W M and W is chosen so that matrix N is non singular.

We can now write the dynamics of y ∆ as

ẏ∆ ≃ A y y ∆ + B y u ∆ , (21)
with A y = N AN -1 and B y = N B.

With these new dynamics, we use a linear quadratic regulator (LQR) to minimize over the control space the following quadratic cost L such that

L = ∞ t0 y ∆⊤ Qy ∆ + u ∆⊤ P u ∆ dt, (22)
where Q and P are the weight matrices for the state and control respectively. The problem then boils down to solving a Riccati Equation, which provides us with the optimal gain matrix K such that u ∆ = -Ky ∆ induces the minimum cost L. The calibration of the weight matrices Q and P allows to modify the behavior of the controller.

C. Integration Within a Multi-Objective Motion Generator

The proposed control is using 6 × n c variables of control. Many robots, especially humanoid ones, are equipped with more dof. and are able to perform other concurrent tasks. To deal with the redundancy and with additional objectives, we propose to use a whole-body motion solver based on a quadratic program (QP). This optimization problem minimizes the tracking error of different weighted objectives. It is convenient to have the acceleration α as a decision variable to take into account dynamical constraints. The QP calculates the optimal reference acceleration αr , by solving

αr = argmin ξ ∥W task (A ob ξ -b ob )∥ 2 , s.t. A eq ξ = b eq , Aξ ≤ b, l b ≤ ξ ≤ u b , (23)
where W task is a positive diagonal matrix made up of diagonal weighting matrices for each objective. The matrices A ob , A eq , and A and the vectors b ob , b eq ,b, l b , and u b contain the corresponding objectives and constraints. Note that our controller naturally fits into this motion solver since it directly provides the desired Cartesian accelerations of the contact bodies and needs only simple Jacobians. For details about how the tasks and constraints for the QP are formulated see [START_REF] Cisneros | Partial Yaw Moment Compensation Using an Optimization-Based Multi-Objective Motion Solver[END_REF] and [START_REF] Cisneros | Robust Humanoid Control Using a QP Solver with Integral Gains[END_REF]. The difference is that in our case, our proposed controller will generate the acceleration for the contacts, which are stored in b ob , and the QP generates the reference joint accelerations only. Figure 1 gives an overall view of the robot and the control structure. The robot's model is reduced and used by the controller to generate the accelerations for the contacts, whose objectives, alongside the posture objective and constraints, are minimized using the QP. The latter generates the reference accelerations for the torque controller to produce joint torques using [START_REF] Escande | Contact planning for acyclic motion with tasks constraints[END_REF]. A passivity based term δτ is added to the joint torques. δτ is calculated using

δτ = (C + λH) (α r -α) , (24)
where α r is the reference velocity vector, obtained by integrating αr , λ is a constant, H is the Inertia matrix, and C is the Coriolis matrix. For more details about the passivityterm, see [START_REF] Cisneros | Robust Humanoid Control Using a QP Solver with Integral Gains[END_REF].

V. SIMULATIONS AND RESULTS

The control is tested in simulation on the biped robot described in [START_REF] Cisneros | Partial Yaw Moment Compensation Using an Optimization-Based Multi-Objective Motion Solver[END_REF], designed using Matlab Simscape Multibody. The robot has 26 dof in total: 6 dof for each leg and 7 for each arm, and a total mass of 77 Kg. The Simscape Multibody Contact Forces Library was used to generate the non-rigid contact models. Each corner of the base of a foot or hand has a virtual sphere attached to it, generating linear and angular stiffness, plus damping. The simulations were run while allowing the simulator to automatically choose the max and min time step size using a variable step solver. This setting provides a simulation with high precision that is able to break down highly dynamical effects and is computationally very costly [START_REF] Li | Reinforcement Learning for Robust Parameterized Locomotion Control of Bipedal Robots[END_REF].

The control is compared to the stabilizer of [START_REF] Murooka | Humanoid Loco-Manipulations Pattern Generation and Stabilization Control[END_REF], where the effect of the compliance in the environment is modeled as a first order lag system on the ZMP. This stabilizer replaces our controller in Figure 1, generating feet damping control. Also, the QP in this case generates the contact forces alongside the reference accelerations, just like in [START_REF] Cisneros | Robust Humanoid Control Using a QP Solver with Integral Gains[END_REF]. We're going to refer to this stabilizer as the feet damping controller.

The tests were done on scenarios with different starting postures, as illustrated in Figure 2. We start the simulation with the robot in contact with the environment, slightly above the equilibrium state, to avoid introducing too much initial energy. Then, the robot is perturbed to assess control robustness. The robot is expected to maintain its balance and to keep the contacts with the environment, which is why we focus on the CoM position of the robot, and on the CoP of each contact (in the contact frame) while respecting the friction constraint. Therefore, we set the weight matrix W such that w f i = 0 0 0.5 and w ti = 0.5 0.5 0 in all scenarios. Furthermore, to challenge the controller, white noise is added to the measured forces, and modeling error is introduced by providing reduced stiffness and damping values to the controller compared to the actual values that we set for the environment. The stiffness is reduced by 20% and the damping by 10%.

A. Two Contacts Case

The robot, in this case, is in a half-sitting configuration while making contacts with its feet on the ground, as shown in Figure 2 (a). The stiffness and damping constant parameters at the contacts are given in Table I for the Right Foot (RF) and Left Foot (LF) under Scenario (a). We choose the LQR control matrices such that

Q = diag( Q com Q c Q c ), with Q com = ( 1e4I 3 I 3 3e4I 3 I 3 ), Q c = ( 3e4I 3 I 3 3e4I 3 I 3 ) and P = diag(I 6
). The robot is hit at its base with a perturbation of 250N at t = 0.5s for 0.1s along the x-axis.

TABLE I CONTACTS STIFFNESS AND DAMPING FOR THE 3 SCENARIOS

Our proposed controller manages to keep the robot's balance, while the feet damping controller struggles until the robot falls down. This is clear when looking at the CoM position in Figure 3 and the CoP of the right foot in Figure 4. The feet damping controller struggles to generate the proper forces in the non-rigid environment, which explains the large fluctuations in the foot CoP and the inability of the controller to maintain the robot's balance after the perturbation. On the other hand, our proposed controller manages to minimize the error of the CoM and to keep the CoP of each foot within its support zone (each foot has a dimension of x : 0.23m × y : 0.13m). The maximum value of the friction ratio at the feet does not exceed 0.176, and no slippage occurs.

To gauge the robustness to modeling errors of our controller, the same scenario is repeated with the same tuning while reducing the linear and angular stiffness of the contacts by 50% and by 70% in the controller. Despite stronger oscillations of the CoM and a bigger steady-state error as seen in Figure 5, the controller is still able to maintain the robot's balance. Reducing the stiffness in the controller by more than 70% might cause the robot to fall due to a CoM error that places the CoM outside of the support region. A finer tuning is perhaps needed to make the controller handle such modeling errors.

B. Four Contacts Case

Here, we show two different scenarios while the robot is making four contacts with the environment, as illustrated in Figure 2: the robot is putting its hands on a table-like surface (scenario (b)), and the robot is making non-coplanar contacts with its hands (scenario (c)). The stiffness and damping constants at the contacts for the Right Foot (RF), Left Foot (LF), Right Hand (RH) and Left Hand (LH) are given in Table I under scenarios (b) and (c). In these multicontact scenarios, we are not going to present the data from the feet damping controller. This is because the robot was already struggling to generate the contact forces through the QP even before any perturbation due to the non-rigidity of the contacts, and also because the fluctuation in the CoP of each contact was very strong. Therefore, the contacts broke quickly and the robot fell. With our proposed approach, the robot is hit with a perturbation of 200N at t = 0.5s and then again with 250N at t = 2.5s for 0.1s. In the case of scenario (b), a bias of 5 • to the orientation of the base is added to the robot at the start of the simulation to test the robustness of the control due to estimation errors. In this case, we set the control matrices such that

Q = diag( Q com Q c Q c Q c Q c ), with Q com = ( 1e4I 3 100I 3 3e4I 3 100I 3 ) , Q c = ( 3e4I 3 100I 3 3e4I 3 100I
3 ) and P = diag(I 6 ). By setting the angular positions and velocities gains at 100, we were able to minimize the CoM error while regaining the default orientation of the base. We noticed that giving too much weight to the orientations results in a worse behavior. In scenario (c), we set Q com = ( 1e5I 3 I 3 1e4I 3 I 3 ), Q c = ( 1e4I 3 I 3 1000I 3 I 3 ) and P = diag(I 6 ).

We see in Figures 6, 7 and 8 that the proposed controller is able to maintain the balance of the robot, minimizing the CoM error and maintaining the CoP in its support zone (each hand has a dimension of x : 0.105m × y : 0.135m). No slippage of contacts occurred, and the maximum value of the friction ratio at each contact in each scenario is given in Table II.

VI. CONCLUSION AND FUTURE WORK

In this paper, we proposed a controller to manage balance and force control in the same loop in task space for multiple flexible-contacts configurations. Using a flexibility model of the contacts with centroidal dynamics, a trade-off between force feedback and the kinematics of the system is established to create an LQR control that generates the required accelerations of the contacts for the QP solver. The performance and robustness of the controller has been demonstrated in simulations.

The main focus of this paper is the definition of the model, the control variables and an outline of the architecture of this control framework. Thus, one limitation of the presented control policy is that it deals only with static references. While adding or removing a contact in a statically stable trajectory is possible using this method, more dynamic cases require an adaptation through a linearization around a timedependent reference trajectory, using the same architecture. Another limitation to mention is that the controller doesn't generate forces that are bound to respect unilaterality constraints nor to stay in friction cones. However, this issue can be solved by adding an explicit constraint to guarantee the feasibility of the control trajectories using a QP, similar to [START_REF] Wieber | Trajectory Free Linear Model Predictive Control for Stable Walking in the Presence of Strong Perturbations[END_REF].

The next step would be to implement this control on a real legged robot, where modeling errors of the stiffness of the environment, estimation errors and noisy force signals could pose a challenge. To overcome these challenges, the trade-off weighting of positions and forces will be crucial, requiring to give higher weight to the most reliable signals. Further extension of this work will include an integration of the controller within a motion planner to achieve dynamic locomotion in the context of compliant multi-contacts.

APPENDIX

The operator ∆ represents the Euclidean difference for positions and velocities in the state vector for instance v ∆ = v -v * . As for the orientations, it is the axis-sine of angle representation of the relative orientation between the orientation matrices, for instance Ω ∆ = Θ(RR * ⊤ ).

The formula used in the linearization process is:

(XY + Z + S(w)) ∆ ≃ X ∆ Y * + X * Y ∆ + Z ∆ + S(w ∆ ) (25)
where X, Y and Z all ∈ R 3×3 , and w ∈ R 3 . (.) ∆ is the error of (.), (.) * is its reference value. S(w) is the skew-symmetric operator of vector w. As we can notice, the operator ∆ functions similarly to the derivative operator. As for the term

Θ R c,i,r R ⊤ c,i = Ω R c,i,r R ⊤ c,i,b R ⊤ ,
the linearization was done in the limits where the estimation sinθ ≈ θ can be used, as the following:

Ω (R M ) ∆ ≃ 1 2 V ec R M -R ⊤ M ∆ = C b Ω ∆ + R * Ω ∆ c,i,b , (26)
where

C b = 1 2 3 i=1 S (e i ) R * M S (e i ), R M = R c,i,r R ⊤ c,i,b R ⊤ , V ec (S (ω)) = ω, e 1 = 1
0 0 , e 2 = 0 1 0 , and e 3 = 0 0 1 .

At the end of the linearization process, we get equation ( 14), with matrices A and B given by

A =        F 0 F 1 F 2 . . . F n 0 D 1 0 . . . 0 0 0 D 2 0 0 . . . . . . . . . . . . . . . 0 0 0 0 D n        , B =        0 0 • • • 0 G 1 0 0 0 0 G 2 0 0 . . . . . . . . . . . .
F i =     0 0 0 0 0 0 0 0 F i31 0 -1 m K f,d,i R * 0 F i41 F i42 F i43 F i44     , D i =     0
with

T i =     R * ⊤ T i12 T i13 -T i12 S R * p * c,i 0 T i22 0 R * ⊤ K -1 t,p,i K t,d,i 0 0 0 0 0 0 0 0     , V i =     V i11 0 V i13 0 0 -R * ⊤ C b,i R 0 V i24 0 0 0 0 0 0 0 0     , T i12 = -R * ⊤ S R * p * c,i -T i13 S(R * ṗ * c,i ) + S (ω * ) S R * p * c,i , T i13 = R * ⊤ K -1 f,p,i K f,d,i , T i22 = -R * ⊤ C b,i + K -1 t,p,i K t,d,i S R * ω * c,i , V i11 = R * ⊤ I + K -1 f,p,i K f,d,i S (ω * ) R * , V i13 = R * ⊤ K -1 f,p,i K f,d,i R * , V i24 = R * ⊤ K -1
t,p,i K t,d,i R * .

Fig. 1 .

 1 Fig. 1. Diagram showing the robot, its reduced model and the control loop. The world, CoM and contact frames are highlighted for the reduced model. The yellow block is our proposed controller.

Fig. 2 .

 2 Fig. 2. Starting postures of scenario (a) (left), (b) (middle) and (c) (right)

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. CoM position error (x component) in scenario (a)

Fig. 6 .Fig. 7 .

 67 Fig. 6. CoM position error (x component) in scenarios (b) and (c)

Fig. 8 .KR

 8 Fig. 8. Right foot and right hand CoP (x component) in scenario (c)

  ,p,i + K f,d,i S (ω * )) R * , F i41 = R I -S f * c,i -S R * p * c,i,b K dp R * , K dp = K f,d,i S (ω * ) + K f,p,i , F i42 = R I K t,p,i C b,i R, F i43 = -R I S R * p * c,i,b K f,d,i R * , F i44 = -R I K t,d,i R * .The matrix M in (19) is written as:

M. Hamze and A. Benallegue are with Université Paris-Saclay, UVSQ, Laboratoire d'Ingénierie des Systèmes de Versailles, France. E-mails: marwan.hamze@uvsq.fr, abdelaziz.benallegue@uvsq.fr

[START_REF] Escande | Contact planning for acyclic motion with tasks constraints[END_REF] M. Benallegue, R. Cisneros and A. Benallegue are with CNRS-AIST JRL (Joint Robotics Laboratory), IRL, National Institute of Advanced Industrial Science and Technology (AIST), Tsukuba, Ibaraki, Japan. Emails: mehdi.benallegue@aist.go.jp, rafael.cisneros@aist.go.jp.


F Ben Abderkader

G Kairo

S Tchamitchian

M Bonnet

M Cousin

N Barbouche

L P Belzunces

J Luc Brunet

Effects of Clothianidin exposure on semen parameters of honey bee drones

Keywords: pesticides, drones, spermatozoa, fertility, metabolic state

Many problems have been reported on honey bees colonies including fertility problems of queens resulted in production failure. Pesticides can be the cause of this failure in connection with the quality of sperm drone. Thus, the aim of this study was to assess the influence of exposure to syrup contaminated with clothianidin at 0.1 µg/L on semen parameters of drones. Results showed a significant decrease of semen volume and sperm concentration and an increase in sperm mortality rate. As for the energetic state, clothianidin increased cell redox potential, the ATP content of spermatozoa as well as the lactate dehydrogenase activity (LDH). It was concluded that exposure to clothianidin during the sexual maturity of drones could affect the semen quality.

Introduction

Populations of honey bees are declining worldwide for more than a decade (Grossman 2013). Colony losses are attributed to numerous causes including exposure to pesticides (DeGrandi-Hoffman et al. 2013). Increasingly, there are concerns about the way by which pesticides may weaken and kill bee colonies at sublethal levels (Frazier et al. 2008). In parallel, the increased use of agrochemicals coincided with increased beekeeper complaints on problems associated with queen performance in the hive (Burley et al. 2008). Many beekeepers have noted problems with supersedure after the introduction of a new queen as well as an inability of colonies to naturally renew queen (Sanford, 2001 in [START_REF] Pettis | A scientific note on varroa destructor resistance to coumaphos in the united states[END_REF]. It is known that the reproductive live of queen depends, among others, on the success of her mating with drones(Connor 2010). Decreasing sperm quality of drones might play a role in the current problems associated with queen performance (Burley et al. 2008). Queens mated with sperm of low viability would become drone layers faster than well-mated queens (Collins 2000). Pesticides are among factors that have an effect on drone semen quality. Burley et al. (2008) report that a decrease in quality of sperm contained in the spermatheca of queens mated with miticide-exposed drones could lead to queen loss and supersedure. Miticides used in honey bee colonies can also reduce drone survival [START_REF] Rinderer | Varroa in the mating yard: 1. The effects of varroa jacobsoni and apistan on drone honey bees[END_REF]) and production [START_REF] Rinderer | Varroa in the mating yard: 1. The effects of varroa jacobsoni and apistan on drone honey bees[END_REF], body and mucus gland weights [START_REF] Rinderer | Varroa in the mating yard: 1. The effects of varroa jacobsoni and apistan on drone honey bees[END_REF]) and spermatozoa production (Fell and Tignor 2001). Besides, many studies have reported adverse effects of pesticides on the reproductive system of vertebrates. Bal et al. (2012a) report reduced sperm production in rats exposed to imidacloprid, a neonicotinoid, at 2 mg/kg bw/day, a dose representing about 1/250 th of the LD50 per day. Exposure to clothianidin at low doses during critical stages of sexual maturation of male rats has moderate detrimental effects on reproductive system but more severe effects are observed at higher dose levels (Bal et al. 2012b). The reproductive system may be more sensitive to exposure to clothianidin early in the development (prenatal and early postnatal) (Bal et al. 2012b). Moreover, it has been reported that some neonicotinoid pesticides induce oxidative stress in the serum of rat [START_REF] Mohany | Immunological and histological effects of exposure to imidacloprid insecticide in male albino rats[END_REF]), liver (El-Gendy et al. 2010) and testis [START_REF] Zhang | Oxidative stress: Role in acetamiprid-induced impairment of the male mice reproductive system[END_REF] of mice. Honey bee colonies are often exposed not only to miticides but also to other insecticides through contaminated nectar and pollen stored and subsequently shared among nestmates including drones (DeGrandi-Hoffman et al. 2013) which could affect their role inside the colony. The objective of this study was to determine the effects of clothianidin on the quality of semen from drones daily exposed up to their sexual maturity through contaminated food gathered in field and brought to the hive by foragers. To this end, analyses were focused on semen volume, concentration, mortality rate, redox potential, energetic state (ATP) of spermatozoaand lactate dehydrogenase (LDH) activity.

Materials and Methods

Experiments were performed from May to July 2013 (two repetitions) in Avignon (France) with Apis mellifera L. colonies carefully monitored for their health status. To stimulate drone production, drone combs were introduced in 10 strong colonies (35,000-45,000 bees/colony; 1 comb/colony) 25 days before the beginning of experiments. These colonies were treated conventionally with Amitraz in September 2012 to control varroa population. The day before emergence, brood frames were placed into an incubator in the dark at 34°C and 70% relative humidity. Queenless and droneless nuclei composed of 5000 honey bees, one brood frame, and four empty frames with no food storage were placed in an outdoor tunnel (32 x 8 m) covered with an insect proof mesh divided 4 compartments. Two nuclei were placed in each compartment (8 x 8 m). To avoid overexposure to the sun, umbrellas were placed above each nucleus and food feeder. No food was supplied inside the nuclei. Three hundred emerging drones were introduced in each droneless nucleus and enabled for free movement in the colony.

Exposure to clothianidin

The honey bee nuclei were introduced in the tunnel three days before drone introduction to train bees to visit a feeder placed at equal distance from the nuclei. This feeder consisted of a yellow cup containing sugar syrup (50% sucrose solution, w/v) on which a piece of cork and filter paper were placed on the surface to enable the honey bees to land. Treated nuclei received sugar syrup with a final concentration of 50% sucrose, 0.1% DMSO (v/v) and 0.1 µg/L clothianidin. Great care was taken to avoid the direct contact of the sunlight with the feeder. Controls received sugar syrup containing 50% sucrose and 0.1% DMSO. Honey bees were fed or exposed to clothianidin from 9:00 to 12:00 a.m. each day for a period of 20 days. Outside this period, honey bees received water and pollen ad libitum. The amounts of syrup and pollen foraged were measured daily in each compartment. Thus, during their sexual maturity, drones were chronically fed with a food gathered to the hive by foragers, as in natural conditions. In total, 22 queenless nuclei were used in the experiment, 10 for treatment and 12 for control (4 and 6 for the first experiment in 2013, 6 and 6 for the second experiment in 2013).

Drone survival rate, maturity rate and semen collection

Drones were collected after 20 days. The survival rate was recorded. Semen was gathered by a manual eversion of the drone endophallus. Briefly, drones were stimulated to ejaculate by pressing on the thorax, which usually resulted in eversion of the endophallus. Semen was collected from the tip of the endophallus with a glass capillary connected to a syringe filled with Kiev solution (36 g/L trisodium citrate, 3.6 g/L sodium bicarbonate, 0.6 g/L potassium chloride, 5 g/L glucose, 3 g/L sulfanilamide, pH 8.5, osmotic pressure = 486 mOs/ml). The average volume of semen was determined. During semen collection, the maturity rate of drones from each colony was assessed by the ability to provide sperm after stimulation. The average semen volume per drone was also determined. A part of the fresh semen from each sample was used for the determination of the spermatozoa concentration, mortality rate and metabolic activity (as assessed by measuring the reducing potential, adenosine triphosphate (ATP) content) and another part of fresh semen from each sample was transferred into a microcentrifuge tube, diluted two-fold in Kiev solution and centrifuged at 4°C for 20 min at 16,000 g. The supernatant (diluted seminal plasma) was kept for LDH assay.

Spermatozoa concentration

The semen was diluted in the Kiev solution (1/1500 th ) and spermatozoa were counted under a phase contrast microscope by using a Neubauer improved Type Petroff cell with 0.02 mm of depth. Five repetitions were done for each sample.

Sperm mortality

Sperm mortality was assessed with the conventional dead-cell stain propidium iodide (PI) from the Live/Dead ® Sperm viability kit (Life Technologie, France). This test was performed in 96-black-well microplate. According to the recommendations of the manufacturer, propidium iodure was diluted in DMSO at a final concentration of 60 µM. Five µl of PI solution were added to each well containing 1.10 7 spermatozoa. The microplate was incubated at 37°C for 5 min before measuring the fluorescence intensity.

Sperm reduction potential

Sperm reduction potential was estimated by the capacity of spermatozoids to reduce resazurin into resorufin, which is fluorescent (λexc = 569 nm, λem = 586 nm) and absorbs at 470 nm (Erb and Ehlers 1950;[START_REF] Rousseau | Apis mellifera (hymenoptera: Apidae) drone sperm quality in relation to age, genetic line, and time of breeding[END_REF]. The test was performed in quadruplicate with a Prestoblue® kit (Invitrogen). Ninety µL of diluted semen samples and 10 µL of Prestoblue® were put in a 96 white-well microplate and kept in darkness for 10 min before measuring the absorbance at 570 nm according to the recommendations of the manufacturer.

ATP content

ATP content was determined with ATPlite® kit (PerkinElmer, Courtaboeuf, France) in the same wells as the viability test. The ATPLite® assay system is based on the production of light caused by the reaction of ATP with added luciferase and D-luciferin. The emitted light is proportional to the ATP concentration within certain limits. As recommended by the supplier, 50 µl of a mammalian cell lysis solution (PerkinElmer®) were added to the 100 µl of diluted semen in each well and the microplate was covered and gently shaken for 5 min with an orbital shaker at 700 rpm. Then, 50 µl of substrate solution was added, the microplate was covered and shaken for 5 min at 700 rpm and then was kept in the darkness for 10 min before measuring the luminescence intensity (LI).

Lactate deshydrogenase

The principle of the measure of LDH activity is based on the reduction of NAD + to NADH/H + during the conversion of lactate to pyruvate by LDH. The assay medium contained diluted semen and 50 mMtriethanolamine buffer, 5 mM EDTA, 0.2 mM reduced β-Nicotinamide adenine dinucleotide (β-NADH) and 2 mM sodium Pyruvate.LDH activity was spectrophotometrically followed at 340 nm for 5 min.

Statistical analysis

Mann-Whitney U test was used to determine differences between groups. Results are presented as mean ± S.E. Values were considered statistically significant if P< 0.05. The SPSS/PC program (Version 16.0) was used for the statistical analysis.

Results and discussion

The nicotinic acetylcholine receptors (nAChRs) of the nervous system of insects are the primary and selective targets for Neonicotinoid insecticides including clothianidin [START_REF] Tomizawa | Selective toxicity of neonicotinoids attributable to specicifty of insect and mammalian nicotinic receptors[END_REF]. As some nAChR subunits are expressed in other organs, any neonicotinoid may have multiple effects outside the central nervous system (Bal et al. 2012b). To our knowledge, no studies have been carried out on the toxicity of Neonicotinoid on drones exposed during the sexual maturity stage at which sperm has completed migration from the testes to the seminal vesicles and when mucus glands are fully developed [START_REF] Rousseau | Apis mellifera (hymenoptera: Apidae) drone sperm quality in relation to age, genetic line, and time of breeding[END_REF]. In this context, drones were exposed to clothianidin via food brought back to the hive by foragers during the period from emergence to drone sexual maturity. Hence, drones directly consuming the syrup with clothianidin at 0.1 μg/L upon the return of foragers. The apparent effects show no significant effects on the foraging behavior of workers during the 20 days of exposure either for contaminated syrup (Figure 1A) or the pollen (Figure 1B).One recent study showed that limited sublethal effects on the tested behavioral functions of winter bees were observed following a chronic exposure to clothianidin at 15 lg/kg(Alkassab and Kirchner 2016). Although many factors are different under field conditions compared to semi-fields conditions, ours results play an important role in assessing the exposure risk. A total of 46.04 ± 25.33% of survived drones was recorded in exposed group (33.88 ± 19.6% for control) (Figure 2). The percentage of mature drones recorder in our study was higher in exposed group (60.97 ± 4.33%) compared to control (56.35 ± 8.56%) (Figure 3A) but statically these results are not significant. The quantity and quality of semen parameters such as ejaculate volume, sperm concentration and sperm mortality are used as biological endpoints to evaluate drone fitness and spermatozoal competitiveness (Gençer and Kahya 2011).In this context, the quality of semen was studied in greater detail by analyzing sperm count and sperm mortality. The mean semen volume per drone as shown in Figure 3B was lower in drones exposed to clothianidin (median of 0.7 µl/drone vs 0.82 in controls, P<0.05). The exposure to clothianidin induced a decrease in total spermatozoa concentration (median of 8.20 × 10 6 spermatozoa vs 10.88 × 10 6 spermatozoa in controls, P< 0.001, Figure 4A) and an increased spermatozoa mortality rate (median of 39.92% vs. 47.57% in controls, P< 0.001, Figure 4B). In animals, a chronic exposure to clothianidin causes a decrease of spermatozoa count and the appearance of deformed spermatozoa (Watts 2011). Bal et al. (2012b) found also that exposure to clothianidin during 90 days decreases significantly epididymal sperm concentration in male rats. Unlike [START_REF] Straub | Neonicotinoid insecticides can serve as inadvertent insect contraceptives[END_REF] who found that exposure drones to clothianidin before emergence did not affect the sperm quantity, our results are in agreement with (Kairo et al. 2016) who find that exposure to Fipronil induced a decrease in total spermatozoa concentration and increase in sperm mortality. This suggests that reproductive effects of clothianidin are not limited to vertebrates and the number of species that could be affected seems much underestimated. At the cellular level, reducing potential corresponds to the ability of cells to metabolically reduce compounds or to passively reduce oxidants resulting from metabolic activity or ionizing rays. Reducing potential was increased by exposure to clothianidin (1.03 absorbance units (AU)) compared to control (0.85 AU, P< 0.001, Figure 5A). The cellular ATP content in spermatozoa produced by mitochondria was also increased (median of 2048.48 luminescence intensity (LI) vs. 1457.28 LI in controls, P<0.05, Figure 5B). The redox potential is important for control of metabolism. Redox potentials are used to infer the direction and free energy cost of reactions involving electron transfer, one of the most ubiquitous and important types of biochemical reactions (Milo and Phillips 2016). According to [START_REF] Trachootham | Redox regulation of cell survival[END_REF], oxidative stress can either enhance cell survival or promote cell death, depending on the magnitude and duration of the stress, the genetic background and redox states of the cells. The same authors reported that redox regulation of the key factors affecting cell death/survival is often bifurcated (i.e., the same protein can either be activated or inhibited by redox alteration). Though ATP is often claimed to be the energy currency of the cell, in fact, for the energetic balance of the cell the carriers of reducing power are themselves no less important (Milo and Phillips 2016). Sperm energy metabolism can be also characterized by the activities of marker enzymes reflecting the capacities of various metabolic pathways such as lactate dehydrogenase (LDH) (Al-Lawati and Bienefeld 2009).The LDH activity in seminal plasma of drones exposed to clothianidin was also significantly higher (1.54 ± 0.41mAU/min/10 6 spz) compared to control (0.97 ± 0.16mAU/min/10 6 spz, P< 0,05, Figure 6). LDH is essential for metabolic processes which provide energy for survival, motility, capacitation and fertility of spermatozoa [START_REF] Sirat | Effect of cryoprotectants on release of various enzymes from buck spermatozoa during freezing[END_REF]. The seminal fluid LDH can be used as good indicator of sperm viability [START_REF] Stamatiadis | Relationship between ldh activity in seminal plasma and viability of ram spermatozoa[END_REF]). It has thus been proposed that increased LDH levels, as shown in treated group in our study, can be used as a good indicator of less of integrity of the plasma membranes (Dube et al. 1982). Moreover, extracellular activity LDH and lactate increases under the condition of oxidative stress, since the cell integrity can be disrupted during the lipid peroxidation process (Jovanovic et al. 2010). Although spermatozoa exhibit antioxidant defenses, an overload of these defenses can lead to cell death, explaining the high mortality rate of spermatozoa in drones exposed to clothianidin.

Conclusion

In conclusion, this study is the first report investigating the effect of exposure to clothianidin, via contaminated food brought to the hive by foragers, on the quality of honey bee drone semen. Clothianidin had negative effects and elicit a decrease of semen volume and sperm concentration, and an increase of sperm mortality rate and other sperm parameters. These results reveal for the first time the toxicity of neonicotinoid insecticides to reproduction in honey bee drone exposed during their sexual maturity. The results were obtained from new investigation methods, and therefore introduce a new perspective on the studies on queen performances and colony development.

Figure 1 .

 1 Figure 1. Effect of chronic clothianidin exposure on the cumulative foraged syrup (A) and pollen (B). The data represent the mean ± the standard deviation of the daily cumulative foraged quantity. For each treatment, the data correspond to the set of values from 2 experiments conducted in 2013 (n =6 for control and n= 5 for clothianidin). Vertical bars represent standard error of the mean

Figure 2 .

 2 Figure 2. Effects of chronic clothianidin exposure on drone survival rate. The data represent the mean values ± standard deviations obtained from drone populations recovered from hives in 2 experiments monitored in 2013 (n=12 for control and n=10 for clothianidin). Vertical bars represent standard error of the mean

Figure 3 .

 3 Figure 3. Effects of chronic clothianidin exposure on drone maturity rate (A) and on the semen volume (B). The data represent the mean values ± standard deviations obtained from drone populations recovered from hives in 2 experiments monitored in 2013 (n=12 for control and n=10 for clothianidin) Vertical bars represent standard error of the mean. *Comparison of control with clothianidin-treated group (P < 0.05)

Figure 4 .

 4 Figure 4. Effects of chronic clothianidin exposure on spermatozoid (SPZ) number (A) and semen mortality rate (B). The data represent the mean values ± standard deviations obtained from drone populations recovered from hives in 2 experiments monitored in 2013 (n=12 for control and n=10 for clothianidin). Vertical bars represent standard error of the mean. *** Comparison of control with clothianidin-treated group (P < 0.001)

Figure 5 .

 5 Figure 5. Effects of chronic clothianidin exposure on Reducing potential (A) of semen expressed in absorbance units (AU) and on ATP content expressed in luminescence intensity (LI). The data represent the mean values ± standard deviations obtained from drone populations recovered from hives in 2 experiments monitored in 2013 (n=12 for control and n=10 for clothianidin). Vertical bars represent standard error of the mean. * and *** Comparison of control with clothianidin-treated group (P < 0.05) and (P < 0.001) respectively

Figure 6 .

 6 Figure 6. Effects of chronic clothianidin exposure on LDH activity expressed in milli absorbance unities (mUA). The data represent the mean values ± standard deviations obtained from drone populations recovered from hives in 2 experiments monitored in 2013 (n=12 for control and n=10 for clothianidin). Vertical bars represent standard error of the mean. *Comparison of control with clothianidin-treated group (P < 0.05)

BEN ABDERKADER et al. (2018) / Journal of new sciences, Agricultureand Biotechnology, 59 (1), 3791-3798

Volume 59 (1). Published November, 01, 2018 www.jnsciences.org E-ISSN 2286-5314 Acknowledgments This work was supported in part by INRA (France), Agence Universitaire de la Francophonie (AUF) and by grants from the European Union FEAGA Beekeeping program managed by the INRA-France AgriMer Agreement 11-42 R 5. References