Fast state transfer in a Λ-system: a shortcut-to-adiabaticity approach to robust and resource optimized control

Precise control of quantum systems is required to realize a number of applications in quantum information and precision quantum measurements [1–3]. The simultaneous fulfillment of constraints on protocol duration, fidelity, robustness against parameter variations, and feasibility to implement the desired interaction, can be handled with optimal control theory [4–8]. Here, a cost functional is defined that quantifies the quality of a solution and penalizes, e.g., extreme values or strong variations in the control fields [9, 10]. This cost functional is then numerically optimized. Building on pioneering efforts on selective excitation of molecular systems and NMR [11], optimal control theory has been successfully applied in various quantum control problems, including manipulation of Bose–Einstein condensates [12, 13] and transport of single atoms in optical tweezers [14]. Benign control problems permit the use of gradient methods on a large number of control parameters. However, for time constrained problems, these methods may converge on sub-optimal solutions, and thus require multi-start approaches with different initial guesses. This process is numerically inefficient and offers no guarantee for the identification of the optimal solution [14]. We recall that the evaluation of a cost function involving the transfer fidelity requires the numerical solution of the Schrödinger equation, which is often time consuming.

In this paper we propose an extension to a method which by construction is guaranteed to reach the desired final state [15–18]. This implies that numerical solution of the Schrödinger equation is not necessary. Thus, we need only assess the constraints on duration, resources and robustness of the protocol, which are all explicitly evaluated for each candidate solution. The starting point of the analysis is to define a time dependent trajectory ϕ₀(t) for our wave function reaching from the initial ϕ_i to the desired final state ϕ_f. We construct ϕ₀(t) as a non-degenerate eigenstate of a Hamiltonian, H₀(t). Subject to a slowly varying H₀(t), the time dependent solution of the Schrödinger equation will adiabatically follow ϕ₀(t) and end up in the final state with certainty [19]. The requirement of adiabaticity, however, only permits solutions of long duration, which may be useless due to the time resource available or due to the effect of decoherence mechanisms.

To speed up adiabatic protocols, we therefore invoke shortcut-to-adiabaticity (STA) [15–18]. STA constructs an explicitly modified Hamiltonian $\tilde{H}(t)$ which suppresses the non-adiabatic transitions and forces a quantum system to follow the eigenstates of H₀(t), thus maintaining perfect state transfer at finite evolution times. For the purpose of optimization, we now exploit the freedom in choosing the trajectory, i.e., H₀(t), such that the cost functional for the STA modified Hamiltonian is minimized. The method is visualized in figure 1, where each point in the H₀-space represents a realization of H₀(t), mapped by STA to a time dependent Hamiltonian $\tilde{H}(t)$ shown as a point in the $\tilde{H}$-space. The cost function is indicated as the vertical dimension in $\tilde{H}$-space. We search in H₀-space for the trajectory, leading to the physical control Hamiltonian in $\tilde{H}$-space that optimizes the costs. This type of approach has been applied in a two-level system, where the state trajectory was parametrized in order to find optimally robust solutions [20]. Previous work has also combined STA and optimal control theory for studying e.g. atom transport [21–24].

Zoom In Zoom Out Reset image size

We shall primarily be concerned with the duration and the energy requirements associated with application of strong control fields [25, 26]. We demonstrate our method by the application to population transfer in a three-level Λ-system. Recently, Du et al [27] reported a successful experimental application of STA to the STIRAP protocol in this system (see also [28]). By extending the pulse parametrization and optimizing in this space we find solutions which are almost twice as fast as the solution reported in [27], while still satisfying the experimental requirements.

The paper is outlined as follows: in section 2 we explain the optimization strategy in detail. In section 3 we apply the strategy to the three-level Λ-system. Finally, in section 4 we conclude the paper.

The goal of STA is to transport an initial state, ϕ_i, to a final state, ϕ_f along a chosen state trajectory ϕ₀(t) such that ϕ₀(t) is the instantaneous eigenstates of a reference Hamiltonian, H₀(t), and ϕ₀(0) = ϕ_i and ϕ₀(T) = ϕ_f. We note that time-evolution with the Hamiltonian H₀(t) does not follow ϕ₀(t) due to non-adiabatic transitions among the instantaneous eigenstates. STA suppresses these transitions by introducing a counter-diabatic term,

