Chattering Phenomenon in Quantum Optimal Control

We present a quantum optimal control problem which exhibits a chattering phenomenon. This is the first instance of such a process in quantum control. Using the Pontryagin Maximum Principle and a general procedure due to V. F. Borisov and M. I. Zelikin, we characterize the local optimal synthesis, which is then globalized by a suitable numerical algorithm. We illustrate the importance of detecting chattering phenomena because of their impact on the efficiency of numerical optimization procedures.


Introduction
Consider the experiment in which a ball bounces up and down on the ground. We assume that the impact with the ground is inelastic and that the ball is only subject to the gravity. In the ideal case in which the ball changes its speed instantaneously at each bounce, an infinite number of bounces is performed in the finite time of the process. Chattering thus refers to an observable (here the speed) having very fast discontinuities, which lead in the mathematical limit to an infinite number of jumps in a finite-time interval [1]. This type of process can also be found in quantum physics. Examples are the quantum Zeno effect and dynamical decoupling in which a repeated observation of the system and a periodic sequence of instantaneous pulses prevent, respectively, its time evolution [2] or its coupling with the environment [3]. The possibility of chattering was also established in Optimal Control Theory (OCT). OCT was founded in the sixties by Pontryagin and his co-workers [4], with a rigorous framework to design control protocols for driving a dynamical system from a given initial state into a desired target state, while minimizing energy or other resources. Chattering was found in this field by A.T. Fuller in a planar system [5,6,7]. It consists of Figure 1: Schematic description of the time scale invariance of the optimal control in the chattering process. The optimal solution of the quantum control problem described in this paper is plotted in the top panel. Near the final time t f , the control switches infinitely many times with an asymptotic invariant structure by time dilation as represented on the two lower panels (successive zooms around t = t f given by the different plot colors). a control that switches infinitely many times over a finite time period [1]. This observation runs counter to common experience for which control is viewed as a piecewise continuous function, while for chattering, the control lies in a larger class of functions [8]. In Fuller's example, the number of switchings accumulates with a geometric progression at the final control time. The control law has a time scale invariance near the final time as schematically represented in Fig. 1. At first Fuller's example was considered a curiosity, but this type of phenomenon is very widespread in optimal control, as rigorously shown few years later by I. Kupka [9]. While optimal trajectories exhibiting chattering have been found in relevant examples from medicine [10] to classical [1] and space dynamics [11], this phenomenon has, to the best of our knowledge, not yet been studied in quantum control [12,13,14,15,16]. The existence, the role and the ubiquity of this process in quantum systems remain an open question. Such control schemes are interesting from a fundamental point of view even if they turn out to not be feasible in experiments. They may also have a rather severe impact on the numerical search of optimal solutions [11,17]. It is therefore important to understand why chattering is occurring and how it can be avoided [18]. We propose in this paper to describe this phenomenon in quantum control by studying a simple but fundamental quantum system. We introduce on this key example a systematic procedure to design the optimal control protocol. The impact on the efficiency of optimization procedures is also described.

The model system
Let us consider the control of a three-level quantum system described by a pure state belonging to a Hilbert space spanned by the states |1 , |2 , and |3 . As in a standard STIRAP process [19], in a suitable rotating frame and in the rotating wave approximation, the dynamics of the system are controlled by two pulses ∆ and u that couple, respectively, states |1 and |2 and states |2 and |3 . We assume that ∆ is a constant and that the only control parameter is u = u(t). The resulting dynamics areẊ where X is a vector of real coordinates (x 1 , x 2 , x 3 ) with the condition x 2 1 + x 2 2 + x 2 3 = 1 (see A). The two skew-symmetric matrices Ω 1 and Ω 3 are given by generating respectively the rotations along the x 1 -and x 3 -directions. In the following, the control cannot exceed a certain physical bound, that is, |u(t)| ≤ u 0 for some u 0 > 0. A time rescaling results in the multiplication of the two parameters ∆ and u 0 of the problem by a positive scalar, which leads to the normalization u 0 = 1. Starting from any state X 0 on the unit sphere, the goal of the control is to steer the system to the state (0, 0, 1), still denoted by |3 , while minimizing the population transfer to the state |1 . This control protocol is interesting in practice if, e.g., the state |1 is the only state of the system subject to an unwanted relaxation process. Notice that this latter is not modeled by Eq. (1). The control protocol can be formalized as an optimal control problem by introducing the cost functional C = t f 0 x 2 1 (t)dt to be minimized, where t f is the control time, which is not fixed. The existence of a solution in finite time is not obvious and it is a consequence of the general theory developed in [1], as explained in D. We give below an argument showing that, once its existence has been established, such a control has a chattering behaviour.

