Layane Abboud
email: layane.abboud@supelec.fr

Andrea Cozza
email: andrea.cozza@supelec.fr

Lionel Pichon
email: lionel.pichon@lgep.supelec.fr

Performance Analysis of the Matched-Pulse-Based Fault Detection

Introduction

The problem of fault detection and location in wire networks has gained an increasing importance in the last few years [START_REF] Furse | Down to the wire[END_REF]. Wired networks are found in all modern systems, and are used to transmit different signals (control, alarm, etc.). That is why the issue of safe and reliable wiring systems is among the primary concerns of researchers and government agencies today [START_REF]Review of federal programs for wire system safety[END_REF].

There are several methods for wire testing, such as visual inspection, impedance testing [START_REF] Furse | Down to the wire[END_REF], and reflectometry methods which are widely used today to help detecting and locating wire faults. These methods send a predefined testing signal down the wire network to be examined. They include Time Domain Reflectometry (TDR), which uses a fast rise time step or pulsed signal as the testing signal, Frequency Domain Reflectometry [START_REF] Furse | Frequency domain reflectometry for on-board testing of aging aircraft wiring[END_REF] which uses multiple sinusoidal signals, sequence TDR [START_REF] Smith | Analysis of spread spectrum time domain reflectometry for wire fault location[END_REF] which uses pseudo noise, etc. Generally, hard faults (open and short circuits) are detectable through standard reflectometry, while soft faults (damaged insulation, etc.) are more critical to detect, especially when dealing with complex wire networks configurations.

In [START_REF] Abboud | Utilization of matched pulses to improve fault detection in wire networks[END_REF], we introduced the Matched-Pulse approach (MP), based on the properties of Time Reversal [START_REF] Fink | Time Reversal of Ultrasonic fields -part 1: basic principles[END_REF], as an improvement of the existing standard TDR. The MP method proposes to adapt the testing signal to the network under test, instead of using a predefined testing signal, as opposed to reflectometry methods. We have shown that this method results in a higher echo energy from the fault to be detected, when compared to TDR.

In this paper, we propose a sensitivity analysis of the MP approach, i.e., to evaluate the effect of the network topology (junctions, discontinuities, etc.) on the performance of this approach when compared to standard TDR. Through topological analysis, we will illustrate the role of different network elements governing the performances of both the TDR and MP. A mathematical analysis then follows in order to compare our approach to the TDR. The whole study will allow us to predict, for a given system configuration, to which extent the MP would be more beneficial compared to TDR, and state the configurations where our approach presents major advantages.

The paper is organized as follows: first we give the general assumptions we considered in our work, and then give a brief reminder about the MP approach, followed by a topological study of the influence of the network elements on the TDR performance. A mathematical analysis following this study allows to establish a comparison criterion of the MP and TDR, and finally simulation results permit a better validation of the discussed ideas.

The wire network

In this paper, we will consider uniform lossless transmission lines. For the sake of simplicity, we also consider resistive parallel loads and do not use reactive components, so that there is no dispersion. We will also consider that the testing signal (i.e., the injected signal) is the input to the system, and the reflected voltage wave is the output. In the presence of an eventual fault, we can either analyze the reflected signal directly or take the difference of the reflected signals, both in the presence and absence of the fault, in order for the echo corresponding to the fault to be more easily detected.

In this paper, provided that the systems we are studying are linear time invariant, we will be considering the difference of the reflected signals. Let the transfer function of the reference system without the fault be noted as H 0 ( f ), and the one corresponding to the system with the fault noted as H F ( f ). Consequently, analyzing the difference signal is equivalent to analyzing the output of a system whose transfer function

H( f ) is H( f ) = H F ( f ) -H 0 ( f ) (1)
We will refer to this system as the difference system.

The MP approach

Before proceeding with our topological study, we remind the basic concept of the MP approach. Instead of injecting a predefined testing signal into the network to be analyzed, as for existing reflectometry methods, the idea is to adapt the injected signal to the system, so that the energy of the echo from the fault is maximized, thus increasing its detection probability. Such a tailor-made signal, or MP, can be synthesized considering the matched filter approach used in signal processing [START_REF] Papoulis | Signal Analysis[END_REF].

According to this theory, a matched filter is obtained by correlating a known signal, or template, with an unknown signal to detect the presence of the template in the unknown signal. This is equivalent to convolving the unknown signal with a timereversed version of the template. The matched filter is the optimal linear filter for maximizing the signal to noise ratio (SNR) in the presence of additive white noise.

To resume, in the MP approach, the testing signal is synthesized as follows: we first inject the TDR testing signal in the network under test. The received echo is next time reversed; this signal is now the MP testing signal. We clearly notice that this signal changes according to the network's architecture, as opposed to the standard TDR.

In this paper, we propose an analysis of the performance of both TDR and MP, in order to better understand and characterize this approach. Through topological analysis and a mathematical study, we will be able to understand the reason why, in some cases, the MP presents major advantages over the TDR, and predict, for a given network architecture, the performance of the TDR and MP.

Topological study

After a brief description of the general assumptions considered in this paper, along with reminding the general idea of the MP approach, we now propose to analyze the impact of different network elements on the performances of both the TDR and MP approaches. But first let us introduce the topological representation we chose to better illustrate the discussed ideas.

Equivalent topological representation

