Designing adiabatic time evolution from high-frequency bichromatic sources

We investigate the quantum dynamics of a two-level system driven by a bichromatic ﬁeld, using a nonpertur-bative analysis. We make special emphasis in the case of two large frequencies, where the Magnus expansion can fail

Introduction: Perturbing a system out-of-equilibrium is at the heart of physics, as it allows to extract information about its properties by just measuring the response to the perturbation.Besides small perturbations, one can also produce non-linear effects of high complexity, and steady states with novel properties such as Floquet topological insulators, skyrmions or time-crystals [1][2][3][4][5][6].Applying periodic perturbations has shown to be a versatile tool to manipulate physical systems.For instance, they allow to control spin qubits in quantum dots [7][8][9][10][11], or to induce new electronic, dynamical and topological properties [12][13][14][15][16].These works typically consider monochromatic driving, although bichromatic fields have been used in a few occasions [17][18][19][20], showing that their potential has not been fully explored.
The periodically driven two-level systems is one of the fundamental models in quantum mechanics.Its physical realization has been successfully implemented in quantum dots [21][22][23][24], complex molecules [25,26], superconducting devices [27], and many other systems [28].Its universality relies on the fact that many quantum mechanical systems, when truncated to their low lying states by lowering the temperature, can reduce to the dynamics between the ground state and the first excited state.
The Hamiltonian describing the unperturbed two-level system usually displays a splitting ∆ z .Then, one chooses this direction as the quantization axis, and performs transitions between the ground state and the excited state to probe the system, which is the guiding principle in techniques such as nuclear and electron spin resonance.The Hamiltonian describing the model can be written as and where V i , ω i and φ i correspond to the different amplitudes, frequencies and phases of the external source, respectively.The dynamics of this simple Hamiltonian can be complicated, even in the monochromatic case, as the three energy scales involved (∆ z , V 1 and ω 1 ) can lead to very different behavior.The standard perturbative analysis in V 1 ≪ ∆ z , ω 1 explains the linear response regime.This is commonly used to probe the system and obtain information about its physical properties [29].On the other hand, one can consider the (high frequency) strongly driven regime V 1 > ω 1 and ω 1 ≫ ∆ z , which produces the spectral changes typically found in Magnus expansions [30], and can be used to dynamically tune the properties of the system [2,5,12,31,32].Finally, the resonant behavior corresponds to the case ω 1 ≃ ∆ z , which produces a transfer of spectral weight from the ground state to the excited state.This is used for state preparation in many experiments [33] or to induce single qubit gates [24].On top of that, environmental degrees of freedom in experiments can also couple to the external field, producing an undesired large signal if they are resonant.
In this work we study a two-level system driven by a bichromatic field and discuss the advantages over the monochromatic case.We analyze the limitations of high frequency expansions, when more than one frequency is present, and introduce a description which captures the full dynamics, including the fine details of the micromotion.This is in contrast with expansions where just the stroboscopic evolution is considered, missing the dynamics within a period, which can be crucial for a topological analysis [34].
We find that when the two frequencies are large, but close to each other, the field amplitudes can be controlled to produce an effective adiabatic evolution.We demonstrate that the adiabatic behavior is robust to noise, and highly tunable.Furthermore, we show that non-adiabatic corrections to the resulting effective Hamiltonian can be beneficial, and used to engineer adiabatic rotations on the Bloch sphere.We also demonstrate that the longtime dynamics is controlled by multi-photon resonances.In addition, our analysis can be extended to the analysis of a multi-level system under multi-chromatic driving in a straightforward manner.
Motivation: As an illustrative example to understand the breakdown of high frequency expansions in a multifrequency case, let us first consider a rather simple trigonometric property of the function g (t) = cos (ωt).It can generally be written as: being ω = ω 1 − ω 2 and their sum arbitrary.If ω i ≫ ω, one can interpret the slowly evolving oscillatory function g (t) as coming from the difference of two large frequencies, whose difference is very small.What in this case makes possible to exactly map the two high frequencies to an adiabatic evolution, is the specific relation between their Fourier components, where only the crossed terms g ±1,∓1 of the two-dimensional Fourier expansion n1θ1+n2θ2) g θ contribute, being g θ the function g (t) re-parametrized according to θ i = ω i t [35].This illustrates a specific case where a system driven by two initially large frequencies, will not give a converging result using a high frequency expansion.The reason is that a system with such a driving term (with ω smaller than all the characteristic energies of the model), obviously requires an adiabatic analysis [36] due to its slow time evolution.Furthermore, it provides some intuition about the requirements to engineer a specific dynamical behavior in a quantum system, by studying its Fourier decomposition.
Bichromatic two-level system: Let us now move to the problem at hand.We choose a simple harmonic protocol for the external drive f i (t, ω i , φ i ) = cos (ω i t + φ i ), although other choices are possible (i labels each different component of the drive).One could also choose each term in Eq.2 coupled to a different, non-commuting, degree of freedom (e.g., to a σ y component).However, this is not necessary for the present analysis and will be discussed below.We consider the specific case of bichromatic drive, where the sum in Eq.2 is restricted to two terms only.Each term is characterized by a different amplitude, frequency and phase.If we parametrize ω i t → θ i to obtain a two-dimensional Fourier representation [35], the only non-vanishing Fourier components are the terms H 0,±1 and H ±1,0 , which characterize the time-evolution operator for small field amplitudes.This indicates that for weak amplitudes, the behavior is dominated by oscillations with frequencies ω 1,2 only.
If instead we perform a non-perturbative analysis of H (t), by applying the transformation U (t) = exp −i V (t) dt we find the following transformed multi-chromatic time-dependent Hamiltonian: where we have defined being this result still valid for an arbitrary number of drive components.The advantage of Eq.4 relies on the fact that, from a Jacobi-Anger expansion [37], one finds non-perturbative expressions in V i .Notice that this type of Hamiltonian is directly obtained in electronic systems via the Peierls substitution [12], indicating that our results will be valid for seemingly different systems, connected by unitary transformations.However in the latter case, the transformation to the interaction picture is not required.
In the monochromatic case, the phase φ 1 in Eq.4 does not affect the spectrum, because it is a gauge degree of freedom that sets the origin of the time-evolution, but it affects the dynamics.However, in the multi-chromatic case the phase differences are relevant and both, the spectrum and the dynamics are affected.This provides an extra degree of freedom in Floquet engineering, absent in the monochromatic case, although we will not make use of it in the present work.
The two-dimensional Fourier decomposition of Eq.4 (we parametrize θ i = ω i t): leads to the following expression for the Fourier components of the Hamiltonian: where α i = V i /ω i and we have used the Jacobi-Anger expansion e iz sin(θ) = ∞ n=−∞ J n (z) e inθ .Eq.6 can be further simplified if the two frequencies are commensurate [14,37], however we will consider the general case.
Notice that Eq.6 contains an infinite number of Fourier components, and in contrast with the unrealistic, but pedagogical case of g (t) above, the adiabatic behavior will only happen if we can enhance the crossed Fourier components H±1,∓1 .Fortunately, due to the non-perturbative expressions in α i this is now possible, and we choose values of α i that maximize J 1 (α i ), while requiring ω − ≪ ∆ z and ω i ≫ ∆ z (we have defined ω ± = ω 1 ± ω 2 and fixed φ i = 0 for simplicity): sin (ω i t) σ y + . . .Noticing that the fast oscillating term in the second line averages to zero in this regime, and that the third line contributes with a small amplitude correction to the time evolution, we obtain the next leading Hamiltonian: Actually as ω i ≫ ∆ z , all terms with a fast frequency dependence, of the order of ω i or larger, tend to zero as ω i → ∞.To demonstrate this we have numerically calculated the time evolution operator Ũ (t) in the interaction picture.Fig. 2 shows the real part of the component Ũ 1,1 (t) over time, which is related with the occupation probability of the excited state.Different colors indicate the time-evolution averaged over a hundred realizations of white noise, following a normal distribution with standard deviation σ.For weak noise, the plot displays adiabatic behavior (notice that T − ∼ 125 while T 1,2 ∼ 0.6 is the period of the driving fields) with small and fast amplitude oscillations around the mean value, coming from the high frequency corrections (see Fig. 2, inset).This is one of the main results of this work: One can combine two high frequency drives to effectively evolve the system adiabatically.Actually, Fig. 2 contains a combination of two different periods, which can be externally controlled: harmonic oscillations due to the static part of the effective Hamiltonian (first term in Eq.8), and a slow frequency modulation due to the adiabatic term (second term in Eq.8).Reducing the difference ω − = ω 1 − ω 2 makes the static part dominate at short times, while its increase makes the non-linear, adiabatic term to take over (these two cases are explicitly shown in the Appendix A, Fig. 4 with an extended discussion about the role of each term in the effective time-dependent Hamiltonian).
In addition, the robustness of our prediction is illustrated with the changes of the time evolution as the noise increases.Fig. 2 shows that the original behavior persists for weak values of noise, until σ becomes of the order of the dominant energy scale, where the oscillations are strongly damped.This indicates that experiments would have time to perform a few adiabatic cycles before the effect of noise takes over.
Corrections to the effective adiabatic Hamiltonian: We have shown that it is possible to drive a system with two different frequencies and produce effective adiabatic behavior.The dominant part of the Hamiltonian, shown in Eq.8, characterizes the slow evolution, but corrections due to oscillatory terms with higher frequency are also present and they can produce transitions between the adiabatic eigenstates (specially relevant are the terms proportional to σ y in Eq.7, because they do not commute with the leading order Hamiltonian).To understand their effect we now introduce a general formalism to study Hamiltonians with two different frequencies, where each drive couples to a different degree of freedom.The general form of the Hamiltonian can be compactly written as: the functions ∆ (t) and Ṽ (t) are, for the moment, general harmonic functions, and ǫ is introduced to organize the perturbative series (we take ǫ → 1 at the end of the calculations).To study the non-perturbative dynamics we will calculate the time-evolution operator using multiplescales analysis.This method deals with the fastest timescales first, and then each correction characterizes new processes taking over at longer times.Importantly, this method includes a renormalization procedure for the secular terms [38][39][40], making the solutions valid at high/low frequency, as well as near a resonance.This approach only neglects processes that take over at longer timescales than the order of the expansion in ǫ.
In the first step we parameterize the time evolution operator U (t) → U ( τ ) in terms of a set of time-scales τ n = ǫ n t, and expand in powers of ǫ the equation of motion for the time-evolution operator, with U ( τ ) = n ǫ n U n ( τ ).The equations to zeroth and first order in ǫ result in: where we have omitted the τ dependence in U n ( τ ).The lowest order solution can be easily obtained from Eq.11 by direct matrix exponentiation: where the matrix u 0 (τ 1 ) comes from the boundary condition and will be determined later on, during the renormalization procedure.Eq.13 can now be inserted in Eq.12 and solved by choosing τ0)dτ0 the solution to the homogeneous equation.The solution results in: where we have defined Eqs.13 and 14 are the formal solutions for the time evolution operator, which now need to be particularized for the case of interest and renormalized, if needed.For our present purpose we fix the specific form of the periodic functions to ∆ (t) = ∆z + µ cos (ωt) and Ṽ (t) = β cos (Ωt).In this case ∆z corresponds to the static part (first term in Eq.8), µ to the dominant time-dependent term (second term in Eq.8), and β to the dominant correction, non-commuting with Eq.8 (the value of the two frequencies ω and Ω is arbitrary for the moment).The unperturbed solution is obtained from Eq.13 and displays the non-linear phase evolution typically obtained in timedependent systems: where u 0 (τ 1 ) still needs to be determined.Similarly, the first order solution in ǫ is obtained from Eq.14.We do not write here the full form of the solution, because of its length and because it is enough to show that is proportional to (details of the calculation in the AppendixB): This correction to U 0 ( τ ) diverges if the denominator in Eq.17 vanishes.This is a common feature of timedependent perturbation theory, indicating the breakdown of the solution, but these resonances can be renormalized in multiple-scales analysis, and produce nonperturbative corrections to U 0 ( τ ).Strictly speaking, the resonance condition can only be fulfilled for commensurate frequencies (which differentiates this case with the one of incommensurate frequencies), and for very specific values of the parameters.However, if the denominator in Eq.17 becomes smaller than β, the perturbative series still diverges and should be renormalized as well, making the difference between incommensurate and commensurate frequencies (with very long total period) merely a mathematical curiosity, for the physically relevant timescales of this setup.Therefore, one can relax the strict relation between the parameters for a resonance, to just the approximate one: nω ± ∆z 2 − Ω 2 β.To renormalize the resonant terms, one needs to separate resonant and off-resonant contributions.The amplitude corrections produced by the off-resonant terms in Eq.17 are of order β, and can be neglected if we focus on U 0 (t) only.However, resonant corrections contribute to leading order and need to be included.We assume that the system is in the regime ω ≪ ∆z ≪ Ω, where ω corresponds to an adiabatic drive and Ω to a high frequency one.This situation is analogous to the one obtained in Eq.7 for the effective adiabatic Hamiltonian.In this situation, several resonances can contribute (i.e., several values of n fulfill nω ± ∆z 2 − Ω 2 β), while in a different regime the analysis would be simpler, because the resonances do not need to be included.Then, in the spirit of multiple-scales analysis, we require that the secular terms produced by the resonances are cancelled by ∂ τ1 u 0 (τ 1 ) in Eq.14.This requirement leads to the following flow equation for u 0 (τ 1 ): where n 0 corresponds to the set of resonances {±n 0 } which fulfill the approximate resonance condition above.This equation allows to determine u 0 (τ 1 ), which encodes the non-perturbative correction to U 0 (t).The lowest order renormalized solution becomes: Jn 0 ( µ ω )σy (19) Notice that the smaller ω is, the larger is the set of resonances ±n 0 that needs to be included, increasing the contribution from the non-perturbative correction.Furthermore, this correction strongly depends on the ratio µ/ω.This indicates that if the system is far from a resonance, or µ is not in the region where J n0 (µ/ω) has a relevant weight, the behavior is similar to that of the unperturbed solution (at least to time-scales of the order of β −2 ).
The multiple-scales analysis can be continued to higher orders in a very systematic way, however the results presented here are quite accurate for the range of parameters under consideration.In Fig. 3 we show a comparison between the numerical and the analytical approximation Comparison between the exact and the approximate dynamics to first order in ǫ for the real part of U 1,1 (t), in the presence of resonances (Eq.19).The short time dynamics is well captured by the unperturbed solution, but the slow oscillations at longer times are obtained from the renormalization of resonances.Parameters: ω/ ∆z = 0.1, Ω/ ∆z = 2, β/ ∆z = 0.2 and µ/ ∆z = 2.In this case the dominant resonance is obtained for n = ±10, and its contribution perfectly captures the slow modulation of the oscillations.
for the time-evolution operator.One can identify three different time scales in this plot: 1.The shortest time-scale is given by slow harmonic oscillations coming from the static part of the unperturbed solution ∆z (as those oscillating between ±1 in Fig. 2).
2. The next time-scale corresponds to the non-linear phase evolution.Proportional to µ, introduces the anharmonic oscillations happening at intermediate times and define the adiabatic period T − .
3. The longest time-scale is produced by the nonperturbative correction produced by H1 (t).It produces the long-time modulation observed in Fig. 3 (in this case T − ∼ 60 and the long-time modulation has period τ long ∼ 600, i.e., one order of magnitude larger).
We also show the dynamics for the off-diagonal component of U (t) in the Appendix B, Fig. 7, to confirm that the long-time behavior is controlled by the nonperturbative correction u 1 (τ 1 ), which produces the rotation proportional to σ y in Eq.19.
Conclusions: We have demonstrated that bichromatic driving provides new possibilities to externally control quantum systems.An interesting one is that two high frequencies ω 1,2 can produce effective adiabatic evolution, with frequency controlled by the difference ω 1 − ω 2 .This effect requires to strongly drive the system beyond the perturbative regime, and provides an example of the breakdown of high frequency expansions.
The effect can be used in experiments where low frequencies are out of reach due to equipment restrictions, or if the slow, monochromatic drive resonantly couples to environmental (or undesired) degrees of freedom [41].This is because the effective Hamiltonian (Eq.6) is a function of the coupling strength between the drive and each degree of freedom (in this case controlled by Bessel functions).Then, as the coupling to the environmental modes is different, their Fourier components will be tuned at a different rate with the field amplitude and generally suppressed, while the one of interest is being enhanced.This would allow to reach the desired monochromatic behavior for the degree of freedom of interest, while reducing the undesired signal from the environment [42].
As the effective adiabatic Hamiltonian (Eq.8) generally contains high frequency corrections, we have also studied their effect.This is equivalent to a bichromatic system with slow and fast frequencies, coupled to different, non-commuting degrees of freedom.In this case, we have shown that the high frequency corrections of the effective Hamiltonian can be used to engineer controlled single-qubit rotations.At short time-scales the adiabatic part dominates, and one can switch between free and adiabatic evolution by adjusting the frequency difference ω − = ω 1 − ω 2 .At longer time-scales the high frequency corrections become relevant and produce adiabatic evolution between the ground and the excited state.This extra adiabatic evolution along a perpendicular direction is controlled by resonances involving the slow and the fast frequencies, and its period depends on the amplitude of both time-dependent terms (µ/ω and β).This provides a highly tunable mechanism to implement single-qubit gates using two off-resonant fields only.
Further applications of our results are the possibility to externally control quantum pumping [43,44] in higher dimensional systems[45], or to describe Floquet topological phases at low frequencies.This is because our approach (multiple-scales analysis) allows a complete characterization of the evolution operator, which is required for the topological analysis [46].In qubits, it would also be interesting to study the competition between geometric and dynamical phases, as the difference between frequencies is tuned.This could be implemented in several experimental setups such as quantum dots, N-V centers, single-ion magnets or superconducting junctions.
This work was supported by the Spanish Ministry of Economy and Competitiveness through Grant No. MAT2017-86717-P and we acknowledge support from CSIC Research Platform PTI-001.Á. G.-L. acknowledges the Juan de la Cierva program.
[1] Lindner NH, Refael G, Galitski V. Floquet topological insulator in semiconductor quantum wells.Na-  where ∆z corresponds to the static part of the effective Hamiltonian (first term in Eq.A13), µ to the dominant periodic modulation with frequency ω (second term in Eq.A13), and β to the transverse oscillating correction with frequency Ω.This way, the unperturbed solution corresponds to: and the first order solution is obtained from the rotated Hamiltonian: To calculate the correction from Eq.B7, one can write the exponentials in terms of Bessel functions: and perform the following integrals: These integrals will produce secular terms if ∆z ± nω ∓ Ω = 0, but they can be cancelled by the term ∂ τ1 u 0 (τ 1 ) in Eq.B7, defining the flow equation.Furthermore, even for case where the resonance condition is approximately fulfilled ∆z ± nω ∓ Ω β only, the perturbative solution would not converge, and the renormalization can be applied.Fig. 5 graphically shows, for a specific case, the harmonics n that require renormalization (n = 10 in this case).Larger ω increases the distance between different harmonics, making more difficult to find a resonance.In this case, beyond ω 5 resonances are no longer possible, signaling the transition to the high frequency regime [47].Once the subset of n values which produces secular terms is identified, the cancellation with ∂ τ1 u 0 (τ 1 ) in Eq.B7 produces the following flow equation: where n 0 is the set of pairs of integers fulfilling the condition n 0 ≃ ± Ω± ∆z ω .Then, the lowest order solution is given by: U 0 (t) ≃ e − i  This non-perturbative correction indicates that the oscillations of each independent energy level are now modulated by a transition between the ground and the excited level with frequency β 4 n0 J n0 µ ω .Fig. 6 compares the exact dynamics with the one generated by U 0 (t) without/with renormalization (left/right).The addition of non-secular contributions produces the plot in the main text (Fig. 3 main text).However, the importance of the resonances is evident from the comparison between the left and right plot for U 0 (t) in Fig. 6.Finally, we plot in Fig. 7 the exact dynamics of U 1,2 (t), to confirm that the resonances control the rotation proportional to σ y and that their contribution is non-perturbative.This is confirmed by noticing that the off-diagonal part initially vanishes, but acquires values of the order of one for times of approximately half the period obtained from the resonances.

