P Shirkovskiy

A Laurin

N Jeger-Madiot

D Chapelle

M Fink

R K Ing

Airborne ultrasound surface motion camera: application to seismocardiography

Keywords: Non-contact ultrasonic imaging, surface motion, non-contact seismocardiography, airborne ultrasound vibrometry, cardio-vascular signals, cardiac mechanics

The recent achievements in the accelerometer-based seismocardiography field indicate a strong potential for this technique to address wide variety of clinical needs. Recordings from different locations on the chest can give a more comprehensive observation and interpretation of wave propagation phenomena than a single-point recording, can validate existing modeling assumptions (such as the representation of the sternum as a single solid body), and provide better identifiability for models using richer recordings. Ultimately, the goal is to advance our physiological understanding of the processes to provide useful data to promote cardiovascular health. Accelerometer-based multichannel system is a contact method and laborious for use in practice, also even ultralight accelerometers can cause non-negligible loading effects. We propose a new contactless ultrasound imaging method to measure thoracic and abdominal surface motions, demonstrating that it is adequate for typical seismocardiogram use. The developed method extends non-contact surface-vibrometry to fast 2D mapping by originally combining multi-element airborne ultrasound arrays, a synthetic aperture implementation and pulsed-waves. Experimental results show the ability of the developed method to obtain 2D seismocardiographic maps of the body surface 30×40 cm 2 in dimension, with a temporal sampling rate of several hundred Hz, using ultrasound waves with the central frequency of 40 kHz. Our implementation was validated in-vivo on eight healthy human participants. The shape and position of the zone of maximal absolute acceleration and velocity during the cardiac cycle were also observed. This technology could potentially be used to obtain more complete cardio-vascular information than single-source SCG in and out of clinical environments, due to enhanced identifiability provided by distributed measurements, and observation of propagation phenomena.

INTRODUCTION

Generally Seismocardiogram (SCG) systems have used accelerometers [START_REF] Inan | Ballistocardiography and Seismocardiography: A Review of Recent Advances[END_REF] or microwave radars [START_REF] Tavakolian | Radar mechanocardiography: A novel analysis of the mechanical behavior of the heart[END_REF] to measure the thoracic surface motion produced by the heart beats. With a few exceptions, these methods have been limited to single-point measurements on the chest surface. While automatic annotation algorithms have been developed [START_REF] Laurin | Accurate and consistent automatic seismocardiogram annotation without concurrent ECG[END_REF][START_REF] Khosrow-Khavar | Automatic Annotation of Seismocardiogram with High Frequency Precordial Accelerations[END_REF], and published studies demonstrate potential uses for clinical applications [START_REF] Salerno | Seismocardiography for monitoring changes in left ventricular function during ischemia[END_REF][START_REF] Salerno | Seismocardiographic changes associated with obstruction of coronary blood flow during balloon angioplasty[END_REF][START_REF] Salerno | Seismocardiography for detection of coronary artery disease[END_REF][START_REF] Tavakolian | Precordial vibrations provide noninvasive detection of early-stage hemorrhage[END_REF][START_REF] Laurin | Seismocardiograms return Valid Heart Rate Variability Indices[END_REF], SCG has so far largely remains confined to scientific laboratories. This is in part due to the lack of physiological and physical understanding of the processes that underlie the phenomena, causing lasting uncertainties as to the exact interpretation of the signal's peaks and valleys, as well as some variability in accelerometer placement.

A 2D time-map of chest wall movements would provide more complete data, and thus help in the physiological and clinical interpretation of SCG. Although some initial work has been done on array-acquisition of SCG [START_REF] Cozic | Development of a cardiac acoustic mapping system[END_REF][START_REF] Munck | Body Surface Mapping of the Mechanical Cardiac Activity[END_REF], physiological interpretation of these data remain either absent, or tenuously linked to individual valve events. This study's proposed system uses airborne ultrasonic waves to detect movement in a roughly 30×40 cm 2 area of interest without direct contact. Some other methods were explored to achieve a similar purpose, namely laser speckle imaging [START_REF] Ramachandran | Cardiac Functional Analysis by Laser Speckle Interferometry[END_REF][START_REF] Singh | Reconstruction of sequential cardiac in-plane displacement patterns on the chest wall by laser speckle interferometry[END_REF], scanning laser vibrometry [START_REF] Munck | Body Surface Mapping of the Mechanical Cardiac Activity[END_REF][START_REF] Tabatabai | Novel Applications of Laser Doppler Vibration Measurements to Medical Imaging[END_REF] and stereoscopic video [START_REF] Poh | Non-contact, automated cardiac pulse measurements using video imaging and blind source separation[END_REF]. These latter methods were dismissed due to their respective limitations: the lower sensitivity for techniques using video cameras [START_REF] Wijenayake | Real-Time external respiratory motion measuring technique using an RGB-D camera and principal component analysis[END_REF], the necessity to synchronize heart beat ECG with laser vibrometer in the sequential scanning mode, the need to have a good reflective surface for all optical based systems. Furthermore, due to implementation difficulties all these techniques have usually been limited to 30 acquisition channels on ~ 25×25 cm 2 area [START_REF] Cozic | Development of a cardiac acoustic mapping system[END_REF][START_REF] Munck | Body Surface Mapping of the Mechanical Cardiac Activity[END_REF].

Airborne ultrasonic waves of low power in measurements have the advantage of being lowrisk, and having lower cost of implementation. The first patents for non-contact measurement of body movement with ultrasound were filed in the 70s [17,[START_REF] Northrup | Ultrasonic respiration/convulsion monitoring apparatus and method for its use[END_REF]. Later studies have described and validated ultrasound-based systems for various monitoring applications [START_REF] Min | Noncontact respiration rate measurement system using an ultrasonic proximity sensor[END_REF][START_REF] Jeger-Madiot | Non-contact and through-clothing measurement of the heart rate using ultrasound vibrocardiography[END_REF]. The acoustical methods required for 2D surface vibrometry with time-and velocity-resolutions adequate for SCG have only recently become feasible [START_REF] Jeger | 3D airborne ultrasound vibrometer for the detection of skin surface heterogeneities[END_REF].

Our goal in this study was to demonstrate a new non-contact method for non-plane surface vibrometry. It consists in an airborne ultrasound surface motion camera (AUSMC) that is able to achieve multi-channel SCG data measurements on the human thorax and abdomen.

A secondary goal was to validate the capability of the AUSMC system to perform 2D SCG measurements and to perform qualitative examination of the 2D velocity map under cardiac forces to develop hypotheses of the thoracic and abdominal systems as a whole. The underlying aim was to explore the developed AUSMC system to complement existing techniques for performing meaningful medical measurements, and to understand the mechanical principles behind cardiothoracic movements.

MATERIAL AND METHODS

The AUSMC system is composed of nine emission panels arranged in a square matrix (Fig. 1a) and one reception panel positioned at the same location as the central emission panel. Each emission panel is composed of 32 piezoelectric transducers (Murata MA40S4S) uniformly distributed over an area of 24 ×24 cm 2 dimension and operating around the central frequency of 40 kHz. The central panel combines both the emitting and receiving elements on the same plane. The reception panel is composed of 16×16 square array of receivers (Knowles microphone FG-23629 equipped with 40dB electronic amplifier). The horizontal and vertical pitch of microphones is 16 mm i.e. nearly twice the wavelength in the air. A camera is included in the experimental setup to overlay a visual representation of the participant to the acoustical images.

This configuration addresses several imaging challenges: 1. The physical separation between emitter and receiver arrays enables simultaneous emission and reception. In this configuration the round-trip time is no longer a constraint, and operation with higher frequency rate of ultrasound illumination is allowed. 2. Multiple panels are used to increase the surface of the emission aperture and to ensure an adequate imaging surface, comparable with the size of receiving aperture. At 40 kHz, the wavelength in the air is roughly equal to 8.6mm. The roughness of the skin is less than 300µm. Because of an acoustic impedance four order of magnitude higher than that of the air, the skin acts as a specular surface, unlike the diffusion regime where the acoustic wave of illumination would be randomly backscattered. To ensure that all the imaged surface reflect acoustic waves toward the whole reception array, an emission aperture nine times larger than receiving is used. The whole system works in synthetic aperture mode [START_REF] Jensen | Synthetic aperture ultrasound imaging[END_REF] where after each complete step of illumination the whole acoustic image is computed from the impulse responses between all the emitting transducers and receiving microphones. However, to get a high framerate in the synthetic aperture mode the number of impulse responses to measure must be maintained low. In this context all transducers at the same position on each of the nine emission panels are electrically connected together (Fig. 1b). Finally, the number of independent emission channels are reduced from 288 to only 32. Each independent channel drives nine transducers located on nine distinct emission panels. Despite the fact that nine transducers (each on its own panel) radiate simultaneously and because of long distance between them and specular reflection, an acoustic wave generated mainly from one single transducer is reflected toward the reception aperture for a given point of the observation surface. Acoustic waves from the eight other transducers are reflected toward the outside region.

This setup is valid if the measurement surface is planar and parallel to the whole aperture. It is less efficient when the measurement surface is curved. However, as the experimental results show, the AUSMC system is still usable on the torso on individuals even if some areas are not sufficiently illuminated.

Emission channels are driven by the 32 electrical outputs of a programmable D/A acquisition card D-TACQ ACQ1001. Each channel is connected to 9 transducers and the voltage amplitude should reach 20 V pk-pk. The emission waveform is a linearly modulated frequency signal that ranges from 35 kHz to 45 kHz. Its time duration is 400 μs. The emission channel is excited in a fixed order. Successive emissions are delayed by 200 µs. The emission sequence enclosing the excitation of the 32 emission channels lasts 6.4 ms. It is repeated continuously during 12 s. The framerate is 1/6.4 ms and then equal to 157 Hz.

The receiver array is connected to a 256-channels A/D acquisition card D-TACQ ACQ196 (16 bits, 320 kSamples/s) and reflected wavefronts are recorded continuously during the entire emission time up to 12 s.

SIGNAL PROCESSING

The AUSMC system works in the synthetic aperture mode. During one illumination step all the 32 emission channels emit successively to acquire the impulse responses from one channel to the entire 256 microphones. Thereafter a beamforming is both processed in transmission and reception modes. For the transmission mode only one of nine transducers is taken into account. The choice depends on the point of the observation surface and reflection law assuming the observation surface to be planar. The beamforming in transmission mode is achieved by linearly combining the signals from the 32 output channels. The beamforming in reception mode is realized thereafter and results in an acoustic signal for each point of observation of the surface. Generally, the observed surface is decomposed into 60 x 80 pixels i.e. 5mm pitch for both horizontal and vertical directions. Each illumination step lasts 6.4ms. The following illumination step leads to the same computation process as described below. The surface motion for each pixel is then determined by cross correlating the acoustic signals computed for two successive illumination steps. The out of plane surface displacement amplitude between the two illumination steps is equal to the half of the product between the time delay between the two acoustic signals and air acoustic velocity. The frame rate of the surface motion video is equal to the inverse of the illumination step duration and is equal to 157 Hz. The surface motion computed corresponds to the differential out of plane displacement amplitude of the surface and therefore corresponds to the out of plane velocity of the surface. Out of plane acceleration is obtained by time differentiation. The important advantage of this method is its ability to realize both beamforming in transmission and reception modes from the impulse responses acquired during each illumination step. Finally, it allows to get a high framerate with a good focusing dynamic -i.e. focusing lobe level vs side lobe level -and a reasonable spatial resolution.

SEISMOCARDIOGRAMS

Investigation of the AUSMC device for SCG monitoring was performed on eight healthy male subjects (23-35 years old). The study was approved by the appropriate legal and ethical authority (Comité de Protection des Personnes (CPP) de Sud-Méditerranée II en date du 12/05/2017). The participants are sitting in front of the measurement system -at a 720 mm distance -with blocked respiration. Acquisition time was 12 seconds. AUSMC, accelerometer-based and laser vibrometer (LV) measurements for some points are performed concurrently for comparison purposes.

RESULTS

We reported here a new method for surface motion monitoring, and tested its capability to acquire 2D vibration maps on human torso. The developed device provided in non-contact mode surface acceleration, velocity or displacement maps over time of the 30×40 cm2 area.

Accelerometer-based, LV and AUSMC SCG signals received from the same places are bandpass filtered at 15 -60 Hz and compared. Both signals from accelerometer and laser vibrometer and signal from the AUSMC device at the same measurement point match well, especially when we consider the temporal waveform (Fig. 4 Top and Middle).

Before mapping the body surface motion, we first characterized the spatial resolution and sensitivity of our system. The spatial resolution was determined using a small disc of 7 mm diameter acting as a point like reflector and located at 720 mm distance. The focusing diameter at -6dB is found to be equal to 30 mm (Fig. 2). This value is equal to 3.5 times the aerial wavelength at 40 kHz for a focal length to aperture width ratio equal to 3. The sensitivity of the AUSMC system is limited by the digital noise of the calculated out of plane velocity images. This noise is of normal distribution with a standard deviation equal to 95µm/s (Fig. 3).

We observed the shape and position of the zone of maximal absolute velocity and acceleration during cardiac cycle (Fig. 5). We also observed a downward-propagating wave on the abdomen. We estimated its velocity by computing the slope of the slanted sections of each systolic wave train on the central-line velocity maps (Fig. 4 Bottom). The mean velocity for the 8 participants was 2.6 ± 0.4 m/s.

DISCUSSION

The study presents preliminary results that support the capability of the AUSMC system to perform SCG measurements in non-contact mode. It is clear, though, that the AUSMC system can potentially go beyond the measurement of SCG timings and amplitudes, by allowing the observation of the motion of the entire thorax and abdomen during heart beats. This multi-channel system can give a greater understanding of the relationship between the biological signals in SCG and the physiological phenomenon creating them.

The slanted sections of each systolic period, a wave propagation phenomenon of ~2.6± 0.4 m/s that is observed in all participants on the abdomen, is much too slow to be the pulse wave of the descending aorta, which travels at ~10 m/s in a healthy normotensive human [START_REF] Mancia | Guidelines for the management of arterial hypertension: The Task Force for the Management of Arterial Hypertension of the European Society of Hypertension (ESH) and of the European Society of Cardiology (ESC)[END_REF]. We have been unable to identify an existing phenomenon observed by clinicians propagating at similar speeds. As such, it presents an intriguing, completely new observation that may have been made possible only with the development of the techniques described in this paper. This wave could potentially be caused by the heart pulling and pushing on the diaphragm, which then propagates in the general viscera, and the velocity could then be related to fat distribution, muscle tone, or other physiologically relevant variables. Studies examining this specific observation have been designed and started by our team.

Results are consistent with the hypothesis that the sternum moves as a solid object under heart forces. These observations are also compatible with the idea that the thoracic cage reacts to cardiac forces as a rather stiff structure, whose movements could be adequately represented with some relatively small set of vibration modes. This is of much value for model validation purposes. Furthermore, distributed measurements enhance identifiability, and therefore we expect more comprehensive and accurate characterization of model properties.

Obtained 2D maps can help for development of new version of the AUSMC system more adapted to the SCG needs and for creation and optimization multisensory accelerometer systems, for example for wearable SCG.

CONCLUSION

A new method was developed for non-contact surface-vibrometry by combining multielement airborne ultrasound arrays, a synthetic aperture implementation and pulsed-waves. Experimental results show the ability of the developed method to obtain 2D surface motion maps with frame rate of several hundred Hz, however, the principle of operation makes it possible to achieve 1 kHz of frame rates using ultrasound waves with the central frequency of 40 kHz.

As a valid method for measuring SCG timings, the AUSMC system could provide a quicker, less invasive alternative to existing accelerometer-based techniques. In a routine visit to the doctor, for example, a device embedded in the wall could perform measurements while the patient is standing on the weighing scale, minimizing effort and time. Also seismocardiographic data collection from multiple locations through multiple channels may help increase the reliability of the method. The AUSMC system has the potential to unlock a new dimension in the study of SCG, although more work is required to establish it firmly. Thus the ultrasound non-contact multichannel SCG measurement system can be a promising technology for 2D cardiac mechanical signal recording, and, in conjunction with ECG, which records electrical cardiac signals, may become an important diagnostic tool in the clinical diagnosis of cardiac diseases, such as valvular heart diseases, heart failure and ventricular dyssynchrony.

Captions for Figures

FIGURE 1 .

 1 FIGURE 1. a) Experimental setup. Data acquisition is performed on sitting participants. b) Schematic top view of the emission/reception panel. Each transmission channel among 32 is used to drive 9 transducers distributed on nine separate sub-panels. The diagram shows that only one transducer is effective, the reflected acoustic waves from the other eight are not correctly directed to the central panel containing the receivers.

