Effect of Phase Errors on a Quantum Control Protocol Using Fast Oscillations

Abstract

It has been recently shown that fast oscillating control fields can be used to speed up an otherwise slow adiabatic process, making the system always follow an instantaneous eigenvector closely. In applying this method though, one typically assumes perfect phase relations among the control fields. In this work, we discuss the effect of potential static phase errors. We show that the latter can in some cases produce higher fidelities, leading to an unexpected improvement of the method. This is shown numerically and explained via a perturbative expansion of the error produced by the control strategy. When high-precision phase control is accessible, the results suggest that the phases of the control field can be used as free parameters whose optimization can be beneficial for the control protocol.

1. Introduction

The capability of preparing and manipulating with high precision the state of a quantum system is essential for the study of quantum phenomena, both for the investigation of fundamental aspects and for the development of technological applications. One of the most interesting control methodologies relies on a fundamental result of Quantum Mechanics, namely the adiabatic theorem [1]. This states that, given a time-dependent Hamiltonian, if the rate of the time evolution is sufficiently slow and there are no level crossings, then a system initially in an instantaneous eigenstate always remains in the corresponding instantaneous eigenstate at future times.

Although the adiabatic method is simple and robust by construction, for practical high-precision applications it often requires too long total evolution timescales, which are ultimately limited by dissipation and decoherence effects. For this reason, many strategies have been recently proposed for speeding up the quantum adiabatic dynamics, which are referred to as “shortcuts to adiabaticity” [2]. These methods can produce the desired evolution exactly in arbitrarily short time, provided one is able to implement sufficiently strong control fields. However, the implementation of the accelerated protocol is often experimentally challenging. This is so since it typically requires time-dependent tuning of matrix elements which are not controllable in the original Hamiltonian [2,3], therefore demanding new resources and modifications of the experimental setup. Indeed, experimental proofs have been limited to few level systems and typically needed adaptations [4,5,6,7]. The problem can be circumvented by working in carefully selected interaction pictures and relaxing the constraint of following the adiabatic path at all times, i.e., by asking the evolution to match the adiabatic eigenstates only at the beginning and at the end of the protocol, see e.g., [8,9,10]. This, however, results in a completely nonadiabatic method, thus arguably dropping the advantage of using adiabatic strategies.

Recently, it was shown that it is possible to achieve both effective instantaneous following of the adiabatic states and no necessity to control new matrix elements at the same time [3]. This can be done by introducing fast oscillations in the parameters of the initial Hamiltonian which mimic the action of a full shortcut-to-adiabaticity Hamiltonian. Therefore, the necessity of new control knobs is traded for a high degree of controllability of those initially available and for an approximate rather than exact driving. This methodology has been theoretically exemplified in a precise circuit quantum electrodynamics (QED) experimental context in Ref. [11], where the experimental realizability is discussed in detail.

An interesting feature recognized in previous work [3,11] is that this scheme requires precise phase relations between oscillations in different Hamiltonian matrix elements, which may not be easily guaranteed experimentally. Nonetheless, it was also noticed that small phase offsets could interestingly lead to an unpredicted improvement of the method. For this reason, in this work we elaborate in more detail on this phenomenon, studying how small biases in the relative phases of the control fields affect the performance of the protocol. We confirm that improvements can be obtained and we explain their origin in relation to the procedure for constructing the accelerating control fields. This new understanding allows us to speculate on the potential use of the phases of the control elements as additional parameters for optimizing the control method, also opening a connection with optimal control theory.

After reviewing the general method in Section 2, we derive general system-independent expressions for the corrections to the final fidelity in Section 3, which are then applied to the case study of a two-level system in Section 4, connecting our results to previous work. In particular, it is explained that, in some cases, small phase shifts can be exploited to improve the control method without adding higher harmonics to the oscillating control functions. Our central result is Equation (26) highlighting the role of the phase biases on the efficiency of the protocol. In Section 5, it is then shown numerically that these mechanisms are present also in multilevel systems, by considering the case of a sped-up adiabatic protocol which produces entanglement between two qubits dispersively coupled to a resonator [11].

2. Effective Counterdiabatic Driving

In this section, the method introduced in Ref. [3] is reviewed. Let us begin by summarizing the main features. The fundamental idea is that of introducing fast oscillations and use them to perturbatively counteract nonadiabatic transitions. Technically, this is done by decomposing the evolution into time steps determined by the small fundamental period of oscillation T of the control functions. Then, the propagator is written down at the end of each time step as a Magnus (exponential) [12] perturbation expansion in small T. The time integration over a full period of the oscillating control functions averages out the oscillations, leaving an effective expression for the propagator which depends only on the oscillation amplitudes and on the phases. These parameters are then used as free parameters to enforce equality of the actual dynamics, to first orders in T, with the desired adiabatic one. In the following, the details are provided by formulating the problem in the instantaneous adiabatic frame and including potential phase shifts among the control functions.

General Setup

Let the Hamiltonian of the finite-dimensional system be $H\left(t\right)={\sum }_{n=0}^N{E}_n\left(t\right)|n_t\rangle \langle n_t|$, where $E_n\left(t\right)$ are the instantaneous eigenvalues, assumed to be never crossing, while $|n_t\rangle$ are the instantaneous eigenvectors. The Hamiltonian $H\left(t\right)$ is assisted by a correcting control Hamiltonian $H_{\mathrm{c}}\left(t\right)$. The time dependent Schrödinger equation governing the evolution of the system is then ($\hslash =1$)

$$i{\partial }_{t}U\left(t\right)=[H\left(t\right)+H_{\mathrm{c}}\left(t\right)]U\left(t\right).$$

(1)

Let us assume that both $H\left(t\right)$ and $H_{\mathrm{c}}\left(t\right)$ can be written as a linear combination of the same M time-independent Hamiltonians ${\left\{H_k\right\}}_{k=1,\dots ,M}$: this assumption guarantees that $H_{\mathrm{c}}$ can be formally absorbed into H, i.e., that $H_{\mathrm{c}}$ does not include initially unavailable matrix elements. The matrices $\left\{H_k\right\}$ are assumed to be linearly independent and they determine the matrix structure of H:

$$H\left(t\right)=\sum _{n=1}^M{h}_n\left(t\right)H_n,H_{\mathrm{c}}\left(t\right)=\sum _{n=1}^M{F}_n\left(t\right)H_n.$$

(2)

The correcting Hamiltonian $H_{\mathrm{c}}$ is also assumed to include only fast oscillating terms with fundamental frequency $\mathrm{\Omega }$ and a limited—small, in principle—number of harmonics $\pm \ell \mathrm{\Omega }$ with $1\le \ell \le L$ and L a natural number. Let us indicate with ${\mathbb{Z}}_L=\{-L,\dots ,-1,1,\dots ,L\}$ the set of integer numbers from $-L$ to L with zero excluded. One can thus write the control functions $F_n\left(t\right)$ in terms of a truncated Fourier series:

