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

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 π-pulse regime to two- or single-cycle pulses, are also derived.


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][11][12][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][20][21], optimal control [22][23][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 where Ω(t) is the Rabi frequency defined as Ω(t) = −µ 01 E (t)/h (considered positive without loss of generality) with the dipole moment µ 01 coupling the two states; ω 0 > 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 φ 0 (t) with respect to the Bohr frequency: and the corresponding detuning ∆: with t i = −T/2 the initial time, t f = −t i the final time with full duration T = t f − t i , and the initial phase φ 0,i . We will restrict our study to the near-resonant caseφ ω 0 for the ease of implementation. The number N of oscillations in the pulse (with respect to the resonant frequency ω 0 in this near-resonant case) is such that The physical implementation of the few-cycle field imposes its zero time-integrated area [37]: which will be imposed for all the derived fields in this work. After application of the resonant (or rotating wave) transformation: we obtain (where we have for convenience shifted the origin of energy such that a traceless matrix is derived) 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 Ω ω 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 and |+, −1 , 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 |+ to state |− (or reciprocally can absorb a photon from state |− to state |+ ).
Taking the simple example of a resonant (∆ = 0) and flat (constant) pulse of amplitude Ω 0 and duration T, the well-known π−pulse condition TΩ 0 = π in (4) gives Ω 0 /ω 0 = 1/(2N ). Thus, the constraint of a few oscillations in the field, typically N < 3, implies that Ω 0 , while still smaller than ω 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.

Few-Cycle-Pulse Adiabatic Floquet Theory
We consider the model (1) written with the time-dependent phase θ + φ(t), where the initial phase θ ≡ −φ 0,i serves for the Floquet representation, and the explicit dependence of the Rabi frequency Ω(t): The quasienergy operator (A19) features the two parameters Ω ≡ Ω(t) and ω ≡ ω(t) that can be designed independently as functions of time. They will both normalized with ω 0 . The resonant transformation (6) is represented in the Floquet picture by the transformation (which dresses the upper state with minus one photon): which leads to Note that we have omitted as before a diagonal term equal to −ω/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 Ω(t) = Ω 0 Λ(t) with Ω 0 the peak value and the shape 0 ≤ Λ(t) ≤ 1, with Λ(±∞) → 0, Λ(0) = 1, and of characteristic width τ, and a chirped frequency ω(t), of characteristic width around the resonance ∆ 0 , requires in general Ω 0 ∼ ∆ 0 1/τ. The chirped frequency induces for the full coupling Ω(t) cos φ(t) in Equation (1) the phase ω(s)ds, which can be rewritten as a function of the detuning; see Equation (2): We consider smooth Gaussian pulses Λ(t) = e −(t/τ) 2 , which allows in principle a transfer exponentially accurate as a function of the pulse duration [38] (see Appendix A.2). Defining τΩ 0 ∼ TΩ 0 = 2κπ (with a larger κ ensuring a better adiabatic passage) leads to Ω 0 /ω 0 = κ/N . This shows that adiabatic passage, for a given κ, breaks more strongly the rotating wave approximation for a smaller number N of cycles, since Ω 0 gets closer to ω. This requires to take into account the full quasienergy operator.
The quasienergies λ ±,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 N . The complete spectral information is contained in one Floquet zone composed of two surfaces λ ±,0 within a given band of energy of widthhω. The other surfaces can be constructed using the periodicity of the spectrum: λ ±,k ≡ λ ±,0 + khω, for any (positive or negative) integer k. In Figure 1 we display a few surfaces. Figure 2 shows cross-sections, one for a constant ω (left frame) and for constant Ω (right frame). Around the one-photon resonance ω ≈ ω 0 (for Ω ω 0 ), one can recognize the surfaces of Figure 1 of Ref. [17]: For Ω = 0, the horizontal line of energy −0.5hω 0 corresponds to the lower state of the system |−, 0 , with the notation |±, k for the bare state |± dressed with k photon, i.e., the state |± ⊗ |k with |k ≡ e ikθ . The crossing line (at ω = ω 0 ) corresponds to the upper state dressed with minus one photon |+, −1 . One can distinguish in the plane Ω = 0 resonances appearing as crossings at ω = ω 0 /(2k + 1) with k = 0, 1, 2 · · · For Ω = 0, they become avoided crossings. For ω = ω 0 /(2k), we have exact crossings for any Ω, 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 Ω), and the valleys to avoided crossings (i.e., to resonances). One can observe that for increasing Ω, the position of the resonances are shifted in the direction of larger ω. This can be interpreted as a Stark shift of the states. This implies that moving along a straight line with ω ≈ ω 0 for growing Ω 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.
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 onephoton resonance correspond thus to the constant distances |λ + − λ − | 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 Ω 0 = ∆(−∞), is given by (2). The initial phase φ 0,i , which does not modify the dynamics in the adiabatic limit, will be chosen to satisfy the zero time-integrated area (5). Considering N oscillations in the pulse, i.e., Tω 0 /2π = N , and Gaussian pulses with the estimate T ∼ 4πτ/3 for the pulse duration gives τω 0 ≈ 3N /2, i.e., This shows that smaller values N require larger ∆ 0 /ω 0 to keep the same τ∆ 0 . We have obtained the most extreme situation with a low number of cycles N ≈ 2.5 allowing ultrahigh-fidelity adiabatic passage with the following parameters. We consider the trajectory shown in Figure 3 around the one-photon resonance ∆ 0 /ω 0 = 0.625, which slightly avoids the three-photon-resonance curve, and the smallest quantity τ∆ 0 = 2.35 satisfying adiabatic passage (see the discussion above). We have determined numerically the value φ 0,i ≈ −0.3032π 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 ∆(−∞) = ∆ 0 ; we have checked that we have the same final result if we start with ∆(−∞) = −∆ 0 (but with a different initial phase φ 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 τ∆ 0 from (13).

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 π-pulse by expanding the Rabi frequency as a Fourier series: with the envelope of each pulsed mode defined as which satisfies 0 ≤ Λ n (t) ≤ 1. This particular Fourier expansion (14) implies that Ω(t) ≥ 0 and , exact resonance, ∆ = 0, and for any initial phase), imposing Ω N = 0 guarantees that the zero area condition (5) is satisfied.
If one considers a single field Ω(t) ≡ Ω 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 Ω 1 (t): The diagonal term can be neglected with respect to the coupling when (considering the condition for the peak Rabi frequencies) Ω 1 In the π-pulse regime, the pulse area of the pulse envelope is then Ω 1 T/2 = π, and the above condition becomes to neglect the diagonal term. This consists precisely to one of the RWA conditions, i.e., a large number of resonant cycles.