FIGURE 2 .

 2 FIGURE 2. Acoustic image of a small disc of 7 mm diameter realized with the AUSMC system.

FIGURE 3 .

 3 FIGURE 3. Experimental noise density function (blue bar) and normal distribution fitting function (red curve) with a standard deviation equal to 95µm/s.

FIGURE 4 .

 4 FIGURE 4. Accelerometer and vibrometer signals compared to the AUSMC signal sampled at points of placement of the sensors (top and middle figures). Mean systolic time: 0.26±0.03 s, mean diastolic time 0.60±0.14 s. B-SCAN along the central line passing through the sternum obtained from the signals of the AUSMC system (bottom figure). Velocity is counted positively when pointing outward of the thorax. The velocity of the downward-propagating surface wave on the abdomen is computed from the slope of the slanted regions.

FIGURE 5 .

 5 FIGURE 5. Body surface motion mapping of the mechanical cardiac activity. Top left image: picture of test subject with highlighted region of interest. Other images: pictures of out of plane surface motion velocity. Note: propagation of surface waves is observed on the body surface during the systole and diastole phases of cardiac cycle. The color scale represents the value of the out of plane surface motion velocity. MC -mitral valve closure, AC -aortic valve closure.

ACKNOWLEDGEMENTS

This work was supported by Chaire Safran -Fonds ESPCI Georges Charpak.


A Biancalani

I Chavdarovski

Z Qiu

A Bottino

D Del Sarto

A Ghizzo

Ö Gürcan

P Morel

I Novikau

Saturation of energetic-particle-driven geodesic acoustic modes due to wave-particle nonlinearity

The nonlinear dynamics of energetic-particle (EP) driven geodesic acoustic modes (EGAM) is investigated here. A numerical analysis with the global gyrokinetic particlein-cell code ORB5 is performed, and the results are interpreted with the analytical theory, in close comparison with the theory of the beam-plasma instability. Only axisymmetric modes are considered, with a nonlinear dynamics determined by waveparticle interaction. Quadratic scalings of the saturated electric field with respect to the linear growth rate are found for the case of interest. As a main result, the formula for the saturation level is provided. Near the saturation, we observe a transition from adiabatic to non-adiabatic dynamics, i.e., the frequency chirping rate becomes comparable to the resonant EP bounce frequency. The numerical analysis is performed here with electrostatic simulations with circular flux surfaces, and kinetic effects of the electrons are neglected.

Introduction

Two main issues related to magnetic confinement fusion are the turbulent transport, and the dynamics of energetic particles (EP), produced by fusion reactions, or injected for heating purposes. Zonal (i.e. axisymmetric) electric fields are observed to interact with turbulence in tokamaks, in the form of zero-frequency zonal flows (ZF) [START_REF] Hasegawa | Nonlinear behavior and turbulence spectra of drift waves and Rossby waves[END_REF]; [START_REF] Rosenbluth | Poloidal Flow Driven by Ion-Temperature-Gradient Turbulence in Tokamaks[END_REF]; [START_REF] Diamond | TOPICAL REVIEW: Zonal flows in plasma -a review[END_REF] and finite-frequency geodesic acoustic modes (GAM) [START_REF] Winsor | Geodesic Acoustic Waves in Hydromagnetic Systems[END_REF]; [START_REF] Zonca | Radial structures and nonlinear excitation of geodesic acoustic modes[END_REF]. Geodesic acoustic modes, due to their finite frequency, can also interact with EP via inverse Landau damping, leading to EP-driven GAM (EGAM) given that EP drive is sufficient to overcome the threshold condition induced by GAM Landau damping and continuum damping [START_REF] Fu | Energetic-Particle-Induced Geodesic Acoustic Mode[END_REF]; [START_REF] Qiu | Nonlocal theory of energetic-particle-induced geodesic acoustic mode[END_REF][START_REF] Qiu | Kinetic theories of geodesic acoustic modes: Radial structure, linear excitation by energetic particles and nonlinear saturation[END_REF][START_REF] Qiu | Geodesic acoustic mode excitation by a spatially broad energetic particle beam[END_REF]; [START_REF] Zarzoso | Fully kinetic description of the linear excitation and nonlinear saturation of fast-ion-driven geodesic acoustic mode instability[END_REF]; [START_REF] Wang | Hole-Clump Pair Creation in the Evolution of Energetic-Particle-Driven Geodesic Acoustic Modes[END_REF]; [START_REF] Girardo | Relation between energetic and standard geodesic acoustic modes[END_REF]; [START_REF] Miki | Finite-Orbit-Width Effects on Energetic-Particle-Induced Geodesic Acoustic Mode[END_REF]; [START_REF] Wang | Simulation study of high-frequency energetic particle driven geodesic acoustic mode[END_REF]; [START_REF] Qu | Energetic Geodesic Acoustic Modes Associated with Two-Stream-like Instabilities in Tokamak Plasmas[END_REF]; [START_REF] Sasaki | A branch of energetic-particle driven geodesic acoustic modes due to magnetic drift resonance[END_REF][START_REF] Sasaki | Toroidal momentum channeling of geodesic acoustic modes driven by fast ions[END_REF]. The mechanism of EGAM excitation can be understood with the analogy to the well known beam plasma instability, where a supra-thermal electron beam interacts with a 1D plasma in a strong axial magnetic field, and the mode is driven unstable by inverse Landau damping [START_REF] Qiu | Kinetic theories of geodesic acoustic modes: Radial structure, linear excitation by energetic particles and nonlinear saturation[END_REF]; [START_REF] Berk | Fast excitation of EGAM by NBI[END_REF]. The similarity lies in the low frequency (ω Ω i , with Ω i being the cyclotron frequency) and toroidally symmetric mode structure of the GAM, which leads to the conservation of µ and P φ , i.e. respectively the magnetic momentum and the toroidal angular momentum. As a result, the EGAM dynamics is quasi-1D in (J, θ) space (with J being here the action conjugate to the poloidal angle θ). This quasi-1D feature of the EGAM determines not only the linear excitation, but also the nonlinear dynamics. Understanding the dynamics of EGAMs is crucial due to their interaction with turbulence, which can modify the turbulent transport [START_REF] Zarzoso | Impact of Energetic-Particle-Driven Geodesic Acoustic Modes on Turbulence[END_REF][START_REF] Zarzoso | Nonlinear interaction between energetic particles and turbulence in gyro-kinetic simulations and impact on turbulence properties[END_REF]. Moreover, a strong nonlinear interaction of EGAM with EP is observed in tokamaks [START_REF] Lauber | ITPA technical meeting on energetic particles[END_REF]; Horváth et al. (2016), leading potentially to a strong redistribution of the EP in phase space.

