John Folkesson

Jacques Leederkerken

Rob Williams

Andrew Patrikalakis

John Leonard

A Feature Based Navigation System for an Autonomous Underwater Robot

We present a system for autonomous underwater navigation as implemented on a Nekton Ranger autonomous underwater vehicle, AUV. This is one of the rst implementations of a practical application for simultaneous localization and mapping on an AUV. Besides being an application of real-time SLAM, the implemtation demonstrates a novel data fusion solution where data from 7 sources are fused at dierent time scales in 5 separate estimators. By modularizing the data fusion problem in this way each estimator can be tuned separately to provide output useful to the end goal of localizing the AUV, on an a priori map. The Ranger AUV is equipped with a BlueView blazed array sonar which is used to detect features in the underwater environment. Underwater testing results are presented. The features in these tests are deployed radar reectors.

Introduction

Navigating underwater requires localization without GPS. Stationary features in the environment can be used to help the robot stay localized [START_REF] Dissanayake | solution to the simultaneous localization and map building (SLAM) problem[END_REF], [START_REF] Newman | \Towards constant-time slam on an autonomous underwater vehicle using synthetic aperture sonar[END_REF], [START_REF] Leonard | \Mapping partially observable features from multiple uncertain vantage points[END_REF], [START_REF] Leonard | \Incorporation of delayed decision making into stochastic mapping[END_REF]. Examples of sonar mapping can be seen in [START_REF] Ribas | using an imaging sonar for partially structured underwater environments[END_REF], [START_REF] Leonard | \Stochastic mapping using forward look sonar[END_REF] and [START_REF] Tena | \Concurrent mapping and localization using sidescan sonar[END_REF]. This is commonly known as simultaneous localization and mapping, SLAM . Some underwater robots have the advantage of rather expensive sensors such as Doppler velocity logs for measuring the movement of the robot. In this application the robot must be expendable which limits the choice of sensors. Instead of a single high quality sensor of the motion we have a number of weaker sensors that must be fused to estimate the motion of the robot.

The problem of fusing data from dissimilar sensors is challenging for a number of reasons [START_REF] Castellanos | \Multisensor fusion for simultaneous localization and map building[END_REF], [START_REF] Durrant-Whyte | \Toward a fully decentralized architecture for multisensor datafusion[END_REF], [START_REF] Moutarlier | \Stochastic multisensory data fusion for mobile robot location and environmental modelling[END_REF], [START_REF] Guivant | \Using absolute non-Gaussian non-white observations in Gaussian SLAM[END_REF], [START_REF] Sukkarieh | high integrity IMU/GPS navigation loop for autonomous land vehicle applications[END_REF]. First there is a temporal issue in that some sensors need to be incorporated with a simple calculation at a high update rate while other information might need intensive calculation but at a much lower update rate. Another issue is observablity. Generally each sensor provides only a partial observation of the total state. The problem is complicated by the fact that some sensors provide measurements in the relative frame of the robot while others are in the absolute earth frame. Measurement errors that are correlated over time due to bias or environmental disturbances can lead to over condent estimates if the errors are treated as independent.

A more practical issue is that of getting the system to work robustly. This requires that it be possible to examine each sensor's contribution to the nal estimate separately. Only then can one hope to understand the system well enough to tune the parameters.

Project Overview

The task of clearing a lane from sea to shore through a mine eld has traditionally been carried out by human divers. This method is associated with very high casualty rates but in many cases there are no alternatives. It is therefore very desirable to come up with methods that do not cost lives.

The system presented here is part of a larger project. The larger mission scenario starts with a target eld of underwater mines along a shoreline. The mine eld will be scouted rst by a REMUS AUV. The data collected from the REMUS will be used to create an a priori map of the mines and other strong sonar reecting objects along the shore. A particular target mine will be selected. The Nekton Ranger will then be sent on a mission to locate and attach to the target. As the Ranger is a less expensive AUV than REMUS it is meant to represent an expendable robot that could be set to detonate an explosive charge after attaching itself to the target.

The a priori map will not be perfect. It will have an absolute accuracy of about 5 meters. The relative accuracy between nearby features will be considerably better. The Ranger will be be able to get to within 70 meters of the target on the surface and to dive to the target depth. The task then becomes one of distinguishing between several nearby targets to nd the designated target. The AUV must then navigate to that target accurately.

The approach taken is to have the Ranger carry out simultaneous localization and mapping, SLAM, were the SLAM state is initialized with the a priori map including its uncertainty. By then matching the features detected by the Ranger to the a priori map the robot pose will become correlated to the a priori map. When this correlation is suciently accurate the target can be acquired with condence.

The a priori map will initially have relatively tight correlations between features that go like the inverse of the distance between the features. The absolute uncertainty of the features will be much looser allowing the map to shift as a whole to agree with where the Ranger nds it. Thus upon matching a feature seen by the Ranger with an a priori map feature both the Ranger pose and the map will shift to reconcile any dierences. The side scan sonar used by the REMUS to collect the data for the a priori map is dierent than the forward looking blazed array sonar on the Ranger. It is therefore likely that the maps will look dierent both in feature positions and numbers.

3 System Overview Figure 1 is a schematic over the navigation system. Its inputs, raw motion data and sonar detections, are written to as data arrives. These writes then trigger the lters downstream to update as data is pushed and pulled through the system. The data from the motion sensors is fed to sensor specic modules called input motion lters in gure 1. These model the individual sensors error characteristics. This sensor data is fused using ve separately tuned estimators. Each estimator has a specic purpose and frequency as summarized in table 1 The outputs of the input motion lters are pushed to the predict lter and to two separate Extended Kalman Filters (EKF). The predict lter uses a model of the forces and torques acting on the vehicle to estimate the motion over short periods of time. Each of the EKF lters has a specic function. The dead-reckoning lter provides the best continuous estimate of the robot pose. This estimate is used as the frame for accumulating data from the sonar. The reference EKF is identical to the dead-reckoning EKF aside from its inputs and outputs. GPS and long base line transponder data (LBL) are the additional inputs. These give measurements of the absolute x and y position of the robot. They also cause the estimated trajectory to have discontinuities in it. These 'jumps' are why the reference estimate is not used for feature tracking. The reference EKF also estimates more parameters than the deadreckoning EKF. The additional parameter estimates such as water current velocity can only be estimated using the GPS or LBL. The corrected parameters and their covariance are then passed to the predict lter to be used in prediction for both EKF lters.

The feature based navigation, FBN, lter consists of two separate estimators. The tracking lter and the global EKF. The tracking lter accumulates information from the sonar detections over short snippets of the robot path, (SLAM). The local map and trajectory from these snippets are then used to form composite measurements. These are then used to update the global EKF. Matching to the prior map is done on the global EKF. [START_REF] Leonard | \Incorporation of delayed decision making into stochastic mapping[END_REF] The Predict Filter A propeller on the aft end of the robot is mounted on a universal joint that can be pan/tilt rotated relative to the robot. An estimate of the pan/tilt angle is used to decompose this propeller thrust into robot pitch/yaw torques and a force forward. The movements of the robot through the water are resisted by forces and torques proportional to the linear and angular velocities. These forces then produce a convergent velocity for a given RPM. The parameters are tuned to give the empirical .75 m/s terminal velocity at 1,200 RPM, Fig. 2. The prediction model includes damping, actuation transfer, water velocity, buoyancy, and wave surge parameters, all modeled as Gaussian.

These various forces and torques are integrated at a rate of 100 Hz using Newton's equations to provide the EKF lters with estimates of the motion between the sensor measurements. The EKF lter can be congured to estimate and pass to the predict lter values for any of the parameters. The predict lter will then used those along with its own defaults for the others when making predictions. The right plot of Fig. 2 shows the improvement in the dead-reckoning from using the reference EKF estimated parameters.

5 Dead-reckoning EKF The dead-reckoning EKF estimates the robot pose, its velocity and angular actuation parameters, a total of 15 dof. As input it has depth, altitude velocity, and orientation from a 3-axis compass. The output is a smooth and accurate estimate of the motion between two sonar pings.

The dead service module does the translation from the EKF outputs at around 30 Hz to two pose estimates at the ping rate (1-2 Hz). The two pose estimates are a 6 dof incremental change estimate P i and a 4 dof earth frame measurement, P e . The two estimated poses are then able to account for the correlations between estimates taken at dierent times. The EKF provides snapshots of the uncertainty at the time of each sonar ping. We must break this uncertainty up into a part relative to the previous pose and a part relative to the earth. This is done by accumulating the predictions used in the EKF. Both poses' mean are the EKF values. The covariance of P i is taken as the accumulated prediction value. This then gives a system of matrix equations that can be solved for the covariance of P e . In this way we were able to have the tracking lter run at a much slower rate than the EKF lters. At the same time we had the exact pose and covariance estimate that the EKF provided.

Reference EKF

The reference EKF adds estimates of velocity damping, linear actuation, and water velocity parameters to the state vector estimated by the dead-reckoning EKF. A total of 27 dof. The purpose of the reference lter is quite dierent than that of the dead-reckoning EKF. The reference EKF is used to smooth the GPS and LBL measurements for use by the global EKF which operates at a much lower frequency. Only the xy parts are actually used. The global EKF will only use the reference xy estimate for an update if the covariance in xy has decreased (ie. when the Reference EKF has had a GPS update).

The reference EKF will only use a GPS measurement if its estimated covariance is less than that of the reference EKF. The reason for this seemingly odd update rule is that the GPS measurements have time correlated errors. Using consecutive GPS measurements and treating them as independent will cause the lter to become over condent very quickly. The rule softens this and essentially ensures that the covariance of the EKF will never get much lower than that of the single best GPS measurement.

The Tracking Filter

The tracking lter has the function of tracking features in the sonar data and using them to correct the motion of the robot over times of 10 to 100 seconds. These times correspond roughly to the time a feature is in the eld of view of the sonar. The output of this lter is a composite measurement that combines the dead-reckoning and feature tracking information over this interval. So a composite measurement gives the positions of any features tracked, the estimated pose change for the robot and the absolute (z,,, ) of the robot at the end of the interval.

The tracking is done in a small scale but highly detailed SLAM lter. This SLAM lter captures all the characteristics of the blazed array sonar. It estimates the entire path of the robot over the interval. Thus it adds 6 dimensions to its state for the robot pose on each new ping. The features are parameterized in the sensor frame of the rst pose they were seen from. In that way the feature will move when this pose is re-estimated. The parameterization in that frame is the range and bearing to the feature. This ensures that the uncertainty in these parameters will be closer to Gaussian than simple Cartesian coordinates would be.

The tracking estimator is a square root lter [START_REF] Dellaert | \Square root sam: Simultaneous location and mapping via square root information smoothing[END_REF]. This is a Gaussian estimator similar to the Extended Information Filter, EIF, which in turn is equivalent to the EKF. It has several important advantages. The most important is that as the square root of the information matrix is used the condition number of this matrix will be smaller. A smaller condition number means that the matrix equations that must be solved are numerically more well-posed. This will virtually always lead to an improvement over EKF or EIF for SLAM. Secondly, one can in our case do the updates in constant time. This is because the short snippets of path we use have no loops in them. Partial covariances can be computed in constant time as well. Estimating the full covariance would take a longer time but is actually never needed.

The details of this lter would take us beyond the scope of this paper and will appear in a later publication. Here we present the system as a whole and concentrate on how the parts work together in a particular application.

The Global EKF

The global EKF is updated with the composite measurements. This begins by augmenting the state with the pose and features of the composite measurement. The EKF state can then be updated using the composite measurement. The older previous pose can then be marginalized out. This implies simple deleting the rows and columns of the older pose. The new features can then be matched and merged with the older features including the a priori map features.