Few-Cycle Resonant Rabi Oscillations
In the resonant situation, ∆ = 0, we consider the number of cycles N = 2 and Ω 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 φ(t) = ω 0 (t − t i ) ± π/2 (i.e., φ 0,i = ∓π/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 + ≈ 0.99 for a peak amplitude Ω 1 ≈ 2π/T = ω 0 /N close to the standard π-pulse condition (with a relative deviation of order 2 × 10 −3 ). It also shows an ultrahigh-fidelity population transfer but for a high and quite unrealistic pulse area TΩ 1 /2 ≈ 6.38π. Figure 6 shows the weak dependance of the population transfer over the initial phase φ 0,i near the π-pulse regime. The maximum is obtained for the initial phase φ 0,i = ±π/2. We have determined that the population transfer depends on φ 0,i much more strongly for higher pulse area regimes (not shown).  (1). The transfer reaches a local maximum (P + ≈ 0.99) close to the standard π-pulse condition but does not reach ultrahigh fidelity (except for an undesirable large pulse area TΩ 1 /2 ≈ 6.38π).  Figure 7 shows the single-cycle Rabi oscillations, N = 1, for the phase of the field φ(t) = ω 0 (t − t i ) ± π/2, which satisfies the zero time-integrated area condition (5). Surprisingly, the maximum transfer (which occurs close to the standard π-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 δ distribution type interaction).  Figure 5 but for a single cycle N = 1. The transfer reaches a good fidelity P + ≈ 0.996 close to the standard π-pulse condition: However, in both situations, N = 1, 2, the transfer is not perfect, and cannot reach ultrahigh fidelity near the π-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 Ω n , 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.

2 N -Resonance Strategy
We consider here the situation N ≥ 2 and φ(t) = ω 0 t, i.e., an exactly resonant problem ∆ = 0, with φ 0,i = −ω 0 t i . We add a single 2N -mode to the main mode: Ω(t) = Ω 1 (t) + Ω 2N (t). The zero area condition (5) is automatically satisfied since Ω N = 0. In this case, the fast oscillating n = 2N -mode term features a resonance with the counter-rotating coupling: If we neglect the oscillating terms [second and third lines of (18)], the area of the resulting coupling is which should be π at the lowest order: According to (A52) and neglecting the faster oscillating term, one concludes that the dominant correction is given by the Ω 1 (t) counter-rotating term [first term of the second line of (18)]. Considering the peak terms, one can evaluate the (diagonal) correction: Its difference should be much smaller than the coupling term, giving the condition of the ratio We implemented a systematic search over the parameter Ω 1 for N = 2 and found an ultrahigh-fidelity transfer P + ≈ 0.9999 for Ω 1 ≈ 0.276ω 0 [leading to Ω 4 ≈ 0.448ω 0 according to (20)], corresponding to an equal distribution of both contributions in (20): . This gives in this situation for the above condition which significantly alleviates condition (17) (by a factor of 3).

Two-Modes-Resonance Strategy
In this strategy, still with N ≥ 2, we consider alternatively two additional nonresonant fields Ω 2N −1 and Ω 2N +1 (and φ(t) = ω 0 t, i.e., ∆(t) = 0) with the simplifying condition Ω 2N −1 = Ω 2N +1 ≡ Ω 0 , giving the total field This field satisfies the zero-area condition (5) as before since Ω N = 0. The Hamiltonian shows that the two fields generate a resonance: The area of the resulting coupling, which should be π, is Neglecting the fastest oscillating term, we obtain for the effective Hamiltonian i.e., for the peak values: The ratio of the (peak) diagonal correction with respect to the (peak) coupling is given by This ratio is minimum when Ω 1 = 0: R = Ω 0 /(3ω 0 ) = π/(3Tω 0 ). On the contrary, when Ω 0 = 0, the ratio is maximum. This solution minimizing the correction for the complete transfer is thus for Ω 1 = 0 and TΩ 0 = π, i.e., The diagonal terms can be neglected, from the above condition, when This shows that the introduction of the two fields in the limit of weak Ω 1 allows one to alleviate again condition (17) by a factor of 3, similarly to the preceding strategy.

Chirped Few-Cycle Pulses: Stark-Shift Compensation Strategy
The effective Hamiltonian (16) with a single-mode field Ω(t) = Ω 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 Ω 1 (t)]: still with a π-pulse: Ω 1 T/2 = π, where we have neglected the second term which is much smaller than the first one. The phase (2) can be then written where the initial phase φ 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 ±π/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 +π/2 has been chosen) with, respectively, ω 0 the mean frequency, which is the frequency that has to be tuned in practice, and φ 0 (t) the modulated phase which is the part of the phase that has to be shaped: Equation (35) shows a small shift of the mean frequency compared to the resonance ω 0 and a small phase modulation amplitude (of approximately 1/8N ). For instance, we obtain a mean frequency ω 0 = 1.0234ω 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 + ≈ 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 i.e., π/2 added to an odd function featuring modulating modes, with ω 0 the mean frequency, in order to satisfy condition (5) with an even single-mode field Ω(t) = Ω 1 (t), for any integer number of cycles N . This will be used for a more systematic optimization in the next section. Figure 8. High-fidelity population transfer for N = 2 with a single π-pulse mode Ω(t) = Ω 1 (t) and the phase (34) featuring a chirped frequency. From top to bottom frames: population dynamics, modulated phase φ 0 (t), detuning ∆(t), and full coupling Ω 1 (t) cos φ(t).

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).