$$F_n\left(t\right)=\sqrt{\mathrm{\Omega }}\sum _{l\in {\mathbb{Z}}_L}{f}_{n,l}{e}^{il\mathrm{\Omega }t+i{\varphi }_{n,l}}.$$

(3)

Most relevant for the present work are the potentially different phases ${\varphi }_{n,l}$ occurring in different control functions. The proportionality with $\mathrm{\Omega }$ is chosen since we have in mind applications to avoided crossing scenarios; as will become clear in Section 4, for this kind of problems this choice allows to counteract nonadiabatic effects to order T. In fact, it permits one to enhance second order (in T) terms in the expansion of the propagator which can then be used to counteract first order undesired transitions.

For studying nonadiabatic effects, it is now convenient to move to a time-dependent frame defined by the instantaneous adiabatic basis, that is the basis of instantaneous eigenvectors multiplied by the dynamical and geometric phase factors [13],

$$A\left(t\right)=\sum _{n=0}^N{e}^{-i{\int }_0^t{E}_n\left(t^{\prime }\right)d{t}^{\prime }-{\int }_0^t{\mathrm{\Theta }}_n\left(t^{\prime }\right)d{t}^{\prime }}|n_t\rangle \langle n_0|,$$

(4)

where ${\mathrm{\Theta }}_n\left(t\right)=\langle n_t|{\partial }_t{n}_t\rangle$ is the effective vector potential which produces the geometric phase [14]. In this frame, the Hamiltonian can be written like $H^{\left(A\right)}\left(t\right)=H_{\mathrm{c}}^{\left(A\right)}\left(t\right)+H_{\mathrm{na}}^{\left(A\right)}\left(t\right)$ and the Schrödinger Equation (1) is thus

$$i{\partial }_t{U}_{\mathrm{A}}=H^{\left(A\right)}{U}_{\mathrm{A}},$$

(5)

where $U_{\mathrm{A}}=A^{†}U$ and, introducing ${\omega }_{nm}\equiv E_n\left(t\right)-E_m\left(t\right)-i[{\mathrm{\Theta }}_n\left(t\right)-{\mathrm{\Theta }}_m\left(t\right)]$ and using Equations (2) and (3),

$$H_{\mathrm{c}}^{\left(A\right)}=A^{†}{H}_{\mathrm{c}}A=\sum _{n,m}^{1,N}\sum _{k=1}^M\sum _{l\in {\mathbb{Z}}_N}{f}_{k,l}{e}^{il\mathrm{\Omega }t+i{\varphi }_{k,l}}{e}^{i{\int }_0^t{\omega }_{nm}\left(t^{\prime }\right)d{t}^{\prime }}\langle n_t|H_k|m_t\rangle |n_0\rangle \langle m_0|,$$

(6)

$$H_{\mathrm{na}}^{\left(A\right)}=i{\partial }_t{A}^{†}A=i\sum _{\underset{n\ne m}{n,m}}^{1,N}{e}^{i{\int }_0^t{\omega }_{nm}\left(t^{\prime }\right)d{t}^{\prime }}\langle {\partial }_t{n}_t|m_t\rangle |n_0\rangle \langle m_0|.$$

(7)

Recall that the controllable matrices $\left\{H_k\right\}$ were introduced in Equation (2). Using Equation (2) and the relation $\langle {\partial }_t{n}_t|m_t\rangle =\frac{\langle n_t|{\partial }_{t}H\left(t\right)|m_t\rangle }{E_{nm}\left(t\right)}$ with $E_{nm}\left(t\right)=E_n\left(t\right)-E_m\left(t\right)$, which can be obtained by taking the time derivative of the eigenvalue equation $H\left(t\right)|n_t\rangle =E_n\left(t\right)|n_t\rangle$, Equation (7) can be further rewritten like

$$H_{\mathrm{na}}^{\left(A\right)}=\sum _{\underset{n\ne m}{n,m}}^{1,N}\sum _{k=1}^M{e}^{i{\int }_0^t{\omega }_{nm}\left(t^{\prime }\right)d{t}^{\prime }}\frac{{\partial }_t{h}_k\left(t\right)}{E_{nm}\left(t\right)}\langle n_t|H_k|m_t\rangle |n_0\rangle \langle m_0|.$$

(8)

Let us now consider an evolution time equal to one period $T=2\pi /\mathrm{\Omega }$. A formal solution of Equation (5) can be written in exponential form as a Magnus expansion [12] (assuming that T is sufficiently small for the expansion to converge): $U_{\mathrm{A}}\left(T\right)=\exp \left[{\sum }_n{\mathcal{M}}^{\left(n\right)}\left(t\right)\right]$. The first terms are

$${\mathcal{M}}^{\left(1\right)}\left(T\right)=-i{\int }_0^T{H}^{\left(A\right)}\left(t_1\right)d{t}_1$$

(9a)

$${\mathcal{M}}^{\left(2\right)}\left(T\right)=\frac{{(-i)}^2}{2}{\int }_0^{T}d{t}_1{\int }_0^{t_1}d{t}_2\left[H^{\left(A\right)}\left(t_1\right),H^{\left(A\right)}\left(t_2\right)\right]$$

(9b)

$${\mathcal{M}}^{\left(3\right)}\left(T\right)=\frac{i}{3!}{\int }_0^{T}d{t}_1{\int }_0^{t_1}d{t}_2{\int }_0^{t_2}d{t}_3\left(\left[H^{\left(A\right)}\left(t_1\right),\left[H^{\left(A\right)}\left(t_2\right),H^{\left(A\right)}\left(t_3\right)\right]\right]\right.$$

(9c)

$$+\left.\left[H^{\left(A\right)}\left(t_3\right),\left[H^{\left(A\right)}\left(t_2\right),H^{\left(A\right)}\left(t_1\right)\right]\right]\right).$$

(9d)

The terms of the Magnus expansion generated by $H_{\mathrm{c}}^{\left(A\right)}$ alone will be denoted as ${\mathcal{M}}_{\mathrm{c}}^{\left(k\right)}$ with k relating to the number of nested commutators involved. Similarly, those generated by $H_{\mathrm{na}}^{\left(A\right)}$ will be denoted as ${\mathcal{M}}_{\mathrm{na}}^{\left(k\right)}$, while terms whose integrands are mixed commutators of $H_{\mathrm{c}}^{\left(A\right)}$ and $H_{\mathrm{na}}^{\left(A\right)}$ will be denoted by ${\mathcal{M}}_{\mathrm{mix}}^{\left(k\right)}$. As an example, ${\mathcal{M}}_{\mathrm{mix}}^{\left(0\right)}=0$ and