The matching of the features is done in three ways. First any features that were tracked in the tracking lter from one composite measurement to another are matched based on a key provided by the tracking lter. Second a threshold is placed on the mahanolobis distance between single features that are relatively well correlated to one another to determine matching features. Third multiple feature matches are search for. All of the pairs of possibly matched features are used a a basis over a hypothesis space. This space is search exhaustively if small enough or selectively if too large. The details of the selective search and match criteria are original and require an entire paper to describe. One part of our criteria is the joint compatibility test [START_REF] Neira | association in stocastic mapping using the joint compatibility test[END_REF]. This forms a distance between sets of points which is essentially the Mahanolobis distance of a pseudo measurement equating the point in each set. d = S 1 :

(

Where: S = C + J ;x C xx J T ;x :

(2)

Here C is the covariance of this pseudo measurement (a parameter) and should reect the nite size of the actual features, C xx is the covariance of the EKF state and J ;x is the Jacobian of the measurement. This is balanced against a threshold that is based on certain statistical measurements accumulated over the path of the robot. The problem of nding candidate combinations to apply this test to is where the diculty lies. The number of possibilities quickly becomes unmanageable as the number of features in the map increases. For this application we deal with a relatively small number of features and a very slow update rate for the maps. However we still need to limit this search at times. We start by forming matchdata ij objects between each pair of features, labeled i and j. This is a N 2 operation where N is the number of features. The match data consists of r ij , the Euclidean distance between the points,c ij , the covariance of rij, and a ij , the angle of the vector between them. rij = jx i x j j;

(3) c ij = J r;x C xx J T r;x (4) aij = arctan(y i y j ; x i x j )

(5)

matchdata ij = (r ij ; c ij ; a ij ) (6)
Remember that we are using a compass and will take full advantage of knowing the orientation of the robot and thus these angles in some cases. The cases that we feel reasonably condent about the angles between two features is when the two features are relatively tightly correlated. This means that they were either seen near each other in time as the robot moved or that they were matched to features that were in turn tightly coupled to one another.

We see that this concept of two features being tightly coupled is important. We dene tightly coupled as c ij < some threshold. Having this we can form N lists of matchdata ij objects where the i th list contains all the j that are tightly coupled to i. These are the local patches that we will try to nd multiple matches between.

Continuing to try and limit the number of possible matches we can note all pairs that are possibly the same point feature by either r ij or rij 2 =c ij being relatively small. These will be stored on a list of possible matches. We can then iterate through this list and search for multiple matches between the local patches associated with each of the possibly matched features. We check if there are other possible matched features between the local patches by checking the list.

Experiments

The experiments consisted of a series of underwater tests. A target eld was deployed consisting of 12 radar reectors anchored 1 meter o the sea oor. These were meant to represent a mine eld obstructing access to a part of the shoreline. The targets were deployed in three rows perpendicular to the shore. Two of the rows had ve reectors and a center row had 2 reectors. This target eld was then surveyed by the REMUS AUV (see Figure 1). This data was then used to form a hand made a priori map for the Ranger AUV to use, Figure 3 left.

The nal map made by the Ranger is seen to be very similar to the a priori map however it is shows some features displaced from the priori positions. The absolute xy of the map is of no signicance here the important point is that the robot pose is correlated to the a priori map. This despite signicant Fig. 3. On the left is the show the SLAM map as black ellipsoids and the a priori map as red crosses. This map made no attempt to match the two maps to one another. On the right is a map made with matching. The robot succeeded in matching on all 6 dives demonstrating the localization of the robot is working. dierences in the GPS location x. The robot here dove six separate times and was able to correctly match features to the map each time. In at least one instance, (the rst return to start location), the ranger FBN Filter gave better estimated positions then the Reference EKF had. [START_REF] Moutarlier | \Stochastic multisensory data fusion for mobile robot location and environmental modelling[END_REF]

Conclusion

The system was able to localize the robot on the a priori map in real time. This will allow an autonomous AUV to home in on a pre-designated target. The approach of fusing the dierent sensor data using ve separate estimators works well. This approach allows us to modularize the system which adds exibility and presumably some robustness. The various estimates give us a more precise view of the overall estimation process than a single integrated estimator could give. We can then identify the source of any errors more quickly. These error then give us a valuable guide to adjustments and renements of the estimators. The GPS was able to be used to improve the dead-reckoning estimate without introducing any discontinuities in the pose estimate. That thanks to a dual EKF implementation.

In the future we plan to work on the nal control of the robot to the designated target. We would also like to add some active perception to the system by having the robot rotate slightly both in pitch and yaw when targets are in view. This will help the convergence of the tracked feature's covariances.

Fig. 1 .

 1 Fig. 1. The Nekton Ranger is shown to the left. The separate modules of the navigation system are shown on the right.

Fig. 2 .

 2 Fig. 2. The tuning of the prediction for the yaw and pitch angle is done by comparing to the 3D compass as shown on the bottom left (yaw) and. upper left (pitch). The gure to the right compares the dead-reckoning and reference estimates of the robot trajectory. The three paths have been oset to separate them for comparison. The path to the right is the reference EKF output. It uses the GPS during the turn arounds on the surface. The middle is the dead-reckoning EKF, sharing the predict lter with the reference EKF. This allows it to benet from the parameter estimate made in the reference EKF. The left is the EKF, without any parameter sharing.

Table 1 .

 1 The ve estimators that fuse the various sensor data are summarized.

	Estimator	Output	Rate	Inputs
	Predict Filter Motion between	100 Hz Propeller Actuation signals, Various
		measurements		parameters
	Deadreck. EKF Smooth pose trajectory 30 Hz 3-Axis Compass, Depth, Altitude,
		for tracking	async	Predict output
	Reference EKF Corrected xy and	30 Hz 3-Axis Compass, Depth, Altitude,
		parameter estimates	async	GPS, LBL, Predict output
	Tracking Filter Composite	1 to 2	Sonar Detections,
	Squareroot SAM measurements	Hz	Dead-reckoning output
	Global EKF	Localization on an	.01 to .1 Composite Measurements and
		a priori map	Hz	Reference output


Paolo Forni
email: paolo.forni@inria.fr

Alain Sarlette
email: alain.sarlette@inria.fr.

Thibault Capelle
email: thibault.capelle@lkb.upmc.fr

Emmanuel Flurin
email: emmanuel.flurin@lkb.upmc.fr

Samuel Deléglise
email: samuel.deleglise@lkb.upmc.fr

Pierre Rouchon
email: pierre.rouchon@mines-paristech.fr.

Adiabatic elimination for multi-partite open quantum systems with non-trivial zero-order dynamics

We provide model reduction formulas for open quantum systems consisting of a target component which weakly interacts with a strongly dissipative environment. The time-scale separation between the uncoupled dynamics and the interaction allows to employ tools from center manifold theory and geometric singular perturbation theory to eliminate the variables associated to the environment (adiabatic elimination) with high-order accuracy. An important specificity is to preserve the quantum structure: reduced dynamics in (positive) Lindblad form and coordinate mappings in Kraus form. We provide formulas of the reduced dynamics. The main contributions of this paper are (i) to show how the decomposition of the environment into K components enables its efficient treatment, avoiding the quantum curse of dimension; and (ii) to extend the results to the case where the target component is subject to Hamiltonian evolution at the fast time-scale. We apply our theory to a microwave superconducting quantum resonator subject to material losses, and we show that our reduced-order model can explain the transmission spectrum observed in a recent pump probe experiment.

Introduction

The evolution of a quantum system interacting with an environment is rigorously described by a Schrödinger equation on the joint Hilbert space. However, the complexity of the environment hampers the study of the system as a whole and one often resorts to the Born-Markov approximation to obtain a Lindblad master equation [START_REF] Breuer | The Theory of Open Quantum Systems[END_REF] describing the target system alone, and the environment's effect summarized by dissipation or "decoherence" operators. Similarly, when a quantum system consists of several interacting components, e.g. a main computing subsystem coupled to an ancillary subsystem expressing a measurement device, one often seeks to analyze a dynamical equation for the main subsystem alone, approximately including the effect of the ancillary subsystem. In this perspective, model reduction methods come to aid to the physicists interested in gaining better physical insights, in running simplified numerical simulations, and in designing the dynamics of a target subsystem by smartly engineering its interaction with other subsystems, as in the case of reservoir engineering [START_REF] Poyatos | Quantum reservoir engineering with laser cooled trapped ions[END_REF].

A classical approach to model reduction for quantum systems makes use of the time-scale separation between a slow subsystem of interest and the fast auxiliary subsystems coupled to it, and eliminates the fast variables in a procedure denominated as adiabatic elimination. In closed quantum systems -where the evolution stays unitary under Hamiltonian dynamics -adiabatic elimination is performed by means of standard perturbation theory techniques [START_REF] Sakurai | Modern Quantum Mechanics[END_REF]. In contrast, the treatment of open quantum systems -including decoherence under Lindbladian dynamics -is more involved. In the literature, adiabatic elimination in the latter case has been addressed for specific examples separately: lambda systems up to second-order [START_REF] Brion | Adiabatic elimination in a lambda system[END_REF], a specific atom-optics example [START_REF] Atkins | Approximate master equations for atom optics[END_REF], systems where excited states decay toward n ground states [START_REF] Mirrahimi | Singular Perturbations and Lindblad-Kossakowski Differential Equations[END_REF][START_REF] Reiter | Effective operator formalism for open quantum systems[END_REF], systems with Gaussian dynamics and subject to continuous measurement [START_REF] Ond Řej Černotík | Adiabatic elimination of gaussian subsystems from quantum dynamics under continuous measurement[END_REF].

However, general approaches to adiabatic elimination of Lindblad systems -and maintaining the positivity-preserving quantum structure, beyond a standard linear systems treatment via singular perturbation theory -have attracted much less attention. In [START_REF] Kessler | Generalized Schrieffer-Wolff formalism for dissipative systems[END_REF], Kessler has developed a generalization of the Schrieffer-Wolff formalism; in [START_REF] Gough | Singular perturbation of quantum stochastic differential equations with coupling through an oscillator mode[END_REF][START_REF] Bouten | Approximation and limit theorems for quantum stochastic models with unbounded coefficients[END_REF], the authors address quantum stochastic differential equations in the limit where the speed of the fast system goes to infinity. A geometric approach to adiabatic elimination has been introduced by [START_REF] Azouit | Adiabatic elimination for open quantum systems with effective lindblad master equations[END_REF][START_REF] Azouit | Towards generic adiabatic elimination for bipartite open quantum systems[END_REF], where the authors explore an asymptotic expansion of the reduced dynamics by a careful application of center manifold techniques [START_REF] Carr | Applications of Centre Manifold Theory[END_REF] and geometric singular perturbation theory [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF]. In order to succesfully retain the physical interpretation, the reduced dynamics is expressed by Lindblad equations and is mapped to the original dynamics via a trace-preserving completely-positive (CPTP) map, also called Krauss map.

The present work builds upon the geometric approach of [START_REF] Azouit | Adiabatic elimination for open quantum systems with effective lindblad master equations[END_REF][START_REF] Azouit | Towards generic adiabatic elimination for bipartite open quantum systems[END_REF] and brings forward two novel features. First, unlike in [START_REF] Azouit | Adiabatic elimination for open quantum systems with effective lindblad master equations[END_REF][START_REF] Azouit | Towards generic adiabatic elimination for bipartite open quantum systems[END_REF] where the target system was assumed to be static in the ideal case, we here develop formulas for the case

dρ dt = L(ρ) = -i [H, ρ] + µ D Lµ (ρ),
where each L µ is a "decoherence" operator on H, H is a Hermitian "Hamiltonian" operator on H, and D is a superoperator defined by:

D Lµ (ρ) := L µ ρL † µ - 1 2 L † µ L µ ρ - 1 2 ρL † µ L µ .
In this paper, we consider the composite Hilbert space H := H A ⊗ H B of a target quantum system on H B and its environment on H A . The dynamics on H satisfies a time scale separation:

dρ dt = L A (ρ) + εL int (ρ) + εL B (ρ) + (-i)[ HB , ρ], (1)
where ε is a small positive parameter; L A and L B are Lindbladian super-operators acting exclusively on H A and H B respectively; L int is a Lindbladian superoperator which captures the interaction between H A and H B . Here we assume that this interaction is Hamiltonian and expressed as:

L int (ρ) := -i A ⊗ B † + A † ⊗ B, ρ ,
where A and B respectively are non-necessarily-Hermitian operators acting on H A and H B only. The resonant interaction from of L int models a wide range of applications; general interactions will be addressed by future works. Finally, HB is a Hamiltonian operator on H B , thus expressing fast unitary dynamics on the target system; its presence is the first novelty in our paper. For a set of interesting situations, the dynamics of typical quantum systems can be expressed in a rotating frame where the term HB would vanish. However, several reasons can justify to keep this term. For instance, in many significant situations the vanishing of HB is not rigorous and involves an additional treatment of appearing fast time-varying parameters in the equation via averaging theory; or, HB can be a term of particular interest like a field to be measured with the quantum device or an actuation towards applying some operation on the target system. As a second novelty, we consider a generalized setting where

H A = k H (k)
A is composed of a non-necessarily-finite number of Hilbert spaces

H (k) A . Each subsystem on H (k)
A is strictly dissipative and interacts with H B only. Then, system (1) reads as:

dρ dt = k L (k) A (ρ) + εL (k) int (ρ) + εL B (ρ) + (-i) HB , ρ (2)
where

L (k) A acts on H (k)
A only and where

L (k) int (ρ) := -i A (k) ⊗ B † + A (k) † ⊗ B, ρ ,
captures the Hamiltonian interaction between H (k)

A and H B , with A (k) non-necessarily-Hermitian operators acting on H (k)

A only. The interaction is here restricted to the case of the same operator B for each subsystem k. While the general case will be the subject of future research, having the same operator B for each interaction still models a wide range of applications.

For ε = 0, the system is uncoupled and the solution trajectories stay separable for all times, namely for ρ(0)

= k ρ (k) A (0) ⊗ ρ B (0) we have ρ(t) = k ρ (k)
A (t) ⊗ ρ B (t) for all times, with each factor in the product following its independent dynamics. To apply adiabatic elimination, we assume that each part of the environment is highly dissipative and relaxes fast to a unique steady state, i.e.: for any initial state ρ 0 on H A ⊗ H B , the solution of the uncoupled system ε = 0 converges to k ρ(k) A ⊗ ρ B (t) where, for each k, ρ(k) A is the unique solution of L (k) A ρ(k) A = 0; and ρ B (t) satisfies ρB = -i[ HB , ρ B ] with ρ B (0) = Tr A (ρ 0 ). For ease of presentation, we will also denote ρA := k ρ(k) A .

Asymptotic expansion

Both in the bi-partite and the K-partite case, for the uncoupled system ε = 0, there exists an asymptotically stable center manifold M 0 of same dimension as H B , on which the dynamics have imaginary eigenvalues. It thus follows from Fenichel's Invariant Manifold Theorem [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF] that, for small enough ε > 0, there exists an invariant and attractive manifold M ε which has the same dimension as M 0 and which is ε-close to it. Furthermore, by virtue of linearity and Carr's result [START_REF] Carr | Applications of Centre Manifold Theory[END_REF], M ε is a vector subspace and its approximation can be computed up to arbitrary precision. The quantum particularity, as explained in [START_REF] Azouit | Towards generic adiabatic elimination for bipartite open quantum systems[END_REF], is that such approximation should retain a physical interpretation by preserving the quantum structure: (i) the mapping from the reduced space to the complete space is a mapping between density operators, and M ε can be parameterized by M ε := {ρ ∈ H : ρ = K(ρ s ), ρ s ∈ H s } for some Hilbert space H s that has same dimension as H B , and where K(•) is a Kraus map1 ; (ii) the reduced dynamics on M ε are Lindbladian, i.e. ρs = L s (ρ s ) for some Lindbladian superoperator L s .

In other words, we aim to find a Kraus map ρ = K(ρ s ) and a Lindbladian L s such that the following invariance equation is satisfied for all ε small enough and for all ρ s :

L A (K(ρ s )) + εL int (K(ρ s )) + εL B (K(ρ s )) + (-i)[ HB , K(ρ s )] = K(L s (ρ s )). (3)
By virtue of Carr's result [START_REF] Carr | Applications of Centre Manifold Theory[END_REF], we first parameterize both the Kraus map and the Lindbladian as infinite series:

K(ρ s ) := +∞ h=0 ε h K h (ρ s ), L s (ρ s ) := +∞ h=0 ε h L s,h (ρ s ); (4)
then, by identifying the terms of the same order of ε in the invariance equation ( 3), we obtain an invariance relation at all orders h. At zero-order, we have:

L A (K 0 (ρ s )) + (-i) HB , K 0 (ρ s ) = K 0 (L s,0 (ρ s )) . (5)
Similarly, the first-order invariance condition reads as:

LA (K1(ρs)) + Lint (K0(ρs)) + LB (K0(ρs))

-i HB, K1(ρs) = K0 (Ls,1(ρs)) + K1 (Ls,0(ρs)) , (6)
whereas the second-order invariance condition reads as:

LA (K2(ρs)) + Lint (K1(ρs)) + LB (K0(ρs)) -i HB, K2(ρs) = K0 (Ls,2(ρs)) + K1 (Ls,1(ρs)) + K2 (Ls,0(ρs)) . (7)
For more details about the asymptotic expansion approach to adiabatic elimination, we refer the reader to [START_REF] Azouit | Towards generic adiabatic elimination for bipartite open quantum systems[END_REF][START_REF] Azouit | Adiabatic elimination for open quantum systems with effective lindblad master equations[END_REF].

Reduced-model formulas

The aim of this Section is to provide explicit solutions to the zero-, first-, and second-order invariance equations ( 5)- [START_REF] Breuer | The Theory of Open Quantum Systems[END_REF] for the case of K-partite systems as introduced in Section 2.1, i.e. for model [START_REF] Atkins | Approximate master equations for atom optics[END_REF]. We immediately observe that the zero-order ( 5) is naturally solved by setting:

L s,0 (ρ s ) := -i HB , ρ s , K 0 (ρ s ) := k ρ(k) A ⊗ ρ s . (8)
At first order, let the Kraus map have the following structure inspired by [START_REF] Azouit | Towards generic adiabatic elimination for bipartite open quantum systems[END_REF]:

K(ρs) = K0(ρs) + εK1(ρs) := (I -iεM ) (ρA ⊗ ρs) I + iεM † + O ε 2 , (9)
where

M := k M (k) , M (k) := F (k) 1 ⊗ B † + F (k) 2
⊗ B for any k. This would immediately imply that:

K1(ρs) = -iM (ρA ⊗ ρs) + i (ρA ⊗ ρs) M . (10)
The following assumption will be instrumental in establishing our main results.

Assumption 1 There exists c B † ∈ R such that:

HB, B † = c B † B † . ( 11
)
Theorem 1 Consider model [START_REF] Atkins | Approximate master equations for atom optics[END_REF]. Let Assumption 1 hold. Then, the first-order invariance equation ( 6) is satisfied by the Lindbladian Ls,1(ρs) = LB(ρs) and by a map K1 of the form [START_REF] Carr | Applications of Centre Manifold Theory[END_REF] where, for each k,

F (k) 1 , F (k) 2
respectively are the unique solutions of:

L (k) A F (k) 1 ρ(k) A + A (k) ρ(k) A -i c B † F (k) 1 ρ(k) A =0, (12a)
L (k) A F (k) 2 ρ(k) A + A (k) † ρ(k) A + i c * B † F (k) 2 ρ(k) A =0. (12b)
Furthermore, K(ρs) = K0(ρs) + εK1(ρs) is a CPTP map up to second-order terms.

Proof 1 see Appendix A.1.

Remark 1

The first-order ( 6) is also satisfied by the Lindbladian Ls,1(ρs) = LB(ρs)-i [Hs,1, ρs] with Hs,1

:= k Tr A (k) ρ(k) A B † + Tr A (k) † ρ(k) A B
, and by a map K1 of the form [START_REF] Carr | Applications of Centre Manifold Theory[END_REF] where

F (k) 1 , F (k) 2
respectively are the unique solutions of:

L (k) A F (k) 1 ρ(k) A + S (k) (A (k) ρ(k) A ) -i c B † F (k) 1 ρ(k) A =0, (13a)
L (k) A F (k) 2 ρ(k) A + S (k) (A (k) † ρ(k) A ) + i c * B † F (k) 2 ρ(k) A =0, (13b)
where, for an operator Q acting on

H (k) A , notation S (k) (Q) denotes Q -Tr(Q)ρ (k)
A . Furthermore, K(ρs) = K0(ρs) + εK1(ρs) is a CPTP map up to second-order terms. The possibility of having alternative solutions to the first-order invariance equation hinges upon having c B † = 0, by means of which a gauge degree of freedom in the selection of the trace of terms

F (k) 1 ρ(k) A and F (k) 2

ρ(k)

A can be tuned so as to cancel out Hs,1, as in Theorem 1. It appears that gauge choices are instrumental for positivity-preservation in the solution of the second-order invariance equation, as we consider next.

Theorem 2 Consider model [START_REF] Atkins | Approximate master equations for atom optics[END_REF]. Let Assumption 1 hold. Assume furthermore that Ls,1 = LB = 0, with K1 selected according to Theorem 1. Then, the second-order invariance equation is satisfied by a Lindbladian:

Ls,2(ρs) = k -i ℑ z (k) 1 BB † , ρs -i ℑ z (k) 2 B † B, ρs + 2 ℜ z (k) 1 D B † (ρs) + 2 ℜ z (k) 2 DB(ρs) + k>k ′ -i δ (k,k ′ ) B, B † , ρs , (14)
with:

δ (k,k ′ ) = -2ℜ z (k) 0 z (k ′ ) * 0 c B † , z (k) 0 = Tr A (k) ρ(k) A , (15a)
z (k) 1 = Tr F (k) 1 ρ(k) A A (k) † , z (k) 2 = Tr F (k) 2 ρ(k) A A (k) , (15b)
and by a map K2, obtained from formulas ( 22)-(29), such that K(ρs) = K0(ρs) + εK1(ρs) + ε 2 K2(ρs) is a CPTP map up to third-order terms.

Proof 2 see Appendix A.2. c where ṽ is in the range from 0 to 10 GHz.

Application

Microwave superconducting resonators are an important component in various quantum devices, and in particular in the quantum electrodynamics circuits [START_REF] Blais | Cavity quantum electrodynamics for superconducting electrical circuits: An architecture for quantum computation[END_REF][START_REF] Rau | Cavity quantum electrodynamics in superconducting circuits: Susceptibility at elevated temperatures[END_REF] that are one of the most promising current technologies towards building a quantum computer [START_REF]platform used by Google, Intel[END_REF]. Losses due to imperfections in amorphous materials constitute a dominant loss channel of such resonators [START_REF] Wang | Improving the coherence time of superconducting coplanar resonators[END_REF][START_REF] Gao | A semiempirical model for two-level system noise in superconducting microresonators[END_REF], and can be represented by a bath of two-level systems (TLSs). In many practical cases, strong microwave tones are applied with significant frequency detuning with respect to the resonance frequency [START_REF] Regal | Measuring nanomechanical motion with a microwave cavity interferometer[END_REF] in order to activate a parametric interaction between the resonator mode and another circuit degree of freedom. Within this framework, the LKB team has performed a pump probe experiment [START_REF] Capelle | Energy relaxation properties of a microwave resonator coupled to a pumped two-level system bath[END_REF] on a microwave resonator: a strong "pump" drive, at a frequency far detuned from the resonator, is applied to essentially scramble the quantum behavior ("saturate") of the TLS bath, whereas a weak probe tone, assumed not to disturb the bath behavior, is used to retrieve the transmission spectrum of the resonator. The latter allows to extract induced detuning and damping rate. Let H

Q and HC respectively be the Hilbert space of the k-th TLS=qubit and the resonator, and

HQ := k H (k) Q . Respectively denote with σ (k) + and σ (k)
-the raising and lowering operator on the k-th qubit, and with σ

(k) x , σ (k)
y , and σ (k) z the Pauli operators on the k-th qubit. Let a and a † be the annihilation and creation operators in the resonator mode. The experimental setup is modeled by the following system in Lindblad form:

d dt ρ = -i[H, ρ] + Γ- k D σ (k) - (ρ), H = ωc a † a + ve iωr t + v * e -iωr t a † + a + k ω (k) q 2 σ (k) z + i g σ (k)
x a †a .

Here ωc, ωr, and ω (k) q are the frequencies of the resonator, the pump drive, and the k-th qubit respectively, v is the amplitude of the pump, g is the coupling strength between the resonator and each qubit, and Γ-is the dissipation rate associated to σ (k) on each qubit. The goal would be to obtain a reduced order model for [START_REF] Kessler | Generalized Schrieffer-Wolff formalism for dissipative systems[END_REF] which matches the transmission spectrum of this experiment, as in Figure 1.

For each k, let ∆

(k) q = ω (k) q
-ωr and ∆c = ωc -ωr. Under the assumption that ω

(k) q , ωc, ωr ≫ |∆ (k ′ ) q
|, |∆c|, g, Γ-, Γ+ for any k, k ′ , we apply the standard rotating-wave approximation (i.e. first-order averaging) with H0 := ωr a † a + k ωr 2 σ (k) z the Hamiltonian corresponding to the rotating change of frame and H1 = H -H0 the remaining Hamiltonian. The first-order RWA yields:

d dt ρ rwa 1 = -i [H rwa 1 , ρ rwa 1 ] + Γ- k D σ (k) - (ρ rwa 1 ), where H rwa 1 = ∆c a † a + (va + v * a † ) + k ∆ (k) q 2 σ (k) z + ig σ (k) -a † -σ (k) + a
is the Jaynes-Cummings interaction Hamiltonian plus drive. We next apply a unitary coordinate change ρ = Ũ ρ rwa 1 Ũ † on the resonator state, to center it around its well-known steady state under off-resonant drive, namely by a complex field amplitude displacement Ũ := exp v * a † -va /∆c . This yields:

d dt ρ = k L (k) Q (ρ) + gL (k) int (ρ) + (-i) ∆c a † a, ρ , (16)
where:

L (k) Q (ρ) := -i ∆ (k) q 2 σ (k) z + g ṽ ∆c σ (k) x , ρ + Γ-D σ (k) - (ρ), L (k) int (ρ) := σ (k) -a † -σ (k) + a, ρ .
The term with v := iṽ for ṽ ∈ R now expresses an effective, indirect drive on the TLSs. System ( 16) is in the form (2) with L (k)

A := L (k) Q , A (k) := iσ (k)
-, B := a, HB := ∆ca † a. The hypothesis of Theorems 1 and 2 are satisfied since LB = 0 and Assumption 1 holds with c B † = ∆c. Despite the fact that the microwave resonator in consideration is an infinite-dimensional system and our theory is fully developed for finite-dimensional ones, we still hope to get useful insights by applying our theory to this setup. We will later explain how infinite dimensionality affects our results. For now, let Hs be a Hilbert space whose dimension matches the dimension of the resonator space HC , and ρs the density operator on Hs. By virtue of Theorem 2, the reduced model is given in Lindbladian form as follows:

d dt ρs = -i ∆c + g 2 k δ (k) a † a, ρs (17)
+ g 2 k Γ (k) a Da(ρs) + g 2 k Γ (k) a † D a † (ρs), δ (k) = ℑ z (k) 1 + z (k) 2 , Γ (k)
a † = 2ℜ z (k) 1 , Γ (k) a = 2ℜ z (k) 2 , z (k) 1 = Tr -iF (k) 1 ρ(k) Q σ (k) + , z (k) 2 = Tr iF (k) 2 ρ(k) Q σ (k) -
and where, for each k, matrices

F (k) 1 , F (k) 2
satisfy equations [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF]. The solution of such equations can be computed directly since, on a qubit space

H (k)
Q , one can always parameterize operators in terms of Pauli matrices. We immediately find: z

(k) 1 = W (k) 1 Z (k) and z (k) 2 = W (k) 2 Z (k) where W (k) 1 = -4g 2 v 2 8ig 2 v 2 + (Γ-+ i∆c)∆c (Γ-+ 2i∆c) 2 + 4∆ 2 q , W (k) 2 =32ig 4 v 4 + 2i(Γ--i∆c)∆ 4 c (Γ2 + 4∆ 2 q )• • (Γ--2i(∆c + ∆q)) -4g 2 v 2 ∆c Γ 3 --5iΓ 2 -∆c+ 4i∆c(∆ 2 c + 2∆c∆q -∆ 2 q ) + 4Γ-(-∆ 2 c + ∆ 2 q ) , Z (k) = 8g 2 v 2 + ∆ 2 c (Γ 2 -+ 4∆ 2 q ) 8g 2 v 2 (iΓ-+ 2∆c)+ ∆ 2 c (iΓ-+ ∆c) (Γ--2i∆c) 2 + 4∆ 2 q . Coefficients g 2 k δ (k) , g 2 k Γ (k)
a , and

g 2 k Γ (k)
a † appearing in our reduced-order model ( 17) can be visualized for different values of pump detuning ∆c and intraresonator photon number N = ṽ2 /(κ 2 + ∆ 2 c ). As depicted in Figure 1, we can compare the frequency shift of the resonator g 2 k δ (k) with experimental findings from the pump probe experiment. We find that, by properly calibrating the values of g and by selecting a proper density function of the TLSs, we are able to match the resonance shift in the trasmission spectrum observed in the pump probe experiment. However, we have observed that for large drive gains ṽ the second-order dissipation gives g 2 k Γ

(k)

a † > g 2 k Γ (k) a
which would imply that we constantly keep adding energy and the resonator state drifts off to infinity. Our finite-dimensional treatment can obviously not be trusted in this case. Anyway, quantitative agreement between data and our model should hence enable to extract characteristics about the TLS bath, pending other experimental features that will have to be calibrated.

Conclusions

We have studied adiabiatic elimination for open quantum systems in Lindblad form composed by a target subsystem weakly interacting with K strongly dissipative subsystems. The key novel features of our approach are: (i) the decomposition of the environment into K separately treated subsystems; (ii) the presence of fast Hamiltonian dynamics on the target system. The timescale separation between the uncoupled dynamics and the interaction allows model reduction via geometric singular perturbation theory. We have provided formulas for the first-and second-order expansion and shown that the asymptotic expansion of the center manifold retains a physical interpretation: (i) the reduced model evolves according to Lindbladian dynamics; (ii) reduced and original model are related via Kraus map. Each strongly dissipative subsystem contributes linearly to the reduced model at first-order, and does the same at second-order if a specific commutation property about the interaction terms holds. We have applied our proposed theory to the model of a microwave superconducting resonator subject to dielectric losses where the shape of the trasmission spectrum of our reduced-order model matches experimental data. Future work will address: (i) the necessary conditions to satisfy the invariance equation at orders higher than two; (ii) a thorough study of the infinite-dimensional case; (iii) a full generalization of the proposed theory by removing the assumption about the commutator between the original Hamiltonian dynamics of the target and the interaction terms.

A Proofs of Theorems

A.1 Proof of Theorem 1

Let ρ[k]

A denote k ′ =k ρ(k) A . By plugging ( 8) and ( 10) into the first-order invariance condition [START_REF] Bouten | Approximation and limit theorems for quantum stochastic models with unbounded coefficients[END_REF] and by making use of Assumption 1, condition (6) reads as:

k ρ[k] A ⊗ -iL (k) A F (k) 1 ρ(k) A -iA (k) ρ(k) A -c B † F (k) 1 ρ(k) A ⊗ B † ρs + -iL (k) A F (k) 2 ρ(k) A -iA (k) † ρ(k) A + c B † F (k) 2 ρ(k) A ⊗ Bρs + herm.conj. + L B (ρs) -L s,1 (ρs) ⊗ ρA = 0 (18)
It can be proved along the lines [3, Lemma 4] that equations (A.1) together with trace condition [START_REF] Poyatos | Quantum reservoir engineering with laser cooled trapped ions[END_REF] (or [START_REF] Rau | Cavity quantum electrodynamics in superconducting circuits: Susceptibility at elevated temperatures[END_REF]) are always solvable for

F (k) 1 , F (k)
2 . Case of Theorem 1. It can be immediately seen from (A.1) that one can select Ls,1 := LB as long as the two round parenthesis in (A.1) are set to zero. Taking the trace on (12a) and (12b) yields:

c B † Tr F (k) 1 ρ(k) A = c B † Tr F (k) 2 ρ(k) A * = -i Tr A (k) ρ(k) A , (19)
which solves the situation with the announced formulas.

Case of Remark 1. By taking the trace on equations ( 13), we observe that, for each k:

Tr F (k) 1 ρ(k) A = Tr F (k) 2 ρ(k) A * = 0. (20)
Then, by taking the partial trace over HA in (A.1), we immediately have Ls,1(ρs) = LB(ρs) -i[Hs,1, ρs] with Hs,1 as in Theorem 1. Now, plugging Ls,1 in (A.1) yields:

- k ρ[k] A ⊗ iL (k) A F (k) 1 ρ(k) A + iS (k) A (k) ρ(k) A + c B † F (k) 1 ρ(k) A ⊗ B † ρs + iL (k) A F (k) 2 ρ(k) A + iS (k) A (k) † ρ(k) A -c * B † F (k) 2 ρ(k) A ⊗ Bρs + herm.conj. = 0 (21)
In order to solve [START_REF] Regal | Measuring nanomechanical motion with a microwave cavity interferometer[END_REF] it is enough to zero the two round parenthesis of ( 21) for each k -see equations [START_REF] Gao | A semiempirical model for two-level system noise in superconducting microresonators[END_REF].

It is immediate to see from ( 9) that K0(ρs) + εK1(ρs) is a completely positive map, as long as one can neglect the terms of order ε 2 . One concludes that it is also trace-preserving at order ǫ by checking that Tr (K1(ρs)) = 0 , thanks to [START_REF] Poyatos | Quantum reservoir engineering with laser cooled trapped ions[END_REF] for the case of Theorem 1 and thanks to [START_REF] Rau | Cavity quantum electrodynamics in superconducting circuits: Susceptibility at elevated temperatures[END_REF] for the case of Remark 1.

A.2 Proof of Theorem 2

We prove Theorem 2 in five steps:

1. we define the mapping K up to third-order terms and, by collecting powers of ε, we obtain the formulations of K0, K1, K2; 2. we decompose K2 in three terms and show that they satisfy the second-order invariance condition (7); 3. by taking the partial trace w.r.t. HA on (7), we compute Ls,2;

4. we show that mapping K, according to our definition, is a completely positive mapping; 5. finally, we show trace-preservation of K by proving Tr[K2(ρs)] = 0 for all ρs.

Let ρ[k,k ′ ] A denote k ′′ =k,k ′ ρ(k ′′ ) A . Let ρ(k,k ′ ) A denote ρ(k) A ⊗ ρ(k ′ )
A . We use the following notation for operators: if the superscript of an operator respectively is (k) or (k, k ′ ), then it respectively applies to H

(k,k ′ ) (Q) denote Q -Tr(Q)ρ (k,k ′ ) A . Let B1 := B † , B2 := B, A (k) 1 := A (k) , A (k) 2 := A (k) † . For any k, k ′ , let {U (k,k ′ ) B † j B h } h,j∈{1,2} be a family of four operators on H (k) A ⊗ H (k ′ ) A
only, which we will define in the following. Define:

M = k j∈{1,2} F (k) j ⊗ Bj N := k,k ′ h,j∈{1,2} U (k,k ′ ) B † j B h ⊗ B † j B h W (k) µ (t) :=[L (k) A,µ , F (k) 1 ] ⊗ B † + e 2ic B † t [L (k) A,µ , F (k) 2 ] ⊗ B w (k) (t) :=B1 + b (k) e -2ic B † t B2 b (k) := - 1 2c B † Tr µ [LA,µ, F (k) 2 ] ρ(k) A [LA,µ, F (k) 1 ] † . ( 22
) Let b (k) 1 := b (k) ⋆ and b (k) 2 := b (k) , Let L (k,k ′ ) A denote the operator L (k) A + L (k ′ )
A . Let δ hj := h -j. Let f1, f2 ≥ 0 two constants which we will define in the following. Now define:

g(ρs) := k π/(2c B † ) 0 w (k) (t) ρs w (k) (t) † dt f(ρs) := j∈{1,2} fj Bj ρsB † j G(ρs) := TrA   π/c B † 0 k,µ W (k) µ (t)(ρA ⊗ ρs)W (k) µ (t) † dt   K Q 2 (ρs) := +∞ 0 e L A (•)t S k,µ W (k) µ (t)(ρA ⊗ ρs)W (k) µ (t) † dt + c B † τ π ρA ⊗ G(ρs) + c B † ρA ⊗ g(ρs) + ρA ⊗ f(ρs) K(ρs) := I -iεM + ε 2 N (ρA ⊗ ρs) I + iεM † + ε 2 N † + ε 2 K Q 2 (ρs). (23)
By collecting powers of ε in [START_REF] Sakurai | Modern Quantum Mechanics[END_REF] and carrying out straightforward computations, we obtain the formulation K(ρs) = (K0 + εK1 + ε 2 K2)(ρs) where:

K0(ρs) =ρA ⊗ ρs, , (24)
K1(ρs) = -i k ′ j∈{1,2} ρ[k ′ ] A ⊗ F (k ′ ) j ρ(k ′ ) A ⊗ Bjρs + herm. conj., (25)
K2(ρs) =K L 2 (ρs) + K E 2 (ρs) + K Q 2 (ρs), (26)
where:

K L 2 (ρs) = k,k ′ h,j∈{1,2} ρ[k,k ′ ] A ⊗ U (k,k ′ ) B † j B h ρ(k,k ′ ) A ⊗ B † j B h ρs + herm. conj., K E 2 (ρs) = k,k ′ h,j∈{1,2} ρ[k,k ′ ] A ⊗ F (k) h ρ(k,k ′ ) A F (k ′ ) † j ⊗ B h ρsB † j K Q 2 (ρs) = k h,j∈{1,2} ρ[k,k ′ ] A ⊗ K(k) hj (ρs) ⊗ B h ρsB † j + j∈{1,2} τ Tr Fjj (ρA) (k) ρA ⊗ Bj ρsB † j + k ρA ⊗ π 2 B † ρsB + |b (k) 2 | 2 BρsB † - h =j iδ hj b (k) h B h ρsB † j + + ρA ⊗ f1B † ρsB + f2BρsB † , (27)
and where:

K(k) hj (ρs) := +∞ 0 e t L (k) A (•) S (k) Fhj (t, ρA ) (k) dt Fhj (t, ρA ) (k) := exp 2ic B † δ hj t µ L (k) A,µ , F (k) h ρ(k) A L (k ′ ) A,µ , F (k ′ ) j Fhj (ρ A ) (k) := Fhj (t, ρA ) (k) f 1 := k,k ′ Tr F (k) 2 ρ(k,k ′ ) A F (k ′ ) 2 + δkk ′ F (k) 22 + π 2 |b (k) | 2 f 2 := k,k ′ Tr F (k) 1 ρ(k,k ′ ) A F (k ′ ) 1 + δkk ′ F (k) 11 + π 2 . ( 28
)
where δkk ′ denotes the Kronecker delta. For any k, k ′ , let the family of operators {U

(k,k ′ ) B † j B h } h,j∈{1,2}
satisfy the following set of equations:

L (k,k ′ ) A U (k,k ′ ) B † j B h ρ(k,k ′ ) A + S (k,k ′ ) 2ic B † δ hj U (k,k ′ ) B † j B h ρ(k,k ′ ) A - -A (k) † j F (k ′ ) h ρ(k,k ′ ) A = 0, ∀h, j ∈ {1, 2}, (29a)
Tr 2ic B † δ hj U (k,k ′ ) B † j B h ρ(k,k ′ ) A = = Tr A (k) † j F (k ′ ) h ρ(k,k ′ ) A , ∀h = j ∈ {1, 2}, (29b)
Tr U (k,k ′ ) BB † ρ(k,k ′ ) A = -Tr F (k) 1 ρ(k,k ′ ) A F (k ′ ) † 1 + δkk ′ τ F11(ρA) (k) + π 2 , ( 29c
)
Tr U (k,k ′ ) B † B ρ(k,k ′ ) A = -Tr F (k) 2 ρ(k,k ′ ) A F (k ′ ) † 2 + δkk ′ τ F22(ρA) (k) + π 2 |b (k) | 2 . ( 29d
)
Equations ( 29) are always solvable, as proved in [3, Lemma 4]. We will then show that our definition of K2 in ( 22)-(29) indeed satisfies the second-order invariance condition [START_REF] Breuer | The Theory of Open Quantum Systems[END_REF]. We start by observing that, thanks to assumptions LB(•) = 0 and Ls,1 = 0, condition [START_REF] Breuer | The Theory of Open Quantum Systems[END_REF]

+ B h ρsB † j ⊗ F (k) h ρ(k,k ′ ) A A (k ′ ) † j + A (k) h ρ(k,k ′ ) A F (k ′ ) † j , = k,k ′ ρ[k,k ′ ] A ⊗ j,h∈{1,2} B † j B h ρs ⊗ -A (k) † j F (k ′ ) h ρ(k,k ′ ) A + herm. conj. + B h ρsB † j ⊗ -2ic B † δ hj F (k) h ρ(k,k ′ ) A F (k ′ ) † j -L (k,k ′ ) A F (k) h ρ(k,k ′ ) A F (k ′ ) † j + δk,k ′ Fhj (ρA) (k) , (34)
where, in the last equality, we first made use of formulas ( 12) and then we used the following formula generalized from [3, Lemma 6]:

F (k) h L (k,k ′ ) A ρ(k,k ′ ) A F (k ′ ) † j + L (k,k ′ ) A F (k) h ρ(k,k ′ ) A F (k ′ ) † j = = L (k,k ′ ) A F (k) h ρ(k,k ′ ) A F (k ′ ) † j
δk,k ′ Fhj (ρA) (k) , ∀h, j, ∀k, k ′ . (36)

Figure 1 :

 1 photon number

  ) A only, possibly leaving the remaining subsystems identical; superscripts (k, k) or [k, k] are used interchangeably with superscripts (k) and [k] respectively. Let S

  reads as:LA (K2(ρs)) + Lint (K1(ρs)) -i [ HB, K2(ρs)] -K2([ HB, ρs]) = ρA ⊗ Ls,2(ρs).

  making use of definitions (23)-(27), properties (31), simplification (34), and Claim 1, we can compute:LA(K2(ρs)) -i [ HB, K2(ρs)] -K2([ HB, ρs]) + Lint(K1(ρs)) = = k,k ′ ρ[k,k ′ ] A ⊗ h,j∈{1,2} B † j B h ρs ⊗ L (k,k ′ ) A U (k,k ′ ) B † j B h ρ(k,k ′ ) A + 2ic B † δ hj U (k,k ′ ) B † j B h ρ(k,k ′ ) A -A (k) † j F (k ′ ) h ρ(k,k ′ ) A + herm. conj. + k ρ[k] A ⊗ h,j∈{1,2} B h ρsB † j ⊗ 2c B † |δ hj |b (k) h ρ(k) A + Tr Fhj (ρA) (k) ρ(k) A =: E (ρs).