In order to predict the importance of the interaction of EGAMs with turbulence and with EP, it is important to understand their nonlinear saturation mechanisms. One of the reasons of saturation for EGAMs is the wave-particle nonlinear interaction. In this case, the saturation occurs due to the EP nonlinear redistribution in phase space, and the consequent decrease of the energy exchange between the EP and the EGAM. The wave-particle nonlinear dynamics of EGAMs has been investigated analytically in the framework of the gyrokinetic (GK) theory in [START_REF] Qiu | Kinetic theories of geodesic acoustic modes: Radial structure, linear excitation by energetic particles and nonlinear saturation[END_REF], where the analogy with the beam-plasma instability has been pointed out, and numerically with the hybrid MHD-GK code MEGA in [START_REF] Wang | Hole-Clump Pair Creation in the Evolution of Energetic-Particle-Driven Geodesic Acoustic Modes[END_REF], where in particular the frequency chirping due to the creation of structures in the EP phase space has been studied. Another possible reason of EGAM saturation is the wave-wave coupling. This can occur for an EGAM interacting with another EGAM and generating zonal side-bands, or for an EGAM with non-zonal modes, e.g., turbulence. The wave-wave nonlinear dynamics of EGAMs has been studied analytically in Ref. [START_REF] Qiu | On zero frequency zonal flow and second harmonic generation by finite amplitude energetic particle induced geodesic acoustic mode[END_REF], focusing in particular on the EGAM self-coupling. Numerical investigations of the dynamics of EGAMs with wave-particle and wave-wave nonlinearities have been carried out with the hybrid MHD-GK code MEGA Ido et al. (2016) , with a reduced 1D model solved with the code COBBLES Lesur et al. (2016) and with the GK semi-lagrangian code GYSELA Zarzoso et al. (2012) in the absence of turbulence, and with GYSELA in the presence of turbulence [START_REF] Zarzoso | Impact of Energetic-Particle-Driven Geodesic Acoustic Modes on Turbulence[END_REF]. Although much has been understood on the linear and nonlinear dynamics of EGAMs, thanks to the above-mentioned studies, nevertheless, to the present day no theory exists to predict the EGAM saturation level. This study, in the regime of a dominant waveparticle nonlinearity, is the goal of this paper.

The nonlinear saturation of EGAMs is investigated here by means of electrostatic simulations with the gyrokinetic particle-in-cell code ORB5 Jolliet et al. (2007); [START_REF] Bottino | Global simulations of tokamak microturbulence: finite-β effects and collisions[END_REF]; [START_REF] Bottino | Monte Carlo particle-in-cell methods for the simulation of the Vlasov-Maxwell gyrokinetic equations[END_REF]. ORB5 has been succesfully verified against analytical theory and benchmarked against other gyrokinetic codes, for the linear dynamics of GAMs [START_REF] Biancalani | Cross-code gyrokinetic verification and benchmark on the linear collisionless dynamics of the geodesic acoustic mode[END_REF]EGAMs Biancalani et al. (2014); [START_REF] Zarzoso | Analytic dispersion relation of energetic particle driven geodesic acoustic modes and simulations with NEMORB[END_REF]. A detailed comparison with the beam-plasma instability (BPI) [START_REF] O'neil | Collisionless Damping of Nonlinear Plasma Oscillations[END_REF][START_REF] O'neil | Transition of the Dispersion Roots from Beam-Type to Landau-Type Solutions[END_REF] in a 1-dimensional uniform plasma is also done, following the scheme anticipated in Ref. [START_REF] Qiu | Kinetic theories of geodesic acoustic modes: Radial structure, linear excitation by energetic particles and nonlinear saturation[END_REF][START_REF] Qiu | Excitation of Kinetic Geodesic Acoustic Modes by Drift Waves in Nonuniform Plasmas[END_REF]. Similarly to the BPI, the saturation level of EGAM is shown to scale quadratically with the linear growth rate. A similar investigation was also previously done for Alfvén modes (see, for example, Ref. [START_REF] Briguglio | Analysis of the nonlinear behavior of shear-Alfvén modes in tokamaks based on Hamiltonian mapping techniques[END_REF]). As a main result of these considerations, we provide here the theoretical prescription for the quantitative prediction of the EGAM saturation level. The EGAM frequency is also studied, and shown to evolve in time when approaching the saturation, like the BPI (see, for example, Ref. [START_REF] Morales | Nonlinear Frequency Shift of an Electron Plasma Wave[END_REF]; [START_REF] Armon | Chirped resonance dynamics in phase space[END_REF]). The chirping rate is observed to get of the order of magnitude of the squared of the EP bounce frequency near the saturation. This denotes a transition from adiabatic to non-adiabatic regimes. The importance of understanding whether the EGAM nonlinear dynamics is adiabatic or non-adiabatic is linked to the underlying mechanims of the saturation, i.e., whether the saturation is accessed through the local flattening of particle distribution in the resonance region due to wave-particle trapping [START_REF] Qiu | Kinetic theories of geodesic acoustic modes: Radial structure, linear excitation by energetic particles and nonlinear saturation[END_REF] or through the avalanches of phase space mode structures [START_REF] Zonca | Nonlinear dynamics of phase space zonal structures and energetic particle physics in fusion plasmas[END_REF] due to the feedbacks of EP distribution and mode structure in phase space.

The paper is organized as follows. The adopted gyrokinetic model is described in Sec. 2, and the equilibrium and case definition in Sec. 3. The linear dynamics is described in Sec. 4. The saturation levels are investigated in Sec. 5. The regimes of different adiabaticity are investigated by means of the analysis of the frequency, in Sec. 6. Finally, a summary of the results is given in Sec. 7.

The model

The main damping mechanism of GAM and EGAM is the Landau damping, which makes the use of a kinetic model necessary. In this paper we use the global gyrokinetic particle-in-cell code ORB5 Jolliet et al. (2007). ORB5 was originally developed for electrostatic ITG turbulence studies, and recently extended to its multi-species electromagnetic version in the framework of the NEMORB project [START_REF] Bottino | Global simulations of tokamak microturbulence: finite-β effects and collisions[END_REF]; [START_REF] Bottino | Monte Carlo particle-in-cell methods for the simulation of the Vlasov-Maxwell gyrokinetic equations[END_REF]. In this paper, collisionless electrostatic simulations are considered.

The model equations of the electrostatic version of ORB5 is made by the trajectories of the markers, and by the gyrokinetic Poisson law for the scalar potential φ. These equations are derived in a Lagrangian formulation [START_REF] Bottino | Monte Carlo particle-in-cell methods for the simulation of the Vlasov-Maxwell gyrokinetic equations[END_REF], based on the gyrokinetic Vlasov-Maxwell equations of Sugama, Brizard and Hahm Sugama (2000); [START_REF] Brizard | Foundations of nonlinear gyrokinetic theory[END_REF]. The equations for the marker trajectories (in the electrostatic version of the code) are [START_REF] Bottino | Monte Carlo particle-in-cell methods for the simulation of the Vlasov-Maxwell gyrokinetic equations[END_REF]:

Ṙ = 1 m s p B * B * + c q s B * b × µ∇B + q s ∇ φ (2.1) ṗ = - B * B * • µ∇B + q s ∇ φ (2.2) μ = 0 (2.3)
The coordinates used for the phase space are (R, p , µ), i.e. respectively the gyrocenter position, canonical parallel momentum p = m s v and magnetic momentum µ = m s v 2 ⊥ /(2B) (with m s and q s being the mass and charge of the species). v and v ⊥ are respectively the parallel and perpendicular component of the particle velocity. The gyroaverage operator is labeled here by the tilde symbol ˜. B * = B + (c/q s )∇ × (b p ), where B and b are the equilibrium magnetic field and magnetic unitary vector.

Kinetic effects of the electrons are neglected. This is done by calculating the electron gyrocenter density directly from the value of the scalar potential as [START_REF] Bottino | Monte Carlo particle-in-cell methods for the simulation of the Vlasov-Maxwell gyrokinetic equations[END_REF]:

n e (R, t) = n e0 + q s n e0 T e0 φ -φ (2.4)
where φ is the flux-surface averaged potential, instead of treating the electrons with markers evolved with Eqs. 2.1, 2.2, 2.3. We focus on the dynamics of zonal perturbations, by filtering out all non-zonal components (this is to avoid interactions of zonal/non-zonal modes). Wave-wave coupling is neglected, by evolving the bulk-ion and electron markers along unperturbed trajectories. This means that, in Eqs. 2.1, 2.2, 2.3 for the bulk ions, the last terms, proportional to the EGAM electric field, are dropped. The nonlinear wave-particle dynamics is studied by evolving the EP markers along the trajectories which include perturbed terms associated with the EGAM electric field. This means that the EP markers are evolved with Eqs. 2.1, 2.2, 2.3 where the terms proportional to the EGAM electric field are retained.

The gyrokinetic Poisson equation is [START_REF] Bottino | Monte Carlo particle-in-cell methods for the simulation of the Vlasov-Maxwell gyrokinetic equations[END_REF]:

-∇ • n 0 m i c 2 B 2 ∇ ⊥ φ = s =e dW q s δ f s -n e (R, t) (2.5)
with n 0 m i being here the total plasma mass density (approximated as the ion mass density). The summation over the species is performed for the bulk ions and for the EP, whereas the electron contribution is given by -n e (R, t). Here δf = f -f 0 is the gyrocenter perturbed distribution function, with f and f 0 being the total and equilibrium (i.e. independent of time, assumed here to be a Maxwellian) gyrocenter distribution functions. The integrals are over the phase space volume, with dW = (2π/m 2 i )B * dp dµ being the velocity-space infinitesimal volume element.

Equilibrium and simulation parameters

The tokamak magnetic equilibrium is defined by a major and minor radii of R 0 = 1 m and a = 0.3125 m, a magnetic field on axis of B 0 = 1.9 T, a flat safety factor radial profile, with q = 2, and circular flux surfaces (with no Grad-Shafranov shift). Flat temperature and density profiles are considered at the equilibrium for all species. The bulk plasma temperature is defined by ρ * = ρ s /a, with ρ s = c s /Ω i , with c s = T e /m i being the sound speed. Three increasing values of bulk plasma temperature are considered, for investigating the dependence of our results on the Landau damping, corresponding to: ρ * 1 = 1/256 = 0.0039, ρ * 2 = 1/128 = 0.0078, and ρ * 3 = 1/64 = 0.0156 (τ e = T e /T i = 1 for all cases described in this paper).