Application of the Pontryagin Maximum Principle
The existence of chattering can be proved from the Pontryagin Maximum Principle (PMP), which gives a first-order necessary condition for optimality [8]. It can be stated by introducing the Pontryagin Hamiltonian where X is a point on the unit sphere, P the adjoint state of coordinates (p 1 , p 2 , p 3 ), and p 0 a constant equal either to 0 or to − 1 2 . If X(t) is an optimal trajectory with corresponding control u(t), then X(t) is extremal, namely, there exist P(t) and p 0 such thatẊ = ∂H P ∂P ,Ṗ = − ∂H P ∂X , and H P (X(t), P(t), u(t), p 0 ) = max v∈[−1,1] H P (X(t), P(t), v, p 0 ) = 0. Moreover, (P(t), p 0 ) = (0, 0) for every t, whereP(t) = P(t) − (P(t) · X(t))X(t). If p 0 = 0 (resp. p 0 = − 1 2 ) the extremal is called abnormal (resp. normal ). The condition (P(t), p 0 ) = (0, 0) (instead of the weaker but more classical one (P(t), p 0 ) = (0, 0)) reflects the fact that we already know that X(t) belong to the two-dimensional unit sphere. Actually, (X(t),P(t)) represents an element of the cotangent bundle of the sphere.
The maximization condition of the PMP is solved with the switching function Φ(t) = P(t) · Ω 1 X(t). If Φ(t) is different from zero then H P is maximal when the control, called bang, is a constant control of maximum amplitude, u(t) = sign[Φ(t)]. If Φ(t) = 0 on a time-interval (a, b), we say that the restriction of the extremal to (a, b) is a bang arc. When Φ is constantly equal to zero on (a, b) then the restriction of the extremal to (a, b) is said to be a singular arc. Singular arcs can be characterized as follows. A simple computation shows thaṫ Φ = ∆ P · Ω 2 X, where Ω 2 is the generator of the rotations along the x 2 -axis. If both Φ andΦ vanish at time t then eitherP(t) = 0 or Ω 1 X(t) and Ω 2 X(t) are linearly dependent, that is X(t) lies on the equator x 3 = 0. IfP is equal to zero on a time interval then p 0 = − 1 2 and, moreover, 0 = H P = − x 2 1 2 , leading to x 1 = 0 on the same interval. Sinceẋ 1 = −∆x 2 , it follows that X stays at |3 .
We now claim that it is enough to consider normal extremals, i.e., that any optimal trajectory corresponding to an abnormal extremal also corresponds to a normal one. Indeed, if an optimal trajectory reaches |3 in a finite time is optimal as well. Moreover, Y corresponds to a vanishing control for all times in [T,t f ]. Let us denote by Q the adjoint state corresponding to Y, whose existence is guaranteed by the PMP. Then the extremal Y is singular on [T,t f ] andQ| [T,t f ] ≡ 0, whereQ is defined in analogy toP. Such an extremal cannot be abnormal (otherwise (Q(t f ), p 0 ) would be zero). We have shown that for any optimal solution reaching |3 in time T it can be assumed without loss of generality that the extremal is normal andP| [T,t f ] ≡ 0 (this corresponds to replacing P by Q| [0,t f ] and denoting the latter by P).
The last step consists in showing that an optimal trajectory cannot be the finite concatenation of bang arcs in [0, T ], where we recall that T is the time at which X reaches |3 . This fact can be deduced by computing the so-called order of the singular trajectory and using general properties of optimal control theory (see, e.g., [10,20,21]). In our specific case we can exhibit an elementary and explicit proof. By contradiction, assume that u is constantly equal to +1 or to −1 in an interval [T − ε, T ] for some ε > 0. Then Φ is smooth on [T − ε, T ] and an explicit computation based on the dynamicṡ and on the relationP . Then Φ(t) has opposite sign with respect to u in a small left neighborhood of T , contradicting the fact that Hence the only option to reach |3 is a chattering process in which the control switches infinitely many times in a finite time-interval. Notice that, as a byproduct, we obtain that there are no optimal abnormal extremals reaching the target. Indeed, a simple computation shows that the control corresponding to an abnormal extremal reaching the target at time T is necessarily bang on an interval [T − ε, T ] for ε > 0 small enough (since Φ(T ) must be zero andΦ(T ) cannot be zero at the same time).
Having established the chattering behaviour, the next goal is to find the position of the switching points, which form the switching curve. This is not an easy task and an exact analytic expression cannot be derived. Notice however that close to the state |3 the system can be described by the two coordinates x 1 and x 2 with the dynamicṡ Taking only the dominant terms, the system can be locally approximated aṡ which is different from the standard linear approximation. In control literature, system (4) is referred to as the nilpotent approximation of (3) [22]. The dynamics of (4) together with the cost dt (t f free) and the origin x 1 = x 2 = 0 as target state, yield the classical Fuller problem, up to the change of coordinates The chattering trajectories of the Fuller model can be described exactly [7] as recalled in B. We then deduce the following properties for system (4). The optimal solution is bang-bang with an infinite number of switchings near the origin. The switching curve is defined by 0.44623. The switching times t k are given by a geometric 4.1301. The relation between the local optimal syntheses of two systems differing only by high-order terms (such as (3) and (4)), when one of the two syntheses exhibits chattering, is a subtle question. To this purpose we apply the results of [1], which permit to conclude that a system has a Fuller-like optimal synthesis provided it differs from the Fuller model by terms that are small while applying suitable non-isotropic dilations. On the basis of this study described in D, we state the following result. Proposition 1. For every point X 0 sufficiently close to |3 , any optimal solution of (1) connecting X 0 to |3 corresponds to a control having infinitely many discontinuities accumulating at the first time T at which |3 is reached (and possibly staying at |3 for larger times). The optimal synthesis is characterized by a switching curve Γ passing through |3 , whose expression in coordinates (x 1 , x 2 ) is of the form where λ 0 and λ 1 are C 1 function satisfying λ 0 (0) = −λ 1 (0) = ∆ξ. In the coordinates (x 1 , x 2 ) and locally near (0, 0) the optimal control is −1 below the curve Γ and +1 above it.

