Keywords: Anesthesia, hypnosis, noise, saturation, observer, control

This paper deals with the maintenance phase control of general anesthesia during surgery, involving saturated input and noisy output. The objective is to maintain the patient to some given depth of hypnosis, measured by the BIS. The control law is an observer-based control, where a dead-zone observer is built in order to mitigate the presence of the output noise. The global exponential stability of the complete closed-loop system in the noise-free case is guaranteed thanks to a linear matrix inequality condition and an input-to-state property in presence of the noise is also proven. A three-steps optimization algorithm is proposed to determine the parameters of the control law and then evaluated on a patient case.

Dead-zone observer-based control for anesthesia subject to noisy BIS measurement S. Tarbouriech * I. Queinnec * G. Garcia * M. Mazerolles * *

INTRODUCTION

Although an old problem (Soltero et al. (1951)), automatic control of anesthesia remains an active research field. Many control approaches have been proposed in the literature but, yet, most of the time, a human-based control is applied during surgery in hospital. Even so, several control approaches have been evaluated in real life [START_REF] Lemos | Robust control of maintenance-phase anesthesia[END_REF] and reference herein). From a control point-of-view, many directions of research may (and have been) investigated, related to patient uncertainty, nonlinear relationship between internal state and output indicators, drug interactions, among others. Both PID-based feedback and more advanced control techniques have been suggested in that framework [START_REF] Ionescu | Robust predictive control strategy applied for propofol dosing using BIS as a controlled variable during anesthesia[END_REF], [START_REF] Nascu | Advanced model based control studies for the induction and maintenance of intravenous anaesthesia[END_REF], [START_REF] Nogueira | Nonlinear controller for bispectral index tracking: Robustness and on-lineretuning[END_REF], [START_REF] Queinnec | Switched control strategy including time optimal control and robust dynamic output feedback for anaesthesia[END_REF]). Difference of dynamics between the fast dynamics of metabolism and circulation of propofol in the central compartment (blood) and at the site effect and the slow dynamics in muscle and fat parts of the body have also been taken into account [START_REF] Queinnec | Switched control strategy including time optimal control and robust dynamic output feedback for anaesthesia[END_REF]). Moreover, as underlined by Bailey and [START_REF] Bailey | Drug dosing control in clinical pharmacology[END_REF], an important issue for research is measurement noise, EEG signals being corrupted by as much as 10% of noise. Experimental data exhibiting noisy effect on BIS (bispectral index) evolution are, for example, given in [START_REF] Beck | Modeling and control of pharmacodynamics[END_REF].

This paper intends to contribute in that direction, with a particular consideration of noise phenomena in the measurement information. The key idea of the paper is to manipulate an adaptive dead-zone mechanism to induce a high-frequency noise rejection in the observation step. Hence, exploiting recent results about this kind

Research funded in part by ANR via project HANDY, number ANR-18-CE40-0010.

of observers [START_REF] Cocetti | High-gain dead-zone observers for linear and nonlinear plants[END_REF]), an observer and associated feedback are first built for the fast subsystem, the output (represented by the BIS) being affected by noise. In a second step, considering that the observer and control gains are given, (1) the global exponential stability of the complete closed-loop system (fast and slow) in the noise-free case is guaranteed; and (2) an input-to-state (ISS) property from the noise to the state of the closedloop system is ensured. The parameters of the dead-zone adaptation are selected solving a suitable Linear Matrix Inequality (LMI), and therefore do not request any manual tuning. The effectiveness of the proposed observer is shown by means of numerical simulations in a patient-type case on the nonlinear model of the BIS.

Notation. The notation throughout the paper is standard. For a vector x or a matrix A, x and A denote the transpose of x and A, respectively. For two symmetric matrices of the same dimensions, A and B, A > B means that A -B is symmetric positive definite. For a matrix A, He(A) = A + A, λ max (A) (resp. λ min (A)) denotes the maximal eigenvalue (resp. the minimal one) of matrix A. I and 0 stand respectively for the identity and the null matrix of appropriate dimensions. For a partitioned matrix, the symbol stands for symmetric blocks.

PROBLEM STATEMENT

Patient model

From a control engineering point of view, compartmental mammillary models are interesting and sufficiently representative models for drug dynamics in the patient body [START_REF] Lemos | Robust control of maintenance-phase anesthesia[END_REF]). A three-compartment model, known as the Pharmacokinetic (PK) model, is used to describe the circulation of drugs in a patient's body [START_REF] Derendorf | Modeling of pharmacokinetic / pharmacodynamic (pk/pd) relationships: Concepts and perspectives[END_REF]). This model is then completed by a first-order dynamics (that relies the concentration of drug in the central compartment to its action at the brain level, denoted the effect site) [START_REF] Beck | Modeling and control of pharmacodynamics[END_REF]). That leads to a four-state system described as follows 1 : ẋan = Ax an + Bu an (1) with

A =    -(k 10 + k 12 + k 13 ) k 21 k 31 0 k 12 -k 21 0 0 k 13 0 -k 31 0 k e0 /v 1 0 0 -k e0    , B =    1 0 0 0   
(2) where x an = [ x an1 x an2 x an3 x an4 ] , x an1 , x an2 , x an3 are the masses in milligrams of the propofol in the different compartments, x an4 is the effect site concentration in mg/l and u an is the infusion rate in mg/min of the anesthetic. The parameters k ij ≥ 0, ∀i = j, i, j = 1, 2, 3, are the transfer rates of the drug between compartments. The parameters k 10 and v 1 represent the rate of elimination from the central compartment and the volume of the central compartment (blood), respectively. Moreover, the fact to take into account the inter-patient variability (i.e., the variability observed between different individuals) and intra-patient variability (i.e., the variability observed within one particular individual) needs to consider that the parameters k ij are uncertain [START_REF] Coppens | An evaluation of using population pharmacokinetic models to estimate pharmacodynamic parameters for propofol and bispectral index in children[END_REF]). In the literature, several models have been described to express inter-patient variability, generally distinguishing patients according to their sex, age, weight and/or size. We use the Schnider model [START_REF] Schnider | The influence of method of administration and covariates on the pharmacokinetics of propofol in adult volunteers[END_REF]) to define a simulated patient used in the numerical evaluations.

Output indicator

The depth of anesthesia indicator communally used by clinicians is the BIS (the bispectral index), allowing to quantify the level of consciousness of a patient from 0 (no cerebral activity) to typically 100 (fully awake patient). The BIS evolution is directly related to the effect site concentration x an4 , and can be described empirically by a decreasing sigmoid function [START_REF] Bailey | Drug dosing control in clinical pharmacology[END_REF]):

y BIS (x an4 ) = y BIS0 (1 - x γ an4 x γ an4 + EC γ 50 ) (3)
y BIS0 is the BIS value of an awake patient typically set to 100. EC 50 corresponds to drug concentration associated with 50% of the maximum effect and γ is a parameter modeling the degree of nonlinearity.

Error model with noise

During a surgery, the objective is to maintain the BIS around the value y BISe = 50, or at least in an interval between 40 and 60. From (3), this value corresponds to the effect site concentration equal to EC 50 . The equilibrium point of the system is then defined (see [START_REF] Zabi | New approach for the control of anesthesia based on dynamics decoupling[END_REF]):

x e1 = x 4e v 1 , x e2 = k 12 k 21 x e1 , x e3 = k 13 k 31 x e1 , x e4 = EC 50
and the associated input and output are given by u e = k 10 x e1 , y BISe = y BIS0 2