2 N -Resonance Strategy
In order to target and validate this strategy, we impose the resonance, φ(t) = ω 0 t, and we implement an optimization using an expansion of Ω(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 Ω 1 ≈ 0.277ω 0 and Ω 4 ≈ 0.452ω 0 , see Figure 9. These parameters are close to the ones determined in Section 4.3. The area of the full envelope Ω(t) gives 1.5π. The area Ω 1 (t) + 1 4 Ω 2N is, as predicted, close to π: T(Ω 1 + 1 2 Ω 4 ) ≈ 1.006π. Adding more modes, as proposed in the next subsection, will allow one to explore a different strategy to achieve the same complete transfer result.

Two-Modes-Resonance Strategy
We now consider the resonant case φ(t) = ω 0 t with an expansion of Ω over the first 5 = 2N + 1 modes still for N = 2. Adding the fifth mode opens a path to the two-modesresonance 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 Ω 3 and Ω 5 are close to each other (see Figure 10), with a value consistent with Equation (30).  Figure 9 but for a five modes decomposition with the obtained non-zero modes Ω 1 = 0.043ω 0 , Ω 3 = 0.2314ω 0 , Ω 5 = 0.227ω 0 (two-modes-resonance strategy).

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 Ω to the first Fourier mode, Ω(t) = Ω 1 (t), and to decompose the phase φ(t) according to Equation (36) with a single mode: We obtain in both cases a complete transfer with smooth and simple phase and detuning. For N = 2, the optimized parameters are Ω 1 = 0.504ω 0 , ω 0 = 1.038ω 0 , and φ 0,1 = −9.74 × 10 −3 , giving a pulse area 1.0087π. The phase can be compared to (34) and (35): It features a larger shift of the mean frequency (though small compared to ω 0 ) and a smaller modulation amplitude (given by the absolute value of φ 0,1 , which is small compared to the angle π), of frequency ω 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π. One can notice a stronger (and negative, i.e.,φ < ω 0 ) shift of the mean frequency (still much smaller than ω 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 π). The modulation frequency of the phase is here ω 0 .

Conclusions
In this paper, we have explored quantum control with a few-cycle pulse implementing and generalizing two standard strategies, adiabatic Floquet theory and π-pulse. We have shown that adiabatic Floquet is obstructed by the presence of avoided crossing in the quasienergy spectrum induced by the counter-rotating term for the number of cycles 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 π-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 2N -resonance and two-modes-resonance strategies. Complete transfer can be achieved in these cases for a number of cycles as low as 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 multilevel 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.

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(θ + φ(t) + 2π) = H(θ + φ(t)) features a 2π periodic phase, φ(t) = ω(t − t i ) with θ the initial phase at time t i , ω the angular frequency. A typical example is for a quantum system of Hamiltonian H 0 (x) driven by a periodic field of amplitude E via the dipole moment µ(x). The Schrödinger equation is defined on a Hilbert space H, which can be of infinite dimension (e.g., the space of squareintegrable functions H = L 2 (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 H = C N ). The initial phase θ appears as a parameter. One can think of Equation (A2) as a family of equations parameterized by the angle θ. We define the solution ψ(t) of a time independent problem: with the time-independent Floquet Hamiltonian (or quasienergy operator) defined in an enlarged Hilbert space where L := L 2 (S 1 , dθ/2π) denotes the space of square integrable functions on the circle S 1 of length 2π, with the scalar product The solution ϕ(t) is related to the time independent problem as with the initial condition in K constructed as the initial condition ϕ(t i ) of the original problem lifted in the enlarged space as: with 1l L ≡ e i(k=0)θ . In the expression ϕ(x; t; θ), on the left-hand side of (A7), θ is a parameter corresponding to the initial phase of the Hamiltonian, while, in ψ(x, θ; t i ), on the righthand side, θ is a dynamical variable on the same footing as x. The procedure dictated by Equation (A7) to obtain the dynamics in the original space H is as follows: (i) first, by the use of Equation (A8), lift the initial condition in the enlarged space K keeping θ undetermined; (ii) next determine ψ(x, θ; t) in K solution of Equation (A3); (iii) finally replace the angle θ by θ + ω(t − t i ) in ψ and fix the angle θ to the initial phase of the Hamiltonian.
Since K is time-independent, the convenient way to compute the solution is to use an eigenfunction expansion in terms of the eigenelements of K (assuming a discrete spectrum) with the eigenvalues named quasienergies, with whereψ ν (x) := S 1 dθ/2π ψ ν (x, θ) is the average of ψ ν (x, θ) 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 ν being associated to two indices, ν ≡ n, k: where the index n refers to the molecule's Hilbert space H (i.e., n = 1, · · · , N if H = 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): The interpretation of this expansion is as follows: the initial condition is expanded into the set of representatives of each family of eigenvectors ψ n,0 (x, θ) through the coefficient c n (θ) (which depends parametrically on θ); during the evolution, the expansion acquires dynamical phases through the eigenvalues λ n,0 and features a dynamical evolution of θ 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: The 2π-periodic phase θ + φ(t) features an instantaneous (or effective) frequency ω(s)ds, allows one to recover the phase dependence of the previous periodic system when the time of the frequency is frozen.
The state ϕ(x; t; θ) solution of the Schrödinger equation is connected to ψ(x, θ; t) the solution of with with the initial condition 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 ∂ ∂θ 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 |ψ n,k , 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 |ψ 0 n,k ≡ |n ⊗ |k with |n a bare state of the quantum system and |k ≡ e ikθ a (relative) state of k photons before the field is on, i.e., before interaction between the matter and the field, to the eigenstate |ψ 0 n ,k ≡ |n ⊗ |k after the field is off, i.e., after the interaction, corresponds to the transfer from |n to |n 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 ≡ s with 1, one can apply the above formulation with the use the eigenelements of K, which depends now on the parameters: on which we derive an adiabatic approximation 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 ψ r(t i ),ω(t i ) n,k(n) (x, θ) at the initial values of the parameters r(t i ), ω(t i ) through the coefficientc n (θ) (which depends parametrically on θ); during the evolution, the expansion acquires dynamical phases through the integration of the eigenvalues λ r(t),ω(t) n,0 via the time-dependence of the parameters r(t), ω(t) and features a dynamical evolution of θ 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 widthhω 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( ). However it can be made in principle exponentially small, O(e −C/ ) 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][40][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 Ω(t) under the constraints given by the RWA π-pulse conditions. The second loop explores a larger space of solutions by removing the RWA π-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 Ω N = 0 in the Fourier expansion of Ω(t). In the third loop, this condition (5) is satisfied by the specific form of the phase (36). Figure A1. Flowchart of the optimization algorithm. The modified Crank-Nicolson scheme can be found in Ref. [44].

Appendix C. Perturbation Theory Formulated with KAM Techniques
We decompose the Hamiltonian with an unperturbed Hamiltonian H 0 and a perturbation εV 1 : We construct a unitary transformation (so-called KAM or contact transformation) T 1 = e εW 1 , with W † 1 = −W 1 such that e −εW 1 He εW 1 where D 1 is a diagonal operator, i.e., satisfying [H 0 , D 1 ] = 0. Thus the perturbation is reduced from order ε to order ε 2 . In order to determine D 1 and W 1 , one expands Equation (A26) in powers of ε and require that the terms of order ε cancel out: e −εW 1 He εW 1 = 1l − εW 1 + 1 2 ε 2 W 2 1 + · · · (H 0 + εV 1 ) 1l + εW 1 + The remaining perturbation of order 2 can be written as We can apply a second contact transformation e 2 W 2 (with respect toK 0 ), W 2 = −i θ dθ (V 2 − V 2 ), which averages this rest (A44) with respect to θ and leads to correction of order 3 . We thus obtain the effective high frequency Hamiltonian H HF (independent of θ) of second order (A45) As an example, one considers an Hamiltonian of the form: i.e., where in the high frequency limit, rewritten with the dimensionless coefficients Ω −2 /2ω → 0, Ω 2 /2ω → 0, Ω 0 /2ω → 0. We obtain V 1 = 0, (Ω 2 −2 − Ω 2 2 ) i.e., at the second order The high-frequency effective Hamiltonian reads then The diagonal term can be neglected when If we consider a time-dependent Ω(t), the transformation T 1 = e W 1 introduces in the Schrödinger equation an additional term of the form: which should be treated perturbatively like V 1 . It introduces an additional term in the diagonal of (A52):