Numerical results
The optimal synthesis can be computed numerically starting from that of the Fuller model. When we are sufficiently close to |3 , we approximate the switching curve of the quantum system by that of its approximation (4). We consider a The switching curve is plotted in red. The parameter ∆ is set to 10. point of the latter curve of coordinates (ξ∆x 2 20 , x 20 ), with x 20 > 0 (notice that the same method could be used for x 20 < 0). The third component is obtained from 10 − x 2 20 and the adjoint state from the condition P·Ω 1 X = 0 together with H P = 0 [8]. The dynamics of the PMP allows to propagate X and P backward in time, uniquely up to an irrelevant component of P along X. Starting from an initial point with x 20 > 0, one integrates the equations taking u = −1. When the corresponding switching function vanishes, the control switches from −1 to +1. Then one goes on by integrating the equations with u = +1 up to the next switching time and so on. The stability of this approach is justified in E. An optimal control law can be obtained for each value of x 20 . Even if the result applied above to characterize the optimal synthesis (see Theorem 1 in D) can be used only for initial points X 0 close to |3 , numerical simulations show that optimal trajectories have the same structure everywhere in the north hemisphere.
The trajectory starting from a given initial state (say (0, 1, 0)) can then be determined by a Newton algorithm to estimate the right parameter x 20 from which the backward propagation arrives at the initial state. The parameter x 20 is not unique because the forward optimal trajectory has infinitely many switching points near the target. In practice, the choice of x 20 is dictated by the required precision on the final state. For ∆ = 10 and a precision of 10 −3 , numerical simulations lead to x 20 = 6.9 × 10 −4 . Figure 2 depicts the optimal trajectory on the unit sphere. The corresponding control u(t) and switching function Φ(t) are displayed in Fig. 3.
The switching curves are reconstructed numerically by varying the small parameter x 20 and collecting all the switching points of the corresponding trajectories obtained by integrating the extremal flow backwards in time. The result is represented in Fig. 2 and 4. We observe in Fig. 4 that the points (±1, 0, 0) belong to the switching curves of the quantum system, independently The geometric analysis gives the optimal control protocol with a very high numerical precision. Due to their complexity, quantum control problems are usually solved by numerical optimization algorithms in which the control is assumed to be a piecewise constant function [12]. An intriguing question is therefore to study to what extent such solutions can approximate the chattering phenomenon. The numerical simulations presented below use a direct optimal method with the software BOCOP [25] with a fixed control time computed using our analysis. In Fig. 5, we observe that the chattering process can only be reproduced approximately by numerical optimization. Additional results can be found in E. The fineness of the time discretization corresponding here to t f /N , where N is the number of time steps, is a key factor to improve the protocol efficiency and to reproduce the control shape. However, for N = 400, we observe an erratic structure of the control. Without a precise understanding of the optimal control strategy, this switching accumulation could be misinterpreted as numerical instabilities or artifacts, while it is due to the very structure of the optimal control. Reasonable efficiencies of the control protocols are achieved for quite small values of N . Such sub-optimal strategies could be used to bypass the problem due to chattering in the numerical optimizations.

