Andreea Grigoriu
email: aglachap@princeton.edu

Implicit Lyapunov control for Schrödinger equations with dipole and polarizability term

We analyze in this paper finite dimensional closed quantum systems in interaction with a laser field. The characteristic of the problem is that the interaction laser-system is modeled by a first order term (dipole coupling) and a second order term (polarisability coupling), that appear in the Hamiltonian of the Schrödinger equation. In order to determine the control, an implicit Lyapunov trajectory tracking procedure is applied when there is no direct coupling between the target state and the eigenvectors of the internal hamiltonian. The method is applied for the difficult case of degenerate systems too. The controlled Lyapunov function is defined by an implicit equation and its existence is shown by a fix point theorem. The convergence is analysed using the LaSalle invariance principle. The performance of the feedbacks is illustrated by numerical simulations.

I. INTRODUCTION

The control of molecules dynamics at the quantum level is a core issue for numerous applications. This is why quantum control constitutes today a very active research field, both from the theoretical and experimental point of view ( [START_REF] Assion | Control of chemical reactions by feedbackoptimized phase-shaped femtosecond laser pulses[END_REF], [START_REF] Brumer | Coherent chemistry: Controlling chemical reactions with lasers[END_REF], [START_REF] Levis | Selective bond dissociation and rearrangement with optimally tailored, strong-field laser pulses[END_REF], [START_REF] Weinacht | Controlling the shape of a quantum wavefunction[END_REF] etc.). Before presenting the issues addressed in this paper, we first introduce the time-dependent Schrödinger equation, which models the time-evolution of a finite dimensional closed quantum system under the influence of a laser field:

i d dt Ψ(t) = (H 0 + ǫ(t)µ 1 + ǫ 2 (t)µ 2 )Ψ(t), (1)
where H 0 , µ 1 and µ 2 are N × N Hermitian matrices with complex entries. The wavefunction Ψ is a vector in C N , satisfying N j=1 |Ψ j | 2 = 1, thus it lives on the unit sphere S N (0, 1) of C N and ǫ(t) ∈ R is the control (for example the laser intensity).

The term H(t) = H 0 +ǫ(t)µ 1 +ǫ 2 (t)µ 2 is the Hamiltonian of the system and H 0 is the internal Hamiltonian. This last one characterizes the system in the absence of the laser. The matrices µ 1 and µ 2 describe the interaction of the quantum system with the laser field. This type of model is used when the first order term, ǫ(t)µ 1 , also called dipole coupling [START_REF] Rabitz | Optimal control of selective vibrational excitation in harmonic linear chain molecules[END_REF] does not have enough influence on the system to reach the control goal; the goal may become accessible only after adding the polarizability term ǫ 2 (t)µ 2 (see e.g. [START_REF] Dion | Two-frequency IR laser orientation of polar molecules. Numerical simulations for HCN[END_REF], [START_REF] Dion | Laser-induced alignment dynamics of HCN: Roles of the permanent dipole moment and the polarizability[END_REF] and related works).

One of the problems that must be clarified for the nonlinear control equation [START_REF] Altafini | Controllability of quantum mechanical systems by root space decomposition of su(n)[END_REF] is to find if a final time T and a laser pulse ǫ(t) exist such that ǫ(t) is able to steer the system from an initial state Ψ 0 to some arbitrary predefined target Ψ(t = T ) = Ψ target . If the answer to this question is positive the systems is called controllable. Positive results of controllability are obtained by applying the Lie algebra criteria [START_REF] Brockett | Lie theory and control systems defined on spheres[END_REF], [START_REF] Sussmann | Controllability of nonlinear systems[END_REF] and especially the specific result in [START_REF] Turinici | Beyond bilinear controllability: applications to quantum control[END_REF]. This study is detailed in [START_REF] Coron | Quantum control design by Lyapunov trajectory tracking for dipole and polarizability coupling[END_REF].

The theoretical results of controllability do not offer automatically a method to determine the laser field. Very often this task is formulated as a cost functional to be minimized. Several techniques have been developed in this direction: iterative critical point methods (e.g. monotonic algorithms [START_REF] Lapert | Monotonically convergent optimal control theory of quantum systems under a nonlinear interaction with the control field[END_REF], [START_REF] Ohtsuki | Monotonically convergent algorithms for solving quantum optimal control problems of a dynamical system nonlinearly interacting with a control[END_REF], [START_REF] Ohtsuki | Monotonically convergent algorithms for solving quantum optimal control problems described by an integrodifferential equation of motion[END_REF], [START_REF] Ohtsuki | Generalized monotonically convergent algorithms for solving quantum optimal control problems[END_REF], [START_REF] Maday | New formulations of monotonically convergent quantum control algorithms[END_REF], [START_REF] Tannor | Control of photochemical branching: Novel procedures for finding optimal pulses and global upper bounds[END_REF], [START_REF] Salomon | A monotonic method for solving nonlinear optimal control problems[END_REF], [START_REF] Schirmer | Efficient algorithm for optimal control of mixed-state quantum systems[END_REF], [START_REF] Zhu | A rapid monotonically convergent iteration algorithm for quantum optimal control over the expectation value of a positive definite operator[END_REF] ), iterative stochastic techniques (e.g., genetic algorithms [START_REF] Turinici | Eficient algorithms for the laboratory discovery of optimal quantum controls[END_REF]), trajectory tracking methods or local control procedures [START_REF] Chen | Competitive tracking of molecular objectives described by quantum mechanics[END_REF], [START_REF] Gross | Inverse quantum-mechanical control: A means for design and a test of intuition[END_REF], [START_REF] Kosloff | Wavepacket dancing: Achieving chemical selectivity by shaping light pulses[END_REF], [START_REF] Lu | Unified formulation for control and inversion of molecular dynamics[END_REF], [START_REF] Sugawara | General formulation of locally designed coherent control theory for quantum systems[END_REF], [START_REF] Ohtsuki | Quantum control of nuclear wave packets by locally designed optimal pulses[END_REF], [START_REF] Sugawara | Control of quantum dynamics by a locally optimized laser eld. application to ring puckering isomerization[END_REF], [START_REF] Mirrahimi | Lyapunov control of bilinear Schrödinger equations[END_REF], [START_REF] Beauchard | Implicit Lyapunov control of finite dimensional Schrödinger equations[END_REF], [START_REF] Mirrahimi | Trajectory tracking for quantum systems: a Lyapounov approach[END_REF].

For control quantum systems evolving according to equation (1), Lyapunov trajectory tracking techniques [START_REF] Altafini | Controllability of quantum mechanical systems by root space decomposition of su(n)[END_REF], [START_REF] Constantinescu | Parametrizing quantum states and channels[END_REF], [START_REF] Ramakrishna | Control of a coupled two-spin system without hard pulses[END_REF] have been applied. One of the advantages of these approaches is that an explicit formula for the laser field is obtained. The method consists in defining a function V , which is nonnegative and vanishes when Ψ = Ψ target :

V (Ψ, t) = Ψ -Ψ target |Ψ -Ψ target = Ψ -Ψ target 2 , (2)
where .|. denotes the Hermitian product and Ψ target a reference trajectory of (1). Imposing the condition dV /dt 0, V has the properties of a Lyapunov function and we obtain an explicit formula for the control.

In order to make sure that the target Ψ target will be reached using the laser thus obtained a convergence analysis must be provided. An initial positive result of asymptotic stability has been proved in [START_REF] Grigoriu | Lyapounov control of Schödinger equations: beyond the dipole approximation[END_REF] under the hypothesis:

H 1 : H 0 is non degenerate i.e.: λ i = λ j , ∀ i, j ∈ {1, ..., N } with i = j; (λ i ) i=1,2,...,N
are the eigenvalues of H 0 H 2 : direct coupling, through µ 1 , between the target state φ (first eigenvector of H 0 associated to the eigenvalue λ ∈ R; H 0 φ = λφ) and all other eigenvectors i.e.: µ 1 φ j |φ = 0, for every j ∈ {2, .., N }; φ, φ 2 , ..., φ N form an orthonormal system of eigenvectors of H 0 , corresponding to the eigenvalues (λ i ) i=1,2,...,N The explicit formula for the control, that proves to be efficient for the above cases is