$${\mathcal{M}}_{\mathrm{mix}}^{\left(1\right)}=-\frac{1}{2}{\int }_0^{T}d{t}_1{\int }_0^{t_1}d{t}_2\left\{\left[H_{\mathrm{c}}^{\left(A\right)}\left(t_1\right),H_{\mathrm{na}}^{\left(A\right)}\left(t_2\right)\right]+\left[H_{\mathrm{na}}^{\left(A\right)}\left(t_1\right),H_{\mathrm{c}}^{\left(A\right)}\left(t_2\right)\right]\right\}.$$

(10)

Let us remark that, since our objective is that of achieving an expansion in T for the propagator, if the solution of the equation $i{\partial }_t{U}_c^{\left(A\right)}=H_{\mathrm{c}}^{\left(A\right)}{U}_c^{\left(A\right)}$ is known in closed form then it would be convenient to further move to the frame defined by $U_c^{\left(A\right)}$; this would remove from Equation (9) the ${\mathcal{M}}_{\mathrm{c}}^{\left(k\right)}$ terms, being re-summed, leaving an expansion in T as desired. However, this is typically not the case in general. The recipe of the control method here discussed is to compute the expressions (9) while keeping terms up to a certain desired order in T. Then, the amplitudes $\left\{f_{n,k}\right\}$ of the oscillations in Equation (3) are to be used to make vanish as many terms as possible for increasing powers of T.

3. System-Independent Expressions

3.1. Magnus Terms

Using the expressions of Equations (6) and (7) into Equation (9), one can compute the general expressions for the first order terms of the expansion in T of the propagator at the end of a time step. Since it is assumed that the oscillations in $H_{\mathrm{c}}$ are much faster than the timescales of the system, we will assume that all the quantities except the exponentials in Equation (6) are constant within the time interval T. This formally amounts to consider an expansion at $T/2$ for all of them, where only the value evaluated at $T/2$ is kept. This approximation allows us to compute the integrals. Recalling that $H_{\mathrm{c}}\propto \sqrt{\mathrm{\Omega }}$ (see Equation (3)) the results up to order ${\mathrm{\Omega }}^{-3/2}$ are of the form