1 The time dependence is omitted to simplify the notation.

As in [START_REF] Queinnec | Switched control strategy including time optimal control and robust dynamic output feedback for anaesthesia[END_REF] and [START_REF] Haddad | Adaptive control for non-negative and compartmental dynamical systems with applications to general anesthesia[END_REF], we use the linearized function of y BIS around the equilibrium target value y BISe = 50:

y = Cx an + k bis0 (4) with C = [0 0 0 -k bis ] and k bis = γy BIS0 4EC50 , k bis0 = (2+γ)y BIS0 4
.

Le us consider the change of variables x = x an -x e , u = u an -u e and y = y BIS -y BISe with x e = [ x e1 x e2 x e3 x e4 ] and y BISe = y BIS0 /2. Consider also that the magnitude of the control input u an is constrained between 0 and u max , then, u is constrained as follows:

-u e ≤ u ≤ u max -u e (5) where u max -u e is the maximum flow rate of the drug that can be administered in practice.

Moreover, in the paper we want to take into account that the measurement is subject to noise, defined as the additional signal ν, which allows to express the variability of the BIS signal. Then, the system under consideration written as the error model is described by: ẋ

= Ax + Bsat u (y c ) y = C x + ν (6)
where y c ∈ R is the output of the controller and the saturation sat u (y c ) is defined as

sat u (y c ) = -u e if y c < -u e y c if -u e ≤ y c ≤ u max -u e u max -u e if y c > u max -u e (7)

Control problem formulation

As previously discussed in [START_REF] Queinnec | Switched control strategy including time optimal control and robust dynamic output feedback for anaesthesia[END_REF], system (6) may be revwritten as follows:

ẋf = A f x f + A f s x s + B f sat u (y c ) (8a) ẋs = A sf x f + A s x s (8b) y = C f x f + ν (8c)
with a fast dynamics state

x f = [ x1 x4 ] ∈ R n f and a slow dynamics state x s = [ x2 x3 ] ∈ R ns ,
and with matrices A f , A s , A f s , A sf , B f and C f issued from the original matrices A, B and C defined in (2) and (4). Note also that the action u directly acts on the output y through the fast subsystem, which means that the control problem focuses on the fast subsystem disturbed by the slow one.

The key idea of the control problem is to design an observer-based state-feedback controller, taking advantage of a particular kind of observer, namely a dead-zone observer such as proposed in [START_REF] Cocetti | High-gain dead-zone observers for linear and nonlinear plants[END_REF], which allows to deal with the high-frequency noise affecting the system output. Both the observer and the state-feedback controller are designed for the fast subsystem.

Then, in a first step we consider that the slow subsystem can be neglected and the dead-zone observer is build for the fast subsystem as defined by: ẋf

= A f xf + B f sat u (y c ) + Ldz √ σ (ŷ -y) ŷ = C f xf (9)
where xf ∈ R n f is the estimated state, ŷ ∈ R p is the estimated output and σ ∈ R p is a vector whose entries are non-negative and define the magnitude of the dead-zone on the corresponding output channel. L ∈ R n f ×p is the classical observer gain. The function dz √ σ is a decentralized vector-valued dead-zone from R p to R p , where

√ σ ∈ R p ≥0
is a component-wise square root. Moreover, σ is a timevarying vector following the adaptation law:

σ = -Λσ +    (ŷ -y) R 1 (ŷ -y) . . . (ŷ -y) R p (ŷ -y)    , σ ∈ R p ≥0 ( 10
)
where Λ ∈ R p×p is a diagonal positive definite matrix, and R 1 , . . . , R p ∈ R p×p are symmetric positive semi-definite matrices. The dynamics (10) ensures the invariance of the orthant R p ≥0 and then the non-negativity of σ. Therefore √ σ is well defined.

In a second step, we want to compute a state-feedback control y c = K xf , with K ∈ R m×n f , such that the fast subsystem is asymptotically stable. Therefore the observer under consideration is described by ẋf

= A f xf + B f sat u (K xf ) + Ldz √ σ (ŷ -y) ŷ = C f xf (11)
with the σ-dynamics defined in (10).

The closed-loop system issued from (8a), ( 8b), ( 8c), ( 11) can be written:

ẋf = A f x f + A f s x s + B f sat u (K xf ) ẋs = A sf x f + A s x s ẋf = A f xf + B f sat u (K xf ) +Ldz √ σ (C f xf -C f x f -ν) (12)
with (10).

By defining φ u (Kx f ) = sat u (K xf ) -K xf and considering the link between the vector-valued saturation and deadzone functions, sat √ σ (w) = w -dz √ σ (w), the closed-loop system can be written as:

ẋ = A 0 x + B 1 φ u (K 1 x) + B 2 sat √ σ (K 2 x -ν) + B 2 ν σ = -Λσ +    (K 2 x -ν) R 1 (K 2 x -ν) . . . (K 2 x -ν) R p (K 2 x -ν)    (13) with x = x f x s x f ∈ R n , n = 2n f + n s , and
A 0 = A f A f s B f K A sf A s 0 -LC f 0 A f + B f K + LC f ; B 1 = B f 0 B f B 2 = 0 0 -L ; K 1 = [ 0 0 K ] ; K 2 = [ -C f 0 C f ] (14)
Remark 1. The objective of the adaptation mechanism (10) is to manage two antagonistic effects: the adaptation speed and the filtering action tuned by Λ and R 1 , . . . , R p , respectively. Moreover, selecting Λ large enough, and R 1 , . . . , R p sufficiently small a classical Luenberger observer could be recovered. By the same way, the matrix K has to be selected such that

A f + B f K is Hurwitz but slower than A f + LC f .
The control problem we intend to solve can then be summarized as follows:

Problem 1. Find the control and observer gains K and L, the adaptation matrices Λ and R 1 , . . . , R p such that:

(1) For ν = 0, the closed-loop system is globally exponentially stable for any initial condition belonging to R n ; (2) For ν = 0, an Input to State (ISS) property from ν to (x, √ σ) is ensured.

Note that the global stability of the closed-loop system can be investigated because the open-loop matrix A and the closed-loop one A 0 are Hurwitz.

THEORETICAL CONDITIONS

In order to address Problem 1 we consider three steps:

(1) The observer gain L is designed by using Theorem 1 in [START_REF] Cocetti | High-gain dead-zone observers for linear and nonlinear plants[END_REF], recalled below;

(2) The state feedback gain is designed such that A f + B f K is Hurwitz;

(3) The complete closed-loop system (13) is analysed.

Observer gain design

Let us first provide the following results directly derived from [START_REF] Cocetti | High-gain dead-zone observers for linear and nonlinear plants[END_REF], consisting in designing the observer gain L for the observer ( 11)-(10). To this aim let us first write the dynamics of the estimation error

e f = xf -x f :          ėf = (A f + LC f )e f -Lν -Lsat √ σ (C f e f -ν) σ = -Λσ +    (C f e f -ν) R 1 (C f e f -ν) . . . (C f e f -ν) R p (C f e f -ν)    (15)
Theorem 1. If there exist a symmetric positive definite matrix

P f ∈ R n f ×n f , a matrix X ∈ R n f ×p , a symmetric semi-definite positive R ∈ R p×p , two diagonal positive definite matrices U ∈ R p×p , Λ ∈ R p×p such that He P f A f + XC f + C f RC f -X U C f -U -Λ < 0 (16) with R := p j=1
R j , then the gain L as