We represented a physical model by means of a topological one, in order to highlight the elements influencing the wave propagation in the network under test, along with having a simple model to understand. We modeled the uniform lossless transmission lines with delay lines, which means that those lines are not physical objects, they only introduce time delays. We also represented the discontinuities in the system (i.e., junctions, branches, mismatched loads, etc.) with nodes, to illustrate the reflection and the transmission of the incident voltage waves. If we observe the example of Figure 1, we notice how the physical elements are modeled; furthermore, we represented, on the topological model, the reflection and transmission phenomena corresponding to the junction: when a voltage wave arrives at this position, a part is reflected back, and another is transmitted. Another example is illustrated in Figure 2, where a network is represented by means of its topological model. A final important point to discuss in this section is the variation of the amplitude of the peaks in the received TDR echo. In fact, let us consider the example of Figure 3, where the fault is modeled with two parallel lines, τ is the transmission coefficient of the junction, and Γ F and Γ F are respectively the reflection coefficients at the fault position and the source position. We propose to illustrate how the amplitude of a certain peak changes in terms of the number of interactions with the system discontinuities. Therefore we chose a simple trajectory of the wave, which is the one from the source, through the first junction, towards the fault and then all the way back. So, if we inject a pulse into the system, its amplitude will be modified by a factor of τ when passing through the junction. This voltage wave will then reflect on the fault and its amplitude will be once more modified by a factor of Γ F . The reflected wave will then follow the reverse path, and when arriving to the source position, the modification in its amplitude compared to the amplitude of the initial wave is τ 2 Γ F . This example illustrates the idea that whenever a wave passes through several junction, its amplitude will decrease each time due to the reflection phenomena. This fact will be useful in the section 5 when analyzing mathematically the performances of the TDR and MP.

Position of the network elements

In the standard TDR case, when using the difference system, the first peak in the TDR echo we obtain corresponds to the first interaction with the fault. This idea is illustrated in the example of Figure 4. The network to be analyzed contains a fault hidden from the source by a discontinuity, at 14 m from the source. We consider that the fault is a parallel resistive load which value is 600 Ω . The echo obtained without the presence of the fault is shown in Figure 5(a), in the blue line. We notice how the different peaks correspond to the interaction of the pulse with the system discontinuities; in particular the first peak we obtain corresponds to the position of the 140 Ω discontinuity (the time value on the x axis is the round-trip time to the discontinuity). Furthermore, when examining the echo in the presence of the fault, we notice how its presence introduced the first variation at about t = 0.9µs from the origin, corresponding to the round-trip time to the position of this fault. Of course the peak corresponding to the 140 Ω discontinuity remains unchanged, while the peaks corresponding to the discontinuities behind the fault (from the source position) are affected by the fault's presence, as we can clearly see when taking the difference signal shown in Figure 5(b). When analyzing this signal, we notice how the first variation corresponds to the fault's position, that is why we look for this peak in the detection process in TDR. Furthermore, we notice that this peak does not have the highest amplitude, as clearly shown in the figure.

Based on this example, we can notice that the network can be divided into two main parts, regarding the fault's position: the first one upstream of the fault (between the source and the fault) and the second downstream of the fault. When examining the reflected signal in the difference system, the peaks corresponding to the first part will eventually disappear because they are the same with and without the fault (as we already noticed in the previous example); as for the peaks corresponding to the discontinuities downstream of the fault, they will be modified by the presence of the fault, thus leaving several peaks in the difference signal. Another important point is that any change in the elements upstream of the fault would affect all the peaks, including the first peak (in the difference signal), whilst any change in the elements downstream of the fault will not affect this peak.

Detection gain

After discussing the influence of different network elements on the performances of both TDR and MP, we need to find a criterion allowing us to compare those two approaches. Therefore, we propose to define a gain G, which we called the detection gain. Let i(t) be the injected pulse in the time domain, and e(t) the received signal. We note that α i is the amplitude of the peak number i in the echo of the difference system, and t i its position. Given that we do not have any dispersion, and also the system we are studying is a cable network, we can write the impulse response of the system as follows

h(t) = ∑ i α i δ (t -t i ) (2)
where δ (t) is the Dirac pulse. Consequently, the reflected echo is

e(t) = i(t) * h(t) = ∑ i α i i(t -t i ) (3)
Knowing that, in the detection process, we are only interested in the peak corresponding to the fault, we propose to define a gain, which we called the detection gain. Let us first calculate the energy of the received signal. It can be written as follows E = e(t)

2 dt = ∑ i α i i(t -t i ) 2 dt (4)
when we calculate this expression, we find that

E = ∑ i α 2 i i 2 (t -t i )dt + ∑ i = j α i α * j i(t -t i )i(t -t j )dt (5)
If ∆t i j = t j -t i , let us denote by ϕ(∆t i j ) the autocorrelation function of the injected signal

ϕ(∆t i j ) = i(t -t i )i(t -t j )dt = i(s)i(s -t j + t i )ds (6)
In the TDR case, we assume that the injected pulse is normalized in terms of its energy. The energy of the echo can thus be written as

E = ∑ i α 2 i + ∑ i = j α i α * j ϕ(∆t i j ) (7)
We know that the coefficients α i are influenced by the discontinuities in the system (as discussed in section 5); their values decrease according to the number of interactions with the system discontinuities (the evolution of these values is a geometric series).

We also note that Equation ( 7) is composed of two terms: the first one is the average energy of the signal ∑ i α 2 i , and the second one is the mutual energy of the peaks of this signal ∑ i = j α i α * j ϕ(∆t i j ). When examining the expression of the mutual energy, we can notice that it contains two important parts: the first one α i α * j which depends on the number of interactions, and the second one is the correlation function ϕ(∆t i j ). If the temporal support of this function is narrow, then this expression can be neglected and the dominant part would be that of the average energy. In fact, we know that the first peaks we observe have the highest amplitude, and when injecting a pulse into the system, those peaks are separated one from another. We also know that the peaks appearing after those first peaks (in time) would have small amplitudes and the mutual energy between them can be neglected. So, we can reasonably assume that

∑ i = j α i α * j ϕ(∆t i j ) ∑ i α 2 i (8)
In this case, the energy of the first peak (the one we consider when detecting the fault) would be

E F TDR = α 2 1 (9)
In the MP case, the same procedure is applied, but here the injected signal is the time reversed version of the TDR echo. The MP echo is written as

e(t) = i(t) * h(t) * h(-t) (10)
Defining g(t) = h(t) * h(-t), we can write

g(t) = ∑ k β k δ (t -t k ) (11)
where β i is the amplitude of the peak number i. In fact, if we calculate g(t) in terms of α i , t j and t i we find

g(t) = ∑ i ∑ j α i α j δ (t + t j -t i ) (12)
Comparing Equations ( 11) and ( 12) we find that