$$\begin{array}{cc}{\mathcal{M}}_{\mathrm{c}}^{\left(0\right)}\left(T\right)=m_{\mathrm{c}}^{\left(0\right)}{\mathrm{\Omega }}^{-3/2}+o\left({\mathrm{\Omega }}^{-5/2}\right);\hfill & \hfill {\mathcal{M}}_{\mathrm{c}}^{\left(1\right)}=m_{\mathrm{c}}^{\left(1\right)}{\mathrm{\Omega }}^{-1}+o\left({\mathrm{\Omega }}^{-2}\right);\\ {\mathcal{M}}_{\mathrm{c}}^{\left(2\right)}\left(T\right)=m_{\mathrm{c}}^{\left(2\right)}{\mathrm{\Omega }}^{-3/2}+o\left({\mathrm{\Omega }}^{(-5/2}\right);\hfill & \hfill {\mathcal{M}}_{\mathrm{na}}^{\left(0\right)}=m_{\mathrm{na}}^{\left(0\right)}{\mathrm{\Omega }}^{-1}+o\left({\mathrm{\Omega }}^{-2}\right),\end{array}$$

where

$$m_{\mathrm{c}}^{\left(0\right)}=-2\pi i\sum _{n,m}^{1,N}\sum _{k=1}^M\sum _{l\in {\mathbb{Z}}_L}\frac{f_{k,l}}{l}{e}^{i{\varphi }_{k,l}}{\omega }_{nm}\langle n_t|H_k|m_t\rangle |n_0\rangle \langle m_0|$$

(11a)

$$m_c^{\left(1\right)}=-2\pi i\sum _{n,m}^{1,N}\sum _{k,k^{\prime }}^{1,M}\sum _{l=1}^L\frac{f_{k,l}{f}_{k^{\prime },-l}}{l}{e}^{i\left({\varphi }_{k,l}+{\varphi }_{k^{\prime },-l}\right)}\langle n_t|[H_k,H_{k^{\prime }}]|m_t\rangle |n_0\rangle \langle m_0|$$

(11b)

$$\begin{array}{cc}\hfill m_c^{\left(2\right)}=& -\frac{2\pi i}{3}{\displaystyle \sum _{n,m}^{1,N}}{\displaystyle \sum _{k,k^{\prime },k^{″}}^{1,M}}\left\{\sum _{l\in {\mathbb{Z}}_L}{\displaystyle \sum _{l^{\prime }=1}^L}\frac{3}{l{l}^{\prime }}{f}_{k,l}{f}_{k^{\prime },l^{\prime }}{f}_{k^{″},-l^{\prime }}{e}^{i({\varphi }_{k^{″},-l^{\prime }}+{\varphi }_{k^{\prime },l^{\prime }}+{\varphi }_{k,l})}\right.\hfill \\ \hfill & \left.+\sum _{l,l^{\prime }}{f}_{k,l}{f}_{k^{\prime },l^{\prime }}{f}_{k^{″},-l-l^{\prime }}{e}^{i({\varphi }_{k,l}+{\varphi }_{k^{\prime },l^{\prime }}+{\varphi }_{k^{″},-l-l^{\prime }})}\frac{1}{l(l+l^{\prime })}\right\}\langle n_t|\left[H_k,\left[H_{k^{\prime }},H_{k^{″}}\right]\right]|m_t\rangle |n_0\rangle \langle m_0|\hfill \end{array}$$

(11c)

$$m_{\mathrm{na}}^{\left(0\right)}=2\pi \sum _{\underset{n\ne m}{n,m}}^{1,N}\sum _{k=1}^M\frac{h_k^{\prime }(T/2)}{E_{nm}}\langle n_t|H_k|m_t\rangle |n_0\rangle \langle m_0|.$$

(11d)

Notice that $U_{\mathrm{A}}$ of Equation (5) can thus be written as

$$U_{\mathrm{A}}=\exp \left[\frac{m_{\mathrm{c}}^{\left(1\right)}+m_{\mathrm{na}}^{\left(0\right)}}{\mathrm{\Omega }}+\frac{m_{\mathrm{c}}^{\left(0\right)}+m_{\mathrm{c}}^{\left(2\right)}}{{\mathrm{\Omega }}^{3/2}}+o\left({\mathrm{\Omega }}^{-2}\right)\right].$$

(12)

Let us remark that the first order in ${\mathcal{M}}_{\mathrm{mix}}^{\left(1\right)}$ is $o\left({\mathrm{\Omega }}^{-5/2}\right)$.

3.2. Infidelity

Since our objective is that of studying the effect of the phases $\left\{{\varphi }_{k,j}\right\}$ on the efficiency of the method, we want to derive an expression which quantifies the deviations from adiabaticity as a function of such phases. In order to do so, it is convenient to study the probability that the system has deviated from the instantaneous eigenstate at the end of the time step T. This is defined by

$${\mathbb{I}}_{\mathrm{T}}=1-{\left|\langle {\psi }_{\mathrm{tg}}|\psi \left(T\right)\rangle \right|}^2,$$

and we call it the infidelity between the system state $|\psi \left(t\right)\rangle$ and target state $|{\psi }_{\mathrm{tg}}\left(t\right)\rangle$ at the end of one time step T. In order to express ${\mathbb{I}}_{\mathrm{T}}$ in our adiabatic frame, let us note that, given the initial system state $|{\psi }_0\rangle$, it holds that $|{\psi }_{\mathrm{tg}}\left(t\right)\rangle =A\left(t\right)|{\psi }_0\rangle$, with $A\left(t\right)$ defined in Equation (4), while for the evolving state $|\psi \left(t\right)\rangle$ it holds $|\psi \left(t\right)\rangle =U\left(t\right)|{\psi }_0\rangle$, with $U\left(t\right)$ defined by Equation (1). Therefore ${\mathbb{I}}_{\mathrm{T}}$ can be written like

$${\mathbb{I}}_{\mathrm{T}}=1-{\left|\langle {\psi }_0|U_{\mathrm{A}}\left(T\right)|{\psi }_0\rangle \right|}^2,$$

(13)

which is consistent with the fact that the adiabatic evolution in the adiabatic frame is by definition time independent.

In order to compute an expansion of Equation (13) for small T, let us consider the propagator $U_{\mathrm{A}}$ written in exponential form, and let us write the exponent as a power series in T with half-integer powers, $k=1,3/2,2,\dots$. The half-integer powers are necessary due to the proportionality of $H_{\mathrm{c}}$ with $\sqrt{\mathrm{\Omega }}$, see Equations (2) and (3). Thus,

$$U_{\mathrm{A}}\left(T\right)=\exp \left(\sum _k{u}_k{T}^k\right).$$

From Equation (12), one can immediately identify

$$u_1=m_{\mathrm{c}}^{\left(1\right)}+m_{\mathrm{na}}^{\left(0\right)};u_{\frac{3}{2}}=m_{\mathrm{c}}^{\left(0\right)}+m_{\mathrm{c}}^{\left(2\right)}.$$

(14)

For easing up the notation, let $\langle ·\rangle$ denote the average value $\langle {\psi }_0|·|{\psi }_0\rangle$. An expansion of $|〈U_{\mathrm{A}}\left(T\right)〉{|}^2$ can be obtained by expanding the matrix exponential for small T, taking the average value termwise and then compute the modulus squared. Let us remark that, since all the Magnus terms are skew-hermitian, the expectation value $〈{\mathcal{M}}^k〉$ is imaginary and thus $\Re 〈{\mathcal{M}}^{\left(k\right)}〉=0$ for all k. This implies in turn that $\Re 〈u_k〉=0$ for all k. Up to order $T^3$, the infidelity can then be written like

$$\begin{array}{cc}\hfill {\mathbb{I}}_{\mathrm{T}}=& -\left[\Re 〈u_1^2〉+{\left|〈u_1〉\right|}^2\right]T^2\hfill \\ \hfill & -\left[2\Re \left({〈u_1〉}^{*}〈u_{\frac{3}{2}}〉\right)+\Re 〈\{u_1,u_{\frac{3}{2}}\}〉\right]T^{5/2}\hfill \\ \hfill & -\left[\Re 〈\{u_1,u_2\}〉+\Re 〈u_{\frac{3}{2}}^2〉+{\left|u_{\frac{3}{2}}\right|}^2+2\Re {〈u_1〉}^{*}〈u_2〉\right]T^3.\hfill \end{array}$$

(15)

Let us now assume that the system is initially in an instantaneous eigenstate, say $|{\psi }_0\rangle =|j_0\rangle$. Then, using Equations (11), one can note that $〈m_{\mathrm{c}}^{\left(1\right)}〉=0$ and $〈m_{\mathrm{na}}^{\left(1\right)}〉=0$, so from Equation (14) it holds that $〈u_1〉=〈m_{\mathrm{c}}^{\left(1\right)}〉$, $〈u_{\frac{3}{2}}〉=〈m_{\mathrm{c}}^{\left(2\right)}〉$.

3.3. Real Hamiltonians

For real Hamiltonians, the expression (15) for the infidelity can be further simplified. Let us indicate with ${\mathrm{ad}}^{\left(k\right)}\left(X\right)$ the nested commutator $[X\left(t_{i_1}\right),[\dots ,[X\left(t_{i_k}\right),X\left(t_{i_{k+1}}\right)]]\dots ]$ involving k commutators, which appear in the Magnus expansion. Denoting with ${\mathcal{M}}_{\mathrm{mix}}^{\left(n\right)}\left(\mathrm{even}\right)$ the commutators in ${\mathcal{M}}_{\mathrm{mix}}^{\left(n\right)}$ which contain an even number of $H_{\mathrm{c}}^{\left(A\right)}$, we can then state

Observation 1.

For real $H\left(t\right)$ and for an initial state $|{\psi }_0\rangle =|j_0\rangle$, it holds that

(i) 

$〈{\mathcal{M}}_{\mathrm{c}}^{(2n-1)}〉=0$ for $n\in \mathbb{N}$;

(ii) 

$〈{\mathcal{M}}_{\mathrm{mix}}^{\left(n\right)}\left(\mathrm{even}\right)〉=0$.

Proof. 

(i) Due to the properties of the commutator under unitary transformations

$$\langle j_0|{\mathrm{ad}}^{(2n-1)}\left[H_{\mathrm{c}}^{\left(A\right)}\right]|j_0\rangle =\langle j_0|A^{†}{\mathrm{ad}}^{(2n-1)}\left[H_{\mathrm{c}}\right]A|j_0\rangle =\langle j_t|{\mathrm{ad}}^{(2n-1)}\left[H_{\mathrm{c}}\right]|j_t\rangle .$$

(16)

Since $H\left(t\right)$ is real, the matrices $\left\{H_k\right\}$ are real symmetric and so is $H_{\mathrm{c}}$. Then, the commutators ${\mathrm{ad}}^{(2k-1)}\left[H_{\mathrm{c}}\right]$ are real skew-symmetric and are thus orthogonal, with respect to the Hilbert–Schmidt inner product, to the whole subspace spanned by the $\left\{H_k\right\}$. From this it follows, since $|j_t\rangle$ is an eigenstate of a combination of $\left\{H_k\right\}$, that $\langle j_t|{\mathrm{ad}}^{(2n-1)}|j_t\rangle =0$ and thus that $〈{\mathcal{M}}_c^{(2n-1)}〉=0$.

(ii) For (real) symmetric $H\left(t\right)$, $A\left(t\right)$ of Equation (4) is a real orthogonal matrix. Therefore $-i{\partial }_{t}A{A}^{†}=A\left(i{\partial }_t{A}^{†}A\right)A^{†}$ is imaginary skew-symmetric and is thus orthogonal to the subspace spanned by the real symmetric matrices $\left\{H_k\right\}$. Then, it also holds that $H_{\mathrm{na}}^{\left(A\right)}=i{\partial }_t{A}^{†}A$ is orthogonal to the subspace spanned by the matrices $\left\{A^{†}{H}_{k}A\right\}$. Since the commutator of two matrices is symmetric only if one of the two matrices is symmetric and the other is skew-symmetric, the commutators in ${\mathcal{M}}_{\mathrm{mix}}^{\left(n\right)}\left(\mathrm{even}\right)$ are of the form $A^{†}XA$ with X skew-symmetric, and they are thus orthogonal to the subspace spanned by the matrices $\left\{A^{†}{H}_{k}A\right\}$. Their average value is then zero. □

Furthermore, we also have $\Re 〈\{u_1,u_{\frac{3}{2}}\}〉=0$ since $u_1$ is purely real skew-symmetric while $u_{3/2}$ is purely imaginary skew-symmetric. In conclusion, the infidelity of Equation (15) reduces to

$${\mathbb{I}}_{\mathrm{T}}=-\left[\Re 〈u_1^2〉\right]T^2-\left[\Re 〈u_{\frac{3}{2}}^2〉+{\Im }^2〈u_{\frac{3}{2}}〉\right]T^3$$

(17)

$$\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=-\left[\Re 〈{\left(m_{\mathrm{c}}^{\left(1\right)}+m_{\mathrm{na}}^{\left(0\right)}\right)}^2〉\right]T^2-\left[\Re 〈{\left(m_{\mathrm{c}}^{\left(0\right)}+m_{\mathrm{c}}^{\left(2\right)}\right)}^2〉+{\Im }^2〈m_{\mathrm{c}}^{\left(2\right)}〉\right]T^3.$$

(18)

4. Effect of Phase Shifts: Two-Level System

Using the general expressions reported in Section 3, we study in this section the effects of phase shifts ${\varphi }_{k,l}$ in the control functions in a simple yet important scenario, namely that of a two-level system. In particular, the case of a two-level avoided crossing problem is considered in detail, being this a fundamental phenomenon where non-adiabaticity mechanisms become of paramount importance [15]. An avoided crossing occurs when two levels whose energy difference is varied would cross if uncoupled, but, if they are coupled instead, they repel each other after reaching a minimal gap. This phenomenon was first described theoretically by Landau [16], Zener [17] and Majorana [18], and has been studied thoroughly in a variety of experiments, see e.g., [19,20,21,22,23]. Many complex nonadiabatic phenomena can also be reformulated by means of multiple single-avoided-crossing scenarios [15,24,25,26,27]. Our aim is to show that slight phase offsets $\left\{{\varphi }_{k,l}\right\}$ can lead to an improvement of the fast adiabatic method based on fast oscillations discussed in the previous sections, and to explain this effect.

Let $\mathit{\sigma }=\{{\widehat{\sigma }}_1,{\widehat{\sigma }}_2,{\widehat{\sigma }}_3\}$ be the Pauli matrices. Let us consider a two-state problem ($N=2$) with initial Hamiltonian $H\left(t\right)=\mathit{H}\left(t\right)·\mathit{\sigma }$ with $\mathit{H}\left(t\right)=\{h_1\left(t\right),0,h_3\left(t\right)\}$ ($M=2$). That is, we assume that H is controllable only along the ${\widehat{\sigma }}_1$ and ${\widehat{\sigma }}_3$ directions, while the ${\widehat{\sigma }}_2$ direction is not directly available due to experimental constraints. For instance, this can happen in planar superconducting circuits such as the transmon [28] where it is not easy to control both real and imaginary parts of the couplings between bare levels [7]. Note that the Hamiltonian $H\left(t\right)$ is real-valued, and we are thus in the case of Section 3.3. The correcting Hamiltonian $H_{\mathrm{c}}$, being a linear combination of the same control matrices ${\widehat{\sigma }}_1$ and ${\widehat{\sigma }}_3$, can be written in the form $H_{\mathrm{c}}\left(t\right)=\mathit{F}\left(t\right)·\mathit{\sigma }$ with $\mathit{F}\left(t\right)=\{F_1\left(t\right),0,F_3\left(t\right)\}$. As it has been suggested in Ref. [3] a first-order cancellation method can be achieved by choosing control functions of the form $F_1\left(t\right)=A\sqrt{\mathrm{\Omega }}\cos \left(\mathrm{\Omega }t\right)$ and $F_3\left(t\right)=B\sqrt{\mathrm{\Omega }}\sin \left(\mathrm{\Omega }t\right)$ and this will be explained in more detail in the following. Here, since we want to study “dephasing” effects between the control fields, we allow for the presence of phase offsets ${\varphi }_1$ and ${\varphi }_3$ by taking

$$F_1\left(t\right)=A\sqrt{\mathrm{\Omega }}\cos (\mathrm{\Omega }t+{\varphi }_1);F_3\left(t\right)=B\sqrt{\mathrm{\Omega }}\sin (\mathrm{\Omega }t+{\varphi }_3).$$

(19)

This choice, relative to Equation (3), implies that only the fundamental frequency is present, $L=1$, and it moreover holds that ${\varphi }_{k,l}=-{\varphi }_{k,-l}\equiv {\varphi }_k$, $f_{1,1}=f_{1,-1}=A/2$ while $f_{3,1}=-f_{3,-1}=B/\left(2i\right)$.

In order to compute an expression for the infidelity of Equation (18) to the first two orders ($T^2$ and $T^3$), one can use Equation (11). Moreover, one can use the exact expressions for the eigenvectors of $H\left(t\right)$, which, since $H\left(t\right)$ is real-valued, can be written in the parametrized form

$$|1_t\rangle =\left\{-\sin \chi ,\cos \chi \right\},|2_t\rangle =\left\{\cos \chi ,\sin \chi \right\},$$

(20)

with $2\chi =\arctan (h_1\left(t\right)/h_3\left(t\right))$ the azimuthal angle between $\mathit{H}\left(t\right)$ and the z axis. In this way, one determines, using Equation (11),

$$m_{\mathrm{c}}^{\left(0\right)}=\frac{4\pi i}{{\mathrm{\Omega }}^{3/2}}\left[A{h}_3\sin {\varphi }_1+B{h}_1\cos {\varphi }_3\right]{\widehat{\mathrm{\Sigma }}}_2,$$

(21)

$$\begin{array}{c}{m}_{\mathrm{c}}^{\left(1\right)}=\frac{2\pi i}{\mathrm{\Omega }}AB\cos \mathrm{\Phi }{\widehat{\mathrm{\Sigma }}}_2,\hfill \\ m_{\mathrm{c}}^{\left(2\right)}=\frac{4\pi i}{{\mathrm{\Omega }}^{3/2}}\frac{AB}{\sqrt{h_1^2+h_3^2}}\left[(A\sin {\varphi }_1{h}_3+B\cos {\varphi }_3{h}_1){\widehat{\mathrm{\Sigma }}}_3+(A\sin {\varphi }_1{h}_1-B\cos {\varphi }_3{h}_3){\widehat{\mathrm{\Sigma }}}_1)\right]\cos \mathrm{\Phi },\hfill \end{array}$$