ǫ = -kI 1 /(1 + kI 2 ) (3)
with k a constant that belongs to ]0, 1 µ2 [, I 1 = Im µ 1 Ψ|φ and I 2 = Im µ 2 Ψ|φ (see [START_REF] Grigoriu | Lyapounov control of Schödinger equations: beyond the dipole approximation[END_REF], [START_REF] Coron | Quantum control design by Lyapunov trajectory tracking for dipole and polarizability coupling[END_REF]). When some of the (direct) coupling is realized by µ 2 instead of µ 1 i.e: H 3 : µ 1 φ j |φ = 0 or µ 2 φ j |φ = 0, for every j ∈ {2, .., N } the previous feedback (3) is not efficient anymore. Numerical tests in [START_REF] Grigoriu | Lyapounov control of Schödinger equations: beyond the dipole approximation[END_REF], [START_REF] Coron | Quantum control design by Lyapunov trajectory tracking for dipole and polarizability coupling[END_REF] show that the target it is not reached, the systems gets trapped into an unknown state. In order to take into consideration the problems that appear, two alternatives have been proposed: discontinuous feedback and time varying feedback (see [START_REF] Coron | Quantum control design by Lyapunov trajectory tracking for dipole and polarizability coupling[END_REF] for more details). These allow to approximately asymptotically stabilize the system around the target φ.

There are systems for which the target φ is not directly connected by µ 1 or µ 2 with all the other eigenvectors and we are not longer in the hypothesis H 2 or H 3 . An example is provided by the following 5-level quantum system [START_REF] Tersigni | On using shaped light pulses to control the selectivity of product formation in a chemical reaction: An application to a multiple level system[END_REF] with the matrices H 0 , µ 1 and µ 2 defined by:

H 0 =    1.0 0 0 0 0 0 1.2 0 0 0 0 0 1.3 0 0 0 0 0 1.4 0 0 0 0 0 2.15    , µ 1 =    0 0 0 1 1 0 0 0 1 1 0 0 0 0 0 1 1 0 0 0 1 1 0 0 0    , µ 2 =    0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 1 0 0 0 0 1 0 0   . (4)
In Fig. 1, we represent the diagram of the coupling achieved by µ 1 and µ 2 between the eigenvectors of the internal Hamiltonian H 0 of the system (4). One can note that the first eigenvector φ 1 = φ = (1, 0, ..., 0) is directly coupled only with φ 4 and φ 5 ( µ 1 φ|φ 4 = 0, µ 1 φ|φ 5 = 0). There is no direct coupling by µ 1 or µ 2 with φ 2 and φ 3 . In this case the target is φ = (1, 0, 0, 0, 0). Another category of systems is the one with degenerate internal Hamiltonian (hypothesis H 1 is not fulfilled). The common characteristic of all these systems is that there are still controllable and this constitutes a very strong motivation to their study. In this context the main goal of this paper is to introduce a Lyapunov trajectory tracking technique that permits to determine efficient laser fields for this cases too. To this end we will adapt the implicit Lyapunov method introduced in [START_REF] Beauchard | Implicit Lyapunov control of finite dimensional Schrödinger equations[END_REF] for bilinear cases

H(t) = H 0 + ǫ(t)µ 1 .
The general idea is to track a mobile target φ β , instead of a fixed target φ, where β implicitly depends on the state of the system Ψ. As we impose a slow convergence of φ β to φ, we stabilize faster the state of the system around the vector φ β applying a Lyapunov method as the one described by relation [START_REF] Assion | Control of chemical reactions by feedbackoptimized phase-shaped femtosecond laser pulses[END_REF]. This should guarantee that the system is stabilized around φ (see Fig. 2) and it does not get trapped into another unknown state. The balance of the paper is as follows: in section II, we present the implicit Lyapunov technique. An existence result for the implicit Lyapunov function is proved followed by the convergence of the stabilization technique. Section III addresses the case of a degenerate target state. Finally, in section IV , we perform some numerical simulations for a five-dimensional test case followed by conclusions.

II. IMPLICIT LYAPUNOV TRAJECTORY TRACKING

A. Lyapunov function

Recall that two wave functions Ψ 1 and Ψ 2 that differ by a phase θ(t) ∈ R, i.e. Ψ 1 = exp(iθ(t))Ψ 2 , describe the same physical state. To take into account the property we add a fictitiouss control ω (see also [START_REF] Mirrahimi | Lyapunov control of bilinear Schrödinger equations[END_REF]). Hence we replace the evolution equation (1) by:

i d dt Ψ(t) = (H 0 + ǫ(t)µ 1 + ǫ 2 (t)µ 2 + ω(t))Ψ(t), (5)
where ω ∈ R is a new control . We can choose it arbitrarily without changing the physical quantities attached to Ψ.

In the following, for a given β ∈ R, we denote by (λ j β ) 1 j N the eigenvalues of the matrix H(β) = H 0 + βµ 1 + β 2 µ 2 , and by (φ j β ) 1 j N the associated normalized eigenvectors:

(H 0 + βµ 1 + β 2 µ 2 )φ j β = λ j β φ j β . (6)
For simplicity we denote

φ β = φ 1 β and λ β = λ 1 β
We introduce the function V (Ψ) defined by (2) with

Ψ target = φ β(Ψ) : V (Ψ, t) = Ψ -φ β(Ψ) |Ψ -φ β(Ψ) = Ψ -φ β(Ψ) 2 (7)
where the function Ψ → β(Ψ) is implicitly defined as below:

β(Ψ) = Γ( Ψ -φ β(Ψ) 2 ). (8)
The properties of the function Γ are formulated in the following lemma in order to assure the existence of function β. Lemma 2.1: Suppose that hypothesis H 1 holds and consider a continuously differentiable function Γ : R + → [0, β * ] satisfying :

Γ(0) = 0, Γ(s) > 0 for every s > 0 Γ ′ ∞ < 1 8C , (9)
with

C = 1 + max dφ β dβ β=β0 ; β 0 ∈ [0, β * ] , (10)
where φ β is the first eigenvector of the matrix H(β) and

dφ β dβ β=β0
is the derivative of φ β at the point β = β 0 .

Then there exists a unique map β ∈ C ∞ (S N (0, 1); [0, β * ]), such that for every Ψ ∈ S N (0, 1) :

β(Ψ) = Γ( Ψ -φ β(Ψ)
2 ) with β(φ) = 0.

(11) Proof: The proof follows the ideas of [START_REF] Beauchard | Implicit Lyapunov control of finite dimensional Schrödinger equations[END_REF] but needs a minor adaptation to take into account the second order control term ǫ 2 (t)µ 2 and the different definition of the Lyapunov function [START_REF] Assion | Control of chemical reactions by feedbackoptimized phase-shaped femtosecond laser pulses[END_REF]. For completeness we give however the full proof below.

Since we made the assumption that all the eigenvalues of the internal Hamiltonian H 0 are different, we are in the settings of Theorem XII.1 and Theorem XII.4 of [START_REF] Reed | Methods of Modern Mathematical Physics IV: Analysis of Operators[END_REF]. This implies that the eigenvalue and the eigenvectors of the matrix H(β) are analytic functions of β and This implies:

C = 1 + max dφ β dβ β=β0 , β 0 ∈ [0, β * ] < ∞. ( 13
)
• Uniqueness of function β. In order to prove the uniqueness of the function β let us consider the function Γ defined by:

β -→ Γ( Ψ -φ β(Ψ) 2 )
with β ∈ [0, β * ]. For every β 1 and β 2 from [0, β * ], there exists β ∈]β 1 , β 2 [ such that:

Γ( Ψ -φ β1(Ψ) 2 ) -Γ( Ψ -φ β2(Ψ) 2 ) = d dβ Γ( Ψ -φ β(Ψ) 2 ) β=β β 1 -β 2 . ( 14
) Since d dβ Γ( Ψ -φ β(Ψ) 2 ) = -2Γ ′ ( Ψ -φ β1(Ψ) 2 )Re dφ β dβ |Ψ -φ β (15)
and |Re(

dφ β dβ |Ψ )| ≤ | dφ β dβ |Ψ | ≤ ( dφ β dβ )| β=β0 , together with Γ ′ < 1
8C we obtain that the function Γ( Ψ-φ β(Ψ)

2 ) is contracting for fixed Ψ ∈ S N (0, 1). Thus, for any fixed Ψ ∈ S N (0, 1) there exists a unique β(Ψ) ∈ [0, β * ] that verifies [START_REF] Coron | Quantum control design by Lyapunov trajectory tracking for dipole and polarizability coupling[END_REF].

• Existence of function β.

In order to prove the existence of β, let us consider:

F (Ψ, β) = β -Γ( Ψ -φ β(Ψ) 2 ). ( 16
)
This application is C ∞ with respect to Ψ and β, and for fixed Ψ ∈ S N (0, 1) we have:

F (Ψ, β(Ψ)) = 0. (17)
Moreover from relation [START_REF] Kosloff | Wavepacket dancing: Achieving chemical selectivity by shaping light pulses[END_REF] we have:

d dβ F (β, Ψ) = 1 + 2Γ ′ ( Ψ -φ β(Ψ) 2 )Re dφ β dβ |Ψ -φ β , (18)
which is non zero since Γ ′ < 1 8C , with C defined by (10) on the interval [0, β * ]. We are in the hypothesis of the implicit function theorem, this implies the existence of the application Ψ → β(Ψ) that belongs to C ∞ (S N (0, 1); [0, β * ]). Now, we can focus on finding feedback controls such that V is a Lyapunov function. To this end we compute the derivative of V along trajectories of (5) and impose the condition dV /dt ≤ 0. We have:

dV dt = 2(ǫ -β)Im( µ 1 Ψ(t)|φ β ) + 2(ǫ 2 -β 2 )Im( µ 2 Ψ(t)|φ β ) + 2(ω + λ β )Im( Ψ(t)|φ β ) - 2 βRe dφ β dβ |Ψ -φ β , (19)
where Im denotes the imaginary part and Re the real part. For convenience we denote:

I β 1 = Im( µ 1 Ψ(t)|φ β ) and I β 2 = Im( µ 2 Ψ(t)|φ β ).
A simple computation leads to:

β = Γ ′ (V ) 2(ǫ -β)I β 1 + 2(ǫ 2 -β 2 )I β 2 + 2(ω + λ β )Im( Ψ(t)|φ β ) - 2 βRe dφ β dβ |Ψ -φ β . ( 20
)
It follows that:

β = Γ ′ (V )((ǫ -β)I β 1 + (ǫ 2 -β 2 )I β 2 1 + 2Γ ′ (V )Re dφ β dβ |Ψ -φ β + 2(ω + λ β )Im( Ψ(t)|φ β )) 1 + 2Γ ′ (V )Re dφ β dβ |Ψ -φ β . ( 21
)
Remark 2.1: Since Γ ′ (V ) < 1/(8C * ) with C * defined by (10) we have:

1 + 2Γ ′ (V )Re dφ β dβ |Ψ -φ β = 0, (22)
moreover

2Γ ′ (V )Re dφ β dβ |Ψ -φ β < 1 2 . ( 23
)
If we replace relation [START_REF] Mirrahimi | Trajectory tracking for quantum systems: a Lyapounov approach[END_REF] in [START_REF] Maday | New formulations of monotonically convergent quantum control algorithms[END_REF] we obtain:

dV dt = 2((ǫ -β)I β 1 + (ǫ 2 -β 2 )I β 2 + 2(ω + λ β )Im( Ψ(t)|φ β )) × 1 - g(t) 1 + g(t) (24)
with

g(t) = 2Γ ′ (V )Re dφ β dβ |Ψ -φ β . (25)
From Remark 2.1 we have 1 -g(t)

1+g(t)
0. If we denote v = ǫβ and we take:

v = -k(I β 1 + 2βI β 2 )/(1 + kI β 2 ) ω = -λ β -cIm( Ψ(t) | φ β ), (26)
with k and c strictly positive parameters, one gets dV /dt ≤ 0, i.e. V is nonincreasing. Thus we obtain the following feedback control:

ǫ = β -k(I β 1 + 2βI β 2 )/(1 + kI β 2 ) ω = -λ β -cIm( Ψ(t) | φ β ). ( 27
)
Remark 2.2: To guarantee that the denominator 1+kI

β 2 > 0, one notes that |I β 2 | ≤ | µ 2 Ψ(t)|Ψ β | ≤ µ 2 ; therefore 1 + kI β
2 > 0 as soon as k < 1 µ2 . From now on, unless otherwise specified, this condition will be assumed.

B. Convergence analysis

In the following we prove the convergence of the trajectories Ψ of the system (5) toward the target φ. The idea is to use that the feedback presented previously [START_REF] Ramakrishna | Control of a coupled two-spin system without hard pulses[END_REF] for Hamiltonian H(t) = H 0 + βµ 1 + β 2 µ 2 assures the convergence towards the set Z β = {φ β } for every β ∈ ]0, β * ]. So we let β tending to zero and by construction the convergence towards Z β must be faster than the of β towards zero.

Theorem 2.1: Assume that the hypothesis H 1 and H 5 hold, where

H 5 : denote J 1 = {j j = 1, µ 1 φ j β |φ β ) = 0} and J 2 = {j j = 1, µ 2 φ j β |φ β ) = 0}, J 1 ∪ J 2 = {2, 3, . . . , n} and J 1 ∩ J 2 = ∅ on [0, β * ].
Consider [START_REF] Brumer | Coherent chemistry: Controlling chemical reactions with lasers[END_REF] with Ψ ∈ S N (0, 1) and an eigenvector φ ∈ S N (0, 1) of H 0 associated to the eigenvalue λ. If we take the feedback [START_REF] Ramakrishna | Control of a coupled two-spin system without hard pulses[END_REF] with k < 1 µ2 and c > 0, then the limit set of Ψ(t) is reduced to ±φ.

