S Livet

M Collet

M Berthillier

P Jean

J M Cote

Turbomachinery blades damping thanks to optimized shunted piezoelectric circuits

Keywords: Piezoelectricity, Semi-passive control, Damping, Turbomachine, Single-crystal, Negative capacitance

Dynamics of gas-turbine blades are particularly aero-elastic coupling sensitive. These aerodynamic limits can be pushed away by adding extra damping to the structure in order to reach even better compressor performance. However nowadays design and manufacturing techniques in aero-mechanics are achieving their maximum of state-of-the-art.

As in many fields active control would solve easily this kind of instability. But the difficulty remains in the needed energy supply for actuators whereas these components are aimed to be bonded on rotating structures. The capacity of different auto-supplied devices using shunted piezoelectric circuits had been studied here to prevent turbomachine bladed from fluttering.

Before realizing the study on complex turbomachine geometries, the presented technique uses a numerical development thanks to a 1D Euler-Bernoulli beam model combining both mechanical and electrical coupling parameters. A second development thanks to a 3D model had been made using a commercial tool, Comsol software. These approximate models are used to optimize electrically the shunted piezoelectric element and its localization. The results, verified experimentally, let suppose that vibrations can be reduced significantly when shunted piezoelectric circuits are mounted on a real structure.

INTRODUCTION

Each new gas turbine generation require more and more demanding design techniques to correspond to the extreme environments and conditions of use. These have to be especially more and more light and efficient. The design methods for these machines imply an augmentation in volume and speed of the airflow to achieve such proportion of thrust/mass with the use of less components or lighter materials. The space between blade stages in compressors are also tighter with a lower number of stages. Consequently new compressors use extremely long profile blades, which rises the aerodynamics stresses.

The vibrations induced by the airflow in turbomachines are the cause of High Cycle Fatigue (HCF) and ask for regular and unexpected maintenance needs. The minimization of the vibration responses toward to the aerodynamics forces can reduce the HCF and by the way reduce significantly the deterioration phenomena and avoid high cost in repairing. The need of new solutions is trivial as most of the traditional solutions of minimization of HCF effects bring negative effect on the new effective designs, adding mass to the structure for instance. Moreover as the elderly concept engines may match the new expectations, after-design solutions have to be suggested and optimized easily without replacing any components.

The piezoelectric material present a new approach to vibration suppression in aeronautical structures [START_REF] Cross | Shunted piezoelectrics for passive control of turbomachine blading flow-induced vibrations[END_REF][START_REF] Agneni | Shunted piezoelectric patches in elastic and aeroelastic vibrations[END_REF] . They can be adapted to the creation of damping or control elements. Indeed when the electric current is applied to a piezoelectric material, a corresponding strain is produced within the material, which produces an electric voltage. Vibration damping of an elastic structure can be obtain by an piezoelectric actuator embedded into it and shunted to an electric resonant circuit adjusted to the mode to damp. The obtained damping depends on the type of electronic shunt and can be expressed directly taking into account the inherent structural damping coefficient or unstable phenomena as fluttering of the studied system.

In this paper, we present a simple Euler-Bernoulli beam model created on the Hamilton Principle [START_REF] Piefort | Finite Element Modeling of Piezoelectric Active Structures[END_REF][START_REF] T I E R S T E N | Linear Piezoelectric Plate Vibration[END_REF] to be able to describe the studied phenomena and explain the role of the piezoelectric component in energy transfer between the schematically modeled blade structure and the electronic circuit. A more complex 3D model thanks to the Comsol software has been also achieved. First optimization guidelines have been highlighted to maximize the induced damping. An experiment has been settled to verify theoretical results on damping with a passive shunt the damping realized thanks to different piezoelectric shunt sticked onto steal beams of different length in order to be able to quantify the capability of passive control and check how much the piezoelectric system can likely dampen passively the structure.

EULER-BERNOULLI BEAM MODELING THANKS TO HAMILTON PRINCIPLE

Energy description

The Hamilton principle 5 is often used to construct Finite Element models. The beam model studied and exploited in this paper has been built on the same criteria. The variational principal can be written as δ t1 to (L + W)dt =0 with L Lagrangian term and W the virtual work of mechanical and electrical forces. L = T + H with T kinetic energy and H the enthalpy which worths by the IEEE definition 6 of the piezoelectricity, H = U -E i D i with U the intern energy of the studied model.

Kinetic Energy expression

The kinetic energy expression can be expressed as T = [START_REF] Cross | Shunted piezoelectrics for passive control of turbomachine blading flow-induced vibrations[END_REF] 2 Ω ρ ui ui dΩ, which leads to:δT = Ω ρδ(u i )u i dΩ

Strain Energy Expression

For a piezoelectric element we can express the strain energy as:

H = 1 2 Ω † {S}{T }dΩ -1 2 Ω † {E}{D}dΩ
The piezoelectric equations of behaviour are:

{T } = c E {S}- † [e] {E} (1)
{D} =[ e] {S} + ε S {E} (2)
where {T } represents the stresses tensor, {S} the strain tensor, {E} the electric field, {D} the electric displacement, [c E ] the elastic coefficients at constant electric field {E}, {ε S } the dielectric coefficients at constant strain {S},a n d[ e] the piezoelectric coupling coefficients.

We replace {T } et {D} by their respective expression in the Hamilton principle. Considering {S} as a symmetric field, we find for H:

t1 to δHdt = Ω t1 to † {δS} c E {S} dtdΩ - Ω t1 to † {δS} † [e] {E} dtdΩ - Ω t1 to † {δE} [e] {S} dtdΩ - Ω t1 to † {δE} ε S {E} dtdΩ( 3 )

Euler-Bernoulli Beam Model

Kinematic hypothesis of Euler-Bernoulli beam model

The initial 3D problem of elasticity is reduced to 1D problem. All functions (displacement, stresses...) depend on only one variable along the length of the beam, here noted x. We suppose also that the thickness and the width are small in comparison with the length and the beam is charged in its plane. The transverse strain is neglected in this model.

Strain expression in the case of a Euler-Bernoulli beam

Some articles show studies 7 with piezoelectric patches facing each other. These sandwich systems usually show a better efficiency of their system while the piezoelectric elements are perfectly bonded in opposite. In the case of aeronautical structure as blades, this position could not happen and that is why all our study is done with piezoelectric elements placed independently. For a beam with a piezoelectric element bond onto the surface the representative drawing of the model is the following:

✻ w θ ■ ✲ ✛ L s ✲ ✛ L p ✻ ❄ t s ✻ ❄ t p ✲ x ✻ z
We consider the beam widths l s et l p respectively the one of the structure and of the piezoelectric. The variables w et θ correspond to the degrees of freedom of the beam.

By definition the termes ij of the strain tensor i in the hypothesis of small linear perturbations can be expressed as S ij = 1 2 dui dxj + duj dxi . Further we will note S ij = S i for i = j and 2S ij = S 9-(i+j) for i = j.

Piezoelectric coefficient correction

The characteristic coefficients of a piezoelectric are affected by the Euler-Bernoulli type modeling. The beam hypothesis implies a reduction of the number of piezoelectric parameters. We can consider the transverse efforts T xz or T 5 neglected. The efforts T 2 , T 3 , T 4 et T 6 are considered equal to zero by the Euler-Bernoulli model.

For a piezoelectric element the equations can be written as:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ T 1 T 2 T 3 T 4 T 5 T 6 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣
c 11 c 12 c 13 000 c 12 c 22 c 13 000 c 13 c 13 c 33 000 000 c 44 00 0000 c 44 0 00000 c 66

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ E ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ S 1 S 2 S 3 S 4 S 5 S 6 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ - ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣
00 e 31 00 e 31 00 e 33 0 e 15 0 e 15 00 000

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎨ ⎩ E 1 E 2 E 3 ⎫ ⎬ ⎭ (4)
Which gives us:

T 1 = c E 11 -c E 33 c E 11 -c E 13 2 -1 c E 12 c E 33 c E 12 -c E 13 2 + c E 13 2 c E 12 -c E 11 S 1 + c E 33 c E 11 -c E 13 2 -1 c E 12 -c E 33 e 31 + c E 13 e 33 + c E 13 -c E 13 e 31 + c E 11 e 33 -e 31 E 3 (5)
By the same way we calculate the electrical equation part:

⎧ ⎨ ⎩ D 1 D 2 D 3 ⎫ ⎬ ⎭ = ⎡
⎣ 0000 e 15 0 000 e 15 00 e 31 e 31 e 33 00 0

⎤ ⎦ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ S 1 S 2 S 3 S 4 S 5 S 6 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ + ⎡ ⎣ ε 11 00 0 ε 11 0 00 ε 33 ⎤ ⎦ S ⎧ ⎨ ⎩ E 1 E 2 E 3 ⎫ ⎬ ⎭ (6)
D 3 = e 31 S 1 + e 31 S 2 + e 33 S 3 + ε 33 E 3 (7) = e 31 S 1 -e 31 c E 33 c E 11 -c E 13 2 -1 c E 33 c E 12 -c E 13 2 S 1 + c E 33 e 31 + c E 13 e 33 E 3 -e 33 c E 13 2 -c E 11 c E 33 -1 c E 13 c E 12 -c E 11 S 1 + -c E 13 e 31 + c E 11 e 33 E 3 + ε 33 E 3 (8)
when we write the system as in the manner developed by Maurini 8 , the new piezoelectric parameters are:

{T 1 } = c E 11 * {S 1 }-e * 31 {E 3 } (9) {D 3 } = e * 31 {S 1 } + ε S 33 * {E 3 } (10) c E 11 * = c E 11 + c E 33 c E 11 -c E 13 2 -1 2c E 12 c E 13 2 -c E 12 2 c E 33 + c E 13 2 c E 11 ( 11
)
e * 31 = e 31 + c E 33 c E 11 -c E 13 2 -1 -c E 12 c E 33 -c E 13 2 e 31 + c E 13 c E 12 -c E 11 e 33 (12)
ε S 33 * = ε S 33 + c E 33 c E 11 -c E 13
In the case of a Euler-Bernoulli type resolution the displacement u i which is defined as u 1 = -z ∂w ∂x , u 2 =0 et u 3 = w(x, t)w i t hw(x, t) displacement of a point distant from the neutral fiber and ∂w ∂x (x, t)=-θ(x, t) corresponding to the rotation of the section at the coordinate x. We obtain the strain

S 11 = S 1 = ∂H1 ∂x = z ∂ 2 w ∂ 2 x
. We replace the obtained expressions into the Hamilton principle:

t1 to δHdt = Ω t1 to c E 11 * z 2 ∂ 2 δw dx 2 ∂ 2 w ∂x 2 dtdΩ - Ω t1 to δE 3 e * 31 z ∂ 2 w ∂x 2 dtdΩ - Ω t1 to ∂ 2 δw ∂x 2 e * 31 zE 3 dtdΩ - Ω t1 to δE 3 ǫ S 33 * E 3 dtdΩ (14)

Electrical model hypotheses

The electrical field can be defined as a potential gradient φ such as:

→ E= - → ∇ φ.
T h e nφ corresponds to the difference of electrical voltage between the two electrodes of the piezoelectric element. We can approximate the electrical field as constant through the thickness, which implies: E 3 = -φ dz (z direction through the thickness of the piezoelectric element). We suppose that the inferior electrode is grounded: φ(x, y, z, t)=φ(z, t)a n dw e define z p the height of the inferior electrode z p = ts 2 . At the superior electrode we find z = z p + t p = ts 2 + t p . and we have φ(z p ,t)=V (t) -V 0 = V (t), which gives us as solution φ = V tp (z -ts 2 )

E 3 = - ∂φ ∂z = - V t p z (15)
In the framework of these hypotheses we find as piezoelectric strain expressions:

t1 to δHdt = Ω t1 to c E 11 * z 2 d 2 δw ∂x 2 ∂ 2 w

FINITE ELEMENT DECOMPOSITION

1D Finite Element model

We replace, in the strain energy expression, the precedent displacement as a polynomial function and split up into two distinct functions independent respectively from time and space parameters w(x, t)=W (x)η(t). We can distinguish three terms in the strain energy expression 9 : a pure mechanical parameter δH m , two electromechanical coupling parameters δH em et δH me and an electrical parameter δH e .

