# Quantum Control by Few-Cycles Pulses: The Two-Level Problem

François Peyraut, Frédéric Holweck, Stéphane Guérin

<figure><img src="image_1.png" /></figure>

► To cite this version:

François Peyraut, Frédéric Holweck, Stéphane Guérin. Quantum Control by Few-Cycles Pulses: The Two-Level Problem. Entropy, 2023, 25 (2), pp.212. 10.3390/e25020212. hal-04047305

HAL Id: hal-04047305

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

Submitted on 2 Jun 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Quantum Control by Few-Cycles Pulses:
The Two-Level Problem

François Peyraut<fnref n="1" />, Frédéric Holweck <fnref n="1" />,<fnref n="2" /> and Stéphane Guérin <fnref n="3" />,*

**Abstract:** We investigate the problem of population transfer in a two-states system driven by an external electromagnetic field featuring a few cycles, until the extreme limit of two or one cycle. Taking the physical constraint of zero-area total field into account, we determine strategies leading to ultrahigh-fidelity population transfer despite the failure of the rotating wave approximation. We specifically implement adiabatic passage based on adiabatic Floquet theory for a number of cycles as low as 2.5 cycles, finding and making the dynamics follow an adiabatic trajectory connecting the initial and targeted states. Nonadiabatic strategies with shaped or chirped pulses, extending the $\pi$-pulse regime to two- or single-cycle pulses, are also derived.

**Keywords:** quantum control; quantum system driven by an external field; adiabatic passage; adiabatic Floquet theory

## 1. Introduction

Quantum (or coherent) control aims at developing methods for manipulating quantum dynamical processes at the atomic or molecular scale by shaped external electromagnetic fields [1,2]. It can be conceptually formulated as the design of the field constituents (phase, amplitude, and polarization) driving the quantum system, such as atom, molecule [3], photonic field [4], nanostructures (quantum dot) and plasmonic field [5], Bose-Einstein condensate [6,7], superconducting circuit [8,9], ..., from an initial state to a target state featuring the desired outcome. The latter can be a single bare state of the system, a coherent superposition of bare states, or the full propagator (i.e., for any initial state) as it is the case in quantum computation for which quantum gates are devised [10–13]. For instance, applications to quantum information target the superposition of a few states (qubit and its generalization) while the control of chemistry deals with wave packets made up of many states, such as, e.g., an aligned or oriented molecule corresponding to a superposition of many rotational states [14,15]. Various techniques have been proposed and implemented: Adiabatic passage [16,17] and adiabatic Floquet theory [18], composite pulses [19–21], optimal control [22–24], robust single-shot [25] and optimal [26] shaped pulses as a variant of shortcut to adiabaticity [27,28], and reinforcement learning [29,30].

Ultrashort laser pulses that last only a few optical cycles have been transformative tools for studying and manipulating light-matter interactions [31]. They allow in particular high nonlinear effects for strong-field processes in the subfemtosecond extreme ultraviolet/X-ray regime [32]. A different single-cycle THz regime has been used to experimentally demonstrate enhanced orientation of molecules [33]. Such single-cycle THz pulses have been investigated for the control of the rotational wavepackets, resulting in long-lasting orientation [34] and orientational quantum revivals [35].

The goal in this paper is twofold: We implement adiabatic Floquet theory [18] and explore its limit for few cycles. We show that it remarkably applies for a number of cycles as low as 2.5 cycles. In particular, one succeeds to realize an *ultrahigh-fidelity* population

<figure><img src="image_3.png" /></figure>

### check for updates

**Citation:** Peyraut, F.; Holweck, F.; Guérin, S. Quantum Control by Few-Cycles Pulses: The Two-Level Problem. *Entropy* **2023**, *25*, 212.
https://doi.org/10.3390/e25020212

Academic Editor: Ronnie Kosloff

Received: 2 December 2022

Revised: 3 January 2023

Accepted: 10 January 2023

Published: 22 January 2023

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2023 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">1 ICB, UMR 6303, CNRS, University Bourgogne Franche-Comté, UTBM, 90010 Belfort, France</footnote>

<footnote n="2">2 Department of Mathematics and Statistics, Auburn University, Auburn, AL 36849, USA</footnote>

<footnote n="3">3 ICB, UMR 6303, CNRS, University of Bourgogne, 21000 Dijon, France</footnote>

<footnote>* Correspondence: sguerin@u-bourgogne.fr</footnote>
transfer (defined by a deviation from the perfect transfer below or of the order 10⁻⁴) by an appropriate design based on parallel adiabatic passage. Its obstruction for a lower number of cycles is due to the presence of nonadiabatic terms emerging from the counter-rotating term beyond the rotating wave approximation (RWA) [36]. We also investigate how one can control complete population transfer when a two-level system is driven by a shaped and chirped few-cycles pulse, typically N < 3, taking the constraint of a zero field area into account. One generalizes the two standard techniques in this context, adiabatic passage and π-pulse method [36]. The latter is of importance as it is known to achieve the so-called quantum speed limit in RWA [24]. The former is known to induce robust transitions [17]. We achieve complete transfer for a single cycle N = 1 pulsed laser featuring an appropriate chirped (time-dependent) frequency around the resonance.

After having defined the model in Section 2, we apply adiabatic Floquet theory in Section 3. In Section 4, we introduce the notion of generalized π-pulse in the context of few cycles. A generalized π-pulse consists in expanding the Rabi frequency as a Fourier series and selecting the appropriate modes leading to resonance and allowing complete transfer. This gives rise to specific strategies depending on the modes chosen in the Fourier expansion and on the chirping of the frequency. Results from numerical optimization are presented and connected to the preceding strategies in Section 6. We finally conclude in Section 7.

## 2. The Model with a Few-Cycle Pulse

We consider the two-level system {|−〉 , |+〉} driven by a few-cycle pulse (in units such that ħ = 1)

$$H = \begin{bmatrix} -\omega_0/2 & \Omega(t) \cos \phi(t) \\ \Omega(t) \cos \phi(t) & \omega_0/2 \end{bmatrix}, \qquad (1)$$

where Ω(t) is the Rabi frequency defined as Ω(t) = −μ₀ε(t)/ħ (considered positive without loss of generality) with the dipole moment μ₀₁ coupling the two states; ω₀ > 0, and φ(t) are the Bohr frequency and the instantaneous phase, respectively. The effective (instantaneous) frequency ω(t) := ˙φ(t) can be modified at will through the phase φ(t).

One can exhibit an instantaneous relative phase φ₀(t) with respect to the Bohr frequency:

$$\phi(t) = \omega_0(t - t_i) - \phi_0(t). \qquad (2)$$

and the corresponding detuning Δ:

$$\Delta(t) := \omega_0 - \omega(t) = \dot{\phi}_0(t), \quad \phi_0(t) = \int_{t_i}^{t} \Delta(s)ds + \phi_{0,i} \qquad (3)$$

with tᵢ = −T/2 the initial time, tᵣ = −tᵢ the final time with full duration T = tᵣ − tᵢ, and the initial phase φ₀ᵢ. We will restrict our study to the near-resonant case ˙φ ≪ ω₀ for the ease of implementation. The number N of oscillations in the pulse (with respect to the resonant frequency ω₀ in this near-resonant case) is such that

$$\mathcal{N} = T \frac{\omega_0}{2\pi}. \qquad (4)$$

The physical implementation of the few-cycle field imposes its zero time-integrated area [37]:

$$\int_{t_i}^{t_f} \Omega(t) \cos(\phi(t)) dt = 0, \qquad (5)$$
which will be imposed for all the derived fields in this work. After application of the resonant (or rotating wave) transformation:

$$R = \begin{bmatrix} 1 & 0 \\ 0 & e^{-i\phi(t)} \end{bmatrix}, \qquad (6)$$

we obtain (where we have for convenience shifted the origin of energy such that a traceless matrix is derived)

$$\tilde{H} = R^{\dagger}HR - iR^{\dagger}\frac{\partial R}{\partial t} = \frac{1}{2} \begin{bmatrix} -\Delta(t) & \Omega(t) \\ \Omega(t) & \Delta(t) \end{bmatrix} + \frac{\Omega(t)}{2} \begin{bmatrix} 0 & e^{-2i\phi(t)} \\ e^{2i\phi(t)} & 0 \end{bmatrix}. \qquad (7)$$

We notice that the rotating wave transformation changes not the amplitude of the solution but its phase. The first matrix corresponds to the standard (resonant) RWA. The second matrix defines the counter-rotating term. The full Hamiltonian (7) does not feature any approximation in the two-state model. It is well known that the counter-rotating Hamiltonian can be neglected in the limit $\Omega \ll \omega_0$, since, in this case, this term features fast oscillations that have a weak effect compared to the resonant, non-oscillatory, RWA Hamiltonian.

In the frame of the Floquet theory (see Appendix A and the next section), the RWA Hamiltonian characterizes the interaction between the states $|-, 0\rangle$ and $|+, -1\rangle$, corresponding to the bare states of the system dressed with 0 photon and $-1$ photon, respectively, which are near degenerate (exactly degenerate on exact resonance). This means that the system can emit a photon when its state is transferred from state $|+\rangle$ to state $|-\rangle$ (or reciprocally can absorb a photon from state $|-\rangle$ to state $|+\rangle$).

Taking the simple example of a resonant ($\Delta = 0$) and flat (constant) pulse of amplitude $\Omega_0$ and duration $T$, the well-known $\pi$-pulse condition $T\Omega_0 = \pi$ in (4) gives $\Omega_0/\omega_0 = 1/(2N)$. Thus, the constraint of a few oscillations in the field, typically $N < 3$, implies that $\Omega_0$, while still smaller than $\omega_0$, becomes of the same order. This condition clearly prevents to apply the resonant approximation. This can be reformulated as follows: a small number of cycles in the pulse corresponds to a broadening of the spectrum of the field instantaneously available and thus a non-negligible influence of the counter-rotating term.

### 3. Few-Cycle-Pulse Adiabatic Floquet Theory

We consider the model (1) written with the time-dependent phase $\theta + \phi(t)$, where the initial phase $\theta \equiv -\phi_{0,i}$ serves for the Floquet representation, and the explicit dependence of the Rabi frequency $\Omega(t)$:

$$H^{\Omega(t)}(\theta + \phi(t), t) = \frac{\hbar\omega_0}{2} \begin{bmatrix} -1 & 0 \\ 0 & 1 \end{bmatrix} + \hbar\Omega(t) \cos(\theta + \phi(t)) \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}. \qquad (8)$$

The quasienergy operator (A19)

$$K \equiv K^{\Omega,\omega} = -i\hbar\omega \frac{\partial}{\partial\theta} + H^{\Omega}(\theta) \qquad (9)$$