$$\begin{eqnarray}&&{H}_{\mathrm{cd}}={\rm{i}}\displaystyle \sum _{n}| {\partial }_{t}{\phi }_{n}\rangle \langle {\phi }_{n}| -\langle {\phi }_{n}| {\partial }_{t}{\phi }_{n}\rangle | {\phi }_{n}\rangle \langle {\phi }_{n}| ,\end{eqnarray} \tag{ 1 }$$

where $| {\phi }_{n}(t)\rangle$ are all the instantaneous eigenstates of H₀(t) with eigenvalues E_n(t). Regardless of the process duration, a system subject to the total, physical Hamiltonian,

$$\begin{eqnarray}&&\tilde{H}(t)={H}_{0}(t)+{H}_{\mathrm{cd}}(t),\end{eqnarray} \tag{ 2 }$$

experiences no transitions among the instantaneous eigenstates. If this Hamiltonian can be realized experimentally, the initial state ϕ_i evolves with certainty into ϕ_f at time T,

$$\begin{eqnarray}&&| \psi (T)\rangle ={{\rm{e}}}^{{\rm{i}}{\xi }_{0}(T)}| {\phi }_{0}(T)\rangle ={{\rm{e}}}^{{\rm{i}}{\xi }_{0}(T)}| {\phi }_{{\rm{f}}}\rangle ,\end{eqnarray} \tag{ 3 }$$

acquiring a dynamical and geometric phase

$$\begin{eqnarray}&&{\xi }_{0}(t)=-{\int }_{0}^{t}{E}_{0}(t){\rm{d}}{t}+{\rm{i}}{\int }_{0}^{t}\langle {\phi }_{0}(t)| {\partial }_{t}{\phi }_{0}(t)\rangle {\rm{d}}{t}.\end{eqnarray} \tag{ 4 }$$

There are many choices for time dependent H₀(t) which solve a given state transfer. We handle the experimental constraints by defining a cost functional on the physical Hamiltonian $\tilde{H}(t)$ and by optimizing H₀(t) such that the cost is minimized (see figure 1). A unique feature of this optimization strategy compared to traditional optimal control theory, is the guarantee of successful state transfer. Thus, the numerical optimization only focuses on the experimental constraints.

We study application of the optimization strategy to population transfer in a three-level Λ-system, with two lower states ($| 1\rangle$ and $| 3\rangle$) and an excited state ($| 2\rangle$), which has a very short life time. The two lower levels of such a system can, for example, be used to represent a qubit, while the excited state mediates the coupling between the qubit states [29]. We follow [27], where the authors realize the Λ-system in two hyperfine levels of the ground state and a short-lived optically exited state in ⁸⁷Rb atoms. The goal is to transfer population from $| 1\rangle$ to $| 3\rangle$, while avoiding population in $| 2\rangle$ due to its short life time. We assume that a direct coupling is unavailable. Instead, a stimulated Raman transition is used by applying two lasers coupling each of the ground states to the excited state. The adiabatic STIRAP protocol has been successfully applied to this problem [29, 30].

In the rotating wave approximation and under the two-photon resonance condition, the system Hamiltonian is given by (ℏ = 1)

$$\begin{eqnarray}{H}_{{\rm{\Lambda }}}(t)=\displaystyle \frac{1}{2}\left(\begin{array}{ccc}0 & {{\rm{\Omega }}}_{P}(t){{\rm{e}}}^{{\rm{i}}{\phi }_{L}} & 0\\ {{\rm{\Omega }}}_{P}(t){{\rm{e}}}^{-{\rm{i}}{\phi }_{L}} & 2{\rm{\Delta }} & {{\rm{\Omega }}}_{S}(t)\\ 0 & {{\rm{\Omega }}}_{S}(t) & 0\end{array}\right),\end{eqnarray} \tag{ 5 }$$

where Ω_P(t) and Ω_S(t) are the real Rabi frequencies of the pump and Stokes lasers, coupling $| 1\rangle$-$| 2\rangle$ and $| 3\rangle$-$| 2\rangle$, respectively. Δ is the one-photon detuning and ϕ_L is the fixed phase difference between the lasers. We consider the Λ-system under large one-photon detuning, Δ ≫ Ω_P(t), Ω_S(t), where the excited state can be adiabatically eliminated, such that the system reduces to an effective two-level system [31]. Here, the effective Hamiltonian is given by