Let i be the index summing the structural elements s and the piezoelectric ones p.

δH m = † {δη(t)} Li x=0 E i I i † [W ′′ (x)][W ′′ (x)] dx {η(t)} ( 17
)
δH me = δH em = -δV Lp x=0 lp(ts+tp) 2 e * 31 [W ′′ (x)] dx {η(t)} and δH e = - lpLp tp ǫ S 33 * δV V
The elementary stiffness and mass matrices for the system {beam+piezoelectric} can be calculated from the strain energy expression.

Stiffness Matrix

The elementary stiffness matrix can be written as equation 18 with I structural inertia developed in Appendix B.

K e = K e mm K e me K e em K e ee = c E 11 * IL e † [W ′′ ][W ′′ ] lp(ts+tp ) 2 e * 31 † [0, -1, 0, 1] lp(ts+tp) 2 e * 31 [0, -1, 0, 1] - lpLp tp ε S 33 * (18)
We note the piezoelectric characteristics: α =

Mass Matrix

The mass matrix M of the system comes from the kinetic expression T previously calculated and using the same methodology as before for the stiffness matrix. δT = -ρ i Ωi δu i üi dΩb u tu 1 ≈ 0,

u 2 = z ∂w ∂x et u 3 = w.
The influence of the inertia is neglected the kinetic energy.

δT = -ρ i S i † {δη(t)} Li x=0 † [W (x)][W (x)] dx {η(t)} (19)
From the expression 19,w ed r a go u tt h ee l e m e n t a r ym a t r i x[ M e ]:

[M e ]=ρ i I i Li x=0 † [W ′ (x)][W ′ (x)] dx (20)
To complete the study the pure mechanical matrices of a Euler-Bernoulli model are recalled in Appendix C.

FE electromechanical system

For a FE system the assembly of the whole system linked to a shunt circuit in dynamic and in the case of a electrical shunt realized between the electrodes, we have for a piezoelectrical mechanical system the following relation written for one element:

[M mm ]{ ẅ} +[G mm ]{ ẇ} +[K mm ]{w} +[K me ]{V } = {F } (21) [K em ]{w}-[K ee ]{V } = -{q} (22)
I nt h ec a s eo fas h u n tw i t haZ impedance, V = -sZq with s Laplace coefficient.

[

M mm ]{ ẅ} +[G mm ]{ ẇ} + [K mm ]+[K me ][K ee ] -1 [K em ] {w} +[K me ][K ee ] -1 {q} = {F } (23) sZ[K ee ]{q} + {q} +[K em ]{w} = 0 ( 24
)
We can notice that in the case of short circuit V = 0 and we return to a purely mechanical system:

[M mm ]{ ẅ} +[G mm ]{ ẇ} +[K mm ]{w} = {F } (25)
and in the case of open circuit q =0:

[M mm ]{ ẅ} +[G mm ]{ ẇ} + [K mm ]+[K me ][K ee ] -1 [K em ] {w} = {F } (26)
The stiffness value increases from a parameter

K ′ mm =[K me ][K ee ] -1 [K em ]
linked to the electromechanical (the coupling α) and electric (the inherent capacity C p ) parameters. This frequency gap between open and short circuit is regularly used as a optimization criteria 10 .

3D modeling and state space system resolution

The 3D model represented in figure 1 is realized with the Comsol myltiphysics software. The FE electromechanical system from equation 3 is firstly reduced by the methodology developed by Collet [START_REF] Collet | Modal synthesis and dynamical condensation methods for accurate piezoelectric systems impedance computation[END_REF] . This methodology condensation allows a more accurate computation of the piezoelectric impedances. The induced matrices are exploited to create the final model in Matlab. To connect the FE model to shunt circuit we write the system as a Linear Time Invariant System described just below:

Ṅ = AN + BE S = CN + DE (27)
We define the observability variable with all the displacements x = φη with φ functions based on the condensation methodology [START_REF] Collet | Modal synthesis and dynamical condensation methods for accurate piezoelectric systems impedance computation[END_REF] :

N = η η ( 28
)
The input of the system:

E = V F (29)
The output with w = C d1 x deflection of the beam tip:

S = q w (30)
We write the initial system as a linear state problem with M, K and G respectively Mass, Stiffness, Damping Matrices:

ẋ ẍ = 0 I -M -1 K -M -1 G x ẋ + 00 M -1 F u1 M -1 F u2 E (31) S = C d1 0 C d2 0 x ẋ +[D d 0] E (32)
In the case of a shunted piezoelectric with a electrical circuit of Z impedance. We have also the relation V = -Zsq with s the Laplace operator. Here we connected to the piezoelectric a resistive R s and negative capacity C s in series and Z = R s s +1/C s and q = C d1 x + D d V as in our experimental work (developed in section 4.3) but the method can be used for many other kinds of electrical circuits.

V = -(R s s +1/C s )(C d1 η + D d V )=-R s C d1 η -R s D d V -C -1 s C d1 η -1/C s D d V (33)
. :