L = P -1 f X (17)
makes (15) globally exponentially stable to the origin for ν = 0 and ISS from ν to (e f , √ σ). Remark 2. Additional pole-placement constraints may be considered to impose some performance on the observer behavior. Typically, a pole-placement in a band {-α L2 , -α L1 } may be added, and/or a pole-placement in a disk centered in -α L3 and of radius α L4 , as additional LMI conditions [START_REF] Furuta | Pole assignment in a specified disk[END_REF]).

Control gain design

A state feedback gain K is designed such that both A f + B f K and A 0 are Hurwitz, the condition on A 0 ensuring the asymptotic stability of the closed-loop system when saturation and noise do not occur. Hence, given the gain L, if there exist symmetric positive definite matrices

W f ∈ R n f ×n f , W f s ∈ R (n f +ns)×(n f +ns) , and a matrix Y ∈ R m×n f such that W f A f + A f W f + B f Y + Y B f ≤ 0 (18) W f s 0 0 W f A 1 + A 1 W f s 0 0 W f + B 1 [0 Y ] + 0 Y B 1 ≤ 0 (19) are satisfied, with A 1 = A f A f s 0 A sf A s 0 -LC f 0 A f + LC f , then the gain K = Y W -1 f ( 20
)
makes A f + B f K and A 0 Hurwitz. Of course, as suggested in Remark 2, one should add pole-placement constraints on the dynamics of the fast closed-loop system, typically such that the poles of

A f + B f K are slower than those of A f + LC f , namely are slower than α L = | max (Re(A f +LC f ))|,
but faster than those of the open-loop A f , considering a band {-α K2 , -α K1 } and/or a pole-placement in a disk centered in -α K3 and of radius α K4 . These constraints are expressed as additional LMIs in the same decision variables.

Main conditions

With the gains L and K in hand we can provide the following result with respect to the complete closed-loop system (13) to solve Problem 1. Theorem 2. Given the gains L and K defined as in ( 17) and ( 20), respectively. If there exist a symmetric positive definite matrix P ∈ R n×n , a symmetric semi-definite positive matrix R ∈ R p×p , three diagonal positive definite matrices

S ∈ R m×m , U ∈ R p×p , Λ ∈ R p×p such that   A 0 P + P A 0 + 2K 2 RK 2 P B 2 + K 2 U P B 1 -K 1 S -2U -2Λ 0 -2S   < 0 (21) then with R = p j=1 R j ,
(1) for ν = 0 the closed-loop system ( 13) is globally exponentially stable, for any intial condition in R n ; (2) for ν = 0, the closed-loop system ( 13) is ISS from ν to (x, √ σ).

Proof: The strict negativity of the matrix in the lefthand side of relation ( 21) (by using the schur complement) implies that there exists a sufficient small positive scalar c 0 such that M ≤ -2c 0 I with

M = He(P A 0 + K 2 RK 2 ) P B 2 + K 2 U P B 1 -K 1 S -2(1 -c 0 )Λ -2U 0 -2S .
Consider the candidate Lyapunov function V (x, σ) = x P x + 21 σ, with P = P > 0 and 1 is the vector defined by 1 = [1 1 . . . 1] ∈ R p and recalling that σ ∈ R p ≥0 . The function V (x, σ) is positive definite and radially unbounded on R n ×R p ≥0 and satisfies the following bounds,

α 1 (x, √ σ) 2 ≤ V (x, σ) ≤ α 2 (x, √ σ) 2
, where α 1 := min{λ min (P ), 2}, and α 2 := max{λ max (P ), 2}.

The time-derivative of the function along the trajectories of the closed-loop system (13) reads:

V (x, σ) = 2x P ẋ + 21 σ = x (A 0 P + P A 0 )x + 2x P B 1 φ u (K 1 x) + 2x P B 2 sat √ σ (K 2 x -ν) + 2x P B 2 ν -21 Λσ + 2(K 2 x - ν) p j=1 R j (K 2 x -ν).
We use several properties for the nonlinearities φ u and sat √ σ . In particular we use Lemma 1.6 and Lemma 1.4 in [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF] with respect to φ u and sat √ σ , respectively, which lead to the following conditions:

• φ u S(φ u + Gx) ≤ 0 for any x ∈ R n and any positive diagonal matrix S ∈ R m×m . Since in the current case, we are studying the global stability, we choose

G = K 1 ; • sat √ σ (K 2 x -ν) U (K 2 x -ν -sat √ σ (K 2 x -ν)) ≥ 0, for any positive diagonal matrix U ∈ R p×p ; • 1 Λσ -sat √ σ (K 2 x -ν) Λsat √ σ (K 2 x -ν) ≥ 0 Then, recalling matrix M previously defined, one gets V (x, σ) ≤ V (x, σ) -2φ u S(φ u + K 1 x) + 2sat √ σ (K 2 x - ν) U (K 2 x -ν -sat √ σ (K 2 x -ν)) ≤ V, with V =   x sat √ σ (K 2 x -ν) φ u (K 1 x)   M   x sat √ σ (K 2 x -ν) φ u (K 1 x)   -2c 0 1 Λσ + 2x P B 2 ν -4x K 2 Rν + 2ν Rν -2sat √ σ (K 2 x -ν) U ν =   x sat √ σ (K 2 x -ν) φ u (K 1 x)   M   x sat √ σ (K 2 x -ν) φ u (K 1 x)   +2   x sat √ σ (K 2 x -ν) φ u (K 1 x)     P B 2 -2K 2 R -U 0   ν -2c 0 1 Λσ + 2ν Rν
To prove the GES when ν = 0, it suffices to put ν = 0 in the expression of V. Then, the GES is obtained if V (x, σ) ≤ V < 0, which leads to

V (x, σ) ≤ -2c 0   x sat √ σ (K 2 x -ν) φ u (K 1 x)   2 -2c 0 Λσ ≤ -2c 0 ( x 2 + sat √ σ (K 2 x -ν) 2 ) -2c 0 λ min (Λ)σ ≤ -2c 0 ( x 2 +σ) -2c 0 λ min (Λ)σ ≤ -2ρ (x, √ σ) 2 ≤ -2 ρ α2 V (x, σ)
, with ρ a positive scalar. Therefore if relation ( 21) holds then the first item of Theorem 2 is satisfied.

To prove the ISS property, we use the expression of V when ν = 0. It follows:

V

(x, σ) ≤ V ≤ -2ρ (x, √ σ) 2 +2 R ν 2 + P B 2 -2K 2 R -U x sat √ σ (K 2 x -ν) ν ≤ -2ρ (x, √ σ) 2 +2 R ν 2 +2ρ 1 (x, √ σ) ν ≤ -2ρ (x, √ σ) 2 +2 R ν 2 +ρ 1 (x, √ σ) 2 + ρ 1 ν 2
where we use the Young inequality to upper-bound the crossed terms in ν and (x, √ σ). By choosing = ρ ρ1 one gets:

V (x, σ) ≤ V ≤ -ρ (x, √ σ) 2 +(2 R + ρ 2 1 ρ ) ν 2 < 0 provided that (x, √ σ) 2 ≥ ρ -1 (2 R + ρ 2 1 ρ ) ν 2 . Therefore the closed-loop system if ISS from ν to (x, √ σ) with the gain ρ -1 (2 R + ρ 2 1 ρ ).
2