All parameters defined so far, are adopted by ORB5 with no consideration of the mass of the bulk ion species. The choice of the bulk ion mass is done only during the post-processing of the results of ORB5, if we need to know the values of equilibrium or perturbed quantities in non-normalized units. In particular, in the case of a hydrogen plasma, we get a value of the ion cyclotron frequency of Ω i = 1.82 • 10 8 rad/s. With this choice of bulk specie, we can calculate the plasma temperature and sound frequency. The three values of plasma temperature are T i1 = 515 eV, T i2 = 2060 eV, and T i3 = 8240 eV. The sound frequency is defined as

ω s = 2 1/2 v ti /R (with v ti = T i /m i , which
for τ e = 1 reads v ti = c s ). We obtain the following three values of the sound velocity: c s1 = 2.22 • 10 5 m/s, c s2 = 4.44 • 10 5 m/s, c s1 = 8.88 • 10 5 m/s. These correspond to the following three values of the sound frequency: ω s1 = 3.14 • 10 5 rad/s, ω s2 = 6.28 • 10 5 rad/s, and ω s3 = 1.25 • 10 6 rad/s.

The energetic particle distribution function is a double bump-on-tail, with two bumps at v = ±v bump (see Fig. 1), like in Ref. [START_REF] Biancalani | Numerical validation of the electromagnetic gyrokinetic code NEMORB on global axisymmetric modes[END_REF]. In this paper, v bump = 4 v ti is chosen. In order to initialize a distribution function which is function of constants of motion only, the modified variable ṽ

= 2(E -µB max )/m/v ti is used instead of v = 2(E -µB(r, θ))/m/v ti (

Linear dynamics

The linear dynamics of EGAMs in an equilibrium similar to the one considered in Sec. 3, has been investigated in Ref. [START_REF] Zarzoso | Analytic dispersion relation of energetic particle driven geodesic acoustic modes and simulations with NEMORB[END_REF] for a bulk plasma temperature given by ρ * = 1/64 = 0.0156, by means of ORB5 simulations and analytical theory. A scan on the EP concentration is performed here, similarly to what was done in Ref. [START_REF] Zarzoso | Analytic dispersion relation of energetic particle driven geodesic acoustic modes and simulations with NEMORB[END_REF], but with three different values of bulk plasma temperature, corresponding to three different values of ρ * , as described in Sec. 3. The dependence of the linear dynamics (frequency and growth rate) on the EP concentration is shown in Fig. 2. Both the frequency and the growth rates are observed to follow the qualitative scalings as described in Ref. [START_REF] Biancalani | Numerical validation of the electromagnetic gyrokinetic code NEMORB on global axisymmetric modes[END_REF] and [START_REF] Zarzoso | Analytic dispersion relation of energetic particle driven geodesic acoustic modes and simulations with NEMORB[END_REF]. Note, in particular, that the growth rate does not grow linearly with n EP .

The dependence of the frequency on the EP concentration is not observed to change with ρ * . Regarding the growth rate, no change is observed when going from ρ * 3 = 0.0039 to ρ * 2 = 0.0078, meaning that the measured growth rate is basically given by the value of the drive, and the Landau damping here is negligible for the chosen values of EP concentration. On the other hand, when further increasing the temperature, and going to ρ * 1 = 0.0156, a smaller value of growth rate is measured, meaning a higher Landau damping. The transit resonance velocity of the EP can be calculated by knowing the EGAM frequency of a specific simulation. Considering a case with n EP /n i = 0.12 as an example, the frequency is measured as: ω L = 1.2 ω s . For comparison, the GAM frequency for these parameters is ω GAM = 1.8 ω s . Then, the transit resonance velocity in the linear phase is calculated as v 0 = qR ω L = 3.4 v ti , with ω L being the EGAM linear frequency.

Scaling of the saturated amplitudes

In this Section, we focus on the value of the saturated electric field δ E r , and we investigate its dependence on the value of the linear growth rate and of the damping. The corresponding scaling of δ E r sheds light on the mechanism which is responsible for the saturation. A one-to-one comparison with the saturation mechanism of the BPI is also described.

The amplitude of the EGAM at saturation has been measured in different simulations performed with ORB5. The wave-particle nonlinearity only has been considered, by pushing on perturbed trajectories the EP population only. Therefore, even if the linear resonance velocity falls near the tail of the bulk ions (see Fig. 1), no nonlinear interaction of those bulk ions with the EGAM is considered in our model. More in general, the nonlinear interaction of the bulk ions with the EGAM is expected to be important: this has been studied in part analytically in Ref. [START_REF] Qiu | On zero frequency zonal flow and second harmonic generation by finite amplitude energetic particle induced geodesic acoustic mode[END_REF], and its study with ORB5 is in progress. Different values of the bulk temperature, and the EP concentration, have been considered. As an example, the radial structure of a nonlinear simulation with ρ * = 0.0156, n EP /n i = 0.176 is depicted in Fig. 3-a. No sensible change in the radial wave-number is observed when going from the linear phase, to the saturation, and after the saturation. This confirms that in this particular configuration, where all equilibrium radial profiles are flat, EGAM can be treated as a 1-dimensional problem where the radial direction does not play an important role.

When varying only the bulk temperature, both the linear frequency and growth rate are observed to scale with the sound frequency, which is a good normalization frequency (consistently with Fig. 2). The saturation level increases with the linear growth rate, similarly to other instabilities like the BPI in a uniform system and the Alfvén instabilities in tokamaks. This is depicted in Fig. 3-b, where non-normalized units are used (in particular, the ion cyclotron frequency is selected as a time unit not depending on the temperature).

The scalings with the energetic particle concentration are also investigated. The results are shown in Fig. 4a. We obtain that, in the considered regime, the saturated level scales as the quadratic power of the linear growth rate. This quadratic scaling is typical for marginally stable bump-on-tail instabilities, as derived by O'Neil O'Neil (1965); O 'Neil & Malmberg (1968).

We can consider the problem to be similar to a monochromatic beam-plasma system, in which particles are moving in the potential well of the perturbed electric field. Depending on the energy, some particles are trapped inside the well and execute bounce motion with frequency ω b in the frame moving with the wave phase velocity. The resonant particles exchange energy with the mode, causing the amplitude to grow and the particles to redistribute in phase space, flattening the velocity distribution in the vicinity of the resonant parallel velocity v = ωqR 0 . The drive is due to the positive slope of the particle distribution in the velocity space at the resonant parallel velocity, which acts as an inverse Landau damping. In the initial stage ω b γ L , the mode grows exponentially with a linear growth rate, making more and more particles to become trapped in the phase space. After some significant particle velocity redistribution, the power exchange between the wave and the particles is balanced, causing the wave amplitude to saturate. Since the initial perturbation is negligible, the saturation level is determined by the exchange of energy between the mode and a band of resonant particles [START_REF] Levin | Contribution to the Nonlinear Theory of Kinetic Instability of an Electron Beam in Plasma[END_REF]. The chirping of the mode seen in Fig. 5 is a strongly non-linear effect that occurs when the amplitude is large enough to have the trapped (and more generally resonant) particles with ω b ∼ γ L drastically change the dynamics of the mode through the modification of the distribution function and non-perturbative fast particle response. Namely, the mode dynamics is determined by all resonant particles that exhibit a continuous oscillation, trapping and detrapping in the potential well of the mode, thus contributing (nonperturbatively) to the non-adiabatic behavior observed in the simulation.

The quadratic scaling of the saturation level with the damping obtained with our simulations is similar to the saturation of the BPI, where it occurs due to wave-particle trapping [START_REF] O'neil | Collisionless Damping of Nonlinear Plasma Oscillations[END_REF]. For the BPI, the original reference of M. B. [START_REF] Levin | Contribution to the Nonlinear Theory of Kinetic Instability of an Electron Beam in Plasma[END_REF] gives a value of ω b = 3.06 γ L at saturation [START_REF] Levin | Contribution to the Nonlinear Theory of Kinetic Instability of an Electron Beam in Plasma[END_REF] (for comparison, note that more recent numerical calculations find ω b = 3.2 γ L [START_REF] Lesur | Fully nonlinear features of the energetic beamdriven instability[END_REF]; Carlevaro et al. (2016a,b)). For EGAMs, the bounce frequency is given by [START_REF] Qiu | Kinetic theories of geodesic acoustic modes: Radial structure, linear excitation by energetic particles and nonlinear saturation[END_REF]:

ω 2 b = α 1 δ E r , with α 1 ≡ e Vdc 2m EP v 0 qR 0 (5.1)
with m EP being the mass of the energetic particle specie, considered equal to the bulk ion mass in this paper, v 0 the velocity matching the resonance condition, and Vdc = m EP v 2 0 /(eBR) the magnetic curvature drift. Therefore we have:

α 1 = v 0 2qR 2 B = ω L 2RB (5.2)
We emphasize that the value of α 1 depends on ω L . This is a main difference with respect to the BPI, where there is only one value ω lin = ω pe , with ω pe being the plasma frequency.

The dependence of the maximum electric field on the linear growth rate can be measured with the results of the numerical simulations. For the simulations shown in Fig. 4, we find: ), the value of β as given in Eq. 5.5, vs the linear frequency, for the same simulations. The black dashed line is the sqruare root interpolation. For a reference, the black star shows the result obtained for the BPI in Ref. [START_REF] Levin | Contribution to the Nonlinear Theory of Kinetic Instability of an Electron Beam in Plasma[END_REF].

δ E r = α 2 γ 2 L (5.3) 10 -1 10 4 10 5 γ L [2 1/2 v ti / R] δ E r [V/m]
The values of α 2 are found to depend on the bulk temperature. For the three chosen increasing values of ρ * , i.e. ρ * = 0.0039, 0.0078 and 0.0156, we have respectively α 2 = 0.47 • 10 7 V/m, 0.9 • 10 7 V/m, and 2.0 • 10 7 V/m. Finally the relationship between the EP bounce frequency and the linear growth rate is obtained:

ω b = β γ L (5.4)
where β is calculated as β = (α 1 α 2 ) 1/2 /ω s , which yields:

β = β 0 ω L ω GAM 1/2 , with β 0 = 1 ω s ω GAM α 2 2RB 1/2
(5.5)

Note that here, β depends on the EGAM frequency, which changes with the intensity of the drive, given here by the EP concentration (see Fig. 2). As a comparison, note that in the problem of the BPI, solved in Ref. [START_REF] O'neil | Collisionless Damping of Nonlinear Plasma Oscillations[END_REF]; O'Neil & Malmberg (1968); [START_REF] Levin | Contribution to the Nonlinear Theory of Kinetic Instability of an Electron Beam in Plasma[END_REF], on the other hand, the mode frequency is assumed to be constant and equal to the frequency measured in the absence of EP. On the other hand, the values of β are found not to depend on the damping, which changes with the three different bulk temperatures: an interpolation can be drawn for all considered simulations, and shown to depend on ω L only (see Fig. 4b). In this sense, the formula given in Eq. 5.5 is universal for the chosen regime, because it does not depend on the bulk plasma temperature.