Proof: Up to a shift on ω and H 0 , we may assume that λ = 0. Since hypothesis H 1 holds we can apply Theorem XII.1 of [START_REF] Reed | Methods of Modern Mathematical Physics IV: Analysis of Operators[END_REF]. This implies that hypothesis H 5 is fulfilled.

LaSalle's principle (see, e.g., [START_REF] Khalil | Nonlinear Systems[END_REF]Theorem 3.4, page 115]) guarantees that the trajectories of the system (5) converge to the largest invariant set contained in dV /dt = 0, this implies

that V is constant, i.e. V (Ψ) = V . Since by definition β = Γ(V ), β(Ψ) is constant i.e. β = β c .
The equation dV /dt = 0 means that:

I β 1 + 2β c I β 2 = 0, Im( Ψ(t)|φ βc ) = 0, (28)
Since the Ω-limit set is also invariant under the flow generated by (5) it follows, taking into account [START_REF] Reed | Methods of Modern Mathematical Physics IV: Analysis of Operators[END_REF], that this set consists in fact of trajectories of the system:

i d dt Ψ = (H 0 + β c µ 1 + β 2 c µ 2 )Ψ. (29)
The solutions of (29) have the form:

Ψ = N j=1
b j e -itλ j βc φ j βc .

The eigenvectors

(φ j βc ) j∈{1,...,N } of H(β c ) = H 0 + β c µ 1 + β 2
c µ 2 can be chosen to form an orthonormal basis. Moreover, the orthonormal eigenvectors are holomorphic functions of β c (see [START_REF] Kato | Perturbation Theory for Linear Operators[END_REF], page 121).

If β c = 0 we have Γ( V ) = 0 which implies V = 0. In this case the limit set of Ψ only contains φ.

If β c = 0, on the contrary we substitute relation (30) into relation [START_REF] Reed | Methods of Modern Mathematical Physics IV: Analysis of Operators[END_REF] and we obtain:

Im( Ψ(t)|φ βc ) = Im(b 1 ) φ βc , φ βc + N j=2 Im(b j φ j βc , φ βc e -itλ j βc ). = 0. (31)
and

I β 1 + 2β c I β 2 = Im(b 1 ) µ 1 φ βc , φ βc + j∈J1 Im(b j µ 1 φ j βc , φ βc e -itλ j βc ) +2β c Im(b 1 ) µ 2 φ βc , φ βc + k∈J2 Im(b j µ 2 φ j βc , φ βc e -itλ j βc ) = 0. (32)
From equation (31), together with φ j βc , φ βc = 0 for all j = 2, ..., N we obtain that Im(b 1 ) = 0. Moreover, equation (32) becomes:

I β 1 + 2β c I β 2 = j∈J1 Im(b j µ 1 φ j βc , φ βc e -itλ j βc ), +2β c k∈J2 Im(b j µ 2 φ j βc , φ βc e -itλ j βc ) = 0. (33)
We use hypothesis H 5 together with β c = 0 and we obtain b j = 0 for every j ∈ J (see [START_REF] Coron | Quantum control design by Lyapunov trajectory tracking for dipole and polarizability coupling[END_REF] for more details ). This implies that the limit set only contains ±φ βc . We let β c tend to zero and conclusion follows.

III. DEGENERATE CASES

For systems with degenerate internal Hamiltonian positive numerical tests have been performed using discontinuous and time varying controls obtained by applying "explicit" Lyapunov tracking techniques (see [START_REF] Coron | Quantum control design by Lyapunov trajectory tracking for dipole and polarizability coupling[END_REF]). Under more restrictive hypothesis than for non degenerate cases convergence might be obtained. The advantage of applying an implicit Lyapunov technique is that adding a small perturbation of the form βµ 1 + β 2 µ 2 , with β ∈ [0, β * ] the degeneracy of H 0 is withdrawn and a asymptotic stability result as the one in Theorem 2.1 is obtained. Moreover we conserve the same type of hypothesis as in the non-degenerate case.

Theorem 3.1: Assume that the hypothesis H 4 and H 5 hold. Consider [START_REF] Brumer | Coherent chemistry: Controlling chemical reactions with lasers[END_REF] with Ψ ∈ S N (0, 1) and an eigenvector φ ∈ S N (0, 1) of H 0 associated to the eigenvalue λ. If we take the feedback [START_REF] Ramakrishna | Control of a coupled two-spin system without hard pulses[END_REF] with k < 1 µ2 and c > 0, then the limit set of Ψ(t) reduces ±φ.

Proof: It follows the same steps as in Theorem 2.1. Before being able to apply Theorem 3.1, some details have to be discussed. In the above section the space generated by any eigenvector of the internal Hamiltonian H 0 is always of dimension one. Therefore we were tracking without loss of generality the first eigenvector φ of H 0 . This is no longer the case for a degenerate situation. We may try to stabilize the system around an arbitrary eigenvector φ k , which can generate a space of dimension larger than 1. As a consequence, first we need to recall a result from the perturbation theory for finite dimensional Hermitian operators ( [START_REF] Kato | Perturbation Theory for Linear Operators[END_REF], page 121).

Lemma 3.1: Let us consider the N × N dimensional complex matrices H 0 , µ 1 , µ 2 and let us take

H(β) = H 0 + βµ 1 + β 2 µ 2 . ( 34
)
For each real β, there exists an orthonormal basis (φ n (β)) n∈{1,...,N } of C N consisting of eigenvectors of H(β). Moreover, the orthonormal eigenvectors can be chosen as holomorphic functions of β.

In order to track any eigenvector φ k of a degenerate matrix H 0 it is enough to consider Ψ target = φ k β in definition (2) of the function V , where φ k β is defined in the above lemma. Since φ k β is continuously differentiable its derivative dφ k β /dβ is bounded on the interval [0, β * ] and the existence and uniqueness of β(Ψ) is guaranteed.

IV. NUMERICAL SIMULATIONS

We consider next the five-dimensional system (see [START_REF] Tersigni | On using shaped light pulses to control the selectivity of product formation in a chemical reaction: An application to a multiple level system[END_REF]) defined by (4). We use the Lyapunov control [START_REF] Ramakrishna | Control of a coupled two-spin system without hard pulses[END_REF] in order to reach the first eigenvector φ = (1, 0, 0, 0, 0, 0) of energy λ = 1, at the final time T . Note that µ 2 = 1. The function β is defined by (8) with Γ(x) = 0.75x. Even if Γ(x) doesn't satisfy the hypothesis Γ ′ ∞ < 1 8C of lemma 2.1, the convergence of β towards zero takes place (see Fig. 4, right image ). This happens because the condition is much stronger than the one required by the numerical simulations. Simulations of Figure 3 describe the evolution of the population of the trajectory Ψ = (Ψ 1 , Ψ 2 , ..., Ψ 5 ), for the initial state Ψ(t = 0) = (1, 1, 1, 1, 1)/ √ 5. We take c = 0.8.

Simulations of

V. CONCLUSIONS

In this paper we study implicit Lyapunov trajectory tracking procedures for closed quantum systems submitted to an external interaction, a laser field. The interaction between the system and the laser is described by a first order term ǫ(t)µ 1 and a second order term consisting in a polarizability term ǫ 2 (t)µ 2 . More precisely the hamiltonian of the Schrodinger equation that models the evolution is equal to H 0 + ǫµ 1 + µ 2 ǫ 2 , with H 0 the internal Hamiltonian.

The goal is to determine efficient controls for two types of controllable systems. For the first type there is not direct coupling between the target state and all the eigenvectors on the internal hamiltonian H 0 This corresponds in the bilinear setting H 0 +ǫµ 1 with the non controllability of the linearized system. The second one is characterized by non degenerate internal Hamiltonian. A description of the method, a convergence result together with numerical simulations are presented.

Fig. 1 .

 1 Fig. 1.Schematic view of the coupling realized by the matrices µ 1 and µ 2 defined by (4) between the first eigenvector φ of the internal Hamiltonian H 0 and all the other eigenvectors φ 2 , ..., φ 5 . The eigenvectors are represented by circles and next to each of them the corresponding energy level (H 0 ) ii , i = 1, ..., 5 are written. The coupling are represented by edges, each of them being labelled with the corresponding matrix that achieves the coupling.

Fig. 2 .

 2 Fig. 2. Schematic view of the slow convergence of φ β toward the target φ and a faster convergence of Ψ towards the mobile target φ β .

H 4 :

 4 λ j β = λ l β for j = l, j, l = 1, ..., N. (12) in the interval [0, β * ]. Since φ β is analytic the derivative dφ β dβ exists for every point β 0 ∈ [0, β * ], and is bounded on [0, β * ].

Fig. 4 ,Fig. 3 . 5 ;Fig. 4 .

 4354 Fig.3.The populations corresponding to system (4) with trajectoryΨ = (Ψ 1 , Ψ 2 , ..., Ψ 5 ); initial condition: Ψ(t = 0) = (1, 1, 1, 1, 1)√ 5; the feedback is defined by (27)( c = 0.8, k = 0.02)

Fig. 5 . 5 ;

 55 Fig. 5. Evolution of control ǫ ; initial condition: Ψ(t = 0) = (1, 1, 1, 1, 1) √ 5; system defined by (4) with feedback (27)( c = 0.8,k = 0.02).

This work was supported by "Agence Nationale de la Recherche" (ANR), Projet Blanc C-QUID number BLAN-3-139579, from PICS NSF-CNRS ,


Andreea Grigoriu

Stability analysis of discontinous quantum control systems with dipole and polarizability coupling

Keywords: Quantum systems, Stabilization, Lyapunov Control function, Tracking, Differential inclusions

Closed quantum systems under the influence of a laser field, whose interaction is modeled by a Schrödinger equation with a coupling control operator containing both a linear (dipole) and a quadratic (polarizability) term are analyzed. Discontinuous feedbacks, obtained by a Lyapunov trajectory tracking procedure, have been recently proposed to control these type of systems. The purpose of this paper is to study the asymptotic stability by considering the solutions in the Filippov sense. The analysis is developed by applying a variant of LaSalle invariance principle for differential inclusions. Numerical simulations are included to illustrate the efficiency of the discontinuous control.

Introduction

Control of quantum systems using laser fields has been subject to significant developments in the last two decades ( [START_REF] Assion | Control of chemical reactions by feedback-optimized phase-shaped femtosecond laser pulses[END_REF][START_REF] Brumer | Coherent chemistry: Controlling chemical reactions with lasers[END_REF][START_REF] Judson | Teaching lasers to control molecules[END_REF][START_REF] Levis | Selective bond dissociation and rearrangement with optimally tailored, strong-field laser pulses[END_REF][START_REF] Weinacht | Controlling the shape of a quantum wavefunction[END_REF] etc.). The increasing interest on this domain is motivated by the effects of the technique: we can create or break chemical bonds, each time with finesse far beyond the usual macroscopic means (temperature, pression, etc.).

