# Invariant-based optimal composite stimulated Raman exact passage

Tian-Niu Xu, Kaipeng Liu, Xi Chen, Stephane Guerin

► To cite this version:

Tian-Niu Xu, Kaipeng Liu, Xi Chen, Stephane Guerin. Invariant-based optimal composite stimulated Raman exact passage. Journal of Physics B: Atomic, Molecular and Optical Physics, 2019, 52 (23), pp.235501. 10.1088/1361-6455/ab49a9. hal-02427337

HAL Id: hal-02427337

https://hal.science/hal-02427337v1

Submitted on 3 Jan 2020

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Invariant-based optimal composite stimulated Raman exact passage

Tian-Niu Xu

International Center of Quantum Artificial Intelligence for Science and Technology (QuArtist) and Department of Physics, Shanghai University, 200444 Shanghai, China

Kaipeng Liu

Laboratoire Interdisciplinaire Carnot de Bourgogne, UMR CNRS 6303, BP 47870, F-21078 Dijon, France

Xi Chen

E-mail: xchen@shu.edu.cn

International Center of Quantum Artificial Intelligence for Science and Technology (QuArtist) and Department of Physics, Shanghai University, 200444 Shanghai, China
Department of Physical Chemistry, University of the Basque Country, 48080 Bilbao, Spain

Stéphane Guérin

E-mail: sguerin@u-bourgogne.fr

Laboratoire Interdisciplinaire Carnot de Bourgogne, UMR CNRS 6303, BP 47870, F-21078 Dijon, France

**Abstract.** Shortcuts to adiabaticity have been put forward for accelerating slow adiabatic passages in various quantum systems with tremendous applications for performing quantum information processing tasks. In this paper, we propose a hybrid protocol to achieve stimulated Raman exact passage (STIREP) by combing invariant-based inverse engineering, optimal control, and composite pulse approaches. We first derive the general solution and their corresponding pulse shapes by invariant-based inverse engineering. Counterintuitive and optimal (intuitive) pulse sequences are formulated in this context and incorporated into composite sequences. Such composite stimulated Raman exact passage not only features robustness against the fluctuation of laser intensity, but also reduces the operation time and energy cost.

# 1. INTRODUCTION

Coherent manipulation and preparation for quantum states with high fidelity are requisite with many applications ranging from quantum information processing to control of chemical interaction [1–5]. Among all techniques, resonant pulses, rapid
adiabatic passage (RAP), stimulated Raman adiabatic passage (STIRAP) and their variants provide the full population transfer in quantum two- or three-level systems with different advantage and disadvantage [6–11]. Resonant pulses are usually fast but sensitive to the parameter variations, while adiabatic passages are in principle robust but slow and inaccurate. To remedy this, different approaches, as shortcuts to adiabaticity (STA), optimal control and composite pulses [12–17], have been proposed to combine the best of the two worlds. Some techniques, sharing the concept of STA [18, 19] include the counter-diabatic driving, fast-forward scaling, inverse engineering based on Lewis-Riesenfeld dynamical invariant, and single-shot shaped pulses [20, 21], also combined with optimal control [22].

Quantum two- or three-level systems are ubiquitous in atom, molecular, optical physics and solid-state devices. Particular systems might require specific techniques for population transfer in presence of different types of inevitable systematic errors or noise. These processes have been widely studied, and there are substantial theoretical discussions and experimental demonstrations of its efficiency for producing complete population transfer. In STIRAP, the Stokes pulse drives the transition between the initially unpopulated levels $|2\rangle$ and $|3\rangle$. It coherently prepares the system and precedes the overlapped pump pulse partially (driving the transition between the initially populated level $|1\rangle$ and the intermediate level $|2\rangle$). Such a sequence is termed counterintuitive. Population transfer from level $|1\rangle$ to level $|3\rangle$ is approximately achieved by adiabatic passage along a single instantaneous dressed (dark) state. Transient population arising in the intermediate state $|2\rangle$ during the adiabatic transfer is small or even negligible. This technique has quickly attracted considerable attention as a promising control tool for quantum information processing (QIP) [23–26]. QIP demands very high fidelity of operations with the permissible error at most $10^{-4}$, which is difficult to achieve with the standard STIRAP. Because, due to its adiabatic nature, it approaches unit efficiency only asymptotically as the temporal pulse areas increase. For usual Gaussian shape pulse, the necessary area for the $10^{-4}$ benchmark is so large that it may break various restrictions in a real experiment. Nevertheless, the adiabatic process is slow and therefore likely to be affected by decoherence, dissipation or systematic error and do not lead to an exact transfer. In order to solve these problems, several scenarios have been proposed to speed up STIRAP with an ultrahigh fidelity [27–32].

Motivated by the current experimental demonstration of composite Stimulated Raman adiabatic passage in a rare-earth doped solid [33], we hybridize the composite method with pulses designed by invariant-based inverse engineering, including exact counterintuitive or optimal (intuitive) sequences, to derive an exact, fast and robust technique. The inverse engineering approach based on the Lewis-Riesenfeld theory allows one to define dynamical invariant and eigenmodes on which the dynamics can be decomposed [34]. We connect the exact single-mode dynamics to the (approximate) adiabatic passage along a dark state resulting from counterintuitive pulse sequence (STIRAP). We derive from a multi-mode decomposition with a parametrization linear in time a dynamics featuring a counterintuitive pulse sequence, which is shown to be
<figure><img src="image_2.png" /><figcaption>**Figure 1.** (a) Schematic view of Λ-type three-level system for resonant STIRAP. (b-d) The Rabi frequencies Ω<sub>p,s</sub> present the pump (dash-dotted blue line) and Stokes (solid red line) pulses for different composite schemes of three pulse pairs: Gaussian, area-optimal and linear multimode pulses. The latter two sequences define examples of composite intuitive and counterintuitive STIREP, respectively.</figcaption></figure>

the opposite sequence with the same pulse shapes of the area-optimal sequence. This
linear multi-mode driving and the area-optimal pulses are implemented in composite
sequences. We find that the invariant-based composite STIREP takes less and less time
compared to composite STIRAP for more and more sequences.