The considered equilibrium has been chosen in order to excite EGAMs out of a GAM, and not out of a Landau pole, as described in Ref. [START_REF] Zarzoso | Analytic dispersion relation of energetic particle driven geodesic acoustic modes and simulations with NEMORB[END_REF]. This choice of the mode is done, in order to have a one-to-one correspondence with the BPI, where the mode which is excited by the energetic electron beam is the Langmuir wave which is an eigenmode of the system in the absence of EP. Following this consideration, we can consider the interpolation of the results shown in Fig. 4, and take the extrapolation to ω L → ω GAM , which is the limit assumed in the resolution of the BPI. In this case, the extrapolation gives a unique value, which defines the EGAM instability, i.e. β 0 = β(ω L /ω GAM = 1) = 2.66. This is to be compared with the value of β obtained for the BPI [START_REF] Levin | Contribution to the Nonlinear Theory of Kinetic Instability of an Electron Beam in Plasma[END_REF]; [START_REF] Lesur | Fully nonlinear features of the energetic beamdriven instability[END_REF]; Carlevaro et al. (2016a,b), i.e. β BP I = 3.2 (originally estimated as β BP I = 3.06 by Levin).

Finally, by using Eqs. 5.1, 5.2, 5.4 and 5.5, we can write the formula for the saturated electric field as a function of the linear characteristics of the mode:

δE r = 2RBβ 2 0 ω GAM γ 2 L (5.6)
with the value of β 0 = 2.66 in the regime considered in this paper.

Frequency

In this section, we show the results of the measurement of the time evolution of the EGAM frequency. In Sec. 5, we have shown that a quadratic scaling of the saturated electric field on the linear growth rate is found. This quadratic scaling, has a one-toone correspondence on the Langmuir wave problem investigated by O'Neil, where the saturation occurs due to wave-particle trapping. The wave-particle trapping mechanism, is usually referred to as adiabatic, meaning that a slowly increasing potential well traps more and more energetic particles. In this adiabatic regime, the mode frequency varies very slowly with respect to the bounce frequency. On the other hand, in the EGAM case considered here, we show that the saturation is not strictly adiabatic, but a transition between adiabatic and non-adiabatic regime occurs at the time of the saturation.

In this section, we take again the EGAM case with ρ * = 0.0078, n EP /n i = 0.12, as a typical case, and we investigate the variation of the frequency in time in comparison with the bounce frequency. For the measurement of the frequency, we use the radial zonal electric field measured at s=0.5. The measurement of the frequency is performed in two different ways: a) as an average of the period between several EGAM oscillation peaks, as shown in Fig. 5-a; b) with a short-time Fourier transform (STFT), as shown in Fig. 5-b.

With the first technique, namely measuring the frequency by inversion of the period between neighbouring peaks, an upward chirping is observed in the nonlinear phase, of the order of 10% of the linear frequency. This means that the resonance condition changes in time, with resonance velocity slightly increasing at the time of the saturation or in the later phase (see Fig. 5-a). A dominantly upward chirping of EGAMs was previously observed and documented in Ref. [START_REF] Berk | Explanation of the JET n = 0 chirping mode[END_REF]; [START_REF] Wang | Hole-Clump Pair Creation in the Evolution of Energetic-Particle-Driven Geodesic Acoustic Modes[END_REF]. In the cases of interest, the frequency changes with a similar time scale as the evolution of the amplitude. In this regime, we refer to this change of the frequency in time as nonlinear frequency shift, i.e. frequency chirping.

The second technique, consists in measuring the frequency with a short-time Fourier transform (STFT) on a Hamming time-window. With this technique, the error bar in frequency is large (due to the small number of oscillations in the nonlinear regime around the saturation), namely of the order of 10-20% (see Fig. 5-b). With such a big error-bar, no clear upward chirping is observed. Near the time of the saturation, i.e. t 2.2 Ω -1 i , only one mode is observed. This is the condition of application of the direct technique of measurement of the frequency described above, where the frequency can be measured as the inverse of the period among peaks.

As mentioned above, the EP bounce frequency ω b depends on the mode amplitude. For this case, the mode amplitude grows during a linear phase, then enters the nonlinear phase and saturates at the time t = 2.3 • 10 4 Ω -1 i . The corresponding value of ω 2 b for the considered simulation is shown in Fig. 6-a. When the EGAM frequency evolves slowly in time with respect to the inverse of the bounce frequency, then the EP can bounce back and forth many times, and this is called adiabatic dynamics. The adiabaticity parameter, defined as ω /ω 2 b , measures the level of adiabaticity of the dynamics. Here, ω is the time derivative of the frequency. The time evolution of the adiabaticiy parameter for the considered simulation is depicted in Fig. 6-b (with ω being the instantaneous mode . Nonlinear evolution of the frequency, measured as a short-time average of the period between the peaks (left, "a") or with a short-time Fourier transform (right, "b"), for ρ * = 0.0078, nEP /ni = 0.12. frequency). A transition from adiabatic to non-adiabatic dynamics near the saturation is observed. In particular, near the saturation, the EP do not have the time to perform many bounces during the nonlinear modification of the wave. From this respect, the EGAM dynamics and the BPI are not in analogy. The details of the EGAM saturation mechanism will be investigated with a power-balance diagnostics in phase space, and discussed in a separate paper, where a comparison with the nonlinear dynamics described in Ref. [START_REF] Berk | Explanation of the JET n = 0 chirping mode[END_REF]; [START_REF] Wang | Hole-Clump Pair Creation in the Evolution of Energetic-Particle-Driven Geodesic Acoustic Modes[END_REF] will be discussed.

Summary and discussion

The importance of understanding the nonlinear dynamics of EGAM, i.e. energetic particle (EP) driven geodesic acoustic modes (GAM), is due to their possible role in interacting with turbulence and with the EP population present in tokamak plasmas as a result of nuclear reactions and/or heating mechanisms. In particular, the level of the nonlinear saturation of EGAM is directly related to their efficiency in regulating turbulence or redistributing EP in phase space.

In this paper, we have investigated the problem of the nonlinear saturation of EGAM with the gyrokinetic particle-in-cell code ORB5, focusing on the wave-particle nonlinearity. Electrostatic collisionless simulations have been considered, with circular flux surfaces magnetic equilibria, and neglecting the kinetic effects of electrons.

The level of the saturated electric field has been shown to scale quadratically with the linear growth rate, similarly to the beam-plasma instability (BPI) in 1D uniform plasmas, and to Alfvén Eigenmodes (AE) near marginal stability. We note that, due to the flat density and temperature profiles assumed for all species at the equilibrium, the EGAM has a strict 1D dynamics similar to the BPI. In fact, the EGAM electric perturbation is low frequency (ω Ω i ) and zonal (m=n=0) and when radial nonuniformities are negligible, like in our case, the direction parallel to the equilibrium magnetic field becomes the only relevant direction for the wave-particle interaction, i.e. for the inverse Landau damping. We also note that, in the case of beta-induced AE, due to the finite radial mode structure compared to resonant particle radial excursion, deviation from the δE r ∝ γ 2 L is observed due to the competition of resonance detuning and radial decoupling [START_REF] Wang | Nonlinear dynamics of beta-induced Alfvén eigenmode driven by energetic particles[END_REF]; [START_REF] Zhang | Nonlinear Frequency Oscillation of Alfvén Eigenmodes in Fusion Plasmas[END_REF]. A similar deviation is anticipated as one further increases the drive of EPs, from the correspondence of EP pitch angle scattering by frequency chirping EGAMs and the radial wave-particle pumping by frequency chirping finite-n BAEs. This will be the content of our next publication.

We have also investigated the relationship of the bounce frequency of the EP in the field of the EGAM, with the linear growth rate, finding a linear dependence, similarly to the BPI. The linear coefficient, in the case of the EGAM problem, is proportional to the square root of the EGAM frequency, and is not a constant like in the beam-plamsa instability. In the limit of ω L /ω GAM → 1, which corresponds to the BPI problem when the mode frequency does not move sensibly from the value of the plasma frequency ω p , a unique value of the linear coefficient can be estimated: β 0 = 2.66, to be compared with the factor β BP I 3.2 of the BPI. This value of β 0 = 2.66 defines the EGAM in the regime of interest. As a main result, the formula for the saturation level has been provided. The investigation of the dependence of β on equilibrium parameters like the safety factor, the magnetic flux surface elongation, the EP distribution function, and the effect of kinetic electrons is left to a future work.

Finally, we have investigated the temporal variation of the EGAM frequency during the saturation phase. The initial adiabatic regime, defined as the regime where the time derivative of the EGAM frequency is much smaller than the squared EP bounce frequency, has been shown to have a transition to a non-adiabatic regime at the saturation. The exact analytical model of this phenomenon is still a work in progress and will be the first expansion of this work. The simplified picture is that the mode will chirp to maintain maximal exchange of power with the energetic particles [START_REF] Chen | Physics of Alfvén waves and energetic particles in burning plasmas[END_REF], while balancing the Landau damping. The mode frequency is at the same time bound by the wave dispersion relation and plasma equilibrium, making their simultaneous examination with the power exchange equation necessary for determining the rate and direction of chirping. The non-adiabatic behavior is enabled by the non-perturbative response of the energetic particles, which consists of significant distortion of the distribution function in the resonant region, as well as of the particle trapping and detrapping [START_REF] Chen | Physics of Alfvén waves and energetic particles in burning plasmas[END_REF]. Differently, the adiabatic chirping is connected with either external driven fluctuations, or slow energetic particle redistribution, provided that the drive is sufficiently weak [START_REF] Breizman | Critical nonlinear phenomena for kinetic instabilities near threshold[END_REF].

  similarly to Ref. Zarzoso et al. (2012); Biancalani et al. (2014); Zarzoso et al. (2014)). Neumann and Dirichlet boundary conditions are imposed to the scalar potential, respectively at the inner and outer boundaries, s = 0 and s = 1.

Figure 1 .

 1 Figure 1. Initial EP distribution function for a simulation with nEP /ni = 0.12, v bump /vti = 4.

Figure 2 .

 2 Figure 2. Frequency (left, "a") and growth rate (right, "b") vs EP concentration, for simulations with ζ = v bump /v th =4.

Figure 3 .

 3 Figure 3. On the left (Fig. "a"), EGAM normalized radial structure for ρ * = 0.0156, nEP /ni = 0.176. On the right (Fig. "b"), absolute value of the electric field measured at the position of the peak, s = 0.6, for three different simulations with respectively ρ * = 0.0039 (blue), ρ * = 0.0078 (red), ρ * = 0.0156 (green). All simulations here have nEP /ni = 0.30. The time is expressed in units of Ω -1 i .

Figure 4 .

 4 Figure 4. On the left (Fig. "a"), maximum value of the EGAM radial electric vs linear growth rate, for the same simulations as in Fig.2. The red, blue and green crosses refer respectively to ρ * = 0.0039, ρ * = 0.0078, ρ * = 0.0156. The dashed lines are the quadratic fitting formulas. On the right (Fig."b"), the value of β as given in Eq. 5.5, vs the linear frequency, for the same simulations. The black dashed line is the sqruare root interpolation. For a reference, the black star shows the result obtained for the BPI in Ref.[START_REF] Levin | Contribution to the Nonlinear Theory of Kinetic Instability of an Electron Beam in Plasma[END_REF].

  Figure5. Nonlinear evolution of the frequency, measured as a short-time average of the period between the peaks (left, "a") or with a short-time Fourier transform (right, "b"), for ρ * = 0.0078, nEP /ni = 0.12.

Figure 6 .

 6 Figure6. Squared bounce frequency (left, "a") and adiabaticity (right, "b") for the EGAM with ρ * = 0.0078, nEP /ni = 0.12.