Conclusion
In summary, the analysis of a simple system has shown that chattering can appear in quantum optimal control. Such a process may play a key role in numerical optimization algorithms, as illustrated in the example by the rapid oscillations that occur in the numerical solutions. In accordance with the existing results on the ubiquity of such phenomena in OCT [9], we expect chattering to appear in many further examples, especially for higher-dimensional problems. Actually, one can give sufficient conditions on the relations between the commutators of the uncontrolled and controlled Hamiltonians to guarantee the existence of chattering solutions (see [9] and [1,Chapter 4]). Such conditions render the chattering phenomenon more and more frequent as the dimension grows, and our example shows that no obstacle to their appearance comes from the quantum structure of the control problem. It should be noticed, however, that one cannot deduce from the general conditions in high dimensions the optimality of the chattering trajectories, unlike the two-dimensional system considered in this paper.  : Schematic representation of the quantum system with the coupling ∆ and u(t) between the states |1 and |2 and |2 and |3 (red arrows). The black and blue dots indicate respectively the initial and the target states. The black arrow represents the relaxation process.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.

A Derivation of the model system
We show in this section how to derive the system used in this study. We consider a three-level quantum system in a Λ configuration whose dynamics are governed by the Schrödinger equation. The system is described by a pure state |ψ(t) which belongs to a three-dimensional Hilbert space spanned by the basis {|1 , |2 , |3 }. The system is subject to a pump and a Stokes pulses coupling, respectively, states |1 and |2 and states |2 and |3 . They are assumed to be on-resonance with the corresponding frequency transitions. There is no direct coupling between levels |1 and |3 . A schematic representation of the control problem is given in Fig. 6. The time evolution of the system is given by the Schrödinger equation where units such that = 1 are used. In the interaction representation and in the rotating-wave approximation, the Hamiltonian of the system can be written as where ∆ and u(t) are half the Rabi frequencies of the two pulses. We denote by c 1 , c 2 , and c 3 the complex coefficients of the state |ψ(t) , and we introduce the real coefficients x 1 , . . . , x 6 defined by Since H(t) is Hermitian, we have that x 2 i = 1 in accordance to the description of a pure state. Straightforward computations from the Schrödinger equation show that the variables x 1 , x 2 , and x 3 are decoupled from x 4 , x 5 , and x 6 . For our purposes it is sufficient to study the dynamics of the first set of variables, which turn out to bė In the control problem, we consider an initial state of the dynamics such that x 4 = x 5 = x 6 = 0. We deduce that x 2 1 + x 2 2 + x 2 3 = 1 at any time t.