In Sec. 2, we define the system, and review the technique to design pulses by the inverse engineering approach based on the Lewis-Riesenfeld dynamical invariant, for completeness. In Sec. 3, we derive explicit solutions for pulses that are optimal with respect to energy (or equivalently to time for bounded pulses) or to pulse area. Section 4 is devoted to the comparison of three composite protocols. Discussion in Sec. IV concludes the paper.

2. INVARIANT DYNAMICS

2.1. *Model and Hamiltonian*

The traditional STIRAP process is controlled by two counterintuitive Gaussian pulses,
as shown in Fig. 1(a). The pump pulse connects the ground initial state |1⟩ and excited
state |2⟩ at frequency ω<sub>p</sub>. The Stokes pulse links the intermediate state |2⟩ and final
target state |3⟩ at frequency ω<sub>s</sub>. The energies ħω<sub>1</sub>, ħω<sub>2</sub> and ħω<sub>3</sub> correspond to the
discrete states |1⟩, |2⟩ and |3⟩, respectively. The corresponding Hamiltonian, in the
rotating-wave approximation at resonance, ω<sub>1</sub> + ω<sub>p</sub> = ω<sub>2</sub> = ω<sub>3</sub> + ω<sub>s</sub>, reads

$$
H_0(t) = \frac{\hbar}{2} \begin{bmatrix} 0 & \Omega_p(t) & 0 \\ \Omega_p(t) & 0 & \Omega_s(t) \\ 0 & \Omega_s(t) & 0 \end{bmatrix}, \qquad (1)
$$
where $\Omega_p(t)$ and $\Omega_s(t)$ denote the pump and Stokes Rabi frequencies, respectively, considered real functions.

## 2.2. The invariant

The Lewis-Riesenfeld theory [34] allows one to expand exactly the solution $|\psi(t)\rangle$ on dynamical modes $|\phi_n(t)\rangle$, which are defined as eigenvectors of the dynamical invariant $I(t)$, $I(t)|\phi_n(t)\rangle = \lambda_n|\phi_n(t)\rangle$, defined as:

$$ \langle\psi(t)|I(t)|\psi(t)\rangle = \text{const.}, \qquad (2) $$

or equivalently as

$$ dI/dt \equiv \partial I(t)/\partial t + (1/i\hbar)[I(t), H(t)] = 0. \qquad (3) $$

Its eigenvalues are real and time independent and the eigenvectors form a complete orthonormal basis, which allows an expansion of the solution as a superposition of the dynamical modes:

$$ |\psi(t)\rangle = \sum_n C_n e^{i\xi_n} |\phi_n(t)\rangle \qquad (4) $$

with time independent coefficients $C_n$ and the Lewis-Riesenfeld phases (determined from the initial time $t_i$):

$$ \xi_n(t) = \frac{1}{\hbar} \int_{t_i}^{t} \langle \phi_n(t') | i\hbar\partial/\partial t' - H(t') | \phi_n(t') \rangle dt'. \qquad (5) $$

The Hamiltonian (1) features SU(2) dynamical symmetry:

$$ H = \frac{\hbar}{2} (\Omega_p(t) K_1 + \Omega_s(t) K_2), \qquad (6) $$

where $K_1$, $K_2$, and $K_3$ are the spin 1 angular-momentum operators:

$$ K_1 = \begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}, K_2 = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{bmatrix}, K_3 = \begin{bmatrix} 0 & 0 & -i \\ 0 & 0 & 0 \\ i & 0 & 0 \end{bmatrix} \qquad (7) $$

satisfying the commutation relations

$$ [K_1, K_2] = iK_3, \quad [K_2, K_3] = iK_1, \quad [K_3, K_1] = iK_2. \qquad (8) $$

The dynamical invariant is also taken as an element of this algebra

$$ I(t) = \sum_i \alpha_i(t) K_i, \qquad (9) $$

from which, using the definition of the invariant and denoting $|\alpha(t)\rangle = [\alpha_1(t), \alpha_2(t), \alpha_3(t)]^T$ (where the superscript T denotes the transpose), we derive

$$ i\hbar \frac{d}{dt} |\alpha(t)\rangle = B(t) |\alpha(t)\rangle \qquad (10) $$

with the Hermitian matrix

$$ B(t) = \frac{\hbar}{2} \begin{bmatrix} 0 & 0 & i\Omega_s(t) \\ 0 & 0 & -i\Omega_p(t) \\ -i\Omega_s(t) & i\Omega_p(t) & 0 \end{bmatrix} \qquad (11) $$
and (from the unitarity of the propagator):

$$ \langle \alpha(t) | \alpha(t) \rangle = \sum_i \alpha_i^2(t) = \text{const} = \langle \alpha(t_i) | \alpha(t_i) \rangle. \quad (12) $$

To parameterize the three angles $\alpha_i$'s of the expansion (9), we need two (time-dependent) angles $\varphi, \theta$ and a constant $\alpha_0$ in order to satisfy (12), $\alpha_1^2 + \alpha_2^2 + \alpha_3^2 = \alpha_0^2$:

$$ \alpha_1 = \alpha_0 \sin \theta \cos \varphi, \quad \alpha_2 = \alpha_0 \cos \theta \cos \varphi, \quad \alpha_3 = \alpha_0 \sin \varphi. \quad (13) $$

Equation (10) allows one to link the pulses to the angles:

$$ \Omega_p/2 = \dot{\varphi} \cos \theta + \dot{\theta} \sin \theta \cotan \varphi, \quad \Omega_s/2 = -\dot{\varphi} \sin \theta + \dot{\theta} \cos \theta \cotan \varphi. \quad (14) $$

The invariant $I(t)$ reads with the parametrization (13):

$$ I(t) = \alpha_0 \begin{bmatrix} 0 & \sin \theta \cos \varphi & -i \sin \varphi \\ \sin \theta \cos \varphi & 0 & \cos \theta \cos \varphi \\ i \sin \varphi & \cos \theta \cos \varphi & 0 \end{bmatrix}, \quad (15) $$

giving the orthonormal dynamical modes and the associated phases (for the eigenvalues $\lambda_0 = 0$ and $\lambda_{\pm} = \pm \alpha_0$, respectively)