features the two parameters $\Omega \equiv \Omega(t)$ and $\omega \equiv \omega(t)$ that can be designed independently as functions of time. They will both normalized with $\omega_0$. The resonant transformation (6) is represented in the Floquet picture by the transformation (which dresses the upper state with minus one photon):

$$R = \begin{bmatrix} 1 & 0 \\ 0 & e^{-i\theta} \end{bmatrix}, \qquad (10)$$

which leads to

$$R^{\dagger}KR = -i\hbar\omega \frac{\partial}{\partial\theta} + \frac{\hbar}{2} \begin{bmatrix} -\Delta(t) & \Omega(t) \\ \Omega(t) & \Delta(t) \end{bmatrix} + \frac{\hbar\Omega(t)}{2} \begin{bmatrix} 0 & e^{-2i\theta} \\ e^{2i\theta} & 0 \end{bmatrix}. \qquad (11)$$
Note that we have omitted as before a diagonal term equal to $-\omega/2$ times identity in order to work with traceless matrices (without loss of generality for population transfer consideration).

Adiabatic passage through the (one-photon) resonance, using a pulsed (symmetric) Rabi frequency $\Omega(t) = \Omega_0\Lambda(t)$ with $\Omega_0$ the peak value and the shape $0 \le \Lambda(t) \le 1$, with $\Lambda(\pm\infty) \to 0$, $\Lambda(0) = 1$, and of characteristic width $\tau$, and a chirped frequency $\omega(t)$, of characteristic width around the resonance $\Delta_0$, requires in general $\Omega_0 \sim \Delta_0 \gg 1/\tau$. The chirped frequency induces for the full coupling $\Omega(t) \cos\phi(t)$ in Equation (1) the phase $\phi(t) = \int_{t_i}^t \omega(s)ds$, which can be rewritten as a function of the detuning; see Equation (2): $\phi(t) = \omega_0(t-t_i) - \int_{t_i}^t \Delta(s)ds - \phi_{0,i}$. We consider smooth Gaussian pulses $\Lambda(t) = e^{-(t/\tau)^2}$, which allows in principle a transfer exponentially accurate as a function of the pulse duration [38] (see Appendix A.2). Defining $\tau\Omega_0 \sim T\Omega_0 = 2\kappa\pi$ (with a larger $\kappa$ ensuring a better adiabatic passage) leads to $\Omega_0/\omega_0 = \kappa/\mathcal{N}$. This shows that adiabatic passage, for a given $\kappa$, breaks more strongly the rotating wave approximation for a smaller number $\mathcal{N}$ of cycles, since $\Omega_0$ gets closer to $\omega$. This requires to take into account the full quasienergy operator.

The quasienergies $\lambda_{\pm,k}$ can be obtained from the numerical diagonalization of the quasienergy operator (9) [or equivalently (11)], apart a change of energy reference], as shown in Figure 1. We emphasize that the quasienergies do not depend on the number of cycles considered $\mathcal{N}$. The complete spectral information is contained in one Floquet zone composed of two surfaces $\lambda_{\pm,0}$ within a given band of energy of width $\hbar\omega$. The other surfaces can be constructed using the periodicity of the spectrum: $\lambda_{\pm,k} \equiv \lambda_{\pm,0} + k\hbar\omega$, for any (positive or negative) integer $k$. In Figure 1 we display a few surfaces. Figure 2 shows cross-sections, one for a constant $\omega$ (left frame) and for constant $\Omega$ (right frame). Around the one-photon resonance $\omega \approx \omega_0$ (for $\Omega \ll \omega_0$), one can recognize the surfaces of Figure 1 of Ref. [17]: For $\Omega = 0$, the horizontal line of energy $-0.5\hbar\omega_0$ corresponds to the lower state of the system $|-,0\rangle$, with the notation $|\pm,k\rangle$ for the bare state $|\pm\rangle$ dressed with $k$ photon, i.e., the state $|\pm\rangle \otimes |k\rangle$ with $|k\rangle \equiv e^{ik\theta}$. The crossing line (at $\omega = \omega_0$) corresponds to the upper state dressed with minus one photon $|+, -1\rangle$. One can distinguish in the plane $\Omega = 0$ resonances appearing as *crossings* at $\omega = \omega_0/(2k+1)$ with $k = 0, 1, 2, \dots$. For $\Omega \ne 0$, they become *avoided crossings*. For $\omega = \omega_0/(2k)$, we have exact crossings for any $\Omega$, due to the particular symmetry of this model (one can see an example for $k=1$). This means that only odd numbers of photons can be absorbed (or emitted) in such a system. The maxima of the upper surface correspond to crossings (for any $\Omega$), and the valleys to avoided crossings (i.e., to resonances). One can observe that for increasing $\Omega$, the position of the resonances are shifted in the direction of larger $\omega$. This can be interpreted as a Stark shift of the states. This implies that moving along a straight line with $\omega \approx \omega_0$ for growing $\Omega$ allows one to cross dynamically the three-photon resonance, next the five-photon resonance, and so on, as shown in Figure 2. The three-photon resonance avoided crossing represents thus the first dynamical obstruction that can prevent the control of population transfer by a standard chirped pulse for a large enough Rabi frequency.

It is important to emphasize that the labeling of the surfaces can be only local due to the multiple crossings. This is this property which, when the resonance is dynamically crossed (i.e with a change of sign of the detuning during the dynamics), will allow a population transfer by adiabatic passage of a single Floquet eigenstate that will connect state $|-,0\rangle$ at early time with state $|+, -1\rangle$ at late time [17]. We denote around the one-photon resonance ($\omega \approx \omega_0$) the upper (lower) surface as $\lambda_+$ ($\lambda_-$), respectively. Near the one-photon resonance, the surface $\lambda_+$ ($\lambda_-$) is connected to the state $|+, -1\rangle$ ($|-,0\rangle$) for $\omega < \omega_0$, i.e., $\Delta > 0$, and to the state $|-,0\rangle$ ($|+, -1\rangle$) for $\omega > \omega_0$, i.e., $\Delta < 0$, respectively.
<figure><img src="image_5.png" /><figcaption>Figure 1. Quasienergy surfaces (in units of $\hbar\omega_0$) as functions of $\omega/\omega_0$ and $\Omega/\omega_0$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 2. Cross section (full lines) of the quasienergy surfaces (in units of $\hbar\omega_0$) of Figure 1 as functions of (i) $\Omega/\omega_0$ for a constant $\omega = 0.8\omega_0$, i.e., $\Delta = 0.2\omega_0$, where the dashed lines are the quasienergies in the resonant approximation (left frame); (ii) $\Delta/\omega_0$ for a constant $\Omega = 0.4\omega_0$ (full line) and $\Omega = 0$ (dashed lines) (right frame). One can distinguish the three-photon resonance as an avoided crossing located around $\Omega = 1.4\omega_0$ (left frame) and identified by a circle (right frame). The five-photon-resonance avoided crossing can be seen around $\Omega = 2.7\omega_0$ (left frame).</figcaption></figure>

To achieve population transfer, we apply the strategy of parallel adiabatic passage, where the adiabatic dynamics follows a level line $\lambda_+ - \lambda_- = \text{const.}$ around the resonance, which is known to be efficient in the context of RWA for a smooth Rabi frequency [39–41]. In the RWA, the detuning crosses the resonance according to

$$ \Delta(t) = \pm \sqrt{\Omega_0^2 - \Omega^2(t)}, \quad \Delta_0 = \Omega_0 > 0, \qquad (12) $$

with the sign $+$ (–) for $t < 0$ ($t > 0$), i.e., $\Delta(-\infty) = \Delta_0$ [$\Delta(+\infty) = -\Delta_0$] the initial (final) detuning, respectively. We remark that, in the RWA, the symmetric situation: the sign
– (+) for $t < 0$ ($t > 0$), i.e., $\Delta(-\infty) = -\Delta_0$ [$\Delta(+\infty) = \Delta_0$] the initial (final) detuning, respectively, leads to the same population transfer. As mentioned above, the quality of adiabatic passage can be evaluated by the value of $\tau\Delta_0$. It has been shown in [39,41] that Gaussian-pulse parallel adiabatic passage mainly features (i) high fidelity (corresponding to population transfer with a deviation from the perfect transfer of the order $10^{-3}$) for the product $\tau\Delta_0$ as low as $\tau\Delta_0 \approx 2.4$, (ii) a narrow perfect transfer for $\tau\Delta_0 \approx 2.58$, and (iii) ultrahigh fidelity for $\tau\Delta_0 \gtrsim 4.5$.

In the situation of small number of cycles $N$, one has additionally to keep a sufficient distance from the neighbouring three-photon resonance. The level lines around the one-photon resonance correspond thus to the constant distances $|\lambda_+ - \lambda_-|$ which are compared to the three-photon resonance, as represented in Figure 3. The driving phase corresponding to a given trajectory of Figure 3, where the pulse amplitude is chosen as a Gaussian pulse of amplitude $\Omega_0 = \Delta(-\infty)$, is given by (2). The initial phase $\phi_{0,i}$, which does not modify the dynamics in the adiabatic limit, will be chosen to satisfy the zero time-integrated area (5).

<figure><img src="image_7.png" /><figcaption>Figure 3. Contour plot of $|\lambda_+ - \lambda_-|$ around the one-photon resonance. The dashed line shows the position of the three-photon resonance, which occurs for a smaller $\Omega/\omega_0$ when one considers a larger detuning $\Delta/\omega_0$. The level line trajectory (full line) slightly avoids the three-photon-resonance curve (for positive detunings).</figcaption></figure>

Considering $N$ oscillations in the pulse, i.e., $T\omega_0/2\pi = N$, and Gaussian pulses with the estimate $T \sim 4\pi\tau/3$ for the pulse duration gives $\tau\omega_0 \approx 3N/2$, i.e.,

$$ \tau\Delta_0 \approx \frac{3}{2} N \frac{\Delta_0}{\omega_0}. \qquad (13) $$

This shows that smaller values $N$ require larger $\Delta_0/\omega_0$ to keep the same $\tau\Delta_0$.

We have obtained the most extreme situation with a low number of cycles $N \approx 2.5$ allowing ultrahigh-fidelity adiabatic passage with the following parameters. We consider the trajectory shown in Figure 3 around the one-photon resonance $\Delta_0/\omega_0 = 0.625$, which slightly avoids the three-photon-resonance curve, and the smallest quantity $\tau\Delta_0 = 2.35$ satisfying adiabatic passage (see the discussion above). We have determined numerically the value $\phi_{0,i} \approx -0.3032\pi$ satisfying the zero time-integrated area (5). Despite the close presence of the three-photon resonance, we have obtained a remarkable ultrahigh fidelity, as shown in Figure 4. Note that we have chosen a dynamics such that $\Delta(-\infty) = \Delta_0$; we have checked that we have the same final result if we start with $\Delta(-\infty) = -\Delta_0$ (but with a different initial phase $\phi_{0,i}$). We notice that, in practice, since the dynamics does not
strictly satisfy the adiabatic conditions, the population dynamics slightly depends on this initial phase.