β k = k= j-i ∑ i, j α i α j (13)
and also

t k = t j -t i (14)
Note that when we time-reverse the TDR echo to obtain the injected signal in the MP case, we are effectively doing a time shifting operation (i.e., g(t) = h(t) * h(t 1 -t)). The value of the peak at t = 0 corresponds to the fault; this peak has an amplitude equal to the average energy of the injected signal, that is the TDR echo in this case. Under the same assumption that the initial pulse (i.e., the injected pulse in the TDR) is normalized in energy, the energy of the MP peak is

E F MP = ∑ i α 2 i 2 (15)
Let E TDR and E MP be respectively the energies of the testing signals in the TDR and MP cases. In order to be able to compare the TDR and MP approaches, we need to have the same injected energy. That is why we propose to normalize according to the injected energy in both cases. The detection gain G can be thus defined as the ratio of the normalized energy of the MP peak to the normalized energy of the TDR peak

G = E F MP /E MP E F TDR /E TDR (16)
Based on the previous analysis, when the different peaks do not interfere one with another, G can be expressed as follows

G = ∑ ∞ i=1 α 2 i α 2 1 ( 17
)
We can clearly notice that G ≥ 1, and consequently the MP approach always presents an advantage over the TDR; this advantage is greater in the case where we have multiple peaks, resulting in a higher gain. In practice, we consider the first N peaks when calculating the gain, so as to have

G ≈ ∑ N i=1 α 2 i α 2 1 ( 18
)
In the next section we propose to illustrate the previously discussed ideas through simulation results. 1 Ω 80 open Fig. 7 The analyzed network in the second case, where the fault is embedded in the system. All lengths are in meters.

limited number of peaks, thus specifying an inferior limit to the gain (referred to as a threshold in the Table 1). In the first case (Figure 6), we chose a value of the fault equal to 600 Ω , corresponding to a soft fault. The inferior limit of G in this case is 1.55 (we considered N=5). When we calculate its exact value, we find G = 1.73. In the second case (Figure 7) when the fault is masked from the source by several discontinuities, we chose several values of the fault, illustrated in Table 1. In this configuration, the fault is separated from the source by a discontinuity at 5 m and a junction at 14 m. The numerical values obtained when calculating the inferior limit of the detection gain, along with its exact values are illustrated in Table 1. When examining those results, we notice that the more the fault is soft (i.e., the reflection coefficient at the source position is small), the more the MP becomes effective when compared to standard TDR. If we also compared the two studied configurations (Figures 6 and7), when the fault's value is 600 Ω , we notice that when the fault is embedded in the system, the gain's value is greater than the case where the fault is directly in front of the source. In fact, in the first case we have a greater number of peaks than in the second case; so as predicted by Equation 18, we have a greater value of the detection gain.

Conclusion

In this paper, we evaluated the impact of the network topology on the performance of the TDR and MP approaches. A topological study allowed us first to state what are the most influencing factors on the effectiveness of the TDR, then a mathematical analysis proved the advantage of the MP method over the standard TDR. The discussed ideas were finally verified through simulation results. This whole study allowed a better understanding of the factors influencing the TDR performance, thus enabling, for any configuration, to predict the effectiveness of the MP approach compared to standard TDR.

Fig. 1 Fig. 2 FFig. 3

 123 Fig. 1 An example of an equivalent topological representation.

Fig. 4

 4 Fig. 4 An example illustrating a wire network. The fault's value is 600 Ω .

Fig. 5 (

 5 Fig. 5 (a) The reflected TDR signals in the presence and absence of the fault. (b) The difference of the two reflected signals. The studied system is illustrated in Figure 4.

Table 1

 1 Numerical results obtained for the configuration of Figure7.

	Fault value (Ω ) Predicted gain (threshold) Calculated gain
	Short	2.33	2.6
	30	3.4	4
	600	7.7	9.22

6 Simulation results

Analyzed configurations

We consider the configurations illustrated in Figures 6 and7, where the fault is in front of the source in the first case, then masked from the source by several discontinuities in the second case. We chose those configurations because they represent two extreme cases; thus it seems interesting to compare the performances of the MP and TDR in both cases. Fig. 6 The analyzed network in the first case, where the fault is in front of the source. All lengths are in meters.

Numerical results

We simulated the voltage propagation in the configurations of Figures 6 and7 using transmission line theory as presented in [START_REF] Paul | Analysis of multiconductor transmission lines[END_REF]. The characteristic impedance of the lines is chosen to be 75 Ω (such as for some coaxial cables). We will verify the obtained values of the detection gain according to the general formula (Equation 16) and the simplified formula (Equation 18). In this last case, we will be examining a


S Bonnabel
email: silvere.bonnabel@ensmp.fr

M Mirrahimi
email: mazyar.mirrahimi@inria.fr

P Rouchon
email: pierre.rouchon@ensmp.fr

Qubit Hamiltonian identification: A symmetry-preserving observer-based approach ⋆

Keywords: Quantum systems, symmetries, asymptotic observers, estimation, nonlinear systems

We consider a pure two-state quantum system illuminated by two lasers. A photo-detector captures the fluorescence of the system. We build an invariant observer which yields a (local) estimation of the wave function (or density matrix) and the two key parameters (laser de-tuning and the atom-laser coupling strength), parameters of the hamiltonian. The design exploits the symmetries of the system and can be interpreted geometrically. The convergence proof is based on averaging arguments. Simulation with noise illustrates the robustness of the obtained estimation algorithm.

INTRODUCTION

In this paper, we consider a quantum system for which we would like to identify some parameters. This estimation problem has already been addressed with interesting theoretical contributions concerning well-posedness and identifiability (see, e.g., Claude Le Bris et al. [2007] and the reference herein) and the development of estimation algorithms robust with respect to the noises and the modeling errors. Here, we consider as a simple quantum system, a two level system described by a Schrödinger dynamics and a time-continuous population measurement. We propose an algorithm based on nonlinear asymptotic observer techniques preserving the symmetries [START_REF] Bonnabel | Invariant asymptotic observers[END_REF]) to estimate the system parameters. As far as we know, such recent techniques have not been applied to tackle this problem and this paper illustrates their potential interest. Since the observer design exploits the physical symmetries (invariance with respect to the frame-change), the filter equation admits a natural physical and geometrical interpretation that can be extended to higher dimensions.