$$ |\phi_0(t)\rangle = \begin{bmatrix} \cos \varphi \cos \theta \\ -i \sin \varphi \\ -\cos \varphi \sin \theta \end{bmatrix}, \quad \dot{\xi}_0 = 0, \quad (16) $$

$$ |\phi_{\pm}(t)\rangle = \frac{1}{\sqrt{2}} \begin{bmatrix} \cos \theta \sin \varphi \pm i \sin \theta \\ i \cos \varphi \\ -\sin \theta \sin \varphi \pm i \cos \theta \end{bmatrix}, \quad \dot{\xi}_{\pm} = \mp \frac{\dot{\theta}}{\sin \varphi}. \quad (17) $$

## 2.3. Single-mode driving: From adiabatic to exact passage

Adiabatic passage assumes a slow evolution of the parameters such that the dynamics projects approximately along a single adiabatic (dark) state,

$$ |\phi_D(t)\rangle = \begin{bmatrix} \cos \theta \\ 0 \\ -\sin \theta \end{bmatrix}, \quad (18) $$

where $\theta$ is the mixing angle, $\tan \theta = \Omega_p/\Omega_s$, which connects the initial with the final states. We notice that $|\phi_D(t)\rangle$ is a particular case of the single-mode state (16) with $\varphi = 0$. The initial state $|1\rangle$ can in principle connect to the dark state in a counterintuitive sequence, i.e. when the Stokes precedes the pump pulse (i.e. $\theta$ going from 0 to $\pi/2$), but the dynamics cannot permanently project exclusively on it since the dark state is connected to the other bright states through the non-adiabatic coupling $\dot{\theta}$ in adiabatic frame. This non-adiabatic coupling can be made small for a slow variation of the parameters, which defines the adiabatic passage. The Lewis-Riesenfeld theory allows one to derive an alternative dynamical basis, on which one can project at all times the dynamics exactly from the initial condition.
One can define a single-invariant-mode driving, or in short a single-mode driving, by imposing at the initial time $\theta(t_i) = 0$ and $\varphi(t_i) = 0$, which then drives the dynamics into the target state, by imposing the final boundary $\theta(t_f) = \pi/2$ and $\varphi(t_f) = 0$ at time $t_f$. These boundaries have to be taken into account in the definition of the pulses (14) to ensure the desired dynamics. Imposing the complete dynamics occurring along the dark state would imply $\varphi = 0$ at all time, which would correspond to infinite pulse amplitudes from Eqs. (14). Single-mode dynamics thus needs $\varphi \neq 0$ during the dynamics. This single-mode dynamics can be interpreted as an exact passage as opposed to the adiabatic (i.e. approximate) passage along a single adiabatic state. In addition to producing an exact passage, such dynamics can also feature certain optimization criteria or specific properties as studied below by an appropriate time parametrization of the angles $\theta(t)$ and $\varphi(t)$.

## 2.4. Multi-mode driving

One can more generally derive multi-mode solutions [29]. We analyze for instance a three-mode solution, still using $\theta(t_i) = 0$, and denoting the initial boundary $\varphi(t_i) = \varepsilon$, according to (4):

$$|\psi(t)\rangle = c_0|\phi_0(t)\rangle + c_+e^{i\xi_+(t)}|\phi_+(t)\rangle + c_-e^{i\xi_-(t)}|\phi_-(t)\rangle \quad (19)$$

with

$$c_0 = \langle\phi_0(t_i)|1\rangle = \cos\varepsilon, \quad (20)$$

$$c_+ = \langle\phi_+(t_i)|1\rangle = \sin\varepsilon/\sqrt{2}, \quad (21)$$

$$c_- = \langle\phi_-(t_i)|1\rangle = \sin\varepsilon/\sqrt{2}, \quad (22)$$

and the Lewis-Riesenfeld phase $\xi_0 = 0$. For simplicity we consider a constant $\varphi$ and the linear parametrization with $t \in [t_i \equiv 0, t_f \equiv T]$ and $\theta(T) = \pi/2$:

$$\varphi(t) = \varepsilon, \quad \theta(t) = \frac{\pi}{2T}t, \quad (23)$$

giving the phase $\xi_\pm(t) = \mp\xi(t)$ with $\xi(t) = \theta(t)/\sin\varepsilon$ and at final time:

$$\langle 3|\psi(T)\rangle = -(\cos^2\varepsilon + \cos\xi(T)\sin^2\varepsilon). \quad (24)$$

The complete population transfer is then finally achieved when $\cos\xi(T) = \cos(\frac{\pi}{2\sin\varepsilon}) = 1$, i.e.

$$\sin\varepsilon = \frac{1}{4k}, \quad k = \pm 1, \pm 2, \pm 3, \dots \quad (25)$$

The resulting pulses, derived from (14),

$$\Omega_p(t) = \frac{\pi}{T}\sqrt{16k^2 - 1} \sin\left(\frac{\pi}{2T}t\right), \quad \Omega_s(t) = \frac{\pi}{T}\sqrt{16k^2 - 1} \cos\left(\frac{\pi}{2T}t\right) \quad (26)$$

feature a counterintuitive sequence. The smallest pulses energy required for the complete transfer in this case is for $k = \pm 1$. We refer this dynamics to as a linear multi-mode driving.
Alternatively, for a given peak Rabi frequency $\Omega_0$, the fields (for $k = \pm 1$) can be written as

$$ \Omega_p(t) = \Omega_0 \sin \left( \frac{\pi t}{2T_{\text{multi}}} \right), \quad \Omega_s(t) = \Omega_0 \cos \left( \frac{\pi t}{2T_{\text{multi}}} \right), \qquad (27) $$

with the transfer time

$$ T_{\text{multi}} = \sqrt{15\pi}/\Omega_0. \qquad (28) $$

One can emphasize that the derived sequence “sin / cos” (for the pump and Stokes, respectively) is counterintuitive, and is the opposite sequence with the same pulse shapes of the optimal pulses (with respect to energy) “cos / sin” derived below in Section 3. This is of interest because counterintuitive sequences are known as featuring a relatively low transient population in the excited state, which is usually a lossy state, and a certain degree of robustness.