Decreasing the number of cycles degrades the fidelity because one cannot satisfy a sufficiently large $\tau\Delta_0$ from (13).

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Dynamics along the level line (a) of Figure 3 with $\Delta(-\infty) = 0.625\omega_0$ and $\tau\omega_0 = 3.75$ for a Gaussian pulse $\Omega(t) = \Omega_0 e^{-(t/\tau)^2}$, giving $N \approx 2.5$. Upper frame: Populations $P_j := |\langle j|\varphi(t)\rangle|^2$. We obtain an ultrahigh fidelity $P_+ \approx 0.9999$. Middle frame: the normalized Rabi frequency (full line) with its envelope (dotted line) and the instantaneous normalized detuning (dashed line). Lower frame: the instantaneous quasienergies (in units of $\hbar\omega_0$). The dynamics follows $\lambda_-$ [connected to $|-, 0\rangle$ ($|+, -1\rangle$) at early (late) times], close to the three-photon resonance quasienergy $\lambda_{+,3}$ at early times.</figcaption></figure>

# 4. Few-Cycle Generalized $\pi$-Pulse and Non-Adiabatic Regimes

## 4.1. Definition

In this section and the next ones, one studies few-cycle pulses driving the dynamics in nonadiabatic regimes and analyzes how one can control the population transfer, in particular with ultrahigh fidelity.

We introduce the notion of few-cycle generalized $\pi$-pulse by expanding the Rabi frequency as a Fourier series:

$$ \Omega(t) = \sum_{n \ge 1} \Omega_n(t), \quad \text{with } \Omega_n(t) = \Omega_n \Lambda_n(t) \text{ and } \Omega_n \ge 0, \qquad (14) $$

with the envelope of each pulsed mode defined as

$$ \Lambda_n(t) = \frac{1}{2} \left[ 1 - (-1)^n \cos\left(\frac{n\omega_0}{N}t\right) \right], \qquad (15) $$

which satisfies $0 \le \Lambda_n(t) \le 1$. This particular Fourier expansion (14) implies that $\Omega(t) \ge 0$ and $\Omega(t_i) = \Omega(t_f) = 0$. Moreover, one can prove that, for $\varphi(t) = \omega_0(t - t_i) - \phi_{0,i}$ (i.e., exact resonance, $\Delta = 0$, and for any initial phase), imposing $\Omega_N = 0$ guarantees that the zero area condition (5) is satisfied.
If one considers a single field $\Omega(t) \equiv \Omega_1(t)$, the high-frequency effective Hamiltonian from the counter-rotating term (11) reads [see Equation (A52)], where we have neglected the terms of Equation (A55) induced by the time-dependence of $\Omega_1(t)$:

$$ \tilde{H} \approx \frac{1}{2} \begin{bmatrix} -\left(\frac{\Omega_1^2(t)}{4\omega_0} + \Delta(t)\right) & \Omega_1(t) \\ \Omega_1(t) & \frac{\Omega_1^2(t)}{4\omega_0} + \Delta(t) \end{bmatrix}. \quad (16) $$

The diagonal term can be neglected with respect to the coupling when (considering the condition for the peak Rabi frequencies) $\Omega_1 \gg \frac{\Omega_1^2}{2\omega_0}$. In the $\pi$-pulse regime, the pulse area of the pulse envelope is then $\Omega_1 T/2 = \pi$, and the above condition becomes

$$ \mathcal{N} \gg 1/2 \quad (17) $$

to neglect the diagonal term. This consists precisely to one of the RWA conditions, i.e., a large number of resonant cycles.

## 4.2. Few-Cycle Resonant Rabi Oscillations

In the resonant situation, $\Delta = 0$, we consider the number of cycles $\mathcal{N} = 2$ and $\Omega_2 = 0$ such that the zero area condition (5) is satisfied. Figure 5 shows the few-cycle Rabi oscillations produced for the phase of the field $\phi(t) = \omega_0(t - t_i) \pm \pi/2$ (i.e., $\phi_{0,i} = \mp \pi/2$) leading to the maximum population transfer (see Figure 6 and the discussion below). They globally significantly deviate from the standard RWA Rabi oscillations; however, they feature a rather good maximum final transfer, $P_+ \approx 0.99$ for a peak amplitude $\Omega_1 \approx 2\pi/T = \omega_0/\mathcal{N}$ close to the standard $\pi$-pulse condition (with a relative deviation of order $2 \times 10^{-3}$). It also shows an ultrahigh-fidelity population transfer but for a high and quite unrealistic pulse area $T\Omega_1/2 \approx 6.38\pi$. Figure 6 shows the weak dependence of the population transfer over the initial phase $\phi_{0,i}$ near the $\pi$-pulse regime. The maximum is obtained for the initial phase $\phi_{0,i} = \pm \pi/2$. We have determined that the population transfer depends on $\phi_{0,i}$ much more strongly for higher pulse area regimes (not shown).

<figure><img src="image_9.png" /><figcaption>Figure 5. Two-cycle resonant Rabi oscillations: Final population transfer probability to the upper state as a function of the pulse area ($\times \pi$) for a single pulse $\Omega_1(t)$ and the number of cycles $\mathcal{N} = 2$ using the full model (1). The transfer reaches a local maximum ($P_+ \approx 0.99$) close to the standard $\pi$-pulse condition but does not reach ultrahigh fidelity (except for an undesirable large pulse area $T\Omega_1/2 \approx 6.38\pi$).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Contour plot of $\log_{10}(1 - P_+)$ with $P_+$ the final population transfer probability around the $\pi$-pulse regime in the condition of Figure 5, but with varying initial phases $\phi_{0,i}$.</figcaption></figure>

Figure 7 shows the single-cycle Rabi oscillations, $N=1$, for the phase of the field $\phi(t) = \omega_0(t - t_i) \pm \pi/2$, which satisfies the zero time-integrated area condition (5). Surprisingly, the maximum transfer (which occurs close to the standard $\pi$-pulse condition) is even better in this case compared to $N=2$ (for the same pulse area), while the system deviates more from the RWA according to condition (17). This can be attributed to the effect of significant population transfer by off-resonant zero-area pulse described in [42], which counterbalances the deviation from RWA. We notice that this effect does not apply efficiently when the field amplitude increase due to the form of the interaction, which is not fully compatible with the hypothesis of [42] in the limit of strong field (as it does not tend to a Dirac $\delta$ distribution type interaction).

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Single-cycle resonant Rabi oscillations: Same as Figure 5 but for a single cycle $N=1$. The transfer reaches a good fidelity $P_+ \approx 0.996$ close to the standard $\pi$-pulse condition: $\Omega_1 T/2 \approx 1.018\pi$.</figcaption></figure>

However, in both situations, $N=1,2$, the transfer is not perfect, and cannot reach ultrahigh fidelity near the $\pi$-pulse regime.
Our goal in the following consists in determining strategies of pulse shaping inducing ultrahigh fidelity in a π-pulse regime. We will show that this will be achieved by adding higher-order terms Ω<sub>n</sub>, that will allow us to alleviate condition (17), or by considering a chirped frequency (next Section). A more systematic optimization procedure will be conducted in Section 6.

4.3. 2 N -Resonance Strategy

We consider here the situation $\mathcal{N} \ge 2$ and $\phi(t) = \omega_0 t$, i.e., an exactly resonant problem $\Delta = 0$, with $\phi_{0,i} = -\omega_0 t_i$. We add a single $2\mathcal{N}$-mode to the main mode: $\Omega(t) = \Omega_1(t) + \Omega_{2\mathcal{N}}(t)$. The zero area condition (5) is automatically satisfied since $\Omega_{\mathcal{N}} = 0$. In this case, the fast oscillating $n = 2\mathcal{N}$-mode term features a resonance with the counter-rotating coupling:

$$
\begin{equation}
\begin{split}
\tilde{H} ={}& \frac{1}{2} \left( \Omega_1(t) + \frac{1}{4} \Omega_{2\mathcal{N}} \right) \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} \\
& + \frac{1}{2} \left( \Omega_1(t) \begin{bmatrix} 0 & e^{-2i\omega_0 t} \\ e^{2i\omega_0 t} & 0 \end{bmatrix} + \frac{1}{4} \Omega_{2\mathcal{N}} \begin{bmatrix} 0 & e^{-2i\omega_0 t} \\ e^{2i\omega_0 t} & 0 \end{bmatrix} - \frac{1}{4} \Omega_{2\mathcal{N}} \begin{bmatrix} 0 & e^{2i\omega_0 t} \\ e^{-2i\omega_0 t} & 0 \end{bmatrix} \right) \\
& - \frac{1}{8} \Omega_{2\mathcal{N}} \begin{bmatrix} 0 & e^{-4i\omega_0 t} \\ e^{4i\omega_0 t} & 0 \end{bmatrix}.
\end{split}
\tag{18}
\end{equation}
$$

If we neglect the oscillating terms [second and third lines of (18)], the area of the resulting
coupling is

$$
\int_{-T/2}^{T/2} \left( \Omega_1(t) + \frac{1}{4} \Omega_{2\mathcal{N}} \right) dt = \frac{T}{2} \left( \Omega_1 + \frac{1}{2} \Omega_{2\mathcal{N}} \right), \quad (19)
$$

which should be $\pi$ at the lowest order:

$$
\Omega_1 + \frac{1}{2} \Omega_{2\mathcal{N}} = \frac{\omega_0}{\mathcal{N}}. \tag{20}
$$

According to (A52) and neglecting the faster oscillating term, one concludes that the
dominant correction is given by the Ω₁(t) counter-rotating term [first term of the second
line of (18)]. Considering the peak terms, one can evaluate the (diagonal) correction:

$$
\tilde{H}_{\max} = \frac{1}{2} \left( \Omega_1 + \frac{1}{4} \Omega_{2\mathcal{N}} \right) \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} + \frac{\Omega_1^2}{8\omega_0} \begin{bmatrix} -1 & 0 \\ 0 & 1 \end{bmatrix}. \quad (21)
$$

Its difference should be much smaller than the coupling term, giving the condition of
the ratio

$$
R = \frac{\Omega_1^2}{2\omega_0 \left(\Omega_1 + \frac{1}{4}\Omega_{2\mathcal{N}}\right)} = \frac{\Omega_1^2}{\omega_0 \left(\Omega_1 + \frac{\omega_0}{\mathcal{N}}\right)} \ll 1. \quad (22)
$$