(22)

$$m_{\mathrm{na}}^{\left(0\right)}\left(T\right)=\frac{2\pi i}{\mathrm{\Omega }}\eta {\widehat{\mathrm{\Sigma }}}_2,$$

(23)

where we have introduced the relative phase $\mathrm{\Phi }={\varphi }_1-{\varphi }_3$, the matrices

$${\widehat{\mathrm{\Sigma }}}_1=|1_0\rangle \langle 2_0|+|2_0\rangle \langle 1_0|,{\widehat{\mathrm{\Sigma }}}_2=-i(|1_0\rangle \langle 2_0|-|2_0\rangle \langle 1_0|),{\widehat{\mathrm{\Sigma }}}_3=|1_0\rangle \langle 1_0|-|2_0\rangle \langle 2_0|$$

and the nonadiabatic coupling $\eta =\langle 1_t|{\partial }_t{2}_t\rangle =\frac{\langle 1_t|{\partial }_{t}H(T/2)|2_t\rangle }{E_1-E_2}$. The latter can be computed from the expressions in Equation (20) to be $\eta =\frac{h_1{\partial }_t{h}_3-h_3{\partial }_t{h}_1}{2(h_3^2+h_1^2)}.$ Thus, the first constraint equation, $u_1\left(T\right)=m_{\mathrm{c}}^{\left(2\right)}+m_{\mathrm{na}}^{\left(1\right)}=0$, is