Hence V = -C d1 D -1 d η -(C -1 s C d1 )(R s D d ) -1 η -(C -1 s D d +1)(R s D d ) -1 V ; Ṅ = AN + BE V = -C d1 D -1 d -(C -1 s C d1 )(R s D d ) -1 N +[ -(C -1 s D d +1)(R s D d -1 0] E (34) Ṅ = AN+ BE V = C ′ N +[ Dd ′ 0] E (35)
We can write the new linear state system taking into account the electronic circuit with a new observation parameter N ′ = N V ,e n t r yF and output

S ′ = V w . Ṅ ′ = AM -1 F u1 C ′ D ′ d N ′ + M -1 F u2 F (36)
S ′ = ⎡ ⎣ C d2 00 00 0 00 1 ⎤ ⎦ N ′ +0 F (37)
3.4 First study of optimal piezoelectric damping with a resistive shunt Many studies [START_REF] Benjeddou | Advances in piezoelectric finite element modeling of adaptive structural elements: a survey[END_REF][START_REF] Lesieutre | Vibration damping and control using shunted piezoelectric materials[END_REF] have detailed some damping optimization with several electronic shunts involving inductive, resistive, capacitive or switching components. In this paper we focused on the capability of passive damping on structure as it is very difficult, if not impossible, to supply the electric devices in aeronautical structures. Therefore we glanced preferentially at the mechanical coupling carried out between the structure and the piezoelectric element and which plays a very high role in passive control. A first study has been realized to observe the evolution of the maximal damping ratio obtained by a purely resistive shunt as a function of the length of the piezoelectric element from the resolution of sections 3. In this particular case we observe that the maximum of damping, calculated at an optimal resistance for each length ratio, is obtained for a piezoelectric length around the middle of the beam (fig. 2). At the same time of the decreasing of the damping ratio, the added mass on the second part of the beam, where the beam strain is the weaker, makes also the eigen frequency of the first mode diminish. The beam is excited by an electro-magnetic generator with a residual magnetic field to be able to excite the metal beam at distance. The damping on the first mode from a white noise is measured by a laser velocimetre and acquired by a PC with Siglab software.

Resistive shunt 4.2.1 Damping influence from the length ratio of piezoelectric patches

We desired to quantify the capability of damping impact thanks to a resistive shunt and by this way evaluate the maximum of coupling between the piezoelectric and the beam. When we make the resistance shunt on the piezoelectric patch vary from 0 (short circuit) to +∞ (open circuit). This phenomena can be explained by the direct impact which the resistance has onto the stiffness of the whole structure {piezo+beam} as explained in section 3.2. Increasing the resistance implies so an increase of the damping ratio at the beginning and a decrease after passing by an optimum, while the frequency keeps increasing.

As u m m a r yt a b l e1 of all the experiments is drew up for the first mode.

Damping impact due to piezoelectric material type

Two different types of piezoelectric material have been also tested during the experimental studies. Keeping the same size of material, we wanted to test the capability of the new material PMN-33PT from TRS Technology in comparison with PZTs (Appendix A) [START_REF] Cao | Elastic, piezoelectric and dielectric properties[END_REF] .T h ec a p a b i l i t i e so fs u c ham a t e r i a lw e r ea l r e a d yp r o v e di nv i b r a t i o n absorption [START_REF] Rusovici | Design of a single-crystal piezoceramic vibration[END_REF] The PMN-33PT is a single crystal material and is a relatively softer piezoelectric element, which induce lower eigen frequencies of the whole structure. It has so less impact on the host behaviour and deflection of the beam. Moreover it has also better coupling coefficients. For a 55mm beam we observe experimentally a damping gain of about 100% with a PMN-33PT in comparison to the used PZT in section 4.2.1 .T h ee x t r a damping ξ on the structure reaches about 2%, which corresponds to an amplitude loss of -25dB (fig. 4). The relative error which can be read in table 2 on the expected theoretically damping is rather small. We pointed out that most of the frequencies are overestimated which can be explained by the fact that the glue is not taken into account in the different models. In consequence we observe that the glue may add some mass on thestructure butseemstohaveapoor impacton thedamping ratio.

Negative capacitance shunt

To increase the efficiency of the damping due to Joule Effect and heat consummation through the resistance, we studied also the possibility to prevent the energy from getting stocked into the inherent capacity. The study The device works as an equivalent negative capacity C eq = -R7+γR8

(1-γ)R8 C 1 for frequencies above 1 R3C1 .I tc a n be tuned thanks to γ through the potentiometer R 8 variable from 0 to 20kΩ, which explains the lines break on the evans diagram (fig. 6). We can then observe the very high degree of damping thanks to this system in figure 6, where resistive damping evans loops are represented for different values of negative capacities. With a simple PZT bond on a 55mm long beam, we gain more than 4% of damping and -30dB in magnitude on the first mode just before the instability, in comparison to 1% of damping and -18dB in magnitude with a simple resitive shunt (section 4.2.1). Indeed the system remains stable for negative capacities below the capacity at contant stress C T p . This instability corresponds to the highest pick visible on the bode diagram on figure 6 traced for a capacity C s just above -C T p .

CONCLUSION

With the arrival of new designs of engine components and the severe requirements which they have to follow, innovation techniques to control blade vibrations in gas turbines have to be considered and studied. The critical vibration modes are associated to the aspect ratio of blades. A high ratio makes the blade considered as a beam with a typical response of low orders as bending or torsion modes. The model presented here corresponds to these systems. However the characteristics of blade vibrations with a low aspect ratio are closer to a shell model and present more complex modes. Induced fluid-structure coupling vibrations present in this last case have a very low inherent structural damping so this technique of damping with a piezoelectric shunt can be very interesting too if applied here.

The potential of such a damping on a turbomachine is very interesting in theory as shown by the theoretical models presented here and validated by the experimental work. Indeed the gain obtained experimentally is as good as expected from a two-degree-of-freedom model (Euler-Bernoulli Beam) suggested by the variational p r i n c i p l ea sw e l la sf r o mt h e3 DC o m s o lm od e l . T h ep u r e l yp a s s i v em e t h ods e e m ss ot oc o n c l u d ee x pe r i m e n t a l l y to a very large damping addition (almost 2% with a PMN-33PT and about -25dB only with a resistance), which is actually globaly efficient enough to prevent instability from aerodynamics coupling to occur. The negative capacity circuit presents better results, however it can present some difficulties to be adjusted correctly. The passive control remains also always stable and can be adapted more easily to a muti-modal damping with several piezoelectric element correctly adjusted and wisely localised on the structure. For a complex blade geometry, the theoretical model is certainly to crude to represent fairly enough an hypothetic experimental model and will be soon enriched, the experimental perspectives being to be able to implement a passive control on fan blades or disk-assemblies [START_REF] Agnes | Piezoelectric coupling of bladed-disk assemblies[END_REF] . The aim of our next studies will be to prevent fluttering phenomena from appearing and to push as far as possible the aerodynamics limits within the material resistance frontiers thanks to the optimisation presmiss guidelines.

APPENDIX A. PIEZOELECTRIC PROPERTIES

APPENDIX C. EULER-BERNOULLI COEFFICIENTS

The Hamilton principle leads naturally to the pure mechanical results for a Euler-Bernoulli beam that we find out developing the term † [W ′′ ][W ′′ ]w i t hW (x) expressed as with {w 1 ,θ 1 } and {w 2 ,θ 2 } respectively degrees of freedom of nodes at the coordinate x 1 and x 2 :

W (x)= 1 - 3 L 2

  inherent capacity at constant strain of the piezoelectric element of length L e .

Figure 1 .

 1 Figure 1. Picture from 3D Comsol model

  2 and 3.3 for an electric impedance Z = R. The piezoelectric is kept at the base of the root of the cantilever beam and we took here a 55mm long-10mm wide-1mm thick beam. The piezoelectric has the same width as the beam and is 0.5mm thick. These values correspond to the values used later in the experimental work (section 4.1). The simulations are realized for a piezoelectric element from Saint Gobain Quartz (P1 94).

Figure 2 .

 2 Figure 2. Comparison of the evolution of the length of the piezoelectric element modeled in 1D and 3D

1 Figure 3 .

 13 Figure 3. Picture of the experimental assembly

Figure 4 .

 4 Figure 4. Bode and Evans diagrams of a 55mm beam with a PMN-33PT at different resistive shunt

Figure 5 .

 5 Figure 5. Electronic negative capacity diagram

Figure 6 .

 6 Figure 6. Bode and Evans diagrams of a 55mm beam with a PZT at different negative capacity and resistive shunt

  -e 31 c E 33 e 31 + c E 13 e 33 + e 33 -c E 13 e 31 + c E 11 e 33

	2 -1

Table 1 .

 1 Passive damping results from experiments

		55mm beam 70mm beam 90mm beam
	Short-circuit frequency	360.3H z	210.0H z	121.4H z
	Open-circuit frequency	366.5H z	212.7H z	122.7H z
	Structural damping (%)	-0.1768	-0.1295	-0.1376
	Maximal damping at Ropt (%)	-0.9764	-0.744	-0.654
	δ damping at Ropt (%)	-0.80	-0.61	-0.52
	δ magnitude damping	-18 dB	-13 dB	-12 dB

Table 2 .

 2 Passive damping results from 1D and 3D theory

	Beam length	55mm	70mm	90mm
	Model type	1D	3D	1D	3D	1D	3D
	Short-circuit frequency	361.0H z 354.9H z 213.9H z 211.0H z 124.0H z 122.7H z
	Experimental relative error (%)	0.194	-1.50	1.86	0.48	2.1	1.07
	Open-circuit frequency	367.0H z	361 Hz	216.9H z 213.9H z 125.4H z 124.1H z
	Structural damping (%)	-0.18	-0.18	-0.13	-0.13	-0.14	-0.14
	Maximal damping at Ropt (%)	-0.995	-1.00	-0.832	-0.811	-0.655	-0.666
	Experimental relative error (%)	1.90	2.4	11.8	9.0	0.153	1.83

Table 3 .

 3 E (10 9 N.m -2 ) c E Modulus parameters

			11	c E 12	c E 13	c E 33	c E 44	c E 66
	P1 94		153.1 104.9 103.9 136.9 37.5 24.1
	P1 91		120.9 76.3	73.1 112.6 33.6 22.3
	PMN-33PT	115	103	102	103	69	66
	e(C.m 2 )	e 31	e 33	e 15		ε(ε 0 )	ε S 11	ε S 33
	P1 94	-12.4 26.2 20.8		P1 94	2710 2230
	P1 91	-4.9 21.4 17.1		P1 91	1820 1461
	PMN-33PT -3.9 20.3 10.1	PMN-33PT 1434 680

c

Table 4 .

 4 Piezoelectric coupling and permittivity coefficientsAPPENDIX B. INERTIA CALCULATIONThe inertia I xy is defined as I xy = z y z 2 dA with dA = dzdy.For the structure I s =

	ts 2 z=-ts 2	ls 2 y=-ls 2	z 2 dA =	l s t 3 s 12	.
	For the piezoelectric element I p =	ts 2 +tp z= ts 2	lp 2 y=-	lp 2	z 2 dA =	l p 3	t p	3t 2 s 4	+	3t p t s 2	+ t 2 p .


Yves Colin De Verdière
email: yves.colin-de-verdiere@univ-grenoble-alpes.fr

Corentin Le Bihan
email: corentin.le-bihan@ens-lyon.fr

Christophe Bavard

Yves Carrière

On essential-selfadjointness of differential operators on closed manifolds

In all this paper, "smooth" means C ∞ X-1, pp 155-157. However the potentials involved are quite complicated near infinity: they do not admit a polynomial asymptotic behaviour. A classical results in this domain is Gaffney's Theorem which states that, if a Riemannian manifold (X, g) is complete, the Laplace operator on it is ESA. For a clear proof, see , pp 151-152. For a more recent work on this aspect, see .

Introduction

Let X be a compact smooth 1 manifold without boundary equipped with a smooth density |dx|. We will often denote by L 2 the Hilbert space L 2 (X, |dx|). Let P be a differential operator of degree 2 on X with smooth coefficients acting on complex valued functions. We assume in what follows that P is symmetric, i.e., for any pair of smooth complex valued functions f, g on X, we have X P f ḡ|dx| = X f P g|dx|.

The adjoint P of P is then defined as follows: the domain D(P ) is the set of distributions f ∈ L 2 so that P f belongs to L 2 and P is the operator P acting on such distributions. Definition 1.1 A symmetric linear differential operator with smooth coefficients P is essentially self-adjoint (denoted ESA in what follows) if the graph of P in L 2 × L 2 is the closure of the graph of P with domain the smooth functions on X. More explicitely, for each v = P u with u, v ∈ L 2 , there exists a sequence (u n , v n ) with u n smooth, v n = P u n and

(u n , v n ) → (u, v) in L 2 × L 2 . 1
The ESA property is also called Quantum completeness because the evolution equation du/dt = iP u with u(t = 0) = f , with f smooth, has then an unique solution defined for t ∈ R and denoted u(t) = exp(itP )f (see [R-S-75], sec. VIII).

On the other hand, P admits a principal symbol and also a sub-principal symbol: if one chooses local coordinates

x = (x 1 , • • • , x n ) so that |dx| = |dx 1 • • • dx n |, and if P = ∂ ∂x k a kl (x) ∂ ∂x l + b k (x) ∂ ∂x k + c(x)
, the principal symbol is p 2 := -a kl (x)ξ k ξ l and the sub-principal symbol is p 1 := i b k (x)ξ k (see [G-K-L-64] sec. 3). Note that, if P is symmetric, p 1 and p 2 are real valued. We denote p = p 2 + p 1 and call it the symbol of P . Note that p is real valued if P is formally symmetric. The symbol p is independent of the choice of local coordinates as soon as interpreted as a function on the cotangent space T X. The cotangent space is a symplectic manifold and one can use p as an Hamiltonian function on it. We say that P is classically complete if the Hamiltonian flow of p is complete: it means that the maximal interval of definition of any integral curve of the Hamiltonian vector field of p is R.

A natural question is then: how are classical and quantum completeness related? The goal of this article is to give very partial answers to this question. We will state a possible answer as the Conjecture 1.1 Let P be a formally self-adjoint differential operator of degree 2 on C ∞ (X) where X is a closed smooth manifold equipped with a smooth density |dx|, classical and quantum completeness are equivalent.

As we will see, this conjecture holds true in the following cases:

1. Differential operators of degree 2 on the circle of the form

P = a(x)d 2 x + • • •
where all zeroes of a are of finite order.

2. Differential operators of degree 1.

3. Generic and conformally flat Lorentzian Laplacians on 2-tori.

Let us describe other known results on this question: a classical result is that classical completeness and quantum completeness are not equivalent in general; examples of Schrödinger operators on R are given in [R-S-75], section

The case of differential operators on a compact interval

Let X := [α, β] be a compact interval. We consider a differential operator P of degree 2 whose coefficients are smooth up to the boundary. Assuming that

P is symmetric on C ∞ c (]α, β[, |dx|
) (it is usually called formally symmetric), then P is given by the equation (1).

We want to describe the Dirichlet boundary conditions. For that, we assume that P is elliptic near the boundary, i.e. that a does not vanish at the points of ∂X. We will take for domain of P the space

D(P ) := C ∞ ( X, C) ∩ {f |f | ∂X = 0}.
Lemma 2.1 The domain of P is then the set of L 2 functions g so that P g ∈ L 2 , where P is acting on distributions defined in the interior of X and g vanishes on the boundary.

Proof.-We get first that P g ∈ L 2 by looking at

X P f ḡ|dx| with f ∈ C ∞ c (]α, β[).
It follows that g is continuous near the boundary. Then we have, if f ∈ D(P ),

X (P f ḡ -f P g)|dx| = [af g] β α
We have to control the righthandside in terms of the L 2 norm of f which is clearly not possible if g does not vanish on the boundary because a does not vanish at on ∂X.

Localization

Let us prove the following localization Lemma 2.2 Let P be a symmetric operator of degree 2 on a the circle. Let Z ⊂ X be the closed set of points where P is not elliptic. We assume that Z is a finite set

Z = {x 1 , • • • , x N } . Let Ω = ∪ N j=1 [α j , β j
] be a neighbourhood of Z so that P is elliptic near the boundary of Ω. Then P is ESA if and only if the Dirichlet restriction P Ω of P to Ω is ESA.

Proof.-

Let us first prove that, if P Ω is ESA, P is ESA: let us take a ρ ∈ C ∞ c (Ω) with ρ ≡ 1 near Z. Then, if P u = v with u, v ∈ L 2 (X), (1-ρ)u belongs to the Sobolev space H 2 (X) by ellipticity of P on the support of 1-ρ. In particular P ((1-ρ)u) ∈ L 2 . There exists (u n , v n = P u n ) a sequence of smooth functions converging to ((1 -ρ)u, P ((1 -ρ)u) in L 2 by density of C ∞ (X) in H 2 (X). We have now P (ρu) = w with ρu, w ∈ L 2 and support(ρu) ⊂ Ω. ESA of P Ω allows to approximate (ρu, w) by smooth functions (u n , v n = P u n ) and we can assume that u n vanishes near the boundary because ρu does. Then

(u n + u n , v n + v n ) are smooth, converge to (u, v) in L 2 × L 2 and P (u n + u n ) = v n + v n .
This allows to conclude that P is ESA.

Let us now prove that if P is ESA, P Ω is ESA: let us start with

P Ω u = v with (u, v) ∈ L 2 (Ω), u ∈ H 2 near the boundary and u(∂Ω) = 0. We choose ρ ∈ C ∞ c (Ω) with ρ = 1 near Z.
Similarly to the previous argument, we decompose u = ρu + (1 -ρ)u. And by ellipticity of P near ∂Ω we get that (1 -ρ)u belongs to the Sobolev space H 2 0 (Ω) of distributions which are in H 2 (Ω) and vanish at the boundary. By density of smooth functions vanishing at the boundary in the Sobolev space H 2 0 , we get an approximating sequence to ((1 -ρ)u, P ((1 -ρ)u). Now we are left with ρu with P (ρu) ∈ L 2 and we use the fact that P is ESA to get an approximating sequence (u n , P u n ) on X. Choosing ρ 1 ∈ C ∞ c (Ω) with ρ 1 = 1 on the support of ρ, we take u n = ρ 1 u n and we get an approximating sequence for (ρu, P (ρu). This allows to conclude.

Note that the previous localization result extends probably to higher dimensional manifolds, but this extension is much less simple.

3 Essential self-adjointness of differential operators of degree 1

The following property goes back to Friedrichs as cited by .

Lemma 3.1 Let P be a differential operator of degree 1 on a closed manifold X and u ∈ L 2 (X) so that P u ∈ L 2 (X). There exists a sequence u j ∈ C ∞ (X) so that u j → u and P u j → P u both in L 2 (X).

If P is symmetric, this implies that the closure in L 2 ⊕ L 2 of the graph of P restricted to smooth functions is the graph of the adjoint of P . Hence P with domain C ∞ (X) is essentially self-adjoint.

Theorem 3.1 Any symmetric differential operator of degree 1 on a closed manifold is essentially self-adjoint.

This holds in particular for differential operators of the form P :

= i(V + 1 2 div |dx| (V ))
where V is a vector field and div |dx| (V ) := d(ι(V )dx)/dx where ι is the inner product.

In the note [L-19], Nicolas Lerner remarks that this property extends to pseudo-differential operator of degree 1. This is related to our problem because then the Hamiltonian flow is complete at infinity: the Hamiltonian vector field of p is bounded by C ξ and the completeness at infinity follows from Gronwall lemma.

4 Sturm-Liouville operators on the circle

Main result

Any symmetric operator on the circle S 1 = R/Z equipped with the Lebesgue measure |dx| can be written as

P = d x a(x)d x -ib(x)d x -i 1 2 b (x) + c(x) (1)
where d x := d/dx, a, b, c are smooth real valued periodic functions of period 1. We assume always in what follows that the zeroes of a are of finite multiplicities. The symbol p of P is

p = -a(x)ξ 2 + b(x)ξ
The term c(x) plays no role in the essential self-adjointeness, so we will forget it in what follows.

Our main result is Theorem 4.1 For operators P of the previous form, classical completeness of the Hamiltonian flow of p is equivalent to quantum completeness of P .

Our proof consists in describing the properties of a and b leading to classical completeness and to study the quantum completeness in the corresponding cases.

Note that the result in the case where the zeroes of a are non degenerate is also proved using some microlocal analysis in .

Classical completeness

We have the Theorem 4.2 Let p := -a(x)ξ 2 + b(x)ξ where the zeroes of a are of finite multiplicity. Then the Hamiltonian flow of p is complete on T S 1 if and only if the zeroes of a are not simple and b vanishes at these zeroes. Moreover, this flow is complete if and only if it is null complete, i.e. complete when restricted to p -1 (0).

Proof.-Recall that the Hamiltonian differential equation writes

dx/dt = -2a(x)ξ + b(x), dξ/dt = a (x)ξ 2 -b (x)ξ.
The function p is constant along the integral curves. We will denote by (x 0 , ξ 0 ) the data at time 0 of the integral curves that we will consider.

The proof splits into three cases:

1. Assume that a(0) = 0 and b(0) > 0. We will show that the flow is not null complete. Let us look at the set p = 0 near x = 0. This set is the union of the disjoint curves {ξ = 0} and C := {a(x)ξ -b(x) = 0}. The curve C ∩ {x > 0} is oriented by the flow so that x is decaying because then dx/dt = -b(x). Let us start on C, with x 0 > 0 small enough and ξ 0 so that -a(x 0 )ξ 0 + b(x 0 ) = 0. Then -a(x(t))ξ(t) + b(x(t)) = 0 for all t. Along C, we have dx/dt = -b(x). Hence, there exists t 0 > 0 so that x(t 0 ) = 0 and ξ(t 0 ) = +∞. The flow is not null complete: the maximal integral curve is only defined up to t - 0 .

2. Assume that 0 is a non degenerate zero of a and b(0) = 0. Let us start with x 0 = 0 and ξ 0 = 0. We have x(t) = 0 for all t and dξ/dt = a (0)ξ 2 -b (0)ξ. The solution of this differential equation is not defined for all t's because a (0) = 0. The flow is not null complete.

3. Assume now that zeroes of a are degenerate and b vanishes on a -1 (0). We want to prove that the flow is complete. We have, by conservation of p, for any integral curve, there exists E so that -a(x)ξ 2 + b(x)ξ ≡ E.

It follows that ξ stays bounded on any compact interval in x disjoint from a -1 (0). We need only to consider what happens when x(t) comes close to a zero of a, says x = 0.

Let us first assume that x 0 = 0, then x(t) ≡ 0 for all t and dξ/dt = -b (0)ξ. The trajectory is complete.

If x 0 = 0 is close to 0, we get dx/dt = ± -4a(x)E -b(x) 2 = O(|x|).
It follows that x(t) does not reach 0 in finite time and hence the integral curve ie defined for all times.

Simple zeroes of a

We will show in this section that if a(0) = 0 is a simple zero of a then P is not ESA.

Let I := [-α, α] with no other zeroes of a inside I. The point 0 is a regular singular point (see Appendix A) of the differential equation (P -i)u = 0.

The indicial equation writes Ar 2 -iBr = 0 with A := a (0), B = b(0). Hence the solutions of this equation near 0 writes, for x > 0, y(x) = f (x) + x iB/A + g(x) if B = 0 and y(x) = f (x) + g(x) log x if B = 0 with f, g smooth up to x = 0 (see Appendix A) and similarly for x < 0 with x -and log(-x).

Let y + be the unique solution of (P -i)y + = 0, y + (α) = 0, y + (α) = 1 on ]0, α], so that y + satisfies the Dirichlet boundary condition at α. And define y -similarly with y -(-α) = 0. If we extend y + by zero for x < 0, we get a Schwartz distribution Y + and (P -i)Y + is supported by the origine. We have

P Y + = d x d x aY + + d x [a, d x ]Y + -ibd x Y + -ib Y + /2. We check that d x aY + is in L 2
loc . So that (P -i)Y + is near 0 in the Sobolev space H -1 , because d x aY + ∈ L 2 . The derivatives δ (0), • • • of the Dirac distribution are not in H -1 . We have the same result for Y -. It follows that (P -i)Y ± = µ ± δ(0).

Hence there is a non zero linear combination Y of Y + and Y -which satisfies (P -i)Y = 0 and the Dirichlet boundary conditions at ±α. This proves that P I is not ESA and hence P is not ESA by Lemma 2.2.

Degenerate zeroes where b(0) vanishes

Let us assume that all zeroes of a are degenerate. Then if I =]c, d[ is an interval between two zeroes of a and assume a > 0 on I, we will show that there is an explicit unitary map from L 2 (I, |dx|) onto L 2 (R, |dy|) sending C ∞ c (I) (the set of operator with compact support on I) into C ∞ c (R) (set of operator with compact support on R) and sending P to Q = d 2 y + V with an explicit V .

First step: a gauge transform. Let us consider P S := e -iS P e iS where S is smooth and real valued. We get, by an easy calculation,

P S = d x ad x -i(b -2aS )d x -i(b /2 -a S -aS") -aS 2 + bS
Choosing S so that S = b/2a, we get

P S = d x ad x + b 2 /4a
Second step: a change of variable. Let us choose x 0 ∈ I. Let us define y = φ(x) = x x 0 a -1 2 (t)dt. The map φ is smooth diffeomorphism in I onto R. Let us introduce the unitary transform Ω : L 2 (I, |dx|) → L 2 (R, |dy|) defined by Ωf (φ(x)) = a(x) 1/4 f (x). We compute P Ω := ΩP S Ω . We get P Ω = d 2 y + V (y) with

V (y) = b 2 2a + a 2 16a - a 4 (φ -1 (y)) (a 2 /16a + a /4) is bounded. The derivative of (b/a 1/2 ) • φ -1 is b - a b 2a • φ -1
which is also bounded. So (b/a 1/2 ) • φ -1 is bounded by C(1 + |y|) and we get that V (y) ≤ C(y 2 + 1). It follows then from the Farine-Lavis Theorem (Theorem X.38 of [R-S-75]) that P Ω is ESA and hence P is ESA on C ∞ c (I). It follows that P is ESA on C ∞ c (S 1 \ a -1 (0)) and a fortiori on C ∞ (S 1 ).

Degenerate zeroes where b(0) does not vanish

Finally, we study the case where all zeroes of a are degenerate and b does not vanish at least at one of these, say x = 0. We will need the following Lemma 4.1 Let us choose a smooth function E on I :=]0, c] so that E = b/a. There exists two independent solutions u 1 and u 2 of (P -i)u = 0 on I such that u 1 is smooth up to 0, u 2 = u 3 e iE with u 3 smooth up to 0.

It follows that the functions a(x)d x u j are in L 2 and that P is not ESA by the same argument than in Section 4.3.

Proof.-(of Lemma) We check first the existence of u 1 in an elementary way by showing the existence of a full Taylor expansion directly: we start with the Ansatz u 1 (x) = 1+a 1 x+a 2 x 2 +• • • . We get b(0)a 1 +(b (0)/2-1) = 0. Hence a 1 . Then, inductively, we get an expression for a k as a function of the a l for l < k. Applying Malgrange's Theorem 7.1 in [Ma-74], we get a smooth solution u 1 with u 1 (0) = 1.

Then we make the Ansatz u 2 = u 1 v and we get the following differential equation for v:

d x + a a + 2 u 1 u 1 -i b a d x v = 0
It follows that, we can choose

d x v = 1 au 2 1 e iE
If k is the order of the zero of a at x = 0, we can choose local coordinates near 0 so that E = 1/y k-1 , we get d y v = A(y)y -k exp(i/y k-1 ). We can integrate by part and we get

v(y) = A 0 (y)e i/y k-1 - 1 y A 1 (y)e i/y k-1
with A 0 and A 1 smooth. We can iterate the integration by part and get a formal solution v(x) x) . Again we can apply Malgrange's Theorem.