The physical setup consists of an ensemble of identically prepared systems undergoing the same dynamics: for instance, dilute mono-atomic gases. They are very simple systems, in the sense their constituents (atoms), are perfectly identical and interact very weakly with each other. Atoms in such gases can be considered as perfect quantum systems, with a sequence of discrete energy states labeled |i , for i ∈ N, with increasing energies E i = hω a,i depending only on the atomic species considered. In order to measure the population of the ground state (i.e, the state of the lower energy), the system is illuminated with coherent light (a first laser) whose frequency is close to the transition frequency corresponding to the energy gap E j -E 0 = h(ω a, j -ω a,0 ) of a very unstable excited state (i.e, having a very short lifetime) to the ground state transition. It generates a transition |0 → | j for a part of the population illuminated from the ground state to the excited state, which spontaneously decays to the ground state emitting a photon.

⋆ This work was supported in part by the "Agence Nationale de la Recherche" (ANR), Projet Blanc CQUID number 06-3-13957.

The measurement of the number of photons emitted is then directly proportional to the population of the ground state |0 . Suppose there is also another laser whose frequency is close to the transition frequency of another excited state |e to the ground state |0 . The lifetime of the transition between the two latter states is supposed to be much longer than the previous one. To a first approximation the dynamics of the twostate system (ground state |0 and excited state |e having the longest lifetime) is described by a Schrödinger equation. We assume that some parameters are not well known: ∆ which is the difference between the second laser frequency and the atomic transition frequency, µ which is the atom-laser coupling strength and characterizes the Rabi frequency. The goal is to identify in real time ∆ and µ, measuring the ground state population thanks to the first laser (which generates a transition only for ground state population) and the photo-detector. The usual modeling of these open-quantum systems via Lindbald type terms in the density matrix dynamics (see, e.g., [START_REF] Haroche | Exploring the Quantum: Atoms, Cavities and Photons[END_REF]) is analyzed in [START_REF] Mirrahimi | Continuous measurement of a statistic quantum ensemble[END_REF], where singular perturbation techniques are applied to justify the adiabatic (quasi-static) approximations usually made by physicists and leading to a dynamical model described (up to higher order terms) by a Schrödinger equation despite a continuous population measurement.

In this paper we build an observer which estimates the wave function (in fact we rather use the density matrix language) and the two parameters ∆ and µ at the same time. We propose an extension of the observer described in [START_REF] Mirrahimi | Observer-based hamiltonian identification for quantum systems[END_REF], where only µ was estimated, and we show that symmetries (invariance by change of frame described by an element of SU(2)) play an important role. We use the standard notations in quantum physics literature (see e.g. [START_REF] Cohen-Tannoudji | Mécanique Quantique, volume I& II[END_REF], and [START_REF] Haroche | Exploring the Quantum: Atoms, Cavities and Photons[END_REF] for a more advanced lesson).

DYNAMICS AND INVARIANCES

Dynamics

The Schrödinger equation for the system writes:

ı d dt Ψ = ∆ 2 σ z + uµ 2 σ x Ψ, Ψ = Ψ 1 Ψ 2 ∈ C 2 (1)
where we let σ x = 0 1 1 0 , σ y = 0 -i i 0 , σ z = 1 0 0 -1 denote the Pauli matrices, ∆ is the difference between the atomic frequency transition (of ground state to excited state ω e0 = ω a,e -ω a,0 ) and the laser frequency ω, µ is the atomlaser coupling constant and u(t) ∈ R is the slowly varying amplitude of the laser. We have the useful foumlas σ 2 x = 1; σ x σ y = ıσ z (with circular permutation), the output is

y =< σ z Ψ, Ψ >= |Ψ 1 | 2 -|Ψ 2 | 2 = 2|Ψ 1 | 2 -1
as the measurement is the ground state population |Ψ 1 | 2 and the conservation of probability implies |Ψ 1 | 2 + |Ψ 2 | 2 = 1. For a justification (in the frame of singular perturbation theory) of this model (applying weak measurement) see [START_REF] Mirrahimi | Continuous measurement of a statistic quantum ensemble[END_REF].

It is convenient to write the dynamics with the density matrix: let ρ = ΨΨ † denote the complex matrix associated to the projector on the state Ψ. Supposing that the system is pure (meaning it is not entangled to its environment) implies both properties

Tr(ρ) = Ψ † 1 Ψ 1 +Ψ † 2 Ψ 2 = 1 and ρ 2 = ΨΨ † ΨΨ † = ρ. Thus rewriting (1) the system becomes ρ = -ı ∆ 2 σ z + uµ 2 σ x , ρ (2)
μ = 0 (3) ∆ = 0 (4) y = Tr(σ z ρ)
(5) where [,] is the commutator. We assume the laser amplitude to be slowly varying compared to the Rabi frequency |uµ| : | u| << |uµ||u| (the Rabi frequency is a characteristic of the absorption-emission cycle of photons for an illuminated atom). We assume, moreover, that that the frequencies of the laser and the frequencies of the atomic transition are close (i.e, smaller than the atomic width): |∆| << |u|µ.

Invariances (symmetries)

The system is invariant under a change of basis for the wave function Ψ → UΨ where U is any unit matrix of the Lie group SU(2). Indeed consider the transformation ϖ = UρU † , and

ζ x = Uσ x U † , ζ y = Uσ y U † , ζ z = Uσ z U † . With the new variables, the dynamics (2)-(5) writes d dt ϖ = -ı[ ∆ 2 ζ z + uµ 2 ζ x , ϖ] y = Tr(ζ z ϖ) ∆

Invariants, Invariant output, invariant vector fields

We are going to explain the form of the (symmetry-preserving) observer ( 6)-( 7)-( 8)-( 9). For simplicity's sake in this section we consider the reduced system (2)-( 5)

