Fast adiabatic evolution by oscillating initial Hamiltonians
Francesco Petiziol Benjamin Dive Florian Mintert Sandro Wimberger


APS/123-QED




 francesco.petiziol@studenti.unipr.it
Dipartimento di Scienze Matematiche, Fisiche e Informatiche, Università di Parma, Parco Area delle Scienze 7/a, 43124 Parma, Italy
INFN, Sezione di Milano Bicocca, Gruppo Collegato di Parma, Parco Area delle Scienze 7/a, 43124 Parma, Italy

Department of Physics, Imperial College, SW7 2AZ London, UK

Department of Physics, Imperial College, SW7 2AZ London, UK

Dipartimento di Scienze Matematiche, Fisiche e Informatiche, Università di Parma, Parco Area delle Scienze 7/a, 43124 Parma, Italy
INFN, Sezione di Milano Bicocca, Gruppo Collegato di Parma, Parco Area delle Scienze 7/a, 43124 Parma, Italy


We propose a method to produce fast transitionless dynamics for finite dimensional quantum systems without requiring additional Hamiltonian components not included in the initial control setup, remaining close to the true adiabatic path at all times. The strategy is based on the introduction of an effective counterdiabatic scheme: a correcting Hamiltonian is constructed which approximatively cancels nonadiabatic effects, inducing an evolution tracking the adiabatic states closely. This can be absorbed into the initial Hamiltonian by adding a fast oscillation in the control parameters. We show that a consistent speed-up can be achieved without requiring strong control Hamiltonians, using it both as a standalone shortcut-to-adiabaticity and as a weak correcting field. A number of examples are treated, dealing with quantum state transfer in avoided-crossing problems and entanglement creation.


Fast adiabatic evolution by oscillating initial Hamiltonians
    Sandro Wimberger
    December 30, 2023
============================================================






§ INTRODUCTION


Quantum adiabatic processes <cit.> are ubiquitous in quantum science and they represent an important resource for quantum control.
The adiabatic theorem states that when the system Hamiltonian H(t) is varied slowly enough in time from an initial configurations H(t_i) to a final configuration H(t_f), then a state initially in an eigenstate of H(t) will remain so during the whole evolution <cit.>. This makes such protocols intrinsically robust against experimental imperfections. On the other hand, the necessity of very long timescales for their implementation dramatically limits the number of operations which can be performed on the system within reasonable coherence times.

Recently, much work has been done towards the design of adiabatic-inspired control strategies which, on the one hand, inherit the robustness of the adiabatic dynamics, while, on the other, avoid the necessity of slow driving.
Among such so-called “shortcuts to adiabaticity” (STAs), a particularly promising method was introduced under the name of counterdiabatic (CD) <cit.> or transitionless <cit.> quantum driving. The basic idea which is put forward is that it is always possible to reverse-engineer a correcting Hamiltonian (t) such that the  total Hamiltonian H(t) + (t) keeps the system in the instantaneous eigenvectors of H(t) without requiring it to change slowly – that is, such that the adiabatic dynamics is an exact solution of the time-dependent Schrödinger equation.

When facing a quantum control problem, one must take into account that only a fairly restricted number of Hamiltonians can be realized and controlled in practice. The crucial drawback of the CD method is that, although a well-defined expression for computing (t) exists, the correcting field typically requires time-dependent control of complex interactions, and more generally of Hamiltonians which do not belong to the available control setup. As a result, the implementation of (t) is often tricky if even possible <cit.>. For this reason, different STAs have been developed which connect adiabatic states at different desired times, but completely deviate from the adiabatic states at intermediate times <cit.>, therefore giving up the benefits of true adiabaticity.

Here, we propose a method for achieving fast adiabatic driving remaining close to the adiabatic path, without needing new unrealizable terms in the Hamiltonian. This works by modulating the original Hamiltonian of the system in time such that it effectively replicates the dynamics induced by  without needing any additional control Hamiltonian on the quantum system.

In order to do so, we first resort to control-theoretic techniques to study how the matrix structure of the correcting field  is related to the initial set of control Hamiltonians which constitute H. This shows that  can always be emulated, to arbitrary precision, by introducing a suitable (fast) time dependence in the control parameters. Second, we identify a class of Hamiltonians, generalizing the set of real ones, for which the matrix structure of  can be discussed on general algebraic grounds and always involves components which are not directly controllable. They can, however, be simulated arbitrarily well using existing protocols.

Building on these results, we describe the construction of an effective counterdiabatic (E–CD) field (t), which is a time-dependent combination of the initially available control Hamiltonians. The control functions in  will be chosen to be oscillating fast with respect to the natural time dependence of H, and  will be enforced to simulate the dynamics  induced by .
As a result, the E–CD evolution tracks the adiabatic path, being arbitrarily close at a set of sampling time points while slightly deviating at intermediate ones. The general idea of producing the  dynamics approximatively by working with a set of available control terms was also pursued, yet with a different strategy, in Ref. <cit.>.

One must take into account that the adiabatic requirement of slow drivings can be recast as the need of strong fields. With this in mind, we characterize the efficiency of the E–CD scheme in terms of final fidelity and duration for different strengths of the correcting Hamiltonian  as compared with the uncorrected H. 
We show that the E–CD field, when acting as an auxiliary weak term in H, realizes a consistent speed up of the adiabatic evolution. 
The E–CD Hamiltonian also works as a stand-alone STA, giving even better results but at the price of completely losing true adiabaticity. The net result is the ability to generate entanglement or perform state transfer significantly faster which, when taking into account noise from an external environment, inexorably leads to higher fidelities.

After recalling the theory of CD driving in Sec. <ref>, the control-theoretic setting in presented in Sec. <ref>, together with the results on the general matrix structure of the correcting Hamiltonian . The effective CD method is introduced in Sec. <ref>, where the derivation of the E–CD field is fully worked through for a single spin system. The results are then exemplified via a number of applications in Sec. <ref>, the first being the Landau-Zener-Majorana model <cit.>. We afterwards show that a two-qubit entangled state can be prepared with 99.9% fidelity ten times faster with respect to purely adiabatic evolution, using an E–CD field with strength comparable to that of the original Hamiltonian. We further discuss the case of a three-level system dynamically undergoing a sequence of avoided crossings in the energy spectrum <cit.>. In this section, the E–CD method is also benchmarked against traditional finite-time adiabatic driving, and the main advantages and limitations are discussed in detail, before proceeding to the conclusion in Sec. <ref>.



§ COUNTERDIABATIC DRIVING
 
In this section, we recall the central elements of the theory of counterdiabatic fields <cit.> which underpins our method.

A unitary evolution U(t) is always a solution of the Schrödinger equation with Hamiltonian H_U(t)=i ħ∂_t U U^†. Equivalently, the same H_U(t) is a reverse-engineered Hamiltonian producing dynamics U(t) exactly. Our aim is to find this H_U(t) for the U(t) which performs perfect adiabatic transfer. By construction, the new Hamiltonian will give the desired dynamics in an arbitrarily short period of time.

Let the initial Hamiltonian of the physical system be 

    H(t)=∑_n=1^N E_n(t) |n(t)⟩⟨n(t)|,

having instantaneous eigenvalues E_n(t) and instantaneous eigenvectors |n(t)⟩. Let (t) = ∑_n e^-i φ_n(t)/ħ|n(t)⟩⟨n(t_0)|, with φ_n arbitrary phases, be a unitary matrix which diagonalizes H(t) at all times. That is,

    ^†(t) H(t) (t) = diag{E_1(t),…,E_N(t)}.

The corresponding reverse-engineered Hamiltonian = i ħ∂_t ^† reads

    (t) = ∑_n ∂_t φ_n(t) |n(t)⟩⟨n(t)| + (t),

where we have introduced the CD field <cit.> 

    (t) = i ħ∑_n |∂_t n(t)⟩⟨n(t)|.