A Kraus map takes the form ρ = K(ρs) := ℓ M ℓ ρsM † ℓ for some operators M ℓ in order to express any completely positive superoperator[START_REF] Choi | Completely positive linear maps on complex matrices[END_REF], and with ℓ M ℓ ρsM † ℓ = I ensuring trace-preservation i.e. Tr (K(ρs)) = Tr(ρs) = 1.

Then, in order to compute the left-hand side of (30), we observe that the computation of term -i[ HB, K2(ρs)] -K2(Ls,0(ρs)) is simplified by the following set of properties directly implied by Assumption 1:

whereas the computation of term LA(K2(ρs)) is simplified by the following Claim.

Claim 1 L (k) A

K(k)

hj (ρs) + 2ic B † δ hj K(k)

hj (ρs) + S (k) Fhj (ρA) (k) = 0.

Proof 3 Case h = j is proved along the lines of [3, Lemma 1 and Lemma 4]. Case (h, j) = (2, 1) is the hermitian conjugate of case (h, j) = (1, 2) which we are now going to prove. Let L ♯ (•) denote the super-operator L (k)

A , we have that:

for any operator X such that Tr H (k)

A

[X] = 0. First, we formulate K(k) 12 as:

Second, we observe that:

We then conclude from (32) and (33) that:

Furthermore, by considering

we have that:

Observe that the definition of b

h + Tr[ Fhj (ρA) (k) ] = 0 whenever h = j. Furthermore, in the last equality of (34), we observed that Tr[F (k) ] for any h ∈ {1, 2}. The latter two observations and trace condition (29b) will be instrumental in the derivation of Ls,2 in (37). Indeed, by recalling (30), Lindblad Ls,2 can be obtained by taking the partial trace over HA of E (ρs) in expression (36), as follows:

If we now apply definitions [START_REF] Ond Řej Černotík | Adiabatic elimination of gaussian subsystems from quantum dynamics under continuous measurement[END_REF] and property [START_REF] Poyatos | Quantum reservoir engineering with laser cooled trapped ions[END_REF], the expression of Ls,2 in (37) simplifies to:

herm.conj.

which immediately reads as [START_REF] Gough | Singular perturbation of quantum stochastic differential equations with coupling through an oscillator mode[END_REF]. Now, by first subtracting ρA ⊗ Ls,2(ρs) from E (ρs) in expression (36) and then using (29a), we conclude that:

which immediately shows that (30) is satisfied as an identify.

We are now going to prove that K(ρs) in our definition ( 23) is indeed a CPTP mapping. Since terms

already retain the Kraus map form, what remains to prove is complete positivity of