$$(iAB\cos \mathrm{\Phi }+i\eta ){\sigma }_2=0⟹AB\cos \mathrm{\Phi }=-\eta .$$

(24)

If one was not aware of the presence of the phase offset $\mathrm{\Phi }$, a solution for $\mathrm{\Phi }=0$ could be taken to be

$$A=\sqrt{\left|\eta \right|};B=-\mathrm{sign}\left(\eta \right)\sqrt{\left|\eta \right|},$$

(25)

as was done in Ref. [3]. This would be a good choice unless $\sqrt{\left|\eta \right(t\left)\right|}$ was not a smooth function, in which case different solutions may be more convenient. Insertion of Equation (22) into Equation (18) gives the infidelity at the end of one time step T. Assuming that the system is initially in the ground state, $|{\psi }_0\rangle =|1_0\rangle$, the result is

$${\mathbb{I}}_{\mathrm{T}}={\eta }^2{(1-\cos \mathrm{\Phi })}^2{T}^2+\frac{2}{\pi }\left\{{(A{h}_3\sin {\varphi }_1+B{h}_1\cos {\varphi }_3)}^2+\frac{{\eta }^2{\cos }^2\mathrm{\Phi }}{h_1^2+h_3^2}{(A{h}_1\sin {\varphi }_1-B{h}_3\cos {\varphi }_3)}^2\right\}{T}^3.$$

(26)

With this expression at hand, the effect of phase shifts in the initial control functions can now be discussed. For very small T, one can focus on the first order, $T^2$, in Equation (26). This determines the general behaviour; the infidelity oscillates for varying $\mathrm{\Phi }$ with a period of $2\pi$, reaching the maximum value for $\mathrm{\Phi }=\pi +2k\pi$, $k\in \mathbb{Z}$—that is, when the oscillating part of the control functions of Equation (19) are in perfect phase match. On the other hand, the first order term is minimum for $\mathrm{\Phi }=0$, i.e., for control function with totally off-phase oscillating part as chosen in Equation (19). Interestingly, the first order term turns out to be perfectly symmetric with respect to sign of $\mathrm{\Phi }$, so, if higher-order terms are neglected, no benefit can come for the presence of phase offsets, since in such a case it always holds that ${\mathbb{I}}_{\mathrm{T}}\left(\mathrm{\Phi }\right)\ge {\mathbb{I}}_{\mathrm{T}}\left(0\right)$.

Although the first term of the expansion of Equation (26) determines the global behaviour of ${\mathbb{I}}_{\mathrm{T}}$, more details can be obtained by inspecting the second term, of order $T^3$. First of all, one can see that this is no more symmetric with respect to the sign of $\mathrm{\Phi }$, in general. In fact, a term $\propto \sin \mathrm{\Phi }$ appears, breaking the symmetry unless ${\varphi }_1=0$. This asymmetry has important consequences on the efficiency of the method; indeed, for certain values of $\mathrm{\Phi }\ne 0$ the infidelity can be become smaller than ${\mathbb{I}}_{\mathrm{T}}\left(0\right)$ due to a compensation of the terms in Equation (26), leading eventually to an improvement of the fidelity.

Secondly, the $o\left(T^3\right)$ term depends explicitly on ${\varphi }_1$ and ${\varphi }_3$, so when this term becomes relevant the infidelity does not depend any more on the relative phase $\mathrm{\Phi }$ only. This means that, even once the relative phase is fixed, one may be able to tune ${\varphi }_1$ and ${\varphi }_3$ so to further minimize ${\mathbb{I}}_{\mathrm{T}}$.

Thirdly, the $o\left(T^3\right)$ term depends also on the choice of A and B as solutions of the constraint Equation (24). This induces a dependence on the choice of the adiabatic sweep function, via the nonadiabatic coupling $\eta$.

As a final remark, it must be noted that, for very small $\mathrm{\Phi }$, no information can in general be obtained from Equation (26) if $\mathrm{\Phi }$ enters a regime in which its magnitude is comparable with that of T. In this case, in fact, Equation (26) would no longer be a valid expansion, and more terms would need to be considered.

In order to corroborate the previous observations, we treat a specific example in the following section.

4.1. Application: Avoided Level Crossing

The standard scenario where nonadiabatic effects become of great importance is when an avoided crossing of energy levels occurs. A simple model to describe this phenomenon is to consider a two-level Hamiltonian of the form

$$H\left(t\right)=h_3\left(t\right){\widehat{\sigma }}_3+h_1{\widehat{\sigma }}_1$$

(27)

with $h_1$ constant in time while $h_3\left(t\right)$ is a sweep function such that $h_3\left(t_{\mathrm{c}}\right)=0$—that is, such that in absence of the coupling $h_1$ the levels would cross at $t=t_{\mathrm{c}}$. Depending on the rate of the sweep, which determines the speed of the evolution, adiabaticity may be satisfied or violated. The simplest choice for $h_3\left(t\right)$ is a linear sweep $h_3\left(t\right)=\alpha t$, as proposed by Landau, Zener and Majorana [15,16,17,18]. For this case, an exact solution of the Schrödinger equation for t going from $-\infty$ to $+\infty$ can be computed, for the system starting in the ground state, giving an asymptotic probability at time $+\infty$ that the system has jumped to the excited state. This is exponentially small in a certain adiabatic parameter which quantifies the speed of the evolution—namely $\exp (-\tau )$ with $\tau =\pi h_1^2/\alpha$. In practical finite-time problems though, more refined choices for $h_3\left(t\right)$ are much preferable if one wants to achieve the best fidelity at the end of the sweep in a certain total time (while always remaining as close as possible to true adiabaticity, i.e., to being in an instantaneous eigenstate). For this reason, here we will work with the specific choice