The Hamiltonian  drives the instantaneous eigenvectors of H(t) exactly, with relative phase factors φ_n(t), e^-i φ_n(t)/ħ|n(t)⟩. In particular, one can choose φ_n(t)=0 for all n, which gives =. Therefore, implementing just  is sufficient if one is only interested in preserving the populations of the instantaneous eigenvectors. Another interesting choice is φ_n(t) = ∫^t E_n(t') dt', in which case φ_n are the adiabatic dynamic phase factors. This gives = H(t) +, so that  can be interpreted as a correcting field acting beside the initial Hamiltonian H(t).

Alternative expressions of  are useful in order to highlight specific properties. For example, from

    = i ħ∑_m n∑_n=1^N |m⟩⟨m|∂_t H |n⟩⟨n|/E_m-E_n,

one can see that , after having absorbed the geometric phases into the φ_n <cit.>, is purely off-diagonal in the basis of instantaneous eigenvectors (adiabatic basis). In terms of the Hilbert–Schmidt inner product (A^† B), this implies that  is orthogonal to H, and in general to all matrices commuting with H. It is also orthogonal to ∂_t H. Moreover, if we assume that the whole evolution takes place starting from an initial time t_i for a total time τ, one can rescale time according to s=(t-t_i)/τ and see that  scales like 1/τ. In other words, the faster the process, the stronger the field  should be. This property will be important in our discussion and it is particularly interesting in the context of quantum speed limits <cit.>.

By taking the derivative of Eq. (<ref>), one can obtain the relation

    ∂_t H(t) = i/ħ [H(t),(t)] + ∂ D(t),

where ∂ D(t) = ∑_n=1^N ∂_t E_n(t) |n(t)⟩⟨n(t)|.
Equation (<ref>) highlights the fact that ∂ D(t) generates the variation of the instantaneous eigenvalues of H(t), while  is responsible for the variation of eigenvectors. This is so since [∂ D(t),H(t)]=0, and then ∂ D makes H(t) “move” inside the set of Hamiltonians which commute with H at time t, without changing the eigenvectors. On the other hand, the term i[H,]/ħ determines the deviation from zero of the off-diagonal elements of H(t+dt) with respect to the instantaneous adiabatic basis at time t.






§ CONTROL FRAMEWORK
 

In order to realize , it is of great interest to understand how the structure (matrix components) of the Hamiltonian  is related to the structure of the initial Hamiltonian H and to the control resources. 
In this section we exploit techniques from control theory to show that, assuming that  can be computed, its action can be approximated at all times, arbitrarily well, by a time-dependent tuning of the initial parameters in the system Hamiltonian. This proves that the realization of our E–CD scheme is actually possible, and the specific construction will be discussed in the next section, Sec. <ref>.
Second, this study permits us to identify a class of Hamiltonians, that generalizes the set of real Hamiltonians, for which all matrix elements of  are not directly accessible without implementing new terms, which are not present in the original Hamiltonian.

Let H(t), of finite dimension N, be realized by tuning some available time-independent control Hamiltonians = {H_1,…, H_M } via a set of continuous control functions u(t) = {u_1(t), …, u_M(t) }. That is, H(t) can be expressed in the form 

    H{u(t) } = u(t)· =  ∑_k=1^M u_k(t) H_k.
 
Let -i be the vector of skew-Hermitian matrices {-i H_1,…, -i H_M }.
The matrices -i generate the so-called dynamical Lie algebra ℒ of the system <cit.>. This is the smallest algebra which contains -i, all possible commutators [-i H_j, -i H_k] of matrices from -i, all possible commutators of commutators, and so on, considering all possibly nested commutators [-i H_l, […,[-i H_j, -i H_k]]… ].

A basis of ℒ can be constructed by calculating, as a first step, the commutator of all possible pairs of matrices drawn from -i. Among the new obtained matrices, one should select those which are linearly independent from themselves and the original set, and compute the commutator of such new matrices with this original set. The procedure is repeated iteratively until no new linearly independent elements are produced. An explicit algorithm can be found in <cit.>. In our context, the linear span  can be though of as the set of all possible matrices attainable by H(t) at different times, for different values of the control functions.
 
Dynamical Lie algebras are of central importance in the study of the controllability of quantum systems <cit.>. This has its origin in the fact <cit.> that the set of reachable states, i.e., the set of unitary matrices that can be obtained as solution of the controlled Schrödinger equation for different choices of the control functions, coincides with the connected Lie group generated by ℒ.
An intuition for this can be given as follows. If an element - i H belongs to ℒ, then the Lie group element e^-i H t that it generates at time t can be realized, to arbitrary precision, by suitably concatenating group elements generated by Hamiltonians in the initial set -i, which belong to the one-parameter subgroups {e^-i H_1 t,…,e^-i H_M t}. From a control-theoretic perspective, this means that the evolution produced by H in a time t can also be obtained by a sequence of evolutions governed by the Hamiltonians , in general in a different total time.

Lie algebraic methods were also used in Ref.s <cit.>, in the context of STA, for designing feasible shortcuts connecting the same initial and final adiabatic states, but following different paths in the Hilbert space.

Assuming that all Hamiltonians involved are made traceless, ℒ can at most be 𝔰𝔲(N), the algebra of skew-Hermitian traceless matrices generating the group of special unitary matrices SU(N), and has dimension N^2-1. When ℒ=(N), the system is said to be (operator) controllable <cit.>. The Lie group generated by ℒ through the exponential map will be denoted by e^ℒ. 

Since -i (t) is a skew-Hermitian matrix, it must belong to (N) at all times. Therefore, if ℒ = (N) then -i(t) is obviously in ℒ. The first result we show is that ∈ℒ is a general property of , even when ℒ is not equal to (N), but is rather a smaller subalgebra ℒ⊂(N) – that is, even when the system is not fully controllable. This is an interesting result for two reasons. First of all, it allows us to restrict the class of Hamiltonians needed to identify and realize , with respect to the full set of traceless Hermitian matrices. Second, and more importantly in the present work, it means that the action of  can be always approximated, in the sense of the action of the Lie group, by working only with the initially available Hamiltonians.
This result is stated in the following theorem, whose proof is given in Appendix <ref>. In all the following results, we will generally assume that  can be computed. This excludes, in general, cases in which evolving eigenvectors can be degenerate <cit.>.

  Let the Hamiltonian of the system be expressible, at all times and for all values of the control functions u(t)={ u_1(t),…,u_M(t)}, as a linear combination H(t) = ∑_k=0^M u_k(t) H_k of time independent control Hamiltonians ={H_1,…,H_M}. Let ℒ be the Lie algebra generated by the matrices -i. Assuming the  exists, then -i (t) belongs to ℒ, for all times t.


When the Lie algebra ℒ has an additional structure, the general form of  can be characterized more accurately.
In order to do so, let us introduce Cartan decompositions, which are important tools in the study of controllability of quantum systems <cit.>. These are decompositions of the algebra ℒ into the direct sum form ℒ = 𝔥⊕𝔭 which satisfy the commutation relations 
    [𝔥,𝔥]⊆,    [,]⊆,    [,]⊆.

We can then state the following result (some of the quantities were defined in Sec <ref>).

 
Given a Cartan decomposition ℒ=⊕, if ∈ and -i ∂ D∈ for all t, then -i (t) ∈ for all t.


If -i∈, then also -i ∂_t H ∈, since it can be written as a linear combination -i ∂_t H = -i ∂_t u(t) ·. If also -i ∂ D ∈ for all t, then from Eq. (<ref>) it holds [H,] ∈. From commutation relations defining the Cartan decomposition, one can then conclude that -i ∈ for all t.



Let us clarify the result of Theorem <ref> by means of two examples.
The first one is the situation in which  includes exclusively real symmetric Hamiltonians. If this is the case, H(t) can always be diagonalized by a (real) special orthogonal matrix O(t). Thus (t) = i∂_t O O^† is necessarily a purely imaginary, skew-symmetric matrix, and -i (t) has thus no component on -i for all t.
When the set of control matrices generates the full algebra (N), this corresponds to a Cartan decomposition of the form 𝔰𝔬(N) ⊕ℑ, where 𝔰𝔬(N) is the algebra spanned by the real matrices in (N), while ℑ=𝔰𝔬(N)^⊥ is the set of purely imaginary matrices in (N).

The second example is the angular momentum algebra 𝔰𝔲(2) with  containing two arbitrary generators J_i,J_j among the three possible {J_1,J_2,J_3}. Explicit diagonalization and computation of  <cit.> shows that then  is proportional to the third generator.
This can be also seen from Eq. (<ref>): interpreting H geometrically as a vector H⃗(t) = {H_1,H_2,H_3} in the three-dimensional algebra, the matrix ∂ D generates the stretching of the Hamiltonian vector along in its own direction, while  generates rotations on the plane spanned by J_i and J_j; see Fig. <ref>. Therefore, recalling the properties of the rotations in 3D,  must lie in the direction of J_k. In this case the Cartan decomposition is unitarily conjugate to 𝔰𝔬(2)⊕ℑ.



§ EFFECTIVE CD FIELD
 

From Theorem <ref> one knows that (i)  can be approximated by a suitable choice of the control functions u(t) for all t and that (ii)  can be expressed as a linear combination of (nested) commutators of elements of . Observation (i) suggests that the problem of finding an effective CD field can be formulated as that of finding a correcting Hamiltonian , of fully controllable form 

    (t) = c(t)· = ∑_k=1^M c_k(t) H_k,

which emulates . More precisely, the Hamiltonian  should produce a dynamics  emulating the one induced by , i.e., .
While the algebraic structure suggests this can be done, it does not constructively specify what the c_k(t) need to be. To find them, from observation (ii), we choose to adopt a representation of the propagators  and  based on a Magnus expansion, see Appendix <ref>, where terms involving commutators of the Hamiltonians appear naturally. The first two terms are given for completeness,

    M^(1)(t)    = -i/ħ∫_0^t H(t_1) dt_1, 
    
    M^(2)(t)     = (-i/ħ)^2 ∫_0^t dt_1 ∫_0^t_1 dt_2 [H(t_1),H(t_2)].

The first step for designing the effective field  is to choose an ansatz for the control functions c(t) involving a certain number of free parameters. The second, considering a small evolution time t, is to ask the first terms of the Magnus expansion of  and  to coincide up to a desired order in t. This requirement will produce constraint equations for the free parameters.  
Since  and  cannot coincide at all times, one needs to choose a discrete temporal grid {t_k} for decomposing the whole evolution, and enforce the approximate equality at such time points. As a result, the approximating dynamics will match the true one at times {t_k}, while deviating at intermediate times.  

The matrix components of  which do not belong to  will not appear in the first Magnus term of H+, which is essentially the time average of H+. They will appear from the commutators in the following terms. Therefore, if one wants  to reproduce  effectively, it must hold that  gives no contribution to the first term, i.e., has vanishing time average on each interval (t_k,t_k+1). Besides, since the mediated effect will be of some order m greater than one in t, the magnitude of the control functions in  must be proportional to t^-X, for some X>0 depending on m, in order to amplify the mediated effect so that it acts at first order. 

Making these general ideas concrete, we choose the control functions in Eq. (<ref>) to have the form of a truncated Fourier series,

    c_k(t) = ω^X ∑_j=1^L [ A_k,jsin(j ω t) + B_k,jcos(j ω t)],

involving at most L harmonics of the fundamental frequency ω.
This way, the needed time grid is naturally defined by the set of stroboscopic times t_k = t_i+k T, k∈ℕ, where T = 2 π/ω is the period of the control functions and t_i the initial time. The constant amplitudes {A_k,j,B_k,j} are the free parameters used to enforce the constraints on the Magnus expansion at the end of each period. As a result, the control functions so determined will be discontinuous between different periods, due to the jumps in the values of the amplitudes. From a practical perspective, this is of course physically inconvenient, and so one would like to come up with an interpolation providing continuous and possibly smooth functions A_k,j(t), B_k,j(t). We will see that this can indeed be done, and often in a natural way.

The smaller the period T, the better the target dynamics  will be sampled, so in general  will need to oscillate fast with respect to the time-dependence of H. Raising the number of harmonics L permits us to introduce more parameters, and thus to obtain and solve constraint equations produced by higher Magnus terms. Ideally, one would like to find a good compromise between high sampling rates of , good approximation at stroboscopic times, while keeping the fundamental frequency and the number of harmonics as small as possible. This will be discussed in detail in Sec. <ref>. 
A strategy based on similar ideas was introduced in Ref. <cit.> for the purpose of producing effective time-independent Hamiltonians via optimal control.

Let us now focus on the class of Hamiltonians identified by Theorem <ref>. Due to the commutation relations from Eq. (<ref>), the Magnus terms of  involving an odd number of commutators, ^(2k), belong to . Those involving an even number of them, ^(2k+1), belong to  instead. Therefore, the general form of the constraint equations is

    ^(2k) =^(k), 
       ^(2k+1) = 0.


A discussion of the accuracy of the method is presented in Appendix <ref>, where an estimate of the expected error in the probability of not being in the target state, at the end of one period, is derived as a function of the number of solved constraint equations.
Let us show the overall procedure by working through a specific case.




 §.§ Worked through case: single spin
 
We concentrate here on the case of the Lie algebra 𝔰𝔲(2), which can physically describe, for example, a single spin driven by magnetic fields. A first application of STA methods via CD driving to this kind of system was studied in Ref. <cit.>. Since for tracking the instantaneous ground state it is sufficient to implement  without H, as discussed in Sec. <ref>, and for simplicity of presentation, we treat the case in which one wants to approximate  alone by means of . The general case, H+ approximated by H+, will be discussed at the end of this section. For clarity, let us work with the Pauli matrices {,,}, having commutation relations [σ_i,σ_j] = 2 i ϵ_ijkσ_k.
We assume that the Hamiltonian of the system is controllable only along two directions, say x and z, and can thus be written as:

    H(t) = u_x(t)  + u_z(t) .

The field , as discussed in Sec. <ref>, will be directed along the unavailable direction. It can be calculated explicitly <cit.> and it has the form (t) = (t), with:

    (t) = -1/2u_x(t) ∂_t u_z(t) - ∂_t  u_x(t) u_z(t)/u_x(t)^2+u_z(t)^2.

To achieve compensation to first order in T at stroboscopic times, we apply the E–CD Hamiltonian = c_x(t) +c_z(t), with control functions of the form (<ref>). The first non-zero term of the Magnus exponent produced by  is the second one, and has the desired matrix structure.   Generalizing to the choice of two arbitrary generators {σ_i,σ_j}, one has

    ^(2) (T) = - iϵ_ijk2π/ω^2-2X∑_n=1^L 1/n[A_i,n B_j,n-B_i,n A_j,n] σ_k .

From Eq. (<ref>) one can see that the highest order in T is produced by the interplay of sin and cos components related to the same harmonic n, belonging to different control functions c_i and c_j. There is no mixing between different harmonics, and no mixing between fields with the same phase. This feature is due to the properties of the integrals of the form ∫_0^T dt_1 sin(n ω t_1) ∫_0^t_1 dt_2 cos(k ω t_2).
Higher-order Magnus terms can be computed analytically, but the expressions get longer and longer with respect to Eq. (<ref>). 

Now, let us find the constraint equation to first order in T. For this purpose, two amplitudes are sufficient and, since we want ^(2) to be of order T, we choose X=1/2 in Eq. (<ref>).
A possible choice of control functions is then

    c_z(t) = A √(ω)sin(ω t);    c_x(t) = B √(ω)cos(ω t).

Using this functions, we can compute the first Magnus terms. Let t_n=t_i+nT, n=1,2,…, and let us indicate with ∫^T the integral over a full period ∫_t_n ^t_n+T. The exponents of the Magnus expansion for the CD and the effective dynamics at the end of each period n T of the control functions are, respectively,

    (T)     = -i ( ∫^T (t) dt ) , 
        = -i [(t_n+T/2) T + o(T^3)],
    (T)     = -i AB T  + o(T^3/2) .

The last equality in Eq. (<ref>) is obtained by formally Taylor-expanding (t) around the midpoint of the integration interval, t_n+T/2, and then integrating.
Equating order-T terms in Eqs. (<ref>) and (<ref>) one obtains the first constraint equation: AB = (t_n+T/2). This can be straightforwardly solved, and a possible solution for the n-th period is

    A = √(|(t_n+T/2)|),    B=[(t_n+T/2)] A,

where (x) indicates the sign of x and takes care of the case in which (t), and thus the product AB, is negative. Different solutions may be more suitable depending on the structure of f_CD.

The fact that the above solution for A and B, Eq. (<ref>), is a simple function of  evaluated at the midpoint of the period suggests a straightforward manner to interpolate the solutions in all intervals, obtaining continuous and smooth amplitude functions A(t) and B(t). It is indeed sufficient to replace (t_n+T/2) →(t) in Eq. (<ref>) and the same accuracy is maintained. As a conclusion, the effective field can be written as

    (t) = √(|(t)| ω)[ {(t)}cos(ω t)  + sin(ω t) ].
 
Two immediate observations can be made from Eq. (<ref>): First,  is proportional to √(ω), so, as expected, a better sampling can be obtained at the cost of having a stronger field. Second, due to the proportionality with √(|(t)|), the effective CD field inherits to some extent the behavior of the exact one. In particular, it vanishes when  does, that is, when the probability of nonadiabatic transitions is zero.

The third constraint equation, ^(3)=0, can be solved by adding a further term -4 A sin(2 ω t) to Eq. (<ref>). 

We conclude this section by discussing the more interesting case in which H+ approximates H+. This scenario is particularly important in a quasiadiabatic regime. When nonadiabatic effects are strong, the correcting fields become dominant with respect to the initial Hamiltonian H, and the method is highly nonadiabatic. In fact, even if the instantaneous eigenvectors of H are tracked during the evolution, the system is not in an instantaneous eigenstate of the total (corrected) Hamiltonian which is actually applied, H+ (or H +). On the other hand, if the nonadiabatic effects are weak, then one can think of being close to true adiabaticity.

Depending on the problem, it might be convenient to move to the time-dependent interaction picture for computing the effective field. In general, though, it becomes more complicated to compute the terms in the Magnus expansion. In any case, the treatment above can be repeated, and to order T^3 in infidelity, one still obtains the same effective field as in Eq. (<ref>).






§ APPLICATIONS
 

In this section, we discuss three applications of the method introduced in Sec. <ref>. The first one is the Landau-Zener-Majorana (LZM) model  <cit.>. This is the standard milestone in the study of nonadiabatic effects and more generally of the physics near avoided crossings in the energy spectrum.
The second one regards the adiabatic preparation of a two-qubit entangled state, a task of central importance for quantum information processing <cit.>.
The last application generalizes the first one and deals with a three-level system whose natural evolution induces the formation of a sequence of avoided crossings <cit.>. Such a model can describe local parts of many-body nontrivial energy spectra.

We also benchmark the efficiency of the E–CD method against standard (finite time) adiabatic driving. This is done with a focus on the LZM case, it being the simplest nonadiabatic scenario. This analysis underpins the usefulness of the method not only for theoretical purposes but, more importantly, for practical applications in the control of quantum systems.

For setting up the comparison between the E–CD and the adiabatic paradigms, a quantification of their respective performance is first needed. Since this is not an obvious task, let us here introduce the criteria which are adopted for this purpose.
The principle figure of merit we are interested in is the infidelity = 1-|⟨ψ(t_f)| gs(t_f)||⟩^2, i.e., the probability that the final state |ψ(t_f)⟩ of the system is not in the instantaneous ground state |gs(t_f)⟩ at the end of the protocol. We are also interested in making the evolution as fast as possible with as few resources as possible, but there is not a unique way to quantify this. The general idea is that one should not only take into account the full duration of the protocol, but also the intensity of the Hamiltonian which is applied to the system. Indeed, the adiabatic theorem <cit.> (see Appendix <ref>) formally states a property of the solutions of the family of Schrödinger equations 

    iħ∂_s U_τ(s) = τ H(s) U_τ (s),

for varying τ. Therefore, the theorem can be interpreted, on one side, as describing a varying-duration behavior, if Eq. (<ref>) is obtained after a rescaling s=t/τ of the physical time for a fixed Hamiltonian. Alternatively, it can be interpreted as describing an intensity-varying behavior, if Eq. (<ref>) is obtained from considering the duration fixed and amplifying the Hamiltonian like H(t)→τ H(t).

These considerations lead us to study the infidelity  both as a function of the protocol duration and of the ratio between the strength of  and that of the initial Hamiltonian H(t). We formalize this by choosing, as a measure of strength 𝒮(·), the maximal Frobenius norm over the whole evolution,

    𝒮(H) ≡max_t H(t)_F .

The Frobenius norm is defined as A_F=(A A^†)  and provides a quantification of the magnitude of the whole Hamiltonian matrix. Different choices are possible as cost functions depending on the resource one is interested in. For instance, one might be practically interested in the maximal matrix element, or the maximal or average amplitude reached by the control functions.

With these tools at hand, let us now proceed to the discussion of the above mentioned applications.



 §.§ Landau-Zener-Majorana model
 



We demonstrate the general approach detailed in Sec. <ref> for the specific case of the Landau-Zener-Majorana model <cit.>. Such a model describes a linear sweep of the gap between the energy levels of a two-level system. As they approach each other,
the presence of a coupling ħβ/2 prevents a net crossing, and an avoided crossing is produced instead with minimal gap ħβ. Assuming that the sweep spans an energy difference ħα in a time interval t_i ≤ t ≤ t_f, and that the anticrossing takes place at the intermediate time t_c=(t_f -t_i)/2, the Hamiltonian can be written in the form

    (t)/ħ =  α/2t-t_c/t_f-t_i + β/2.



In order to study the adiabatic properties of the system, it is convenient to rescale the quantities in terms of β, which defines a fundamental frequency scale of the problem. Furthermore, let us parametrize the time according to s = t-t_i/t_f-t_i, with 0 ≤ s ≤ 1. In terms of the dimensionless quantities s, ε = α/β and τ=β (t_f-t_i), the Schrödinger equation becomes 


    i ∂ U(s)/∂ s = τ/2[ε(s-1/2)  + ] U(s).


The parameter ε determines the distance, in units of β, from the avoided crossing at the beginning and at the end of the protocol. It also rules the decay rate of the LZM oscillations in the transition probability, i.e., of nonadiabatic transition, after the anticrossing <cit.> (see also Appendix A in <cit.>). Assuming that the system starts in the ground state, the transition probability to the excited state (i.e., of a nonadiabatic transition) is shown in Fig. <ref>. The asymptotic value, for large times, is given by the Landau-Zener formula <cit.>,

    P_LZ = exp(- πτ/2 ε).

The CD field for this problem, satisfying i ∂_s U = τ H_CD(s) U, is H_CD(s) =  f_CD(s) σ_y, with

    f_CD(s) = -1/2 τε/ε^2(s-1/2)^2+1.

This CD protocol was experimentally studied in Refs. <cit.>.
From the results of the previous Sec. <ref>, Eq. (<ref>), the dimensionless effective CD field, satisfying i∂_s U = τ H_E(s) U,  is 

    H_E(s) = √( |f_CD(s)| ω)[-cos(ω s τ) σ_x + sin(ω s τ)σ_z ],

with ω = 2 π/(β T).
The negative sign in front of the cosine is due to the fact that  is always negative.

The dynamics given by H alone can be compared with the one produced by  from Figs. <ref> and <ref>, for β = 1. As expected, the E–CD dynamics oscillates fast, remaining close to the target one. In the vicinity of the anticrossing, where nonadiabatic effects are strong, the deviation between consecutive periods is larger.

Let us now study the comparison between the E–CD and finite-time adiabatic methods more thoroughly. In particular, we show that the E–CD strategy can provide a consistent speed up of adiabatic protocols, characterizing the important regimes of parameters in terms of duration of the protocol and strength of the involved Hamiltonians. 
 
 
 
Specifically, we compute the infidelity obtained with the E–CD Hamiltonian for different total durations τ, after having selected the maximal frequency ω such that, for a certain factor k, the inequality

    𝒮() ≤ k 𝒮(H),

with 𝒮 defined in Eq. (<ref>), holds.

Let us remark that nonadiabatic transitions are exponentially weak in the adiabatic parameter <cit.>. On the other hand, one expects that the efficiency of the E–CD method increases polynomially when raising the frequency ω, since it is based on a perturbative argument. Therefore, one can safely predict that the adiabatic method, for very large τ and for a fixed maximal strength of the control fields, performs better. Nonetheless, what we are interested in is the intermediate range of durations, very far from the asymptotic limit in τ, which is the regime where practical protocols need to work. 




The results for the LMZ model are shown in Figs. <ref> and <ref>. In Fig. <ref>, the case in which only , as a stand-alone shortcut to adiabaticity, is applied is reported. As expected from the above discussion, the infidelity under effective CD driving scales like a power-law for large durations τ, and eventually the adiabatic exponential decay takes over. However, in the intermediate range of τ, the effective method does permit us to achieve an important speed-up of the adiabatic process, for fixed infidelity. The range of improvement gets wider and wider as one allows a stronger correcting field, while still remaining below the maximal norm of H. For example, one can reach a fidelity of ∼ 99.9% twenty times faster using an E–CD field of the same strength of H(t) [k=1 in Eq. (<ref>)].
Figure <ref> shows the same results for the case in which the E–CD Hamiltonian is implemented beside H(t), and thus acts as a weak correcting field to the major dynamics induced by H(t). The final infidelity oscillates strongly for different values of τ (thin gray line), so the points indicate an average of 20 surrounding points. The oscillations can bring the infidelity to much lower values, but never much higher than those represented by the points. This method performs worse in numbers than the previous case (Fig. <ref>) for small τ, but its main advantage is that it allows one to preserve true adiabaticity to a certain extent. That is, the driven states are close to being instantaneous eigenstates of the full Hamiltonian H+, since  is weaker than H and averages to zero in time. One can see that there is an initial regime in which the protocol is too fast for the correction to give a significant improvement, and a second regime in which the infidelity follows essentially the same behavior as in Fig. <ref>.
As an example, an infidelity of at least ∼ 10^-4 can be achieved with a speed up of around ∼ 6.5 times, using a maximal strength of the correcting field satisfying Eq. (<ref>) with k=1/2. A speed up of 2.2 times is obtained in the case k=1/4 for the same infidelity.
All these results confirm that the E–CD method represents an efficient strategy for obtaining fast adiabatic evolution, without requiring strong drivings.

The second criterium for comparing adiabatic and effective CD driving is the following: we consider the total duration to be fixed, 0 ≤ s ≤ 1, and we compare the dynamics produced by the amplified Hamiltonians H_τ(s) with the one produced by  for increasing values of the frequency ω. Due to the intertwined relation between duration and amplitude, we choose as a basis for comparison the integral of the norm over the whole evolution

    ∫_t_i^t_f dt H_X(t)_F .

The results are shown in Fig. <ref>. As for the previous criterium adopted, they are rather promising, confirming the validity of the E–CD method. Indeed, the E–CD Hamiltonian turns out to always produce a better infidelity, for a given total integral norm of the Hamiltonian, in the range of values studied.


The robustness of the method against static errors in the parameters, namely amplitude and relative phase of the fields, is discussed in Appendix <ref>. The results show that the method is not very sensitive to errors: it behaves linearly with respect to amplitude noise and quadratically with respect to phase noise.



 §.§ Two-qubit entanglement creation
 

The ability to reliably produce and manipulate entanglement is at the core of quantum information processing <cit.>. High-fidelity entanglement preparing protocols are thus a central touchstone for the development of quantum technologies. These are typically difficult to realize, especially due to the necessity of strong interactions and the ability of controlling many of them. 

Here, we use the E–CD method for the high-fidelity adiabatic preparation of a two-qubit entangled Bell state, requiring time-dependent control of local terms and one two-body interaction. 
Specifically, let us consider the Hamiltonian

    H(s)     = B(s) H_1 + g(s) H_2 ,
        = - B(s) (^(1)+^(2)) - g(s) (^(1)^(2) + ^(1)^(2)) .

The Hamiltonian H_1 is local, while H_2 has entangled eigenstates. Varying the field B(s) from high values to very small values, while keeping g fixed, the adiabatic path of the ground state connects a separable state |00⟩ to an entangled Bell state |ψ_+⟩=1/√(2)(|00⟩+|11⟩). In particular, we fix g=1, choose the simple time dependence B(s) = ε(1-s), and consider a total duration τ of the protocol.

The dynamical Lie algebra ℒ of the system has dimension four and is isomorphic to 𝔲(2), so the system is not fully controllable, since ℒ⊂(4). However, from Theorem <ref>, we already know that the Hamiltonian  will be inside ℒ, and will thus be a combination of the four basis elements. Since we know from Sec. <ref> that  is orthogonal, with respect to the Hilbert-Schmidt inner product, both to H(t) and its time derivative, we also know that it will not have components along H_1 and H_2. Indeed,  can be computed analytically (Appendix <ref>) and turns out to have the form

    τ(s) = 1/2ε/4 ε^2 (1-s)^2+1 H_3 = τ(s) H_3,

where H_3= ^(1)^(2)+^(1)^(2) is a hopping term. Therefore,  would require the time-dependent control of an extra two-body interaction. 

In order to avoid the additional implementation of H_3, we can use the results of Sec. <ref>, since H_3 ∝ [H_1,H_2]. We thus apply an E–CD field of the form (s) = c_1(s) H_1 + c_2(s) H_2, with control functions

    c_1 (s) = -√(|(s)|ω)cos(ω s τ),
    
    c_2(s) = √(|(s)|ω)sin(ω s τ).

The outcomes of a numerical simulation with parameters ε =5, τ = 5, N_T = ωτ/2 π = 10 are reported in Fig. <ref>. A final infidelity of 1.7 × 10^-3 is produced with a speed up of ∼ 12 times, for  having maximal strength 𝒮() ≤1/2𝒮(H) [see Eq. (<ref>)]. These results are extremely favorable, showing that the E–CD method provides a concrete advantage as a quantum control tool for a difficult and important task such as entanglement creation.



 §.§ Three-level system


We apply now the E–CD scheme to the case of a three-level system undergoing a sequence of LZM avoided crossing. 
The calculation and application of the exact CD field was discussed in detail in Ref. <cit.>. The Hamiltonian, in terms of dimensionless quantities, reads

    H(s) = [ d + ε(s -1/2)             1             0;             1          -2 d             1;             0             1 d - ε(s -1/2) ].




This Hamiltonian can physically describe an effective spin-1. For example, it could describe a local part of the energy spectrum of a molecular nanomagnet <cit.>, which is subject to a constant magnetic field along the x direction, a linearly sweeping magnetic field along the z direction, and an axial zero-field splitting <cit.>.
The Hamiltonian is real, so we are in the situation of Theorem <ref>.
The CD field  cannot be computed analytically for all times, but has the general form 

    = [        0 -i ^(12)  -i^(13);  i ^(12)        0 -i ^(23);  i ^(13)  i ^(23)        0 ],

with ^(12) and ^(23) always negative, while ^(13) assumes also positive values.
The effective field is chosen of the form 

    = √(ω)[  c_3(t)  c_1(t)       0;  c_1(t)       0  c_2(t);       0  c_2(t) -c_3(t) ].


The control functions c(t) are chosen such that the constraint equations given by ^(1) = ^(2) can be easily solved analytically, ensuring precision in infidelity at least to order T^3 at the end of one period. This is done by recalling, as observed in Sec. <ref>, that only sines and cosines of the same harmonics contribute to the first set of constraint equations. They read

    c_1(t)  =     A cos( ω t) - B cos(2 ω t) ,
    
    c_2(t) =     C sin( ω t) - D cos(3 ω t) ,
    
    c_3(t)  =    Bsin(2 ω t) + D sin(3 ω t).
 
The constraint equations are
 

    1/4 T B^2  = - ∫^T ^(12)(t) dt ,


    1/6 T D^2  = - ∫^T ^(23)(t) dt ,


    1/2 T A C = -∫^T ^(13)(t) dt .


Approximating the integrals with the value of the integrated functions multiplied by T, as explained in detail in Sec. <ref>, an interpolated solution for Eqs.
(<ref>) is given by

    A(t)    = - [^(13)] √(2 |^(13)|),     B(t) = 2 √(|^(12)|), 
    
    C(t)    = √(2|^(13)|),    D(t) = √(6 |^(23)|).

The results with parameters τ=25, ε=40, d = 5/2 are shown in Fig. <ref>. A speed up of 2.5 times is obtained, for equal strength 𝒮() = 𝒮(H) of the Hamiltonians. This shows that the E–CD paradigm behaves well also in a situation where more complex nonadiabatic phenomena are present in the dynamics of the quantum system.






§ CONCLUSION


We have introduced and discussed a method for speeding-up adiabatic quantum state transfer without requiring the introduction of new Hamiltonian components, while remaining always close to the true adiabatic path. 
This is achieved by introducing a control Hamiltonian which simulates in real time the effect of a counterdiabatic field <cit.>, but can be adsorbed into the initial Hamiltonian H(t).
As a trade-off, complete control of H(t) is required. However, the algebraic results of Sec. <ref> can be exploited to limit the number of independent control Hamiltonians needed in H(t) to implement the control protocol. For instance, we have shown in Sec. <ref> that, by encoding a two-level problem into a two-qubit framework, a two-qubit entangled state can be prepared with high fidelity in short time, requiring control of a single two-body interaction and a local field. 

Effective counterdiabatic Hamiltonians have been computed for the prototype applications studied, by choosing control functions which produce simple constraint equations for the amplitudes; see Secs. <ref> and <ref>. Remarkably, this allowed us to solve the equations analytically, and to obtain satisfying results without resorting to heavy numerical methods. Nonetheless, the structure of the E–CD problem still leaves room for optimization of the parameters, so further improvement can be predicted by the hybridization with optimal control strategies <cit.>. Moreover, 
one could combine the E–CD method with a preselection of an optimal adiabatic protocol. For instance, one might first conceive a “local” adiabatic driving <cit.>, in which the evolution is accelerated where nonadiabatic effects are small and then slowed down in the vicinity of avoided crossings. This would provide an a priori improvement of the basic adiabatic protocol with respect, for example, to a simple LZM sweep. Such a protocol could then be ulteriorly sped up by means of E–CD corrections.

In particular, the simplicity and versatility of the E–CD method makes it a strong strategy for the control 
of adiabatic processes. Therefore, we expect it to be of practical interest for many applications, ranging from quantum state preparation to adiabatic quantum computing <cit.>, quantum thermodynamics <cit.>, and probing of quantum critical dynamics <cit.>.

The method has been developed for finite-dimensional quantum systems. Still, the scheme could in principle be applied to infinite-dimensional systems with discrete levels, or restricted to well-separated nonpathological parts of the spectrum. Pathological here refers to possible degeneracies or continuous spectra. The first case, which is problematic also for finite dimensions, would require addressing the decoupling of instantaneous eigensubspaces, rather then single instantaneous eigenvectors. On the other hand, continuous spectra would require a revision of the CD framework, starting from the adiabatic theorem itself  <cit.>. Concerning the more mathematical results, noncompactness of the dynamical Lie group of the system would break some of the assumptions used in our control-theoretic setting <cit.>, which should then be investigated more thoroughly.



F.P. is grateful for hospitality at Imperial College London and acknowledges partial funding from the Doctoral International Mobility Program of the University of Parma. B.D. acknowledges funding from the Engineering and Physical Sciences Research Council (UK) administered by Imperial College London via the Postdoctoral Prize Fellowship program. The authors are grateful to Riccardo Mannella for a critical reading of the manuscript.








§ PROOF OF THEOREM <REF>
 

Before proceeding with the proof of Theorem <ref>, and to fix notation, let us recall the basic formulation of the adiabatic theorem.

Considering the Schrödinger equation, i ħ∂_t U(t)=H(t) U(t), from initial time t_i, let us parametrize the physical time t according to t-t_i=s τ, 0≤ s ≤ 1, where τ=t_f - t_i is the total evolution time. The equation becomes 

    i ħ∂_s U_τ(s)=τ H (s) U_τ(s) .

(Adiabatic theorem <cit.>)
Let U_τ(s) be the solution of (<ref>). Then,

    U_τ(s)  - ∑_n e^-i/ħτ∫_0^s E_n(t(s')) ds' |n(t(s))⟩⟨n(0)|  = O(τ^-1),


In the adiabatic limit τ→∞, and for s=1, it holds exactly

    U_τ(1) τ→∞⟶ U_ad(t_i,t) ≡∑_n e^-i/ħ∫_t_i^t E_n(t') dt' |n(t)⟩⟨n(0)|.


Under the control-theoretic setup introduced in Sec. <ref>, we proceed with the proof of Theorem <ref>. 

 Let Ĥ_CD = ( t̂ ) be the time-independent matrix obtained by evaluating (t) at the time instant t=t̂.
In order to prove that ∈ℒ, one can equivalently show that either (i)  can be written at all times as a linear combination of elements of a basis of ℒ, or (ii) that the group elements generated by the matrices Ĥ_CD for all time instants t̂, e^-i Ĥ_CD t, can be written as a concatenation of elements of the group e^ℒ. We follow route (ii). 
With the initial time t_i set to zero and ħ set to one for simplicity of notation, let (t) be the solution of the Schrödinger equation 

    i  ∂_t (t) = (t) (t) = [(t)-H(t)] (t),

where the Hamiltonian H_corr was introduced in Sec. <ref>. Let us assume that a Magnus expansion (see Appendix <ref>) (t)=e^(t) can be formally written, where

    (t)     = -i ∫_0^t ((t_1)-H(t_1)) d t_1 
       -1/2∫_0^t dt_1 ∫_0^t_1 d t_2 [(t_1), (t_2)]+ …

Now, calling (t) and M(t) the Magnus exponents generated by  and H respectively, the terms in (<ref>) can be rearranged as to give

    (t)     = (t) - M(t) + R(t),

with

    R(t) =     -1/2∫_0^t dt_1 ∫_0^t_1 dt_2 [(t_1), H(t_2)] 
         -1/2∫_0^t dt_1 ∫_0^t_1 dt_2 [H(t_1), (t_2)] + …

By expanding in Taylor series around the midpoint t^*=t/2 of the (small) integration interval (0,t), see Eq. (<ref>) in Appendix <ref>, one can write 

    R(t) = -1/12( [^(1),H^(0)]-[^(0),H^(1)] ) t^3 + o(t^4),

where we have used the notation 

    H_X^(k) = 1/k!.∂^k H_X(t)/∂^k t|_t^*.

Therefore R(t)=o(t^3).

Now, considering small times t, and using repeatedly the Zassenhaus formula <cit.> to factorize the exponential exp{(t) },  can be decomposed like

    (t) = (t) U^†(t) e^-1/2[M_CD(t),M(t)] e^o(t^3).

The unitary matrix (t), introduced in Sec. <ref>, is the solution of i ∂_t  = H_corr. 
Proceeding as above, one can estimate -1/2[M_CD(t),M(t)]=o(t^2).
The adiabatic theorem, Eqs. (<ref>) and (<ref>), then states that (t) can be approximated to arbitrary precision through an adiabatic process of long duration τ,

    (t) =(t) U^†(t) e^o(t^2)  + o(τ^-1).

We are then ready to prove the theorem. First of all, let us observe that, for sufficiently small t, one can rewrite

    Ĥ_CD = ∂/∂t̂∫_0^t̂(t_1) dt_1 = 1/t∫_t̂^t̂+t(t_1) dt_1 + o(t).

Applying the exponential map to Eq. (<ref>) and using Eq. (<ref>) we eventually obtain

    e^-i Ĥ_CD t    = (t̂,t̂+t)e^o(t^2), 
        = U_ad(t̂,t̂+t)U^† (t̂,t̂+t)e^o(t^2)+o(τ^-1).

Since U_ad and U can be realized via time dependent control of H(t), they belong to the group e^ℒ.
Equation (<ref>) then means that, for all times t̂, the group element exp(-i Ĥ_CD t) generated by -i Ĥ_CD can be approximated by a concatenation of elements of the group e^ℒ, arbitrarily well for sufficiently small t and large τ. That is, Ĥ_CD∈ℒ for all t̂.




§ MAGNUS EXPANSION
 


 §.§ Basics and useful bounds


The Magnus expansion <cit.> is a representation of the evolution operator U(t), solution of the time-dependent Schrödinger equation i ħ∂_t U(t) = H(t) U(t), in exponential form U(t) = exp{M(t)}. The exponent M(t) is an infinite sum M(t) = M^(1)(t)+M^(2)(t) + … whose first terms read

    M^(1)(t) =     -i/ħ∫_0^t H(t_1) dt_1, 
    
     M^(2)(t) =   (-i/ħ)^2 ∫_0^t d t_1∫_0^t_1 d t_2  [H(t_1),H(t_2)], 
    
     M^(3)(t)  =    (-i/ħ)^3 ∫_0^t d t_1∫_0^t_1 d t_2 
       ∫_0^t_3 dt_3 ( [H(t_1),[H(t_2),H(t_3)]] 
        +[H(t_3),[H(t_2),H(t_1)]]).


We will in general denote with M_X^(k) the k-th Magnus term for the Hamiltonian H_X.
The behavior of the Magnus terms for small integration time t can be conveniently studied by Taylor-expanding the Hamiltonian around the midpoint of the integration interval, and computing the integrals afterwards <cit.>. It turns out that, due to the time symmetry of the expansion, all terms are odd functions of t. In particular, one has

    M^(2 m)(t) = o(t^2m+1);    M^(2m+1)=o(t^2m+3).

Due to its usefulness in Appendix <ref>, for the proof of Theorem <ref>, we do the explicit calculation for the following “generalized” Magnus term, which involves the commutator of two different matrices,

    Ω_2(t)  = 1/2∫_0^t dt_1 ∫_0^t_1 d t_2 [A(t_1),B(t_2)].

Let us denote with t^*=t/2 the midpoint, and let us introduce the notation 

    a_k = 1/k!.∂^k A(x)/∂^k x|_t^*, b_k = 1/k!.∂^k B(x)/∂^k x|_t^*,

Inserting the Taylor expansions A(t) = ∑_k a_k (t-t^*)^k and B(t) = ∑_k b_k (t-t^*)^k into Eq. (<ref>) one obtains

    Ω_2(t)     = 1/2∑_k,n^0,∞1/n+1{1-(-1)^k+n+2/k+n+2
        - (-1)^n+1[1-(-1)^k+1]/k+1} [a_k,b_k] (t/2)^n+k+2, 
        = 1/4 [a_0,b_0] t^2 + 1/24( [a_1,b_0]-[a_0,b_1] ) t^3 
        + 1/48([a_0,b_2]+[a_2,b_0] ) t^ 4 + o(t^5).




 §.§ Error in the infidelity for the E–CD method
 

It is useful in our study to estimate the error when the expansion is truncated. In particular, we will focus on the error in the infidelity at the end of one period T of the oscillating control functions, for the problem of approximating  by means of , and for the class of Hamiltonians characterized by Theorem <ref> (see the end of Sec. <ref>). We will assume that the system starts in the ground state |gs(t_i)⟩ at the initial time t_i. Writing the states in terms of the evolution operators, and denoting with ⟨·⟩ the expectation value ⟨gs(t_i)|·|gs(t_i)⟩, the infidelity can be written in the form 

    = 1- |⟨^†⟩ |^2.
 
We then write the propagators in terms of their Magnus expansions, ^†(T)=e^-(T) and (T) = e^(T), and use the Baker-Campbell-Hausdorff formula for computing the product ^†. Assuming that one has solved the constraint equations, discussed in Sec. <ref>, involving the first 2m Magnus terms of , then it holds 

    ∑_k=1^m ^(k) = ∑_k=1^2m^(k).

We can then write

    ⟨^†⟩ =   ⟨exp{^(2m+1) + ^(2m+2)-^(m+1). . 
        -1/2[^(1),^(2m+1)] 
       . . -1/2[^(m+1),^(2)] + …}⟩ .

Let us remark that, since the oscillations in  are much faster than the timescales of , the behavior in Eq. (<ref>) does hold for (T), that is
^(2 m)(T) = o(T^2m+1) and ^(2m+1)(T)= o(T^2m+3).
For the dynamics , though, the integration is over a full period and the same is thus not true. Recalling the prefactor √(ω) in , the E–CD Magnus terms are of order ^(m)(T) = o(T^m/2).  Equation (<ref>) can then be expressed like

    ⟨^†⟩ = ⟨exp{^(2m+1) + o(T^m+1)}⟩.

One can now expand the matrix exponential in Taylor series and take the expectation value termwise. Taking the modulus squared, for a generic Hermitian matrix X,  it holds

    |⟨ e^i X t⟩|^2 = 1 +  t^2 ( Re⟨ X^2 ⟩ + |⟨ X ⟩|^2 ) + o(t^3),

Applying the same reasoning to Eq. (<ref>) and inserting into Eq. (<ref>) one finally obtains

    (T)     = - Re⟨( ^(2m+1))^2 ⟩ - |⟨^(2m+1)⟩|^2 + o(T^2m+3/2), 
        = o(T^2m+1).




§ ROBUSTNESS


We test the robustness of the E–CD method against possible experimental imperfection for the LZM case discussed in Sec. <ref>. We take into account possible errors in the amplitude and relative phase of the driving fields. This is done by adding a static relative offset δ. More specifically, with reference to Eq. (<ref>),

    sin(ω s τ)    → (1+δ) sin(ω s τ), 
    sin(ω s τ)    →sin(ω s τ + 2 πδ),

with -1/2 ≤δ≤ 1/2. The quantity F = |⟨gs(t_f) |ψ(t_f)||⟩^2 is the final fidelity between system state |ψ(t_f)⟩ and target (ground) state |gs(t_f)⟩, while F_0 is the fidelity in the absence of errors. The behavior of the relative error |1-F/F_0| is shown in Fig. <ref> as a function of the offset δ for the two error situations described by Eq. (<ref>). Both positive and negative values of δ are considered. The E–CD method results in being mostly linearly stable against amplitude errors, while it turns out to be quadratically sensitive to phase errors. It should be remarked that in some cases small errors in amplitude can incidentally lead to a small improvement of the method: this can be seen from the two regimes observable for positive δ. The negative peak around 10^-1 indicates the transition between a region δ<10^-1 where the protocol improves and a region for greater values where it performs worse instead. The behavior is symmetric with respect to the sign of δ for phase noise.








§ CD FIELD FOR THE TWO-QUBIT PROBLEM
 
In this appendix, we report the explicit computation of the CD field for the two-qubit problem described in Sec. <ref>. The initial Hamiltonian in dimensionless units, see Eq. (<ref>), is

    H(s)/ħ = -ε(1-s) [^(1)+^(2)]-[^(1)^(2)+^(1)^(2)].

As a first step, we need to diagonalize H(s). First of all, H(s) can be partially diagonalized by writing it in the (time-independent) basis

    ℬ_1 = {|00⟩; |01⟩+|10⟩/√(2); |10⟩-|01⟩/√(2); |11⟩}.

Let Q be the change-of-basis matrix from basis ℬ_0 ={|00⟩, |01⟩,|10⟩,|11⟩} to ℬ_1, thus having the kets in ℬ_1 as columns. 
The Hamiltonian reads

    [ -1+2 ε(s-1)           0           0          -1;                       0           0           0;           0           0           2           0;          -1           0           0  -1-2ε(s-1) ].

We thus see that two levels are actually decoupled among themselves and from the rest of the spectrum.
What remains to be diagonalized is a two-by-two real symmetric matrix (formed by the four corner elements). For this we can use the usual convenient trigonometric parametrization of two-by-two unitary matrices <cit.>, and the matrix diagonalizing (<ref>) is

    P(s)=[ -sin[θ(s)]          0          0  cos[θ(s)];          0          1          0          0;          0          0          1          0;  cos[θ(s)]          0          0  sin[θ(s)] ],

with θ(s) = 1/2arctan[1/2 ε (1-s)].
Eventually, the full diagonalizing matrix is U(s)=P(s) Q, and thus  the CD field can be readily computed, recalling that s=(t-t_i)/τ,

    (s)     = i/τ∂ U(s)/∂ s U(s)^†, 
    
        = 1/2 τε/1+4 ε^2 (s-1)^2 [^(1)^(2)+^(1)^(2)].










51
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Messiah(1961)]messiah1961qm
  author author A. Messiah, https://books.google.co.uk/books?id=mwssSDXzkNcC title Quantum Mechanics (publisher Dover
  Publications, New York, year 1961)NoStop
[Kato(1950)]kato1950
  author author T. Kato, 10.1143/JPSJ.5.435 journal
  journal Journal of the Physical Society of Japan volume 5, pages 435 (year
  1950)NoStop
[Demirplak and Rice(2003)]rice1
  author author M. Demirplak and author S. A. Rice, 10.1021/jp030708a journal journal J. Phys. Chem. A volume 107, pages 9937 (year 2003)NoStop
[Berry(2009)]berry1
  author author M. V. Berry, http://stacks.iop.org/1751-8121/42/i=36/a=365303
  journal journal J. Phys. A volume 42, pages 365303 (year
  2009)NoStop
[Bason et al.(2012)Bason,
  Viteau, Malossi, Huillery,
  Arimondo, Ciampini, Fazio,
  Giovannetti, Mannella, and Morsch]bason1
  author author M. G. Bason, author M. Viteau,
  author N. Malossi, author P. Huillery, author
  E. Arimondo, author
  D. Ciampini, author
  R. Fazio, author V. Giovannetti, author R. Mannella,  and author O. Morsch, http://www.nature.com/nphys/journal/v8/n2/full/nphys2170.html journal journal Nature Phys. volume
  8 (year 2012)NoStop
[Jarzynski(2013)]jarz1
  author author C. Jarzynski, 10.1103/PhysRevA.88.040101 journal journal Phys. Rev. A volume
  88, pages 040101 (year 2013)NoStop
[Demirplak and Rice(2008)]rice2
  author author M. Demirplak and author S. A. Rice, 10.1063/1.2992152 journal journal J. Chem. Phys. volume 129, pages 154111 (year 2008)NoStop
[Ibáñez et al.(2012)Ibáñez, Chen, Torrontegui,
  Muga, and Ruschhaupt]ibanez2012
  author author S. Ibáñez, author X. Chen,
  author E. Torrontegui, author J. G. Muga,  and author A. Ruschhaupt, 10.1103/PhysRevLett.109.100403 journal journal
  Phys. Rev. Lett. volume 109, pages
  100403 (year 2012)NoStop
[Torrontegui et al.(2014)Torrontegui, Martínez-Garaot, and Muga]torrontegui2014
  author author E. Torrontegui, author S. Martínez-Garaot,  and author J. G.  Muga, 10.1103/PhysRevA.89.043408
  journal journal Phys. Rev. A volume 89, pages 043408 (year
  2014)NoStop
[Martínez-Garaot et al.(2014)Martínez-Garaot, Torrontegui, Chen, and Muga]martinez2014
  author author S. Martínez-Garaot, author E. Torrontegui, author X. Chen,   and author J. G. Muga, 10.1103/PhysRevA.89.053408 journal journal Phys. Rev. A volume 89, pages 053408 (year 2014)NoStop
[Deffner et al.(2014)Deffner, Jarzynski, and del
  Campo]deffner2014
  author author S. Deffner, author C. Jarzynski,
   and author A. del Campo,  10.1103/PhysRevX.4.021013 journal journal Phys. Rev. X volume 4, pages
  021013 (year 2014)NoStop
[Li and Chen(2016)]li2016
  author author Y.-C. Li and author X. Chen,  10.1103/PhysRevA.94.063411 journal journal Phys. Rev. A volume 94, pages 063411 (year 2016)NoStop
[Baksic et al.(2016)Baksic,
  Ribeiro, and Clerk]Baksic2016
  author author A. Baksic, author H. Ribeiro,   and author A. A. Clerk, 10.1103/PhysRevLett.116.230503 journal journal Phys. Rev. Lett. volume 116, pages 230503 (year 2016)NoStop
[Opatrný and Mølmer(2014)]opatrny201
  author author T. Opatrný and author K. Mølmer, http://stacks.iop.org/1367-2630/16/i=1/a=015025
  journal journal New Journal of Physics volume 16, pages 015025 (year
  2014)NoStop
[Landau(1932)]landau
  author author L. D. Landau, @noop  journal journal Phys.
  Z. Sowjetunion volume 2, pages 46 (year
  1932)NoStop
[Zener(1932)]zener
  author author C. Zener, 10.1098/rspa.1932.0165 journal
  journal Proc. R. Soc. A volume 137,  pages 696 (year 1932)NoStop
[Majorana(1932)]majorana
  author author E. Majorana, 10.1007/BF02960953 journal
  journal Nuovo Cimento volume 9,  pages 43 (year 1932)NoStop
[Shevchenko et al.(2010)Shevchenko, Ashhab, and Nori]nori
  author author S. Shevchenko, author S. Ashhab,
   and author F. Nori, https://doi.org/10.1016/j.physrep.2010.03.002 journal
  journal Physics Reports volume 492,  pages 1  (year 2010)NoStop
[Theisen et al.(2017)Theisen, Petiziol, Carretta, Santini, and Wimberger]Theisen2017
  author author M. Theisen, author F. Petiziol,
  author S. Carretta, author P. Santini,  and author S. Wimberger, 10.1103/PhysRevA.96.013431 journal journal Phys.
  Rev. A volume 96, pages 013431
  (year 2017)NoStop
[Campbell and Deffner(2017)]Campbell2017
  author author S. Campbell and author S. Deffner, 10.1103/PhysRevLett.118.100601 journal journal Phys. Rev. Lett. volume 118, pages 100601 (year
  2017)NoStop
[Deffner and Campbell(2017)]DeffnerCampbell2017
  author author S. Deffner and author S. Campbell, http://stacks.iop.org/1751-8121/50/i=45/a=453001
  journal journal Journal of Physics A:
  Mathematical and Theoretical volume 50, pages 453001 (year 2017)NoStop
[D'Alessandro(2007)]dalessandro2007
  author author D. D'Alessandro, https://books.google.co.uk/books?id=e5M0id5enzQC
  title Introduction to Quantum Control and Dynamics,  Chapman & Hall/CRC Applied Mathematics & Nonlinear Science (publisher CRC Press, Boca Raton, year 2007)NoStop
[Schirmer et al.(2001)Schirmer, Fu, and Solomon]Schirmer2001
  author author S. G. Schirmer, author H. Fu,  and  author A. I. Solomon, 10.1103/PhysRevA.63.063410 journal journal Phys. Rev. A volume 63, pages 063410 (year 2001)NoStop
[Dirr and Helmke(2008)]Dirr
  author author G. Dirr and author U. Helmke,  10.1002/gamm.200890003 journal journal GAMM-Mitteilungen volume 31, pages 59 (year 2008)NoStop
[Glaser et al.(2015)Glaser,
  Boscain, Calarco, Koch,
  Köckenberger, Kosloff, Kuprov, Luy, Schirmer, Schulte-Herbrüggen, Sugny, and Wilhelm]Glaser2015
  author author S. J. Glaser, author U. Boscain,
  author T. Calarco, author C. P. Koch, author
  W. Köckenberger, author
  R. Kosloff, author I. Kuprov, author B. Luy, author S. Schirmer, author T. Schulte-Herbrüggen, author D. Sugny,  and author F. K. Wilhelm, 10.1140/epjd/e2015-60464-1 journal journal The European Physical Journal D volume 69, pages 279 (year
  2015)NoStop
[Jurdjevic and Sussmann(1972)]Jurdjevic
  author author V. Jurdjevic and author H. J. Sussmann, http://www.sciencedirect.com/science/article/pii/0022039672900356
  journal journal Journal of Differential
  Equations volume 12, pages 313 
  (year 1972)NoStop
[Khaneja and Glaser(2001)]KHANEJA2001
  author author N. Khaneja and author S. J. Glaser, http://www.sciencedirect.com/science/article/pii/S0301010401003184
  journal journal Chemical Physics volume 267, pages 11  (year
  2001)NoStop
[Verdeny et al.(2014)Verdeny, Rudnicki, Müller, and  Mintert]Verdeny2014
  author author A. Verdeny, author L. Rudnicki,
  author C. A. Müller,  and  author F. Mintert, 10.1103/PhysRevLett.113.010501 journal journal Phys. Rev. Lett. volume 113, pages 010501 (year 2014)NoStop
[Chen et al.(2010)Chen,
  Lizuain, Ruschhaupt, Guéry-Odelin, and Muga]Chen2010
  author author X. Chen, author I. Lizuain,
  author A. Ruschhaupt, author D. Guéry-Odelin,  and author J. G. Muga, 10.1103/PhysRevLett.105.123003 journal journal
  Phys. Rev. Lett. volume 105, pages
  123003 (year 2010)NoStop
[Nielsen and Chuang(2010)]NielsenChuang
  author author M. A. Nielsen and author I. L. Chuang, 10.1017/CBO9780511976667 title Quantum Computation and Quantum Information (publisher Cambridge University Press, Cambridge, UK, year
  2010)NoStop
[Jansen et al.(2007)Jansen,
  Ruskai, and Seiler]Jansen2007
  author author S. Jansen, author M.-B. Ruskai,
   and author R. Seiler, 10.1063/1.2798382 journal journal
  Journal of Mathematical Physics volume 48,  pages 102111 (year 2007)NoStop
[Vitanov(1999)]vitanov1
  author author N. V. Vitanov, 10.1103/PhysRevA.59.988 journal
  journal Phys. Rev. A volume 59,  pages 988 (year 1999)NoStop
[Vitanov and Garraway(1996)]vitanovgarraway
  author author N. V. Vitanov and author B. M. Garraway, 10.1103/PhysRevA.53.4288 journal journal Phys. Rev. A volume
  53, pages 4288 (year 1996)NoStop
[Lim and Berry(1991)]limberry
  author author R. Lim and author M. V. Berry, http://stacks.iop.org/0305-4470/24/i=14/a=014 journal journal J. Phys. A volume
  24, pages 3255 (year 1991)NoStop
[Malossi et al.(2013)Malossi, Bason, Viteau, Arimondo, Mannella, Morsch, and  Ciampini]malossi1
  author author N. Malossi, author M. G. Bason,
  author M. Viteau, author E. Arimondo, author
  R. Mannella, author
  O. Morsch,  and author
  D. Ciampini, 10.1103/PhysRevA.87.012116 journal journal Phys.
  Rev. A volume 87, pages 012116
  (year 2013)NoStop
[Dykhne(1962)]dykhne1962
  author author A. M. Dykhne, http://www.jetp.ac.ru/cgi-bin/e/index/e/14/4/p941?a=list
  journal journal Sov. Phys. JETP volume 14, pages 941 (year
  1962)NoStop
[Davis and Pechukas(1976)]davis1976
  author author J. P. Davis and author P. Pechukas, 10.1063/1.432648 journal
  journal The Journal of Chemical Physics volume 64, pages 3129 (year
  1976)NoStop
[Collison et al.(2003)Collison, Murrie, Oganesyan, Piligkos, Poolton, Rajaraman, Smith, Thomson, Timko, Wernsdorfer, Winpenny, and McInnes]collison1
  author author D. Collison, author M. Murrie,
  author V. S. Oganesyan, author S. Piligkos, author
  N. R. J. Poolton, author
  G. Rajaraman, author
  G. M. Smith, author
  A. J. Thomson, author
  G. A. Timko, author
  W. Wernsdorfer, author
  R. E. P. Winpenny,  and author
  E. J. L. McInnes, 10.1021/ic034541b journal journal Inorg. Chem.  volume 42, pages 5293 (year
  2003)NoStop
[Roland and Cerf(2002)]Roland2002
  author author J. Roland and author N. J. Cerf, 10.1103/PhysRevA.65.042308 journal
  journal Phys. Rev. A volume 65,  pages 042308 (year 2002)NoStop
[E. Farhi and Sipser(2000)]fahri2000
  author author  E. Farhi,  J. Goldstone, S. Gutmann and author
  M. M. Sipser, https://arxiv.org/abs/quant-ph/0001106 journal journal arXiv:quant-ph/0001106v1  (year 2000)NoStop
[Albash and Lidar(2018)]albash2016
  author author T. Albash and author D. A. Lidar, 10.1103/RevModPhys.90.015002 journal journal Rev. Mod. Phys. volume
  90, pages 015002 (year 2018)NoStop
[Alicki and Kosloff(2018)]Alicki2018
  author author R. Alicki and author R. Kosloff, https://arxiv.org/abs/1801.08314 journal
  journal arXiv:1801.08314v2  (year
  2018)NoStop
[del Campo et al.(2014)del
  Campo, Goold, and Paternostro]Campo2014
  author author A. del
  Campo, author J. Goold,  and  author M. Paternostro, http://dx.doi.org/10.1038/srep06208 http://10.0.4.14/srep06208
  https://www.nature.com/articles/srep06208#supplementary-information
  journal journal Scientific Reports volume 4, pages 6208 (year
  2014)NoStop
[Greiner et al.(2002)Greiner, Mandel, Esslinger, Hänsch, and Bloch]Greiner2002
  author author M. Greiner, author O. Mandel,
  author T. Esslinger, author T. W. Hänsch,  and  author I. Bloch, http://dx.doi.org/10.1038/415039a http://10.0.4.14/415039a journal journal Nature volume 415,  pages 39 (year 2002)NoStop
[Polkovnikov(2005)]Polkovnikov2005
  author author A. Polkovnikov, 10.1103/PhysRevB.72.161201 journal journal Phys. Rev. B volume
  72, pages 161201 (year 2005)NoStop
[del Campo et al.(2012)del
  Campo, Rams, and Zurek]delcampo2012
  author author A. del
  Campo, author M. M. Rams,   and author W. H. Zurek, 10.1103/PhysRevLett.109.115703 journal journal Phys. Rev. Lett. volume 109, pages 115703 (year 2012)NoStop
[Saberi et al.(2014)Saberi,
  Opatrný, Mølmer, and del Campo]manybody2014
  author author H. Saberi, author T. Opatrný,
  author K. Mølmer,  and  author A. del Campo, 10.1103/PhysRevA.90.060301 journal journal Phys. Rev. A volume 90, pages 060301 (year 2014)NoStop
[Maamache and Saadi(2008)]Maamache2008
  author author M. Maamache and author Y. Saadi, 10.1103/PhysRevLett.101.150407 journal journal Phys. Rev. Lett. volume 101, pages 150407 (year
  2008)NoStop
[Genoni et al.(2012)Genoni,
  Serafini, Kim, and Burgarth]Genoni2012
  author author M. G. Genoni, author A. Serafini,
  author M. S. Kim,  and  author D. Burgarth, 10.1103/PhysRevLett.108.150501 journal journal Phys. Rev. Lett. volume 108, pages 150501 (year 2012)NoStop
[Magnus(1954)]Magnus1954
  author author W. Magnus, 10.1002/cpa.3160070404 journal
  journal Communications on Pure and Applied Mathematics  volume 7, pages 649 (year
  1954)NoStop
[Blanes et al.(2009)Blanes,
  Casas, Oteo, and Ros]Blanes2009
  author author S. Blanes, author F. Casas,
  author J. Oteo,  and author J. Ros, https://doi.org/10.1016/j.physrep.2008.11.001 journal
  journal Physics Reports volume 470,  pages 151  (year 2009)NoStop