Remark 3. If the condition ( 21) is not feasible, the local case could be studied by deriving conditions inspired from [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF], and therefore by using Lemma 1.6 as mentioned in the proof but with G = K 1 being a new decision variable, and considering symmetric bounds on u in (7). In this case, a set S 0 in which the local exponential stability is ensured has to be characterized. Moreover, the local context should be addressed in the case where we would like to take into account the inter-patient variability (a wide range of adult patients uncertainty), that is the presence of uncertainty in the system and the analysis of the robustness of the closed loop.

         conditions (18), (19) pole-placement : α K1 = 1 2 α L , α K4 = 1 α K3 = 1.1 1 2 α L -α K4 Keep K = Y W -1 f . •
Step 3. Analysis of the closed-loop system (13). An observer-based control has been proposed, where a dead-zone observer is built in order to reduce the effect of the output noise. The dead-zone observer and the feedback gain are designed for the fast subsystem. Then, the global exponential stability of the complete closedloop system in the noise-free case is guaranteed thanks to a linear matrix inequality condition and an input-tostate property in presence of the noise is also proven. A three-steps optimization algorithm is proposed to select the parameters of the control law without any manual tuning, and then evaluated on a patient case.

These preliminary results pave the way for future work.

In particular, the reference tracking problem should be considered by adding an integrator in the control loop. Furthermore, in order to be closer to real-life anesthesia, sampled output or even event-triggered output controller could be investigated and delays should be considered in the control problem. Patient variability should also be considered as part of the control problem.

  4. NUMERICAL ILLUSTRATIONLet us consider a patient whose characteristics are: male, 53 years, 177 cm and 77 kg. The associated pharmacokinetic parameters are computed with the Schnider model[START_REF] Schnider | The influence of method of administration and covariates on the pharmacokinetics of propofol in adult volunteers[END_REF]), allowing to define: associated to BIS = 50 given by:x e1 = 14.51mg, x e2 = 64.26mg, x e3 = 809.2mg, x e4 = 3.4mg/l, u e = 6.08mg/min The parameters used in the output equation (3) are EC 50 = 3.4mg/l and γ = 3 (Bailey and Haddad (2005)), allowing to determine those of the linearized output equation (4): k bis = 22.06 and k bis0 = 125. To solve Problem 1, one considers the following algorithm, with pole-placement constraints as discussed in Remark 2 and in Section 3.2: • Step 1. Design of the observer Find P f , X, Λ, R and U solution to the feasibility problem issued from Theorem 1 given by:    condition (16) Λ ≤ 100, R ≤ Λ pole-placement: α L1 = 6, α L3 = 6, α L4 = 4 Keep the observer gain L = P -1 f X, and compute α L . • Step 2. Design of the controller gain K Find W f , W f s , Y solution to the feasibility problem:

  Find P , S, U , R and Λ solution to the optimization problem: max trace(R) under conditions (21) and Λ ≤ 100 • Step 4. The observer and controller gain are L and K solution to Step 1 and Step 2, respectively. Matrices of the adaptation law Λ and R are solution to Step 3. Note that as in Cocetti et al. (2019), we impose a bound on Λ (namely 100) to reduce the time scale among the observer dynamics and the adaptation dynamics. Applying the algorithm on the patient model data above given, one obtains the following solution: system, we consider an open-loop induction phase such as described in Queinnec et al. (2019) before to close the loop. It means that one injects 100 mg/min of propofol during 30 seconds followed by 1 minute at rest (u an = 0). Then at time t = 1.5 minutes, σ and xf are initialized to 0 2 and the loop is closed (vertical dashed line in the figure). Excepted for the initial induction bolus, a bound on the flow rate of propofol (between 0 and 20 mg/min) is considered. The time evolution of the noisy BIS output and control input u an are shown in Figures 1 and 2, respectively. The evolution of the deadzone bound √ σ is plotted in Figure 3.

Fig. 1 .

 1 Fig. 1. Time evolution of the BIS for the system involving the dead-zone observer

Fig. 4 .

 4 Fig. 4. Time evolution of u an for the system involving a Luenberger observer without dead-zone attenuation 5. CONCLUSION Control of general anesthesia has been addressed in this paper, considering both input saturation and output noise.An observer-based control has been proposed, where a dead-zone observer is built in order to reduce the effect of the output noise. The dead-zone observer and the feedback gain are designed for the fast subsystem. Then, the global exponential stability of the complete closedloop system in the noise-free case is guaranteed thanks to a linear matrix inequality condition and an input-tostate property in presence of the noise is also proven. A three-steps optimization algorithm is proposed to select the parameters of the control law without any manual tuning, and then evaluated on a patient case.

Following the initial bolus, a more accurate initial estimate could be possibly considered.


D Ben Ali

T Badr

T Brézillon

R Dubessy

H Perrin

A Perrin

Detailed study of a transverse field Zeeman slower

We present a thorough analysis of a Zeeman slower for sodium atoms made of permanent magnets in a Halbach configuration. Due to the orientation of the magnetic field, the polarisation of the slowing laser beam cannot be purely circular leading to optical leakages into dark states. To circumvent this effect, we propose an atomic state preparation stage able to significantly increase the performances of the Zeeman slower. After a careful theoretical analysis of the problem, we experimentally implement an optical pumping stage leading to an increase of the magneto-optical trap loading rate by 3.5. Such method is easy to set up and could be extended to other Zeeman slower architectures.

I. INTRODUCTION

Zeeman slowers are a popular technology widely used to decelerate thermal beams of atoms to velocities where they can be laser cooled and trapped. The various designs that have been proposed and implemented since its first realisation [1] all rely on the radiation pressure from a contra-propagating laser beam. In order to compensate for the Doppler frequency shift and maintain the laser beam at resonance with an electronic transition of the atom, a static and spatially varying magnetic field is used to induce an opposite Zeeman frequency shift. The shape of this magnetic field is chosen to fit the velocity change of the atoms along their trajectories.

The orientation of the magnetic field inside the apparatus defines a preferred quantization axis which in turn determines the polarisation of laser beams. In the different Zeeman slower implementations, the magnetic field is either colinear [1][2][3][4][5][6][7][8][9] or orthogonal to the atom propagation axis [9][10][11][12][13]. In the first case, pure circular polarisation of the laser beam is achievable allowing for an atom-light interaction limited to a cycling transition, within an effective two-level system. In the other case, a pure circular polarisation of the laser beam is not possible and only the fraction of the laser field which is properly polarised contributes to a cycling transition. The interaction of the atom with other polarisation components of the laser beam may lead to optical leakages into dark states. Appropriate strategies are required to avoid this effect and preserve the performances of the Zeeman slower.

Here we present a thorough analysis of a transverse magnetic field Zeeman slower based on permanent magnets in a Halbach configuration [14]. Such design has the advantage to be compact and easy to install and remove without the need for electric power nor water cooling. It also allows for very smooth magnetic field profiles with excellent transverse homogeneity and low stray magnetic fields. It has originally been introduced for 87 Rb atoms [12] and we show how to adapt it to 23 Na atoms taking into account their specific level structure. Thanks to numerical simulations based on the resolution of optical Bloch equations, we are able to precisely identify the origin of the leakage into dark states, due to the orientation of the magnetic field in this configuration. This allows us to propose an adequate response based on a proper atomic state preparation upstream from the Zeeman slower. We finally present experimental measurements demonstrating the performances of our implementation and conclude with general considerations about Zeeman slower technologies.