ρ = -ı ∆ 2 σ z + uµ 2 σ x , ρ y = Tr(σ z ρ)
since the group only acts on ρ and not on ∆ and µ. To find the form of observers which preserve the symmetries (i.e, are invariant under the action of SU( 2), see section 7), we follow the method of [START_REF] Mirrahimi | Observer-based hamiltonian identification for quantum systems[END_REF]. Even if the method does not fully apply to the two-state quantum system under study (all the results are local and the dimension of the group r must be strictly smaller than the dimension of the state space n), it gives guidelines which help us build a symmetry-preserving observer. We recall some basic definitions and results of this latter paper in the appendix (section 7). In particular, we need to build an invariant output error, an invariant frame, and scalar invariants of the group action (see eq (24) of section 7).

Invariant output error:

The output y = Tr(σ z ρ) is a scalar. It is invariant under the group action since for any

U ∈ SU(2) if we let ϖ = UρU † and ζ z = Uσ z U † we have y = Tr(ζ z ϖ) = Tr(Uσ z ρU † ) = Tr(σ z ρ).
Thus an output error is (see the definition 5 of section 7) : ŷy = Tr(σ z ( ρ -ρ)).

Invariant vector fields: The system (2) is invariant thus the second member of ( 2) is made of invariant vector fields (in the sense of the definition 3 of section 7). Inspiring from the form of the dynamics, let us take as invariant vector fields

-ı[σ x , •], -ı[σ y , •] and -ı[σ z , •].
They provide a global parameterization of the tangent bundle of the space of 2 × 2 projector matrices with trace 1. Notice that, they are functionally dependant.

Scalar invariants: A complete set of scalar invariants is a full rank function ( ρ, σ x , σ y , σ z ) → I( ρ, σ x , σ y , σ z ) ∈ R n+m-r which is invariant under the group action (where n is the dimension of the state space, m is the dimension on the inputs on which (ψ g ) g∈G act, and r is the dimension of the group G). Here n = 2, m = 3 and r = 3. Locally there are n + mr = 2 + 3 -3 = 2 scalar invariants (see [START_REF] Olver | Equivalence, Invariants and Symmetry[END_REF]). In fact we take the 3 invariants: Tr(σ x ρ), Tr(σ y ρ), Tr(σ z ρ). They are functionally dependent since Tr 2 (σ x ρ) + Tr 2 (σ y ρ) + Tr 2 (σ z ρ) = 1. (see section 3.2).

A SYMMETRY-PRESERVING OBSERVER

The observer

Consider the observer

d dt ρ = -ı ∆ 2 σ z + u μ 2 σ x , ρ (6)
-K ρ (Tr (σ z ρ) -y) (σ z ρ + ρσ z -2Tr(σ z ρ) ρ) (7) d dt μ = -uK µ Tr (σ y ρ) (Tr (σ z ρ) -y) (8) d dt ∆ = -uK ∆ Tr (σ x ρ) (Tr (σ z ρ) -y) (9)
where K ρ , K µ and K ∆ are positive scalars. It preserves the symmetries (see eq (24) of section 7) since ( 7) can be written (see section 3.2 for the proof)

d dt ρ = -ı ∆ 2 σ z + u μ 2 σ x , ρ + ıK ρ (Tr (σ z ρ) -y) (Tr (σ y ρ) [σ x , ρ] -Tr(σ x ρ) [σ y , ρ])
As we did for the true system, let us suppose that u is constant, and | ∆| ≈ ε|uµ|. To be able to apply the standard perturbation techniques for this type of physical system we choose the gains

K ρ = 4k ρ ε|u|µ, K µ = 2k µ ε 2 µ 2 , K ∆ = 2k ∆ ε 2 |u|µ 2 with ε > 0 small (ε ≪ 1), and k ρ , k µ , k ∆ ∼ 1.

Geometrical interpretation with the Bloch sphere

The Bloch sphere is a geometrical representation of the pure state space of a two-level quantum mechanical system. An important property is that any density matrix ρ can be written

ρ = 1 + Xσ x + Y σ y + Zσ z 2 , with ζ = X Y Z ∈ S 2
where 1 denotes the identity 2 × 2 matrix. We have Tr (σ x ρ) = X, Tr (σ y ρ) = Y and Tr (σ z ρ) = Z. Note that, the coordinate Z in this section is the output y. The commutation operation -ı[σ x , ρ] corresponds to the wedge product (1, 0, 0) T with (X,Y, Z) T (circular permutations allow to complete the correspondences). The dynamics of ζ is

d dt ζ = ( uµ 2 , 0, ∆ 2 ) T ∧ ζ
The dynamics does not depend on the choice of the orthonormal frame (invariance under the action of SO( 3)).

The correction term in (7) writes

-K ρ ( Ẑ -Z)[-X Ẑσ x - Ŷ Ẑσ y + (1 -Ẑ2 )σ z ]
and corresponds on the Bloch sphere to

-K ρ ( Ẑ -Z)( Ŷ , -X, 0) T ∧ ζ , where ζ = ( X, Ŷ , Ẑ) T .
It respects the symmetries (invariance by change of orthonormal frame) since it bounds ζ to remain on S 2 (meaning physically that the system is pure).

Choice of the invariant correction terms

This paragraph explains the geometrical motivations of the choice of the gains. More rigorous proofs are given in the next section. The correction term of ( 7) is such that it is a tangent vector to the sphere (so ζ remains in S 2 ) and its direction is such that it tends to make output error Ẑ -Z decrease. Indeed the dynamics of ζ is:

d dt ζ = 1 2   u μ 0 ∆   × ζ -K ρ ( Ẑ -Z) ζ × 0 0 1 × ζ But < ζ × 0 0 1 × ζ , 0 0 1 >=< ζ × 0 0 1 , ζ × 0 0 1 > ≥ 0
is always positive (where <, > denotes the usual scalar product). The correction term is a vector which is in fact always pointing towards "north" (i.e, Z = 1) if Ẑ < Z and towards "south" (i.e, Z = -1) if not.