Since the first successful laboratory experiments obtained at the beginning of the 90s [START_REF] Assion | Control of chemical reactions by feedback-optimized phase-shaped femtosecond laser pulses[END_REF][START_REF] Judson | Teaching lasers to control molecules[END_REF] many applications of this method have been developed: designing logical gates in future quantum computers, investigations of imaging by nuclear magnetic resonance -NMR, study of protein dynamics, molecular detection, molecular orientation and alignment, construction of ultra-short laser etc.. From the beginning, the complexity of chemical phenomena that arise during the interaction laser-quantum system has required the introduction of theoretical methods as an important step to experimental phase. This type of analysis can reveal the set of objectives that can be achieved, and the nature of the laser pulse that can be used. In this context, we consider the time dependent Schrödinger equation, that models the evolution of Email address: andreea.lachapelle@gmail.com (Andreea Grigoriu).

a quantum system:

i d dt Ψ(t) = H(t)Ψ(t) (1)
where H(t) is an Hermitian operator, called the Hamiltonian and Ψ a complex function called the wavefunction. When the system is controlled by selecting a convenient laser intensity ǫ(t), the interaction between the laser and the system is described by an operator µ 1 , also called dipole coupling [START_REF] Rabitz | Optimal control of selective vibrational excitation in harmonic linear chain molecules[END_REF]. Thus, we recover a billinear form of the Schrödinger equation, formally written:

i d dt Ψ(t) = (H 0 + ǫ(t)µ 1 )Ψ(t). (2)
In this case

H(t) = H 0 + ǫ(t)µ 1
, where H 0 is the internal Hamiltonian operator, that characterizes the system when the laser is shut down (ǫ(t) = 0). In the limit of small laser intensities the first order term ǫ(t)µ 1 may be enough to adequately describe the interaction, however, situations exist where the dipole coupling does not have enough influence on the system to reach the control goal; the goal may become accessible only after taking into account terms of higher order in the expansion of H(t), for example a polarizability term ǫ 2 (t)µ 2 (see e.g. [START_REF] Dion | Two-frequency IR laser orientation of polar molecules. Numerical simulations for HCN[END_REF][START_REF] Dion | Laser-induced alignment dynamics of HCN: Roles of the permanent dipole moment and the polarizability[END_REF] and related works).

In the following, we focus on the case where a second order them is added in the expansion of the Hamiltonian:

H(t) = H 0 + ǫ(t)µ 1 + ǫ 2 (t)µ 2 . (3)
For numerical reasons a finite dimensional setting is considered. The operators will be restrained to a linear space spanned by a N dimensional set D. This set can contain for example the first N eigenvalues of the infinite dimensional internal hamiltonian H 0 . For simplicity we conserve the same notations, i.e. we denote by H 0 , µ 1 and µ 2 , N × N Hermitian matrices with complex coefficients and by Ψ a N dimensional complex vector. One important problem is to determine efficient laser fields to control quantum systems whose Hamiltonian are defined by [START_REF] Brockett | Lie theory and control systems defined on spheres[END_REF]. For this purpose an analysis of the controllability has to be pursued, i.e. ask if any admissible quantum state can be attained with some admissible laser field. This can be studied via the general accessibility criteria [START_REF] Brockett | Lie theory and control systems defined on spheres[END_REF][START_REF] Sussmann | Controllability of nonlinear systems[END_REF] based on Lie brackets; more specific results can be found in [START_REF] Turinici | Beyond bilinear controllability: applications to quantum control[END_REF]. A detailed presentation has been made in [START_REF] Coron | Quantum control design by Lyapunov trajectory tracking for dipole and polarizability coupling[END_REF].

Even if positive results of controllability for systems, with Hamiltonian defined by (3), have been obtained, finding efficient numerical algorithms to determine the control field remains a very difficult task. A solution is to present the problem as a minimization of a cost functional, that describes the goal to be achieved, and eventually some other constraints. This approach leaded to procedures such as stochastic iterative approaches (e.g., genetic algorithms) [START_REF] Li | Optimal dynamic discrimination of similar molecules through quantum learning control[END_REF], iterative critical point methods (monotonic algorithms) [START_REF] Maday | New formulations of monotonically convergent quantum control algorithms[END_REF][START_REF] Tannor | Control of photochemical branching: Novel procedures for finding optimal pulses and global upper bounds[END_REF][START_REF] Zhu | Uniform rapidly convergent algorithm for quantum optimal control of objectives with a positive semi-definite Hessian matrix[END_REF], trajectory tracking or local control procedures ( [START_REF] Beauchard | Implicit Lyapunov control of finite dimensional Schrödinger equations[END_REF][START_REF] Chen | Competitive tracking of molecular objectives described by quantum mechanics[END_REF][START_REF] Ferrante | Control of quantum systems using model-based feedback strategies[END_REF][START_REF] Grivopoulos | Lyapunov-based control of quantum systems[END_REF][START_REF] Mirrahimi | Lyapunov control of bilinear Schrödinger equations[END_REF][START_REF] Rabitz | Quantum control design via adaptive tracking[END_REF][START_REF] Sugawara | General formulation of locally designed coherent control theory for quantum systems[END_REF] etc.). One advantage of this class of methods is that we obtain explicit control fields. Another one is that few propagations in time are required to approach the solution of the time-dependent Schrödinger equation (TDSE). This is an important aspect when larger systems are considered.

Lyapunov trajectory tracking techniques have been applied for systems with Hamiltonian (3) in order to determine the control ǫ. A first positive result has been obtained by adapting the analysis presented in [START_REF] Jurdjevic | Controllability and stability[END_REF][START_REF] Mirrahimi | Lyapunov control of bilinear Schrödinger equations[END_REF], that deals with bilinear quantum systems H 0 + ǫ(t)µ 1 .

The success of the feedback control depends on whether there exists (non-zero) direct coupling, through µ 1 , between the target state and all other eigenstates. When the same property holds for Hamiltonian H(t) = H 0 + ǫ(t)µ 1 +ǫ 2 (t)µ 2 the same type of feedback formulas hold. When some of the (direct) coupling is realized through µ 2 instead of µ 1 , the previous feedback formulas do not hold any more and two alternatives have been proposed (see [START_REF] Coron | Quantum control design by Lyapunov trajectory tracking for dipole and polarizability coupling[END_REF] for more details): discontinuous feedback and time varying feedback.

Only approximative asymptotic stability results have been proved for this last two situations. This paper focuses on the case of discontinuous feedback obtained for quantum systems with Hamiltonian defined by (3). The goal is to prove stability results, and especially asymptotic stability considering the solutions of the quantum system (1), with Hamiltonian H given by ( 3), in the Filippov sense.

The balance of the paper is as follows: in Section 2 we introduce the main notations, the Lyapunov tracking procedure followed by the construction of the discontinuous feedback. Then, we study the existence of solutions in the Filippov sense. In Section 3 we prove a first stability result followed by an asymptotic stability result. The last two sections are dedicated to numerical simulations and conclusions.

2 Lyapunov trajectory tracking

Lyapunov function

We consider equation ( 1), with Hamiltonian H(t) given by ( 3), that describes the evolution of a N -level quantum system submitted to an external action:

i d dt Ψ(t) = (H 0 + ǫ(t)µ 1 + ǫ 2 (t)µ 2 )Ψ(t). (4)
The wave function

Ψ = (Ψ j ) N j=1 is a vector in C N , veri- fying N j=1 |Ψ j | 2 = 1, i.e.
Ψ belongs to the unit sphere S N (0, 1) of C N . The function Ψ represents a complete physical description of the state of the quantum system at every instant t.

Recall that two wave functions Ψ 1 and Ψ 2 that differ by a phase θ(t) ∈ R, i.e. Ψ 1 = exp(iθ(t))Ψ 2 , describe the same physical state. To take into account the property we add a fictitious control ω (see also [START_REF] Mirrahimi | Lyapunov control of bilinear Schrödinger equations[END_REF]). Hence we will replace the evolution equation ( 4) by:

i d dt Ψ(t) = (H 0 + ǫ(t)µ 1 + ǫ 2 µ 2 + ω(t))Ψ(t), (5)
where ω ∈ R is a new control. We can choose it arbitrarily without changing the physical quantities attached to Ψ. We assume in the sequel that the state space is S N (0, 1) and the dynamics given by ( 5) admits two independent controls ǫ and ω.

In order to obtain an explicit formula for the laser field ǫ(t), we apply a Lyapunov trajectory tracking technique. The method consists in introducing a time varying function V (Ψ(t)):

V (Ψ(t)) = Ψ -φ|Ψ -φ = Ψ -φ 2 , (6)
with Ψ a smooth solution of (5) and φ an eigenvector of H 0 associated to the eigenvalue λ.

The function V is nonnegative for all t > 0 and all Ψ ∈ S N (0, 1) and vanishes when Ψ = φ. We search for feedback controls such that V is a Lyapunov function. To do that we compute formally the derivative of V along the trajectories of (5):

dV dt = 2ǫIm( µ 1 Ψ(t)|φ ) + 2ǫ 2 Im( µ 2 Ψ(t)|φ ) +2(ω + λ)Im( Ψ(t) | φ ), (7)
where Im denotes the imaginary part. For convenience we denote:

I 1 = Im( µ 1 Ψ(t)|φ ) and I 2 = Im( µ 2 Ψ(t)|φ ).
Then note that if, for example, one takes

ǫ(I 1 , I 2 ) = -kI 1 /(1 + kI 2 ) ω = -λ -cIm( Ψ(t)|φ), (8)
with k and c strictly positive parameters, one gets

dV /dt = -2k(I 1 /(1 + kI 2 )) 2 -2c(Im( Ψ(t) | φ )) 2 ≤ 0,
and thus V is nonincreasing. However, even if the feedback is chosen such that is V monotonically decreasing, this does not automatically imply that the minimum value will be reached. A convergence analysis is required.

Discontinuous feedback

The theoretical result (see Theorem 2.1) in [START_REF] Grigoriu | Lyapounov control of Schödinger equations: beyond the dipole approximation[END_REF] shows that tracking to φ works well when all eigenstates of H 0 , φ 2 , ....φ N , other than φ are coupled to φ by µ 1 , i.e. φ j , µ 1 φ = 0, j = 2, . . . , N . For the important case when some of the coupling are realized by µ 2 instead of µ 1 formulas (8) are ineffective. Discontinuous and time varying feedback have been proposed to stabilize the system (see [START_REF] Coron | Quantum control design by Lyapunov trajectory tracking for dipole and polarizability coupling[END_REF]).

The introduction of discontinuous feedback laws is motivated by the formula of the derivative of V with respect to time. We remark that dV /dt reads as the sum of 2(ω + λ)Im( Ψ(t) | φ) and a function U (ǫ):

dV dt = 2U (ǫ) + 2(ω + λ)Im( Ψ(t) | φ). (9)
Here U (ǫ) = ǫ 2 I 2 + ǫI 1 is a second order function of ǫ, with coefficients depending on I 1 and I 2 . Consequently, the condition dV /dt 0 depends on the sign of a second order function. Thus, the idea is to divide the space defined by I 1 and I 2 into disjoint regions. To each region we assign different formulas for the control ǫ(I 1 , I 2 ) such that U (ǫ) 0 in any point (I 1 , I 2 ).

To this goal we consider the regions (see Fig. 1):