We implemented a systematic search over the parameter Ω₁ for $\mathcal{N} = 2$ and found an ultrahigh-fidelity transfer $P_+ \approx 0.9999$ for $\Omega_1 \approx 0.276\omega_0$ [leading to $\Omega_4 \approx 0.448\omega_0$ according to (20)], corresponding to an equal distribution of both contributions in (20): $\Omega_1 \sim \frac{1}{2}\Omega_{2\mathcal{N}} \sim \omega_0/(2\mathcal{N})$. This gives in this situation for the above condition

$$
\mathcal{N} \gg 1/6, \tag{23}
$$

which significantly alleviates condition (17) (by a factor of 3).
### 4.4. Two-Modes-Resonance Strategy

In this strategy, still with $\mathcal{N} \ge 2$, we consider alternatively two additional nonresonant fields $\Omega_{2\mathcal{N}-1}$ and $\Omega_{2\mathcal{N}+1}$ (and $\phi(t) = \omega_0 t$, i.e., $\Delta(t) = 0$) with the simplifying condition $\Omega_{2\mathcal{N}-1} = \Omega_{2\mathcal{N}+1} \equiv \Omega_0$, giving the total field

$$ \Omega(t) = \frac{\Omega_1}{2} \left[ 1 + \cos\left(\frac{\omega_0}{\mathcal{N}} t\right) \right] + \Omega_0 \left[ 1 + \cos\left(\frac{\omega_0}{\mathcal{N}} t\right) \cos(2\omega_0 t) \right]. \quad (24) $$

This field satisfies the zero-area condition (5) as before since $\Omega_{\mathcal{N}} = 0$. The Hamiltonian shows that the two fields generate a resonance:

$$ \begin{aligned} \tilde{H} ={}& \frac{1}{2} \left[ \Omega_1(t) + \Omega_0 + \frac{1}{2} \Omega_0 \cos\left(\frac{\omega_0}{\mathcal{N}} t\right) \right] \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} \\ & + \frac{1}{4} \Omega_0 \cos\left(\frac{\omega_0}{\mathcal{N}} t\right) \left( \begin{bmatrix} 0 & e^{2i\omega_0 t} \\ e^{2i\omega_0 t} & 0 \end{bmatrix} + \begin{bmatrix} 0 & e^{-2i\omega_0 t} \\ e^{-2i\omega_0 t} & 0 \end{bmatrix} \right) \\ & + \frac{1}{2} [\Omega_1(t) + \Omega_0] \begin{bmatrix} 0 & e^{-2i\omega_0 t} \\ e^{2i\omega_0 t} & 0 \end{bmatrix} + \frac{1}{4} \Omega_0 \cos\left(\frac{\omega_0}{\mathcal{N}} t\right) \begin{bmatrix} 0 & e^{-4i\omega_0 t} \\ e^{4i\omega_0 t} & 0 \end{bmatrix}. \end{aligned} \quad (25) $$

The area of the resulting coupling, which should be $\pi$, is

$$ \int_{-T/2}^{T/2} \Omega(t) dt = T\left(\frac{\Omega_1}{2} + \Omega_0\right) = \pi. \quad (26) $$

Neglecting the fastest oscillating term, we obtain for the effective Hamiltonian

$$ \tilde{H} \approx \frac{1}{2} \left[ \Omega_1(t) + \Omega_0 + \frac{1}{2} \Omega_0 \cos\left(\frac{\omega_0}{\mathcal{N}} t\right) \right] \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} + \frac{1}{8\omega_0} [\Omega_1(t) + \Omega_0]^2 \begin{bmatrix} -1 & 0 \\ 0 & 1 \end{bmatrix}, \quad (27) $$

i.e., for the peak values:

$$ \tilde{H}_{\max} \approx \frac{1}{2} \left(\Omega_1 + \frac{3}{2}\Omega_0\right) \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} + \frac{(\Omega_1 + \Omega_0)^2}{8\omega_0} \begin{bmatrix} -1 & 0 \\ 0 & 1 \end{bmatrix}. \quad (28) $$

The ratio of the (peak) diagonal correction with respect to the (peak) coupling is given by

$$ R = \frac{(\Omega_1 + \Omega_0)^2}{2\omega_0(\Omega_1 + \frac{3}{2}\Omega_0)} \ll 1. \quad (29) $$

This ratio is minimum when $\Omega_1 = 0$: $R = \Omega_0/(3\omega_0) = \pi/(3T\omega_0)$. On the contrary, when $\Omega_0 = 0$, the ratio is maximum. This solution minimizing the correction for the complete transfer is thus for $\Omega_1 = 0$ and $T\Omega_0 = \pi$, i.e.,

$$ \Omega_0 = \frac{\omega_0}{2\mathcal{N}}. \quad (30) $$

The diagonal terms can be neglected, from the above condition, when

$$ \mathcal{N} \gg 1/6. \quad (31) $$

This shows that the introduction of the two fields in the limit of weak $\Omega_1$ allows one to alleviate again condition (17) by a factor of 3, similarly to the preceding strategy.

### 5. Chirped Few-Cycle Pulses: Stark-Shift Compensation Strategy

The effective Hamiltonian (16) with a single-mode field $\Omega(t) = \Omega_1(t)$ indicates that the lowest order perturbative correction induces a dynamical (diagonal) Stark shift. It
can be compensated by a detuning properly shaped [according to (A55) if one considers additionally the time dependence of $\Omega_1(t)$]:

$$ \Delta(t) = -\frac{\Omega_1^2(t)}{4\omega_0} + \frac{\Omega_1^2}{64N^2\omega_0} \sin^2\left(\frac{\omega_0}{N}t\right) \approx -\frac{\omega_0}{32N^2}\left[3 + 4\cos\left(\frac{\omega_0}{N}t\right) + \cos\left(\frac{2\omega_0}{N}t\right)\right], \quad (32) $$

still with a $\pi$-pulse: $\Omega_1 T/2 = \pi$, where we have neglected the second term which is much smaller than the first one. The phase (2) can be then written

$$ \phi(t) = \left(1 + \frac{3}{32N^2}\right)\omega_0(t - t_i) + \frac{1}{8N}\left[\sin\left(\frac{\omega_0}{N}t\right) - \sin\left(\frac{\omega_0}{N}t_i\right)\right] \\ + \frac{1}{64N}\left[\sin\left(\frac{2\omega_0}{N}t\right) - \sin\left(\frac{2\omega_0}{N}t_i\right)\right] - \phi_{0,i}, \quad (33) $$

where the initial phase $\phi_{0,i}$ has to be chosen to satisfy the zero-area condition (5). We have checked that the dynamics weakly depends on its precise value (with a final population transfer oscillating between 0.9996 and 0.9998 as a function of the initial phase). One can prove that the constant part of the phase has to be $\pm \pi/2$ in order to satisfy (5) in this case of a single-mode field. The phase takes then the form (where the constant part of the phase $+\pi/2$ has been chosen)

$$ \phi(t) = \frac{\pi}{2} + \tilde{\omega}_0 t - \tilde{\phi}_0(t) \quad (34) $$

with, respectively, $\tilde{\omega}_0$ the mean frequency, which is the frequency that has to be tuned in practice, and $\tilde{\phi}_0(t)$ the modulated phase which is the part of the phase that has to be shaped:

$$ \tilde{\omega}_0 = \left(1 + \frac{3}{32N^2}\right)\omega_0, \quad \tilde{\phi}_0(t) = -\frac{1}{8N}\sin\left(\frac{\omega_0}{N}t\right)\left[1 + \frac{1}{4}\cos\left(\frac{\omega_0}{N}t\right)\right]. \quad (35) $$

Equation (35) shows a small shift of the mean frequency compared to the resonance $\omega_0$ and a small phase modulation amplitude (of approximately $1/8N$). For instance, we obtain a mean frequency $\tilde{\omega}_0 = 1.0234\omega_0$ and an approximate phase modulation amplitude of 0.0625 for $N=2$. Figure 8 shows the phase shaping and the dynamics in this case resulting into the transfer $P_+ \approx 0.9997$, close to an ultrahigh-fidelity transfer (and more than 30 times more accurate than the optimal nonchirped pulse shown in Figure 5).

More generally, we can show that the phase can take the following form

$$ \phi(t) = \frac{\pi}{2} + \tilde{\omega}_0 t - \sum_{n \ge 1} \tilde{\phi}_{0,n} \sin\left(\frac{n\omega_0}{N}t\right), \quad (36) $$

i.e., $\pi/2$ added to an odd function featuring modulating modes, with $\tilde{\omega}_0$ the mean frequency, in order to satisfy condition (5) with an even single-mode field $\Omega(t) = \Omega_1(t)$, for any integer number of cycles $N$. This will be used for a more systematic optimization in the next section.
<figure><img src="image_12.png" /><figcaption>**Figure 8.** High-fidelity population transfer for $N = 2$ with a single $\pi$-pulse mode $\Omega(t) = \Omega_1(t)$ and the phase (34) featuring a chirped frequency. From top to bottom frames: population dynamics, modulated phase $\tilde{\phi}_0(t)$, detuning $\Delta(t)$, and full coupling $\Omega_1(t) \cos \phi(t)$.</figcaption></figure>

# 6. Numerical Optimization

We provide the numerical results of a more systematic optimization procedure for the population transfer problem as described in Appendix B (see the flowchart of Figure A1). Since the parameter landscape is a priori large, we orient the search around the three strategies investigated in Sections 4 and 5 and check the convergence of our algorithm in each situation. We consider the number of cycles $N = 2$ for the strategies where the detuning is zero, which necessitates the cancellation of the $n = N$ mode. On the other hand, we investigate the single-cycle limit $N = 1$ for the more flexible situation with a time-dependent detuning (i.e., a chirped frequency).

## 6.1. 2N-Resonance Strategy

In order to target and validate this strategy, we impose the resonance, $\phi(t) = \omega_0 t$, and we implement an optimization using an expansion of $\Omega(t)$ limited to the first $4 = 2N$ modes for $N = 2$. The mode $n = 2$ is set to zero in order to satisfy Equation (5). Note that, without imposing it, the optimization procedure leads to the vanishing of the third mode $n = 3$, as predicted in Section 4.3. With this strategy, one obtains a complete transfer for $\Omega_1 \approx 0.277\omega_0$ and $\Omega_4 \approx 0.452\omega_0$, see Figure 9. These parameters are close to the ones determined in Section 4.3. The area of the full envelope $\Omega(t)$ gives $1.5\pi$. The area $\Omega_1(t) + \frac{1}{4}\Omega_2N$ is, as predicted, close to $\pi: T(\Omega_1 + \frac{1}{2}\Omega_4) \approx 1.006\pi$.
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Complete population transfer for $\mathcal{N} = 2$ and a four-mode Fourier decomposition with the nonzero modes $\Omega_1 = 0.277\omega_0$, $\Omega_4 = 0.452\omega_0$ (2N-resonance strategy) resulting from the optimization procedure. Upper frame: population dynamics. Middle frame: pulse envelope $\Omega(t)$. Lower frame: Full coupling $\Omega(t) \cos \omega_0 t$.</figcaption></figure>