As shown below numerically, such multi-mode driving provides a fast and robust protocol. It can be explained by a combined resonant and adiabatic process, see Appendix A.

### 3. OPTIMAL PULSES

In the resonantly driven three-level system (1), the optimal pulses leading to the transfer have been already addressed [35]. Here we rederive them adding explicit formula. We denote $u_p \equiv \Omega_p/2$, $u_s \equiv \Omega_s/2$, $t_i \equiv 0$ and $t_f \equiv T$, the total duration of the process, i.e. the transfer time. Optimality can be analyzed in term of pulse energy (defined through $u_p$ and $u_s$ for simplicity)

$$ \mathcal{E} = \hbar \int_0^T [u_p^2(t) + u_s^2(t)] dt, \qquad (29) $$

or in term of generalized pulse area

$$ \mathcal{A} = \int_0^T \sqrt{u_p^2(t) + u_s^2(t)} dt. \qquad (30) $$

By using Pontryagin's maximum principle, see Appendix B, we obtain for the energy-optimal pulses

$$ \Omega_p(t) = \frac{\sqrt{3\pi}}{T} \cos \left( \frac{\pi t}{2T} \right), \quad \Omega_s(t) = \frac{\sqrt{3\pi}}{T} \sin \left( \frac{\pi t}{2T} \right), \qquad (31) $$

satisfying the boundary $\Omega_p^2 + \Omega_s^2 = \Omega_0^2$, with transfer time

$$ \mathcal{E} = \hbar \frac{3\pi^2}{4T}. \qquad (32) $$

One can remark that minimizing the cost $\mathcal{E}$ is equivalent to minimizing the transfer time for bounded controls, $u_p^2 + u_s^2 \le \text{const}$. Denoting $\Omega_0$ the peak of the pump and Stokes Rabi frequency,

$$ \Omega_p(t) = \Omega_0 \cos \left( \frac{\pi t}{2T_{\text{min}}} \right), \quad \Omega_s(t) = \Omega_0 \sin \left( \frac{\pi t}{2T_{\text{min}}} \right), \qquad (33) $$
we get for the optimal time:

$$T_{\min} = \sqrt{3\pi}/\Omega_0, \qquad (34)$$

which is $\sqrt{5}$ times smaller than the linear multi-mode transfer time. These pulses (33) represent time-optimal pulses.

As emphasized in Sec 2.4, the derived intuitive sequence cos / sin is the opposite sequence with the same pulse shapes of the linear multi-mode driving.

If one considers the minimization of the effective pulse area (30), defining the area-optimal pulses, one can remark that this quantity is time-reparametrization invariant as it leads to the same trajectory in the $(\theta, \varphi)$ space for any time-reparametrization, which is given by Eq. (B.15). The energy-optimal pulses (B.25) feature an example of area-optimal pulse with which one can determine the resulting optimal area:

$$\mathcal{A} = \frac{\sqrt{3\pi}}{2}. \qquad (35)$$

Another convenient parametrization, symmetric around $t = T/2$ (i.e. such that $\varphi(T - t) = \varphi(t)$), such that $\varphi(T/2) = \pi/3$, is given by

$$\varphi(t) = \begin{cases} 0 & \text{if } t \le 0 \text{ or } t \ge T \\ \frac{\pi e^{-\alpha t(t-T)/T^2} - 1}{3 e^{\alpha/4} - 1} & \text{otherwise.} \end{cases} \qquad (36)$$

Using the (dimensionless) quantity $\alpha$ such that $\alpha \gg 1$ allows to start and end the process with pulses close to zero. One can choose for instance $\alpha = 36$. An alternative parametrization, leading to smooth pulses, which start and end exactly at zero, is given by

$$\varphi(t) = \begin{cases} 0 & \text{if } t \le 0 \text{ or } t \ge T \\ \frac{\pi}{3} e^{4\beta} \exp\left(\frac{\beta T^2}{t(t-T)}\right) & \text{otherwise.} \end{cases} \qquad (37)$$

The (dimensionless) quantity $\beta$ measures the width of the pulse. The value $\beta = 2$ leads to pulses close to the ones obtained with parametrization (36) and $\alpha = 36$.

For comparison, the completely overlapping fields $\Omega_p(t) = \Omega_s(t) = \Omega_0/\sqrt{2}$, also satisfying the boundary $\Omega_p^2 + \Omega_s^2 = \Omega_0^2$, leads to complete population transfer, for the duration

$$T_{\pi} = 2\pi/\Omega_{0}, \qquad (38)$$

with the pulse energy $\mathcal{E} = \hbar\pi^2/T_{\pi}$ and pulse area $\mathcal{A} = \pi$. Alternatively, the successive pulses (without overlap), $\Omega_p(t) = \Omega_0/\sqrt{2}$ (for $0 \le t \le T/2$, 0 otherwise) and $\Omega_s(t) = \Omega_0/\sqrt{2}$ (for $T/2 \le t \le T$, 0 otherwise) also leads to the complete transfer for the same pulse area $\mathcal{A} = \pi$, with the transfer time $T = 2\sqrt{2\pi}/\Omega_0$. As a consequence, one can remark that all these pulses of small pulse area ($\pi$ or less for the optimal one) feature intuitive sequences (or fully overlapping pulses), with a relatively large transient population in the excited state
## 4. COMPOSITE PROTOCOLS

The implementation of the previously defined optimal pulses induces a dynamics sensitive to systematic errors. Composite pulses are designed to solve this issue. They use a sequence of pulses with appropriate phases, $\Omega_p(t) \to \Omega_p(t)e^{i\mu}$ and $\Omega_s(t) \to \Omega_s(t)e^{i\nu}$, to control the dynamics in a robust way [16, 17]. These relative phases convert the propagator $U(t)$ of the STIRAP dynamics associated to the Hamiltonian (1) into

$$ U_{\mu,\nu} = \psi U \psi^*, \psi = \begin{bmatrix} e^{i\mu} & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & e^{-i\nu} \end{bmatrix}. \qquad (39) $$

The composite STIRAP derived in [17] implements a high order fidelity, which is robust to variation of systematic errors, with universal phases, i.e. independent of the specific pulse shape, delay, and area. The dynamics with a sequence of N STIRAP's (with N an odd number) is represented by the propagator