A = {Ψ| I 2 (Ψ) < -|I 1 |}, B = {Ψ||I 2 (Ψ) > |I 1 |}, C = {Ψ| -|I 1 | I 2 (Ψ) |I 1 |} (10)
and we define the control as follows:

ǫ(I 1 (Ψ), I 2 (Ψ)) =        k 1 I 2 , in A 0, in B -k 2 I 1 /(1 + k 2 I 2 ), in C ω = -λ -cIm( Ψ(t) | φ ). ( 11
)
with k 1 , k 2 , c > 0. Remark 2.1 Under some restrictions for k 1 and k 2 that will be introduced later on, the condition U (ǫ) 0 is fulfilled on the region A and C. On region B we have U (ǫ) = 0. Remark 2.2 In order to guarantee 1+k 2 I 2 > 0 in equation [START_REF] Grigoriu | Lyapounov control of Schödinger equations: beyond the dipole approximation[END_REF], one notes that

|I 2 | ≤ | µ 2 Ψ(t)|φ | ≤ µ 2 ; therefore 1 + k 2 I 2 > 0 as soon as k 2 < 1
µ2 . From now on, unless otherwise specified, this condition will be supposed satisfied.

We replace the feedback (11) into equation ( 5) and we obtain a discontinuous right side equation:

i d dt Ψ(t) =                      H 0 + k 1 I 2 µ 1 + (k 1 I 2 ) 2 µ 2 - λ -cIm( Ψ(t)|φ ) Ψ, in A H 0 -λ -cIm( Ψ(t)|φ Ψ, in B H 0 -k 2 I1 1+k2I2 µ 1 + (k 2 I1 1+k2I2 ) 2 µ 2 -λ -cIm( Ψ(t)|φ ) Ψ, in C (12)
3

Existence of solutions in Filippov sense

Before discussing the stability of the system defined by [START_REF] Grivopoulos | Lyapunov-based control of quantum systems[END_REF] we need to study the existence of solutions. One idea is to consider solutions in the Filippov sense.

Definition 2.1 Let us consider the equation

ẋ(t) = f (x(t)) (13)
with piecewise discontinuous function f : D → R d , where d is the dimension of the space and D ⊆ R d a compact set. We recall that a solution in the Filippov sense of ( 13) is a locally absolutely continuous map such that:

ẋ ∈ F(x(t)) (14)
with

F := δ>0 µ(S)=0 conv(f (x + δB) \ N )) ( 15
)
where µ is the Lebesgue measure, S is a arbitrary set of measure zero, conv(A) is the smallest closed convex set containing A, B is the unit ball of R d and f is a discontinuous function.

In our case a solution in the Filippov sense, Ψ of ( 12) is a locally absolutely continuous map such that:

d dt Ψ ∈ F(Ψ(t)) (16)
with F defined by [START_REF] Levis | Selective bond dissociation and rearrangement with optimally tailored, strong-field laser pulses[END_REF] and f given by

f (Ψ) =                      -i H 0 + k 1 I 2 µ 1 + (k 1 I 2 ) 2 µ 2 - λ -cIm( Ψ(t)|φ ) Ψ, in A -i H 0 -λ -cIm( Ψ(t)|φ Ψ, in B -i H 0 -k 2 I1 1+k2I2 µ 1 + (k 2 I1 1+k2I2 ) 2 µ 2 -λ -cIm( Ψ(t)|φ ) Ψ, in C (17)
The way F (Ψ) is defined implies that it is a nonempty, bounded, closed, convex set. In the same time it is an upper semicontinuous function of Ψ. Thus, we are in the hypothesis of Theorem 1(page 70) in [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF]) and the existence of an absolutely continuous function solution of ( 16) is assured.

(Ψ(t) -Ψ) = 0. ( 19
)
We can give now a first stability result:

Theorem 3.1 Consider (12) with Ψ ∈ S N (0, 1) a solution in the Filippov sense and an eigenstate φ ∈ S N (0, 1) of H 0 associated to the eigenvalue λ. Take the constants k 1 > 1, k 2 < 1 µ2 and c > 0. The solution Ψ = φ of the inclusion ( 16) is stable.

Proof 3.1 Up to a shift on ω and H 0 , we can assume that λ = 0. Since the function V (Ψ) is C 1 with respect to Ψ, we can define the upper derivate by:

V * = dV dt * = sup y∈F (Ψ) (∇V • y). ( 20
)
For almost all t the derivative Ψ exists and satisfies the differential inclusion [START_REF] Li | Optimal dynamic discrimination of similar molecules through quantum learning control[END_REF]. For these t there exists:

V = d dt V (Ψ(t)) = ∇V • Ψ. ( 21
)
Theorem 1, page 153 in [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF] says that if V * 0 then φ is a stable point. In order to verify the fulfillment of this condition it is sufficient to make sure that dV /dt = ∂V • f 0 only on the domains of continuity of the function f defined by [START_REF] Maday | New formulations of monotonically convergent quantum control algorithms[END_REF]. In this domains we have F (Ψ) = f (Ψ).

On the discontinuity points of the function f the set F is defined as the closure of a convex set. This operation does not increase the upper boundary of the expression ∇V • f (see [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF] for more details).

On the interior of the region

B = {Ψ||I 2 (Ψ) > |I 1 |} the control ǫ(t) is zero, this implies dV dt = -2c(Im( Ψ(t)|φ )) 2 0. ( 22
)
We have the same property on the interior of the region

C = {Ψ| -|I 1 | I 2 (Ψ) |I 1 |}, since the control ǫ(t), is chosen such that: dV dt = 2 -k 2 I 1 1 + kI 2 I 1 + 2 k 2 2 I 2 1 (1 + kI 2 ) 2 I 2 -2c(Im( Ψ(t)|φ )) 2 = -2k 2 I 2 1 (1 + k 2 I 2 ) 2 -2c(Im( Ψ(t)|φ )) 2 0, ( 23
)
same conclusion on the interior of the region A = {Ψ| I 2 (Ψ) < -|I 1 |} since by hypothesis k 1 > 1. The condition V * 0 is fulfilled, thus we can apply here Theorem 1, page 153 in [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF], and the conclusion follows.

Asymptotic stability analysis

In the following we prove an asymptotic stability result for the system defined by [START_REF] Grivopoulos | Lyapunov-based control of quantum systems[END_REF] around the target φ. We apply a LaSalle type result for differential inclusions introduced in [START_REF] Ryan | An integral invariance principle for differential inclusions with applications in adaptive control[END_REF].

Theorem 3.2 Consider (12) with Ψ ∈ S N (0, 1) a solution in the Filippov sense and an eigenstate φ ∈ S N (0, 1) of H 0 associated to the eigenvalue λ. Take the feedback [START_REF] Grigoriu | Lyapounov control of Schödinger equations: beyond the dipole approximation[END_REF]

with k 1 > 1, k 2 < 1
µ2 and c > 0. Under the hypothesis:

(1) λ j = λ l for j = l, (2) for any j = 2, .., N : µ 1 φ j |φ = 0 or µ 2 φ j |φ = 0, where φ 1 , . . . , φ N is an orthogonal system of eigenvectors of H 0 corresponding to the eigenvalues (λ i ) i=1,...,N , the ω limit set of Ψ(t) reduces to ±φ.

Proof 3.2 Up to a shift in ω and H 0 , we may assume that λ = 0. Since we consider the solutions of the system (12) in the Filippov sense, the stability analysis will be made for the system defined by (16).

Theorem 2.11 in [START_REF] Ryan | An integral invariance principle for differential inclusions with applications in adaptive control[END_REF] says that the trajectories of the system (16) converge to the largest weekly invariant set contained in

E = {Ψ ∈ S N (0, 1)|0 ∈ V }, where V (Ψ) = {∇V (Ψ) • u, u ∈ F(Ψ)}, with V defined by (6).
Let us first compute the differential inclusion F defined by [START_REF] Levis | Selective bond dissociation and rearrangement with optimally tailored, strong-field laser pulses[END_REF], associated to the discontinous function f defined by [START_REF] Grivopoulos | Lyapunov-based control of quantum systems[END_REF].

If I 2 (Ψ) > |I 1 | the function f is continuous then: F = -i(H 0 -λI)Ψ. (24)
If -|I 1 | < I 2 (Ψ) < |I 1 | the function f is continuous thus: F = -i H 0 -k 2 I 1 1 + k 2 I 2 µ 1 + (k 2 I 1 1 + k 2 I 2 ) 2 µ 2 - λ -cIm( Ψ(t)|φ ) Ψ. ( 25
)
If I 2 (Ψ) < -|I 1 | the function f has the same property and:

F = -i H 0 + k 1 I 2 µ 1 + (k 1 I 2 ) 2 µ 2 -λ -cIm( Ψ(t)|φ ) Ψ. ( 26
)
On the contrary on the set {Ψ|I 2 (Ψ) = |I 1 (Ψ)|} the function f is discontinuous, hence

F (Ψ) = [b, c 1 ], if b c 1 [c 1 , b], if c 1 b (27)
Same property on the set {Ψ|I 2 (Ψ) = -|I 1 (Ψ)|}:

F (Ψ) = [a, c 1 ], if a c 1 [c 1 , a], if c 1 a, ( 28
)
where we have denoted:

a = -i H 0 + k 1 I 2 µ 1 + (k 1 I 2 ) 2 µ 2 -λ -cIm( Ψ(t)|φ ) Ψ, b = -i H 0 -λ -cIm( Ψ(t)|φ Ψ c 1 = -i H 0 -k 2 I 1 1 + k 2 I 2 µ 1 + (k 2 I 1 1 + k 2 I 2 ) 2 µ 2 - λ -cIm( Ψ(t)|φ ) Ψ. ( 29
)
If we take in consideration the computations made in the proof of Theorem 3.1 we have:

V =                  -k 1 I 2 (I 1 + k 1 I 2 2 ) -2c(Im( Ψ(t)|φ )) 2 , if I 2 < -|I 1 | -2k 2 I 2 1 (1+k2I2) 2 -2c(Im( Ψ(t)|φ )) 2 , if -|I 1 | < I 2 < |I 1 | -2c(Im( Ψ(t)|φ )) 2 , if I 2 > |I 1 |
(30) On the contrary on the discontinuity set {Ψ|I 2 (Ψ) = |I 1 (Ψ)|} we use relation [START_REF] Weinacht | Controlling the shape of a quantum wavefunction[END_REF] and we obtain:

V = [f, g], if f g [g, f ], if g f (31)
In the same way, on the set {Ψ|I 2 (Ψ) = -|I 1 (Ψ)|}, considering relation [START_REF] Zhu | Uniform rapidly convergent algorithm for quantum optimal control of objectives with a positive semi-definite Hessian matrix[END_REF] and we have:

V = [g, h], if g h [h, g], if h g, (32)
where:

f = -k 1 I 2 (I 1 + k 1 I 2 2 ) -2c(Im( Ψ(t)|φ )) 2 , g = -2c(Im( Ψ(t)|φ )) 2 h = -2k 2 I 2 1 (1 + k 2 I 2 ) 2 -2c(Im( Ψ(t)|φ )) 2 . ( 33
)
Since (I 1 + k 1 I 2 2 ) > 0, it follows that the limit set is characterized by:

I 1 = 0, I 2 = 0, Im( Ψ(t)|φ ) = 0, ( 34
)
and therefore ǫ = 0. This implies that the set E consists in fact of trajectories of the uncontrolled system:

i d dt Ψ = H 0 Ψ. ( 35
)
with solutions of the form:

Ψ = N j=1 b j e -iλj t φ j . (36)
We substitute (36) in (34) and we obtain:

Im( Ψ(t)|φ ) = Im(b 1 ) φ, φ + N j=2
Im(b j φ j , φ e -iλj t ) = 0. (37)

I 1 (Ψ) = Im(b 1 ) µ 1 φ, φ + j∈J1 Im(b j µ 1 φ j , φ e -iλj t ) = 0. ( 38
) I 2 (Ψ) = Im(b 1 ) µ 2 φ, φ + j∈J2 Im(b j µ 2 φ j , φ e -iλj t ) = 0. ( 39
)
Without loss of generality we take φ = φ 1 . From equation (34) and (37), together with φ j , φ = 0 for all j = 2, ..., N we obtain Im(b 1 ) = 0. Since along the trajectories Ψ in Ω(Ψ), I 1 (Ψ) ≡ 0, we have j∈J1 Im(b j µ 1 φ j , φ e -iλj t ) = j∈J1 B ′ j sin(λ j t + θ j ) = 0. The functions sin(λ j t + θ j ) are linearly independent as λ j are all different, hence the sum can only vanish if all coefficients B ′ j vanish. Observe now that B ′ j = 0, j ∈ J 1 if and only if b j = 0, j ∈ J 1 . Using Im(b 1 ) = 0 we have:

I 2 (Ψ) = j∈J2 Im(b j µ 2 φ j , φ e -iλj t ) = j∈J2 B j sin(λ j t + θ j ). ( 40
)
Since I 2 (Ψ) = 0 following the same arguments as above b j = 0 for j = 2, ..., N . Together with equality (36) this leaves only Ψ = b 1 e -iλt φ = b 1 φ (we assumed λ = 0). Since Im(b 1 ) = 0 the only case remained is Ψ = ±φ. This concludes the proof of Theorem 3.2.

Numerical simulations

We consider next the five-dimensional system (see [START_REF] Tersigni | On using shaped light pulses to control the selectivity of product formation in a chemical reaction: An application to a multiple level system[END_REF]) defined by:

H 0 =         1.0 0 0 0 0 0 1.2 0 0 0 0 0 1.3 0 0 0 0 0 1.4 0 0 0 0 0 2.15         , µ 1 =         0 0 1 1 1 0 0 1 1 1 1 1 0 0 0 1 1 0 0 0 1 1 0 0 0         , µ 2 =         0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0         . (41)
In Fig. 2 simulations describe the evolution of the population of the trajectory Ψ = (Ψ 1 , Ψ 2 , ..., Ψ 5 ), for the initial state Ψ(t = 0) = (0, 1/ √ 4, 1/ √ 4, 1/ √ 4, 1/ √ 4). We take k 1 = 1.1, k 2 = c = 0.8. We remark that the discontinuous laser field [START_REF] Grigoriu | Lyapounov control of Schödinger equations: beyond the dipole approximation[END_REF] is efficient to reach the first eigenstate φ = (1, 0, 0, 0, 0, 0) of energy λ = 1, at the final time T . Note that here µ 2 = 1. The Fig. 3 describes the evolution of the Lyapunov function defined by [START_REF] Coron | Quantum control design by Lyapunov trajectory tracking for dipole and polarizability coupling[END_REF].

Conclusions

We study in this paper the control of Schrödinger equation. The particularity of the problem is that the interaction between the system and the laser is not described just by a first order term ǫ(t)µ 1 , but also by a second order, polarizability term ǫ 2 (t)µ 2 . In a previous work discontinuous feedback with memory terms were introduced in order to exploit the polarizability coupling. The present paper studies this discontinuous case and focuses on obtaining an asymptotic stability result. Related numerical simulations are also presented. ; system defined by (41) with feedback [START_REF] Grigoriu | Lyapounov control of Schödinger equations: beyond the dipole approximation[END_REF]. We take k1 = 1.1, k2 = c = 0.8.

Fig. 1 .

 1 Fig. 1. Schematic view of the regions A, B,C. We consider arbitrary units for I1 and I2.

Fig. 2 . 4

 24 Fig. 2. The population of the system (41) with trajectory Ψ = (Ψ1, Ψ2, ..., Ψ5); initial condition:Ψ(t = 0) = (0, 1/ √ 4, 1/ √ 4, 1/ √ 4, 1/ √ 4); the feedback is defined by[START_REF] Grigoriu | Lyapounov control of Schödinger equations: beyond the dipole approximation[END_REF], with k1 = 1.1, k2 = c = 0.8.

Acknowledgements

This work was supported by "Agence Nationale de la Recherche" (ANR), Projet Blanc C-QUID number BLAN-3-139579 and Department of Energy. The author is also grateful to professor Jean-Michel Coron for useful discussions and comments.


Augustin Manecy
email: augustin.manecy@gipsa-lab.fr

Nicolas Marchand
email: nicolas.marchand@gipsa-lab.fr

member IEEE Stephane Viollet
email: stephane.viollet@univ-amu.fr

Bio-Inspired Hovering Control for a Robot Equipped with only a Decoupled Eye and a Rate Gyro

Keywords: Rotational Vestibulo-Ocular Reflex, tVOR Translational Vestibulo-Ocular Reflex, VFR Visual Fixation Reflex

This work provides an hovering control strategy for a sighted robot, the eye of which being decoupled from the body and controlled by means of a tiny rotative piezo motor. The main purpose of this paper is to show the effectiveness and the efficiency of this fundamental bio-inspired mechanical decoupling. Indeed, it exhibits several benefits:

• it enables to stabilize the robot's gaze on the basis of three bio-inspired oculomotor reflexes (ORs) : a visual fixation reflex (VFR), a translational and rotational vestibuloocular reflexes (tVOR and rVOR), • the eye can better, quickly and accurately compensate for sudden, untoward disturbances caused by the vagaries of the supporting head or body, • it yields a reference visual signal that can be used to unbias the rate gyro used to implement the VORs and to stabilize the hovering robot, • it increases the tracking accuracy with moving targets compared to without OR, This paper shows also that lateral disturbances are rejected 2 times faster with the decoupled eye robot, and roll perturbations induce a retinal error 20 times smaller. The occulomotor reflexes enables to cancel retinal error 6 times faster with 5 times lower retinal error picks. The conclusion of the paper is that decoupled eye must be considered as an efficient autonomous flight solution.

ACRONYMS FOV

Field Of View. rVOR

the quad-rotor. This system offers high frequency and good resolution, but is not mobile, and the robot is not totally autonomous. The same idea is exploited by [START_REF] Kim | Accurate modeling and robust hovering control for a quad-rotor vtol aircraft[END_REF] but using a CCD camera. An other strategy is to use markers disposed in the environment (one under the robot and one in front), as in [START_REF] Mori | Vision-based guidance control of a small-scale unmanned helicopter[END_REF] to reconstruct its attitude and position. In the same spirit, the use of embedded camera and different geometrical markers (five) in [START_REF] Zhang | Autonomous hovering of a vision/imu guided quadrotor[END_REF], provides an efficient estimate of both position and orientation of the robot. The same kind of task is too achieved in [START_REF] Karl | Low-cost visual tracking of a landing place and hovering flight control with a microcontroller[END_REF], with markers and but with the coarse visual sensor of the Wii remote. Finally, others technique using optical flow are proposed as in [START_REF] Rondon | Vision-based altitude, position and speed regulation of a quadrotor rotorcraft[END_REF] to estimate the altitude, position and speed above a specific geometric road [START_REF] Schafroth | Modeling, system identification and robust control of a coaxial micro helicopter[END_REF]. In all these approaches, the position follows from visual sensors when the attitude is either obtained thanks to the vision as in [START_REF] Gurdan | Energy-efficient autonomous four-rotor flying robot controlled at 1 khz[END_REF] or with embedded IMU's. In this last case, the gyrometer's bias is classically compensated using the accelerometer.

Contrary to the previously cited papers, in this paper, the robot is assumed to have an "eye" with free movements w.r.t. the "body" of the robot. That way, the visual sensor is able to rotate ant hence to change the gaze (eye's line of sight). This additional degree of freedom mimicks the mechanical decoupling between eye and body that is so characteristic of animals such as for instance the hoverfly. However, we will consider the robots eye as a sensitive and accurate visual position-sensing device (PSD), able to sense the position of an edge (or a bar) within only its small field of view FOV (here, FOV = 5 in opposition to FOV upper than 50 in the previous references). This sensors performance in the task consisting of locating an edge is a 40-fold improvement in resolution versus the interphotodiode angular resolution [43]. It can, therefore, be said to be endowed with hyperacuity [44]. For further details about the performance (i.e., accuracy and calibration) of this hyperacute visual PSD, see [40] and [43]. Gaze stabilization is a difficult task because the eye control system must compensate both quickly and accurately for any sudden, untoward disturbances caused by the vagaries of the supporting head or body. This finely adapted mechanism is way beyond what can be achieved in the field of present-day robotics. In addition with this controlled visual sensor, we assume that only a rate-gyro with an unknown bias embarked on the robot. This is more restrictive w.r.t. other contributions where an entire IMU is used.

The aim of this paper is to emphasize the important advantages of this added degree of freedom, even in the case when less MEMS sensors are embarked. For this, a robot with a decoupled eye is presented in the next section as well as its nonlinear dynamical model. In Section III, an observer is proposed to estimate position, speed, angles and rate-gyro bias. Based on this observer, the eye control and then the whole robot control are given. Section IV gives an extensive comparison of the behaviour of the robot with a decoupled eye and without. These results are then discussed in the last section.

II. SYSTEM OVERVIEW

A. Description of the twin-engine hovering robot

As shown in the CAD of figure 1, the robot will consist of a twin-engine aerial robot with three degrees of freedom (rotation around the horizontal axis: θ r , right and left translation: X and an eye rotation θ er with respect to the robot's body). Thanks to a mechanical decoupling between the eye and its mechanical support (the head), the eye can rotate freely in the robot's frame. In addition, the eye's orientation could be finely controlled by means of a fast and accurate piezo motor (PCB motor) featuring an extreme compactness (diameter of 20mm) and a very low mass (1 gramm). The mass of the overall robot will be about 100 grams and the robot will be completely autonomous in terms embedded computational resource and power supply. A wireless link based on a Bluetooth module will be implemented for the monitoring of several parameters and the sending of input reference signals to the robot from a personal computer. As shown in figure 3, the UAV roll angle θ r can be controlled by applying a differential rotational speed on the propellers. We assumed in this work that the robot flies at a constant altitude H. As the flying robot is under-actuated, its position X along the horizontal axis was controlled by adjusting its attitude around the roll axis. As a consequence, the robot remained still in space (hovering condition) when θ r = 0 • . In this work, the robot was considered to hover above a target placed onto the ground. In this paper, we deliberately placed our hovering control strategy in a bio-inspired minimalistic framework where the objective was to stabilize an underactuated hovering robot by means of only a drifting rate-gyro and an eye with a limited field-of-view. Figure 2 shows that the future realization of our robot will share many similarities with the fly :

cm

Propeller

• A rate-gyro : the fly has gyroscopic haltere organ measuring its bodys angular speed around the three degree of rotation (pitch, roll and yaw) (citer Hengstenber) whereas the robot is equipped with a classical MEMS rate-gyro.

• An optical position sensing device : the compound eye of the fly is able to locate a contrasting target placed in a small fontal part of the visual field (citer Collett and Boedecker) and the robot is equipped with an eye endowed with hyperacuity (citer Vodka).

• a neck : the fly has non-less than 23 pairs of muscles to control its heads orientation (citer Straussfeld). The robot has a decoupled eye actuated by means of a tiny position servomotor (rotative piezo motor).

• A proprioceptive sensor in the neck : the fly has prosternal organs which consists of a pair of mechanosensitive hair fields located in the neck region (citer Preuss Hengstenberg) and the robot is equipped with a contactless magnetoresistive sensor measuring the orientation of the eye relative to the head.

• A gaze stabilization : in the freely flying sandwasp, active gaze-stabilization mechanisms prevent the incoming visual information from being affected by disturbances, such as large bodys rotations around the roll axis (Citer Zeil). and a fly (a). These two dynamic under-actuated systems are able to measure their bodys rotational speed Ωr by means of a rate-gyro (for the robot) and halteres (for the fly) and to locate a contrasting target θt placed in a small part of their FOV. The fly has non-less than 23 pairs of muscles in its neck to stabilize its gaze θg whereas the robot can control the angular position of its eye θer by means of miniature rotative piezo motor.

To facilitate the following considerations, we introduced some notation.

• θ g : angular position of the gaze in the inertial frame θ g = θ er + θ r .

• ǫ r : retinal error defined by ǫ r = θ tθ g .

• X: position of the robot along the horizontal axis in the inertial frame.

• V x : speed of the robot along horizontal axis in the inertial frame.

• Y : position of the robot along the vertical axis in the inertial frame. Y was supposed to be constant (Y (t) = H ∀t). This assumption is not too strong because the control of height could be decoupled of the other control (citer ...) and realized independently.

• Ω r : roll rotational speed of the robot. Estimated values are denoted by an additional hat (e.g., θ). Reference values are denoted with a star (e.g., θ ⋆ ) and measured values are denoted with a bar (e.g., θ)

B. Non-linear model of the hovering robot

A classical non-linear dynamic model was assumed for the robot in the inertial frame:

                 Vx = -(F 1 + F 2 ) sin(θ r ) m -K vx V x Ẋ = V x Ωr = L (F 2 -F 1 ) I z θr = Ω r (1)
Where L is the distance between the robots center and a propeller, I z is the inertia momentum around the roll axis, F 1 and F 2 are respectively the thrust generated by the propeller 1 and 2 and K vx is the flapping coefficient which is supposed to be constant.

An inner speed feedback loop makes the robots angular speed (Ω r ) follow faithfully the rotational speed set points (Ω ⋆ r ) yielded by the attitude controller (see III-E). To control the rotational speed (Ω r ), we assumed that the propellers are controlled directly by adjusting the thrust value. Therefore the control input signal for a propeller is composed of a nominal thrust and a differential thrust. The nominal thrust (T ⋆ ) counteracts the gravity and the differential thrust (δ ⋆ ) generates the torque responsible of a roll rotation. The propeller control input signals are defined by

F ⋆ 1 (p) = T ⋆ (p) -δ ⋆ (p) and F ⋆ 2 (p) = T ⋆ (p) + δ ⋆ (p).
We consider the dynamic of the propellers as a first order system with a time constant equal to τ mot (see table ??).

Figure 3 shows the complete model of the robot including the model of the propellers and the model of the rate-gyro, the visual sensor and the angle sensor (used for measuring the eye in robot angle θ er ). For the dynamic simulations of the robot, computed under Matlab/Simulink environment, we used a set of parameters given in table I (see Appendix).

III. HOVERING BY GAZING

This work aims at describing a new control strategy for implementing a robust, accurate and fast hovering flight. In our approach where the gaze control takes a key role ( [START_REF] Kerhuel | A sighted aerial robot with fast gaze and heading stabilization[END_REF] and [START_REF] Kerhuel | Steering by gazing: An efficient biomimetic control strategy for visually guided micro aerial vehicles[END_REF]), the "hovering-by-gazing", which is presented here, amounts to maintaining the gaze automatically oriented toward a stationary (or moving) target. Contrary to many previous approaches (citer...) where the robots attitude is estimated from an inertial measurement unit (IMU) including rate gyro, accelerometers and magnetormeters, the robots roll θ r and its position X (with respect to the target) are here estimated only from the robots angular speed Ω r and the eye-in-robot orientation θ er . Our control strategy makes the robot minimize its retinal error signal thanks to different occulomotor reflexes and hover accurately above a target by estimating its attitude and position. Here, we will describe the eye control system and the hovering control system and explain how they interact.

A. The non-linear observer

As shown in figure 4, the state observer is the corner stone of our hovering control strategy. This observer delivers at its output non-less than 4 estimated states on the basis of only two measurements: the eye-in-robot orientation θer and the robots angular speed Ωr . Indeed, it provides an estimation of linear speed V x , position X, roll angle θ r and rate-gyro bias ∆ g :

                           Vx Ẋ θr ∆g       =     -T m sin( θr ) -K vx Vx + L 1 ( θer -θer ) Vx + L 2 ( θer -θer ) Ωr -∆g + L 3 ( θer -θer ) L 4 ( θer -θer )     θer = -arctan X Y 0 -θr (2)
We decided to implement a non linear observer because of the strongly non linear equations giving the evolution of the linear speed V x and the position X of the robot (see equation ( 2)). Non linearities in Vx provide better estimation during transient, and non linearities in θer steady state error in position estimation. Observer gain L = L 1 L 2 L 3 L 4 T was tuned using the classical LQG method using the linearised system around the origin.

B. An unbiased rate-gyro

To increase the efficiency of the observer, the rate-gyro bias is estimated, assuming that this bias is slowly varying. But, instead of using a classical unbias method based on the accelerometer measurements, we choose to estimate the rate-gyro bias only thanks to the eye system. And more particularly, thanks to the eye-in-robot angle which provides us a good means to avoid drift in position estimation. As a consequence, the rate-gyro bias is observable, only if the position estimation is accurate. This accuracy is achieved thanks to the non-linearity of the estimated output θe r. Where Ω r is the actual rotation speed, ∆ g is the rate gyro bias, and µ is an unknown noise (supposed centred).

Ωr = Ω r + ∆ g + µ (3) ∆g = 0 (4
Remark 1: In system (2), for the estimation of θ r , it is known from (3) that the actual rotation speed is Ω r = Ωr -∆ gµ. So the term Ωr -∆g is just the "model" term allowing to take care of rate gyro bias and L 3 ( θerθer ) is the innovation term. As supposed in (4), ∆g contains only an innovation term.

C. Eye controller

The control of the eyes orientation allows to keep the gaze locked onto the target placed onto the ground. In our bioinspired approach, the hovering control strategy consists of merging three complementary oculomotor reflexes :

• A rotational vestibulo ocular reflex, called rVOR, yielding at its out the signal θ ⋆ er θr which is simply equal to the opposite of the estimated roll angle θr • A translational vestibulo ocular reflex, called tVOR, which relies on the estimation of the robots linear position X, assuming the altitude of the robot is known. This reflex compensates for any translation disturbance applied to the robots body by yielding the output signal θ ⋆ er X which contributes to keep the eye locked onto the target.

• A visual fixation reflex where the visual feedback loop acts to cancel the retinal signal error ǫ r by controlling the eyes orientation θ er through the control input signal θ ⋆ er V F R (see figure 5). A PI controller acts to annul the retinal error ǫr . The ZSL "Zero-Setting Limiter" is here to prevent runaway of the eye if the target is lost. Gopt is the static gain of the visual sensor and d is a pure delay due to the image processing.

Visual Sensor

It is worth to note that the ZSL function shown in figure 5 clamps the retinal error back to zero whenever the latter tends to become higher (or lower) than a specified positive (or negative) level. This ZSL used in previous studies ( [START_REF] Viollet | A high speed gaze control system based on the vestibulo-ocular reflex[END_REF], [START_REF] Kerhuel | Steering by gazing: An efficient biomimetic control strategy for visually guided micro aerial vehicles[END_REF]) serves the same function as the limiter block used to model the inhibition of the smooth pursuit reflex whenever the position error goes beyond a fixed threshold (Citer Robinson and Stark and Young).

1) Visual Fixation Reflex control: This is the most important part of the control, which allows our robot to hover in windy conditions. This control a PI controller keeping the retinal error close to 0. The retinal error loop control could be summarized as it is shown by figure 5. To design the PI controller the pure delay due to image processing and the sampling is approximated by a first order Pade approximation.

The PI-controller gives an angle reference (θ ⋆ er V F R ) to compensate disturbance. In this part, to prevent modelisation error, the PI-controller is designed to obtain good stability margin (M φ = 65, M M ≥ 0.5).

2) Vestibulo-Ocular Reflexes control: These reflexes allows to counteract the two kind of robot's movements, to keep the target into the FOV. The rVOR is implemented here by a feedforward, using the estimation of the roll angle θ r to compensate rotation of the body. And the tVOR minimize the impact of lateral displacement onto the retinal error ǫ r .

3) Final eye control: To summarize, the reference angle θ ⋆ er (see figure 4) results from the contribution the three reflexes (rVOR, tVOR and VFR) as follows:

θ ⋆ er = θ ⋆ er θr rV OR + θ ⋆ er X tV OR + θ ⋆ er V F R V F R (5)
So, it is easy to see in figure 6, the contribution of the different reflexes. And in figure 7 and ?? the advantage provided by the VOR (implemented by the feedforward: θ ⋆ er θr and θ ⋆ er X ). To conclude, when an unknown disturbance occurs, the VFR reacts first, followed by the VORs to compensate movements engaged to reject the wind. And it is the contrary when a wished movement is accomplished: the VORs react first, and the VFR compensates model errors. The retinal error 5 times lower with the ORs on: less than 7 while it reaches more than 35 when ORs are off. And the retinal error is vanished 6 times faster, in only 0.25 seconds with ORs activated and more than 1.5 seconds without ORs.

D. Rotation speed controller

The speed rotation controller enables to track the speed rotation set points (Ω ⋆ r ) given by the attitude controller (see III-E). To control the speed rotation it is assumed the propellers are commanded directly in thrust (see II-B).

A LQG controller-observer is chosen for tracking rather than a simple PID controller because the measure Ωr also contains the rate-gyro dynamics. This choice allows to take care of it and increase the tracking quality. So it gives the following system: Final command: Different reflexes for different situations. When a voluntary displacement appears (here the robot moves 1 meter from the target at time t = 5s), the VORs react first, and generates a rotation to compensate robot movements (the movements are estimated thanks to the observer and the control Ω ⋆ r ). Then, during the displacement, the VFR appears to compensate model errors of the observer, and keep the target in FOV. Contrarily, when a lateral perturbation looms (at time t = 7s), the optic controller is the first to react, and generates a rotation to keep the target into the FOV. And only after few milliseconds, the VORs appear (when the robot begins to move to counteract the displacement generated by the disturbance).

Ωr Ωr = -2L Izτmot 0 1 τgyro -1 τgyro Ω r Ωr + 2L Izτmot 0 δ ⋆ (6)
3 $ er 3 $ er X (feedforward) 3 $ er 3r (feedforward) 3 $ er V F R (feedback)
Where Ω r is the actual roll rotation speed, Ωr is the rotation speed measured by the rate-gyro. Then an integral effect is added to reject static error in the Ω r pursuit. Finally, the closed loop time response for Ω r is less than 20ms, and the command noise is less than 2%.

E. Position and attitude controller

The position and attitude controller is implemented by a pseudo LQG controller with the non linear observer. The dynamic of the speed rotation loop, previously described (see III-D) is assumed to be very fast, and could be neglected. So the control to track the position reference (on X axis) and stabilize the attitude is the speed rotation of the robot (Ω r ). This control signals is then obtained by a simple LQG estimated-state-feedback which the estimated states (V x , X and θ r ) are provided by the non-linear observer. To obtain the LQG state feedback the system (1) is linearised around the origin with the equilibrium input Ω req = 0. To cancel the steady state error, an integral effect is added on the state X.

F. Control strategy for non decoupled eye system

The control strategy applied to the non decoupled eye robot is exactly the same than for the decoupled eye robot. The difference being that the eye is fixed (see remark below), so the rVOR, tVOR and VFR are disabled. It is assumed that the measure θer is directly replaced by ZSL(ǫ r ) (retinal 7. When the ORs are off, the retinal error reaches more than 35, and the eye have a lot of difficulties to lock the target because the robot movements are too quick. Luckily, when the robot come back to a roll angle close to 0, the eye can locked the target. Note the role of the ZSL which is here very important to avoid run away of the eye. Contrarily, when the ORs are ON, the retinal error keeps less than 7, and the target comes back to the FOV in less than 0.2s. The ORs allow to recover efficiency the target, even if asked movements are fast.

error and the anti run away), where:

ZSL : R → R, ǫ r → ǫ r if ǫ r F OV 0 if ǫ r > F OV (7)
So the observer is fed by ZSL(ǫ r ) which corresponds to the angle θ t + θ r , while the target is in the FOV. Remark 2: As for the decoupled eye system, the different angles for the non decoupled eye system are linked by equation ?? with θ er = 0.

IV. SIMULATIONS & COMPARISON OF BOTH BEHAVIOUR

A. Simulations conditions

For simulations, sensors have the characterization described in appendix by tabular II. Conditions have chosen to be as close as possible to real world. So noises and default are implemented. All parameters used for the control are voluntary worse initialised and estimated, it is mean that all parameters differs by 10% compared to the model parameters.

B. Comparison of performances for the two systems

In figure 8 and9, is applied to both systems, a same lateral disturbance and roll disturbance. So, a lateral disturbance, involving a half FOV displacement, is applied at time t = 2.5s and a roll impulse wind is applied at t = 13s.

1) Translation disturbances: Even if the same lateral disturbances is applied to both systems, it is easy to remark that the D-EYE system present a faster rejection in both cases. As it is shown in figure [START_REF] Kerhuel | A sighted aerial robot with fast gaze and heading stabilization[END_REF], the F-EYE system can not adopt big roll angle because of it necessarily implies a retinal error increase. Indeed it is visible in figure 9 that the limit of FOV is reached, for F-EYE system, even if the roll angle is small. In conclusion, the decoupled eye allows to be more responsive to lateral disturbances without taking the risk to loose the target from FOV.