Adding more modes, as proposed in the next subsection, will allow one to explore a different strategy to achieve the same complete transfer result.

## 6.2. Two-Modes-Resonance Strategy

We now consider the resonant case $\phi(t) = \omega_0 t$ with an expansion of $\Omega$ over the first 5 = 2N + 1 modes still for $\mathcal{N} = 2$. Adding the fifth mode opens a path to the two-modes-resonance strategy described in Section 4.4. The mode 2 is still set to zero in order to satisfy Equation (5) and, without imposing it, the optimization procedure leads to the vanishing of the fourth mode $n = 4$. One also obtains in this case a complete population transfer resulting from the combination of the $2N - 1 = 3$ and $2N + 1 = 5$ modes. We notice that, as predicted in Section 4.4, the optimized amplitudes $\Omega_3$ and $\Omega_5$ are close to each other (see Figure 10), with a value consistent with Equation (30).
<figure><img src="image_14.png" /><figcaption>Figure 10. Same as Figure 9 but for a five modes decomposition with the obtained non-zero modes $\Omega_1 = 0.043\omega_0$, $\Omega_3 = 0.2314\omega_0$, $\Omega_5 = 0.227\omega_0$ (two-modes-resonance strategy).</figcaption></figure>

### 6.3. Stark-Shift Compensation Strategy

One finally considers the situation with a nonzero detuning and a low number of cycles, $N = 2$ or $N = 1$. In both cases, it is sufficient to restrict $\Omega$ to the first Fourier mode, $\Omega(t) = \Omega_1(t)$, and to decompose the phase $\phi(t)$ according to Equation (36) with a single mode:

$$ \phi(t) = \frac{\pi}{2} + \tilde{\omega}_0 t - \tilde{\phi}_{0,1} \sin\left(\frac{\omega_0}{N} t\right). \qquad (37) $$

We obtain in both cases a complete transfer with smooth and simple phase and detuning.

For $N = 2$, the optimized parameters are $\Omega_1 = 0.504\omega_0$, $\tilde{\omega}_0 = 1.038\omega_0$, and $\tilde{\phi}_{0,1} = -9.74 \times 10^{-3}$, giving a pulse area $1.0087\pi$. The phase can be compared to (34) and (35): It features a larger shift of the mean frequency (though small compared to $\omega_0$) and a smaller modulation amplitude (given by the absolute value of $\tilde{\phi}_{0,1}$, which is small compared to the angle $\pi$), of frequency $\omega_0/2$. The resulting dynamics (not shown) is similar to the one obtained in Figure 8.

Optimization for the case $N = 1$ is displayed in Figure 11. It gives a pulse area $1.067\pi$. One can notice a stronger (and negative, i.e., $\dot{\phi} < \omega_0$) shift of the mean frequency (still much smaller than $\omega_0$) and a larger modulation amplitude compared to the above case $N = 2$ (but still smaller than the compensation shown in Figure 8), which is still small (compared to the angle $\pi$). The modulation frequency of the phase is here $\omega_0$.
In both cases, $\mathcal{N} = 1, 2$, the resulting phase is smooth, it features small deviations and small modulation amplitude at the frequency close to the resonant frequency $\omega_0$ and leads to an instantaneous frequency which is also close to $\omega_0$. All these ingredients indicate a realistic experimental implementation in principle.

<figure><img src="image_15.png" /><figcaption>**Figure 11.** Complete population transfer for $\mathcal{N} = 1$, a single-mode for the amplitude: $\Omega(t) = \Omega_1(t)$ and the phase: $\phi(t) = \pi/2 + \tilde{\omega}_0 t - \tilde{\phi}_{0,1} \sin(\omega_0 t)$ with the non-zero coefficients $\Omega_1 = 1.067\omega_0$, $\tilde{\omega}_0 = 0.889\omega_0$, $\tilde{\phi}_{0,1} = 3.56 \times 10^{-2}$. From top to bottom frames: population dynamics, modulated phase $\tilde{\phi}_0(t)$, detuning $\Delta(t)$, full coupling $\Omega_1(t) \cos \phi(t)$.</figcaption></figure>

7. Conclusions

In this paper, we have explored quantum control with a few-cycle pulse implementing and generalizing two standard strategies, adiabatic Floquet theory and $\pi$-pulse. We have shown that adiabatic Floquet is obstructed by the presence of avoided crossing in the quasi-energy spectrum induced by the counter-rotating term for the number of cycles $\mathcal{N} < 2.5$. This shows the relevance of using the tool of adiabatic Floquet theory for analyzing and controlling the dynamics even for few cycles.

We have implemented a generalization of $\pi$-pulse transfer by expanding the pulse shape in a Fourier expansion. We have shown two particular strategies on resonance with a few modes involved, named $2\mathcal{N}$-resonance and two-modes-resonance strategies. Complete transfer can be achieved in these cases for a number of cycles as low as $\mathcal{N} = 2$. We have also explored a chirped frequency associated to a simple pulsed-shape field. We have
shown that an appropriate chirping allows an ultrahigh-fidelity transfer for a number of cycles as low as $N = 1$ (single-cycle pulse).

In a future work, we will explore the control in a multilevel system since the broadening of the spectrum induced by the small number of cycles is expected to impact significantly excited states. The present study will serve as an important exploratory work that will guide the strategies with the clear advantage of limiting the parameters landscape compared to blind optimal control strategies. More specifically, we will use the few-cycle-pulse adiabatic Floquet theory, including multistate effects, such as rotational and vibrational states in molecules (similarly to Ref. [43]), but beyond the RWA (similarly to Figure 1), which will allow one to find adiabatic trajectories as functions of the phase and pulse amplitudes (when they exist), connecting the initial and targeted states. In addition, we have constructed specific parametrizations in amplitude (14) and (15) and in phase (36) that satisfy the physical constraint of zero time-integrated field area (5) with few parameters to be optimized. Nonadiabatic few-cycle regimes with shaped and chirped pulses in multi-level systems will be then explored on this basis. This will find for instance applications in the fine control of angular wavepackets towards the control of the orientation and more generally of the rotation of molecules [15] by single-cycle THz pulses.

**Author Contributions:** F.P. and S.G. carried out the numerical simulations. S.G. proposed the initial idea about the investigation of quantum control by few-cycles pulses. All authors worked on the design of the research project. All authors participated in the analysis of the results and in the paper redaction. All authors have read and agreed to the published version of the manuscript.

**Funding:** We acknowledge support from 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).

**Conflicts of Interest:** The authors declare no conflict of interest.

## Appendix A. Adiabatic Floquet theory

### Appendix A.1. Floquet Theory for Periodic Systems

The Floquet theory allows the treatment of a periodic time dependent system by transforming it into a time-independent problem via the use of an enlarged space as follows [18]. A periodic system of Hamiltonian $H(\theta + \phi(t) + 2\pi) = H(\theta + \phi(t))$ features a $2\pi$ periodic phase, $\phi(t) = \omega(t - t_i)$ with $\theta$ the initial phase at time $t_i$, $\omega$ the angular frequency. A typical example is

$$H(x; \theta + \omega(t - t_i)) = H_0(x) - \mu(x)\mathcal{E}\cos(\theta + \omega(t - t_i)) \quad (A1)$$

for a quantum system of Hamiltonian $H_0(x)$ driven by a periodic field of amplitude $\mathcal{E}$ via the dipole moment $\mu(x)$. The Schrödinger equation

$$i\hbar \frac{\partial}{\partial t} \varphi = H(x; \theta + \omega(t - t_i)) \varphi, \quad \varphi \in \mathcal{H} \quad (A2)$$

is defined on a Hilbert space $\mathcal{H}$, which can be of infinite dimension (e.g., the space of square-integrable functions $\mathcal{H} = L_2(\mathbb{R}^n, d^n x)$, where $n$ is the number of the degrees of freedom of the molecule) or of finite dimension (e.g., in $N$-level models $\mathcal{H} = \mathbb{C}^N$). The initial phase $\theta$ appears as a parameter. One can think of Equation (A2) as a family of equations parameterized by the angle $\theta$. We define the solution $\psi(t)$ of a time independent problem:

$$i\hbar \frac{\partial}{\partial t} \psi = K\psi, \quad \psi \in \mathcal{K} \quad (A3)$$

with the time-independent Floquet Hamiltonian (or quasienergy operator)

$$K = -i\hbar\omega \frac{\partial}{\partial \theta} + H(\theta), \quad (A4)$$
defined in an enlarged Hilbert space

$$ \mathcal{K} := \mathcal{H} \otimes \mathcal{L}, \qquad (A5) $$

where $\mathcal{L} := L_2(\mathbb{S}^1, d\theta/2\pi)$ denotes the space of square integrable functions on the circle $\mathbb{S}^1$ of length $2\pi$, with the scalar product

$$ \langle \xi_1 | \xi_2 \rangle_{\mathcal{L}} := \int_{\mathbb{S}^1} \frac{d\theta}{2\pi} \xi_1^*(\theta) \xi_2(\theta). \qquad (A6) $$

The solution $\varphi(t)$ is related to the time independent problem as

$$ \varphi(x; t; \theta) = \psi(x, \theta + \omega(t - t_i); t) \qquad (A7) $$

with the initial condition in $\mathcal{K}$ constructed as the initial condition $\varphi(t_i)$ of the original problem lifted in the enlarged space as:

$$ \psi(x, \theta; t_i) = \varphi(x; t_i; \theta) \otimes \mathbb{1}_{\mathcal{L}} \qquad (A8) $$

with $\mathbb{1}_{\mathcal{L}} \equiv e^{i(k=0)\theta}$. In the expression $\varphi(x; t; \theta)$, on the left-hand side of (A7), $\theta$ is a parameter corresponding to the initial phase of the Hamiltonian, while, in $\psi(x, \theta; t_i)$, on the right-hand side, $\theta$ is a dynamical variable on the same footing as $x$. The procedure dictated by Equation (A7) to obtain the dynamics in the original space $\mathcal{H}$ is as follows:

(i) first, by the use of Equation (A8), lift the initial condition in the enlarged space $\mathcal{K}$ keeping $\theta$ undetermined;

(ii) next determine $\psi(x, \theta; t)$ in $\mathcal{K}$ solution of Equation (A3);

(iii) finally replace the angle $\theta$ by $\theta + \omega(t - t_i)$ in $\psi$ and fix the angle $\theta$ to the initial phase of the Hamiltonian.

Since $\mathcal{K}$ is time-independent, the convenient way to compute the solution is to use an eigenfunction expansion