ρs).

Claim 2 There exists τ > 0 such that K Q 2 (ρs) is a completely positive mapping.

Proof 4 Consider an Hilbert basis {|n } 1≤n≤d of HA ⊗ HB. Let H be any Hilbert space of finite dimension. Let:

For each n and ν, select any |φn , |ψν ∈ H and define:

We are then going to prove that there exists τ > 0 such that Ψ| K(|Φ Φ|) |Ψ ≥ 0. Standard computations give:

where zn,µ := φn|ψν and

(t) is a completely-positive superoperator on X, the d 2 × d 2 Hermitian matrix r n ′ ,ν ′ ,n,ν (t) is non-negative, and therefore, for any vector z ∈ C d 2 :

Now take a d 2 vector z such that (42) is satisfied. We then have:

Since the propagator e tL A (•) is a completely-positive mapping of the form e tL A (•) (X) = θ H θ (t)XH(t) † θ for some operators H θ (t), we then have [3, Lemma 1]:

≥ 0, and we thus conclude that z † Rz = 0 implies z † M z ≥ 0. Assume now that z † Rz = 0 and z † M z = 0. Inequality (44) then implies that for any t ≥ 0

and therefore M z = 0. Then, by virtue of [3, Lemma 2], we conclude that there exists τ > 0 such that Kτ is completely positive.