$$ U^{(N)}(t) = U_{k_N, q_N} \tilde{U}_{\mu_{N-1}, \nu_{N-1}} \cdots U_{\mu_3, \nu_3} \tilde{U}_{\mu_2, \nu_2} U_{\mu_1, \nu_1}, \qquad (40) $$

where $\tilde{U}$ is a backward STIRAP propagator from state $|3\rangle$ to state $|1\rangle$, i.e. with the exchange of the order of the pump and Stokes pulses. In this section we compare the implementation of the composite STIRAP for three different pulse shapes: the Gaussian, area-optimal and linear multimode pulses, as illustrated in Fig. 1 (b-d).

**Protocol 1.** The Gaussian pulses are the most widely used form in STIRAP, see Fig. 1(b). They can be written as

$$ \Omega_p(t) = \Omega_0 e^{-(t-\frac{\tau}{2})^2/\sigma^2}, \quad \Omega_s(t) = \Omega_0 e^{-(t+\frac{\tau}{2})^2/\sigma^2}, \qquad (41) $$

where $\tau$ is the delay, the $\sigma$ is the pulse width and $\Omega_0$ is the peak amplitude of both Rabi frequencies. A positive delay defines a so-called counterintuitive sequence, since the Stokes precedes the pump pulse. We emphasize that such counterintuitive sequence is of importance compared to an intuitive sequence when the upper state is lossy, since it leads to much less transient population in this state. The fidelities corresponding to the change of the pulse amplitude from $\Omega_0$ to $\Omega_0(1 + \eta)$, via a systematic error $\eta$, is shown in Fig. 2 (solid red lines). We can see that the fidelities get better with respect to the systematic error $\eta$ for a larger sequence number $N$. However, the Gaussian pulses feature a large area and a large energy compared to the other invariant-based and optimal pulses as analyzed below.

**Protocol 2.** We use the parametrization (36) with the geodesic (B.20), (B.23), designed to minimize the pulse area (30). Setting the same peak Rabi frequency as for the Gaussian pulses of Protocol 1, we obtain the transfer time $T_{\text{area}} \approx 2.8\mu$s for the peak Rabi frequency $\Omega_0 = 4$MHz. With the area-optimal pulses shown in Fig. 1(c), the fidelities versus the variation of the systematic errors for different $N$ are shown in Fig. 2 (dotted black line).

**Protocol 3.** In this protocol, we consider the counterintuitive linear multimode driving pulse sequence (26) with $k=1$, i.e. of minimum energy, see Fig. 1(d). The
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Fidelity versus the variation of pump and Stokes pulses from $\Omega_{s,p}$ to $\Omega_{s,p}(1 + \eta)$ for Gaussian (solid red line), area-optimal (dotted black line) and linear-multimode (dashed blue line) pulses, for a single pulse pair ($N = 1$) and the sequences $N = 3, N = 5$. Parameters are as follows for: $\tau = 2\mu s$, $\sigma = 2\mu s$ (leading to a transfer time for Gaussian pulses that we can estimate as $T_{Gaussian} = 16\mu s$), $T_{area} = 2.8\mu s$, $T_{multi} = 3\mu s$, $\Omega_0 = 4$ MHz.</figcaption></figure>

energy associated to these pulses is $\mathcal{E} = 15\pi^2\hbar/(4T)$, for a given interaction time $T$. We obtain $T_{multi} \approx 3\mu s$ for $\Omega_0 = 4$MHz. The fidelities versus variation of systematic errors are shown in Fig. 2 (dashed blue line). Even though, as expected, the robustness of area-optimal pulses is worse than the other two protocols for a single pulse, it becomes comparable already for a sequence number $N = 5$.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** (a) Time transfer and (b) energy (denoted $\bar{E} \equiv \mathcal{E}$) for Gaussian (solid red line), area-optimal (dotted black line) and linear multimode pulses (dashed blue line) with ($N = 1, 2, \dots, 7$).</figcaption></figure>

**5. DISCUSSION AND CONCLUSION**

We can compare the transfer time and the energy of the studied protocols for area-optimal, Gaussian and linear multimode pulses at the same peak Rabi frequency (see Fig. 3). They are all robust to systematic errors and feature all a high fidelity that can satisfy the quantum information processing constraint [30].

The transfer times of composite STIRAP and composite STIREP are shown in Fig. 3(a) for different sequence numbers. We can conclude that the area-optimal pulses and linear multimode pulses, even if they are not designed for optimizing the time transfer, are much shorter than the Gaussian shape pulses. *The more sequence we use in the dynamics, the more time we can save.* The linear multimode pulses are only a little longer than the area-optimal pulses. We remark that, obviously, we could save even more time by using the time-optimal pulses (33). We can notice that the area-optimal pulses need less energy than the linear multimode pulses, see Fig. 3(b), as expected for intuitive compared to counterintuitive sequences. It means we can combine the area-optimal pulses and composite pulses into a hybrid protocol which features pulses starting and ending at zero, a short transfer time and a low energy cost (but at the price of high transient population in the upper state, see below).

STIRAP was initially developed for controlling atomic and molecular systems, in
which the dark state (18) plays a major role to transfer population over a dissipative intermediate state, rendering the dynamics insensitive to the loss from it (in the adiabatic limit). In order to quantify it, we calculate the time area of transient population in the excited state as $P_e = \int_0^T \sin^2 \varphi(t) dt$ for a single pulse pair. Such area $P_e \approx 0.264$ for the (counterintuitive) linear multimode driving pulse is twice larger than the one, 0.138, for (counterintuitive) Gaussian pulse, but twice lower than the one, 0.494, for the (intuitive) area-optimal protocol. Therefore, the linear multimode pulses should find applications in quantum information processing when the upper state is relatively weakly lossy or the process times are sufficiently short.

We have used single- and two-photon resonances in the derived hybrid protocols since they are based on the search of optimal pulse sequences already for a single shot. Detuned excitation corresponding to lower coupling strength would not lead to optimal sequences. In a system featuring large inhomogeneous broadening as in rare-earth doped solids, the use of composite detuned sequences has been shown to be more efficient and more robust [33] since resonant sequence is sensitive to single-photon inhomogeneous broadening. It would be of interest to reconsider the problem of optimality in $\Lambda$ systems including this additional effect (through a Lindblad equation with the density matrix formalism) and a single-photon detuning.