$$ \varphi(x; t; \theta) = \sum_{\nu} c_{\nu} e^{-i\lambda_{\nu}(t-t_i)/\hbar} \psi_{\nu}(x, \theta + \omega(t-t_i)) \qquad (A9) $$

in terms of the eigenelements of $\mathcal{K}$ (assuming a discrete spectrum)

$$ K \psi_{\nu}(x, \theta) = \lambda_{\nu} \psi_{\nu}(x, \theta), \qquad (A10) $$

with the eigenvalues named quasienergies, with

$$ c_{\nu} = \langle \psi_{\nu} | \varphi(t_i) \otimes 1 \rangle_{\mathcal{K}} = \langle \tilde{\psi}_{\nu} | \varphi(t_i) \otimes 1 \rangle_{\mathcal{H}}, \qquad (A11) $$

where $\tilde{\psi}_{\nu}(x) := \int_{\mathbb{S}^1} d\theta/2\pi \psi_{\nu}(x, \theta)$ is the average of $\psi_{\nu}(x, \theta)$ over the phase, or equivalently, its constant Fourier component.

The Floquet spectrum contains an infinite number of elements, but its eigenelements feature a periodic structure with the index $\nu$ being associated to two indices, $\nu \equiv n, k$:

$$ \psi_{n,k}(x, \theta) = \psi_{n,0}(x, \theta) e^{ik\theta}, \quad \lambda_{n,k} = \lambda_{n,0} + k\hbar\omega, \qquad (A12) $$

where the index $n$ refers to the molecule's Hilbert space $\mathcal{H}$ (i.e., $n = 1, \dots, N$ if $\mathcal{H} = \mathbb{C}^N$), and $k$ are all positive or negative integers. This allows one to classify the Floquet eigenstates in families labeled by $n$. The individual members within one family are distinguished by the index $k$. The eigenfunction expansion can be simplified using only one representative of each family (e.g., the one with $k = 0$):

$$ \varphi(x; t; \theta) = \sum_{n} \tilde{c}_{n}(\theta) e^{-i\lambda_{n,0}(t-t_i)/\hbar} \psi_{n,0}(x, \theta + \omega(t-t_i)) \qquad (A13) $$
with

$$ \tilde{c}_n(\theta) := \langle \psi_{n,0}(\theta) | \varphi(t_i) \rangle_{\mathcal{H}}. \qquad (A14) $$

The interpretation of this expansion is as follows: the initial condition is expanded into the set of representatives of each family of eigenvectors $\psi_{n,0}(x, \theta)$ through the coefficient $\tilde{c}_n(\theta)$ (which depends parametrically on $\theta$); during the evolution, the expansion acquires dynamical phases through the eigenvalues $\lambda_{n,0}$ and features a dynamical evolution of $\theta$ in the eigenvectors.

## Appendix A.2. Adiabatic Floquet Theory

The adiabatic Floquet theory extends the preceding formulation when the Hamiltonian features an additional non-periodic timescale, which is typically slow with respect to the periodic one. These parameters are gathered under the function $r(t)$ (for instance the time-dependent amplitude of the coupling for the case of a pulse-shaped field and its time-dependent polarization if an interacting polarized laser field is considered), while the non-periodic time-dependence of the phase itself (referred to as chirped pulse) is considered specifically:

$$ H^{\mathbf{r}(t)}(\theta + \varphi(t) + 2\pi) = H^{\mathbf{r}(t)}(\theta + \varphi(t)). \qquad (A15) $$

The $2\pi$-periodic phase $\theta + \varphi(t)$ features an instantaneous (or effective) frequency $\omega(t) = \dot{\varphi}(t)$, such that $\varphi(t) = \int_{t_i}^t \omega(s)ds$, allows one to recover the phase dependence of the previous periodic system when the time of the frequency is frozen.

The state $\varphi(x; t; \theta)$ solution of the Schrödinger equation

$$ i\hbar \frac{\partial}{\partial t} \varphi = H^{\mathbf{r}(t)}(\theta + \varphi(t)) \varphi \qquad (A16) $$

is connected to $\psi(x, \theta; t)$ the solution of

$$ i\hbar \frac{\partial}{\partial t} \psi = K^{\mathbf{r}(t), \omega(t)} \psi, \qquad (A17) $$

with

$$ \omega(t) = \dot{\varphi}(t), \quad \varphi(t) = \int_{t_i}^{t} \omega(s)ds \qquad (A18) $$

and

$$ K^{\mathbf{r}(t), \omega(t)} = H^{\mathbf{r}(t)}(\theta) - i\hbar\omega(t) \frac{\partial}{\partial\theta} \qquad (A19) $$

through

$$ \varphi(x; t; \theta) = \psi(x, \theta + \varphi(t); t), \qquad (A20) $$

with the initial condition

$$ \psi(x, \theta; t_i) = \varphi(x; t_i; \theta) \otimes \mathbb{1}_{\mathcal{L}}. \qquad (A21) $$

We remark that, at this stage, this correspondence between the original and Floquet dynamics with additional non-periodic time dependent parameters is exact. It does not require a specific number of oscillations of the field in the pulse, nor specific (slow) variation of these parameters.

The physical interpretation of the operator $-i\hbar \frac{\partial}{\partial \theta}$ is the relative photon number operator in the limit of large number of photons present in the interacting field. One can thus assign two labels to the eigenstates $|\psi_{n,k}\rangle$, with the label $n$ coming from the matter and the label $k$ interpreted as the relative photon number which *dresses* the state. A transfer from the initial eigenstate $|\psi_{n,k}^0\rangle \equiv |n\rangle \otimes |k\rangle$ with $|n\rangle$ a bare state of the quantum system and $|k\rangle \equiv e^{ik\theta}$ a (relative) state of $k$ photons before the field is on, i.e., before interaction between the matter and the field, to the eigenstate $|\psi_{n',k'}^0\rangle \equiv |n'\rangle \otimes |k'\rangle$ after the field is off, i.e., after the interaction, corresponds to the transfer from $|n\rangle$ to $|n'\rangle$ with the absorption (emission) of
$k' - k$ photons for negative (positive) $k' - k$. When the field is on, the Floquet eigenstates cannot be simply decomposed as a tensor product of a bare matter state and a field state, but are rather an entangled state, as a superposition of such tensored states.

In the limit of slow variations, $t \equiv \epsilon s$ with $\epsilon \ll 1$, one can apply the above formulation with the use the eigenelements of $K$, which depends now on the parameters:

$$K^{\mathbf{r}(t),\omega(t)} \psi_{\mathbf{v}}^{\mathbf{r}(t),\omega(t)}(x, \theta) = \lambda_{\mathbf{v}}^{\mathbf{r}(t),\omega(t)} \psi_{\mathbf{v}}^{\mathbf{r}(t),\omega(t)}(x, \theta), \quad (A22)$$

on which we derive an adiabatic approximation

$$\varphi(x; t; \theta) \approx \sum_n \tilde{c}_n(\theta) e^{-i \int_{t_i}^t ds \lambda_{n,k(n)}^{(s),\omega(s)}/\hbar} \psi_{n,k(n)}^{\mathbf{r}(t),\omega(t)}(x, \theta + \varphi(t)) \quad (A23)$$

with

$$\tilde{c}_n(\theta) := \langle \psi_{n,k(n)}^{\mathbf{r}(t_i),\omega(t_i)}(\theta) | \varphi(t_i) \rangle_{\mathcal{H}}. \quad (A24)$$

Similarly as before, this expansion means that the initial condition is expanded into the set composed of one representative of all the families of eigenvectors $\psi_{n,k(n)}^{\mathbf{r}(t_i),\omega(t_i)}(x, \theta)$ at the initial values of the parameters $\mathbf{r}(t_i), \omega(t_i)$ through the coefficient $\tilde{c}_n(\theta)$ (which depends parametrically on $\theta$); during the evolution, the expansion acquires dynamical phases through the integration of the eigenvalues $\lambda_{n,0}^{\mathbf{r}(t),\omega(t)}$ via the time-dependence of the parameters $\mathbf{r}(t), \omega(t)$ and features a dynamical evolution of $\theta$ and of the parameters in the eigenvectors. The choice of the representatives labelled by $k(n)$ is such that their corresponding eigenvalues are localized in a band of energy of width $\hbar\omega$ defining a Floquet zone; the coefficients $k(n)$ depends thus in general on $n$.

A remarkable property is when a single eigenstate is followed by the dynamics, since, then, the resulting single dynamical phase is an irrelevant global phase. Moreover, the control of the dynamics reduces to an analysis of the trajectories of the eigenvalues and, more precisely, to the connection of a trajectory between an initial and a final state [17]. This means that if a state is targeted at a specific time, for instance at the end of the interaction, one just has to drive the dynamics along an adiabatic trajectory, which connects the initial and the desired targeted states; the exact form of the trajectory being unimportant. The dynamics can be thus selective and robust with respect to the precise instantaneous values of the parameters.

The (local) correction of the eigenfunction expansion (A23) is in general $O(\epsilon)$. However it can be made in principle exponentially small, $O(\epsilon^{-C/\epsilon})$ with C a positive constant, with the use of smooth controls and a global analysis of the dynamics [38]. Obstruction of adiabatic passage occurs when the eigenvalues get closer, which can induce non-adiabatic transitions between the eigenstates, as typically given by Landau-Zener avoided crossing. Adiabatic passage can be made more efficient when the eigenvalues are parallel to each other during the dynamics [39–41].

## Appendix B. Algorithms

Our numerical simulations were implemented under Matlab using the built-in function `fmincon`. This command is dedicated to finding the minimum of a constrained non-linear multivariable function.

Our script is based on three loops as shown in the workflow in Figure A1. The first loop is dedicated to the initialization of the Fourier coefficients of $\Omega(t)$ under the constraints given by the RWA $\pi$-pulse conditions. The second loop explores a larger space of solutions by removing the RWA $\pi$-pulse constraint. Finally the last loop realizes the optimization by adding the detuning.

Note that the condition (5) is imposed in the first two loops by choosing $\Omega_N = 0$ in the Fourier expansion of $\Omega(t)$. In the third loop, this condition (5) is satisfied by the specific form of the phase (36).
<figure><img src="image_16.png" /><figcaption>Figure A1. Flowchart of the optimization algorithm. The modified Crank-Nicolson scheme can be found in Ref. [44].</figcaption></figure>

## Appendix C. Perturbation Theory Formulated with KAM Techniques

We decompose the Hamiltonian with an unperturbed Hamiltonian $H_0$ and a perturbation $\epsilon V_1$:

$$H = H_0 + \epsilon V_1. \qquad (A25)$$
We construct a unitary transformation (so-called KAM or contact transformation) $T_1 = e^{\epsilon W_1}$, with $W_1^\dagger = -W_1$ such that