The organisation of the paper is as follows. In Section II we present our experimental design of a Zeeman slower in a Halbach configuration dedicated to sodium atoms. This is the opportunity to introduce a Matlab routine [15] able to optimise the design of any Zeeman slower in a Halbach configuration. In Section III we present results of numerical simulations of the atom-light interaction upstream from the Zeeman slower. Depending on the orientation of the magnetic field and on the laser beam polarisation, we show how the atoms tend to accumulate in the dark hyperfine state of the sodium ground state. We then propose a method to circumvent this issue relying on an optical pumping stage with two laser beams. As we show, an atom optically pumped into the right magnetic substate doesn't suffer from significant leakages into dark states once inside the Zeeman slower. In Section IV we present the experimental performances of our Zeeman slower design and stress how the optical pumping stage significantly enhances the flux of decelerated atoms feeding a magneto-optical trap (MOT). We finally discuss our results and give general conclusions in section V.

II. EXPERIMENTAL DESIGN

Like all alkali atoms, the sodium ground state presents an hyperfine structure [START_REF] Steck | Sodium D Line Data[END_REF] and cycling transitions are only available on the D 2 line. Our Zeeman slower design follows the same architecture than the one described in [12] adapting its shape to fit sodium atoms specifici-ties. As in [12], we work with an increasing magnetic field configuration. As a result, the laser beam addresses the |F = 2, m F = -2 to |F = 3, m F = -3 transition, where F and F represent respectively the total angular momentum of the ground and excited states and m F and m F their projections onto the quantization axis. In the rest of the paper we label |S = |F = 2, m F = -2 the ground state of this slowing transition.

A. Magnetic field profile within a simple model

Modelling the cycling transition as a two level system and using a semi-classical approach, the force induced by the radiation pressure onto the atoms writes [START_REF] Metcalf | Laser Cooling and Trapping[END_REF]:

F = -η Γ 2 hku z ( 1
)
where h is the Planck constant, k the light wave vector, Γ the natural line width of the transition, u z the unit vector along the Zeeman slower axis z and

η = s 0 1 + s 0 + 4δ 2 /Γ 2 .
(

) 2
Here δ = ω laser -ω atom is the laser beam detuning to the atomic transition, s 0 = I/I sat the resonant saturation parameter, I the laser beam intensity and I sat the saturation intensity of the transition. Assuming that the atom propagates along the z axis against the laser beam, in the presence of a magnetic field B and taking into account the Doppler frequency shift due to the atom velocity v = vu z , the detuning writes

δ = δ 0 + kv -µ B B h (g F m F -g F m F ) (3)
δ 0 + kv + µ B B h ( 4
)
where δ 0 is the detuning for an atom at rest and in the absence of magnetic field, µ B the Bohr magneton and g F and g F the Landé factors of the hyperfine states. If δ remains constant along the atom propagation, the radiation pressure force is also constant and energy conservation implies the following relation for the velocity

v(z) at position z v(z) 2 = v(0) 2 -η hk m Γ z. (5)
Combining Eq. ( 4) and ( 5) sets the model profile of the magnetic field along the atom propagation [1]:

B(z) = B(0) + ∆B 1 -1 - z L (6)
with B(0) = h(δ -δ 0 -kv(0))/µ B , ∆B = hkv(0)/µ B and L = mv(0) 2 /(ηhkΓ ).

B. Choice of the experimental parameters

According to Eq. ( 6), in order to set the magnetic field shape for experimental implementation, the values of ∆B, B(0) and L -or equivalently η -have to be defined.

The choice of ∆B determines the initial velocity v(0) which corresponds to the largest velocity of the atoms that will be decelerated in the Zeeman slower. Ideally, in order to maximize the flux of decelerated atoms, the latter has to be as large as possible compared to the peak value of the velocity distribution in the thermal effusive beam of atoms originating from the oven, which writes [START_REF] Ramsey | Molecular Beams[END_REF]

P (v) = m k B T 2 v 3 2 exp - mv 2 2k B T . (7)
In the experiment, the oven has a temperature of T = 570 K. Then, imposing v(0) = 950 m/s allows for the slowing of about 64% of the atomic flux. This leads to ∆B = 1152 G which is on the order of the highest amplitude of the magnetic field achievable with rare-earth magnets in a Halbach configuration taking into account typical experimental constraints. For a given value of the laser beam intensity I, the parameter η only depends on |δ| and is bounded from above by s 0 /(1 + s 0 ) < 1. This sets a lower limit on the value of L. To compensate for possible local imperfections of the magnetic field achieved experimentally, it is safer to design its profile with a strictly positive value for |δ|. Setting η = 0.52 leads to L = 0.96 m. This in turn imposes a minimum value for the laser intensity I tot > I min 13.6 mW/cm 2 . Here, the laser polarisation is chosen linear and orthogonal to the direction of the magnetic field, which maximises the projection onto the relevant σ -polarisation addressing the slowing transition, i.e I = I tot /2.

Once ∆B and L are set, the value of B(0) is arbitrary but must remain positive to avoid a zero magnetic field region which cannot be easily designed using permanent magnets. This in turn determines the detuning of the laser beam δ 0 since |δ| δ 0 . A value of B(0) = 300 G leads to δ 0 /(2π)

-2 GHz which allows to almost suppress the effect of the laser beam onto the magneto-optical trap where the atoms are captured after deceleration. Moreover, with such a value for B(0), electronic level crossings within the excited state of the atoms are also avoided inside the Zeeman slower [START_REF] Steck | Sodium D Line Data[END_REF], as illustrated in figure 1. Above 300 G, there is a single isolated σ -transition, which justifies the use of an effective two-level model to describe the deceleration process [START_REF] Joffe | [END_REF]. The residual leakage into dark states due to the other, far-detuned, transitions is studied quantitatively in section III.

In order to define the output velocity of the Zeeman slower, one could in principle slightly shorten the length of the magnetic field profile. Since for z L the slope of the magnetic field profile is almost vertical, see Eq. (6) is actually more efficient experimentally to slightly adjust the detuning δ 0 to optimize it.

C. Experimental implementation

The magnetic field of a Zeeman slower should ideally be transversally homogeneous across the atomic beam and should not leak outside the apparatus. Only few technological designs are actually able to produce magnetic fields respecting these two constraints. The most common one relies on current-carrying wires wound around the atomic beam. It readily produces a transversely homogeneous magnetic field pointing along the axis of the solenoid. However it often represents a large amount of metal that needs to be carefully wound in order to produce a smooth magnetic field profile. Moreover, once built it cannot be easily disassembled and limits baking possibilities for the vacuum chamber. Finally, in order to reach high magnetic field values, it requires a large current circulating in the coils [20], representing hundreds of watts of power that needs to be dissipated.

Alternatively it is possible to use permanent magnets in order to mimic a Halbach configuration [14]. Here, the orientation of the magnetization of the magnets is such that the magnetic field is homogenous within the Halbach ring and quickly falls to zero outside [12]. Fig. 2 (a) shows an example involving eight NdFeB parallelepipedic magnets [21] of size 6 × 6 × 128 mm 3 and remanent field B R = 10.8 kG. The resulting magnetic field, deduced from numerical calculations, points along the y axis and is relatively homogeneous transversely as shown in Fig. 2 (b),(c) and (d). In order to tune the amplitude of the magnetic field along the z-axis, one can slightly tilt the magnets individually with respect to the symmetry axis of the ring. Such assembly can be seen as a building block for Zeeman slowers.