≡ (v 0 +v 1 x+• • • )e iE(
5 Lorentzian Laplacians on surfaces

General facts on Lorentzian tori

We will consider for X a 2-torus with a smooth Lorentzian metric g. Recall that a Lorentzian metric on a surface is a smooth non degenerate symmetric 2-form of signature (1, 1). There is, as in the Riemannian case, an associated geodesic flow (the Hamiltonian flow of the dual metric), a canonical volume form and a Laplace operator, which is an hyperbolic operator.

The null curves: a smooth curve γ is said to be null if, at every point x of γ, and any tangent vector V to γ at the point x, we have g x (V, V ) = 0, i.e. the tangent spaces to γ are isotropic. Locally the null curves are the leaves of two transverse foliations. This is not always true globally.

A closed null leaf γ is a simple closed curve which is null. There is then a neighborhood of γ with two null foliations: F + is close to the tangent space to γ and F -is transversal to γ. We can then define a Poincaré map P γ as follows. Take a point x 0 on γ and a germ of leaf of F -, C at x 0 . Then P γ is a germ of diffeomorphism (C, x 0 ) into itself obtained by following the null leaves of F + . The map P γ is uniquely defined modulo smooth conjugation by germs of diffeomorphisms. Note that P γ is orientation preserving because X is orientable.

Examples of non geodesically complete Lorentzian surfaces

It is known that Lorentzian metrics on the 2-torus are not always geodesically complete. It is the case for example for the Clifton-Pohl torus: Let T be the quotient of R 2 \ 0 by the group generated by the homothety of ratio 2. On T , the Clifton-Pohl Lorentzian metric is g := dxdy/(x 2 + y 2 ). The associated Laplacian

g = (x 2 + y 2 )∂ 2 /∂x∂y is formally self-adjoint on L 2 (T, |dxdy|/(x 2 + y 2 )).
There is also a much simpler example, namely the quotient on (R + x × R y , dxdy) by the group generated by (x, y) → (2x, y/2). The manifold is not closed, but non completeness sits already in a compact region.

It is known these metric are not geodesically complete. What about ESA of g ?

Some results

We will prove a rather general result: Theorem 5.1 1) If the metric g admits a closed null leaf for which the Poincaré section is not tangent to infinite order to the identity, then g is not geodesically complete. Under the same assumptions, g is not ESA.

2) If g is conformal to a flat metric with a smooth conformal factor on a 2-torus, then g is ESA.

Remark 5.1 In the first case, the proof of the null-incompleteness of the geodesic flow is due to Yves Carrière and Luc Rozoy [C-R-94]. We will reprove it.

Note that the conformal class of a Lorentzian metric is determined by the null foliations; hence ESA is a property of these foliations.

For the proof of Theorem 5.1, we will need two lemmas: Lemma 5.1 The null-geodesic completeness is invariant by conformal change.

Proof.-If g = e φ g 0 , the dual metric satisfy g = e -φ g 0 and hence the geodesic flow restricted to g = 0 are conformal with a bounded ratio.