† Email address for correspondence: www2.ipp.mpg.de/∼biancala

Acknowledgements

Interesting discussions with N. Carlevaro, G. Montani, F. Zonca on the nonlinear dynamics of the energetic particles are acknowledged. Interesting discussions with A. Mishchenko and B. Finney McMillan on the treatment of EGAMs with ORB5, and with the whole ORB5 team, are also acknowledged. Part of this work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053, within the framework of the Nonlinear energetic particle dynamics (NLED) European Enabling Research Project, WP 15-ER-01/ENEA-03 and Nonlinear interaction of Alfvénic and Turbulent fluctuations in burning plasmas (NAT) European Enabling Research Project, CfP-AWP17-ENR-MPG-01. The views and opinions expressed herein do not necessarily reflect those of the European Commission. Part of this work has been funded by the Centre de Coopération Universitaire Franco-Bavarois -Bayerisch-Französisches Hochschulzentrum, for the grant FK 30 15. Simulations were performed on the Marconi supercomputer within the framework of the OrbZONE and ORBFAST projects. Part of this work was done while two of the authors (A. Biancalani and I. Novikau) were visiting LPP-Palaiseau, whose team is acknowledged for the hospitality. One of the authors (I.C.) wants to thank the Max Planck Princeton Center for the support during this work. This paper is dedicated to Prof. Francesco Pegoraro.


G Zheng
email: gang.zheng@inrialpes.fr

D Boutat
email: driss.boutat@

J P Barbot
email: barbot@ensea.fr

SINGLE OUTPUT DEPENDENT OBSERVABILITY NORMAL FORM

This paper gives the sufficient and necessary conditions which guarantee the existence of a diffeomorphism in order to transform a nonlinear system without inputs into a canonical normal form depending on its output. Moreover we extend our results to a class of systems with inputs.

1. Introduction. Since Luenberger's work [START_REF] Luenberger | An introduction to observers[END_REF], the design of an observer for observable linear systems with linear outputs has been a well-known concept. In order to use the same observer for nonlinear systems, the so-called observability linearization problem for nonlinear systems was born. The sufficient and necessary conditions which guarantee the existence of a diffeomorphism and of an output injection to transform a single output nonlinear system without inputs into a linear one with an output injection were firstly addressed in [START_REF] Krener | Linearization by output injection and nonlinear observer[END_REF]. Then, for a multi-output nonlinear system without inputs, the linearization problem was partially solved in [START_REF] Krener | Nonlinear observer with linearizable error dynamics in[END_REF]. The complete solution to the linearization problem was given in [START_REF] Xia | Nonlinear observer design by observer error linearization in[END_REF]. Another approach was introduced for the analytical systems in [START_REF] Kazantzis | Nonlinear observer design using Lyapunov's auxiliary theorem[END_REF] by assuming that the spectrum of the linear part must lie in the Poincaré domain and it was generalized in [START_REF] Krener | Nonlinear observer design in the Siegel domain through coordinate changes[END_REF] by assuming that the spectrum of the linear part must lie in the Siegel domain. These assumptions are not in generally fulfilled. Other approaches using quadratic normal forms were given in [START_REF] Boutat-Baddas | Quadratic Observability normal form[END_REF] and [START_REF] Chabraoui | Observability quadratic characteristic numbers[END_REF]. All these approaches enable us to design an observer for a larger class of nonlinear systems.

Meanwhile, other researchers worked on designing nonlinear observers directly, such as high-gain observers [START_REF] Hammouri | Bilinearization up to the output injection[END_REF], [START_REF] Gauthier | Deterministic observation theory and applications[END_REF], [START_REF] Hammouri | Observer synthesis for state-affine systems[END_REF]. Nevertheless, even if the conditions which guarantee the linearization method to design an observer were not generically fulfilled, this method still remained important for the nonlinear observer design first because it works well for non-analytic systems, and second because it could be used in the adaptive theory and also be useful for the observation of systems with unknown inputs. All these reasons explain why researchers carry on investigating this matter.

In [START_REF] Respondek | Time scaling for observer design with linearization error dynamics[END_REF], the author gave the sufficient and necessary geometrical conditions to transform a nonlinear system into a so-called output-dependent time scaling linear canonical form. While [START_REF] Guay | Observer linearization by output diffeomorphism and output-dependent time-scale transformations[END_REF] gave the dual geometrical conditions of [START_REF] Respondek | Time scaling for observer design with linearization error dynamics[END_REF].

In this paper, as an extension of [START_REF] Zheng | Output Dependent Observability Linear Normal Form[END_REF], we will propose a method to deduce the geometrical conditions which are sufficient and necessary to guarantee the existence of a local diffeomorphism z = φ (x) which transforms the locally observable dynamical system

ẋ = f (x), y = h(x), (1.1)
where

x ∈ U ⊂ IR n , f : U ⊂ IR n → IR n and h : U ⊂ IR n → IR are sufficiently smooth, into the following form ż = A(y)z + β(y), y = z n = Cz, (1.2)
where

A(y) =        0 • • • 0 0 0 α 1 (y) • • • 0 0 0 . . . . . . . . . • • • . . . 0 • • • α n-2 (y) 0 0 0 • • • 0 α n-1 (y) 0        , β(y) =        β 1 (y) β 2 (y) . . . β n-1 (y) β n (y)       
, and α i (y) = 0 for y ∈ ]-a, a[ and a > 0. This kind of linearization is called Single Output Dependent Observability normal form (SODO normal form). For dynamical systems in the form (1.2), we may for example apply the following high gain observer [START_REF] Busawon | A simple observer for a class of nonlinear systems[END_REF]:

ż = A(y)ẑ + β(y) -Γ -1 (y) R -1 ρ C T (C ẑ -y), 0 = -ρR ρ -ĀT R ρ -R ρ Ā + C T C. (1.3) where Γ (y) is the n × n diagonal matrix Γ (y) = diag n-1 i=1 α i (y), n-1 i=2 α i (y), • • • , α n-1 (y), 1 and Ā is the n × n matrix defined as follows Ā =      0 • • • 0 0 1 • • • 0 0 . . . . . . . . . . . . 0 • • • 1 0      .
Indeed, here the output of system (1.2) is considered as an input of (1.3). Setting e = z -ẑ, the observation error can be obtained as follows:

ė = A (y) -Γ -1 (y) R -1 ρ C T C e,
And the convergence of such observer is proved in [START_REF] Busawon | A simple observer for a class of nonlinear systems[END_REF], thus in section 4 we simply highlight the design of such observer for systems in the form (1.2). Moreover, we generalize our result to a class of systems with inputs. Then, some useful corollaries are discussed in order to deal with affine systems and the so-called left invertibility problem.

This paper is organized as follows. The next section addresses notations and technical results which are key to prove our main result. In section 3, we present our method to deduce the geometrical conditions for a nonlinear system without inputs in order to transform it into a SODO normal form. Section 4 is devoted to the generalization of our results to a class of systems with inputs. In the same section, some practical particular cases are studied, including the state affine systems and the left invertibility problem. Throughout this paper, examples are discussed in order to highlight our theoretical results.

2. Notations and technical results. Throughout this article, L i-1 f h for 1 ≤ i ≤ n denotes the (i -1) th Lie derivative of output h in the direction of f, and set θ i = dL i-1 f h as its differential. Assume that system (1.1) is locally observable, thus

θ = θ 1 , • • • , θ n
T is a basis of the cotangent bundle T * U of U . Then, we also consider the vector field τ 1 defined in [START_REF] Krener | Linearization by output injection and nonlinear observer[END_REF] as follows

θ i (τ 1 ) = 0, f or 1 ≤ i ≤ n -1, θ n (τ 1 ) = 1, (2.1)
and by induction we define

τ k = (-1) k-1 ad k-1 f (τ 1 ) , f or 2 ≤ k ≤ n. (2.2) It is clear that {τ 1 , • • • , τ n } is a basis of the tangent bundle T U of U .
Let us recall a famous result from [START_REF] Krener | Linearization by output injection and nonlinear observer[END_REF]. Theorem 2.1. The following conditions are equivalent i) There exist a diffeomorphism and an output injection which transform system (1.1) into normal form (1.2) with α k (y) = 1 for 1 ≤ k ≤ n -1.

ii) [τ i , τ j ] = 0 for 1 ≤ i, j ≤ n.

If for some 1 ≤ k ≤ n -1 the functions α k (y) in the form (1.2) are not constant, then ii) of Theorem 2.1 is not fulfilled. Consequently, the rest of this section is devoted to use [τ i , τ n ] in order to determine all the functions α i (y) for 1 ≤ i ≤ n -1.

Lemma 2.2. For a system in the form (1.2) we have for 1 ≤ k ≤ n -1,

τ k = 1 π k ∂ ∂z k + A k k-1 (z n ) z n-1 + η k k-1 (z n ) ∂ ∂z k-1 + k-2 i=1 A k i (z n ) z n-k+i + n-1 j=n-k+i+1 n-1 l=j T k j,l (z n ) z j z l ∂ ∂zi + k-2 i=1 n j=n-k+i+1 η k i (z n ) z j + O [3] zn (z n-k+i+1 , • • • , z n-1 ) ∂ ∂zi , (2.3)
where

π n = 1 and π k-1 = π k α k-1 for 2 ≤ k ≤ n, η k i (z n ) and T k j,l (z n ) are some smooth functions of z n , O [3] zn (z n-k+2 , • • • , z n-1
) represents the residue higher than order 2 with coefficient which is function of z n and

A k i (z n ) = (-1) k-i+1   S k k-i,1 π ′ i π 2 i + k-1 m=k-i+1 S k k-i,m-k+i+1 π ′ k-m π 2 k-m   m j=k-i+1 α k-j     π n-k+i , (2.4)
where S k k-i,1 and S k k-i,m-k+i+1 are defined as follows

S k j,1 = 1, S k j,l = S k-1 j-1,l + S k-1 j,l-1 , f or 2 ≤ k ≤ n, 1 ≤ j ≤ k -1 and 1 ≤ l ≤ k -j. (2.5)
and S k 0,l = S k i,0 = 0.

Proof. For a system in the (1.2) form, equation (2.1) gives

τ 1 = 1 π1 ∂ ∂z1 , then, we use equation (2.2) to obtain τ 2 = 1 π2 ∂ ∂z2 + π ′ 1 π 2 1 π n-1 z n-1 + π ′ 1 π 2 1 β n ∂ ∂z1
, and

τ 3 = 1 π 3 ∂ ∂z 3 + π ′ 1 π 2 1 α 1 + π ′ 2 π 2 2 π n-1 z n-1 + π ′ 1 π 2 1 α 1 + π ′ 2 π 2 2 β n ∂ ∂z 2 - π ′ 1 π 2 1 π n-2 z n-2 + π ′ 1 π 2 1 π n-1 ′ π n-1 z 2 n-1 ∂ ∂z 1 - π ′ 1 π 2 1 π n-1 ′ β n + π n-1 β ′ n z n-1 + π ′ 1 π 2 1 π n-1 β n-1 + β n β ′ n ∂ ∂z 1 .
Then by an induction, for 3 < k ≤ n, we get