K ρ is much bigger than K µ and K ∆ , meaning that the time scale of the convergence of Ẑ towards Z is faster. Let us thus suppose Ẑ -Z is close to 0. If u μ > uµ, then ζ tends to rotate faster around the X-axis than ζ , and since d dt Z = uµY , the difference Ẑ -Z tends to increase if Y > 0 and decrease if not. That explains (8) which also writes d dt μ = -uK µ Ŷ ( Ẑ -Z). The design (9) can be explained in the same way.

CONVERGENCE ISSUE

In this section we prove the existence of a Lyapunov function for the reduced system (after having performed an averaging of the periodic perturbations). To apply the rotating wave approximation one writes the system in the interaction frame, i.e, one makes the time dependant change of variables:

ρ = e -ı uµtσx 2 ξ e ı uµtσx 2 , ρ = e -ı uµtσx 2 ξ e ı uµtσx 2 .
We have thus

d dt ξ = [ı uµσ x 2 , ξ ] + e ı uµtσx 2 ( d dt ρ)e -ı uµtσx 2
and one can derive a similar formula for d dt ξ . Thus

d dt ξ = -ı ∆ 2 e ıuµtσ x σ z , ξ (10)
which shows the interest of the interaction frame. Note that σ z σ 2k x = σ 2k x σ z and σ z σ 2k+1 x = -σ 2k+1

x σ z for any k, and these relations are also true when σ z is replaced by σ y thus σ z e ı uµtσx 2 = e -ı uµtσx 2 σ z , σ y e ı uµtσx 2 = e -ı uµtσx 2 σ y . So the observer, in the interaction frame, reads:

d dt ξ = -ı ∆ 2 e ıuµtσ x σ z + u( μ -µ) 2 σ x , ξ -K ρ Tr e ıuµtσ x σ z ( ξ -ξ ) × e ıuµtσ x σ z ξ + ξ e ıuµtσ x σ z -2Tr e ıuµtσ x σ z ξ ξ d dt μ = -uK µ Tr e ıuµtσ x σ y ξ Tr e ıuµtσ x σ z ( ξ -ξ ) d dt ∆ = -uK ∆ Tr σ x ξ Tr e ıuµtσ x σ z ( ξ -ξ ) . ( 11
)
Secular first-order approximation Also called quasi-resonant, it consists in neglecting the terms rotating with high frequencies uµ and 2uµ, by averaging their influence on the evolution of ρ. The true dynamics consists of small oscillations around the solution of the averaged system. In order to compute the averaged system we use for a > 0

e ıaσ x σ z = σ z + aσ y - a 2 2 σ z -• • • = cos(a)σ z + sin(a)σ y (12)
and similarly e ıaσ x σ y = cos(a)σ y -sin(a)σ z .

(13) We write the cosine and sine as linear combinations of exponentials and we use interference formulas of the type:

Tr e ıuµtσ x σ z ( ξ -ξ ) e ıuµtσ x σ z ξ = 1 4 2Tr σ z ( ξ -ξ ) σ z ξ + 0 + 2Tr σ y ( ξ -ξ ) σ y ξ + 0
+ oscillating terms with mean 0.

We therefore get the following autonomous system for the averaged system/observer:

d dt ξ = 0 , d dt µ = 0 , d dt ∆ = 0 d dt ξ = -ı u( μ -µ) 2 σ x , ξ - K ρ 2 Tr σ y ( ξ -ξ ) σ y ξ + ξ σ y -2Tr σ y ξ ξ - K ρ 2 Tr σ z ( ξ -ξ ) σ z ξ + ξ σ z -2Tr σ z ξ ξ d dt μ = - uK µ 2 [Tr σ y ξ Tr σ z ( ξ -ξ ) -Tr σ z ξ Tr σ y ( ξ -ξ ) ] d dt ∆ = 0
Lemma 1. For all K ρ , K µ > 0, ξ and μ converge locally towards ξ and µ (c.f. [START_REF] Mirrahimi | Observer-based hamiltonian identification for quantum systems[END_REF]).

Proof. We consider the Lyapunov function

V = 1 2 Tr 2 (σ y ( ξ -ξ )) + 1 2 Tr 2 (σ z ( ξ -ξ )) + 1 K µ ( μ -µ) 2 Note that Tr d dt σ y ( ξ -ξ ) = -u( μ-µ) 2 Tr 2σ z ξ + • • • and Tr d dt σ z ( ξ -ξ ) = u( μ-µ) 2 Tr 2σ y ξ + • • • . When developing d dt V
the terms having u as a factor compensate each other and there only remain the terms with K ρ as a factor. Let δ ξ denote ξξ . We have

dV dt = -K ρ Tr 2 (σ y δ ξ )(1 -Tr 2 (σ y ξ )) -K ρ Tr 2 (σ z δ ξ )(1 -Tr 2 (σ z ξ )) + 2K ρ Tr (σ y δ ξ ) Tr (σ z δ ξ ) Tr σ y ξ Tr σ z ξ ≤ -K ρ Tr(σ y δ ξ )Tr(σ z ξ ) -Tr(σ z δ ξ )Tr(σ y ξ ) 2 ≤0.
Here we have used the fact that Tr 2 (σ x ξ ) + Tr 2 (σ y ξ ) + Tr 2 (σ z ξ ) = 1. See [START_REF] Mirrahimi | Observer-based hamiltonian identification for quantum systems[END_REF] for the end of the proof which is a standard application of the Lasalle principle.

Second order secular approximation Since the secular (nonoscillating) first-order terms vanish when computing d dt ξ and d dt ∆, one can not prove the convergence of the observer. One needs to compute the second-order approximation only for d dt ξ and d dt ∆. We apply the Kapitsa method described in e.g, Landau and Lifshitz [1982]. Suppose uµ = ν is large. ξ obeys a differential equation with a high frequency source term d dt ξ = f (ξ , νt). We proved that the mean of ξ over a period is constant. Integrating high frequency terms yields high frequency terms with same frequency and smaller amplitude: we seek a solution of the type