Lemma 5.2 The ESA property is invariant by conformal change.

Proof.-If g u = v, we have also g 0 u = e φ v and e φ v is in L 2 as soon as v is. Hence, if g 0 is ESA, there exists a sequence (u n , w n = g 0 u n ) n∈N converging in L 2 to (u, e φ v) and e -φ v n converges to v.

This proves part 2 of Theorem 5.1.

Normal forms

It is well known and due to Sternberg that a smooth germ of map (R, 0) → (R, 0) whose differential at the origin is in ]0, 1[∪]1, +∞[ is smoothly conjugated to y → λy and hence is the time 1 flow of the vector field µy∂ y with λ = e µ . A similar result hold for more degenerate diffeomorphisms: we assume that g admits a closed null-leaf so that the Poincaré map is of the form P (y) = y + y k R(y) where R(0) = 0 and k ≥ 2. It is proved in (see Theorem 4) Theorem 5.2 Any such map is the flow at time 1 of a vector field V = A(y)∂ y with A ∼ A 0 y k .

Let γ be closed null-leaf of g and U a neighbourhood of γ so that we have the two null foliations F + and F -. We have the: Theorem 5.3 Let γ be a closed leaf whose Poincaré map is P = Id + R with R of order k. There exists coordinates near γ so that the metric g is conformal to g 0 = dx(dy -a(y)dx) with a(y) ∼ a 0 y k , a 0 = 0.

Proof.-Let us parametrize the closed leaf γ by x ∈ R/Z and extend the coordinate x in some neighbourhood U of γ so that the null foliation F - is given by dx = 0. Choose then for y any coordinate in U so that y = 0 on γ. We introduce the differential equation dy/dx = A(x, y) associated to the foliation F + close to γ. Note that A(x, 0) = 0. Let φ x (y) be the flow of this differential equation. The map y → φ 1 (y) is the Poincaré map of γ. By Theorem 5.2, we can choose a vector field a(y)∂y so that the time one flow is the same Poincaré map; and denote by (φ 0 ) x (y) this flow. Let us consider the germ of diffeomorphism near γ defined by

F : (x, y) → x, y = (φ 0 ) x • φ -1
x (y) . The map F sends the integral curves of dy -bdx onto the integral curves of dy -adx and is periodic of period 1 because the time 1 flows are the same. Hence the two null foliation are given respectively by dx = 0 and dy -a(y )dx = 0. The Theorem follows.

Proof of Theorem 5.1, part 1

The idea is to use the normal form which, being invariant by translation in x, allows a separation of variables and hence application of Theorem 4.1.

Let us first prove the null incompleteness. Using the normal form and the conformal invariance of null completeness, we have to study near y = 0 the Hamiltonian h = -η(a(y)η + ξ). The function ξ is a constant of the motion. Let us take initial conditions with y 0 > 0, ξ 0 > 0, and a(y 0 )η 0 + ξ 0 = 0. We have, using that a(y)η + ξ 0 stays at 0, dy/dt = 2a(y)η = -2ξ 0 . Hence y(t) vanishes for a finite time t 0 and, we have then η(t 0 ) = ∞. Null incompleteness follows.

The Lorentzian Laplacian associated to g = dx(dy -a(y)dy) is given by

= ∂ y a(y)∂ y + ∂ 2 xy
Let us look at solutions of u = v of the form u(x, y) = e 2πix v(y) with v compactly supported near 0. We have u(x, y) = e 2πix (∂ y a(y)∂ y + 2iπ∂ y ) v(y)

The operator P := ∂ y a(y)∂ y + 2iπ∂ y is a Sturm-Liouville operator already studied in section 3. P is not ESA. It follows then that there exists v compactly supported near 0 and L 2 so that P v = w ∈ L 2 and there is no sequences (v n , w n = P v n ) converging in L 2 × L 2 to (v, w). The result follows.

Genericity

The goal of this section is to show that, for a generic Lorentzian metric on the 2-torus, there exists at least one null closed curve γ whose Poincaré map is hyperbolic, i.e. the differential of P γ at the point x 0 of γ is not tangent to the identity. It follows that, for a generic metric on the 2-torus, the geodesic flow is not complete and the Lorentzian Laplacian is not ESA. A C ∞ -generic property is a property which holds for an open dense subset of the metrics in the C ∞ topology.

We have the Proposition 5.1 The existence of a closed null hyperbolic curve is a C ∞ -generic property of Lorentzian metrics on 2-tori.

The following argument is due to Etienne Ghys.

Proof.-The openness of the set of metric with a closed null hyperbolic geodesic is evident.

We say that the metric g splits if the null leaves belongs to two distincts foliations F + and F -. We say that the metric g is orientable if there is a smooth non vanishing vector field V on X so that g(V, V ) is strictly positive everywhere. Any orientable metric splits: the two foliations are the boundaries of the connected component C + of the cone g > 0 containing V , indeed using an orientation of X, we choose F + so that the frame generated by F + and V is positively oriented. We now study the two different cases. 1. Case where g splits: the genericity then follows from the fact that having an hyperbolic closed leaf is a generic property for a foliation of a torus (see ), here for F + .

2. Case where g do not split: we introduce in this case a two-fold cover Y of X for which the lift G of the metric g is orientable. This cover Y is equipped with an involution J exchanging the two null foliations. Let us take a null closed curve γ of one of these foliations. Then J(γ) is a null closed curve of the other. They cannot cross: they have the same rotation number, because J is homotopic to the identity. Moreover all intersections have the same sign because both foliations as well as Y are orientable. It follows that the projection of γ onto X is simple. Still having a closed hyperbolic leaf for the foliation F + of G is generic and F -= J(F + ).

Further questions

There are still several open problems in this setting; we see at least five of them:

1. Prove our conjecture 1.1.

2. Describe the self-adjoint extensions in the case of Lorentzian tori in a geometrical way.

3. If we choose a self-adjoint extension, are there interesting spectral asymptotics?

4. Extend to higher dimensional Lorentzian manifolds 5. Extend to pseudo-differential operators of principal type.

A Appendix: Regular singular points of linear differential equations of order two

For this section, one can look at [Co-Le-55] and .

We consider a linear differential equqation We assume that a(0) = 0 and 0 is a zero of finite order k of a. The singular point x = 0 of P is said to be regular if b (resp. c) vanishes at order at least k -1 (resp. k -2) at x = 0. Otherwise 0 is an irregular singular point. If x = 0 is a regular singular point, we introduce the indicial equation: a (k) (0)r(r -1) + kb (k-1) (0)r + k(k -1)c (k-2) , r ∈ C

We call r 1 , r 2 the two roots of the indicial equation. Then the following holds:

• If Im(r 1 -r 2 ) / ∈ Z, there exist two independent solutions of P u = 0 on a small interval ]0, c[ of the form u j = x r j + v j (x)

• If Im(r 2 -r 1 ) ∈ N, we have u 1 = x r 1 + v 1 (x) and u 2 = x r 2 + (v 2 (x) log x + v 3 (x))

where the functions v j are smooth on [0, c[ and v 1 (0) = v 2 (0) = 1.

  P u := (a(x)d 2 x + b(x)d x + c(x))u = 0


Yves Colin De Verdière
email: yves.colin-de-verdiere@univ-grenoble-alpes.fr

A proof of a trace formula by Richard Melrose

The goal of this note is to give a new proof of the wave trace formula proved by Richard Melrose in the impressive paper . This trace formula is an extension of the Chazarain-Duistermaat-Guillemin trace formula (denoted "CDG trace formula" in this paper) to the case of a sub-Riemannian ("sR") Laplacian on a 3D contact closed manifold. The proof uses a normal form constructed in the papers [CHT-18, CHT-21], following the pioneering work , in order to reduce to the case of the invariant Laplacian on the 3D-Heisenberg group. We need also the propagation of singularities results of Victor Ivrii, Bernard Lascar and Richard Melrose .

1 The CDG trace formula

The following result was proved by and refined by (see also and Appendix F for the history of the trace formulae):

Theorem 1.1 Let (M, g) be a closed connected smooth Riemannian manifold and λ 1 = 0 < λ 2 ≤ • • • the spectrum of the Laplace operator, then we have the following equality of Schwartz distributions:

∞ j=1 e it √ λ j = T 0 (t) + γ∈P T γ (t) mod C ∞
where P is the set of periodic geodesics, SingSupp(T 0 ) = {0} and, for γ a periodic geodesic, SingSupp(T γ ) ⊂ {L γ } ∪ {-L γ }, where L γ is the length of γ. Moreover, if γ is non degenerate, T γ (t) = L 0 e iπm(γ)/2 πdet(Id -P γ )

1 2

(t + i0 -L γ ) -1 1 + ∞ j=1 a j (t -L γ ) j
where • L 0 is the length of the primitive geodesic associated to γ

• m(γ) is the Morse index of γ

• P γ is the linearized Poincaré map.

This result gives a nicer proof of the main result of my thesis saying that, in the generic case, the Length spectrum, i.e. the set of lengths of closed geodesics, is a spectral invariant. This formula holds for any elliptic selfadjoint pseudo-differential operator P of degree 1 by replacing the periodic geodesics by the periodic orbits of the Hamiltonian flow of the principal symbol of P and the Morse index by a Maslov index. Our goal is to prove that the same statement holds for sR Laplacians on a closed contact 3-manifold.

What we mean here by 'Melrose's trace formula' is the formulae of theorem 1, interpreted in this contact sR setting. Richard Melrose gave a proof which I found difficult (see Theorem 6.4 in ).

Review of basic facts and notations

For more details on this section, one can look at the paper .

Contact 3D sR manifolds

In what follows, M is a closed (compact without boundary) connected manifold of dimension 3 equipped with a smooth volume form |dq|. We consider also an oriented contact distribution globally defined as the kernel of a non vanishing real valued 1-form α so that α ∧ dα is a volume form. Let also g be a metric on the distribution D = ker α. The "co-metric" g is defined by g (q, p) := p |Dq 2 where the norm is the dual norm of g(q). To such a set of data is associated

• A geodesic flow denoted by G t , t ∈ R: the Hamiltonian flow of

√

g . The geodesics are projections of the orbits of that flow onto M and are everywhere tangent to D. We will often prefer to consider the geodesic flow as the restriction to g = 1 of the Hamiltonian flow of1 2 g . • A Laplacian which is locally given by ∆ = X 1 X 1 +X 2 X 2 where (X 1 , X 2 ) is an orthonormal frame of D and the adjoint is taken with respect to the measure |dq|.

• A canonical choice of a 1-form α g defining D by assuming that dα g restricted to D is the oriented g-volume form on D.

From a famous theorem of Hörmander, we know the Laplacian is subelliptic and hence has a compact resolvent and a discrete spectrum

λ 1 = 0 < λ 2 ≤ • • • with smooth eigenfunctions. It follows from the sR Weyl law #{j|λ j ≤ λ} ∼ Cλ 2 that Trace(e it √ ∆ ) = ∞ j=1 e it √ λ j
is a well defined Schwartz distribution often called the wave trace (see Appendix A for the link with the wave equation). Our goal is to extend the CDG formula to this case.

The form α g defines a Reeb vector field R by the equations α g ( R) = 1 and ι( R)dα g = 0. This vector field admits the following Hamiltonian interpretation: the cone Σ = D ⊥ is a symplectic sub-cone of T M \ 0 1 . We define the Hamiltonian ρ : Σ → R by ρ(α) = α/α g where α ∈ Σ is a covector vanishing on D. The Hamiltonian vector field of ρ is homogeneous of degree 0 and the projection of this field onto M is the Reeb vector field R (see section 2.4 of ).

The 3D Heisenberg group H 3

For this section, the reader could look at Section 3.1 of . We identify the Heisenberg group H 3 with R 3

x,y,z with the group law

(x, y, z) (x , y , z ) = (x + x , y + y , z + z + 1 2 (xy -yx ))
and the Lie algebra is generated by X, Y, Z with

X = ∂ x + 1 2 y∂ z , Y = ∂ y - 1 2 x∂ z , Z = ∂ z
We choose (D, g) by asking that (X, Y ) is an oriented orthonormal basis of D for g and |dq| = |dxdydz|. We have [X, Y ] = -Z and the Reeb vector field is Z. The Laplacian is ∆ 3 = -(X 2 + Y 2 ) and can be rewritten as

∆ 3 = -|Z| X √ Z 2 + Y √ Z 2 (1)
on the complement of the kernel of Z. We write this as

∆ 3 = |Z|Ω
where Ω is an harmonic oscillator with spectrum {2l + 1|l = 0, 1, • • • } (see , prop. 3.1). We will need the following "confining" result (see ):

Lemma 2.1 Given T > 0 , σ 0 ∈ Σ \ 0 and U a conic neighbourhood of σ 0 , there exists a conic neighbourhood V ⊂⊂ U of σ 0 so that ∀t ∈ [-T, T ], G t (V ) ⊂ U .

Speed of propagation

First, we have the following Theorem 3.1 If u is a solution of an sR wave equation,

∀t ∈ R, Support(u(t)) ⊂ B(Support(u(t = 0)) ∪ Support(du/dt(t = 0)), |t|)
where B(A, r) is the closed sR neighbourhood of radius r of A.

This result follows from the Riemannian case by passing to the limits (see Section 3 of ). We will also use the following Theorem due to Victor Ivrii, Bernard Lascar and Richard and revisited by Cyril ): Theorem 3.2 If e(t, q, q ) is the wave kernel of a sR Laplacian whose characteristic manifold Σ is symplectic, i.e. e is the Schwartz kernel of cos(t √ ∆), then W F (e) ⊂ {(q, p, q , p , t, τ )|τ = ± √ g , (q, p) = G ±t (q , p )} ∪ {(q, p, q, p, t, 0)} where G t is the geodesic flow.

The local wave trace for the Heisenberg group

As a preparation, we will prove the local trace formula for the 3D-Heisenberg group H 3 . The Laplacian ∆ 3 commutes with Z. We can hence use a partial Fourier decomposition of L 2 (H 3 ) identifying it with the Hilbert integral

L 2 (R 3 ) = ⊕ R H ζ dζ
where

H ζ := {f |f (x, y, z + a) = e iaζ f (x, y, z)} which is identified to L 2 (R 2
) by looking at the value of f at z = 0. In what follows, we will omit the space H 0 which corresponds to a flat 2D-Euclidian Laplacian. Using this decomposition, the Laplacian rewrites as follows:

∆ 3 = ∞ l=0 (2l + 1) ⊕ R |ζ|K l ζ dζ
where the operator K l ζ is the projector on the l-th Landau level of ∆ ζ , the restriction of ∆ 3 to H ζ , which is a magnetic Schrödinger operator

H ζ := -(∂ x + iζy/2) 2 + (∂ y -iζx/2) 2 on R 2 with magnetic field ζdx ∧ dy. The Schwartz kernel of K l ζ satisfies K l ζ (m, m) = |ζ| 2π (see Appendix B).
Hence the half-wave operator writes

e it √ ∆ 3 = ∞ l=0 ⊕ R e it √ (2l+1)|ζ| K l ζ dζ
and the local distributional trace is given by

Trace e it √ ∆ 3 f = 1 2π R 3 f (q)|dq| ∞ l=0 R e it √ (2l+1)|ζ| |ζ|dζ for f ∈ C ∞ 0 (H 3 )
. This trace can be explicitly computed by using the distributional Fourier transform

∞ 0 e iτ u u 3 du = 6(τ + i0) -4 . We get, for t = 0, Trace e it √ ∆ 3 f = 6 πt 4 ∞ l=0 1 (2l + 1) 2 H 3 f |dq|.
In particular, this trace is smooth outside t = 0 which is consistent with the fact that there is no periodic geodesic in H 3 .

The same result holds for the trace formula microlocalized near Σ:

Proposition 4.1 Let P be a pseudo-differential operator of degree 0, which is compactly supported and so that W F (P )

∩ {ζ = 0} = ∅. Then Trace e it √ ∆ 3 P is smooth outside t = 0. Proof.-We will first prove that Trace e it √ ∆ 3 χ ∆ 3 |Z| 2 , q , with χ ∈ C ∞ 0 (R × H 3 ), which rewrites as 1 2π ∞ l=0 H 3 |dq| R e it √ (2l+1)|ζ| χ 2l + 1 |ζ| , q |ζ|dζ , is smooth outside t = 0. Let us define I l (t) = R e it √ (2l+1)|ζ| χ 2l + 1 |ζ| , q |ζ|dζ .
We split the integral into ζ > 0 and ζ < 0. Looking at the first integral and using the new variable s = ζ/2l + 1, we get

I l (t) = 2(2l + 1) 2 ∞ 0 e it(2l+1)s χ 1 s 2 , q s 3 ds
Observing that the function χ 1 s 2 , q s 3 is a smooth classical symbol supported away of 0, integrations by part give that I l (t) = O (l -∞ ) as well as its derivatives.

Then we introduce

P := 1 2π 2π 0 e itΩ P e -itΩ dt
which is again a pseudo-differential operator. We check that P -P = [Q, Ω] for a pseudo-differential operator Q: if

S t := t 0 e isΩ (P -P )e -isΩ dt,
this holds true with

Q = 1 2πi 2π 0 e -itΩ S t e itΩ dt .
The corresponding part of the trace vanishes, because Ω commutes with ∆ 3 . Then the full symbol of P can be written as

∞ j=0 |ζ| -j p j I ζ , q
with p j compactly supported. This allows to reduce to the first case.

5 The trace formula for compact quotients of H 3

This section can be skipped. It contains an example with a direct derivation of Melrose's trace formula.

Let us give a co-compact subgroup Γ of H 3 . We will prove Melrose's trace formula for M := Γ\H 3 with Laplacian ∆ M which is ∆ 3 , defined in Section 2.2, restricted to Γ-periodic functions. We fix some time T > 0 and will look at the trace formula for |t| ≤ T . We choose χ D a smoothed fundamental domain, ie χ D ∈ C ∞ 0 (H 3 ) with γ∈Γ χ D (γq) = 1. We will denote by

D • • • the integral H 3 χ D • • • .
We start with e M (t, q, q ) = γ∈Γ e 3 (t, q, γq )

where e 3 is the half-wave kernel in H 3 . Because of the finite speed of propagation and the fact that Γ is discrete, we have only to consider a finite sum for the trace:

Trace e it √ ∆ M = D e 3 (t, q, q)|dq| + γ∈Γ\Id, minq d(q,γq)≤T D e 3 (t, q, γq)|dq|
The first term is smooth outside t = 0 while the second one is given by the CDG trace formula.

With more details: let c > 0 be small enough so that the cone C c = {g < cζ 2 } has the property (P) that there is no γ ∈ Γ \ Id with a geodesic from q to γq of length smaller than T starting with Cauchy data in this cone. This is possible thanks to Lemma 2.1. We can hence split the integrals I γ (t) = D e 3 (t, q, γq)|dq| into two pieces ), is equal to 1 near 0 and is supported in ] -c, c[. The first term is smooth by Theorem 3.2 and property (P). The second term corresponds to the elliptic region and hence we use the parametrix for the wave equation given by "FIO"s as given in the CDG trace formula. We get then that the singularities of the wave trace locate on the length spectrum.

I γ (t) = Trace e it √ ∆ 3 τ γ χ D P + Trace e it √ ∆ 3 τ γ χ D (Id -P ) with τ γ (f ) = f •γ -1 and P = ψ(∆ 3 /|Z| 2 ) where ψ belongs to C ∞ o (R
Note that the "heat trace" can be computed from the explicit expression of the spectrum. This is worked out in Appendix E.

Normal forms

In what follows, M is a closed 3D sR manifold of contact type equipped with a smooth volume. We denote by ∆ the associated Laplacian. The proof of the Melrose formula will be done by using a normal form allowing a reduction to the case of Heisenberg.

Classical normal form

Theorem 6.1 Let Σ be the characteristic manifold, i.e. the orthogonal of the distribution with respect to the duality, and let σ 0 ∈ Σ \ 0, then there exists a conical neighbourhood U of σ 0 and an homogeneous symplectic diffeomorphism χ of U onto a conical neighbourhood of (0, 0, 0; 0, 0, 1) in T H 3 , so that g

H 3 • χ = g M .
We use first [Me-84] (Prop. 2.3) (or [CHT-21] (Theorem 2.1)) to reduce to ρI where ρ is the Reeb Hamiltonian and I the harmonic oscillator Hamiltonian: this means that there is an homogeneous canononical transformation χ from a conic neighborhood of σ 0 into Σ σ × R 2 u,v so that ρI • χ = g with I = u 2 + v 2 . Let us denote by ζ the principal symbol of Z in Equation ( 1). Then we use the normal form of (Prop. 6.1.3) to reduce ρ to |ζ| by a canonical transformation. We get then the normal form |ζ|I which is the canonical decomposition of g H 3 used in [CHT-18] (see the principal symbols in Equation ( 1)).

Quantum normal form

This is a 3-step reduction working in some conical neighbourhood C of a point of Σ.

1. Using FIO's associated to χ, we first reduce the Laplacian to a pseudodifferential operator of the form |Z|Ω + R 0 where R 0 is a pseudodifferential operator of degree 0. This step is worked out in Theorem 5.2 of [CHT-18].

2. We can improve the previous normal form so that R 0 commutes with Ω: the cohomological equations {|ζ|I, a} = b (on T H 3 ) where b, vanishing on Σ and homogeneous of degree j can be solved as shown in Appendix D which is an improvement of what is proved in . It follows that we get a normal form ∆ 3 + R 0 with R 0 commuting with Ω: the full symbol of R 0 is independent of the angular part of the (u, v) variables.

3. Using the spectral decomposition of Ω, we get a decomposition ∆ ≡ ⊕ ∞ l=0 (2l + 1)∆ l Π l where the ∆ l 's are pseudo-differential operators of the form

∆ l = |Z| + 1 2l + 1 R 0
and Π l is the projector on the eigenspace of eigenvalue 2l + 1 of Ω.

We can then use a reduction of the pseudo-differential operators |Z| + 1 2l+1 R 0 to |Z| by conjugating by elliptic pseudo-differential operators A l depending smoothly of ε = 1/(2l + 1) and commuting with Ω as in [DH-72], proposition 6.1.4. We have

A -1 l |Z| + 1 2l + 1 R 0 A l ≡ |Z|
where ≡ means modulo smoothing operators in C.

Proof of the Melrose trace formula

Let us fix T > 0 and try to prove Melrose's trace formula for times t ∈ J := [-T, T ]. Let us fix, for each σ ∈ Σ, a conical neighbourhood U σ of σ as in section 6.2. We then take W σ V σ U σ so that, for any z ∈ V σ and any t ∈ J, G t (z) ∈ U σ where G t is the geodesic flow. This is clearly possible using the classical normal form and the Lemma 2.1. We then take a finite cover of Σ by open cones W α := W σα and a finite pseudo-differential partition of unity (χ 0 , χ α (α ∈ B)) so that W F (χ 0 ) ∩ Σ = ∅, χ α = Id in W α and W F (χ α ) ⊂ V α . We have then to compute the traces of cos t √ ∆ χ 0 and cos t √ ∆ χ α . We will prefer to use the wave equation now because the

operator √ ∆ is not a pseudo-differential operator! We know from Theorem 3.2 that, for t ∈ J, W F cos(t √ ∆)χ α is a subset of {(z, z, t, 0)|z ∈ V α } ∪ {(z, G ±t (z), t, τ = ±g (z))|z ∈ V α }.
If u(t) = cos(t √ ∆)χ α u 0 , we have u tt + ∆u = 0, u(0) = χ α u 0 , u t (0) = 0. We can hence use the normal form and denote by ≡ the equality "modulo smooth functions of t ∈ J" to get

Z α (t) := Trace(cos(t √ ∆)χ α ) ≡ Trace(cos(t ∆ 3 + R 0 ) χ α ),
where χ α is the PDO obtained by Egorov theorem when we take the normal form. Then, because R 0 commutes with Ω,

Z α (t) ≡ ∞ l=0 Trace cos t (2l + 1) |Z| + 1 2l + 1 R 0 Π l χ α and Z α (t) ≡ ∞ l=0 Trace A -1 l cos t (2l + 1)|Z| A l Π l χ α Z α (t) ≡ ∞ l=0 Trace cos t (2l + 1)|Z| A l Π l χ α A -1 l .
We can assume that A l is invertible on W F (χ α ) and put

χ l α = A l χ α A -1 l . We get Z α (t) ≡ ∞ l=0 Trace cos t (2l + 1)|Z| A l Π l A -1 l χ l α
Using the fact A l commutes with Ω and hence with Π l , we get finally

Z α (t) ≡ ∞ l=0 Trace cos t (2l + 1)|Z| Π l χ l α Π l
We can then apply a variant of the Proposition 4.1, more precisely of its proof, where P is replaced by ⊕ ∞ l=0 Π l χ l α Π l and using the fact that the A l 's and hence the χ l α too are uniformly bounded pseudo-differential operators. It remains to study the part Z 0 (t) = Trace cos(t √ ∆)χ 0 which involves the elliptic part of the dynamics for which we can use the FIO parametrix as in . More precisely, for t ∈ J, the geodesic flow maps the microsupport of χ 0 away of Σ, therefore there exists a parametrix for U (t)χ 0 given by Fourier integral operators as in the paper [DG-75] and the calculation of the trace therefore follows the same path.

Appendices A Wave and half-wave equations

Let ∆ be a self-adjoint positive sub-elliptic operator on a closed manifold. The wave equation is

∂ 2 u ∂t 2 + ∆u = 0, u(t = 0) = u 0 , ∂u ∂t (t = 0) = v 0
This gives a one parameter group

U (t) = (U 0 (t), U 1 (t)) on L 2 ×L 2 . The trace of U 0 (t) is Z 0 (t) = Trace cos t √ ∆ = ∞ j=1 cos t λ j
One can introduce also the half-wave equation

∂u ∂t = i √ ∆u, u(t = 0) = u 0 . The trace of the half-wave group is Z(t) = ∞ j=1 e it √ λ j . We have the relation Z = 2H(Z 0 )
where H is the L 2 -projector multiplying the Fourier transform by the Heaviside function. It follows that the singularities of both distributions are easily related.

In the elliptic case, one can work directly with the half-wave group because √ ∆ is still an elliptic pseudo-differential operator ). This is no longer the case for sub-elliptic operators.

B The value of K l

ζ (m, m)

Recall that K l ζ is the orthogonal projector on the l-th Landau level with a magnetic field in R 2 equal to ζdx ∧ dy. An easy rescaling shows that

K l ζ (m, m) = |ζ|K l 1 (m, m).
We know from the Mehler formula (see , p. 168) that the heat kernel of the magnetic Schrödinger operator with constant magnetic field equal to 1 is given on the diagonal by e(t, m, m) = 1 4π sinh t On the other hand, we have

e(t, m, m) = ∞ l=0 e -(2l+1)t K l 1 (m, m)
and

1 4π sinh t = 1 2π ∞ l=0 e -(2l+1)t
Identifying both sums as Taylor series in x = e -t gives

K l 1 (m, m) = 1 2π .

C Toeplitz operators

Let Σ be a symplectic cone with a compact basis. Louis Boutet de Montvel and Victor Guillemin associate in to such a cone an Hilbert space and an algebra of operators called the Toeplitz operators with the same properties as the classical pseudo-differential operators. The latter case corresponds to the cone which is a cotangent cone. For an introduction, one can look at .

Two examples are implicitly present in this paper:

1. Harmonic oscillator: the harmonic oscillator Ω = -d 2 x +x 2 is an elliptic self-adjoint Toeplitz operator; the cone Σ is R 2 u,v \0 with the symplectic form du ∧ dv and the dilations λ.(u, v) = ( √ λu, √ λv). The symbol of Ω is u 2 + v 2 .

2. Quantization of the Reeb flow: if Σ ⊂ T X \ 0 is the characteristic cone of our sR Laplacian, one can quantize the Reeb Hamiltonian ρ as a first order elliptic Toeplitz operator of degree 1.

D A cohomological equation

The following proposition is a global formulation of the formal cohomological equations discussed in [CHT-18] (section 5.1 and Appendix C) with a simple proof:

Proposition D.1 We consider the cohomological equation

{|ζ|I, A} = B (2)
where A, B are smooth homogeneous functions in the cone C := {I < c|ζ|} with compact support in q ∈ H 3 . If B is homogeneous of degree j and vanishes on Σ := {I = 0}, Equation (2) admits a solution A homogeneous of degree j -1.

Restricting to ζ = 1, reduces to prove the following Lemma. Any smooth function f in some disk in C admits a Fourier expansion

f (w) = ∞ n=0 f n (|w| 2 )w n + ∞ n=1 g n (|w| 2 ) wn
where the f n 's and the g n 's are smooth2 . We can use this expansion with f = a(z, .) and g = b(z, .). We consider only the sum of powers of w n . The second can be worked out in a similar way. We then put a = ∞ We can solve it, for I = 0, by

a n (z, I) = 2 I 0 -∞ b n (z + u, I)e in2u/I du = 0 -∞ b n (z + Is/2, I)e ins ds
and a n (z, 0) = b n (z, 0)/in. We need to prove that a n is smooth: we check that a n is continuous and then the derivative are given by the same kind of integrals with derivatives of b n as shown by the last expression of a n . The continuity of a n follows from an integrations by parts in the first expression of a n . Now we want to add up the series n a n w n . I was not able to do that directly and will proceed as follows: the sum n a n w n is convergent as a formal series along Σ, because a n w n = O(I n/2 ). Using Borel procedure, we need only to solve our cohomological equation with a flat righthandside. This follows clearly from the expression

a(z, w) = 0 -∞ b(z + It/2, e it w)dt .
In fact only the behaviour as I → 0 could be a problem, but we can divide b by any power of I.

E The heat trace for a compact quotient of

H 3 Let us consider the discrete subgroup Γ = √ 2πZ 2 ×πZ of H 3 identified with R 3
x,y,z as in Section 2.2. The spectrum of the sub-Laplacian defined in section 2.2 on M = H 3 /Γ is the union of the spectrum of the flat torus R 2 / √ 2πZ 2 and the eigenvalues 2m(2l + 1), m ≥ 1, l ≥ 0, with multiplicities 2m. The corresponding part of the complexified heat trace is hence

Z o (z) = ∞ m=1 2m ∞ l=0 e -2m(2l+1)z
with (z) > 0. Summing with respect to l gives

Z o (z) = ∞ m=1 m sinh 2mz
which we rewrite as

Z o (z) = 1 4z ∞ m∈Z 2mz sinh 2mz - 1 4z
The Fourier transform of x sinh x is π 2 1+cosh πξ . Applying Poisson summation formula to the last expression of Z o , we get

Z o (z) = π 2 16z 2 - 1 4z + π 2 4z 2 ∞ n=1 1 1 + cosh π 2 n/z
The first term gives the Weyl law. Each term in the sum w.r. to n is equivalent to π 2 2z 2 e -π 2 n/z We observe that the lengths of the periodic geodesic of M are the numbers 2π √ n. Hence we recover also the length spectrum giving contributions of the order of exp(-L 2 /4z) as in the Riemannian case as proved in . It could be nice to derive an exact formula for the wave trace from our expression of the heat trace. Similarly, the heat trace for the Riemannian Laplacian on M was computed by Hubert Pesce .

F A short history of the trace formulae

The trace formulae were first discovered independently by two groups of physicists: Martin Gutzwiller for a semi-classical Schrödinger operator and Roger Balian & Claude Bloch in a very impressive series of papers for Laplacians in Euclidean domains . In , page 154, the authors suggested already a possible application to the inverse spectral problems3 , an industry which just started at the end of the sixties. From the point of view of mathematics, the Poisson summation formula can be interpreted as a trace formula for the Euclidian Laplacian on flat tori. Similarly, the famous Selberg trace formula [Se-56] (see also ) is a trace formula for the Laplacian on hyperbolic surfaces. Then my thesis , inspired by the work of Balian and Bloch and the Selberg trace formula, uses the complex heat equation for general closed Riemannian manifold. The definitive version, the CDG formula, using wave equation, was discovered by Jacques Chazarain and the tandem Hans . They use the power of the Fourier Integral Operators calculus . See [CdV-07] for a review paper. Later results cover the cases of manifolds with boundaries and semi-classical versions.

Author references: no conflist of interest, no funding.

  w) ∈ R × {|w| < c}, w = |w|e iθ , b smooth, compactly supported in z and I = |w| 2 . We assume that b(z, 0) = 0.Then Equation (3) admits a smooth solution a depending smoothly of b.

  n=1 a n (z, I)w n + a 0 (z, I) and b = ∞ n=1 b n (z, I)w n + Ic 0 (z, I). The factorization b 0 (z, I) = Ic 0 (z, I) follows from the assumption on b. We can take a 0 (z, I) = 2 z -∞ b 0 (s, I)ds. The equation for a n , with n ≥ 1, writes ina n