$$\begin{eqnarray}{H}_{0}(t)=\displaystyle \frac{-1}{2}\left(\begin{array}{cc}{{\rm{\Delta }}}_{\mathrm{eff}} & {{\rm{\Omega }}}_{\mathrm{eff}}\\ {{\rm{\Omega }}}_{\mathrm{eff}} & -{{\rm{\Delta }}}_{\mathrm{eff}}\end{array}\right),\end{eqnarray} \tag{ 6 }$$

with

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{eff}}=\displaystyle \frac{{{\rm{\Omega }}}_{P}{{\rm{\Omega }}}_{S}}{2{\rm{\Delta }}},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\mathrm{eff}}=\displaystyle \frac{{{\rm{\Omega }}}_{P}^{2}-{{\rm{\Omega }}}_{S}^{2}}{4{\rm{\Delta }}}.\end{eqnarray} \tag{ 8 }$$

The eigenstates of the effective Hamiltonian are given by

$$\begin{eqnarray}&&| {a}_{0}\rangle =\cos \theta | 1\rangle -\sin \theta | 3\rangle ,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&| {a}_{-}\rangle =\sin \theta | 1\rangle +\cos \theta | 3\rangle ,\end{eqnarray} \tag{ 10 }$$

where the mixing angle is defined as $\theta =\arctan \tfrac{{{\rm{\Omega }}}_{P}}{{{\rm{\Omega }}}_{S}}$. Imposing the boundary conditions

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Omega }}}_{P}(0)}{{{\rm{\Omega }}}_{S}(0)}=0,\ \ \ \ \displaystyle \frac{{{\rm{\Omega }}}_{S}(T)}{{{\rm{\Omega }}}_{P}(T)}=0,\end{eqnarray} \tag{ 11 }$$

ensures that $| {a}_{0}(0)\rangle =| 1\rangle$ and $| {a}_{0}(T)\rangle =-| 3\rangle$. The reference pulses Ω_P(t) and Ω_S(t) thus define the trajectory ϕ_o(t). If the process duration is large enough, by the adiabatic theorem, the system follows the eigenstate, $| {a}_{0}(t)\rangle$, and the STIRAP protocol therefore realizes the state transfer.

3.1. STA for a STIRAP trajectory

We can use the STA formalism to follow the eigenstates of H₀(t) even when the system is driven more rapidly. To do this we calculate the counter-diabatic Hamiltonian, equation (1). The result is

$$\begin{eqnarray}{H}_{\mathrm{cd}}({t})=\left(\begin{array}{cc}0 & -{\rm{i}}{{\rm{\Omega }}}_{a}\\ -{\rm{i}}{{\rm{\Omega }}}_{a} & 0\end{array}\right),\end{eqnarray} \tag{ 12 }$$

where ${{\rm{\Omega }}}_{a}=\tfrac{{\dot{{\rm{\Omega }}}}_{P}{{\rm{\Omega }}}_{S}-{\dot{{\rm{\Omega }}}}_{S}{{\rm{\Omega }}}_{P}}{{{\rm{\Omega }}}_{P}^{2}+{{\rm{\Omega }}}_{S}^{2}}$.

Implementing $H(t)={H}_{0}(t)+{H}_{\mathrm{cd}}(t)$ requires temporal control of the relative phase between the Stokes and pump lasers, which can be circumvented by transforming to a frame, defined by the unitary transformation

$$\begin{eqnarray}U=\left(\begin{array}{cc}{{\rm{e}}}^{-{\rm{i}}\gamma (t)/2} & 0\\ 0 & {{\rm{e}}}^{{\rm{i}}\gamma (t)/2}\end{array}\right),\end{eqnarray} \tag{ 13 }$$

with $\gamma =\arctan \left(\tfrac{{{\rm{\Omega }}}_{a}}{{{\rm{\Omega }}}_{\mathrm{eff}}}\right)+{\phi }_{L}$ [27]. The resulting Hamiltonian