2) Rotation disturbances: In the case of the roll disturbance, the advantage of the decoupled eye system is more impressive. While the disturbance implies a big roll angle movement and a big retinal error, for no-decoupled eye system. The fast dynamic of the eye, provides to the decoupled eye system, the ability to counteract immediately the roll disturbance (retinal error and roll movements are very small.) To conclude, the decoupled eye provides a rejection, for lateral disturbances, 1.5 times faster than a classical visual robot. And a roll disturbance has induces practically any lateral displacement and less than 3 roll angle for the biomimetic robot, while the classical robot undergoes a 3 centimetres displacement and more than 15 roll angle.

V. DISCUSSION OF A DECOUPLED EYE ADVANTAGES

Here are balanced advantages and disadvantages to have a decoupled eye rather than no decoupled eye with a large FOV. The first disadvantage is naturally the additional electro-mechanical complexity introduced by the new degree of freedom. As the same, we supposed that the image treatment is not too long. Here we summarize some advantages to have a decoupled eye rather than a fixed eye with large FOV:

• Better peripheral acuity: Technically, it is difficult to obtain a sensor with a large FOV and a good peripheral vision. The peripheral image is often blurred, and worse mainly if the robot is moving. So the measured error between the target direction and the gaze ǫ r would be less accurate in the periphery. A thin FOV, yields a better measurement for the retinal error.