Last but not least, recent experiments on STA for speeding up STIRAP demonstrate that our hybrid protocols should be implementable for fast and robust state preparation and also for producing quantum gates with various systems such as cold atoms, superconducting circuits and solid-state spins [37–39].

## Acknowledgments

This work was supported by NSFC (11474193), SMSTC (18010500400 and 18ZR1415500), and the Program for Eastern Scholar. X.C. also thanks the Ramón y Cajal grant (RYC-2017-22482). L.K. acknowledges support from the CSC (China Scholarship Council). S.G. acknowledges additional support from the Investissements d'Avenir programs, project ISITE-BFC /IQUINS (ANR-15-IDEX-03), the EUR-EIPHI Graduate School (17-EURE-0002) and from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 765075 (LIMQUET).

## Appendix A. Connection between multi-mode driving and resonant pulses

In this appendix, one interprets the *linear multi-mode driving*, which provides fast and robust state transfer, in term of a resonantly driven spin-1 particle. To this end, we transform the Hamiltonian (1) into the adiabatic frame

$$H_{ad} = A^{\dagger}(t)H(t)A(t) - i\hbar A^{\dagger}(t)\frac{\partial}{\partial t}A(t), \qquad (A.1)$$
where

$$A(t) = \begin{pmatrix} \sin \theta / \sqrt{2} & \cos \theta & \sin \theta / \sqrt{2} \\ 1 / \sqrt{2} & 0 & -1 / \sqrt{2} \\ \cos \theta / \sqrt{2} & -\sin \theta & \cos \theta / \sqrt{2} \end{pmatrix}.$$

The transformed Hamiltonian becomes

$$H_{ad} = \Omega J_z + \dot{\theta} J_y, \qquad (A.2)$$

where $\Omega = \sqrt{\Omega_s^2(t) + \Omega_p^2(t)}/2$, and the matrix elements for spin-1 particle are

$$J_x = \frac{\hbar}{\sqrt{2}} \begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 0 \end{pmatrix}, \quad J_y = \frac{\hbar}{\sqrt{2}} \begin{pmatrix} 0 & -i & 0 \\ i & 0 & -i \\ 0 & i & 0 \end{pmatrix}, \quad J_z = \hbar \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & -1 \end{pmatrix}.$$

The time evolution in the adiabatic frame is analogous to the rotation of a spin-1 system around an effective magnetic field given by $\hat{\Omega}_0 = \Omega J_z + \dot{\theta} J_y$, with the amplitude $\Omega_0 = \sqrt{\Omega^2 + \dot{\theta}^2}$. The state vector is initially oriented along the $J_z$, corresponding to the dark state $|\phi_D(t)\rangle$. In the adiabatic limit $\dot{\theta} \ll 1$, the precession axis is $\hat{\Omega}_0 = J_z$, so the state vector does not precess, and the state evolution will be along the dark state without any transition. Regarding the multi-mode driving, we substitute the pulses (26) and have

$$\Omega = \frac{\pi}{2T} \sqrt{16k^2 - 1}, \quad \dot{\theta} = \frac{\pi}{2T}. \qquad (A.3)$$

In this case, the state vector starts to precess when $\dot{\theta} > 0$ and $\hat{\Omega}_0 \neq J_z$, resulting in the nonadiabatic evolution. The resonant condition, $\Omega_0 T = 2k\pi$, which is fulfilled, allows the state vector to be driven back to the target dark state, that is, the ideal target state $|\phi_D(T)\rangle$. Therefore, the linear multi-mode driving can be interpreted by a resonant driving in the adiabatic frame, thus producing a fast and robust state transfer.

## Appendix B. Optimal Control Theory

In order to formulate the problem, one uses the fact that this resonant three-state problem is of SU(2) symmetry. One can separate the dynamics into two independent sub-dynamics each of them being associated to the real and imaginary parts of the initial condition, here taken as real $\langle 1|\psi(0)\rangle = 1$, leading to the system of equations

$$\dot{\mathbf{x}} = (u_p T_1 + u_s T_2)\mathbf{x}, \qquad (B.1)$$

with

$$T_1 = \begin{bmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}, \quad T_2 = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & -1 \\ 0 & 1 & 0 \end{bmatrix}, \quad T_3 = \begin{bmatrix} 0 & 0 & 1 \\ 0 & 0 & 0 \\ -1 & 0 & 0 \end{bmatrix}, \quad (B.2)$$

$[T_1, T_2] = T_3$, and $\mathbf{x} = [x_1 \equiv \langle 1|\psi \rangle, x_2 \equiv i\langle 2|\psi \rangle, x_3 \equiv -\langle 3|\psi \rangle]^T$, where the superscript T denotes the transpose.
To determine optimal control fields $\mathbf{u}(t)$ of a dynamical system $\dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}(t); \mathbf{u}(t))$ (of dimension $N$) with respect to the minimization of a given cost

$$J(\mathbf{u}(t)) = \int_{0}^{T} g(\mathbf{x}(t), \mathbf{u}(t))dt, \qquad (B.3)$$

we use Pontryagin's maximum principle, which provides necessary conditions for optimality [36]. The maximum principle states that the trajectories of the extremal vector $\mathbf{x}(t)$ and of the corresponding adjoint state $\mathbf{p}(t)$ formed by the Lagrange multipliers, $\mathbf{p}(t) \equiv [p_1(t), \dots, p_N(t)]$, fulfill Hamilton's equations

$$\dot{\mathbf{x}} = \frac{\partial H_c}{\partial \mathbf{p}}, \quad \dot{\mathbf{p}} = -\frac{\partial H_c}{\partial \mathbf{x}}, \qquad (B.4)$$

associated to the control pseudo-Hamiltonian

$$H_c(\mathbf{p}(t), \mathbf{x}(t); \mathbf{u}(t)) = \mathbf{p} \cdot \mathbf{f}(\mathbf{x}(t); \mathbf{u}(t)) - p_0 g(\mathbf{x}(t), \mathbf{u}(t)), \qquad (B.5)$$