$$h_3\left(t\right)=h_1\tan \left[c_1\left(1-2\frac{t}{t_f}\right)\right],$$

(28)

with $c_1=\arctan (\Delta /h_1)$ where $2\Delta$ is the initial energy gap, $t_f$ is the total duration time of the sweep, and $t\in [0,t_f]$. This kind of sweep function was introduced in Ref. [11], following general recipes from Refs. [29,30], and guarantees that the instantaneous rate of the evolution satisfies an instantaneous adiabatic condition for all times.

Using the expressions (19) and (25), a correcting Hamiltonian for this problem is

$$H_{\mathrm{c}}\left(t\right)=\sqrt{\mathrm{\Omega }|{\eta }_t|}[-\mathrm{sign}\left({\eta }_t\right)\cos \left(\mathrm{\Omega }t\right){\widehat{\sigma }}_1+\sin \left(\mathrm{\Omega }t\right){\widehat{\sigma }}_3],$$

(29)

with ${\eta }_t=-h_1{\partial }_t{h}_3/2(h_1^2+h_3^2)$.

Introducing the dimensionless quantities $\tau =t_f{h}_1,\omega =\mathrm{\Omega }/h_1,\xi =t_f\eta$ and rescaling the physical time like $s=t/t_f$, $s\in [0,1]$, the Schrödinger equation can be written in the dimensionless form:

$$i{\partial }_{s}U\left(s\right)=\tau [H\left(s\right)+H_c\left(s\right)]U\left(s\right),$$

(30a)

$$H\left(s\right)=\tan \left[c_1(1-2s)\right]{\widehat{\sigma }}_3+{\widehat{\sigma }}_1,$$

(30b)

$$H_c\left(s\right)=\sqrt{\frac{\omega \left|\xi \right|}{\tau }}[-\mathrm{sign}\left(\xi \right)\cos \left(\omega s\tau \right){\widehat{\sigma }}_1+\sin \left(\omega s\tau \right){\widehat{\sigma }}_3].$$

(30c)

The numerical results discussed in detail in the following are obtained with the values of the dimensionless parameters $\Delta /h_1=10,\tau =0.1$. The instantaneous energy levels and the shape of the sweep function of Equation (28) for these parameters are plotted in Figure 1a, while the evolution given by the control method, for $\omega /2\pi =15$ is compared with the adiabatic dynamics in Figure 1b.

Figure 2 and Figure 3 show the dependence of the infidelity on the phase offset $\mathrm{\Phi }$ for a sweep function of the form of Equation (28) and for fixed ${\varphi }_3=0$ (so $\mathrm{\Phi }={\varphi }_1$). Numerical data are indicated with symbols. In Figure 2, the general behaviour for an ample range of phase errors $-\pi \le \mathrm{\Phi }\le \pi$ is shown. It is appreciable how this is very well described by the oscillating term $\propto {(1-\cos \mathrm{\Phi })}^2$ of Equation (26). Nonetheless, such a description becomes ineffective for values of $\mathrm{\Phi }$ closer to zero. This can be seen from Figure 3a, where the same general behavior for $-\pi \le \mathrm{\Phi }\le \pi$ is shown for different values of T, in semilogarithmic scale; the blue dashed line represents the $o\left(T^2\right)$ term in Equation (26), which is a good description only when $\mathrm{\Phi }$ is not small. The behavior near zero is captured by taking into account also the $o\left(T^3\right)$ term (red solid line). The latter introduces an asymmetry around zero with respect to the sign of $\mathrm{\Phi }$ which becomes evident for small $\mathrm{\Phi }$: this is shown in Figure 3c for different values of T, with the blue line produced by considering both terms in Equation (26). As T increases, so does the asymmetry, due to the increasing importance of higher-order terms. From Figure 3c, it is also interesting to note that the asymmetry can produce a minimum which can be smaller than the value for $\mathrm{\Phi }=0$, witnessing again the fact that small phase offsets can be beneficial for the control method. Once the relative phase $\mathrm{\Phi }$ is fixed, one has still freedom to choose the phase ${\varphi }_1={\varphi }_3$; although the first term is invariant with respect to this choice, the second one is not. This leaves room for optimization of ${\varphi }_1$ for further minimizing the infidelity in each time step. As an example, in Figure 3b, the infidelity is shown for $\mathrm{\Phi }=0$ and varying ${\varphi }_1$. It is evident that some choices of ${\varphi }_1$ can strongly reduce the infidelity with respect to a null ${\varphi }_1$. Numerical simulations show that performing such an optimization in each time step can decrease the infidelity at the end of the whole protocol almost by one order of magnitude.

5. Beyond Two Levels: Two Atoms in a Cavity

In Ref. [11] an effective counterdiabatic scheme was proposed for speeding up an adiabatic entangling protocol between two superconducting qubits coupled to a transmission line resonator. Due to the complexity of the full multi-level problem, exact analytical expressions in the general case are hardly accessible. Therefore, we discuss the data produced by numerical simulation in the light of the results of the previous section, showing that an insight into the dependence on the phase shift can be obtained also in this scenario. The Hamiltonian of the problem in the rotating wave approximation is of Jaynes–Cummings form:

$$H={\omega }_r{\widehat{a}}^{†}\widehat{a}+\sum _{k=1}^2\frac{{\omega }_k}{2}{\widehat{\sigma }}_3^{\left(k\right)}+\sum _{k=1}^2{g}_k[{\widehat{\sigma }}_{+}^{\left(k\right)}\widehat{a}+{\widehat{\sigma }}_{-}^{\left(k\right)}{\widehat{a}}^{†}],$$

(31)