• Reduce rotational optic flow: In the same, a decoupled eye allows to significantly reduce the rotational optic flow. And, as explained in [START_REF] Hateren | Blowfly flight and optic flow. ii. head movements during flight[END_REF], the rotational optical flow gives any informations of attitude contrary to the translatory optic flow, dominating with a decoupled eye. Then, the decoupled eye provides better informations for attitude estimation.

• Faster disturbance rejection: The decoupled eye system, enables more stiff roll angle to counteract sudden, untoward disturbances, without loosing the target from FOV. Indeed, for a same lateral disturbance, the decoupled eye system rejects 10 times faster the retinal error and 1.5 times faster the lateral offset.

• Lower calculation cost: A small FOV implies reduced image treatment.

VI. CONCLUSION

In this paper, a new way to stabilize an aerial robot and reject wind disturbances during hovering, thanks to a decoupled eye system, is developed. This decoupled eye allows to be more efficient to reject sudden untoward disturbances thanks to occulomotor reflexes (rVOR, tVOR and VFR) possible because of the low eye's inertia. These reflexes are inspired by the insects behaviour, and provide to the robot the ability to counteract more efficiency sudden untoward disturbances. In spite of have a small field of view renders more difficult the tracking of a target, the new degree of freedom provided by the neck compensate largely this handicap. Indeed, the fast dynamics of the eye, allowing in a first hand, to keep the target in FOV thanks to ORs. And in the other hand, to provide (thanks to the proprioceptive angle sensor) a measurement of the eye-in-robot which enables to unbias the rate-gyro while estimate the current attitude and position. In this way, a new unbias approach for rate gyro is proposed, using the visual loop to estimate rate-gyro drift.

The success of the presented method is possible because the eye dynamics are very fast in comparison with the robots body dynamics. The proprioceptive sensor, providing the eye-in-robot measurement θer , is fundamental, because it enables to estimates the attitude and position of the robot in spite of the decoupling between the eye and the body.

In a future work, this methodology will be implemented in an autonomous hovering robot, equipped with a very light and compact piezo motor to actuate the eye. And an expansion to a total autonomous aerial quadri-rotor with 6 degrees of freedom will be considered. In order to control the robot only with the measure of its gaze orientation, and the rate-gyro's datas, while rejecting sharply wind disturbances.

VII. APPENDIX:SENSORS AND ACTUATORS

The robot is equipped with two motors and propellers allowing it to control roll angle and consequently its position on X axis. The height of the robot is fixed, it is mean that the robot evolves along an horizontal rail without friction. And the eye system is actuated by a servo motor allowing it to control the gaze.

A summary of the different sensors characteristics is given in the table II.

Fig. 1 .

 1 Fig. 1.CAD view of the 100-grams micro-air vehicle in which the miniature rotary piezomotor (PCB Motor) controls the orientation of the eye relative to the body 'angle θer). The orientation of the robot around the roll axis (angle θr) is controlled by applying a differential rotation speed onto the propellers. The robot itself is mounted at the tip of a rotating arm allowing an horizontal translation of the robot.

Fig. 2 .

 2 Fig.2. Similarities between the hovering robot with a decoupled eye (b) and a fly (a). These two dynamic under-actuated systems are able to measure their bodys rotational speed Ωr by means of a rate-gyro (for the robot) and halteres (for the fly) and to locate a contrasting target θt placed in a small part of their FOV. The fly has non-less than 23 pairs of muscles in its neck to stabilize its gaze θg whereas the robot can control the angular position of its eye θer by means of miniature rotative piezo motor.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Block diagramm of the complete system. Green variables represent control signal, blues correspond to the different physical variables and red to different measures.

Fig. 5 .

 5 Fig.5. Block diagram of the Visual Feedback Loop used to implement the visual fixation reflex (VFR). A PI controller acts to annul the retinal error ǫr . The ZSL "Zero-Setting Limiter" is here to prevent runaway of the eye if the target is lost. Gopt is the static gain of the visual sensor and d is a pure delay due to the image processing.

Fig.

  Fig.Different reflexes for different situations. When a voluntary displacement appears (here the robot moves 1 meter from the target at time t = 5s), the VORs react first, and generates a rotation to compensate robot movements (the movements are estimated thanks to the observer and the control Ω ⋆ r ). Then, during the displacement, the VFR appears to compensate model errors of the observer, and keep the target in FOV. Contrarily, when a lateral perturbation looms (at time t = 7s), the optic controller is the first to react, and generates a rotation to keep the target into the FOV. And only after few milliseconds, the VORs appear (when the robot begins to move to counteract the displacement generated by the disturbance).

  Fig.7. When the ORs are off, the retinal error reaches more than 35, and the eye have a lot of difficulties to lock the target because the robot movements are too quick. Luckily, when the robot come back to a roll angle close to 0, the eye can locked the target. Note the role of the ZSL which is here very important to avoid run away of the eye. Contrarily, when the ORs are ON, the retinal error keeps less than 7, and the target comes back to the FOV in less than 0.2s. The ORs allow to recover efficiency the target, even if asked movements are fast.

Fig. 8 .

 8 Fig. 8. Evolution of X position for the same lateral and roll disturbances. At the top, it is shown the evolution of the position X, and the evolution of roll angle at the bottom.

Fig. 9 .

 9 Fig.9. Evolution of the retinal error for the two systems for the two disturbances. Retinal error for the decoupled eye system at the top and for the no-decoupled eye system at the bottom.

• θ t : angular position of the target in the inertial frame.• θ r : roll angle of the robot.• θ er : angle between the eye and the robot in the robot frame. This angle is mechanically constrained to a maximum angle: |θ er | < θ er M AX .

ACKNOWLEDGMENT

The authors would like to thank Robert Mahony, who have imagined with Stephane Viollet this system.