Combining eight successive Halbach rings as shown in Fig. 3 (a) and (c) allows to reproduce the model profile of Eq. ( 6) taking into account the parameters introduced in section II B (see Fig. 3 (b) and(d)). Since the shape of the model profile is quite steep toward the end, an additional ring made of 10 mm-side cubic magnets of N35 grade (B R = 11.7 kG) is added at the end of the structure [12]. In order to adjust the geometry of the Zeeman slower shown in Fig. 3 (a), we have optimised by hand the values of the diameters of the different Halbach rings, their position along the z-axis and the angle of the magnets with respect to the symmetry axis of the ring. Since a design where the angles of the magnets are all different would require to be able to hold them individually, we have opted for a simpler mechanical solution close to Ref. [12] implementation where all the magnets apart from the last ring share the same angle in a guide. This leads to Fig. 3 (a) and (b) where our experimental measurements show very good agreement with the theoretical expectation.

It is worth noting that since the magnetic field produced by each magnet decreases within its length on each end (see also Fig. 2 (b)), the total length of the magnets assembly is larger than the extension of the magnetic field profile. Moreover, as the model profile becomes extremely steep toward the output of the Zeeman slower, reproducing experimentally its shape with accuracy, and in particular its maximum, is beyond reach. In our experimental realisation, the measured value of the magnetic field maximum is 100 G lower than the theoretical one. This implies that the experimental values for the capture velocity and the output velocity will be smaller than the theoretical ones. It is not an issue as they can be readily fine-tuned by changing the laser detuning. The experimental optimum for the MOT loading rate is actually δ 0 /(2π) = -1885 MHz.

Adding a magnetic shield made of a 1 mm-thick soft iron sheet around the permanent magnets assembly improves the fall off of the magnetic field outside the Zeeman slowing region (see inset Fig. 3 (b)). At the position of the MOT, the residual magnetic field drops from 1 G to 0.5 G thanks to the shielding. The complete assembly is also fairly easy to install and remove on the experimental setup and the resulting magnetic field is much smoother than what is achievable with wound wires.

Setting as constraints that the different Halbach rings must lie next to each other along the z-axis and that the diameter of the rings must decrease continuously at the exception of the last ring, and letting the angles of each magnet be free parameters, we have been able to check using a Matlab optimisation routine [15] that it is possible to almost perfectly reproduce the model profile of the magnetic field as shown in Fig. 3 (d). The routine is versatile and once the number of Halbach rings is fixed, allows for the optimisation of any parameter of the mag- nets geometry: their sizes, the diameter of the different rings and their angle with respect to the symmetry axis of the Zeeman slower. It is also possible to optimize the parameters of the model profile, for example η, which can be interesting if the length of the magnets is fixed since it controls the length of the Zeeman slower. We have made a few tests with different model profiles and the optimisation routine converges typically after a few tenths of seconds. The field resulting from the optimisation is in good agreement with the model profile. The relative error is typically smaller than 0.1% when the angle of each magnet is optimised independently (see Fig. 3 (c)).

(b) z [m] B [G] Bx By Bz (a) x [m] z [m] y [m] (c) x [m] y [m] By [G] (d) x [m]

III. NUMERICAL SIMULATIONS A. Atom trajectories

Since the experimental implementation of the Zeeman slower doesn't produce a magnetic field profile which perfectly reproduces the one deduced from Eq. ( 6), it is interesting to simulate the classical trajectories of atoms belonging to different velocity classes to check how the efficiency of the assembly is affected. From Eq. (1) we deduce the acceleration of each atom, which depends on its velocity and on the amplitude of the magnetic field at the position of the atom. We can then solve the equation of motion and infer the trajectory of the atoms in velocity space within the two-level model. Taking into account the experimental magnetic field profile shown in Fig. 3 (b) and setting the laser intensity to I tot = 87 mW/cm 2 and the detuning to δ 0 /(2π) = -1885 MHz as in the experiment, see section IV, we obtain the results presented in Fig. 4 (a). Up to a velocity slightly smaller than 900 m/s, the atoms decelerate inside the Zeeman slower to reach a velocity close to 13 m/s. Atoms starting with an initial velocity smaller than this latter value are ultimately decelerated to a zero velocity before being pushed back. Atoms with an initial velocity close to 1150 m/s start at resonance with the laser beam and are therefore slightly slowed down before entering the Zeeman slower. All the other velocity classes are unaffected. From Fig. 4 (b), we observe that the detuning δ of all the decelerated atoms increases up to a value on the order of -Γ where it locks, the force remaining constant for the rest of the trajectory. Taking into account the initial velocity distribution of the atoms at the output of the oven following Eq. ( 7), depicted in Fig. 4 (c), we obtain in Fig. 4 (d) the final velocity distribution of the atoms at the output of the Zeeman slower. About 60% of the initial atoms are hence slowed down inside the magnets assembly.