We are now going to prove trace-preservation of K(ρs). Since Tr[K1(ρs)] = 1 and Tr[K1(ρs)] = 0 for all ρs, what remains to prove is Tr[K2(ρs)] = 0 for all ρs. First, by a subsequent application of formulas (29b), (35), and (12), we have:

Secondly, it is straightforward to prove from (29c) and (29d

Finally, thanks to observations (45), (47) and the definition of f1, f2 in (28), we compute the trace of K2 from (27) as follows:

for all ρs. We bring to the attention of the reader that the gauge choices in K2 have been carefully designed so as to yield crucial quantum properties:

• τ is selected according to Claim 2 to yield complete positivity of K2;

• the trace of operators U

B † B † are selected according to (29b) to yield cancellation of the BBρs and B † B † ρs terms and their hermitian conjugates in the invariance equation, thus ensuring the Lindblad form of Ls,2;

• constants b (k) are selected according to [START_REF] Reiter | Effective operator formalism for open quantum systems[END_REF] to yield cancellation of the BρsB and B † ρsB † terms in the invariance equation, thus ensuring the Lindblad form of Ls,2;

• the definition of f1, f2 and the trace of operators U

BB † are selected according to (28) and (29c)-(29d) respectively in order to cancel out the positive terms in the trace of K2 corresponding to quadratic terms, as shown in (47), thus ensuring trace-preservation of K.


Gerardo Cardona

Alain Sarlette
email: alain.sarlette@inria.fr

Pierre Rouchon

Exponential stochastic stabilization of a two-level quantum system via strict Lyapunov control

This article provides a novel continuous-time state feedback control strategy to stabilize an eigenstate of the Hermitian measurement operator of a two-level quantum system. In open loop, such system converges stochastically to one of the eigenstates of the measurement operator. Previous work has proposed state feedback that destabilizes the undesired eigenstates and relies on a probabilistic analysis to prove convergence. In contrast, we here associate the state observer to an adaptive version of so-called Markovian feedback (essentially, proportional control) and we show that this leads to a global exponential convergence property with a strict Lyapunov function. Furthermore, besides the instantaneous measurement output, our controller only depends on the single coordinate along the measurement axis, which opens the way to replacing the full state observer by lower-complexity filters in the future.

Introduction

After the technological developments of the last decade, we are now at a stage where measuring and controlling quantum systems is experimentally realizable, as in the first real time quantum feedback experiment [START_REF] Sayrin | Real-time quantum feedback prepares and stabilizes photon number states[END_REF]. From this pioneering work, involving essentially discrete-time logic, the physics community has been moving towards engineering designs that are truly promising for quantum computation [START_REF] Campagne-Ibarcq | Using spontaneous emission of a qubit as a resource for feedback control[END_REF][START_REF] Devoret | Superconducting circuits for quantum information: an outlook[END_REF][START_REF] Liu | Comparing and combining measurementbased and driven-dissipative entanglement stabilization[END_REF] and which pose control problems in terms of continuoustime stochastic systems.

A basic element in quantum control is the use of quantum non-demolition (QND) measurements [START_REF] Haroche | Exploring the Quantum: atoms, cavities and photons[END_REF]. This is essentially a continuous-time version of the projection postulate, where performing a continuous measurement makes the quantum state progressively converge to a random eigenstate of the measurement operator. Each such eigenstate is a steady-state of the dynamics, and hence remains unperturbed under the backaction associated to this quantum measurement.

QND eigenstates are thus natural equilibria of a measured quantum system, and stabilizing such a system towards one target eigenstate thanks to an appropriate feedback law would be a basic building block towards more involved control procedures. Several papers have indeed considered ways to stabilize a target QND eigenstate, see e.g. [START_REF] Van Handel | Feedback control of quantum state reduction[END_REF][START_REF] Yamamoto | Feedback control of quantum entanglement in a twospin system[END_REF][START_REF] Mirrahimi | Stabilizing feedback controls for quantum systems[END_REF][START_REF] Tsumura | Global stabilization of ndimensional quantum spin systems via continuous feedback[END_REF][START_REF]Global stabilization at arbitrary eigenstates of n-dimensional quantum spin systems via continuous feedback[END_REF]. They prove global asymptotic convergence with probabilistic arguments, and based on a state feedback controller assuming that a state observer (often called the quantum Bayes filter) efficiently captures the evolution of the quantum state.

In our sense these results leave room for improvement in two directions. First, the feedback laws proposed in existing work are quite complicated, especially their implementation based on a quantum state observer would scale poorly with increasing system dimension as will be needed in advanced quantum technology. Second, the associated convergence analysis is rather involved and none of these feedback schemes are shown to ensure exponential convergence. It is wellknown that exponential convergence is an indication of robustness when the system would interact with other subsystems or would perform subject to some perturbations. In other basic quantum settings, exponential convergence does hold with rather direct proofs. This is the case for the stabilization of target states that are not eigenstates of the measurement operator, via so-called Markovian feedback [START_REF] Wiseman | Quantum theory of continuous feedback[END_REF][START_REF] Wiseman | Quantum Measurement and Control[END_REF]. This essentially comes down to proportional control, and in [START_REF] Wiseman | Bayesian feedback versus markovian feedback in a two-level atom[END_REF] it was shown how the closed-loop equation can advantageously be reformulated as having changed the dissipation operators to have the target as steady state. A drawback of this approach is that, under non-ideal measurement conditions, a steady state far from the QND eigenstate would be subject to significant noise; a state close to a QND eigenstate in contrast would have little noise, but the Markovian feedback gain and the associated convergence rate towards such state approaches zero as the target approaches a QND eigenstate. Nonwithstanding, the convergence of the openloop system under QND measurement towards the set of its steady states can also be shown to be exponential. The absence of a proven similar property for the selection of one target QND steady state thus appears as an avoidable gap.

In the present paper, we identify a solution to the latter problem, i.e. exponential stabilization of a QND eigenstate at least for a qubit system, and we pave the way towards addressing the first issue, i.e. the real-time implementation complexity associated to a full state observer. Our solution in fact combines the state feedback ideas of [START_REF] Van Handel | Feedback control of quantum state reduction[END_REF][START_REF] Yamamoto | Feedback control of quantum entanglement in a twospin system[END_REF][START_REF] Mirrahimi | Stabilizing feedback controls for quantum systems[END_REF][START_REF] Tsumura | Global stabilization of ndimensional quantum spin systems via continuous feedback[END_REF][START_REF]Global stabilization at arbitrary eigenstates of n-dimensional quantum spin systems via continuous feedback[END_REF] with the Markovian feedback in the sense of Wiseman [START_REF] Wiseman | Quantum theory of continuous feedback[END_REF]. It should be noted that the idea of using this kind of dynamical feedback is also considered in [START_REF] Zhang | Locally optimal measurement-based quantum feedback with application to multi-qubit entanglement generation[END_REF][START_REF] Xu | Measurement of the topological chern number by continuous probing of a qubit subject to a slowly varying hamiltonian[END_REF], but for different goals. As far as we know our proposal is the first one that provides a convergence proof ensuring almost-sure exponential asymptotic stability with a strict Lyapunov function in the stochastic sense. The interpretation of the control law comes down to an adaptive version of the proportional feedback gains, depending on the distance of the measured QND coordinate from its target value. This coordinate is still assumed to be estimated by a perfectly converging quantum filter for now, but in the future it seems natural to replace it by lower-order filters.

The structure of the paper is as follows. In section 2 we set up some preliminaries on quantum systems under continuous measurements in finite dimensions. We recall the asymptotic behavior for the open loop system and provide an exponential convergence rate for a general QND system, via an original Lyapunov function, under a non-degeneracy condition. In section 3 we review static output feedback (Markovian feedback). We explain how to tune the feedback gain as a function of the target steady state. We provide a Lyapunov function that proves the exponential convergence towards any chosen target, except that the convergence rate tends to 0 for a QND eigenstate as target. This highlights the necessity of another approach for QND eigenstates. In section 4 we consider state feedback control, we review some of the existing results on global stabilization for QND eigenstates and discuss the main difficulties related to this kind of approach. In section 5 we present our main result, i.e. the novel state feedback controller and its exponential convergence proof.

QND diffusive measurements

We consider the basic open-loop model for a quantum system under Quantum Non Demolition (QND) continuous measurements [START_REF] Barchielli | Quantum trajectories and measurements in continuous time: the diffusive case[END_REF] in finite dimension 1 < n < ∞ :

dρ = D(L, ρ)dt + √ ηM(L, ρ)dW, (1)
dY = √ η Tr(Lρ + ρL † )dt + dW, (2)
where ρ belongs to the set of density matrices S = {ρ ∈ C n×n : ρ ≥ 0, ρ = ρ † , Tr(ρ) = 1}, Tr( • ) denotes the trace, W is a standard Wiener process and dY is the measurement process, L = L † is a Hermitian matrix in C n×n . The general properties of such stochastic differential equations (1) are well known: existence and uniqueness of solutions in the state-space S. We have used the following usual notations:

D(L, ρ) = LρL † - 1 2 L † Lρ - 1 2 ρL † L, (3)
M(L, ρ) = Lρ + ρL † -Tr((L + L † )ρ)ρ. (4)
Throughout the paper we assume that L = L † admits a non degenerate spectrum: all its eigenvalues (λ ℓ ) 1≤ℓ≤n are distinct and its normalized eigenvectors (| ψℓ ) 1≤ℓ≤n , with the bra-ket notation (see, e.g., [START_REF] Nielsen | Quantum computation and quantum information[END_REF]), form an orthonormal frame. Such QND systems behave in a particular way: each eigenvector of L, thus ρℓ | ψℓ ψℓ |, is a steady-state of Eq. (1); using Itô calculus and taking the expectation we have

d dt E(Tr(ρρ ℓ )) = 0, (5)
where Tr(ρρ ℓ ) = ψℓ |ρ| ψℓ corresponds to the so-called trace fidelity, reaching its maximum 1 only when ρ = ρℓ . Thus Tr(ρρ ℓ ) is a martingale for (1). Moreover as shown in [START_REF] Van Handel | Feedback control of quantum state reduction[END_REF][START_REF] Mirrahimi | Stabilizing feedback controls for quantum systems[END_REF][START_REF] Benoist | Large time behavior and convergence rate for quantum filters under standard non demolition conditions[END_REF], all trajectories converge to one of the above steady states, with Tr(ρ t ρℓ ) giving the probability to end up on ρℓ conditioned on knowing ρ t (in particular, for t = 0). We propose in the following lemma a proof of this asymptotic result where we estimate the convergence rate via an original exponential Lyapunov function.

L| ψℓ = λ ℓ | ψℓ , 1 ≤ ℓ ≤ n}. Let V (ρ) = n ℓ ′ =1 ℓ<ℓ ′ Tr(ρρ ℓ ) Tr(ρρ ℓ ′ ). Then E(V (ρ t ) | ρ 0 ) ≤ e -rt V (ρ 0 )
for all t ≥ 0, with rate r ≥ 0 given by

r = 1 2 η min ℓ ′ ,ℓ =ℓ ′ (λ ℓ -λ ℓ ′ ) 2 . Proof. Let ξ ℓ := Tr(ρρ ℓ ) for 1 ≤ ℓ ≤ n. It satisfies dξ ℓ = -1 2 η(λ ℓ -̟(ξ))) 2 ξ ℓ dt + √ η(λ ℓ -̟(ξ)))ξ ℓ dW (6) with ξ = (ξ ℓ ′ ) 1≤ℓ ′ ≤n and ̟(ξ) = ℓ ′ λ ℓ ′ ξ 2 ℓ ′ . With the ξ coordinates, V becomes a C 2 function of ξ, V (ξ) = ℓ ′ , ℓ<ℓ ′ ξ ℓ ξ ℓ ′
, with ξ obeying the system of n stochastic differential equations ( 6) driven by a common scalar Wiener process W . We can directly compute the Markov generator A of the ξ-system on V using the formula (15) recalled in appendix. Standard algebraic computations yield