where the constant $p_0 > 0$ can be chosen for convenience since it amounts to multiply the cost function by a constant. For almost all $0 \le t \le T$ the function $H_c(\mathbf{p}(t), \mathbf{x}(t); \mathbf{u}(t))$ is maximum at certain controls $\mathbf{v}(t) = \mathbf{u}(t)$.

We use the angle coordinates, which defines the dynamical mode $|\phi_0(t)\rangle$ (16), i.e. with $\theta(0) = 0, \theta(T) = \pi/2, \varphi(0) = \varphi(T) = 0$. The equations of the dynamics (14)

$$\dot{\varphi} = u_p \cos \theta - u_s \sin \theta, \quad \dot{\theta} = (u_p \sin \theta + u_s \cos \theta) \tan \varphi \qquad (B.6)$$

can be simplified as

$$\dot{\varphi} = v_p, \quad \dot{\theta} = -v_s \tan \varphi, \qquad (B.7)$$

using a rotation on the control fields

$$\begin{bmatrix} v_p \\ v_s \end{bmatrix} = \begin{bmatrix} \cos \theta & -\sin \theta \\ -\sin \theta & -\cos \theta \end{bmatrix} \begin{bmatrix} u_p \\ u_s \end{bmatrix}. \qquad (B.8)$$

The minimization of the energy (29), $\mathcal{E} = \int_0^T [u_p^2(t) + u_s^2(t)] dt = \int_0^T [v_p^2(t) + v_s^2(t)] dt$, leads to the control pseudo-Hamiltonian

$$H_c = \lambda_{\varphi} v_p - \lambda_{\theta} v_s \tan \varphi - \frac{1}{2}(v_p^2 + v_s^2) \qquad (B.9)$$

with $\lambda = [\lambda_{\varphi}, \lambda_{\theta}]^T$ the costate. The Hamilton equations lead to the equation of motion (B.7) and to

$$\dot{\lambda}_{\varphi} = \frac{v}{\cos^2 \varphi} \lambda_{\theta}, \quad \dot{\lambda}_{\theta} = 0. \qquad (B.10)$$

This implies that $\lambda_{\theta}$ is a constant of motion (taken positive). We apply the maximum principle:

$$\frac{\partial H_c}{\partial v_p} = \lambda_{\varphi} - v_p = 0, \qquad (B.11)$$

$$\frac{\partial H_c}{\partial v_s} = -\lambda_{\theta} \tan \varphi - v_s = 0, \qquad (B.12)$$
which leads to $\lambda_{\varphi} = v_p$ and $\lambda_{\theta} \tan \varphi = -v_s$, i.e.

$$H_c = \frac{1}{2}(v_p^2 + v_s^2) = \frac{1}{2}(\lambda_{\varphi}^2 + \lambda_{\theta}^2 \tan^2 \varphi). \qquad (B.13)$$

We determine the optimal trajectory by quadrature, using $H_c = C/2 = \text{const.}$, since $H_c$ features an effective autonomous system:

$$\dot{\varphi} = \pm \lambda_{\theta} \tan \varphi_m \sqrt{1 - \frac{\tan^2 \varphi}{\tan^2 \varphi_m}}, \qquad (B.14)$$

$$\frac{d\varphi}{d\theta} = \frac{\dot{\varphi}}{\dot{\theta}} = \pm \frac{\tan \varphi_m}{\tan^2 \varphi} \sqrt{1 - \frac{\tan^2 \varphi}{\tan^2 \varphi_m}}, \qquad (B.15)$$

$$\tan \varphi_m = \frac{\sqrt{C}}{\lambda_{\theta}}, \qquad (B.16)$$

where we have assumed $\varphi \ge 0$ (which is satisfied for $u_p$ and $u_s$ both positive) and $\dot{\varphi}(t = T/2) = 0$ by symmetry, leading to $\varphi(T/2) = \varphi_m$ maximum at $t = T/2$ and positive (negative) branch of (B.14),(B.15) for $t \in [0, T/2]$, $\varphi$ increasing from 0 to $\varphi_m$ ($t \in [T/2, T]$, $\varphi$ decreasing from $\varphi_m$ to 0). Taking into account the boundary conditions, we determine the solution of (B.14):

$$\sin \varphi = \sin \varphi_m \sin (\pi t/T) \qquad (B.17)$$

with $\lambda_{\theta}T = \pi \cos \varphi_m$, and the one of (B.15), $\theta_+(\varphi_+)$, for the positive branch of the geodesic (where $\varphi$ increases from 0 to $\varphi_m$, denoted $\varphi_+$):

$$\theta_+ = \operatorname{atan}\left(\frac{\sin \varphi_+}{\sqrt{\tan^2 \varphi_m - \frac{\sin^2 \varphi_+}{\cos^2 \varphi_m}}}\right) - \cos \varphi_m \operatorname{asin}\left(\frac{\sin \varphi_+}{\sin \varphi_m}\right). \qquad (B.18)$$

By imposing by symmetry that $\theta(\varphi_m) = \pi/4$, we get

$$\varphi_m = \pi/3, \qquad (B.19)$$

which leads to

$$\theta_+ = \operatorname{atan}\left(\frac{\sin \varphi_+}{\sqrt{3 - 4 \sin^2 \varphi_+}}\right) - \frac{1}{2}\operatorname{asin}\left(\frac{2}{\sqrt{3}} \sin \varphi_+\right). \qquad (B.20)$$

This also leads to the time dependence of $\varphi$:

$$\sin \varphi = \frac{\sqrt{3}}{2} \sin (\pi t/T) \qquad (B.21)$$

and consequently to the one of $\theta$:

$$\theta = \operatorname{atan}\left(\frac{1}{2} \tan (\pi t/T)\right) - \frac{\pi t}{2T}. \qquad (B.22)$$

For the negative branch (i.e. $\varphi$ decreasing from $\varphi_m$ to 0), Eq. (B.21) still holds. The negative branch $\theta_-(\varphi_-)$ of the geodesic reads:

$$\theta_- = \frac{\pi}{2} - \operatorname{atan}\left(\frac{\sin \varphi_-}{\sqrt{3 - 4 \sin^2 \varphi_-}}\right) + \frac{1}{2} \operatorname{asin}\left(\frac{2}{\sqrt{3}} \sin \varphi_-\right). \qquad (B.23)$$
We note that $\arcsin(\sin \varphi_{-} / \sin \varphi_{m}) = \pi - \pi t / T$ in this case since $\varphi_{m}$ is maximum. Equation (B.22) hold as well if we assume the definition

$$ \operatorname{atan}(x) \in \begin{cases} [0, \pi/2[ & \text{for } x \ge 0, \\ ]\pi/2, \pi[ & \text{otherwise.} \end{cases} \qquad (B.24) $$

Inserting these relations into the definition of the pulses (14), we obtain for the optimal pulses with respect to the energy as follows,

$$ u_p(t) = \frac{\sqrt{3\pi}}{2T} \cos \left( \frac{\pi t}{2T} \right), \quad u_s(t) = \frac{\sqrt{3\pi}}{2T} \sin \left( \frac{\pi t}{2T} \right). \qquad (B.25) $$

## References

[1] Bergmann K, Theuer H and Shore B 1998 *Rev. Mod. Phys.* **70** 1003

[2] Guérin S and Jauslin H R 2003 *Adv. Chem. Phys.* **125** 147

[3] Král P, Thanopulos I and Shapiro M 2007 *Rev. Mod. Phys.* **79** 53

[4] Saffman M, Walker T G and Mølmer K 2010 *Rev. Mod. Phys.* **82** 2313

[5] Vitanov N V, Rangelov A A, Shore B W and Bergmann K 2017 *Rev. Mod. Phys.* **89** 015006

[6] Schiemann S, Kuhn A, Steuerwald S and Bergmann K 1993 *Phys. Rev. Lett.* **71** 3637

[7] Hioe F T and Eberly J H 1981 *Phys. Rev. Lett.* **47** 838

[8] Oreg J, Hioe F T and Eberly J H 1984 *Phys. Rev. A* **29** 690

[9] Caroll C E and Hioe F T 1988 *J. Opt. Soc. Am. B* **5** 1335

[10] Oreg J, Bergmann K, Shore B W and Rosenwaks S 1992 *Phys. Rev. A* **45**, 4888

[11] Vitanov N V and Stenholm S 1997 *Phys. Rev. A* **55** 648

[12] Levitt M H and Freeman R 1979 *J. Magn. Reson.* **33** 473

[13] Freeman R, Kempsell S P and Levitt M H 1980 *J. Magn. Reson.* **38** 453

[14] Levitt M H 1986 *Prog. NMR Spectrosc.* **18** 61

[15] Freeman R 1977 Spin Choreography (Spektrum, Oxford).

[16] Torosov B T, Guérin S and Vitanov N V 2011 *Phys. Rev. Lett.* **106** 233001

[17] Torosov B T and Vitanov N V 2013 *Phys. Rev. A* **87** 043418

[18] Torrontegui E, Ibàñez S, Martínez-Garaot S, Modugno M, del Campo A, Guéry-Odelin D, Ruschhaupt A, Chen X and Muga J G 2013 *Adv. At. Mol. Opt. Phys.* **62** 117

[19] Guéry-Odelin D, Ruschhaupt A, Kiely A, Torrontegui E, Martínez-Garaot S and Muga J G, arXiv:1904.08448.

[20] Torrontegui E, Chen X, Modugno M, Schmidt S, Ruschhaupt A and Muga J G 2012 *New J. Phys.* **14** 013031

[21] Daems D, Ruschhaupt A, Sugny D and Guérin S 2013 *Phys. Rev. Lett.* **111** 050404

[22] Van-Damme L, Schraft D, Genov G T, Sugny D, Halfmann T and Guérin S 2017 *Phys. Rev. A* **96** 022309

[23] Hennrich M, Legero T, Kuhn A and Rempe G 2000 *Phys. Rev. Lett.* **85** 4872

[24] Kuhn A, Hennrich M and Rempe G 2002 *Phys. Rev. Lett.* **89** 067901

[25] Wilk T, Webster S C, Specht H P, Rempe G and Kuhn A 2007 *Phys. Rev. Lett.* **98** 063601

[26] Sørensen J L, Møller D, Iversen T, Thomsen J B, Jensen F, Staanum P, Voigt D and Drewsen M 2006 *New J. Phys.* **8** 261

[27] Dridi G, Guérin S, Hakobyan V, Jauslin H R and Eleuch H 2009 *Phys. Rev. A* **80** 043408

[28] Chen X, Lizuain I, Ruschhaupt A, Guéry-Odelin D and Muga J G 2010 *Phys. Rev. Lett.* **105** 123003

[29] Chen X and Muga J G 2012 *Phys. Rev. A* **86** 033405

[30] Torosov B T and Vitanov N V 2013 *Phys. Rev. A* **87** 043418

[31] Li Y-C and X Chen 2016 *Phys. Rev A* **94** 063411
[32] Mortensen H L, Sørensen J J W H, Mølmer K and Sherson J F 2018 *New J. Phys.* **20** 025009

[33] Bruns A, Genov G T, Hain M, Vitanov N V and Halfmann T 2018 *Phys. Rev. A* **98** 053413

[34] Lewis H R and Riesenfeld W B 1969 *J. Math. Phys.* **10** 1458

[35] Boscain U, Charlot G, Gauthier J-P, Guérin S and Jauslin H-R 2002 *J. Math. Phys.* **43** 2107

[36] D’Alessandro D 2007 *Introduction to quantum control and dynamics*, Taylor & Francis Ltd

[37] Du Y-X, Liang Z-T, Li Y-C, Yue X-X, Lv Q-X, Huang W, Chen X, Yan H and Zhu S-L 2016 *Nat. Commun.* **7** 12479

[38] Zhou B B, Baksic A, Ribeiro H, Yale C G, Heremans F J, Jerger P C, Auer A, Burkard G, Clerk A A and Awschalom D D 2017 *Nat. Phys.* **13** 330

[39] Vepsäläinen A, Danilin S and Paraoanu G S 2019 *Sci. Adv.* **5** eaau5999