τ k = 1 π k ∂ ∂z k + A k k-1 (z n ) z n-1 + η k k-1 (z n ) ∂ ∂z k-1 + k-2 i=1   A k i (z n ) z n-k+i + n-1 j=n-k+i+1 n-1 l=j T k j,l (z n ) z j z l   ∂ ∂z i + k-2 i=1   n j=n-k+i+1 η k i (z n ) z j + O [3] zn (z n-k+i+1 , • • • , z n-1 )   ∂ ∂z i ,
where

A k i (z n ) = (-1) k-i+1   S k k-i,1 π ′ i π 2 i + k-1 m=k-i+1 S k k-i,m-k+i+1 π ′ k-m π 2 k-m   m j=k-i+1 α k-j     π n-k+i ,
with the coefficients S k i given by the rule (2.5).

In order to determine the α i (y) for 1 ≤ i ≤ n -1, we impose that

∂ ∂z i h • φ -1 = 0, for 1 ≤ i ≤ n -1, 1, when i = n.
Now, we are ready to state a set of differential equations which enables us to compute functions

α i for 1 ≤ i ≤ n -1.
Proposition 2.3. If there exists a diffeomorphism which transforms system (1.1) into form (1.2) then

[τ k , τ n ] = λ k (y)τ k + G [1] n + R, f or 1 ≤ i ≤ n -1,
where

G [1] n = k-2 i=1 1 π k T k k,n-k+i z n-k+i ∂ ∂z i + 1 π k T k k,k z k ∂ ∂z 2k-n ,
and

R = k-2 i=1   n j=n-k+i+1 ηk i (z n ) + O [2] zn (z n-k+i+1 , • • • , z n-1 )   ∂ ∂z i and λ k (y) = diag{δ k 1 (y), • • • , δ k i (y), • • • , δ k k (y), 0, • • • , 0}, f or 1 ≤ i ≤ k -1, (2.6)
where

δ k k = A n k + π ′ k π k and δ k i = A n i -A n n-k+i - (A k i ) ′ A k i for 1 ≤ i ≤ k -1, and
A k i is given in (2.4). Proof. According to equation (2.3), for 1 ≤ k ≤ n -1 we have [τ k , τ n ] = A n k + π ′ k π k 1 π k ∂ ∂z k + k-2 i=1 A n i -A n n-k+i - A k i ′ A k i A k i z n-k+i + 1 π k T k k,n-k+i z n-k+i ∂ ∂z i + 1 π k T k k,k z k ∂ ∂z 2k-n + k-2 i=1   n j=n-k+i+1 ηk i (z n ) + O [2] zn (z n-k+i+1 , • • • , z n-1 )   ∂ ∂z i . Set λ k (y) = diag{δ k 1 (y), • • • , δ k i (y), • • • , δ k k (y), 0, • • • , 0}, where δ k k = A n k + π ′ k π k and δ k i = A n i -A n n-k+i - (A k i ) ′ A k i for 1 ≤ i ≤ k -1, then [τ k , τ n ] = λ k (y)τ k + G [1] n + R. (2.7)
Remark 1. In equation (2.7), λ k (y) could be uniquely determined since G [START_REF] Boutat-Baddas | Quadratic Observability normal form[END_REF] n might be separated according to the coefficients of second-order terms in τ n .

Finally, the following result enables us to determine all the functions α i (y) for all 1 ≤ i ≤ n -1. Proposition 2.4. If there exists a diffeomorphism which transforms system (1.1) into form (1.2), then

α i = πi πi+1 for 1 ≤ i ≤ n -2, and α n-1 = π n-1 , where    π i = c i exp exp δ i i -δ n-1 i -δ i+1 i+1 dy -Bn-1 i dy , f or 1 ≤ i ≤ n -2, π n-1 = c n-1 exp δ n-1 n-1 -Ān n-1 2 dy , (2.8)
with Bk 1 = 0 and for

1 ≤ i, k ≤ n -1 and 1 ≤ i ≤ n -1 Bk i = k-1 m=k-i+1 S k k-i,m-k+i+1 π ′ k-m π k-m .
(2.9)

Proof. Define

B k i = π ′ i π i + Bk i .
(2.10)

According to equation (2.4), for 1 ≤ i, k ≤ n -1, A k i ′ A k i = B k i ′ B k i - π ′ i π i + π ′ n-k+i π n-k+i . As δ k k = A n k + π ′ k π k , hence δ n-1 i = A n i -A n 1+i - B n-1 i B n-1 i + π ′ i π i - π ′ 1+i π 1+i = δ i i -δ 1+i 1+i - π ′ i π i + Bn-1 i ′ / π ′ i π i + Bn-1 i . which yields i = c i exp exp δ i i -δ n-1 i -δ 1+i 1+i dy -Bn-1 i dy , f or 1 ≤ i ≤ n -2,
where Bn-1 i is defined in (2.9) and c i ∈ R, c i = 0.

As

δ n-1 n-1 = 2 π ′ n-1 πn-1 + Ān n-1 , where Ān n-1 = n-1 m=2 S n 1,m π ′ n-m πn-m , then π n-1 = c n-1 exp δ n-1 n-1 -Ān n-1 2 dy .
Remark 2. For system (1.2), if we set

α i = s(y) for 1 ≤ i ≤ n -1, then δ n-1 n-1 = A n n-1 + π ′ n-1 π n-1 = 2 π ′ n-1 π n-1 + n-1 i=2 π ′ n-i π n-i .

By the definition of π

i for 1 ≤ i ≤ n -1, we have π k = s n-k for 1 ≤ k ≤ n -1, therefore δ n-1 n-1 = 2 s ′ s + n-1 i=2 i s ′ s = l n s ′ s ,
where l n = n(n-1) 2 + 1. In such a way, we obtain the same result as the one stated in [START_REF] Respondek | Time scaling for observer design with linearization error dynamics[END_REF].

3. Main result. If there exists a diffeomorphism which transforms system (1.1) into form (1.2), then equation (2.8) of Proposition 2.4 gives all α i for 1 ≤ i ≤ n -1. Therefore, let us consider a new family of vector fields defined as follows:

τ 1 = π 1 τ 1 and τ i+1 = 1 α i [ τ i , f ], f or 1 ≤ i ≤ n -1. (3.1) Set θ( τ 1 , • • • , τ n ) =          0 0 • • • 0 1 0 . . . • • • π n-1 l 2,n . . . • • • . . . • • • . . . . . . π 2 • • • • • • . . . π 1 l n,2 • • • • • • l n,n          := Λ,
where

l k,j = θ k ( τ j ) f or 2 ≤ k ≤ n and n -k + 2 ≤ j ≤ n.
Consider the following R n -valued form ω

ω = Λ -1 θ := (ω 1 , ω 2 , • • • , ω n ) T , (3.2)
for 1 ≤ s ≤ n, we have

ω s = n θ m . (3.3)
Then, the following algorithm gives all the components of ω.

Algorithm 1.

f or 1 ≤ j ≤ n, r n,j = • • • = r n-j+2,j = 0 and r n-j+1,j = 1.

f or 2 ≤ k ≤ n -1 and 1 ≤ j ≤ n, r n-k,j = - k i=2 l k,n-k+i-(j-1) r n-k+i-(j-1),j ,

2).

2) There exists a family of functions α i (y) for 1 ≤ i ≤ n -1 such that the family of vector fields τ i for 1 ≤ i ≤ n defined in (3.1) satisfies the following commutativity conditions

[ τ i , τ j ] = 0, f or 1 ≤ i, j ≤ n.
(3.4)

3) There exists a family of functions α i (y) for 1 ≤ i ≤ n -1 such that the R n -valued form ω defined in (3.2) satisfies the following condition dω = 0.

(3.5)

Proof. Assume that there exists a diffeomorphism which transforms system (1.1) into form (1.2), then we compute α i (y) for 1 ≤ i ≤ n -1 from equation (2.8) in Proposition 2.4. Thus, it is easy to show that τ 1 = 1 π1 ∂ ∂z1 which yields that τ 1 = ∂ ∂z1 and then, by construction we obtain τ i = ∂ ∂zi for 2 ≤ i ≤ n. Consequently, we have [ τ i , τ j ] = 0 for 1 ≤ i, j ≤ n.

Reciprocally, assume that there exist α i > 0 for 1 ≤ i ≤ n -1 such that [ τ i , τ j ] = 0 for 1 ≤ i, j ≤ n, then it is well-known ( [START_REF] Isidori | Nonlinear control systems, 2nd edtion[END_REF], [START_REF] Nijmeijer | Nonlinear Dynamical Control Systems Springer[END_REF]) that we can find a local diffeomorphism φ = z such that

φ * ( τ i ) = ∂ ∂z i . As φ * ( τ i ) = ∂ ∂zi is constant, hence ∂ ∂z i φ * (f ) = φ * ([ τ i , f ]) = α i φ * ( τ i+1 ) = α i ∂ ∂z i+1 , thus ∂ ∂zi φ * (f ) = α i ∂ ∂zi+1 for 1 ≤ i ≤ n -1.
Consequently, by integration we obtain: φ * (f ) = A(y)z + β(y). Moreover, as dh

• τ i = 0 for 1 ≤ i ≤ n -1 and dh • τ n = 1, we obtain h • φ -1 = z n .
Finally, in order to prove that in Theorem (3.1) Condition 2) is equivalent to Condition 3), it is sufficient to prove that equation (3.4) is equivalent to equation (3.5).

Recall that for any two vector fields X, Y, we have

dω(X, Y ) = L X (ω(Y )) -L Y (ω(X)) -ω([X, Y ]).
Setting X = τ i and Y = τ j , we obtain

dω( τ i , τ j ) = L A τi ω( τ j ) -L A τj ω( τ i ) -ω([ τ i , τ j ]).
As ω( τ j ) and ω( τ i ) are constant, then we have

dω( τ i , τ j ) = -ω([ τ i , τ j ]).
Because ω is an isomorphism and ( τ i ) 1≤i≤n is a basis of T U, then equation (3.4) is equivalent to equation (3.5). Remarks 1. i) The R n -valued form ω can be viewed as an isomorphism T U n → U × IR n which brings each τ i to the canonical vector basis ∂ ∂z I . Moreover, dω = 0 means that there is a local diffeomorphism φ : U → U such that ω is the tangent map of φ.

ii) The diffeomorphism φ(x) = z is determined by ω = φ * (x), which can be given locally as follows

z i = φ i (x) = γ ω i + φ i (0) f or 1 ≤ i ≤ n,
where γ is a smooth path from 0 to x lying in a neighborhood V 0 ⊆ U of 0.

The following simple example is studied in order to illustrate Theorem 3.1. Example 1. Let us consider the following system

             ẋ1 = γ(y) 1+x4 x 1 x 3 , ẋ2 = β(y) 1+x4 x 1 , ẋ3 = µ(y)x 2 , ẋ4 = γ(y)x 3 , y = x 4 . (3.6) which gives            θ 1 = dx 4 , θ 2 = γdx 3 + γ ′ x 3 dx 4 , θ 3 = γµdx 2 + 2γ ′ γx 3 dx 3 + (γµ) ′ x 2 + (γ ′ γ) ′ x 2 3 dx 4 , θ 4 = γµ β 1+x4 dx 1 + 2γ ′ µ + (γµ) ′ γx 3 dx 2 + 2γ ′ γµx 2 + γ (γµ) ′ x 2 + 3γ (γ ′ γ) ′ x 2 3 dx 3 + O [2] (x 1 , x 2 , x 3 ) θ 1 .