AV = -η 2 n ℓ ′ =1 ℓ<ℓ ′ (λ ℓ -λ ℓ ′ ) 2 ξ ℓ ξ ℓ ′ .
Since each component of ξ(t) remains non-negative for all t, we have AV ≤ -η 2 min ℓ ′ ,ℓ =ℓ ′ (λ ℓ -λ ℓ ′ ) 2 V and we conclude with Gronwall's inequality.

Since V (ρ) = 0 is equivalent to {ρ = ρℓ for some ℓ}, standard results on stochastic stability imply that, for any initial state ρ 0 ∈ S, the solution ρ t of (1) converges almost surely to one of the states ρℓ , as proved independently in [START_REF] Van Handel | Feedback control of quantum state reduction[END_REF][START_REF] Mirrahimi | Stabilizing feedback controls for quantum systems[END_REF].

As explained in [START_REF] Stockton | Deterministic dicke-state preparation with continuous measurement and control[END_REF], such QND diffusive measurement can be viewed as preparing a pure state among the | ψℓ in a non-deterministic manner. As recalled in our introduction, several papers have proposed feedback controllers to actually stabilize one target steady state (see e.g. [START_REF] Van Handel | Feedback control of quantum state reduction[END_REF][START_REF] Yamamoto | Feedback control of quantum entanglement in a twospin system[END_REF][START_REF] Mirrahimi | Stabilizing feedback controls for quantum systems[END_REF][START_REF] Tsumura | Global stabilization of ndimensional quantum spin systems via continuous feedback[END_REF]), establishing global convergence with probabilistic arguments. In section 5 we propose a new controller, based on adapting so-called quantum Markovian feedback [START_REF] Wiseman | Quantum theory of continuous feedback[END_REF], for which we can prove exponential convergence with a Lyapunov function like in the open-loop case. In the next section we therefore first recall the stabilization properties of such Markovian feedback, or static output feedback, which are in fact complementary to our goal.

Static output feedback for QND systems

When adding a control input to the system, we have:

dρ = exp(-iF udt) (ρ + D(L, ρ)dt + √ ηM(L, ρ)dW ) exp(iF udt) -ρ . ( 7
)
Here u is a (real) control signal and F is a Hermitian matrix. We use the propagator notation to emphasize that when u dt contains Wiener processes, as we will do, exp(-iF udt) must be expanded to second order in a causal Itō calculus. In analogy to classical systems, consider the simplest control scheme, a proportional output feedback of the form

udt = f dt + κdY, (8)
where κ is a constant gain and f a constant bias. Following [START_REF] Wiseman | Quantum theory of continuous feedback[END_REF] and Itō rules, the static controller (8) leads to the following closed-loop stochastic master equation:

dρ = -if [F, ρ]dt - iκ √ η 2 [F L + LF, ρ]dt + D(L -iκ √ ηF, ρ)dt + D(iκ 1 -ηF, ρ)dt + √ ηM(L, ρ)dW -iκ[F, ρ]dW. ( 9
)
The next lemma shows that the controller (8) cannot be tuned to achieve global asymptotic stability towards any QND eigenstate ρℓ .

Lemma 3.1. There exists no combination of F , and constants f and κ such that the proportional output feedback (8) would yield a closed-loop dynamics (9) which globally converges towards a target QND eigenstate ρℓ for some chosen ℓ ∈ {1, . . . , n}.

Proof. We can rewrite Eq. ( 9) as

dρ = -if [F, ρ]dt -iκ √ η[F, Lρ + ρL]dt + D(L, ρ)dt + κ 2 D(F, ρ)dt + √ ηM(L, ρ)dW -iκ[F, ρ]dW.
Then if ρℓ is a steady-state of ( 9) and since f, κ are constants, then This obstruction seems to only appear for QND eigenstates ρℓ . The next lemma indeed shows that, at least for a two-level system, any other pure state can be exponentially stabilized by a static output feedback with a fixed measurement operator L and detection efficiency η = 1. We use the standard Pauli matrix and Bloch sphere notation for the qubit system. Lemma 3.2. Consider (9) with n = 2, η = 1, and

-if [F, ρℓ ]-iκ √ η[F, L ρℓ + ρℓ L]+D(L, ρℓ )+ κ 2 D(F, ρℓ ) =
L = Γ 2 σ z = Γ 2 1 0 0 -1 , F = σ y = 0 -i i 0 . Take θ / ∈ {kπ : k ∈ Z} and set f = -Γ 2 sin( θ) cos( θ), κ = Γ 2 sin( θ).
Then the closed-loop system exponentially stabilizes the pure state:

ρ = 1 2 1 + cos( θ) sin( θ) sin( θ) 1 -cos( θ) , as the Lyapunov function V (ρ) = 1 -Tr(ρρ) decreases according to E(V (ρ t )) = e -rt V (ρ 0 ) ∀t ≥ 0 with rate r = Γ(sin θ) 2 .
The proof, left to the reader, is based on standard matrix manipulations showing dE(V ) = -Γ(sin θ) 2 E(V ) dt. By modifying the actuation Hamiltonian as F = U ᾱσ y U † ᾱ, where U ᾱ = exp(-i ᾱσ z /2) is a rotation of angle ᾱ around the z axes of the Bloch associated to the qubit, the same feedback gains stabilize U ᾱ ρU † ᾱ. With ᾱ ∈ [0, 2π] and θ ∈]0, π[, any pure state different from the two eigenstates of σ z are thus obtained via U ᾱ ρU † ᾱ.

For more complex systems, with n > 2, we can use the algebraic criterion elaborated in [START_REF] Ticozzi | Quantum markovian dynamics: Invariant subsystems, attractors, and control[END_REF][START_REF]Analysis and synthesis of attractive quantum markovian dynamics[END_REF][START_REF] Ticozzi | Stabilizing quantum states by constructive design of open quantum dynamics[END_REF] to analyze global convergence towards a unique pure state under such static output feedback. Thanks to positivity of ρ, this can be done by looking at the ensemble average closed-loop dynamics

dρ/dt = -if [F, ρ] -iκ 2 [F L + LF, ρ] + D(L -iκF, ρ
) and its convergence towards pure states. We think that at least when the goal pure state ρ is such that Tr(ρρ ℓ ) > 0 for all ℓ ∈ {1, . . . , n}, then a control Hamiltonian F can be chosen among a "reasonable" set and f and κ exist which ensure global exponential convergence of (9) towards ρ. At this point we leave the question of feedback constrained to typical available F for n > 2, and we move on to the stabilization of the elusive targets, that is the QND eigenstates ρℓ .

u dt = f (ρ) dt (10)
and no direct feedback of dY . Although they all start with Lyapunov control techniques and the open-loop martingale V (ρ) = 1 -Tr(ρρ ℓ ), they also involve specific stochastic arguments combining support of the closed-loop trajectories with Doob inequalities in order to establish the global asymptotic almost sure convergence for the closed-loop system [START_REF] Mirrahimi | Stabilizing feedback controls for quantum systems[END_REF][START_REF] Tsumura | Global stabilization of ndimensional quantum spin systems via continuous feedback[END_REF]. Let us explain the basic reason why it is difficult, as already highlighted in [START_REF] Van Handel | Feedback control of quantum state reduction[END_REF], to construct quantum state feedback and a strict closed-loop Lyapunov function in this case. With ( 7), [START_REF]Global stabilization at arbitrary eigenstates of n-dimensional quantum spin systems via continuous feedback[END_REF], we have E(V t+dt |ρ t , u t ) = V t + u t Tr(i[F, ρ t ]ρ ℓ )dt. The quantum state feedback u t = -Tr(i[F, ρ t ]ρ ℓ ) makes this expectation decreasing. However any QND state ρℓ ′ with ℓ ′ = ℓ is a closed-loop steady state since Tr([F, ρℓ ′ ]ρ ℓ ) = 0, for any choice of F . In a stochastic convergence setting, such steady states are harder to treat, as significant portions of the state can even converge to maxima of a stochastic Lyapunov function. Researchers have proposed several ways to solve this problem.

In all pure Lyapunov feedback approaches like [START_REF] Van Handel | Feedback control of quantum state reduction[END_REF][START_REF] Yamamoto | Feedback control of quantum entanglement in a twospin system[END_REF], the presence of several steady states remains, and an additional proof element e.g. based on sending undesired equilibria to infinite values of the Lyapunov function is used to obtain almost global stability. In [START_REF] Mirrahimi | Stabilizing feedback controls for quantum systems[END_REF] and in [START_REF] Tsumura | Global stabilization of ndimensional quantum spin systems via continuous feedback[END_REF], the controllers are explicitly designed to apply perturbations which remove the undesired steady states. In [START_REF] Mirrahimi | Stabilizing feedback controls for quantum systems[END_REF] this is done with a discontinuous control that gives a constant input whenever ρ is too close to the bad QND eigenstates ρℓ ′ and then switching to the above Lyapunov feedback once ρ is close enough to the goal state ρℓ . In [START_REF] Tsumura | Global stabilization of ndimensional quantum spin systems via continuous feedback[END_REF] a continuous approach is proposed where the control is perturbed by the addition of a term that only vanishes on the target state. More precisely, for the two-level system of lemma 3.2, a smooth feedback law

u = -α Tr(i[σ y , ρ]σ z ) + β(1 -Tr(ρσ z )), (11)
with α, β > 0 such that β 2 8αη < 1, is applied to stabilize the excited state ρ1 = (1 + σ z )/2. The convergence proof then involves a combination of Lyapunov function on part of the state space, and probabilistic arguments on the other part. With the feedback [START_REF] Wiseman | Quantum theory of continuous feedback[END_REF], the set {ρ ∈ S : V (ρ) ≤ β 2 8αη } satisfies dE(V ) ≤ 0 defining a kind of attraction set for the control Lyapunov function. Thus trajectories that do not exit this region of attraction will converge to the target. Outside the region of attraction, both feedback schemes rely on ensuring that trajectories will reach again this region of attraction. We refer to [START_REF] Mirrahimi | Stabilizing feedback controls for quantum systems[END_REF][START_REF] Tsumura | Global stabilization of ndimensional quantum spin systems via continuous feedback[END_REF][START_REF]Global stabilization at arbitrary eigenstates of n-dimensional quantum spin systems via continuous feedback[END_REF] for more details.

All those methods only establish asymptotic convergence. Having a stronger global convergence result like exponential convergence is important for robustness and estimation of convergence speed. The robustness is practically important in particular towards unmodeled dynamics. This could concern other system elements like actuators, but at very least there is the quantum filter which estimates the state ρ from the measurement inputs with a finite convergence speed. The convergence speed is particularly important in applications where feedback control is applied to protect the fragile quantum systems from perturbing effects, like decoherence (see simulations section). While our paper was under review, a parallel paper [START_REF] Liang | On exponential stabilization of spin-1/2 systems[END_REF] has proposed an adaptation of the feedback and of several stochastic techniques used in [START_REF] Tsumura | Global stabilization of ndimensional quantum spin systems via continuous feedback[END_REF]. They obtain global asymptotic convergence, which is further proved to be asymptotically exponential, but without giving guarantees at finite time. As far as we know, a feedback scheme ensuring global exponential convergence via simple Lyapunov argument has remained an open issue.

Exponential feedback stabilization of two-level QND systems

We present now our main result. We consider [START_REF] Yamamoto | Feedback control of quantum entanglement in a twospin system[END_REF] with

L = Γ 2 σ z , F = σ y
, where σ z , σ y are as in Lemma 3.2. When u = 0, Lemma 2.1 shows that ρ t converges as t → ∞ to the set of states {|e e|, |g g|} corresponding to the eigenvalues 1 and -1 respectively of σ z . We propose a feedback law combining quantum-state and output feedback, in the form:

u t dt = f (ρ t )dt + κ(ρ t )dY ( 12
)
where f and κ are regular functions of ρ such that the closed-loop dynamics

dρ t = -if (ρ)[σ y , ρ]dt+D Γ 2 σ z -i √ ηκ(ρ)σ y , ρ dt + D i (1 -η)κ(ρ)σ y , ρ dt + ηΓ 2 M(σ z , ρ)dW -iκ(ρ)[σ y , ρ]dW (13)
is well-posed. For any realization ρ t , if ρ 0 belongs to the compact set S of density operators, then ρ t remains in S for all t > 0. We will target the stabilization of the state |e e|; the treatment for target |g g| is exactly similar. Theorem 5.1 below shows that f and κ can be chosen in order to have an explicit strict closedloop Lyapunov function, converging in average exponentially to zero, and implying the global convergence of the state towards |e e|.

Theorem 5.1. Consider equation [START_REF] Wiseman | Bayesian feedback versus markovian feedback in a two-level atom[END_REF] with

f (ρ) = -ηΓ 2 (1-Tr(ρσ z ) 2 ) and κ(ρ) = ηΓ 2 (1-Tr(ρσ z )). Set V (ρ) = 1 -Tr(ρσ z ).
Then |e e| is globally exponentially stable (in the sense of definition 7.1 in appendix with p = 1) and

E(V (ρ t ) | ρ 0 ) ≤ e -rt V (ρ 0 ) for all t ≥ 0 and ρ 0 ∈ S, with convergence rate r = ηΓ 4 . Proof. Let z = Tr(ρσ z ), x = Tr(ρσ x ), y = Tr(ρσ y ) and set ξ = √ 1 -z. We express f = -ηΓ 2 (2 -ξ 2 )ξ 2 , κ = ηΓ 2 ξ 2 .
After the change of coordinates (x, y, z) → (x, y, ξ), system (13) satisfies:

dx = ηΓ(2-ξ 2 )ξ 2 (1-ξ 2 ) + 2ηΓξ 2 -ηΓξ 4 x dt + 2Γη((1 -ξ 2 )(ξ 2 -x))dW, dy = -Γydt -2Γη(1 -ξ 2 )ydW, dξ = ηΓ 2 (2 -ξ 2 )x -ξ 2 (1 -ξ 2 ) -1 2 ((2 -ξ 2 ) -x) 2 ξdt -ηΓ 2 ((2 -ξ 2 ) -x)ξdW.
We now compute the Markov generator A of the Lyapunov function V (x, y, ξ) = ξ which yields

AV = ηΓ 2 -ξ 2 (1 -ξ 2 ) -1 2 ((2 -ξ 2 ) 2 + x 2 ) ξ = -ηΓ 2 ( 1 2 + 3 2 z 2 + 1 2 x 2 )ξ ≤ -ηΓ 4 ξ. Since dE(V t ) = E(AV (ξ))dt, by Gronwall's inequality E(V t | V 0 ) ≤ e -rt V 0 .
To establish convergence of the state, we use theorem 7.1 in the appendix. The Lyapunov function V is non-negative and it is a C 2 function of (x, y, ξ). It is obviously upper bounded by V = ξ < k 2 ( x 2 + y 2 + ξ 2 ) p with k 2 = p = 1. For a lower bound, we must take into account that the system evolves in the set S = {x, y, ξ :

x 2 + y 2 + (1 -ξ 2 ) 2 ≤ 1}; a direct development yields x 2 + y 2 ≤ ξ 2 (2 -ξ 2 ) ≤ 2ξ 2 , such that V = ξ ≥ ( (x 2 + y 2 + ξ 2 )/3) p = k 1 ( x 2 + y 2 + ξ 2 ) p with p = 1 and k 1 = 1/ √ 3.
The latter also implies that AV (ξ) ≤ -k 3 ( x 2 + y 2 + ξ 2 ) p with p = 1 and k 3 = rk 1 . Thus all the assumptions of theorem 7.1 are satisfied and we conclude that (x, y, ξ) = (0, 0, 0) is exponentially stable in the sense of definition 7.1 in appendix with p = 1. With Lemma 3.2 and Theorem 5.1 we have achieved full exponential stabilization of a two-level system. The static output feedback of Lemma 3.2 can stabilize all states of a qubit except a QND eigenstate; the feedback of Theorem 5.1 stabilizes the QND eigenstates. The overhead cost in Theorem 5.1 with respect to static output feedback is that an estimator is needed for the state component Tr(ρσ z ). Interestingly, this is much less than the full quantum state ρ, and it corresponds to the component Tr(Lρ + ρL † )dt that is directly present in the measurement output dY with L ∝ σ z . This shows that besides allowing a more standard exponential convergence proof, our controller also possibly opens the door to simpler feedback schemes than fullstate-feedback.

Regarding convergence speed, the result of Theorem 5.1 is at most 4 times slower than the open loop convergence towards a random eigenstate (Lemma 2.1). This is much faster than the slow convergence of Lemma 3.2 for target states close to |e e|, and one may envision to speed up the latter by adapting the strategy of Theorem 5.1.

Numerical test

We test the results of Theorem 5.1 through numerical simulations. We fix η = 0.5 and consider 100 trajectories to sample the stochastic behavior. In fig. 1 we show the behavior of the closed loop system (13) with initial condition ρ 0 = I/2. The simulation appears to confirm the exponential convergence. In fig. 2, we show the behavior when we add a standard perturbation to the nominal behavior, namely a decoherence term of the form Γ 2 D(|g e|, ρ) that expresses energy loss of the excited state |e e| by spontaneous photon emission. By varying Γ ∈ {0, Γ/100, Γ/10}, we see that the Lyapunov function goes from about 0 to about 0.1 and 0.5, which corresponds to z = 1, z = .99 and z = .75 respectively.

Conclusion

A quantum system under quantum non-destructive (QND) measurements is a particular stochastic system, where any eigenstate of the measurement operator is a steady state of the stochastic dynamics. While stabilizing a particular eigenstate under continuous-time QND measurement appears as a basic problem of quantum control, it is also remarkably different from standard control problems. Indeed, we have first considered a proportional output feedback (see Lemma 3.2), with which we recall that we can stabilize all target states except the QND steady states. This indicates that for the QND eigenstates as targets, one must make a more involved controller in order to bias the stochastic evolution towards the desired extremum. Previous work like [START_REF] Van Handel | Feedback control of quantum state reduction[END_REF][START_REF] Yamamoto | Feedback control of quantum entanglement in a twospin system[END_REF][START_REF] Mirrahimi | Stabilizing feedback controls for quantum systems[END_REF][START_REF] Tsumura | Global stabilization of ndimensional quantum spin systems via continuous feedback[END_REF] has turned to specific stochastic controller designs to achieve this goal, on the basis of a quantum state observer. In this paper we have shown that just adapting the proportional output feedback gain, as a function of estimated quantum state, is sufficient to obtain global exponential stabilization, something that previous state feedback controls do not achieve or at least could not prove. In our case, the proof is a direct consequence of standard stochastic convergence theorems with a strict Lyapunov function. Interestingly, this controller architecture can be re-interpreted as indirectly engineering state-dependent dissipation operators.

While this provides a simple and well-proved state feedback law, the complexity of its verbatim implementation would still be dominated by the quantum state observer, which is needed to estimate the quantum state on the basis of the measurement outputs. However, our controller depends on the qubit state through only the z coordinate, which is the measured coordinate. This paves the way towards using simpler filters to extract accurate z information from the mea-surement output, and leading to more scalable controllers. The extension of this technique to higherdimensional systems is indeed a subject of ongoing work, with promising preliminary results.

Appendix

We recall some results from stochastic stability and refer the reader to [START_REF] Kushner | Stochastic stability[END_REF][START_REF]Stochastic Stability and Control[END_REF][START_REF] Khasminskii | Stochastic stability of differential equations[END_REF] for the proof of these results and for further reference. We consider diffusion processes x t on R n . They correspond to solution of Itō stochastic differential equations of the form dx = µ(x)dt + σ(x)dW, [START_REF] Zhang | Locally optimal measurement-based quantum feedback with application to multi-qubit entanglement generation[END_REF] where the µ, σ are regular functions of x that satisfy the usual conditions for existence and uniqueness of SDE's [START_REF] Arnold | Stochastic Differential Equations: Theory and Applications[END_REF]. For a C 2 function V the Markov generator associated with ( 14) is

AV = i µ i ∂ ∂x i V + 1 2 i,j σ i σ j ∂ 2 ∂x i x j V, (15)
and E(V (x t ) |V (x 0 )) = V (x 0 ) + E t 0 AV (x s )ds .

In this appendix we assume that S, a compact subset of R n , is positively invariant. We use the following specific asymptotic stability definition.

Definition 7.1 (Khasminskii [START_REF] Khasminskii | Stochastic stability of differential equations[END_REF]). Consider the diffusion process on S governed by [START_REF] Zhang | Locally optimal measurement-based quantum feedback with application to multi-qubit entanglement generation[END_REF] with 0 ∈ S, µ(0) = σ(0) = 0 and p > 0. Consider the Euclidean norm on R n : x = i x 2 i . The equilibrium solution x t = 0 is said exponentially p-stable if, for some constants C > 0 and r > 0 such that

E( x t p | x 0 ) ≤ C x 0 p e -rt .
The stochastic counterpart of Lyapunov's second method provides a sufficient criterion for proving such asymptotic stability in terms of the expectation of a scalar function V as shown by the following theorem. Theorem 7.1 (Khasminskii [START_REF] Khasminskii | Stochastic stability of differential equations[END_REF]). Let V (x) be a nonnegative real-valued twice continuously differentiable function with respect to x ∈ S everywhere except possibly at x = 0. Assume that exist strictly positive constants p, k 1 , k 2 and k 3 such that for all x ∈ S we have 1. AV (x) ≤ -k 3 x p , 2. k 1 x p ≤ V (x) ≤ k 2 x p . Then the equilibrium solution x t = 0 is exponentially p-stable in the sense of definition 7.1.

We note that exponential convergence on a compact set S is stronger than almost sure convergence, and thus of all other weaker notions of convergence. This follows from dominated convergence argument to interchange the order of the limit and the expectation, see e.g. [8, Lemma 4.9] for technical details.

Lemma 2 . 1 (

 21 Exponential stability of QND systems). Consider equation (1) with L a non-degenerate Hermitian matrix and denote as ρℓ the states { ρℓ := | ψℓ ψℓ | :

  0 and √ ηM(L, ρℓ ) -iκ[F, ρℓ ] = 0. Since L ρℓ = λ ℓ ρℓ , the second condition implies that κ[F, ρℓ ] = 0. By plugging this into the other steadystate constraint, we get f [F, ρℓ ] = 0 and κ 2 D(F, ρℓ ) = -κ 2 2 [F, [F, ρℓ ]] = 0. If f = κ = 0, we are in open-loop and ρℓ is not globally asymptotically stable. When κ or f are not zero, we must have [F, ρℓ ] = 0. But then we get d (Tr(ρ t ρℓ )) = 0, thus Tr(ρ t ρℓ ) = Tr(ρ 0 ρℓ ) is time invariant and we cannot have global convergence towards ρℓ .

Figure 1 :

 1 Figure 1:In blue: Ensemble average of V (ρ) = 1 -Tr(ρσ z ) over 100 trajectories (in gray) of the closed-loop system[START_REF] Wiseman | Bayesian feedback versus markovian feedback in a two-level atom[END_REF] with initial condition ρ 0 = I/2. In green: Bound obtained from theorem 5.1.

Figure 2 :

 2 Figure 2: Comparison of the ensemble average of V (ρ) = 1 -Tr(ρσ z ) over 100 trajectories when adding to (13) a decoherence term Γ 2 D(|g e|, ρ) with Γ ∈ {0, Γ/100, Γ/10}.