$$\begin{eqnarray}\tilde{H}(t)=\displaystyle \frac{-1}{2}\left(\begin{array}{cc}{\tilde{{\rm{\Delta }}}}_{\mathrm{eff}} & {\tilde{{\rm{\Omega }}}}_{\mathrm{eff}}\\ {\tilde{{\rm{\Omega }}}}_{\mathrm{eff}} & -{\tilde{{\rm{\Delta }}}}_{\mathrm{eff}}\end{array}\right),\end{eqnarray} \tag{ 14 }$$

can be implemented with a real Rabi frequency ${\tilde{{\rm{\Omega }}}}_{\mathrm{eff}}=\sqrt{{{\rm{\Omega }}}_{\mathrm{eff}}^{2}+{{\rm{\Omega }}}_{a}^{2}}$ and ${\tilde{{\rm{\Delta }}}}_{\mathrm{eff}}={{\rm{\Delta }}}_{\mathrm{eff}}+\dot{\gamma }$, and since the basis state populations are unaffected by this transformation, the new $\tilde{H}(t)$ also realizes the population transfer. We now look for the physical pulses, ${\tilde{{\rm{\Omega }}}}_{P}(t)$ and ${\tilde{{\rm{\Omega }}}}_{S}(t)$, applied to the original three-level system, that realize $\tilde{H}(t)$. We thus solve for the values of ${\tilde{{\rm{\Omega }}}}_{P}$ and ${\tilde{{\rm{\Omega }}}}_{S}$ that yield equations (7) and (8) with Ω_eff, Δ_eff replaced by ${\tilde{{\rm{\Omega }}}}_{\mathrm{eff}}$ and ${\tilde{{\rm{\Delta }}}}_{\mathrm{eff}}$. The result is

$$\begin{eqnarray}&&{\tilde{{\rm{\Omega }}}}_{P}(t)=\sqrt{2{\rm{\Delta }}(\sqrt{{\tilde{{\rm{\Delta }}}}_{\mathrm{eff}}^{2}+{\tilde{{\rm{\Omega }}}}_{\mathrm{eff}}^{2}}+{\tilde{{\rm{\Delta }}}}_{\mathrm{eff}})},\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\tilde{{\rm{\Omega }}}}_{S}(t)=\sqrt{2{\rm{\Delta }}(\sqrt{{\tilde{{\rm{\Delta }}}}_{\mathrm{eff}}^{2}+{\tilde{{\rm{\Omega }}}}_{\mathrm{eff}}^{2}}-{\tilde{{\rm{\Delta }}}}_{\mathrm{eff}})}.\end{eqnarray} \tag{ 16 }$$

For any choice of the reference pulses, Ω_P(t) and Ω_S(t), that fulfill the conditions, equation (11), we can calculate ${\tilde{{\rm{\Omega }}}}_{\mathrm{eff}}(t)$ and ${\tilde{{\rm{\Delta }}}}_{\mathrm{eff}}(t)$ and therefore the physical pulses ${\tilde{{\rm{\Omega }}}}_{P}(t)$ and ${\tilde{{\rm{\Omega }}}}_{S}(t)$. As long as the elimination of the excited state remains valid, subjecting the three-level atom to these physical pulses will yield the perfect transfer between the ground states in any finite time interval. Figure 1 illustrates how equations (15) and (16) represents the mapping between H₀-space and $\tilde{H}$-space,

3.2. Cost functional

For any choice of the reference pulses that fulfill the boundary conditions, STA provides the corresponding physical pulses through equations (15) and (16). However, the physical pulses might violate constraints set by the experiment. The constraints considered here include the peak intensity of the lasers and robustness against a scaling of the control parameters and are based on the experiment reported in [27].

The peak intensity of the laser can be included in the cost functional as the dimensionless quantity

$$\begin{eqnarray}&&{\tilde{{\rm{\Omega }}}}_{\mathrm{peak}}=\max \{{\tilde{{\rm{\Omega }}}}_{S/P}(t)\}/{{\rm{\Omega }}}_{0},\end{eqnarray} \tag{ 17 }$$

where we choose Ω₀ = 2π · 5 MHz to define the scale. Minimizing ${\tilde{{\rm{\Omega }}}}_{\mathrm{peak}}$ is equivalent to minimizing the peak intensity since ${I}_{S/P}\propto {\tilde{{\rm{\Omega }}}}_{S/P}^{2}$.