ξ = ζ + g 1 (ζ ,t) ν + g 2 (ζ ,t) ν 2 + • • •
where ζ is the mean of ξ over a period (recall

ζ = 0 + O(1/ν)).
Let us compute the first-order term, g 1 (ζ ,t), and neglect the second-order terms. On one hand, we have:

d dt ξ = 0 + g ′ 1 (ζ ,t) ν where g ′ 1 =
∂ 2 g 1 is the partial derivative of g 1 with respect to its second variable. But using (10) and neglecting third order terms:

d dt ξ = -ı ∆ 2 e iνtσ x σ z , ζ + g 1 (ζ ,t) ν . (14)
Gathering the two latter equations

g ′ 1 (ζ ,t) = -ν ı[ ∆ 2 e iνtσ x σ z , ζ ] + • • •
where we have not written the small terms compared to ν.

Integrating with respect to time t the last equation yields:

g 1 (ζ ,t) = ı[ ∆ 2 e iνtσ x σ y , ξ ] + 1 ν • • • .
Thus ( 14) can be re-written as

d dt ξ = -ı[ ∆ 2 e ıνtσ x σ z , ζ + g 1 (ζ ,t)/ν] + • • • = -ı[ ∆ 2 e ıνtσ x σ z , ζ ] + ∆ 2 4ν [e ıνtσ x σ z , [e iνtσ x σ y , ξ ]] + • • • .
Now let us compute the temporal mean (over a period) and only keep the secular terms. We apply ( 12) and ( 13). Applying the Jacobi identity

d dt ξ = -ı ∆ 2 2uµ [σ x , ξ ] + • • •
where we have not written the oscillating terms of 0 mean nor the terms of order O( 1 ν 2 ). Note that, here, we find the term ∆ 2 2uµ corresponding to the standard Bloch-Siegert shift.

We also need to develop ∆ up to second order terms (coefficients of 1 ν ) since it has constant mean. Its time derivative verifies

d dt ∆ = -uK ∆ Tr σ x ξ Tr e ıuµtσ x σ z ( ξ -ξ ) (15)
Note that, ξ and ξ are themselves solutions of differential equations depending on t via oscillating terms of frequency ν. they have an oscillating part of frequency ν and another one of frequency 2ν, since their time derivatives are sum and products of at most two oscillating terms of frequency ν. When considering up to the second order, only the terms of frequency ν in ξ and ξ can have a secular effect, as they multiply the ν-frequency term e ıuµtσ x σ z . Thus in this second order secular approximation the oscillating terms of frequency 2ν can be neglected.

But, the ν-frequency part of ξ is due to the integration of -ı ∆ 2 e ıuµtσ x σ z , ξ and so is for ξ . Thus

ξ = ζ + ı ∆ 2uµ e ıuµtσ x σ y , ξ + • • • ξ = ζ + ı∆ 2uµ e ıuµtσ x σ y , ξ + • • • (16)
where ζ (resp ζ ) is a solution of the averaged equation for ξ (resp ξ ) and the non-written terms are either 2ν-frequency terms or are of order O( 1 ν 2 ). Up to the second order, the secular terms of (15) can be calculated as the sum of two parts: 1. replacing ξξ by its ν-frequency part in order to compute the secular terms of Tr e ıuµtσ x σ z ( ξξ ) ; this will be then multiplied by -uK ∆ Tr σ x ξ ; 2. replacing ξ in Tr σ x ξ by its ν-frequency terms, leaving ξξ as it is and calculating the secular terms of Tr σ x ξ Tr e ıuµtσ x σ z ( ξξ ) .

On one hand, we have

Tr e ıνtσ x σ z ( ξ -ξ ) = 1 ν ∆Tr σ x ξ -∆Tr(σ x ξ ) + • • •
where we have not written the oscillating terms of 0 mean. Here, we have applied ( 16) and the following relation:

ı ∆ 2ν e ıνtσ x σ z [e ıνtσ x σ y , ξ ] = ı ∆ 2ν σ z [σ y , ξ ] -ıσ y [σ y , ξ ] + • • • .
On the other hand, for the ν-frequency part of Tr σ x ξ , using (16), we have:

Tr σ x [e ıνtσ x σ y , ξ ] = Tr cos νt σ x [σ y , ξ ] -sin νt σ x [σ z , ξ ]
= 2ıTr cos νt σ z ξ + sin νt σ y ξ .

Therefore, up to the second order, we have

d dt ∆ = - K ∆ µ Tr σ x ξ 2 ∆ -Tr σ x ξ Tr (σ x ξ ) ∆ + K ∆ ∆ 2µ [Tr σ y ξ Tr σ y ( ξ -ξ ) -Tr σ z ξ Tr σ z ( ξ -ξ ) ].
We have obtained the following, locally convergent, triangular system:

d dt ξ order 2 = -ı ∆ 2 2uµ [σ x , ξ ] (17) d dt ξ order 1 = -ı u( μ -µ) 2 σ x , ξ - K ρ 2 Tr σ y ( ξ -ξ ) σ y ξ + ξ σ y -2Tr σ y ξ ξ - K ρ 2 Tr σ z ( ξ -ξ ) σ z ξ + ξ σ z -2Tr σ z ξ ξ d dt μ order 1 = - uK µ 2 [Tr σ y ξ Tr σ z ( ξ -ξ ) -Tr σ z ξ Tr σ y ( ξ -ξ ) ] d dt ∆ order 2 = - K ∆ µ Tr σ x ξ 2 ∆ -Tr σ x ξ Tr (σ x ξ ) ∆ + K ∆ ∆ 2µ [Tr σ y ξ Tr σ y ( ξ -ξ ) -Tr σ z ξ Tr σ z ( ξ -ξ ) ].
Tuning the gains for the linearized system The lemma 1 proves (after using averaging arguments) that μ -µ and ξξ converge (locally) to 0 for any K ρ , K µ > 0. The last equation of (17) shows that once ξξ has converged to 0, ∆ -∆ converges to 0 for K ∆ > 0. Nevertheless in order to get a robust observer the tuning of the gains must respect the time scales. To choose appropriate gains, we consider the first order approximation of (17) around particular equilibrium points which are such that the linearized system writes simply around these points.