B The Fuller model
This paragraph briefly describes the main results that can be established for the classical Fuller model. The interested reader will find the proofs of the different statements in textbooks of mathematical control theory [1,7]. The Fuller model is an optimal control problem in R 2 . The dynamics of the state (x, y) are governed by the differential equationṡ The optimal trajectories have the discrete symmetry (x(t), y(t), u(t)) → (−x(t), −y(t), −u(t)), and a scaling symmetry defined by the family of transformations where λ is a positive parameter. This means that if (x(t), y(t), u(t)) is an optimal solution then (x λ (t), y λ (t), u λ (t)) is also a solution, with initial state (λ 2 x 0 , λy 0 ) and cost J λ = λ 5 J . We deduce that if (x,ȳ) is a switching point for an optimal trajectory, i.e., the corresponding control goes from ±1 to ∓1 when the trajectory crosses (x,ȳ), then the optimal synthesis has a switching curve of equation x = −ξsign[y]y 2 , where ξ is a positive constant such that the curve passes through (x,ȳ). The second step of the analysis consists in applying the Pontryagin Maximum Principle, which is a necessary condition for optimality [4]. The Pontryagin Hamiltonian can be expressed as where (p x , p y ) is the adjoint state and p 0 is a constant multiplier equal either to 0 or to − 1 2 . If (x(t), y(t)) is an optimal trajectory with corresponding control u(t) then there exist (p x (t), p y (t)) and p 0 such that (p x (t), p y (t), p 0 ) = (0, 0, 0), x = ∂H P ∂px ,ẏ = ∂H P ∂py ,ṗ x = − ∂H P ∂x ,ṗ y = − ∂H P ∂y , and The switching function Φ is given by Φ(t) = p y (t). As in the case studied in Section 3, one can restrict the analysis to normal extremals. Integrating system (5) from the state (−ξy 2 0 , y 0 ) with y 0 > 0 and u = +1, and imposing that the corresponding switching function vanishes at a point (ξy 2 1 , y 1 ) with y 1 > 0, we obtain that ξ is solution of a polynomial equation of order 4, given by The optimal trajectory has the following symmetries: denoting by t k the k- Finally, starting from the state (−ξy 2 0 , y 0 ), y 0 > 0, it can be shown that t Ful f = 1+α α−1 y 0 . An example of optimal trajectory is plotted in Fig. 7. Figure 8 compares the results obtained for the quantum system and its approximation (Eq. (4) in the main text), which coincides with the Fuller model up to a linear rescaling of the first coordinate. In Fig. 8, we project the optimal trajectory of the quantum system onto the (x 1 , x 2 )-plane. This comparison highlights that the two solutions are very close to each other near the target state.
C Properties of the switching function for the three-level quantum system We focus in this paragraph on the switching function Φ for the three-level quantum system. Using the Hamiltonian equations for the adjoint state in the normal   case (Eq. (2) in the main text), one computes that Φ and its time derivatives can be expressed as Moreover, on a segment where Φ = 0 and, therefore, u is constantly equal to +1 or −1, one has Note that, if x 2 = 0 and Φ = 0, then p 3 = x3 x2 p 2 and hencė Notice also that, since H P = 0 and using again that Φ = 0, we have Hence, for x 2 x 3 = 0 and x 1 = 0, the derivativeΦ has the same sign as x 2 x 3 , which means that switches only occur from −1 to +1 if x 2 x 3 > 0 and from +1 to −1 if x 2 x 3 < 0. If x 2 = 0 and x 1 = 0 thenΦ = 0, which implies that p 1 = 0. ThereforeΦ = 0 and It follows that everywhere on the set {x | x 2 x 3 > 0} switches only occur from −1 to +1, while on {x | x 2 x 3 < 0} they only occur from +1 to −1. This result proves a statement of the main text, i.e. if the switching curve goes out of the north hemisphere then it passes through the points (±1, 0, 0). This property does not depend on the parameter ∆.