$$e^{-\epsilon W_1} H e^{\epsilon W_1} = H_0 + \epsilon D_1 + \epsilon^2 V_2, \qquad (A26)$$

where $D_1$ is a diagonal operator, i.e., satisfying $[H_0, D_1] = 0$. Thus the perturbation is reduced from order $\epsilon$ to order $\epsilon^2$.

In order to determine $D_1$ and $W_1$, one expands Equation (A26) in powers of $\epsilon$ and require that the terms of order $\epsilon$ cancel out:

$$ \begin{aligned} e^{-\epsilon W_1} H e^{\epsilon W_1} &= (\mathbf{1} - \epsilon W_1 + \frac{1}{2}\epsilon^2 W_1^2 + \dots)(H_0 + \epsilon V_1)(\mathbf{1} + \epsilon W_1 + \frac{1}{2}\epsilon^2 W_1^2 + \dots) \\ &= H_0 + \epsilon([H_0, W_1] + V_1) + \frac{1}{2}\epsilon^2([H_0, W_1], W_1) + \epsilon^2[V_1, W_1] + \dots \\ &= H_0 + \epsilon D_1 + \frac{1}{2}\epsilon^2[V_1, W_1] + \frac{1}{2}\epsilon^2[D_1, W_1] + \dots \end{aligned} \qquad (A27) $$

This leads to the two equations

$$[H_0, W_1] + V_1 - D_1 = 0, \qquad (A28a)$$

$$[H_0, D_1] = 0, \qquad (A28b)$$

and the first terms of the resulting second order perturbation read:

$$\epsilon^2 V_2 = \frac{\epsilon^2}{2} [V_1, W_1] + \frac{\epsilon^2}{2!} [D_1, W_1] + \dots \qquad (A29)$$

The solution of Equation (A28) can be given using the eigenvalues and eigenvectors of $H_0$ which we will denote by $\lambda_\nu^0$ and $|\nu, j\rangle$ (we use an index $\nu$ that labels the different eigenvalues, and $j$ distinguishes different basis vectors corresponding to a degenerate eigenvalue). We define a projection operator $\Pi_{H_0}$ that extracts from the perturbation $V_1$ the diagonal component with respect to the eigenbasis of $H_0$:

$$\Pi_{H_0} V_1 \equiv \bar{V}_1 = \sum_{\nu, j, j'} |\nu, j\rangle \langle \nu, j| V_1 |\nu, j'\rangle \langle \nu, j'|. \qquad (A30)$$

With this notation a solution of (A28) can be written as

$$D_1 = \bar{V}_1 = \text{diagonal part of } V_1, \qquad (A31a)$$

$$W_1 = - \sum_{\nu, j, j', \nu' \neq \nu} \frac{|\nu, j\rangle \langle \nu, j| V_1 |\nu', j'\rangle \langle \nu', j'|}{\lambda_{\nu}^0 - \lambda_{\nu'}^0}. \qquad (A31b)$$

$D_1$ is interpreted as the averaging of $V_1$ with respect to $H_0$. The perturbation theory does not apply when a resonance appears, which is detected by a degeneracy $\lambda_\nu^0 = \lambda_{\nu'}^0$, leading to a zero denominator of $W_1$.

The KAM transformations can be iterated. One usually applies another KAM transformation consisting in extracting the averaging of $V_2$, and the effective Hamiltonian of second order reads

$$H_{\text{eff}} = H_0 + \epsilon \bar{V}_1 + \epsilon^2 \bar{V}_2. \qquad (A32)$$
## Appendix D. High-Frequency Perturbation Theory

High frequency perturbation theory applies to situations in which the frequency $\omega$ of the perturbation is high with respect to the internal frequencies of the considered system. We formulate with a Floquet Hamiltonian of the form

$$K = -i\hbar\omega \frac{\partial}{\partial\theta} + H_0 + V_1(\theta). \qquad (A33)$$

Since we are interested in the limit $\hbar\omega \to \infty$, we define a small parameter $\epsilon := 1/(\hbar\omega)$ and we rewrite

$$K = \hbar\omega\hat{K} \qquad (A34)$$

with

$$\hat{K} = -i \frac{\partial}{\partial \theta} + \epsilon(H_0 + V_1). \qquad (A35)$$

The difference with the preceding section is that here $H_0$ and $V$ are both of order $\epsilon$. Thus we take as the unperturbed Floquet Hamiltonian just

$$\hat{K}_0 := -i \frac{\partial}{\partial \theta}. \qquad (A36)$$

If the frequency is large compared with the frequencies of the system, there will not be any resonances. We can thus proceed with the perturbative KAM transformation: $e^{\epsilon W_1(\theta)}$, with $W_1^\dagger = -W_1$ such that

$$e^{-\epsilon W_1} \hat{K} e^{\epsilon W_1} = \hat{K}_0 + \epsilon D_1 + \epsilon^2 V_2, \qquad (A37)$$

where $D_1$ is a $\theta$-independent operator such that $[\hat{K}_0, D_1] = 0$. The generator $W_1$ of the contact transformation is determined by the equations

$$[K_0, W_1] + H_0 + V_1 - D_1 = 0, \qquad (A38a)$$

$$[K_0, D_1] = 0. \qquad (A38b)$$

Equation (A38a) can be written as

$$-i \frac{\partial W_1}{\partial \theta} + H_0 + V_1 - D_1 = 0, \qquad (A39)$$

whose general solution is given by

$$W_1 = -i \int^{\theta} d\theta (H_0 + V_1 - D_1), \qquad (A40)$$

where we have removed an arbitrary constant. Since $W_1(\theta)$ is a multiplication operator acting on functions of the angle $\theta$, it must be necessarily $2\pi$-periodic. This condition determines $D_1$ uniquely in terms of averaging

$$\bar{V}_1 := \int_0^{2\pi} \frac{d\theta}{2\pi} V_1(x, \theta) \qquad (A41)$$

as

$$D_1 = H_0 + \bar{V}_1. \qquad (A42)$$

Thus we obtain

$$W_1(\theta) = -i \int^{\theta} d\theta (V_1(\theta) - \bar{V}_1). \qquad (A43)$$
The remaining perturbation of order $\epsilon^2$ can be written as

$$ \epsilon^2 V_2 = \frac{\epsilon^2}{2!} [V_1 + \bar{V}_1 + 2H_0, W_1] + \epsilon^3 \dots, \qquad (A44) $$

We can apply a second contact transformation $e^{\epsilon^2 W_2}$ (with respect to $\hat{K}_0$), $W_2 = -i \int^\theta d\theta (V_2 - \bar{V}_2)$, which averages this rest (A44) with respect to $\theta$ and leads to correction of order $\epsilon^3$. We thus obtain the effective high frequency Hamiltonian $H^{\text{HF}}$ (independent of $\theta$) of second order

$$ H^{\text{HF}} = H_0 + \bar{V}_1 + \epsilon \bar{V}_2. \qquad (A45) $$

As an example, one considers an Hamiltonian of the form:

$$ H = \frac{1}{2}\Omega_0 \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} + \frac{1}{2}\Omega_{-2} \begin{bmatrix} 0 & e^{-2i\omega t} \\ e^{2i\omega t} & 0 \end{bmatrix} + \frac{1}{2}\Omega_2 \begin{bmatrix} 0 & e^{2i\omega t} \\ e^{-2i\omega t} & 0 \end{bmatrix}, \qquad (A46) $$

i.e.,

$$ \hat{K} = \frac{K}{\omega} = -i \frac{\partial}{\partial \theta} + \frac{1}{2} \frac{\Omega_0}{\omega} \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} + \frac{1}{2} \frac{\Omega_{-2}}{\omega} \begin{bmatrix} 0 & e^{-2i\theta} \\ e^{2i\theta} & 0 \end{bmatrix} + \frac{1}{2} \frac{\Omega_2}{\omega} \begin{bmatrix} 0 & e^{2i\theta} \\ e^{-2i\theta} & 0 \end{bmatrix}, \qquad (A47) $$

where

$$ \epsilon H_0 \equiv \frac{1}{2} \frac{\Omega_0}{\omega} \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}, \quad \epsilon V_1 \equiv \frac{1}{2} \frac{\Omega_{-2}}{\omega} \begin{bmatrix} 0 & e^{-2i\theta} \\ e^{2i\theta} & 0 \end{bmatrix} + \frac{1}{2} \frac{\Omega_2}{\omega} \begin{bmatrix} 0 & e^{2i\theta} \\ e^{-2i\theta} & 0 \end{bmatrix}, \qquad (A48) $$

in the high frequency limit, rewritten with the dimensionless coefficients $\Omega_{-2}/2\omega \to 0$, $\Omega_2/2\omega \to 0$, $\Omega_0/2\omega \to 0$. We obtain $\bar{V}_1 = 0$,

$$ W_1(\theta) = \frac{1}{4}\Omega_{-2} \begin{bmatrix} 0 & e^{-2i\theta} \\ -e^{2i\theta} & 0 \end{bmatrix} + \frac{1}{4}\Omega_2 \begin{bmatrix} 0 & -e^{2i\theta} \\ e^{-2i\theta} & 0 \end{bmatrix}, \qquad (A49) $$

and $\bar{V}_2 = \frac{1}{2}[\bar{V}_1, W_1] + [H_0, W_1] + \dots$, where

$$ [H_0, W_1] = 0, \quad [\bar{V}_1, W_1] = \frac{1}{4}(\Omega_{-2}^2 - \Omega_2^2) \begin{bmatrix} -1 & 0 \\ 0 & 1 \end{bmatrix}, \qquad (A50) $$

i.e., at the second order

$$ \bar{V}_2 = \frac{1}{8}(\Omega_{-2}^2 - \Omega_2^2) \begin{bmatrix} -1 & 0 \\ 0 & 1 \end{bmatrix}. \qquad (A51) $$

The high-frequency effective Hamiltonian reads then

$$ H^{\text{HF}} = \frac{1}{2} \begin{bmatrix} -\frac{\Omega_{-2}^2 - \Omega_2^2}{4\omega} & \Omega_0 \\ \Omega_0 & \frac{\Omega_{-2}^2 - \Omega_2^2}{4\omega} \end{bmatrix}. \qquad (A52) $$

The diagonal term can be neglected when

$$ \Omega_0 \gg \frac{1}{2} \frac{|\Omega_{-2}^2 - \Omega_2^2|}{\omega}. \qquad (A53) $$

If we consider a time-dependent $\Omega(t)$, the transformation $T_1 = e^{\epsilon W_1}$ introduces in the Schrödinger equation an additional term of the form:

$$ -i\epsilon T^\dagger T = -i\epsilon^2 \hat{W}_1 = -i \frac{\dot{\Omega}}{4\omega^2} \begin{bmatrix} 0 & e^{-2i\theta} \\ -e^{2i\theta} & 0 \end{bmatrix}, \qquad (A54) $$
which should be treated perturbatively like $\epsilon V_1$. It introduces an additional term in the diagonal of (A52):

$$H^{\text{HF}} = \frac{1}{2} \left[ -\frac{\Omega_2^2 - \Omega_2^2}{4\omega} + \frac{\Omega_2^2 - \Omega_2^2}{16\omega^3} \frac{\Omega_0}{\Omega_0^2 - \Omega_2^2} - \frac{\Omega_2^2 - \Omega_2^2}{16\omega^3} \right]. \quad (A55)$$

## References

1. Shapiro, M.; Brumer, P. *Quantum Control of Molecular Processes*; Wiley-VCH: Weinheim, Germany, 2012.

2. Shore, B.W. *Manipulating Quantum Structure Using Laser Pulses*; Cambridge University Press: Cambridge, UK, 2011.

3. Bergmann, K.; Theuer, H.; Shore, B.W. Coherent population transfer among quantum states of atoms and molecules. *Rev. Mod. Phys.* **1998**, *70*, 1003. [CrossRef]

4. Nisbet-Jones, P.B.R.; Dilley, J.; Ljunggren, D.; Kuhn, A. Highly efficient source for indistinguishable single photons of controlled shape. *New J. Phys.* **2011**, *13*, 103036. [CrossRef]

5. Rousseaux, B.; Dzsotjan, D.; Francs, G.C.d.; Jauslin, H.R.; Couteau, C.; Guérin, S. Adiabatic passage mediated by plasmons: A route towards a decoherence-free quantum plasmonic platform. *Phys. Rev. B* **2016**, *93*, 045422; Erratum in: *Phys. Rev. B* **2016**, *94*, 199902. [CrossRef]

6. Dorier, V.; Gevorgyan, M.; Ishkhanyan, A.; Leroy, C.; Jauslin, H.R.; Guérin, S. Non-linear stimulated Raman exact passage by resonance-locked inverse engineering. *Phys. Rev. Lett.* **2017**, *119*, 243902. [CrossRef]

7. Dupont, N.; Chatelain, G.; Gabardos, L.; Arnal, M.; Billy, J.; Peaudecerf, B.; Sugny, D.; Guéry- Odelin, D. Quantum State Control of a Bose-Einstein Condensate in Optical Lattice. *PRX Quantum* **2021**, *2*, 040303. [CrossRef]

8. Motzoi, F.; Gambetta, J.; Rebentrost, P.; Wilhelm, F.K. Simple pulses for elimination of leakage in weakly nonlinear qubits. *Phys. Rev. Lett.* **2009**, *103*, 110501. [CrossRef] [PubMed]

9. Vepsäläinen, A.; Danilin, S.; Paraoanu, G.S. Superadiabatic population transfer in a three-level superconducting circuit. *Sci. Adv.* **2019**, *5*, eaau5999. [CrossRef]

10. Nielson, M.A.; Chuang, I.L. *Quantum Computation and Quantum Information*; Cambridge University Press: Cambridge, UK, 2000.

11. Sangouard, N.; Lacour, X.; Guérin, S.; Jauslin, H.R. Fast SWAP gate by adiabatic passage. *Phys. Rev. A* **2005**, *72*, 062309. [CrossRef]

12. Jones, J.A. Designing short robust NOT gates for quantum computation. *Phys. Rev. A* **2013**, *87*, 052317. [CrossRef]

13. Rousseaux, B.; Guérin, S.; Vitanov, N.V. Arbitrary qudit gates by adiabatic passage. *Phys. Rev. A* **2013**, *87*, 032328. [CrossRef]

14. Koch, C.P.; Lemeshko, M.; Sugny, D. Quantum control of molecular rotation. *Rev. Mod. Phys.* **2019**, *91*, 035005. [CrossRef]

15. Lapert, M.; Guérin, S.; Sugny, D. Field-free quantum cogwheel by shaping of rotational wave packets. *Phys. Rev. A* **2011**, *83*, 013403. [CrossRef]

16. Vitanov, N.V.; Rangelov, A.A.; Shore, B.W.; Bergmann, K. Stimulated Raman adiabatic passage in physics, chemistry, and beyond. *Rev. Mod. Phys.* **2017**, *89*, 015006. [CrossRef]

17. Yatsenko, L.P.; Guérin, S.; Jauslin, H.R. Topology of adiabatic passage. *Phys. Rev. A* **2002**, *65*, 043407. [CrossRef]

18. Guérin, S.; Jauslin, H.R. Control of Quantum Dynamics by Laser Pulses: Adiabatic Floquet Theory. *Adv. Chem. Phys.* **2003**, *125*, 147.

19. Levitt, M.H. *Spin Dynamics: Basics of Nuclear Magnetic Resonance*; John Wiley and Sons: New York, NY, USA, 2008.

20. Torosov, B.T.; Guérin, S.; Vitanov, N.V. High-fidelity adiabatic passage by composite sequences of chirped pulses. *Phys. Rev. Lett.* **2011**, *106*, 233001. [CrossRef]

21. Genov, G.T.; Schraft, D.; Halfmann, T.; Vitanov, N.V. Correction of Arbitrary Field Errors in Population Inversion of Quantum Systems by Universal Composite Pulses. *Phys. Rev. Lett.* **2014**, *113*, 043001. [CrossRef]

22. Koch, C.P.; Boscain, U.; Calarco, T.; Dirr, G.; Filipp, S.; Glaser, S.J.; Kosloff, R.; Montangero, S.; Schulte-Herbrüggen, T.; Sugny, D.; et al. Quantum optimal control in quantum technologies. Strategic report on current status, visions and goals for research in Europe. *EPJ Quantum Technol.* **2022**, *9*, 19. [CrossRef]

23. Boscain, U.; Sigalotti, M.; Sugny, D. Introduction to the Pontryagin Maximum Principle for Quantum Optimal Control. *PRX Quantum* **2021**, *2*, 030203. [CrossRef]

24. Boscain, U.; Charlot, G.; Gauthier, J.-P.; Guérin, S.; Jauslin, H.R. Optimal control in laser-induced population transfer for two- and three-level quantum systems. *J. Math. Phys.* **2002**, *43*, 2107. [CrossRef]

25. Daems, D.; Ruschhaupt, A.; Sugny, D.; Guérin, S. Robust Quantum Control by a Single-Shot Shaped Pulse. *Phys. Rev. Lett.* **2013**, *111*, 050404. [CrossRef] [PubMed]

26. Dridi, G.; Liu, K.; Guérin, S. Optimal robust quantum control by inverse geometric optimization. *Phys. Rev. Lett.* **2020**, *125*, 250403. [CrossRef] [PubMed]

27. Chen, X.; Lizuain, I.; Ruschhaupt, A.; Guéry-Odelin, D.; Muga, J.G. Shortcut to Adiabatic Passage in Two- and Three-Level Atoms. *Phys. Rev. Lett.* **2010**, *105*, 123003. [CrossRef]

28. Guéry-Odelin, D.; Ruschhaupt, A.; Kiely, A.; Torrontegui, E.; Martínez-Garaot, S.; Muga, J.G. Shortcuts to adiabaticity: Concepts, methods, and applications. *Rev. Mod. Phys.* **2019**, *91*, 045001. [CrossRef]
29. Ding, Y.-C.; Ban, Y.; Martín-Guerrero, J.D.; Solano, E.; Casanova, J.; Chen, X. Breaking adiabatic quantum control with deep learning. *Phys. Rev. A* **2021**, *103*, L040401. [CrossRef]

30. Giannelli, L.; Sgroi, P.; Brown, J.; Paraoanu, G.S.; Paternostro, M.; Paladino, E.; Falci, G. A tutorial on optimal control and reinforcement learning methods for quantum technologies. *Phys. Lett. A* **2022**, *434*, 128054. [CrossRef]

31. Kärtner, F.X. *Few-Cycle Laser Pulse Generation and Its Applications*; Springer: Berlin/Heidelberg, Germany, 2004.

32. Brabec, T.; Krausz, F. Intense few-cycle laser fields: Frontiers of nonlinear optics. *Rev. Mod. Phys.* **2000**, *72*, 545. [CrossRef]

33. Egodapitiya, K.N.; Li, S.; Jones, R.R. Terahertz-Induced Field-Free Orientation of Rotationally Excited Molecules. *Phys. Rev. Lett.* **2014**, *112*, 103002. [CrossRef]

34. Xu, L.; Tutunnikov, I.; Gershnabel, E.; Prior, Y.; Averbukh, I.S. Long-Lasting Molecular Orientation Induced by a Single Terahertz Pulse. *Phys. Rev. Lett.* **2020**, *125*, 013201. [CrossRef]

35. Shu, C.-C.; Hong, Q.-Q.; Guo, Y.; Henriksen, N.E. Orientational quantum revivals induced by a single-cycle terahertz pulse. *Phys. Rev. A* **2020**, *102*, 063124. [CrossRef]

36. Shore, B.W. *The Theory of Coherent Atomic Excitation*; Wiley: New York, NY, USA, 1990.

37. Pollock, C.R. Ultrafast optical pulses. In *Progress in Optics*; Wolf, E., Ed.; Elsevier: Amsterdam, The Netherlands, 2008; pp. 211–249.

38. Berry, M.V. Quantum phase corrections from adiabatic iteration. *Proc. R. Soc. Lond. A* **1987**, *414*, 31.

39. Guérin, S.; Thomas, S.; Jauslin, H.R. Optimization of population transfer by adiabatic passage. *Phys. Rev. A* **2002**, *65*, 023409. [CrossRef]

40. Dridi, G.; Guérin, S.; Hakobyan, V.; Jauslin, H.R.; Eleuch, H. Ultrafast stimulated Raman parallel adiabatic passage by shaped pulses. *Phys. Rev. A* **2009**, *80*, 043408. [CrossRef]

41. Guérin, S.; Hakobyan, V.; Jauslin, H.R. Optimal adiabatic passage by shaped pulses: Efficiency and robustness. *Phys. Rev. A* **2011**, *84*, 01343. [CrossRef]

42. Vasilev, G.S.; Vitanov, N.V. Complete population transfer by a zero-area pulse. *Phys. Rev. A* **2006**, *73*, 023416. [CrossRef]

43. Thomas, S.; Guérin, S.; Jauslin, H.R. State-selective chirped adiabatic passage on dynamically laser-aligned molecules. *Phys. Rev. A* **2005**, *71*, 013402. [CrossRef]

44. Hong, J.; Liu, Y.; Munthe-Kaas, H.; Zanna, A. Globally conservative properties and error estimation of a multi-symplectic scheme for Schrödinger equations with variable coefficients. *Appl. Numer. Math.* **2006**, *56*, 814–843. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.