Consider the linearized error variables

ξ -ξ = 1 + Xσ x + Ỹ σ y + Zσ z 2 , μ = μ -µ, ∆ = ∆ -∆
where we have Tr (σ x ξ ) = X, Tr σ x ξ = X and the similar formulas for Y, Z, Ŷ , Ẑ. Consider the linearized system around ξ = 1-σ z 2 (i.e, Z = -1). Up to second order terms, we have

d dt ( ξ -ξ ) = (-u μ -K ρ Ỹ )σ y .
This can be written as

d dt X = 0 d dt Ỹ = -u μ -K ρ Ỹ d dt Z = 0 (18)
We also have Let us now write the linearized system around the other equi-

d dt μ = -uK µ (0 -Ỹ )/2 = uK µ Ỹ /2 (19)
librium point ρ = 1-σ x 2 (i.e, X = -1) d dt ( ξ -ξ ) = - K ρ 2 ( Ỹ σ y + Zσ z ), d dt μ = 0 and d dt ∆ = - K ∆ µ ∆ (20)
The interesting equations for the tuning are ( 18), ( 19) and (20).

To respect the time scales, set 0 < ε ≪ 1 and choose the gain K ρ ≈ ε|uµ|. In this case the observer filters the high frequencies and the average of ξξ tends to decrease. Choose a slower characteristic time of convergence for the parameters μ and ∆ with respect to the characteristic time of convergence of ξ

K ρ = 2k ρ ε|u|µ, K µ = 2ε 2 µ 2 , K ∆ = k ∆ ε 2 |u|µ 2
where k ρ , k ∆ > 0 are any scalar of order 1. The tuning of the gains is only made for the linearized system around particular equilibrium points. Nevertheless we think the nonlinear structure of the observer, based on the symmetries (and thus very close to the structure of the system), allows a good global behavior, as the simulations show.

SIMULATIONS

We take for the initial conditions:

CONCLUSION

One could have written directly the system on the Bloch sphere with the variables (X,Y, Z). But the averaging computations in the interaction frame are easier to write with Pauli matrices. Moreover, as it is proved in [START_REF] Mirrahimi | Observer-based hamiltonian identification for quantum systems[END_REF] when ∆ = 0, our method can be extended to more general systems for which the interpretation on the Bloch sphere is not possible anymore.

  and µ are unchanged by the transformation and we still have ∆ = μ = 0. ζ x , ζ y , ζ z respect the commutation relations of the Pauli matrices. Thus the system is invariant under the action of the transformation group SU(2) (see definition 2 of section 7).

Fig. 1 .

 1 Fig. 1. Measured output, output error, estimations of the parameters µ and ∆ without noise.

  ρ 0 = 1 + cos π 5 σ x + sin π 5 cos π 1.4 σ y + sin π 5 sin π 1.4 σ z 2 µ = 1 , ∆ = 1 5 , ρ0 = σ x ρ 0 σ xWe choose for the control u and the gains:u = 1, K ρ = 2ε|u|µ, K µ = 2ε 2 µ 2 et K ∆ = 2ε 2 |u|µ 2 with ε = 1 5 .The results are given by fig 1. In fig 2, The measured signals were added a white gaussian noise of amplitude 20% .

Fig. 2 .

 2 Fig. 2. Measured output, output error, estimations of the parameters µ and ∆ with noisy measurement. 7. APPENDIX : SYMMETRY-PRESERVING OBSERVERS

In this section we recall the basic definitions and results of [START_REF] Bonnabel | Invariant asymptotic observers[END_REF]. Consider the smooth system

where ϕ g and ψ g are local diffeomorphisms depending on a parameter g which is an element of a Lie group G such that

• ϕ e (ξ ) = ξ for all ξ ∈ X • ϕ g 2 φ g 1 (ξ ) = ϕ g 2 g 1 (ξ ) for all g 1 , g 2 ∈ G, ξ ∈ X . and ψ g verifies similar conditions. In section 2.3 we have G = SU(2) and for any

The property also reads d dt X = f (X,U), i.e., the system remains unchanged under the transformation (23). Definition 3. A vector field w on X is said to be G-invariant if the system d dt x = w(x) is invariant. This means w(ϕ g (x)) = Dϕ g (x) • w(x) for all g, x. Definition 4. An invariant frame (w 1 , ..., w n ) on X is a set of n linearly point-wise independent G-invariant vector fields, i.e (w 1 (x), ..., w n (x)) is a basis of the tangent space to X at x. Definition 5. The smooth map ( x, u, y) → E( x, u, y) ∈ R p is an invariant output error if

• the map y → E( x, u, y) is invertible for all x, u • E x, u, h( x, u) = 0 for all x, u • E ϕ g ( x), ψ g (u), h(ϕ g (x), ψ g (u)) = E( x, u, y) for all x, u, y Definition 6. (pre-observer). The system d dt x = F( x, u, y) is a pre-observer of ( 21)-( 22) if for all x, u F x, u, h(x, u) = f (x, u).

The definition does not deal with convergence; if moreover x(t) → x(t) as t → +∞ for every (close) initial conditions, the pre-observer is an (asymptotic) observer.

Definition 7. The pre-observer u, y).

The property also reads d dt X = F( X,U, h(X,U)), with X = ϕ g (x), U = ψ g (u). Assume that the output map is G-equivariant (see ?]), which is the case for the quantum mechanical system considered in this paper since the output y is a scalar invariant. Then a sufficient condition for the system d dt x = F( x, u, y) to be a G-invariant pre-observer for the G-invariant system d dt x = f (x, u) is:

where E is an invariant output error, ( x, u) → I( x, u) ∈ R n+m-r is a full-rank invariant function, the L i 's are smooth functions such that for all x, L i I( x, u), 0 = 0, and (w 1 , ..., w n ) is an invariant frame. This result is a consequence of the theorem 2 of [START_REF] Bonnabel | Invariant asymptotic observers[END_REF]. I is called a complete set of scalar invariants and verifies I(ϕ g ( x), ψ g (u)) = I( x, u) for any g ∈ G.