Robust optimal control of two-level quantum systems

TL;DR

This paper develops optimal control strategies for two-level quantum systems that are robust against uncertainties, using Pontryagin's principle to derive minimal energy and time solutions and analyzing the control landscape complexity.

Contribution

It introduces a method to find global optimal pulses for robustness in quantum control, providing insights into the control landscape's structure and complexity.

Findings

Derived global optimal pulses for robustness against uncertainties.

Established bounds on the control landscape dimension.

Provided estimates of control complexity for robust quantum control.

Abstract

We investigate the time and the energy minimum optimal solutions for the robust control of two-level quantum systems against offset or control field uncertainties. Using the Pontryagin Maximum Principle, we derive the global optimal pulses for the first robustness orders. We show that the dimension of the control landscape is lower or equal to 2N for a field robust to the N th order, which leads to an estimate of its complexity.

Tables

Table 1. Table 1: Robustness of order n 𝑛 n with respect to the offset terms (type δ 𝛿 \delta ) or to the control field inhomogeneities (type α 𝛼 \alpha ) for the local measure. The pulse area and the minimum time are given for the energy and time minimum control problems, respectively.

Type	Cost	Orders	Area or time ($\times \pi$)
$\delta$	time	1/2/3	2.0/2.44/3.54
$\alpha$	time	1/2/3	1.86/2.71/3.56
$\delta$	energy	1/2/3	1.45/1.81/2.11

Equations

$$\dot{\vec{q}}(t)=\left(\begin{array}{ccc}0& \delta & -(1+\alpha )u_y\\ -\delta & 0& (1+\alpha )u_x\\ (1+\alpha )u_y& -(1+\alpha )u_x& 0\end{array}\right)\vec{q}(t),$$

$$\vec{q}(t)={\vec{q}}_0(t)+\delta {\vec{q}}_1(t)+\cdots +{\delta }^N{\vec{q}}_N(t)+O({\delta }^{N+1}),$$

$${\vec{q}}_0(0)={}^t(0,0,1)\mapsto {\vec{q}}_0(t_f)={}^t(0,0,-1),$$

$${\vec{q}}_i(0)={}^t(0,0,0)\mapsto {\vec{q}}_i(t_f)={}^t(0,0,0),\text{for}i\in \{1,\cdots ,N\}.$$

$$\frac{d}{dt}\left(\begin{array}{c}{\vec{q}}_0\\ {\vec{q}}_1\\ {\vec{q}}_2\\ ⋮\\ {\vec{q}}_N\end{array}\right)=\left(\begin{array}{ccccc}{H}_0& 0& 0& \cdots & 0\\ {\partial }_{\delta }H& H_0& 0& & 0\\ 0& {\partial }_{\delta }H& H_0& & 0\\ ⋮& & \ddots & \ddots & ⋮\\ 0& \cdots & 0& {\partial }_{\delta }H& H_0\end{array}\right)\left(\begin{array}{c}{\vec{q}}_0\\ {\vec{q}}_1\\ {\vec{q}}_2\\ ⋮\\ {\vec{q}}_N\end{array}\right),$$

$$H_0=\left(\begin{array}{ccc}0& 0& -u_y\\ 0& 0& u_x\\ u_y& -u_x& 0\end{array}\right),{\partial }_{\delta }H=\left(\begin{array}{ccc}0& 1& 0\\ -1& 0& 0\\ 0& 0& 0\end{array}\right).$$

$$H_P=\sum _{i=0}^N{\vec{p}}_i\cdot {\dot{\vec{q}}}_i+p^0{f}^0.$$

$$H_P={\vec{p}}_0\cdot ({\vec{q}}_0\times \vec{u})+\sum _{k=1}^n{\vec{p}}_k\cdot ({\vec{q}}_k\times \vec{u}+{\vec{q}}_{k-1}\times {\vec{e}}_z)+p^0{f}^0,$$

$$C_E={\int }_0^{t_f}{f}^0(u_x,u_y)dt={\int }_0^{t_f}[u_x^2+u_y^2]dt,$$

$$C_t={\int }_0^{t_f}dt=t_f$$

$$H_P=\vec{u}\cdot {\vec{\Omega }}_0+{\vec{e}}_z\cdot {\vec{\Omega }}_1+p^0{f}^0.$$

$$\dot{\vec{q}}=\frac{\partial H_P}{\partial \vec{p}},\dot{\vec{p}}=-\frac{\partial H_P}{\partial \vec{q}},$$

$$H(\vec{x},\vec{p})=\underset{(u_x,u_y)\in U}{\max }{H}_P(\vec{q},\vec{p},u_x,u_y).$$