D A sufficient condition for chattering
In this section we present an adaptation of the results presented in Chapter 3 of the book [1] by Zelikin and Borisov, concerning sufficient conditions for the appearance of Fuller-like chattering in a two-dimensional optimal synthesis. We consider a control system of the form We also suppose that φ x i , φ y i are smooth and small in the following sense: denoting by g κ the anisotropic dilatation g κ (x, y) = (κ 2 x, κy), one has lim sup for every (x, y) ∈ R 2 . Then, the optimal control problem satisfies the following properties.
Theorem 1. For every (x 0 , y 0 ) in a sufficiently small neighbourhood of (0, 0), Problem (9) admits a solution. Such a solution reaches the fixed point (0, 0) in finite time and its corresponding control has infinitely many discontinuities accumulating at the final time. Moreover, there is no other solution of (9) up to prolongation by a constant trajectory at (0, 0). The optimal synthesis is characterized by a switching curve of the form where λ 0 and λ 1 are C 1 function satisfying λ 0 (0) = −λ 1 (0) = ξ∆, where ξ is as in (6). The optimal control is −1 above Γ and +1 below it.
All the proof material is in the book [1] by M. I. Zelikin and V. F. Borisov. Nevertheless, since this theorem is not explicitly stated as such, we propose to retrace the ingredients of the proof and point out the relevant statements in [1].
Proof. First, let us highlight the fact that the existence of a solution of (9) is not obvious. Indeed, it is a priori possible that the cost decreases by reaching the target later in time. As a consequence we are lacking compacteness for existence of the solution of the problem with free final time. Existence of optimal trajectories can be obtained once an extremal synthesis has been constructed, using a field-of-extremals argument ([1, Theorem 3.3]).
E Numerical optimization procedure E.1 Backward stability of the optimal synthesis For the quantum system, we point out that the optimal trajectory can be computed numerically starting from that of the Fuller model. When we are sufficiently close to |3 , we approximate the switching curve of the quantum system by that of its approximation. As described in the main text, we then integrate backwards in time the dynamics of the state and of the adjoint state. We stress that it makes sense to integrate backwards the extremal flow since, according to the results in [1], the map Υ associating with a switching point the subsequent one is hyperbolic and has Γ as stable manifold towards |3 . As a result, the k-th iteration of the inverse of Υ tends towards Γ as k grows.

E.2 A direct optimization method
We apply in this paragraph a direct optimization method to design numerically the solution of the optimal control problem. We use the open access optimal solver BOCOP [25]. A direct optimization approach is a procedure in which the state and the time are discretized in time, transforming the initial optimal control problem into a nonlinear constrained optimization problem. In this optimization protocol, the optimal control problem is slightly modified. The goal is to minimize the cost C = t f 0 x 2 1 (t)dt, while reaching a state as close as possible of the target |3 in a fixed time t f (note that the algorithm does not converge if the time is free). In order to have the fairest comparison possible between the two approaches, we estimate as follows the time of the optimal solution presented in the main text. We assume that the control time is the sum of the time τ taken by the quantum system to go from (0, 1, 0) to a point with first two components equal to (ξ∆x 2 20 , x 20 ), x 20 > 0 small, and the time t Ful f corresponding to the initial condition (ξ∆x 2 20 , x 20 ) for the Fuller problem. For x 20 = 6.9 × 10 −4 we obtain τ = 2.5890 and t Ful f = 0.0011, which leads to t f = 2.5901. We set t f to 2.59 in the numerical optimization procedure. The time subdivision is regular and given by the number of time steps N , going from N = 50 to N = 400. The other optimization parameters are set to their default or recommended values in BOCOP.
In addition of the controls displayed in Fig. 5 of the main text, we report in Fig. 9 additional numerical results about the distance to the target and the cost C with respect to N . When N becomes larger, the numerical optimal process seems to converge towards the optimal solution both in terms of final state and cost. For quite small values of N , we observe that the efficiency of the control is already reasonable. We point out that such sub-optimal strategy could be a possible option to bypass the chattering phenomenon in the optimization procedure.