(a) x [m] z [m] y [m] (b) z [m] By [G] (c) x [m] z [m] y [m] (d

B. Optical pumping

As mentioned above it is not possible to define a purely circular polarisation for the slowing laser beam in a Zeeman slower using permanent magnets in a Halbach configuration, because the resulting magnetic field is perpendicular to the beam propagation axis. The best choice, a linear polarisation orthogonal to the magnetic field, can be written as the sum with equal weights of σ -and σ + polarisations. Hence an atom interacting with the slowing beam can be optically pumped into dark states, which cannot be slowed down. In fact, for the increasing field Zeeman slower we built, only the ground state |S can be efficiently slowed. In order to evaluate the fraction of slowed atoms, we perform a numerical calculation to compute precisely the fraction of atoms that remain in this state during their propagation from the output of the oven to the magneto-optical trap chamber. We start with the low magnetic field section, i.e. the first 50 centimetres of the ballistic flight, from the output of the oven to the input of the Zeeman slower, where significant optical pumping occurs. We then show how this depumping effect can be mitigated by introducing a transverse repumping section. Finally we focus on the dynamics inside the Zeeman slower where optical pumping is very limited thanks to the large Zeeman splitting, see Fig. 1.

We proceed as follows. We use the framework of the optical Bloch equations (OBE) to describe the interaction of the internal atomic degrees of freedom with the laser beams [START_REF] Steck | Sodium D Line Data[END_REF][START_REF] Cohen-Tannoudji | Atom-Photon Interactions[END_REF]. Because we are mainly interested in optical pumping phenomena that occur at long time scales and result in population transfer, we reduce the full OBE system to rate equations describing only the time evolution of the atomic populations. We perform this reduction by considering that the atomic coherences follow adiabatically the populations. This approximation is not appropriate to describe the short time dynamics but we checked that it leads to an accurate determination of the long term dynamics and correctly recovers the steady state. Moreover it is quite simple to implement numerically and results in a dramatic speed up of the computation besides reducing the memory requirements (576 × 576 matrices to describe the full sodium D 2 line within the OBE, 24 × 24 matrices for the rate equation model).

(a) z [m] v [m/s] (b) z [m] δ [GHz] (c) v [m/s] P (v) (d) v [m/s] P (v)
Starting with an atom with a longitudinal velocity v at the output of the oven, we compute the evolution of its internal state populations due to the interaction with the detuned Zeeman slower beam, starting from a thermal equilibrium equipartition between the ground state sublevels. Between the output of the oven and the input of the Zeeman slower itself the atom propagates for 50 cm in the presence of a low magnetic field, mostly due to the Earth magnetic field. Therefore the atom absorbs on average very few photons and we can neglect the radiation pressure force, which in turn allows to simplify the numerical treatment of the dynamics. However the depumping effects are not negligible, as shown below. We typically follow the evolution of the internal atomic states by monitoring the populations every centimetre along the atom trajectory.

Figure 5 shows the computed evolution of the population in level |S along the propagation for different initial velocities close to the output of the oven, with the same values of the detuning δ 0 and the total laser intensity I tot as those used in section III A. For velocity classes close to the resonance, a strong optical pumping towards states of the F = 1 hyperfine manifold occurs. As the state |S belongs to the F = 2 hyperfine manifold at low magnetic fields, only a small fraction of the atoms can be slowed by the Zeeman slower. For the theoretical velocity distribution of Fig. 4 (c), we compute that only 3.7 % of the atoms enter the Zeeman slower in the |S state.

In order to mitigate this optical pumping effect we propose to pump the atoms back into the |S state using a transverse beam, with two frequencies tuned to the |F = 1 → |F = 2 and to the |F = 2 → |F = 2 transitions with purely circular polarisation. Because the lasers are perpendicular to the atomic beam and can be resonant for all velocity classes, this process can be made very efficient, as can be seen in figure 6 at position z = -0.15 m. Unfortunately, due to geometrical constraints this optical pumping stage cannot be done immediately at the entrance of the Zeeman slower and the atoms propagate again in a low magnetic field for Finally we consider the evolution of the atomic populations inside the Zeeman slower, in the presence of the increasing magnetic field, under the assumption that the atomic populations follow adiabatically the magnetic field changes. More specifically we model the continuous magnetic field inside the Zeeman slower by discrete steps of constant magnetic field, with an increment of 3 G, and compute the atomic population evolution for this particular field value, over the appropriate distance (determined mainly by the local magnetic field gradient).

The output populations are then fed to the next segment. As the magnetic field increases, the matrix elements of the electric dipole between the different atomic states are modified as the system progressively enters the Paschen-Back regime and the two circular components of the linearly polarised Zeeman slower laser beam become more and more asymmetrically detuned with respect to the allowed transitions, see figure 1. At sufficiently high magnetic field these two effects ensure that the slowing transition is almost closed such that the atoms are slowed, provided they reach the 'slowing area' in the |S internal state. Here we numerically investigate to which extent this affirmation is true and in particular we quantify the effect of the intermediate field area where level crossings could favour optical pumping into non slowed states and reduce the efficiency of the Zeeman slower.

Figure 6 shows the computed evolution of the population in level |S along the propagation for different initial velocities, from the oven to the magneto-optical trap chamber, including the repumping section using a transverse laser beam at position z = -0.15 cm.

After an initial depumping section, as in figure 5, the atoms are repumped into the level |S by the transverse beams, resulting in 31 % of the atoms entering the Zeeman slower in state |S . The populations subsequently evolve in the presence of the detuned Zeeman laser beam and the increasing magnetic field, resulting in a more complex behaviour. First, above the Zeeman slower capture velocity v c 900 m/s, the population in level |S is depleted inside the Zeeman slower. As the Zeeman slower is not designed to slow these atoms this is not an issue. Then, for velocities smaller than v c , the depumping mechanism becomes less efficient as the field increases and about 28 % of the atoms end in the slowable state at the position where the Zeeman laser becomes resonant (indicated by the vertically hatched pattern). Therefore the additional depumping due to this section is negligible in our setup. For the lowest velocities the depumping is significant inside the Zeeman slower, even at high magnetic fields, due to the long interaction times.

Our model predicts an increase of the slow atom flux at the output of the Zeeman slower, with about four times more atoms when the repumping section is present. This figure is computed by comparing the number of atoms in the |S state for each velocity class at the position inside the Zeeman slower where it is resonant with the slowing laser, convoluted by the velocity distribution, respectively with and without the transverse repumping section described above.

IV. EXPERIMENTAL PERFORMANCES

In order to set up an atomic state preparation on the experimental apparatus, we have installed a laser beam propagating horizontally along the horizontal x-axis orthogonal to the atomic beam of sodium, roughly 15 cm before the entrance of the Zeeman slower. It corresponds to the closest position where a window is available on the set-up to let the beam pass through the vacuum chamber. The laser is circularly polarised and its frequency is set close to the |F = 2, m F = -1 → |F = 2, m F = -2 transition. The total power in the beam is 25 mW for a 1/e 2 -diameter of 8 mm. We use an electro-optical modulator to produce sidebands at 1.771 GHz from the carrier frequency. They both represent about 20% of the total laser power. In order to split the Zeeman sublevels of the atoms in this region where only the Earth magnetic field is present, we apply a homogeneous magnetic field pointing in the direction of the propagation axis of the laser and produced by two 20-turn coils in a Helmholtz configuration. Its amplitude is on the order of 1 G. While the carrier optically pumps the |F = 2 atoms into the |F = 2, m F = -2 state, the blue sideband is used to repump the atoms from the |F = 1 to the |F = 2 states.

The slowing beam itself possesses a waist of 0.5 cm for an intensity of I tot = 87 mW/cm 2 . Its detuning is δ 0 /(2π) = -1885 MHz and it is linearly polarised along the x axis which means that only half of the total power of the beam contributes to the slowing transition.

Loading a magneto-optical trap (MOT) from the slowed atomic beam for 3 s allows to test the performances of the atomic state preparation. The dependence of the captured atom number with the laser detuning is shown in figure 7[23]. The pump laser intensity has been optimized here such that increasing it further does not change the final populations anymore. When the pumping laser is resonant with the |F = 2, m F = -2 → |F = 3, m F = -3 transition, the atoms end up cycling over the transition, exchanging a large number of photons with the laser which increases significantly the divergence and the deflection of the atomic beam, degrading in turn the performances of the MOT loading. This is no more the case when the laser is resonant with the |F = 2, m F = -1 → |F = 2, m F = -2 transition since the |F = 2, m F = -2 state is in this case a dark state and all the atoms get optically pumped into it. It is worth noting that since the magnetic field direction is different in the optical pumping area (along the x axis) compare to the slowing region inside the Zeeman slower (along the y axis), we assume the spin of the atoms to adiabatically follow the local direction of the magnetic field along their propagation.

In order to compare the efficiency of the atomic state preparation with the simulations of section III B, we have measured the dependence of the MOT atom number with its loading time. Comparing the slope of the atom number with and without the atomic state preparation, we measure a factor of 3.5 increase in the atomic loading. This result agrees well with the factor of four deduced from the simulations. We attribute the small discrepancy to the probable over-estimation of the efficiency of the repumping section in the computation as we did not include the spatially varying intensity of the laser in the numerics and modelled the real gaussian shaped laser beams by top-hat beam profiles. We have also measured the peak velocity of the slowed beam by monitoring the fluorescence signal S induced by a probe beam crossing the atomic beam at ±45 • [1]. This signal depends on the velocity of the atoms through Doppler effect. For a single atom with longitudinal velocity v, the probe beam detuning δ p writes

|F = 2, mF = -1 → |F ′ = 2, m F ′ = -2 |F = 2, mF = -2 → |F ′ = 3, m F ′ = -3 Detuning [MHz]
δ p = δ p,0 -k p • v = δ p,0 ± k p v √ 2 (8)
with δ p,0 its detuning for an atom at rest and k p the modulus of its wavevector. We present in Fig. 8 (a) the dependence of S with δ p,0 and the detuning δ 0 of the slowing laser beam. For a given δ 0 , the observed peak in the fluorescence signal S corresponds to the peak velocity v peak of the velocity distribution of the slowed atoms : δ p,0 = ∓k p v peak / √ 2. The observed velocity width is limited by the natural width of the transition. As shown in Fig. 8 (b), the peak velocity of the slowed atomic beam depends linearly with δ 0 as expected from Eq. ( 4) if we assume that the velocity of the atoms remains constant after exiting the Zeeman slower. This is a fair assumption from the results of the simulations shown in Fig. 4 (a). For a detuning of -1885 MHz we find v peak = 28(3) m/s, slightly higher than the estimation of section III A. This discrepancy can be explained by a 1 % uncertainty on the magnetic field measurement near the end of the Zeeman slower, where the magnetic field is high and varies rapidly, see figure 3.

V. DISCUSSION

We now turn to more general considerations on the available methods to mitigate the optical pumping fects in the transverse field Zeeman slower considered in this paper. This discussion is enabled by the accurate modelling of the optical pumping dynamics that we have developed. We will focus on the increasing field Zeeman slower for sodium described in the present work, for which we have shown that a transverse repumping strategy is efficient, resulting in up to 28 % of the total atomic flux being slowed.

Let us first consider the case of a perfectly polarised slowing light, which could be achieved by a more complex arrangement of permanent magnets resulting in a longitudinal magnetic field (see for example Ref. [3]). We compute that 32 % of the flux would be slowed, due to optical depumping occurring during the initial population redistribution in the |F = 2 manifold, before the entrance of the Zeeman slower. Again, a transverse repumping section is needed to reach the theoretical limit of 60 % of the flux set by the estimated velocity capture.

The authors of Ref. [12] have shown that another solution was to use a longitudinal repumping laser, copropagating with the slowing light, with a frequency either fixed or rapidly modulated over the whole Doppler width. We only computed the former case as the latter would be much more involved. Briefly we extended our OBE description by adding the contribution of the second laser and instead of performing a rotating wave approximation we performed a Floquet expansion truncated at second order. This allows to consider single photon absorption as well as two photon transitions but neglects higher order phenomena. We find that for a repumping laser tuned to the |F = 1 → |F = 2 transition (at low magnetic field) about 19 % of the total flux is slowed.

From our computations it seems that it is particularly efficient inside the Zeeman slower, where it pumps back to the |S state a set of velocity classes locally resonant at a given magnetic field. We may expect that the result could be even better with a modulated repumping frequency. However this strategy is very demanding in terms of laser sources and probably necessitates an independent laser system for the repumping beam, which can be very costly, especially for sodium atoms.

To conclude, we have presented an efficient optical pumping method able to significantly increase the performance of a transverse field Zeeman slower made of permanent magnets in a Halbach configuration. We have experimentally demonstrated an increase by a factor of 3.5 of the atomic loading rate in the MOT. This figure could be even improved if the distance between the optical pumping stage and the entrance of the Zeeman slower could be further reduced. We believe that this technique could be extended to many different Zeeman slower architectures in order to optimise their efficiencies and in turn reduce the MOT loading time or limit the total atomic flux to improve the background vacuum pressure.

FIG. 2 .

 2 FIG. 2. (a) Three-dimensional representation of a Halbach ring made of 8 parallelepipedic magnets of size 6 × 6 × 128 mm and remanent field BR = 10.8 kG. The mean diameter of the ring is 5.4 cm and its angle with respect to the symmetry axis of the assembly is 0.97 • (b) Amplitude of the three components of the magnetic field B along the symmetry axis of the Halbach ring. The results are based on the exact formula for the magnetic field induced by a cuboid magnet presented in the supplementary information of Ref. [12]. The vertical dotted lines indicate the length of one magnet of the assembly. (c) Amplitude of By in the x -y plane for z = 0. The orientation of the magnetization of each magnet of the Halbach ring is indicated with a white arrow.(d) Zoom on the central part of (c).

FIG. 3 .

 3 FIG. 3. (a) Three-dimensional representation of the experimental implementation of our Zeeman slower. It is fabricated with consecutive Halbach rings made of NdFeB permanent magnets. The first eight rings are built with parallelepipedic 30SH grade magnets of size 6 × 6 × 128 mm 3 . The diameter of the first ring is 80 mm and the angle of the magnets with respect to the symmetry axis of the assembly is 1.25 • . The last ring is made of 10 mm-side cubic magnets of N35 grade. Its diameter is 61 mm. (b) Longitudinal profile of the y-component of the magnetic field deduced from Eq. (6) (red dashed line), obtained from exact calculations taking into account the geometry of the magnets shown in (a) (black line) or from experimental measurements with (blue crosses) or without (red circles) a magnetic shield surrounding the magnets assembly. The inset shows in log scale a zoom of the last centimetres of the magnetic field profile. (c) Three-dimensional representation of a Zeeman slower where the magnets in the different rings can have different angles with respect to the symmetry axis of the assembly. (d) Comparison of the model profile of the magnetic field (red dashed line) and the one obtained from exact calculations of the field produced by magnets assembly shown in (c) (black solid line).

FIG. 4 .

 4 FIG. 4. (a) Simulation of the trajectories in velocity space for atoms with different initial velocities inside the actual Zeeman slower implemented in the experiment (see Fig. 3 (b) for the magnetic field profile). (b) Behaviour of the detuning δ for the atomic trajectories shown in (a). (c) Initial velocity distribution at the output of the oven, following Eq. (7). (d) Velocity distribution at the output of the oven (blue bars) and at the output of the Zeeman slower (red bars). The bin size for histograms (c) and (d) is 25 m/s.

FIG. 5 .FIG. 6 .

 56 FIG. 5. Population in state |S for different initial velocities as a function of the distance travelled by the atoms from the output of the oven at z = -0.5 m to the entrance of the Zeeman slower at z = 0 m. Velocity dependent depumping along the propagation is due to the presence of the counterpropagating linearly polarised Zeeman slower laser beam, see text for detail.

8 FIG. 7 .

 87 FIG. 7. Dependence with the optical pumping laser detuning of the number of atoms captured in the magneto-optical trap after 3 s of loading with the atomic state preparation (black line, see text) or without (red line). The detuning is zero when the laser is resonant with the |F = 2, mF = -2 → |F = 3, m F = -3 transition. The error bars come from the standard deviation computed over two realisations of the experiment.

FIG. 8 .

 8 FIG. 8. (a) Fluorescence signal S in arbitrary units as a fonction of δp,0 and δ0 induced by a probe beam crossing the atomic beam at -45 • . (b) Dependence of v peak as a fonction of δ0 for a probe beam crossing the atomic beam at respectively -45 • (red crosses), 45 • (blue plus signs). The red line is a linear fit of the results. The shaded area around the line represents the uncertainty of the fit.

ACKNOWLEDGMENTS

This work has been supported by the Region Ile-de-France in the framework of DIM IFRAF (Institut Francilien de Recherche sur les Atomes Froids) and by the ANR Project No. 11-PDOC-021-01. LPL is UMR 7538 of CNRS and Paris 13 University.