$$\{\begin{array}{l}{\dot{\vec{\Omega }}}_0={\vec{\Omega }}_1\times {\vec{e}}_z,\\ {\dot{\vec{\Omega }}}_k={\vec{\Omega }}_k\times {\vec{\Omega }}_0+{\vec{\Omega }}_{k+1}\times {\vec{e}}_z,k\in \{1,\cdots ,N-1\}\\ {\dot{\vec{\Omega }}}_N={\vec{\Omega }}_N\times {\vec{\Omega }}_0,\end{array}$$

$$\{\begin{array}{l}{\dot{\vec{\Omega }}}_0={\vec{\Omega }}_1\times {\vec{e}}_z,\\ {\dot{\vec{\Omega }}}_k=\frac{1}{r}{\vec{\Omega }}_k\times {\vec{\Omega }}_0+{\vec{\Omega }}_{k+1}\times {\vec{e}}_z,k\in \{1,\cdots ,N-1\}\\ {\dot{\vec{\Omega }}}_N=\frac{1}{r}{\vec{\Omega }}_N\times {\vec{\Omega }}_0,\end{array}$$

$$\vec{q}(t)={\vec{q}}_0(t)+\alpha {\vec{q}}_1(t)+\cdots +{\alpha }^N{\vec{q}}_N(t)+O({\alpha }^{N+1}).$$

$$\frac{d}{dt}\left(\begin{array}{c}{\vec{q}}_0\\ {\vec{q}}_1\\ {\vec{q}}_2\\ ⋮\\ {\vec{q}}_N\end{array}\right)=\left(\begin{array}{ccccc}{H}_0& 0& 0& \cdots & 0\\ H_0& H_0& 0& & 0\\ 0& H_0& H_0& & 0\\ ⋮& & \ddots & \ddots & ⋮\\ 0& \cdots & 0& H_0& H_0\end{array}\right)\left(\begin{array}{c}{\vec{q}}_0\\ {\vec{q}}_1\\ {\vec{q}}_2\\ ⋮\\ {\vec{q}}_N\end{array}\right),$$

$$H_0=\left(\begin{array}{ccc}0& 0& -u_y\\ 0& 0& u_x\\ u_y& -u_x& 0\end{array}\right).$$

$${\vec{\Omega }}_k=\sum _{i=k}^N{\vec{\ell }}_{i,i-k}+\sum _{i=k+1}^N{\vec{\ell }}_{i,i-k-1},k\in \{1,\cdots ,N-1\},$$

$${\vec{\Omega }}_N={\vec{\ell }}_{N,0},$$

$$\frac{d}{dt}\left(\begin{array}{c}{\vec{\Omega }}_0\\ {\vec{\Omega }}_1\\ {\vec{\Omega }}_2\\ ⋮\\ {\vec{\Omega }}_N\end{array}\right)=\left(\begin{array}{ccccc}{H}_0& H_0& 0& \cdots & 0\\ 0& H_0& H_0& & ⋮\\ 0& 0& H_0& \ddots & 0\\ ⋮& & & \ddots & H_0\\ 0& 0& 0& \cdots & H_0\end{array}\right)\left(\begin{array}{c}{\vec{\Omega }}_0\\ {\vec{\Omega }}_1\\ {\vec{\Omega }}_2\\ ⋮\\ {\vec{\Omega }}_N\end{array}\right).$$

$$H_P=\vec{u}\cdot {\vec{\Omega }}_0+p^0{f}^0,$$

$$\tilde{H}=\sqrt{{\Omega }_{0x}^2+{\Omega }_{0y}^2}.$$

$$H={\Omega }_{0x}^2+{\Omega }_{0y}^2,$$

$$\{\begin{array}{l}{\dot{\vec{\Omega }}}_0=({\vec{\Omega }}_0+{\vec{\Omega }}_1)\times \vec{u},\\ {\dot{\vec{\Omega }}}_1=({\vec{\Omega }}_1+{\vec{\Omega }}_2)\times \vec{u},\\ {\dot{\vec{\Omega }}}_2=({\vec{\Omega }}_2+{\vec{\Omega }}_3)\times \vec{u},\\ ⋮\\ {\dot{\vec{\Omega }}}_N={\vec{\Omega }}_N\times \vec{u}.\end{array}$$

$${\dot{\vec{q}}}_k(t)=\left(\begin{array}{ccc}0& {\Delta }_k& -u_y\\ -{\Delta }_k& 0& u_x\\ u_y& -u_x& 0\end{array}\right){\vec{q}}_k(t),$$

$${\dot{\vec{\ell }}}_k={\vec{\ell }}_k\times \vec{u}+{\Delta }_k{\vec{\ell }}_k\times {\vec{e}}_z,k\in \{1,\cdots ,N\},$$

$$G_{\text{NOT}}=\left(\begin{array}{ccc}1& 0& 0\\ 0& -1& 0\\ 0& 0& -1\end{array}\right).$$

$$\dot{R}=\left(\begin{array}{ccc}0& \delta & -u_y\\ -\delta & 0& u_x\\ u_y& -u_x& 0\end{array}\right)R.$$

$$R(t)=R_0(t)+\delta R_1(t)+\cdots +{\delta }^N{R}_N(t)+O({\delta }^{N+1}).$$

$$\frac{d}{dt}\left(\begin{array}{c}{R}_0\\ R_1\\ R_2\\ ⋮\\ R_N\end{array}\right)=\left(\begin{array}{ccccc}{H}_0& 0& 0& \cdots & 0\\ {\partial }_{\delta }H& H_0& 0& & 0\\ 0& {\partial }_{\delta }H& H_0& & 0\\ ⋮& & \ddots & \ddots & ⋮\\ 0& \cdots & 0& {\partial }_{\delta }H& H_0\end{array}\right)\left(\begin{array}{c}{R}_0\\ R_1\\ R_2\\ ⋮\\ R_N\end{array}\right),$$

Full text

1 Introduction

Quantum control techniques are nowadays at the core of emergent quantum technologies in a multitude of domains extending from molecular and solid state physics to Nuclear Magnetic Resonance (NMR) and Magnetic Resonance Imaging [1, 2, 3, 4, 5]. One of the main obstructions to the experimental realization of open-loop control processes is their high sensitivity to experimental imperfections and model uncertainties. Since the start of quantum control, this question has motivated the development of pulse design methods addressing such robustness issues [6, 7, 8]. Adiabatic quantum control techniques were first applied with success in some examples, but these protocols have a limited efficiency in a general setting due to the requirement of high energy and long duration fields [9, 10]. Composite pulses [11, 12, 13] and shortcut to adiabaticity techniques [14, 18, 19] have also been proposed, but they cannot reach the physical limits of the dynamical process in terms of time or efficiency. We will show that Optimal Control Theory (OCT) can be a perfect tool to overcome these difficulties [8, 20, 21, 22]. OCT is a general approach allowing to manipulate the system dynamics by determining the control field that minimizes a cost functional which can be, e.g., the control duration or its energy [23, 24, 25]. However, optimal control fields are not robust by construction and this issue is still at the center of a vivid debate [1, 2, 26, 27]. Different numerical approaches ranging from the simultaneous control of an inhomogeneous ensemble of quantum systems [28, 29, 30, 31, 32, 33] to pseudo-spectral methods [34, 35] have been proposed. Only local optimal solutions are obtained, with no certitude about the global optimality of the control process. In this work, we show how this fundamental question can be solved through the Pontryagin Maximum Principle (PMP) [20, 21, 22]. This approach has already been used with success in different optimal quantum problems [36, 37, 38, 39, 40, 41, 42, 43, 44, 45] and we propose here to extend its range of application to the design of robust control protocols. The PMP transforms the optimal control problem into a generalized Hamiltonian system subject to a maximization condition and some boundary constraints. In this framework, the goal consists in finding the Hamiltonian trajectory reaching the target state, while minimizing the cost functional which defines the optimization procedure. A key advantage of the PMP is that it reduces the initial infinite-dimensional control landscape [2, 46] to a finite space of low dimension. As shown below, this property is crucial in the search for globally optimal controls.

In this paper, we establish time and energy minimum optimal control strategies leading to a robust and precise state to state transfer of two-level quantum systems. The implementation of quantum gates is also analyzed. The measure of the robustness is given by the deviation of the control fidelity against offset or field inhomogeneities. This description reproduces the standard experimental uncertainties that can be encountered in quantum information processing, in NMR and in atomic or molecular physics [1, 15, 16, 17]. The robustness is defined either locally by expanding the state of the system order by order with respect to the unknown parameters [18, 19] or globally by considering a discrete inhomogeneous ensemble of quantum systems [28, 29]. A precise definition will be given later. Ultra-precise or broadband excitation profiles are respectively realized based on the first or second measure, which will be called local or broadband robustness below. The two definitions will be considered in the different examples. Note that the controllability of the different systems is assumed and not discussed in this work [47, 48].

The paper is organized as follows. Section 2 introduces the model we study. The optimal solutions for the energy and time-minimum inversion robust against offset inhomogeneities are presented in Sec. 3. Section 4 focuses on the robustness with respect to control field imperfections, while Sec. 5 is dedicated to the broadband control of an ensemble of spins. A comparison with the results obtained with a numerical optimization algorithm is made in Sec. 6. The method is generalized in Sec. 7 to the robust implementation of one-qubit gates. Conclusion and prospective views are given in Sec. 8. Technical computations are reported in the Appendices A, B and C.

2 The model system

We consider the Bloch representation of a two-level quantum system whose dynamics is governed by the Bloch equation. The Bloch vector $\vec{q}(t)={{}^{t}(x,y,z)}$ satisfies the following differential system:

[TABLE]

where $u_{x}$ and $u_{y}$ are the two control fields. The parameters $\delta$ and $\alpha$ represent respectively the offset and control field inhomogeneities. We first consider the case where $\alpha=0$. We assume that the solution of the Bloch equation can be written as a perturbative expansion in $\delta$ up to a given order:

[TABLE]

with $\vec{q}_{i}={{}^{t}(x_{i},y_{i},z_{i})}$. The vector $\vec{q}_{0}$ is the homogeneous part of the solution and $\vec{q}_{i}$ the inhomogeneous contribution due to the $i$th- order term of the expansion. We investigate the robust control of the inversion of the Bloch vector, i.e. the goal is to bring in a time $t_{f}$ the state $\vec{q}(t)$ from the north pole to the south pole of the Bloch sphere. Other initial or target states can be analyzed in the same way. The control process can be expressed in the perturbative expansion as:

[TABLE]

Note that the target states of the inhomogeneous contributions ensure that the offset term does not modify the final state of the system up to the order $N$ in $\delta$, improving thus the robustness of the control protocol. Plugging Eq. (2) into Eq. (1), it is straightforward to show that the differential system governing the dynamics of each vector $\vec{q}_{i}$ is given by:

[TABLE]

where:

[TABLE]

The originality of the method consists in directly solving the PMP applied to Eq. (4). We show that the optimal control fields must satisfy a certain differential system depending upon a finite number of parameters for a robustness at order $N$. This number can be interpreted as the dimension of the control landscape. A careful investigation of this landscape allows us to detect the global optimal solution of the process, at least for low robustness orders.

3 Energy and time minimum inversion of two-level quantum systems robust against offset inhomogeneities

We consider in this paragraph the robust optimal inversion with respect to offset uncertainties. In this case, the PMP is formulated from the pseudo-Hamiltonian $H_{P}$ which can be written as follows [20, 21, 22]:

[TABLE]

This leads to:

[TABLE]

where $\vec{p}_{i}$ is the adjoint state of $\vec{q}_{i}$ and $p_{0}$ a negative constant, which is set to $-1/2$ and to $-1$ for the energy and time optimal control problems, respectively [20, 21, 22]. $f^{0}$ is a function of $u_{x}$ and $u_{y}$ whose integral over time gives the associated cost functional $C$ to minimize. We have:

[TABLE]

for the energy and

[TABLE]

for the time, where the control duration $t_{f}$ is not fixed. In Eq. (7), the vectors $\vec{u}$ and $\vec{e}_{z}$ have the coordinates $(u_{x},u_{y},0)$ and $(0,0,1)$ and $\times$ denotes the vector product of two three-dimensional vectors.

We introduce the angular momenta $\vec{\ell}_{a,b}=\vec{p}_{a}\times\vec{q}_{b}$ and the partial sums $\vec{\Omega}_{n}=(\Omega_{nx},\Omega_{ny},\Omega_{nz})=\sum_{i=k}^{n}\vec{\ell}_{i,i-k}$. Using the properties of the scalar triple product, we arrive at:

[TABLE]

The PMP states that the coordinates of the Bloch vector $\vec{q}$ and of the corresponding adjoint state $\vec{p}$ fulfill the Hamiltonian’s equations associated with $H_{P}$:

[TABLE]

the control fields being given by the maximization condition [20, 21, 22]:

[TABLE]

The set $U$, which defines the constraint on the pulses, is given by $U=\mathbb{R}^{2}$ and by $u_{x}^{2}+u_{y}^{2}\leq 1$ for the energy and time minimization problems, respectively. Note that the constraint $u_{x}^{2}+u_{y}^{2}\leq 1$ avoids the occurrence of very intense control fields which are not relevant experimentally. We obtain $H=\frac{\Omega_{0x}^{2}+\Omega_{0y}^{2}}{2}+\Omega_{1z}$, with $u_{x}=\Omega_{0x}$ and $u_{y}=\Omega_{0y}$, for the energy-minimum and $H=r+\Omega_{1z}$, with $u_{x}=\Omega_{0x}/r$, $u_{y}=\Omega_{0y}/r$ and $r=\sqrt{\Omega_{0x}^{2}+\Omega_{0y}^{2}}$, for the time-minimum problem.

Using the maximization condition of the PMP, a straightforward computation then leads to:

[TABLE]

in the energy case and to

[TABLE]

for the time-minimum problem. In the two situations, we deduce that $\Omega_{0z}(t)=cst$. This constraint becomes $\Omega_{0z}(t)=0$ when the initial point is the north pole of the Bloch sphere. The differential systems (8) and (9) can be interpreted as the conditions to satisfy by the control fields to realize the control process. Moreover, the initial phase of the control fields is irrelevant, which means that we can set $\Omega_{0y}(t=0)=0$. Another point is related to the fact that the norm $|\vec{\Omega}_{N}|$ is constant in time, and can be set to $1$ without loss of generality (this is equivalent to a time rescaling). We obtain that $\vec{\Omega}_{N}(t=0)$ only depends on one angle, i.e. $\vec{\Omega}_{N}(0)=(\cos\vartheta,\sin\vartheta,0)$. Finally, the field depends on $2N$ parameters at order $N$: $\Omega_{0x}(0)$, $\Omega_{kx}(0)$ and $\Omega_{ky}(0)$ for $k\in\{1,\cdots N-1\}$, and the angle $\vartheta$. This number of parameters is also the dimension of the control landscape. The last point to solve is to adjust these parameters to realize the transfer (3). Many different solutions can exist. For $N$ small enough, the low dimension of the control landscape allows us to find the global optimal solution minimizing $C_{E}$ or $C_{t}$.

We now describe the results for the energy-minimum problem. We consider here only one control field along $\vec{e}_{x}=(1,0,0)$, i.e. $u_{y}(t)=0$. We compute the global optimal solutions for the first, second and third robustness orders. The analytical expression of the control field is derived at first order, while numerical optimization techniques are used for the second and third order. The computations are detailed in the Appendix A. Figure 1 displays the control fields and the homogeneous contribution $\vec{q}_{0}$ of the Bloch vector.

We observe the peculiar dynamics of the homogeneous part of the Bloch vector, which oscillates in the $(y_{0},z_{0})$ plane like a damped oscillator, with a number of oscillations increasing with the robustness order in $\delta$. The efficiency of the control protocol is shown in Fig. 2. As could be expected, a better robustness is achieved when higher orders are nullified.

In the time-optimal case, it can be shown that the minimum time to cancel the first robustness order is associated with a pulse of constant intensity along the $x$- direction of the Bloch sphere. This pulse switches between the values 1 and -1 at time $t=3\pi/2$ and it has a total duration of $t_{f}=2\pi$. The pulse structure is not the same at second and third orders. The robust inversion is realized by smooth fields of durations $t_{f}=2.44\pi$ and $t_{f}=3.54\pi$ for the second and third order robust control fields. The details of the computations are given in Appendix B. Figure 3 displays the fields $u_{x}$ and $u_{y}$ and the homogeneous contribution $\vec{q}_{0}(t)$ of the Bloch vector. We recover with this geometrical analysis the transition found numerically in [28, 29] from a square signal to a smooth field.

The different numerical results are collected in Tab. I.

4 Robustness against control field inhomogeneities

The same method can be applied to the field inhomogeneities and the $\alpha$- parameter. We consider the on-resonance case where $\delta=0$. We assume that the solution can be written as a perturbative expansion of the form:

[TABLE]

We obtain the following differential system:

[TABLE]

with:

[TABLE]

The pseudo-Hamiltonian of this system is given by $H_{p}=\sum_{i=0}^{N}\vec{p}_{i}\cdot\dot{\vec{q}}_{i}+p^{0}f^{0}$, where $p^{0}$ and $f^{0}$ depend on the cost functional to minimize. Introducing the angular momenta $\vec{\ell}_{ij}=\vec{p}_{i}\times\vec{q}_{i}$, and the vectors

[TABLE]

we can show, using the Hamilton’s equations, that:

[TABLE]

Pontryagin’s Hamiltonian can be written as follows:

[TABLE]

where $\vec{u}=(u_{x},u_{y},0)$. In this case, the PMP leads to the same condition both for the energy-minimum and time-minimum control problems. For the minimization of the time, the maximization condition gives $u_{x}=\frac{\Omega_{0x}}{\sqrt{\Omega_{0x}^{2}+\Omega_{0y}^{2}}}$ and $u_{y}=\frac{\Omega_{0y}}{\sqrt{\Omega_{0x}^{2}+\Omega_{0y}^{2}}}$. We deduce that the Pontryagin’s Hamiltonian reads (after absorbing the constant $p^{0}f^{0}=-1$ in the definition of the pseudo-Hamiltonian):

[TABLE]

If $\tilde{H}\neq 0$, we can consider the Hamiltonian:

[TABLE]

which is set to $H=1$ without loss of generality. The control fields are then given by $u_{x}=\Omega_{0x}$ and $u_{y}=\Omega_{0y}$. The application of the PMP in the energy minimum case leads to the same expression. The first consequence of the PMP is that $u_{x}^{2}+u_{y}^{2}=1$ for any time $t$. The differential system can be written as:

[TABLE]

At time $t=0$, we have $\Omega_{kz}(0)=0$ since the initial state of the Bloch vector is the north pole of the Bloch sphere. Moreover, the initial phase of the control field is irrelevant for the inversion transfer, leading to $u_{y}(0)=\Omega_{0y}(0)=0$ and $u_{x}(0)=\Omega_{0x}(0)=1$. The dimension of the control landscape is $2N$, and this space can be parameterized by $(\Omega_{kx}(0),\Omega_{ky}(0))$ for $k\in\{1,\cdots,2N\}$. The results for the orders 1, 2 and 3 are displayed in Fig. 4. The analytical expression of the first-order robust solution is given in the appendix C.

The minimum times are given in Tab. I. Note the linear evolution of these times as a function of the robustness order.

5 Broadband robust optimal control

The robustness can also be defined through an inhomogeneous ensemble of quantum systems featured by different parameters. A finite number of systems belonging to this set is considered [28, 29, 31, 32, 33]. Our approach also works in this situation and leads to the same control landscape complexity. For an ensemble of quantum systems with different offsets, the dynamics of each element is of the form:

[TABLE]

where $\Delta_{k}$ is the offset of the system $k$. We define the individual angular momentum $\vec{\ell}_{k}=\vec{p}_{k}\times\vec{q}_{k}$, where $\vec{p}_{k}$ is the adjoint state of the system $k$. The optimal control fields which invert simultaneously all the elements of the set satisfy:

[TABLE]

with $u_{x}=\sum_{k}\ell_{k,x}$, $u_{y}=\sum_{k}\ell_{k,y}$ for the energy-minimum problem, and $u_{x}=\frac{1}{r}\sum_{k}\ell_{k,x}$, $u_{y}=\frac{1}{r}\sum_{k}\ell_{k,y}$ with $r=\sqrt{\left(\sum_{k}\ell_{k,x}\right)^{2}+\left(\sum_{k}\ell_{k,y}\right)^{2}}$ for the time-minimum case. Since the initial point of each system is the north pole of the Bloch sphere, we deduce that $\ell_{k,z}=0$. A trajectory depends on $2N$ parameters but this number can be reduced to $2N-2$. Indeed, since the initial phase of the control is arbitrary, we can choose $u_{y}(0)=0$ which leads to $\sum_{i=1}^{N}\ell_{iy}(0)=0$. Moreover, one of the first integrals can be set to $1$ by rescaling the time. Here we choose $H=1$, which leads to $r(0)=\sum_{i=1}^{N}\ell_{ix}(0)=1$. Thus, we have $\ell_{Nx}(0)=1-\sum_{i=1}^{N-1}\ell_{ix}(0)$ and $\ell_{Ny}(0)=-\sum_{i=1}^{N-1}\ell_{iy}(0)$. Finally, we deduce that the total dimension of the control landscape is $2N-2$.

As displayed in Fig. 5, the time-optimal solutions for two, three and four quantum systems are very similar to the trajectories of Fig. 3. For $N=2$, we recover the results derived in Ref. [49]. The right panel of Fig. 2 displays the inversion profile against the offset $\delta$. Note the very large robustness obtained although only four quantum systems are considered.

6 Comparison with the GRAPE algorithm

The broadband inversion has been investigated numerically in Ref. [28, 29]. The authors used the GRAPE algorithm [23], which is based on the PMP. However, the pulses were not optimized in the same way, in the sense that this algorithm aims at controlling a large number of spins ($\simeq$ 100) over a certain range of offset inhomogeneities and for different pulse durations. Here, the goal is to compare this method to the time-optimal control of two, three and four spins. We recall that the control fields are so that $u_{x}=\cos\Phi(t)$ and $u_{y}=\sin\Phi(t)$. As a consequence, we use the GRAPE algorithm to optimize only the phase of the control fields (see Fig. 5 of Ref. [28]).

We proceed as follows. We have determined the optimal time for controlling two spins with $\Delta_{i}\in\{-1/2,1/2\}$, three spins with $\Delta_{i}\in\{-1/2,0,1/2\}$ and four spins with $\Delta_{i}\in\{-1/2,-1/6,1/6,1/2\}$ (see Fig. 5). We optimize a pulse of the same duration with GRAPE for 100 spins so that $\Delta\in[-1/2,1/2]$. A guess field is chosen for the optimization algorithm in order to recover the results of Fig. 5 of Ref. [28]. The result is presented in Fig. 6.

We observe a strong similarity between the GRAPE solution and the global one with two, three and four spins. This result was not obvious since the GRAPE algorithm optimizes a large number of spins. Another interesting point is the robustness profile achieved with each pulse. We compute the fidelity $-z(t_{f})$ obtained by integrating Eq. (19) for 1000 spins with an offset $\Delta\in[-0.6,0.6]$. The result is displayed in Fig. 7.

A remarkable point is the fact that the number of peaks of the robustness profiles is the same if a small or a large number of spins is considered. This suggests a way to limit the computational time of robust control sequences.

7 Implementation of one-qubit quantum gates

The method presented in this work can be adapted to the implementation of quantum gates. A quantum gate can be written as a rotation matrix belonging to $SO(3)$. As an example, a NOT- gate is associated with a rotation of angle $\pi$ about the axis $\vec{e}_{x}$ (or $\vec{e}_{y}$) of the Bloch sphere. This matrix is given by:

[TABLE]

The Bloch equation can be derived for this rotation matrix, which is equivalent to consider the dynamics of three orthogonal Bloch vectors. This matrix will be called below, the Bloch matrix. In the case of offset inhomogeneities, the evolution of the Bloch matrix is governed by the following equation:

[TABLE]

We assume that the solution of the system (22) can be expressed as a perturbative expansion in terms of $\delta$:

[TABLE]

The matrix $R_{0}$ is the homogeneous part of the solution. It can be shown that the dynamics of the system is of the form:

[TABLE]

with:

[TABLE]

At time $t=0$, the matrix $R_{0}(0)$ is the identity matrix and the matrices $R_{k\geq 1}$ are zero. A robust gate $G$ is a transfer of the form:

[TABLE]

Each matrix $R_{k}$ has nine components that we denote by $a^{(k)}_{ij}$, i.e.:

[TABLE]

We introduce the vector $\vec{q}_{k}$ (with nine-elements) defined as the concatenation of the three columns of the matrix $R_{k}$, i.e. $\vec{q}_{k}={{}^{t}(a_{11}^{(k)},a_{21}^{(k)},\cdots,a_{33}^{(k)})}$. The pseudo-Hamiltonian $H_{P}$ of the system is given by:

[TABLE]

where $p^{0}f^{0}$ depends on the cost functional of the system. The components of the adjoint state $\vec{p}$ are given by: $\vec{p}_{k}={}^{t}(b_{11}^{(k)},b_{21}^{(k)},\cdots,b_{33}^{(k)})$. We also define some three-dimensional angular momentum vectors $\vec{\ell}^{(i,j)}$ as:

[TABLE]

and we introduce the vectors $\vec{\Omega}_{k}$ defined as:

[TABLE]

Finally, we can show that the pseudo-Hamiltonian can be expressed in terms of the vectors $\vec{\Omega}_{k}$ as follows:

[TABLE]

with $\vec{u}=(u_{x},u_{y},0)$. Using the Hamilton’s equations, we can deduce the dynamics of each $\vec{\Omega}_{k}$. We find exactly the same equation as Eq. (8), that is:

[TABLE]

We can then apply the PMP. For the time-minimum problem, the maximization of the pseudo-Hamiltonian leads to $u_{x}=\tfrac{\Omega_{0x}}{r}$ and $u_{y}=\tfrac{\Omega_{0y}}{r}$, with $r=\sqrt{\Omega_{0x}^{2}+\Omega_{0y}^{2}}$. Substituting the expressions of $u_{x}$ and $u_{y}$ into Eq. (33), it is then straightforward to show that $\Omega_{0z}$ is constant, but this constant can be different from zero. We set $\Omega_{0z}=I$. The time-optimal control field is solution of the following system of equations:

[TABLE]

Pontryagin’s Hamiltonian is constant and is given by:

[TABLE]

A difference from a state to state control problem is that, in general, $\Omega_{kz}\neq 0$. The dimension of the control landscape can only be reduced by one by noting that $|\vec{\Omega}_{N}|$ is a constant, which can be set to $1$ without loss of generality. Thus, a solution of the system (34) depends on $3N+2$ parameters which are $\Omega_{0x}(0)$, $\Omega_{0y}(0)$, $\Omega_{0z}=I$, $\Omega_{kx}(0)$, $\Omega_{ky}(0)$, $\Omega_{kz}(0)$, $\vartheta$ and $\varphi$, with $\Omega_{Nx}(0)=\sin\vartheta\cos\varphi$, $\Omega_{Ny}(0)=\sin\vartheta\sin\varphi$ and $\Omega_{Nz}(0)=\cos\vartheta$. Figure 8 displays the control fields robust to first and second order in $\delta$. Further work will be necessary to improve this first solution.

8 Conclusion

Using the PMP, we have derived the global robust optimal control strategies (for the lowest orders) for the inversion transfer in the energy and time- minimum cases. The derived pulses have an explicit and relatively simple form which is easily implementable experimentally. The analytical expression of some control fields is also obtained. We stress that the global optimality of the solutions of this work is in sharp contrast with the fields designed by numerical methods, which correspond to local maxima.

These results can be viewed as a first step towards a complete answer of the robustness issue, which is a long-standing problem in quantum control. They also pave the way to other studies using the same approach, such as a transfer robust with respect to the two inhomogeneous parameters or the design of robust propagators, which will be interesting in quantum computing. A first step in this direction is made in Sec. 7 with the derivation of a NOT gate robust to the second order in $\delta$. The main obstacle in these two generalizations will be the dimension of the control landscape, which makes the search for global optimal controls difficult. Another interesting issue would be to generalize this approach to the robustness against noise for which the experimental uncertainties are not constant in time.

**ACKNOWLEDGMENT

**S.J. Glaser acknowledges support from the DFG (Gl 203/7-2). D. Sugny and S. J. Glaser acknowledge support from the ANR-DFG research program Explosys (ANR-14-CE35-0013-01). D. Sugny acknowledges support from the PICS program and from the ANR-DFG research program COQS (ANR-15-CE30-0023-01). The work of D. Sugny has been done with the support of the Technische Universität München – Institute for Advanced Study, funded by the German Excellence Initiative and the European Union Seventh Framework Programme under grant agreement 291763.

Appendix A Application of the PMP in the energy minimum case

We consider the control of the population inversion in the energy minimum case. We show how to determine the control field robust at order one, two and three. The general method consists in first computing a solution of the system (8), and then using the control fields $u_{x}=\Omega_{0x}$ and $u_{y}=\Omega_{0y}$ in the system (1) in order to realize a robust inversion. We add the constraint to use only one control field by setting $u_{y}=0$. Note that better results could be achieved at second order (but not at first order) when two fields are considered.

A.1 Analytical derivation at first order

For $N=1$, the system of Eq. (8) becomes:

[TABLE]

Setting $u_{y}=\Omega_{0y}=0$, we can show that the system simplifies to:

[TABLE]

It has two constants of motion:

[TABLE]

We can set $\ell=1$ without loss of generality. These two constants are associated with two surfaces in the $(\Omega_{0x},\Omega_{1y},\Omega_{1z})$- space. The conservation of $H$ corresponds to a parabolic plane, and the one of $\ell$ to a cylinder of radius $\ell=1$ along the $\Omega_{0x}$- direction. These two surfaces are represented in Fig. 9. The solution of the system (37) belongs to the intersection of these two surfaces. The plot of this intersection for different values of $H$ leads to the phase portrait of the system (37) depicted on the cylinder. We get the same phase portrait as the one of a planar pendulum, with three families of solution [50]. The rotating solutions occur for $H>1$, the oscillating solutions for $-1\leq H\leq 1$, and the separatrix for $H=1$. The trajectories for which $-1\leq H<0$ are not solutions of the problem, since we have in this case $\Omega_{1z}<0$, which is not possible since $\Omega_{1z}(0)=0$.

The solution of the system (37) can be expressed in terms of the Jacobi amplitude function $\operatorname*{am}(u,m)$ [51], which is an angle defined as the inverse of the Incomplete Elliptic Integral of the first kind $\operatorname*{F}(u,m)=\int_{0}^{u}\frac{dt}{\sqrt{1-m\sin^{2}t}}$. Note the relation $d\operatorname*{am}(u,m)/du=\sqrt{1-m\sin^{2}(\operatorname*{am}(u,m))}$ which is useful for integrating the following equations. We only consider the oscillating trajectories which contain the global optimum of the problem. The explicit solution of the system (37) can be expressed in the oscillating case as $\Omega_{0x}=2\sqrt{m}\cos\nu$, $\Omega_{1y}=-2\sqrt{m}\sin\nu\sqrt{1-m\sin^{2}\nu}$ and $\Omega_{1z}=2m\sin^{2}\nu-1$, where $\nu(t)=\operatorname*{am}(t+\rho,m)$, $\rho=\pm\operatorname*{F}\left(\arcsin\left(\tfrac{1}{\sqrt{2m}}\right),m\right)$ and $m=\tfrac{1+H}{2}$. The initial phase $\rho$ is computed so as $\Omega_{10z}=0$. Note that the initial conditions are so that $\Omega_{0x}(0)=\sqrt{2H}$ and $\Omega_{1y}(0)=\pm 1$, the sign being the same as the sign of $-\rho$.

Each trajectory of the phase portrait is associated with a control $\Omega_{0x}(t)$ which is a candidate to realize a robust population inversion. Since a control field depends only on $H$ (and eventually on the sign of $\rho$), the goal is to find the parameter $H$ which inverts the population with a minimum of energy. For that purpose, we integrate the coordinates $\vec{q}_{i}(t)$. At first order in $\delta$, the differential system governing the dynamics of $\vec{q}_{0}=(x_{0},y_{0},z_{0})$ and $\vec{q}_{1}=(x_{1},y_{1},z_{1})$ is given by:

[TABLE]

Since we have $\Omega_{0y}(t)=\Omega_{0z}(t)=0$ and $\vec{q}_{1}(0)=(0,0,0)$, the system simplifies into:

[TABLE]

All the other coordinates are equal to zero. Introducing the angle $\theta(t)$ so that $\sin\theta=y_{0}$ and $\cos\theta=z_{0}$, the problem consists in computing the parameter $H$ which realizes the transfer:

[TABLE]

The differential system is given in these coordinates by:

[TABLE]

The solution of this system can be written as:

[TABLE]

where the sign is the same as the one of $\rho$. The solution $x_{1}(t)$ could be expressed as the sum of a linear term and an Elliptic Integral of the second kind [51]. Figure 10 shows the trajectories $x_{1}=f(\theta)$ starting from $(0,0)$ for different values of $H$.

Note that Fig. 10 represents some trajectories starting from the north pole of the Bloch sphere driven by the control fields solution of the PMP. Optimal robust transfers can be found from Fig. 10. As an example, it can be seen that the point $(\theta=-\pi/2,x_{1}=0)$ belongs to some of the trajectories, which corresponds to an optimal robust excitation transfer.

For the inversion, the optimal time $t^{*}$ is computed so that $\theta(t^{*})=\pm\pi$. We find:

[TABLE]

where $\operatorname*{K}$ is a complete elliptic integral of the first kind [51]. A robust control is achieved by finding the parameter $H$ for which $\int_{0}^{t^{*}}\sin\theta dt=0$, in order to cancel the contribution to the first order term in $\delta$. We get:

[TABLE]

for the two possible signs of $\rho$. The corresponding trajectories in the $(\theta,x_{1})$- plane are represented in Fig. 10, for $\rho>0$ and $\rho<0$. The two optimal solutions are equivalent and give the global optimum. The corresponding pulse $u_{x}=\Omega_{0x}$ is plotted in Fig. 1. Its area is given by:

[TABLE]

Note that other local optima exist and can be used to realize the population inversion.

A.2 Derivation of the control field at second order

We compute in this paragraph the solution at second order with one control field $u_{x}$. Better solutions can be achieved if a second field is considered. The second order solution satisfies the following differential system:

[TABLE]

This system has six constants of the motion given by:

[TABLE]

In the case $u_{y}=\Omega_{0y}=0$, we obtain that $K=I=0$ and that $\Omega_{1x}=\Omega_{2y}=\Omega_{3y}=0$, which also means that $\Omega_{2x}=-1$ (note that the sign could be set to $+1$, but the result is equivalent). The system then becomes:

[TABLE]

It has two constants of motion:

[TABLE]

The conservation of $H$ is associated with the same parabolic plane as the one represented in Fig. 9. The conservation of $J$ is described by a symmetric paraboloid of axis $\Omega_{0x}$ as shown in Fig. 11. The solution belongs to the intersection of these two surfaces, and the phase portrait is obtained by drawing it for every pair $(H,J)$. We plot the phase portrait for two different values of $J$ in Fig. 11.

The system (49) can be integrated by using the fact that:

[TABLE]

The solution $\Omega_{0x}$ can be expressed in terms of Jacobi’s elliptic functions [52]. However, we use here a numerical analysis to find the optimal pulse.

A solution of the system (49) depends on two parameters which are the initial conditions $\Omega_{0x}(0)$ and $\Omega_{1y}(0)$ ($\Omega_{1z}(0)$ is equal to zero by construction). These parameters can be related to $H$ and $J$ through the relations $\Omega_{0x}=\sqrt{2H}$ and $\Omega_{0y}(0)=\sqrt{2J+2\sqrt{2H}}$. In other words, the control landscape is two-dimensional and is parameterized by $H$ and $J$.

We introduce the following time-dependent fidelity:

[TABLE]

where $\vec{q}_{0T}=(0,0,-1)$ and $\vec{q}_{kT,k>0}=(0,0,0)$ are the target states for a robust population inversion. The general method can be described as follows. For a set of parameters $(H,J)$, we integrate numerically the system (49) until an arbitrary time $t_{f}$. We then use the solution $\Omega_{0x}(t)$ as a control field in the system (1). We integrate this system until the time $t=t_{f}$ and we compute the fidelity $F(t)$ of Eq. (52). This leads to the time $t^{*}$ for which the fidelity $F(t)$ is maximum, and we denote by $F(t^{*})=F^{*}$ the corresponding fidelity. We also obtain the area of the control field $A^{*}=\int_{0}^{t^{*}}|\Omega_{0x}(t)|dt$. Since the control landscape is a two-dimensional space, the different quantities $F^{*}$, $t^{*}$ and $A^{*}$ can be determined for every couple $(H,J)$. It is then straightforward to find the global optimal solution of the control problem. This approach is shown in Fig. 12.

It can be seen in Fig. 12 that $F^{*}=0$ is satisfied for many pairs $(H,J)$. One solution is associated with the smallest energy (and with a very small area $A^{*}$). This solution is the global optimum of the control problem, and is given by:

[TABLE]

From a numerical point of view, the target is reached with a precision of $F^{*}=-3.4\times 10^{-6}$, and the area of the pulse is $A^{*}=1.81\pi$. The control field $u_{x}=\Omega_{0x}$ and the homogeneous solution $\vec{q}_{0}$ are represented in Fig. 1.

A.3 Derivation of the control field at third order

With only one control field $u_{x}$, we can show that the system governing the optimal control field is such that $\Omega_{0y}=\Omega_{1x}=\Omega_{2y}=\Omega_{2z}=\Omega_{3x}=0$. The other coordinates satisfy:

[TABLE]

The system has the four following first integrals:

[TABLE]

Since the Bloch vector starts on the north pole of the sphere at time $t=0$, we have $\Omega_{kz}(0)=0$ for all $k$. We thus have $\Omega_{3y}(0)=-1$ (we could also choose $+1$), and a solution of the system (54) depends on the three following parameters:

[TABLE]

which involves that the control landscape is $3$-dimensional. An analytical solution of this system is difficult to compute. The corresponding parameters for a robust inversion are determined with a numerical gradient algorithm. We find: $\Omega_{0x}(0)=1.2384$, $\Omega_{1y}(0)=2.9848$, $\Omega_{2x}(0)=-2.8019$ with and $t^{*}=2.43\pi$ and $F^{*}=-9\times 10^{-8}$. The area of the control field is:

[TABLE]

The control fields and the homogenous solution $\vec{q}_{0}(t)$ are depicted in Fig. 1.

Appendix B Application of the PMP to the time-minimum case

In this section, we apply the PMP to the time-minimum case.

For the optimal control problem at first order, the differential system can be expressed as:

[TABLE]

In addition to the Pontryagin’s Hamiltonian, this system has two first integrals of the form $|\vec{\Omega}_{1}|\equiv 1$ and $I=\vec{\Omega}_{0}\cdot\vec{\Omega}_{1}$. The control landscape of the system (58) is parameterized by $\Omega_{0x}(0)$ and $\vartheta$ so that $\Omega_{1x}(0)=\cos\vartheta$ and $\Omega_{1y}(0)=\sin\vartheta$. A numerical analysis shows that the global optimum occurs for $\vartheta=\pi/2$. Since the quantity $I=\vec{\Omega}_{0}\cdot\vec{\Omega}_{1}$ is constant, and is equal to zero in this optimal situation, we deduce that $\vec{\Omega}_{1}$ is perpendicular to $\vec{\Omega}_{0}$ for any time $t$. This constraint can be fulfilled if and only if $\Omega_{0y}(t)=0$, that is $u_{y}(t)=0$. The control is therefore of the form $u_{x}=\operatorname*{sign}(\Omega_{0x})$.

The robust time-optimal solution at first order is thus a pulse along the axis $\vec{e}_{x}$ of constant amplitude. Some switches between the maximum and the minimum values of the field can exist. In this case, the differential system becomes:

[TABLE]

and the first integrals are given by:

[TABLE]

Note that the conservation of $H$ is associated with two semi-planes intersecting along the line of equation $(\Omega_{0x},\Omega_{0z})=(0,H)$, as shown in Fig. 13.

If $H>1$, the sign of $\Omega_{0x}$ remains the same for any time $t$ which means that no switch appears. The control field $u_{x}$ is a constant pulse and it cannot realize a robust inversion. We now analyze the case $H\in[0,1]$, and we explicitly derive a solution starting at time $t=0$ from the point $\Omega_{0x}(0)=H$ and $\Omega_{0y}(0)=1$. We consider one period, that is one cycle of a trajectory. For one period, we can see that the sign of $\Omega_{0x}$ changes two times. We denote by $T_{1}$ and $T_{2}$ the times of the first and second switches. The solutions can be expressed as follows:

[TABLE]

for $t\in[0,T_{1}]$,

[TABLE]

for $t\in[T_{1},T_{2}]$ and

[TABLE]

for $t\in[T_{2},T]$. Note that the solutions of the system (59) are given for any initial condition so that $\Omega_{0x}(0)=H$ and $\Omega_{0y}(0)=1$. The times $T_{1}$ and $T_{2}$ correspond to the switches of $\Omega_{0x}$. $T$ is the period of the solution. We have $T_{1}=\pi+\arctan\left(\tfrac{H}{\sqrt{1-H^{2}}}\right)$, $T_{2}=3T_{1}-\pi$ and $T=4T_{1}-2\pi$. The evolution of $\Omega_{0x}(t)$ is displayed in Fig. 14.

To summarize, the control field is a pulse of constant amplitude with switches at times $T_{1}$ and $T_{2}$ given in terms of $H$. The final goal is then to determine the value of $H$ which allows us to realize a robust transfer in minimum time. We start by integrating the dynamics of the states $\vec{q}_{0}(t)$ and $\vec{q}_{1}(t)$ in the system (1) with a control field given by $u_{x}=\operatorname*{sign}(\Omega_{0x})$. The system is of the form:

[TABLE]

Introducing the angle $\theta=\arctan(y_{0}/z_{0})$, we can show that the system becomes:

[TABLE]

In these coordinates, a robust transfer is of the form $\theta(0)=0\rightarrow\theta(t_{f})=\pi$ and $x_{1}(0)=0\rightarrow x_{1}(t_{f})=0$. The solutions of this equation are:

[TABLE]

for $t\in[0,T_{1}]$,

[TABLE]

for $t\in[T_{1},T_{2}]$ and

[TABLE]

for $t\in[T_{2},T]$.

We plot in Fig. 15 the solutions in the plane $(\theta,x_{1})$ for different values of $H$ in order to find the optimal solution.

We obtain that the global optimum corresponds to a very particular case, which occurs in the limit $H\to 1$. The corresponding control field is associated with a trajectory infinitely close to the separatrix of Fig. 13. Note that $H$ cannot be exactly equal to $1$, because, in this case, there is no switch. The solution satisfies $T_{1}\to 3\pi/2$ and $t^{*}\to 2\pi$. This point is shown in Fig. 3 of the main text. The duration $t^{*}=2\pi$ is thus the physical minimum time to make a first-order robust inversion with a control field bounded by $1$. The computation of the second and third order robust parameters is made with a numerical gradient algorithm.

Appendix C Robustness against control field inhomogeneities

In this paragraph, we derive the robust optimal field at first order. The differential system to solve is given by:

[TABLE]

Note that the control landscape is parameterized by the two initial values of $\Omega_{1x}(0)$ and $\Omega_{1y}(0)$. In addition to the Pontryagin’s Hamiltonian, this system has $5$ first integrals which are $\Omega_{0z}-\Omega_{1z}=0$, $I_{x}=\Omega_{0x}-2\Omega_{1x}$, $I_{y}=\Omega_{0y}-2\Omega_{1y}$, $J=\vec{\Omega}_{0}\cdot\vec{\Omega}_{1}$ and $M=|\vec{\Omega}_{1}|$. The relation $\Omega_{0z}-\Omega_{1z}=0$ is due to the initial conditions $\Omega_{kz}(0)=0$ for all $k$. Using these constants, the system becomes:

[TABLE]

Since $\Omega_{0x}(0)=1$ and $\Omega_{0y}(0)=\Omega_{0z}(0)=0$, a solution depends only on the two parameters $I_{x}$ and $I_{y}$ and we can show that $2J=1-I_{x}$. This system can be integrated in terms of Jacobi’s elliptic functions. We introduce the parameters $\omega$ and $m$ so that:

[TABLE]

and the Jacobi’s amplitude function:

[TABLE]

We also introduce the phase of the control field $\phi$ so that $\Omega_{0x}=\cos\phi$ and $\Omega_{0y}=\sin\phi$. It can be checked that the following functions are solutions of the problem:

[TABLE]

At time $t=t_{f}$, we must have $\Omega_{0z}(t_{f})=0$ since the Bloch vector reaches the south pole of the Bloch sphere. Furthermore, $\Omega_{0z}$ goes back to zero when the time $t$ is a multiple of $2\operatorname*{K}(m)/\omega$.

In order to locate the global optimum, we use a similar method as in Sec. A.2. We integrate the system (11) with the control fields (73) and we compute the corresponding fidelity for every pair $(I_{x},I_{y})$. We obtain that the optimal control field verifies:

[TABLE]

The result is presented in the left panels of Fig. 4. We use a numerical algorithm to compute the solutions robust at second and third order.