In experiments with many atoms, the spatial laser profile causes different atoms to experience different laser powers depending on their location. Effectively, this corresponds to a random scaling of the Rabi frequency, $\tilde{{\rm{\Omega }}}\to \epsilon \tilde{{\rm{\Omega }}}$, where  ≈ 1. In the adiabatic limit, this scaling does not alter the state transfer, but the values of the time dependent Rabi frequencies of the physical pulses are important when we apply the STA, and the scaling reduces the transfer fidelity, as shown in figure 2. We thus seek solutions which are robust against this perturbation. The sensitivity towards amplitude scaling can be quantified by perturbation theory [20]. To second order in  − 1, the correction to the transfer fidelity is found to be

$$\begin{eqnarray}&&{ \mathcal F }\approx 1-4q{(\epsilon -1)}^{2},\end{eqnarray} \tag{ 18 }$$

where we define the sensitivity,

$$\begin{eqnarray}&&q={\left|{\int }_{0}^{T}{{\rm{e}}}^{{\rm{i}}[{\xi }_{0}(t)-{\xi }_{-}(t)]}\langle {a}_{-}(t)| {U}^{\dagger }(t)\tilde{H}(t)U(t)| {a}_{0}(t)\rangle {\rm{d}}{t}\right|}^{2},\end{eqnarray} \tag{ 19 }$$

where U(t) is given by equation (13) and ξ_i is given by equation (4). By minimizing q we minimize the sensitivity to variations in intensity.

Zoom In Zoom Out Reset image size

To penalize solutions with large peak Rabi frequency, ${\tilde{{\rm{\Omega }}}}_{\mathrm{peak}}$, and large values of the sensitivity, q, we introduce a cost functional. We define our goal based on the peak value and sensitivity found in [27]. Here, two Gaussians are used for the reference pulses at a process duration of T = 0.4 ms, and for such Gaussians we have ${\tilde{{\rm{\Omega }}}}_{\mathrm{peak}}=1.14$ and q = 1.59. We seek to match these values at the lowest possible duration. We heuristically find that a cost functional defined as

$$\begin{eqnarray}&&C=\exp [10({\tilde{{\rm{\Omega }}}}_{\mathrm{peak}}-1.14)]+\exp [2(q-1.59)]\end{eqnarray} \tag{ 20 }$$

represents a balanced minimum of ${\tilde{{\rm{\Omega }}}}_{\mathrm{peak}}$ and q in accordance with our goal.

3.3. Parametrization of a family of reference pulses

As any two functions fulfilling the boundary conditions, equation (11), can be chosen for the reference pulses, it is difficult to search the entire H₀-space. In [27] the authors parameterize the reference pulses as partially overlapping Gaussians. We hence restrict the search to smooth and time symmetric solutions Ω_P(t) = Ω_S(T−t). We choose to define

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{P}(t)={{\rm{\Omega }}}_{0}f(t-T/2-T/10),\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{S}(t)={{\rm{\Omega }}}_{0}f(-t+T/2-T/10).\end{eqnarray} \tag{ 22 }$$

The parametrization in [27] is extended by choosing the parametrization function, f(t), as a sum of Gaussians,

$$\begin{eqnarray}&&f(t)=A\left({{\rm{e}}}^{-{t}^{2}/{(T/6)}^{2}}+\displaystyle \sum _{n=1}^{N}{a}_{n}{{\rm{e}}}^{-{(t-{t}_{n}^{0})}^{2}/{({w}_{n})}^{2}}\right),\end{eqnarray} \tag{ 23 }$$

where A is chosen such that max{f(t)} = 1. Here, the amplitude, offsets and widths ({a_n, t_n⁰, w_n}) are the control parameters. We found N = 4 to offer good solutions.

At this point the problem is reduced to finding the set of control parameters, {a_n, t_n⁰, w_n}, that minimizes the cost functional, C. This is done by locally optimizing several initial guesses, or seeds, for the parameters. The seeds are constructed by choosing random values in a suitable interval. An optimization routine is then employed to iteratively update the control parameters until C is locally minimized, as shown in figure 2(b). The optimal physical pulses are then constructed from the optimized control parameters. It is necessary to generate and optimize multiple seeds, since the optimization landscape often contains multiple local optima, as illustrated in figure 1.