if Σ = {(q, tα(q)|t ∈ R\0, q ∈ M }, the symplectic form on T M restricts to tdα+dt∧α whose square is the volume form -tdt ∧ α ∧ dα

This is proved first in formal power series in(w, w). This reduces the problem to functions which are flat at the origin. A Fourier expansion in θ allows to finish the proof.

They wrote "The analysis of the eigenvalue density as a sum of oscillating terms gives a new insight into the problem of "hearing the shape of a drum" [Kac paper]. . . . It is convenient, for the discussion to start from the fact that the knowledge of eigenvalues determines uniquely the path generating function . . . Thus . . . the lengths of the closed stationary polygons are determined".


Ramin Mehrabifard

Two-dimensional simulation of Argon dielectric barrier discharge (DBD) in plasma actuator structure with COMSOL Multiphysics 1 Two-dimensional simulation of Argon dielectric barrier discharge (DBD) in plasma actuator structure with COMSOL Multiphysics

Keywords: Dielectric Barrier Discharge, Argon Plasma, Plasma Actuator, Body Force, COMSOL Multiphysics, Gas Discharge

Dielectric barrier discharge (DBD) plasma is used for various applications. DBD is also one of the most efficient and low-cost methods for active fluid flow control. In this study, a detailed physical model of DBD in atmospheric pressure at 1kV DC voltage is developed with COMSOL Multiphysics software. Argon gas is also used as a background gas and electrodes are assumed to be copper. Plasma parameters such as electron and ion density, electric field, potential, and temperature for different distances of electrodes (1.0mm, 0.9mm, 0.8mm) have been investigated. Moreover, the effect of dielectric type (Quartz, Silica Glass, Mica) on these key parameters is investigated. The results of the simulation show that the longitudinal distance of the buried electrodes to the exposed electrodes has a direct influence on parameters such as electron temperature, and electron and ion density which are the main factors of fluid flow control. These parameters have the maximum value when mica is used as a dielectric and the lowest value when silica glass is utilized.

Introduction

Due to its exceptional features, cold atmospheric pressure plasma (CAP) has been applied in a wide range of scientific fields [START_REF] Da Ponte | Plasma deposition of PEO-like coatings with aerosol-assisted dielectric barrier discharges[END_REF][START_REF] Hati | Nonthermal plasma technology and its potential applications against foodborne microorganisms[END_REF][START_REF] Lukes | Aqueous-phase chemistry and bactericidal effects from an air discharge plasma in contact with water: evidence for the formation of peroxynitrite through a pseudo-secondorder post-discharge reaction of H 2 O 2 and HNO 2[END_REF][START_REF] Mehrabifard | Improving Cold Atmospheric Pressure Plasma Efficacy on Breast Cancer Cells Control-Ability and Mortality Using Vitamin C and Static Magnetic Field[END_REF][START_REF] Zimmermann | Effects of cold atmospheric plasmas on adenoviruses in solution[END_REF][START_REF] Zimmermann | Effects of cold atmospheric plasmas on adenoviruses in solution[END_REF][START_REF] Zimmermann | Effects of cold atmospheric plasmas on adenoviruses in solution[END_REF][START_REF] Zimmermann | Effects of cold atmospheric plasmas on adenoviruses in solution[END_REF][START_REF] Zimmermann | Effects of cold atmospheric plasmas on adenoviruses in solution[END_REF][START_REF] Zimmermann | Effects of cold atmospheric plasmas on adenoviruses in solution[END_REF][START_REF] Zimmermann | Effects of cold atmospheric plasmas on adenoviruses in solution[END_REF][START_REF] Zimmermann | Effects of cold atmospheric plasmas on adenoviruses in solution[END_REF]. CAPs can involve partially ionized gas, many active species (positive/negative ions, radicals) UV radiation and, transient electric field [START_REF] Graves | The emerging role of reactive oxygen and nitrogen species in redox biology and some implications for plasma applications to medicine and biology[END_REF][START_REF] Mehrabifard | Effect of non-thermal atmospheric pressure plasma on MDA-MB-231 breast cancer cells[END_REF][START_REF] Weltmann | Atmospheric-pressure plasma sources: Prospective tools for plasma medicine[END_REF]. Controlling the Preprint of: Ramin Mehrabifard, (2023). Two-dimensional simulation of Argon dielectric barrier discharge (DBD) in plasma actuator structure with COMSOL Multiphysics 2 separation of flow with the use of atmospheric cold plasma is considered very promising [START_REF] Jayaraman | Low-Reynolds Number Flow Control Using Dielectric Barrier Discharge-Based Actuators AIAA 2007-3974 Low-Reynolds Number Flow Control Using Dielectric Barrier Discharge Actuators[END_REF][START_REF] Shang | Surface plasma actuators modeling for flow control[END_REF][START_REF] Sohbatzadeh | Theoretical investigation of supersonic flow control by nonthermal DC discharge[END_REF]. The separation of the flow from the surface will occur at high angles of attack in a variety of cases, including diffusers, cars, feathers of turbines and airfoils, and will always be accompanied by losses, which will reduce significantly the performance parameters and efficiency. The cold plasma can control this flow by making the body force produced by the difference in charge particle and the electric field in the system [START_REF] Boeuf | Electrohydrodynamic force in dielectric barrier discharge plasma actuators[END_REF]. Several methods are used to obtain the body force created by the non-equilibrium plasma [START_REF] Abdollahzadeh | Numerical modelling of boundary layer control using dielectric barrier discharge[END_REF]. Different models of atmospheric pressure discharge are also used to create body force including direct discharge, microwave discharge, corona discharge, spark and dielectric barrier discharge [START_REF] Boeuf | Electrohydrodynamic force in dielectric barrier discharge plasma actuators[END_REF][START_REF] Georghiou | Numerical modelling of atmospheric pressure gas discharges leading to plasma production[END_REF][START_REF] Shang | Modeling of ac dielectric barrier discharge[END_REF]. Therefore, the dielectric barrier discharge is one of the most common methods for creating body force [START_REF] Jayaraman | Modeling of Dielectric Barrier Discharge Plasma Actuator[END_REF]. In the plasma actuator, the structure of the electrodes is non-symmetrical that are located two sides of the dielectric. Many researchers have used numerical codes to simulate plasma. Studies have shown that in previous work, due to the use of numerical codes, the number of reactions considered for discharge and surface reactions are limited [START_REF] Soloviev | Surface barrier discharge modelling for aerodynamic applications[END_REF].

In this research, all the characteristics of an electric discharge of argon gas in a twodimensional structure used in a plasma actuator have been evaluated. And the effect of the distances of the electrodes from each other and the dielectric material in which the grounded electrode is buried on it, on the body force main parameters (net charge and electric field magnitude) has been investigated.

Plasma model 2.1. Geometry and Description

The general schematic of the actuator is shown in Fig. 1. Copper electrodes dimensions are ( 0.2 3 mm mm 

) and the dielectric dimensions are ( 0.6 14 mm mm 

). The gas used in this structure is argon, which contains neutral particles, electrons, negative ions, and positive ions. The total space around the structure involves argon. The input voltage leads to the production or consumption of each component. Argon Preprint of: Ramin Mehrabifard, (2023). Two-dimensional simulation of Argon dielectric barrier discharge (DBD) in plasma actuator structure with COMSOL Multiphysics 3 reacts in different forms on the surface, which is among the seven probable reactions that are shown in Table 1.

In most of the studies that have been done in this field, pure gases have been used to produce smaller components and simplify computations. Increasing reactions also result in longer computing time. In addition to the seven reactions mentioned above, two surface reactions are considered, as shown in Table 2. Table1. The fluid model is used in this study. By solving the drift-diffusion equations, the density of electrons and energy will be obtained. The governing equations for electric discharge with the drift-diffusion approximation are as follows:

) 1 ( ) 2 ( . ( . ) e e e e n R u n t           ( .
) .

) 3 ( . . ( . ) n E R u n t                    ( . ). . E n D n                This expression . E    
represents the amount of energy obtained form an electron by electric field. e R is the energy derived from non-elastic collisions calculated by the following equation:

) 4 ( gen en Q Q R S q     e S
is power dissipation, gen Q is the thermal source and q is electron charge. e D is electron diffusion coefficient,   energy mobility, and D  energy distribution coefficient.

)