where ${\omega }_r$ and ${\omega }_k$—$k=1,2$—are the transition frequencies of the resonator and of the two qubits, respectively, while $g_k$ denotes the coupling between qubit k and the resonator. The idea of the entangling protocol is that of exploiting an avoided crossing which emerges between states $|n\uparrow \downarrow \rangle$ and $|n\downarrow \uparrow \rangle$, where n is the number of photons while the $|\uparrow \rangle$ denote the states of the qubits, when the two qubits get in resonance with each other far in frequency from ${\omega }_r$. This avoided crossing results in the dispersive regime ${\omega }_k-{\omega }_r\gg g_k$ from the emergence of an effective exchange interaction mediated by the exchange of virtual photons through the resonator. If the system starts in the state $|n\uparrow \downarrow \rangle$, for instance, driving the system adiabatically into qubit resonance produces the entangled singlet state $(|n\uparrow \downarrow \rangle -|n\downarrow \uparrow \rangle /\sqrt{2}$. In order to bring the two qubits into resonance, one drives their transition frequencies, so that the sweep function is $f\left(t\right)=({\omega }_1\left(t\right)-{\omega }_2\left(t\right))/2$. For the system starting in $|1\uparrow \downarrow \rangle$ and for $g_1=g_2\equiv g$, the correcting Hamiltonian $H_{\mathrm{c}}$, derived in Ref. [11], requires a time-dependent driving of the couplings g. It reads

$$H_{\mathrm{c}}=\sqrt{2\mathrm{\Omega }{\eta }_{2,3}}\left\{\sin \left(\mathrm{\Omega }t\right)[{\widehat{\sigma }}_{+}^{\left(1\right)}\widehat{a}+{\widehat{\sigma }}_{-}^{\left(1\right)}{\widehat{a}}^{†}]+\cos \left(\mathrm{\Omega }t\right)[{\widehat{\sigma }}_{+}^{\left(2\right)}\widehat{a}+{\widehat{\sigma }}_{-}^{\left(2\right)}{\widehat{a}}^{†}]\right\},$$

(32)

where ${\eta }_{2,3}=\left|\langle 2_t|{\partial }_t{3}_t\rangle \right|$ and $|2_t\rangle$ and $|3_t\rangle$ are the second and third excited instantaneous eigenvectors of the whole system. As opposed to the two-level case, the nonadiabatic coupling ${\eta }_{2,3}\left(t\right)$ requires numerical evaluation. This correcting Hamiltonian counteracts unwanted transitions between the two levels involved in the crossing, neglecting effects in the rest of the spectrum, since these are the major obstacle to the desired state transfer.

Let us now suppose that a phase shift is introduced such that, in the above Equation (32), $\cos \left(\mathrm{\Omega }t\right)⟶\cos (\mathrm{\Omega }t+\mathrm{\Phi })$. Figure 4 shows the dependence of the infidelity at the end of one time step T, in units of the unperturbed infidelity ${\mathbb{I}}_{\mathrm{T}}\left(0\right)$, as a function of the phase error $\mathrm{\Phi }$. The blue symbols are the results of numerical simulations with parameters ${\omega }_1\left(0\right)/2\pi =6.01$ GHz, $f\left(0\right)/2\pi =10$ MHz, $\mathrm{\Omega }/2\pi =7$ GHz, $g_1/2\pi =g_2/2\pi \equiv g/2\pi =10$ MHz and a sweep function of the form

$$f\left(t\right)=g_0\tan \left[\alpha (1-t/t_f)\right],$$

with $\alpha =\arctan [f\left(0\right)/g_0]$ and $2g_0$ the minimal width of the anticrossing extrapolated from numerical diagonalization. This is similar to the one in Equation (28), but it stops ($f\left(t_f\right)=0$) at the crossing. One can recognize from Figure 4 the same general oscillatory behaviour found in Section 4, which is well captured, far from its nodes, by the black solid line: this line represents a fit with functional dependence as in the first term of Equation (26), namely $k{[1-\cos \left(\mathrm{\Phi }\right)]}^2$ with k free parameter. In the vicinity of the values $\mathrm{\Phi }=\pm 2n\pi$ though, higher order effects again come into play. Focusing in the neighborhood of zero, an asymmetry with respect to the sign of $\mathrm{\Phi }$ is present, so that for small positive $\mathrm{\Phi }$ the infidelity becomes lower than the unperturbed value. We thus see once again that small phase offsets can lead to an accidental improvement of the control protocol. The systematics of this phenomenon further suggest that it can be exploited for optimizing the method.

In conclusion of this analysis, the parameter ranges should be discussed in which these effects are appreciably detectable. This is strictly related to the presence of other secondary phenomena which compete in slightly affecting the value of the infidelity on the same order of magnitude. For the specific framework treated in this section, the essential phenomenon which must be particularly taken into account is that of transitions towards levels other than the two involved in the crossing. We refer to Ref. [11] for a detailed discussion. The strength of this transition channel depends on the total duration of the protocol and on how good the dispersive-regime approximation is. The latter in turn depends on how large is the ratio $({\omega }_k-{\omega }_r)/g$. The effects due to the phase shift are related instead to the strength of second-highest-order term in the expansion of the infidelity in the time step T, see Equation (26), which in turn depends on the system parameters and especially on how small T is. Assuming that T is sufficiently small for the method to converge at the end of each time step, these effects become visible for larger T, while they are reduced as T diminishes as discussed in Section 4.1. In Figure 5 we show how the dependence of the infidelity on the phase offset $\mathrm{\Phi }$ behaves as the parameter g is varied, thus affecting the dispersive-regime approximation. For facilitating the comparison, all the curves are renormalized in such a way that the highest value reached in the interval considered lies at ordinate 1. It is immediately clear the general $\sim {(1-\cos \mathrm{\Phi })}^2$ oscillating behaviour is robust even if g becomes large (larger values would challenge the validity of the Hamiltonian (31) in the experimental platform considered here, due to effects beyond the strong-coupling regime between qubits and resonator [31]). For small $g\sim 3$ MHz, transitions to other levels can be neglected and the time step T is sufficiently small for the $\sim {(1-\cos \mathrm{\Phi })}^2$ behaviour to be completely dominant. For larger values of g, asymmetries in the vicinity of zero become more visible. For larger values of g the whole pattern start to shift with respect to $\mathrm{\Phi }=0$, indicating the emergence of further effects.

6. Discussion

We have discussed the effect of phase errors in the control functions for a control method which performs effective counterdiabatic driving by using fast oscillations. We have derived general expressions for the first order terms of the error produced by the protocol. These have been applied to the case of an avoided crossing problem, well explaining the behaviour found in numerical simulations. The results show that, on the scale of large values of the phase error $\mathrm{\Phi }$ (as compared to $\pi$), the fidelity oscillates as a function of $\mathrm{\Phi }$ and the protocol typically becomes sensibly less efficient. Nonetheless, it turns out to be rather stable against small values of $\mathrm{\Phi }$. Remarkably, small values of $\mathrm{\Phi }$ can also produce higher fidelities. This is due to the small phase error leading to an accidental compensation of the first terms in an expansion of the infidelity. As a result, nonadiabatic transitions are further suppressed with respect to the basic protocol and the method turns out to be more efficient. This mechanism, well understood in the case of a two-level system, is confirmed numerically also for a more complex multilevel system. Also in this scenario, in fact, the infidelity oscillates for large $\mathrm{\Phi }$, while it can diminish for small $\mathrm{\Phi }$.

These results are a good evidence that a fine tuning of the phases of the control field can be used for an optimization of the control protocol, without needing to include more expensive resources, such as new control fields or higher harmonics.