Figure 1 .
Figure 1.Schematic representation where a system driven by a high frequency bichromatic field behaves as adiabatically driven by a single frequency.

Figure 2 .
Figure 2. Exact dynamics during an adiabatic period T− = 2π/ω− for the real part of the upper diagonal element of the time-evolution operator Ũ (t).The presence of white noise following a normal distribution with a standard deviation σ has been considered.The colors indicate different noise strength σ, in units of the dominant energy scale ∆ = ∆zJ 2 1 (α).Parameters: ω1/∆z = 10, ω2/∆z = 10.05,α ≃ 1.8 (first maximum of J1 (α)) and φi = 0.The inset shows the short time dynamics, where the fast oscillations from the fast applied drive (with period T1,2 ≃ 2π/ω1,2) are more evident.All the other components of the evolution operator also evolve adiabatically.

Figure 4 .
Figure 4.In red(blue), the dynamics upper diagonal component of the time evolution operator for ω1/∆z = 10, ω2/∆z = 10.05 (10.01), αi ≃ 1.8 (this is the first maximum of J1 (αi)) and φi = 0.The inset shows the effect of the high frequency corrections, as small oscillations with the frequency ωi of the original drive.For ω− = 0.01 (blue) the static part of the effective Hamiltonian dominates at short time, displaying harmonic oscillations with non-linearities taking over at later times (not shown).For ω− = 0.05 (red) the non-linear part produces an earlier adiabatic frequency modulation.

Figure 5 . 2 − Ω 2
Figure 5. Curves nω ± ∆z 2 − Ω 2 for different values of n (red and blue dots) and a region of width β (green, see Eq.B10).Points in the green area fulfill the resonance condition, and require renormalization.Increasing ω from adiabatic to diabatic values produces a transition from several resonances to none.Parameters for the plot: Ω/ ∆z = 4, ω/ ∆z = 0.3 and β/ ∆z = 0.4.

Figure 6 .
Figure 6.Comparison between the exact dynamics (red solid) and the one generated by U0 (t) (blue dashed) without/with renormalization of the resonances (left/right).We have plotted the real part of the U 1,1 (t) component, but the agreement is valid for all the other components as well.The addition of small corrections of order β to the renormalized U0 (t) leads to an even better agreement, as shown in Fig.3of the main text.Parameters: ω/ ∆z = 0.1, Ω/ ∆z = 2, β/ ∆z = 0.2 and µ/ ∆z = 2.