The crab, grape and Krotov algorithms are widely used as local optimizers on multiple seeds in traditional quantum optimization approaches [5–7, 32], and have also been used to solve control problems in three-level systems [33]. However, for these algorithms the cost functional involves the transfer fidelity in addition to the terms from the experimental constraints. This makes the evaluation of the cost functional computationally expensive, as it requires the Schrödinger equation to be numerically solved. This is not required in our approach, as unit transfer fidelity is guaranteed by the STA formalism. Hence, the cost functional, C, can be quickly evaluated and optimized by a standard optimization routine.

3.4. Results

The result of the optimization is presented in figure 3(a). The triangles and squares mark ${\tilde{{\rm{\Omega }}}}_{\mathrm{peak}}$ and q, respectively, obtained by minimizing the cost functional at the given process duration. The dashed–dotted lines mark the target values of ${\tilde{{\rm{\Omega }}}}_{\mathrm{peak}}=1.14$ and q = 1.59, which are the values from [27] obtained using Gaussian reference pulses, $f(t)=\exp (-{t}^{2}/{(T/6)}^{2})$, with T = 0.4 ms. We search for the lowest process duration where we can find equivalent values of ${\tilde{{\rm{\Omega }}}}_{\mathrm{peak}}$ and q. Both ${\tilde{{\rm{\Omega }}}}_{\mathrm{peak}}$ and q decrease as the process duration increases and the lowest duration where our target is met is found to be T_eqv = 0.25 ms. That is, we find solutions that are nearly twice as fast compared to [27] without compromising energy consumption or robustness. The thick and thick-dashed lines show ${\tilde{{\rm{\Omega }}}}_{\mathrm{peak}}$ and q obtained using Gaussian reference pulses for process durations lower than T = 0.4 ms. The optimized pulses perform significantly better at all durations compared to the Gaussian pulse.

Zoom In Zoom Out Reset image size

In figure 3(b) we plot cost values before and after optimization for 3000 seeds. The figure shows that seeds with low-cost yield good results more often than seeds with high cost. A seed is taken to be successful if the cost after optimization is below the criterion value, log C_success = 0.8. The probability for finding a successful seed for each bin suggests that only seeds with low cost should be optimized. Such low-cost seeds are rare (see insert in figure 3(b)). However, generating and evaluating a seed is computationally inexpensive, and to produce the data points in figure 3(a) we chose to generate 5 million seeds and optimized only the 1000 of them with lowest cost.

We have proposed an optimization strategy for solving state-to-state quantum control problems. Our strategy combines the STA formalism and minimization of a cost functional incorporating resource requirements and a perturbative expression for the robustness. Unlike traditional quantum optimal control algorithms our cost functional does not include the transfer fidelity and is therefore computationally inexpensive to evaluate and optimize. We have demonstrated the capability of our strategy on a control problem in the three-level Λ-system. Here we find solutions that are almost twice as fast as the solution reported in [27], while still obeying experimental constraints. The calculation leading to the perturbative expression for the robustness, equation (18), can be carried out for any perturbation, and, in principle, to any order [20]. This makes our approach especially well-suited for finding solutions that are robust against perturbations. Traditional optimal control algorithms do not immediately offer these solutions, since unit fidelity is required to construct such perturbative expressions.

Our strategy can be applied and is efficient for any system where STA is applicable, i.e., in systems where the adiabatic eigenstates and the counter-diabatic driving Hamiltonian can be easily evaluated and applied. A related method makes use of the so-called Lewis–Riesenfeld invariant to determine the time dependence of Hamiltonian parameters that lead with certainty to a desired final state [17]. In analogy with the mapping in our figure 1, we may use the parametrization of the invariant I(t) as our starting point in an 'I-space' and as long as the associated H(t) is easy to obtain, we may efficiently try several candidate I(t) and optimize for robustness and energy costs. Examples of such systems include two- and three-level systems [21, 22, 27, 34], atoms in harmonic traps [17, 23, 35], quantum many-body systems [36, 37] and quantum heat engines [38–40].

We believe that our approach is a valuable addition to the arsenal of quantum optimal control algorithms, and especially for control problems that require robust solutions.

This work was supported by the European Research Council, the Lundbeck Foundation and the Villum Foundation.