5 ( 5 3 e e e e D T D T          
We used the Townsend coefficients of the electron source that is calculated by the following equation:

) 6 ( 1 M e j j n e j R x a N     M is number of reactions, j
x the molar fraction of the target species for the j reaction, j a the Townsend coefficient for the reaction j, and n N is the total number of neutral particles. Considering the number p of non-elastic electron collisions, we will have:

) 7 ( 1 p j j n e j j R x a N        Which j  
is the energy dissipation of j reaction. In non-electron species, the above equations are used for mass fraction:

Preprint of: Ramin Mehrabifard, (2023). Two-dimensional simulation of Argon dielectric barrier discharge (DBD) in plasma actuator structure with COMSOL Multiphysics 6 ) 8 (

( . )

.

k k k k w u w j R t            k w
is the ionic density, k j is the energy flux of the ions. The electrostatic field is obtained by the following equation:

) 9 ( 0 .( ) r E      0 
is the permittivity of vacuum, and r  is a relative dielectric constant.

Boundary condition

With respect to the boundary conditions for the electron flux and energy flux, the following relation is established:

) 11 ( ) 12 ( 1 ˆ. ( ) (
                    
The right-hand side of equation ( 11) shows the electron production by the secondary electron and  is the secondary electron coefficient. On the surface of electrodes, ions and excited species are neutralized by the surface reaction. Surface interactions on the electrode are indicated by the j  coefficient, which indicates the probability of the function of the j species. The continuity equation for ions is defined as follows :

) 13 (

.( ) .( )

i i i i i i i i n n u n q D n S t               And i
S indicate the electrostatic potential and the electron density change rate.

Result and discussion

The evolution of plasma characteristics by buried electrode displacement are investigated in this study. The transverse distance to the dielectric edge is 0.2 mm and the longitudinal distance from the electrode are considered 0.8mm, 0.9mm and 1 mm. The boundary conditions for free space are Neumann and, on the electrode, surface are Dirichlet conditions (fig 1). Figure 3 shows the electrical potential in the presence of plasma at different distances of the electrodes. The power electrode voltage is 1000 V (DC) and the buried electrode is at zero potential (ground). The process of ion changes in this structure is more significant with the displacement of electrodes. The maximum electron density in the edge of exposed electrode is Variation in density is more significant for distances of 0.8 mm to 0.9 mm.

Fig6. Electron and ion density on the upper space of the wall.

Figure 7 shows the electron temperature distribution at 10 ns for different distances of the electrodes. The electron temperature near the exposed electrode is the maximum power and by moving away from it to the dielectric surface reaches less than 1 electron volt. The increase in the distance between the electrodes also reduces the electron temperature in the discharge zone.

Fig7. Electron temperature distribution at discharge zone at 10ns (a) 0.8mm (b) 0.9mm (c) 1mm.

One of the effective parameters in creating body force is the electric field. Increasing the electric field apart from the direction is the main factor in increasing the body force. In Fig. 8, the density of the field lines and the field size in the space between the two electrodes and above the dielectric surface are shown. Preprint of: Ramin Mehrabifard, (2023). Two-dimensional simulation of Argon dielectric barrier discharge (DBD) in plasma actuator structure with COMSOL Multiphysics 15

The dielectric material plays an important role in the rate of species production. Here are three different materials for the dielectric. Then the effect of these materials on electron and ion density has been investigated. The horizontal distance between the two electrodes is considered to be 1mm here, and the particle density at the moment of 10 nanoseconds is shown in Figure 10. Since the magnitude of the positive and negative charge difference has a direct effect on the magnitude of the body force, this value has been calculated for different dielectrics Figure 11. The highest charge difference is for mica dielectric and the lowest is for silica glass.

Conclusion

In this study, the electric discharge of Argon gas in the two-dimensional structure of plasma actuators has been studied. The density of species at different intervals of the electrodes was studied. Also, in the discharge of argon gas, seven possible reactions and two surface reactions on the dielectric are considered. The magnitude of the field and net electric charge in the discharge zone are the main factor body force production. The results of the simulation of electrical discharge with Comsol show the significant effect of the buried electrode distance on the body force. Longitudinal displacement of 0.1mm (1.0-0.8mm) distance from the exposed electrode has a significant effect on plasma parameters, including ionic and electron and electric field. Also, as a dielectric, Mica can create the largest charge difference compared to Quartz and Glass.

References.

Fig 1 .

 1 Fig 1. Plasma actuator structure and boundary condition for potential

  expresses the electron continuity equation. e n is the electron density, e D is electron diffusion coefficient, and e  electron flux, and u is average velocity of the species and e R is the electron production rate. Equation (2) is the equation of electron flux, which consists of two parts: drift and diffusion. The electron energy density is calculated by the following equation:

Fig 3 .

 3 Fig 3. Electric potential distribution on the cross-section of the actuator for different distances of electrodes (a) 0.8mm (b) 0.9mm (c) 1mm.

Figures 4

 4 Figures 4 and 5 show the effect of buried electrode displacement on electron and ion density at 10ns. Increasing electrode distance decrease the number of electron and ion on the right side of the exposed electrode. The discharge starts at the right side of the exposed electrode and propagate to the wall over time. The maximum electron density in the edge of the exposed electrode is 18 3 1.4 10 (1/ ) m  .

  Fig5. The ion density at 1kV DC voltage over 10ns (a) 0.8mm (b) 0.9mm (c) 1mm. For more investigations of the electron and ion density change due to the displacement of the buried electrode, electron and ion densities are shown on the upper side of the wall (fig 6). Fig 6 shows ion and electron density though a virtual line from points(5, 1.1) to (9 and 1.1) for different displacements of electrodes. Variation in density is more significant for distances of 0.8 mm to 0.9 mm.

Fig 8 .

 8 Fig 8. Density and magnitude of the electric field in plasma actuator (a) 0.8mm (b) 0.9mm (c) 1mm.

Fig 9 .

 9 Fig 9. Electric field arrow made by electrode potential and charges (a) 0.8mm (b) 0.9mm (c) 1mm.

  fig10. Ion (up) and electron (down) density in presence of different dielectric materials.

  Table of collisions and reactions modeled.

	Reactions	Formula			Type	 	( ) eV
	1	e Ar e Ar   		Elastic	0
	2	e Ar e Ars   		Excitation	11.5
	3	e Ars e Ar   		Superelastic	-11.5
	4	e Ar 	2   e Ar 		Ionization	15.8
	5	e Ars 	2   e Ar 		Ionization	4.24
	6	Ars Ars e Ar Ar     	Penning ionization	-
	7	Ars Ar Ar Ar   	Metastable quenching	-
	Table2. Table of surface reactions
	Reaction	Formula	Sticking Coefficient
	1		Ars Ar 	1
			Ar	 	Ar	1

Preprint of:Ramin Mehrabifard, (2023). Two-dimensional simulation of Argon dielectric barrier discharge (DBD) in plasma actuator structure with COMSOL Multiphysics