Then we have

τ 1 = 1+x4 γµβ ∂ ∂x1 . Consequently we obtain              τ 2 = 1 γµ ∂ ∂x2 + (1 + x 4 ) γ (γµβ) ′ (γµβ) 2 x 3 ∂ ∂x1 , τ 3 = 1 γ ∂ ∂x3 -γµ (1 + x 4 ) (γµβ) ′ (γµβ) 2 x 3 ∂ ∂x2 + (γµ) ′ (γµ) 2 + β (γµβ) ′ (γµβ) 2 γx 2 ∂ ∂x1 + R 1,3 τ 1 , τ 4 = ∂ ∂x4 + γ ′ γ + (γµ) ′ (γµ) + (γµβ) ′ (γµβ) x 3 ∂ ∂x3 -(γµ) ′ (γµ) + 2 (γµβ) ′ (γµβ) x 2 ∂ ∂x2 + 1 1+x4 + (γµβ) ′ γµβ x 1 ∂ ∂x1 + R 1,4 (z 3 , z 2 ) τ 1 + R 2,3 (z 2 3 )τ 2 .
A straightforward computation gives

δ 1 1 = 2 (γµβ) ′ γµβ , δ 2 2 = -2 (γµβ) ′ γµβ , δ 3 3 = 2 γ ′ γ + (γµ) ′ γµ + (γµβ) ′ γµβ , δ 3 1 = 4 (γµβ) ′ (γµβ) - (γµβ) ′ (γµβ) ′ / (γµβ) ′ (γµβ) , δ 3 2 = -2 γ ′ γ + (γµ) ′ γµ + 3 (γµβ) ′ γµβ - (γµ)
′ γµ + (γµβ) ′ γµβ ′ / (γµ) ′ γµ + (γµβ) ′ γµβ .
According to equation (2.8) in Proposition 2.4, we obtain

       π 1 = c 1 exp exp δ 1 1 -δ 3 1 -δ 2 2 dy dy = c 1 γµβ, π 2 = c 2 exp exp δ 2 2 -δ 3 2 -δ 3 3 dy - π ′ 1 π1 dy = c 2 γµ, π 3 = c 3 exp 1 2 δ 3 3 - π ′ 1 π1 - π ′ 2 π2 dy = c 3 γ. Thus α 1 = π1 π2 = c1 c2 β, α 2 = π2 π3 = c2
c3 µ and α 3 = π3 π4 = c 3 γ, so the new vector fields are

τ 1 = c 1 (1 + x 4 ) ∂ ∂x 1 , τ 2 = c 2 ∂ ∂x 2 , τ 3 = c 3 ∂ ∂x 3 , τ 4 = ∂ ∂x 4 + x 1 1 + x 4 ∂ ∂x 1 .
It is clear that [ τ i , τ j ] = 0 for all 1 ≤ i, j ≤ 4. Therefore, according to Theorem 3.1, system (3.6) can be transformed into SODO normal form (1.2).

Moreover as

Λ =     0 0 0 1 0 0 γ γ ′ x 3 0 γµ 2γ ′ γx 3 (γµ) ′ x 2 + 2 (γ ′ γ) ′ x 2 3 γµβ 2γ ′ µ + (γµ) ′ γx 3 2γ ′ γµx 2 + γ (γµ) ′ x 2 + 6γ (γ ′ γ) ′ x 2 3 γ x1 (1+x4) 2 µβ + R     , where R = O [2]
x4 (x 1 , x 2 , x 3 ), a straightforward computation gives

ω = Λ -1 θ = d x1 c1(1+x4) , d x2 c2 , d x3 c3 , dx 4 T .
As ω = dφ, thus the diffeomorphism which transforms system (3.6) into SODO normal form (1.2) is

φ(x) = z = ( x 1 c 1 (1 + x 4 ) , x 2 c 2 , x 3 c 3 , x 4 ) T .
with which system (3.6) could be transformed into

       ż1 = 0, ż2 = c1 c2 β(y)z 1 , ż3 = c2 c3 µ(y)z 2 , ż4 = c 3 γ(y)z 3 .
So far, in this paper, we have only considered systems without inputs. The next section is devoted to systems that are also driven by an input term.

4. Extension to systems with inputs. Consider a system with inputs in the following form ẋ = f (x) + g(x, u), y = h(x), (

where x ∈ U ⊂ IR n , f : U ⊂ IR n → IR n , g : U × IR m → IR n , h : U ⊂ IR n → IR are analytic functions and for x ∈ U, g(x, 0) = 0.

For system (4.1), the SODO normal form along its output trajectory y(t) is as follows

ż = A(y)z + β(y) + η(y, u), y = z n = Cz, (4.2)
where A(y) and β(y) are given in (1.2) and η(y, u) = η 1 (y, u), η 2 (y, u),

ii) [g, τ i ] = 0 for 1 ≤ i ≤ n -1.
Proof. From Theorem 3.1, we can state that there exists a diffeomorphism φ such that φ * (f ) = A(y)z + β(y).

For 1 ≤ i ≤ n -1, because φ * ( τ i ) = ∂ ∂zi is constant, hence we have ∂ ∂z i φ * (g) = φ * ([g, τ i ]) = 0.
Therefore φ * (g) = η(y, u). Thus, we obtain the form (4.2). Let us now study some special cases of the term η(y, u). Corollary 4.2. Assume that conditions i) and ii) of Theorem 4.1 are fulfilled,

Remark 3. If g(x, u) = g 1 (x)u 1 + • • • + g m (x)u m ,
a) if [g, τ n ] = 0, then η(y, u) = η(u). b) if g(x, u) = g 1 (x)u 1 + • • • + g m (x)u m and [g k , τ i ] = 0, f or 1 ≤ i ≤ n and 1 ≤ k ≤ m, then η(y, u) = B 1 u 1 + • • • + B m u m ,
where B i are constant vector fields.

Example 2. Let us consider the following system A straightforward computation gives:

         ẋ1 = γ(y) 1+x3 x 1 x 2 + x1 1+x3 u, ẋ2 = µ(y) 1+x3 x 1 , ẋ3 = γ(y)x 2 + u, y = x 3 .
τ 1 = 1 + x 3 γµ ∂ ∂x 1 , τ 2 = 1 γ ∂ ∂x 2 + (1 + x 3 ) (γµ) ′ γµ 2 x 2 ∂ ∂x 1 , τ 3 = ∂ ∂x 3 + (µγ) ′ (µγ) + γ ′ γ x 2 ∂ ∂x 2 + 1 1 + x 3 - (γµ) ′ γµ x 1 ∂ ∂x 1 .
Then, we obtain

δ 1 1 = 0, δ 2 2 = 2 γ ′ γ + (µγ) ′ (µγ) , δ 2 1 = - (µγ) ′ (µγ) -2 γ ′ γ - (µγ)
′ (µγ) ′ / (µγ)
′ (µγ) .

Then, according to (2.8), we have

π 1 = c 1 exp exp δ 1 1 -δ 2 1 -δ 2 2 dy dy = c 1 γµ, π 2 = c 2 exp 1 2 δ 2 2 - π ′ 1 π1 dy = c 2 γ.
which yields α 1 (y) = c1 c2 µ (y) and α 2 (y) = c 2 γ (y) . Therefore, we obtain

τ 1 = c 1 (1 + x 3 ) ∂ ∂x1 , τ 2 = c 2 ∂ ∂x2 and τ 3 = ∂ ∂x3 + x1 1+x3 ∂ ∂x1 . As g = x1 1+x3 ∂ ∂x1 + ∂ ∂x3 = τ 3 then [g, τ 1 ] = [g, τ 2 ] = 0 and system (4.3) is transformed into        ż1 = 0, ż2 = c1 c2 µ (y) z 1 , ż3 = c 2 γ (y) z 2 + u, y = z 3 . (4.4) by the following diffeomorphism φ(x) = z = ( x 1 c 1 (1 + x 3 ) , x 2 c 2 , x 3 ) T .
Following the proposed high gain observer in the form (1.3), the corresponding observer for the system (4.4) can be designed as follows: In addition, in order to solve the left invertibility problem, the Observability Matching Condition (OMC) for system (4.1) with m = 1 is as follows

     . ẑ1 = -ρ 3 γµ (ẑ 3 -z 3 ) , . ẑ2 = c1 c2 µ (y) ẑ1 -3 ρ 2 γ (ẑ 3 -z 3 ) , . ẑ3 = c 2 γ (y) ẑ2 -3ρ (ẑ 3 -z 3 ) + u, 0 0.5 1 1.5 2
L g L i-1 f h = 0, ∀x ∈ U, 1 ≤ i ≤ n -1, L g L n-1 f h = 0.

Conclusion.

In this paper, we have put forward the geometrical conditions which allow us to a nonlinear system can be transformed locally into the SODO normal form by means of a diffeomorphism and of an output injection. In our main result we state two equivalent ways to check these conditions. In the first one, we have used Lie brackets commutativity and the second one was based on the Moreover, an extension of our results is stated for a class of nonlinear systems with inputs.

  and then, equation (3.3) becomes: ω s = n-s+1 m=1 r s,m θ m . Theorem 3.1. The following conditions are equivalent 1) There exists a diffeomorphism which transforms system (1.1) into a SODO normal form (1.

  and also both conditions i) and ii) of Theorem 4.1 are fulfilled, thenη(y, u) = B 1 (y)u 1 + • • • + B m (y)u m .

1 Fig. 4 . 1 .

 141 Fig. 4.1. Observation error between z 1 and ẑ1

2 Fig. 4 . 2 . 3 Fig. 4 . 3 .

 242343 Fig. 4.2. Observation error between z 2 and ẑ2

Corollary 4 . 3 .Remark 4 . 3 . 3 .

 43433 Assume conditions i) and ii) of Theorem 4.1 are fulfilled and the OMC is verified thenη(y, u) = η 1 (y, u), 0, • • • , 0T . The OMC for system (4.1) with m = 1 is equivalent to g ∈ span{ τ 1 }. We give another example in order to highlight Corollary 4.Example Consider the following system       ẋ1 = u, ẋ2 = µ(y)x 1 + µ(y)x 2 1 + x2 1+x1 u, ẋ3 = γ(y) x21+x1 , y = x 3 .

. 1 ∂∂x1 + c 1

 11 From equation(2.8)in Proposition 2.4, we can determine α 1 (y) = c1 c2 µ (y) and α 2 (y) = c 2 γ (y). Thus, we haveτ 1 = c x2 1+x1 ∂ ∂x2 , τ 2 = c 2 (1 + x 1 ) ∂∂x2 and τ 3 = ∂ ∂x3 . As g ∈ span{ τ 1 }, then the OMC condition is fulfilled, therefore system (4.5) could be transformed by the following diffeomorphismφ(x) = z = c2 µ (y) z 1 , ż3 = c 2 γ (y) z 2 , y = z 3 .

  • • • , η n (y, u) T .

	Theorem 4.1. System (4.1) can be transformed into SODO normal form (4.2) by a diffeomorphism if
	and only if
	i) one of conditions in Theorem 3.1 is fulfilled.

Acknowledgments. The authors are deeply grateful to the anonymous reviewers for valuable comments and helpful suggestions that enable us to improve the presentation of this paper. We are also deeply grateful to Professor Vincent Maki for many corrections.