Accelerating adiabatic protocols for entangling two qubits in circuit QED

TL;DR

This paper presents a method to accelerate adiabatic entanglement protocols in circuit QED, achieving faster operation while maintaining high fidelity and robustness against errors and decoherence.

Contribution

It introduces a technique that adds fast oscillations to adiabatic sweeps, significantly speeding up entanglement protocols in circuit QED systems.

Findings

Substantial speedup over existing adiabatic methods.

Effective beyond the rotating-wave approximation.

Robustness against static biases and decoherence.

Abstract

We introduce a method to speed up adiabatic protocols for creating entanglement between two qubits dispersively coupled to a transmission line, while keeping fidelities high and maintaining robustness to control errors. The method takes genuinely adiabatic sweeps, ranging from a simple Landau-Zener drive to boundary cancellation methods and local adiabatic drivings, and adds fast oscillations to speed up the protocol while canceling unwanted transitions. We compare our protocol with existing adiabatic methods in a state-of-the-art parameter range and show substantial gains. Numerical simulations emphasize that this strategy is efficient also beyond the rotating-wave approximation and that the method is robust against random static biases in the control parameters and with respect to damping and decoherence effects.

Figures

Equations

$$\mathcal{H}=$$

$$+\sum _{k=1,2}{g}_k(b_k+b_k^{†})(a+a^{†}).$$

$$H_{\text{JC}}={\omega }_r{a}^{†}a+\sum _{k=1}^2\frac{{\omega }_k}{2}{\sigma }_z^{(k)}+\sum _{k=1}^2{g}_k[{\sigma }_{+}^{(k)}a+{\sigma }_{-}^{(k)}{a}^{†}],$$

$$H=\left(\begin{array}{cccc}\frac{{\delta }_1+{\delta }_2}{2}& g_2& g_1& 0\\ g_2& \frac{{\delta }_1-{\delta }_2}{2}& 0& \sqrt{2}{g}_1\\ g_1& 0& \frac{{\delta }_2-{\delta }_1}{2}& \sqrt{2}{g}_2\\ 0& \sqrt{2}{g}_1& \sqrt{2}{g}_2& -\frac{{\delta }_1+{\delta }_2}{2}\end{array}\right).$$

$$f(t)=({\delta }_1-{\delta }_2)/2g$$

$$i{\partial }_{s}U(s)=\tau H(s)U(s)$$

$$H(s)=\left(\begin{array}{cccc}{\delta }_g& 1& 1& 0\\ 1& f(s)& 0& \sqrt{2}\\ 1& 0& -f(s)& \sqrt{2}\\ 0& \sqrt{2}& \sqrt{2}& -{\delta }_g\end{array}\right).$$

$$H(s)={\delta }_g{D}_0+f(s)D_1+C_1+C_2.$$

$$f_{\text{LZ}}(s)=f_0[1-s].$$

$$f_{\text{PL}}(s)=f_0[1-35s^4+84s^5-70s^6+20s^7].$$

$$f_{\beta }(s)=f_0[1-{\Theta }_k(s)]$$

$${\Theta }_k(x)=\frac{B_s(1+k,1+k)}{B_1(1+k,1+k)},$$

$$f_{\text{RC}}(s)=\frac{g_0{f}_0(1-s)}{g\sqrt{(g_0/g{)}^2+f_0^{2}s(2-s)}},$$

$$f(s)=\frac{g_0}{g}\tan [\alpha (1-s)],$$

$${\mathbb{I}}_f=1-|⟨\psi (t_f)|e(t_f)⟩{|}^2.$$

$$H(t)=\sum _{n=1}^N{E}_n(t)|n_t⟩⟨n_t|.$$

$$U_{\text{d}}^{†}(t)H(t)U_{\text{d}}(t)=\text{diag}\{E_1(t),\dots ,E_N(t)\},$$

$$H_{\text{CD}}(t)=i{\partial }_t{U}_{\text{d}}(t)U_{\text{d}}^{†}(t).$$

$$H_{\text{CD}}=i\sum _{m\ne n}\sum _{n=1}^N\frac{|m_t⟩⟨m_t|{\partial }_{t}H|n_t⟩⟨n_t|}{E_n(t)-E_m(t)}.$$

$$H_{\text{CD}}^{(i,j)}(s)=\text{Im}⟨i|H_{\text{CD}}(s)|j⟩$$

$$H_{\text{p}}(s)=\frac{H_{\text{CD}}^{(2,3)}}{2}[{\sigma }_x^{(1)}{\sigma }_y^{(2)}-{\sigma }_y^{(1)}{\sigma }_x^{(2)}],$$

$$H_{\text{e}}=d_0(t)D_0+d_1(t)D_1+c_0(t)C_0+c_1(t)C_1.$$

$$\sum _k\left[{\mathcal{A}}_k\sin (k\omega t)+{\mathcal{B}}_k\cos (k\omega t)\right];$$

$$2i[C_1,C_2]={\sigma }_y^{(1)}{\sigma }_x^{(2)}-{\sigma }_x^{(1)}{\sigma }_y^{(2)}.$$

$$H_{\text{e}}=c_1(t)C_1+c_2(t)C_2$$

$$c_1(t)=\mathcal{A}\sqrt{\omega }\sin (\omega t);c_2(t)=\mathcal{B}\sqrt{\omega }\cos (\omega t).$$

$$M_{\text{e}}(T)=-i\frac{\mathcal{A}\mathcal{B}}{4}T[{\sigma }_x^{(1)}{\sigma }_y^{(2)}-{\sigma }_y^{(1)}{\sigma }_x^{(2)}]+o(T^{3/2}),$$

$$M_{\text{p}}(t)$$

$$=-i\frac{H_{\text{CD}}^{(2,3)}(T/2)}{2}T[{\sigma }_x^{(1)}{\sigma }_y^{(2)}-{\sigma }_y^{(1)}{\sigma }_x^{(2)}]+o(T^3).$$

$$\mathcal{A}(t)=\mathcal{B}(t)=\sqrt{2H_{\text{CD}}^{(2,3)}(t)}.$$

$$H_{\text{e}}(t)=\sqrt{2\omega H_{\text{CD}}^{(2,3)}(t)}\left[\sin (\omega t)C_1+\cos (\omega t)C_2\right].$$

Full text

I Introduction

The development of quantum technologies depends crucially on the capability of realizing many high-quality logical operations on a quantum system within the coherence times. There are thus two objectives to be realized simultaneously, time and fidelity, which are often the opposite sides of the same coin: Doing a protocol slowly often reduces the associated gate errors, but limits the number of operations that can be performed within the natural lifetime of the system. In addition to this, there is the constraint that a quantum control problem is always subject to limitations on the control resources, that is, on the knobs which one can work with in order to influence the system’s behavior.

If the time and control restrictions are not severe, adiabatic methods are a very interesting option due to their intrinsic robustness. They require a slow driving of the system Hamiltonian in time such that the system, as predicted by the adiabatic theorem of quantum mechanics Messiah (1961), always remains in an instantaneous eigenvector. However, adiabatic methods may be too slow to be a viable strategy.

On a different side of the quantum control spectrum, the most used protocols exploit sequences of evolutions under time-independent effective Hamiltonians. These are typically faster and require less control resources, but they are often unstable with respect to biases in control and system parameters.

In this work we contribute to fill the gap between these two scenarios, by discussing a control protocol producing high fidelities in intermediate-range timescales, while still partially inheriting the robustness of adiabatic strategies and without requiring the addition of new “control knobs.” The protocol is contextualized in the problem of producing entanglement between two qubits in a standard circuit QED (cQED) architecture, namely constituted of two superconducting qubits dispersively coupled to a planar transmission line resonator Blais et al. (2007); Gu et al. (2017).

We first compare different finite-time adiabatic sweeps, ranging from the standard linear Landau-Zener Shevchenko et al. (2010) to boundary cancellation methods (BCMs) Lidar et al. (2009) and local adiabatic drivings (LADs) Roland and Cerf (2002); Rezakhani et al. (2009). These methods provide satisfactory final fidelities, for the experimental parameters considered, only for total durations exceeding 1 $\mu$s, that is, for a total time which is rather long for performing a single operation in the experimental platform considered. For this reason, we then study how the protocol can be sped up by allowing the presence of a degree of nonadiabaticity in the process which is dynamically counteracted by an additional correcting control Hamiltonian. The latter can be absorbed into the original Hamiltonian by adding quickly oscillations in the initial parameters. This strategy was introduced in Ref. Petiziol et al. (2018) and is built upon the theory of counterdiabatic (CD) or transitionless quantum driving Demirplak and Rice (2003); Berry (2009), but, contrary to this one, it does not require the introduction of new Hamiltonian components. It is realizable instead, in the setup under analysis, by allowing the qubit-bus couplings to oscillate quickly. For brevity, we will refer to this method as effective counterdiabatic (eCD) method. Similar ideas in the use of fast oscillating controls for accelerating adiabatic dynamics were also recently pursued in Refs. Boyers *et al.*Boyers, Pandey, Campbell, Polkovnikov, Sels, and Sushkov ; Bukov et al. (2018), based on Floquet-engineering techniques Bukov et al. (2015).

We show that, by adding eCD corrections to standard adiabatic drivings, high fidelities can be attained much more quickly. In particular, by integrating the eCD scheme with a LAD, timescales of practical interest can be achieved while remaining within state-of-the-art values of the system parameters. At the same time, robustness against systematic errors of the latter protocol is proven by showing that small random imperfections in the control parameters only weakly affect the final fidelity. High fidelities are maintained also in the presence of relaxation and decoherence of the two qubits and photon leakage from the resonator, as shown by numerical solution of a proper Lindblad-type master equation.

The physical system under analysis is introduced in Sec. II, where also the basics of the adiabatic control protocol are presented, and different choices for the adiabatic sweep function are discussed. In Sec. III, the CD method is reviewed before describing the eCD method. The method is then applied to speed up the adiabatic process and the results are discussed in Sec IV, where it is shown that best performance is achieved by combining eCD corrections with a local adiabatic driving. The robustness of these results against imperfections in the control parameters, and against dissipation and decoherence effects is analyzed in Sec. V, followed by the conclusion in Sec. VI.

II The model

II.1 General setup

The general model we will consider is that of two (artificial) atoms whose interaction is mediated by a single bosonic mode acting as quantum bus. The Hamiltonian $\mathcal{H}$ thus includes two anharmonic ladders representing the atoms, described by field operators $b_{k}$ and $b_{k}^{\dagger}$ with $k=1,2$, and a single-mode harmonic oscillator described by field operators $a$ and $a^{\dagger}$ Blais et al. (2007): In units of $\hbar$, it reads

[TABLE]

The parameter $\omega_{k}$ denotes the transition frequency of the two lowest levels of atom $k$, $\omega_{r}$ is the oscillator transition frequency, $\alpha_{k}$ denotes the anharmonicity of the atoms’ energy levels, and $g_{k}$ is the coupling between atom $k$ and the resonator. Performing a two-level approximation for the two atoms and assuming the rotating-wave approximation (RWA), the Hamiltonian of Eq. (1) can be rewritten in Jaynes-Cummings form Mandel and Wolf (1995); Wendin (2017); Gu et al. (2017),

[TABLE]

where $\{\sigma_{x}^{(k)},\sigma_{y}^{(k)},\sigma_{z}^{(k)}\}$ are the Pauli matrices acting on the Hilbert space of qubit $k$ and $\sigma_{\pm}^{(k)}=[\sigma_{x}^{(k)}\pm i\sigma_{y}^{(k)}]/2$ are the corresponding raising and lowering operators. This Hamiltonian can physically describe, for instance, atoms in a cavity Haroche and Raimond (2006), atoms Häffner et al. (2008) and trapped ions Saffman et al. (2010) under external laser fields, or superconducting qubits coupled to transmission line resonators Blais et al. (2007); Majer et al. (2007); Devoret and Schoelkopf (2013); Wendin (2017). The latter is the specific framework which we will take into account in this work. In fact, these platforms have proved to be extremely promising candidates for building quantum technology, especially because of their scalability as solid-state devices and their relatively long coherence times Devoret and Schoelkopf (2013); Clarke and Wilhelm (2008); Schoelkopf and Girvin (2008); Gambetta et al. (2017). They have been a warhorse for the experimental investigation of many complex quantum phenomena, such as basic quantum algorithms DiCarlo et al. (2009) and quantum error correction Reed et al. (2012); DiCarlo et al. (2010), quantum teleportation Steffen et al. (2013), quantum simulation Langford et al. (2017), multi-qubit entanglement Barends et al. (2014); DiCarlo et al. (2010) and Jaynes-Cummings physics Fink et al. (2008).

We will assume tunability both of the internal transition frequency of the qubits, $\omega_{k}=\omega_{k}(t)$, and of the qubit-resonator couplings, $g_{k}=g_{k}(t)$. The former is standard for typical superconducting qubits, such as the transmon Koch et al. (2007). This is a Cooper-pair box charge qubit Clarke and Wilhelm (2008) shunted by an additional high capacity which allows large ratios $E_{J}/E_{C}$ between Josephson energy $E_{J}$ and charging energy $E_{C}$, as compared to a standard Copper-pair box. This makes it rather insensitive with respect to low-frequency charge noise. As a result, it exhibits long coherence times with respect to standard charge, flux and phase qubits Koch et al. (2007). The transmon can be loosely though of as a LC electric circuit with the inductance replaced by a Josephson junction which induces the anharmonicity of the energy levels.

Tunability of the coupling can be achieved by replacing, for instance, the standard transmon with a so-called tunable coupling qubit (TCQ) proposed in Refs. Gambetta et al. (2011); Srinivasan et al. (2011). The TCQ is constituted of three superconducting islands connected by Superconducting Quantum Interference Devices (SQUIDs). As suggested in Ref. Srinivasan et al. (2011), it can be thought of as a system composed of two transmonic qubits capacitively coupled. By changing the energies of the islands, which is done by varying the energy and external magnetic flux applied to the SQUIDs, one can change the dipole moments of the two transmons. These in turn can combine in a parallel or antiparallel way, resulting globally in a dipole or quadrupole moment, respectively. In the parallel case, the TCQ effectively couples strongly to the resonator, while in the antiparallel case it does not. As a result, tuning between the two configurations allows one to tune the TCQ-resonator coupling.

II.2 Specific framework

Starting from the Hamiltonian of Eq. (2), let us restrict the discussion to the two-excitation manifold: The Hamiltonian is then a four-by-four matrix. Let us move to a reference frame rotating at the resonator frequency $\omega_{r}$ and introduce the qubit-resonator detunings $\delta_{k}=\omega_{k}-\omega_{r}$. In the dressed basis $\{\ket{0\uparrow\uparrow},\ket{1\uparrow\downarrow},\ket{1\downarrow\uparrow},\ket{2\downarrow\downarrow}\}$, where the first quantum number indicates the number of field excitations while the others refer to the two qubits bare states, the Hamiltonian explicitly reads

[TABLE]

In the dispersive regime $g_{k}/\delta_{k}\ll 1$, an avoided crossing is formed when the two qubits are brought into resonance far from the resonator’s transition frequency Blais et al. (2004, 2007); Majer et al. (2007). This effect shows the evidence of a bus-mediated coupling between the two qubits, produced by the exchange of a virtual photon. Such a process is well understood by second order perturbation theory in $g_{k}/\delta_{k}$ Blais et al. (2007). As shown in Appendix C, by applying the transformation $U=\exp\left[\sum_{k}\frac{g_{k}}{\delta_{k}}(a\sigma_{+}^{(k)}-a^{\dagger}\sigma_{-}^{(k)})\right]$ to (2), the Hamiltonian $UHU^{\dagger}$ to second order exhibits a qubit interaction term $\propto\sigma_{+}^{(1)}\sigma_{-}^{(2)}+\sigma_{-}^{(1)}\sigma_{+}^{(2)}$ Blais et al. (2004); Majer et al. (2007). The latter is responsible for the formation of the anticrossing, and can produce entanglement between the qubits. For example, it can be exploited to realize a $\sqrt{i\text{SWAP}}$ gate Majer et al. (2007). We denote the width of the anticrossing by $2g_{0}$, with $g_{0}$ being positive.

Here, our purpose is to take advantage of the above described avoided crossing for the adiabatic preparation of an entangled Bell state. Rather than kicking the qubits fast into resonance with each other, and then letting them evolve freely under the time-independent Hamiltonian, the idea is to escort them slowly, so that the system always remains in an instantaneous eigenstate of the evolving Hamiltonian as predicted by the adiabatic theorem. If this is the case, the system, initially assumed to be in the first excited state, which is $\sim\ket{1\downarrow\uparrow}$ if $\omega_{1}>\omega_{2}$, will end up at resonance in the entangled (instantaneous eigen)state $\sim(\ket{1\uparrow\downarrow}+\ket{1\downarrow\uparrow})/\sqrt{2}$. Alternatively, starting from the second excited $\sim\ket{1\uparrow\downarrow}$ would produce the singlet state $(\ket{1\uparrow\downarrow}-\ket{1\downarrow\uparrow})/\sqrt{2}$. Entanglement preparation via adiabatic methods in cQED was also studied in different frameworks and architectures; see, for instance, Refs. DiCarlo et al. (2009); Chen et al. (2014); Martinis and Geller (2014); Xu et al. (2016).

For achieving satisfactory final fidelities, slowly means typically too slow though with respect to reasonable coherence times. For this reason, a shortcut-to-adiabaticity method first introduced in Ref. Petiziol et al. (2018) will be used which tracks to the true adiabatic path closely without requiring the introduction of the new terms in the Hamiltonian (3). It requires, in exchange for that, fast oscillation of the available terms.

The setup is the following. Let the protocol have a total duration $t_{f}$. Let the two qubit-bus couplings have the same value, $g_{1}=g_{2}\equiv g$, and let all quantities be rescaled by $g$. The whole protocol proposed here can be formulated also for the case in which $g_{1}\neq g_{2}$ by choosing a different final configuration of the system. Specifically, one should determine the configuration such that the final instantaneous eigenvector is the desired Bell state, which for $g_{1}\neq g_{2}$ occurs when the qubits are not exactly in resonance with each other. We introduce a sweep function

[TABLE]

which drives the (rescaled) qubit detuning $(\delta_{1}-\delta_{2})/g$ from an initial (positive) value $2f_{0}$ to zero. When we rescale the physical time like $s=t/t_{f}$ and introduce $\tau=gt_{f}$ and $\delta_{g}=[\delta_{1}(t)+\delta_{2}(t)]/2g$, the dimensionless Schrödinger equation for the evolution operator reads

[TABLE]

with

[TABLE]

Using the matrices $\{D_{0},D_{1},C_{1},C_{2}\}$ introduced in Appendix A, this Hamiltonian can be expressed more compactly like

[TABLE]

A sketch of the whole setup with the relevant quantities is represented in Fig. 1. An experimental setup very similar to the one considered here was used in Ref. Quintana et al. (2013) for studying Landau-Zener-Stückelberg effects at the avoided crossing. In the following, the parameters which will be used for all the simulations are taken to be $\omega_{r}/2\pi=8.2$ GHz, $g/2\pi=50$ MHz, and the qubit transition frequencies $\omega_{1}/2\pi$, $\omega_{2}/2\pi$ sweep from $6.01$ and $5.99$ GHz, respectively, to $6$ GHz at resonance producing an avoided crossing of width $2g_{0}\sim 14$ MHz. These values are chosen among typical parameter ranges used for cQED systems Majer et al. (2007); Schoelkopf and Girvin (2008); DiCarlo et al. (2009). Translated into the dimensional quantities appearing in Eq. (5), they give $\delta_{g}=-44$, $f_{0}=0.2$.

II.3 Sweep functions

Here we discuss different choices for the finite-time adiabatic sweep function $f(s)$ of Eq. (4) which brings the two qubits into resonance starting from an initial rescaled gap $(\delta_{2}-\delta_{1})/g=2f_{0}$, so that $f(0)=f_{0}$ and $f(1)=0$. We compare their performance in terms of final fidelity with respect to the total duration of the sweep. Specifically, the possibilities that are taken into account are the following, which are represented in Fig. 2:

1. (1)

Linear (Landau-Zener-like):

[TABLE] 2. (2)

Polynomial, with three derivatives vanishing at the beginning and at the end of the sweep.

[TABLE]

Polynomials with higher and lower numbers of zero derivatives were also considered. The lower number case performed worse, while there was negligible improvement in the higher number case. 3. (3)

Regularized Beta function:

[TABLE]

with

[TABLE]

where $B_{s}(a,b)=\int_{0}^{s}y^{a-1}(1-y)^{b-1}dy$, with $\text{Re}(a)>0$, $\text{Re}(b)>0$, $|s|\leq 1$ is the incomplete Beta function Arfken and Weber (2005). This function was proposed in Ref. Rezakhani et al. (2010) as an adiabatic sweep function and has $k$ derivatives vanishing at the boundary. We will show the results obtained with $k=8$, since larger $k$ did not result to improve the performance in our cases. 4. (4)

Roland-Cerf (RC):

[TABLE]

where $2g_{0}$, introduced in Sec. II.2, is the minimum gap between the first and second excited states in the final configuration – that is for $f(s)=0$. The value of $g_{0}$ can be approximately determined by perturbative arguments, see Appendix C and Ref. Blais et al. (2004). Here the value is extrapolated from numerical diagonalization. More details on the derivation of the functional dependence in Eq. (7) are in Appendix B. 5. (5)

Tangent:

[TABLE]

with $\alpha=\arctan(gf_{0}/g_{0})$ and $g_{0}$ is as in point (4). The derivation of this function in connection to local adiabaticity is given in Appendix B.

Cases (2) and (3) belong to the class of so-called boundary cancellation methods (BCMs) Lidar et al. (2009). The latter predict that, once a duration $t_{f}$ sufficient to satisfy the global adiabatic condition (see Appendix B) is chosen, greater (arbitrary, in principle) precision can be achieved by setting to zero an always larger number of derivatives at the beginning and at the end Rezakhani et al. (2009); Lidar et al. (2009). Cases (4) and (5) represent local adiabatic drivings (LADs) instead, in which the instantaneous velocity of the evolution is adapted such as to satisfy a local adiabatic condition Roland and Cerf (2002); Rezakhani et al. (2009) (see again Appendix B). The RC method was proven to be optimal for a one-dimensional (1D) quantum search problem with respect to a time-optimal adiabatic variational principle in Ref. Rezakhani et al. (2009). A RC function for a Landau-Zener avoided crossing problem was also studied, both theoretically and experimentally, in Refs. Bason et al. (2012); Malossi et al. (2013).

The results produced by all the sweep functions considered here can be compared in Fig. 3, where the final infidelity, defined as

[TABLE]

with $\ket{\psi(t_{f})}$ the final system state and $\ket{e(t_{f})}$ the final first excited state that we wish to prepare, is shown as a function of the total duration of the protocol. Figure 3 is produced by solving numerically the Schrödinger equation for different sweep functions and different total time $t_{f}$. The use of the functions $f_{\text{PL}}$ and $f_{\beta}$ (red squares and blue dashed line, respectively) produce almost equivalent performance, and they prove to dramatically increase the precision with respect to a simple Landau-Zener (LZ) sweep at fixed protocol duration. For example, an infidelity of $10^{-3}$, requiring $10$ $\mu$s if a LZ sweep is used, can be attained in $\sim 2-3$ $\mu$s using $f_{PL}$ or $f_{\beta}$. Further improvement is obtained for shorter times using the RC or tangent sweeps $f_{\text{RC}}$ and $f_{\text{tan}}$ (green-solid and purple-circles, respectively). For instance, an infidelity of $10^{-2}$ can be attained in $\sim 0.5$ $\mu$s , while the same precision would require $\sim 1-2$ $\mu$s and $\sim 2.5$ $\mu$s using BCMs and a LZ sweep, respectively. Therefore, if one can accept timescales around 1 $\mu s$, global adiabatic drivings are the best choice, while LADs are to be used for ulteriorly speeding up the process. One can see that the RC sweep produces fidelities highly fluctuating with respect to the final duration. In principle, this effect could be exploited by selecting a duration close to a local minimum of the infidelity, but the high sensitivity suggests that this would be very hard in practice, and the protocol would still be quite unstable with respect to control imperfections.

III Counterdiabatic control fields

In this section, we discuss the use of shortcuts to adiabaticity Torrontegui et al. (2013) for speeding up the desired adiabatic state transfer and entanglement production. The first case which we treat is that of the CD quantum driving Demirplak and Rice (2003); Berry (2009) theory. This provides a nonadiabatic method which guarantees unit fidelity in arbitrary time, but at the cost of introducing many new time-dependent control components and, for very fast times, strong correcting fields Theisen et al. (2017); Petiziol et al. (2018); del Campo (2013). Implementation of the full protocol turns out to be too demanding in the experimental platform considered here, so we will address the case in which only a partial, dominant correction is taken into account. Since this is in turn hard to realize in practice, we will make use of the recently proposed method of effective counterdiabatic driving Petiziol et al. (2018). This allows one to produce the needed correction terms in an effective manner by introducing fast oscillations in the already available parameters only.

III.1 Exact counterdiabatic corrections

Let us derive the CD correcting field for the setup described in Sec. II: This will be the starting ingredient for determining the final eCD Hamiltonian.

Let the original Hamiltonian of the system be

[TABLE]

The theory of counterdiabatic (CD) fields Demirplak and Rice (2003); Berry (2009) provides a recipe for finding a Hamiltonian $H_{\text{CD}}(t)$ which drives the eigenstates of the original Hamiltonian $H(t)$ with no transitions in a desired time.

Given a matrix $U_{\text{d}}(t)$ which diagonalizes $H(t)$ at time $t$,

[TABLE]

the CD field has the form

[TABLE]

A more convenient expression, in terms of the instantaneous spectral properties of $H(t)$, is Berry (2009)

[TABLE]

Exact CD driving has been investigated experimentally, for instance, in Ref. Bason et al. (2012) and, using superconducting qubits, in Ref. Vepsäläinen et al. (2018); Vepsäläinen *et al.*Vepsäläinen, Danilin, and Paraoanu ; Vepsäläinen et al. (2019).

Let us now focus on the system of Eq. (5). For the rest on this subsection, the reference adiabatic sweep we will work with is $f_{\text{PL}}$ as defined in Eq. (6). The CD field cannot be computed analytically, since the Hamiltonian matrix cannot be generally diagonalized by exact methods. Therefore, inspection of its structure is done numerically. From general results on CD fields Petiziol et al. (2018), one already knows that $H_{\text{CD}}$ will have only imaginary matrix entries, since $H(t)$ is real. Therefore, its components are hard to be realized within the specific experimental setup considered here. New control elements should be introduced, or one should radically change the whole architecture. This is the reason why we will treat the problem via the eCD method in the next section. The (imaginary part of the) matrix elements

[TABLE]

are plotted in Fig. 4. One immediately sees that the largely dominant component is a coupling between the bare states $\ket{1\uparrow\downarrow}$ and $\ket{1\downarrow\uparrow}$, $H_{\text{CD}}^{(2,3)}$. Taking also into account that the corresponding two levels are those with minimal gap, it is easy to recognize that this component gives the by far dominant contribution to nonadiabatic errors. Such a component is enhanced, while the others are further suppressed, if the detuning $\delta_{1,2}$ between qubits and the resonator increases.

In the following, we will thus concentrate on the effect of $H_{\text{CD}}^{(2,3)}$. In particular, this will be compensated for via the eCD method, while the others will be neglected. This corresponds to taking a two-level approximation focusing on the subspace spanned by $\ket{1\uparrow\downarrow}$ and $\ket{1\downarrow\uparrow}$, while neglecting transitions outside this space. However, we still work with the numerical values obtained with the full matrix, rather than using an adiabatically reduced two-level matrix. For the moment, let us define the “partially” correcting CD field $H_{\text{p}}$, defined as

[TABLE]

which equals $H_{\text{CD}}$ with $H_{\text{CD}}^{(i,j)}=0$ for $(i,j)\neq(2,3)$ or $(3,2)$.

The evolution of the system under the action of $H(s)+H_{\text{p}}(s)$ gives at least order-$10^{-5}$ infidelity for $t_{f}>1$ ns, and still improves as the total duration increases. The deviation of the infidelity from zero quantifies the error due to neglecting all nonadiabatic couplings apart from the one connecting first and second instantaneous excited states. The CD method would thus give great results if one was able to implement at least the $H_{\text{CD}}^{(2,3)}$ matrix elements. Still, the total duration affects the strength of the corrective field: The faster the protocol, the stronger $H_{\text{CD}}$ is, since $H_{\text{CD}}$ scales like $1/t_{f}$ Berry (2009). A measure of the magnitude of the matrix elements of $H_{\text{CD}}$ for our parameters can be estimated from Fig. 4: It is the value represented in the figure divided by $t_{f}$, in units of frequency. For example, achieving a fidelity greater than $99.999\%$ in 50 ns would require a $H_{\text{CD}}^{(2,3)}$ component with maximal amplitude around 10 MHz. In the architecture under consideration thought, the necessary strong coupling matrix elements are not available, and for this reason we introduce the eCD method in the next section.

III.2 eCD corrections

The eCD method which will be used now has been introduced in Ref. Petiziol et al. (2018), for realizing a CD Hamiltonian in an approximate manner, by working only with the initially available control Hamiltonians. This kind of problem was also addressed in Ref. Opatrný and Mølmer (2014), using different methodologies. A summary of the general procedure for constructing the correcting Hamiltonian is the following:

1. (1)

Write an ansatz Hamiltonian $H_{\text{e}}$ constituted of the available control Hamiltonians $\{D_{0},D_{1},C_{0},C_{1}\}$ given in Appendix A, controlled by means of control functions $\{d_{0}(t),d_{1}(t),c_{0}(t),c_{1}(t)\}$, respectively,

[TABLE] 2. (2)

Choose periodic control functions, of period $\omega=2\pi/T$, written in the form of a truncated Fourier expansion

[TABLE]

1. (3)

Ask the first terms of the Magnus expansion Blanes et al. (2009) generated by $H_{\text{e}}$ to match those of the desired evolution at the end of each period $T$, by enforcing equality through the amplitudes $\{\mathcal{A}_{k},\mathcal{B}_{k}\}$ in the control functions. 2. (4)

Interpolate solutions in different periods such as to obtain smooth control functions. This is straightforward in two-state problems for eCD schemes at first order in $T$.

With the results of the previous section, let us compute the eCD Hamiltonian. One can verify that the matrix necessary to realize $H_{\text{p}}$, $\sigma_{x}^{(1)}\sigma_{y}^{(2)}-\sigma_{y}^{(1)}\sigma_{x}^{(2)}$, can be obtained as the commutator of the coupling matrices

[TABLE]

This suggests choosing an eCD Hamiltonian of the form

[TABLE]

with control functions $c_{1}(t)$ and $c_{2}(t)$ of the form in Eq. (13). Fast oscillations in $c_{1}$ and $c_{2}$ are then tailored to effectively emulate the dynamics induced by $H_{\text{p}}$ of Eq. (12). We make the simple choice

[TABLE]

The proportionality of $c_{1}(t)$ and $c_{2}(t)$ to $\sqrt{\omega}$ is necessary for having the second term in the Magnus expansion to be of first order in $T$. The exponent $M_{\text{e}}$ of the Magnus expansion for the dynamics $U_{\text{e}}$ induced by $H_{\text{e}}$, computed at the end of one period, is

[TABLE]

while the exponent $M_{\text{p}}$ for the Magnus expansion of $U_{\text{p}}$ induced by $H_{\text{p}}$ is, using Eq. (12),

[TABLE]

Let us remark that $H_{\text{CD}}^{(2,3)}$ is positive for all times if the initial adiabatic sweep function is $f_{\text{PL}}$; see Fig. 4. Moreover, for small $T$, it holds that $H_{\text{CD}}^{(2,3)}(T/2)=H_{\text{CD}}^{(2,3)}(t)+o(T^{3})$. Therefore, a solution of the constraint equation $M_{\text{e}}(T)\simeq M_{\text{p}}(T)$ to order $T$ can be obtained from Eqs. (14) and (15) as

[TABLE]

The complete eCD Hamiltonian then reads

[TABLE]

To clarify the whole procedure we illustrate in Fig. 5 how the correcting Hamiltonians affect the matrix structure of the original Hamiltonian $H(t)$.

IV Results

The results given by the proposed quantum control method are reported in Figs. 6 and 7. In Fig. 6, the eCD field is applied to assist an adiabatic sweep driven by the sweep function $f_{\text{PL}}$ of Eq. (6), which proved in Sec. II.3 to be the best global adiabatic sweep; see Fig. 3. In Fig. 7, the reference sweep is instead a local adiabatic one, $f_{\text{tan}}$ of Eq. (8). We favored this one over $f_{\text{RC}}$ of Eq. (7) due to its stability with respect to changes in the total duration, see again Fig. 3.

Before discussing the results, let us explain how they are obtained and represented. The parameters of the eCD corrections — that is, amplitude and frequency — are chosen such as to satisfy reasonable realizability constraints. The latter regard principally the strength of the time-dependent qubit-resonator coupling. The frequency must also be consistent with the RWA used in deriving the Hamiltonian of Eq. (5); namely, it should not be as high as that of neglected terms and should not be resonant with unwanted transitions. For these reasons, in producing the colored solid lines in Figs. 6 and 7, two ceilings on the amplitude and frequency of the eCD field of Eq. (17) are imposed. For the amplitude, we impose that the maximal allowed value of $\sqrt{\omega}\mathcal{A}(t)$ of Eq. (16) does not surpass a factor $k$ times the initial qubit-bus coupling $g$ — that is, $\sqrt{\omega}\mathcal{A}(t)=kg$. Using Eq. (17), this translates into an upper bound on $\omega$ for a given total duration $t_{f}$:

[TABLE]

Note that this expression is reminiscent of an adiabaticity condition; see Eq. (23) in Appendix B. Together with Eq. (18), also the RWA must be taken into account. For this reason, even if Eq. (18) is satisfied for a given $k$ and $t_{f}$, we impose a second ceiling on the frequency: The main (solid colored) curves are produced by requiring it not to surpass $\omega/2\pi=7$ GHz. This threshold is an arbitrary choice, made by considering that in this way the eCD controls are not resonant with any transition and that terms neglected by the RWA are of the order of $(\omega_{r}+\omega_{1,2})/2\pi\simeq 14$ GHz with our working parameters. Nonetheless, in Figs. 6 and 7, one can easily extrapolate the results for different frequency ceilings by making use of the dotted, dashed, and dot-dashed black lines, which represent the eCD results for fixed $\omega/2\pi$ ranging from 4 to 8 GHz. More details about the RWA will be discussed in Sec. IV.3. Summarizing, the solid colored eCD curves in Fig. 6 and 7 are obtained by setting the eCD frequency $\omega$ to the largest value smaller or equal to $2\pi\times 7$ GHz that satisfies Eq. (18). The color palette then indicates the resulting maximal amplitude of the eCD control fields in units of $g$.

IV.1 Boundary cancellation method and eCD

With the previous caveats in mind, one can see from Figs. 6 and 7 that the eCD method helps to achieve an important speed up of the assisted adiabatic evolutions. Let us first discuss Fig. 6, where the result of adding eCD corrections to a boundary cancellation method (BCM) is reported. The use of an eCD correcting field with maximal strength equal to $g$ [$k=1$ in Eq. (18)], giving an overall qubit-bus coupling oscillating between 0 and 100 MHz, produces fidelities above 99.9% with $\sim 30\%$ speed up with respect to the unassisted sweep. The total timescale in this case still remains above $2\mu$s though, so stronger correcting fields are needed in order to access timescales of practical interest. If one allows the eCD field to reach maximal amplitudes of $2g$ [$k=2$ and 3 curves in Eq. (18)] with frequency around $7$ GHz, the same fidelities can be reached, for instance, starting from $\sim 700-800$ ns.

Let us note that, from Fig. 6, two general behaviors can be observed. Initially, condition (18) on $\omega$ is stricter than the ceiling at $7$ GHz: In this regime, the infidelity goes down rapidly with respect to the total time, even if following a polynomial scaling (linear in the logscale in Fig. 6). When condition (18) starts to give frequencies larger than $7$ GHz, the $7$-GHz ceiling takes over, enforcing a slower, but exponential, decay of the infidelity. One can thus see that being able to reach higher ratios $k$ is mostly relevant for faster timescales.

IV.2 Local adiabatic driving and eCD

Let us now turn to Fig. 7. In this case, the eCD corrections are applied to assist a LAD, namely $f_{\text{tan}}$ of Eq. (8). One can immediately see that fidelities between 99% and 99.99% are achievable in shorter duration with respect to the BCM+eCD case of Fig. 6. For higher fidelities, the BCM+eCD strategy gives better results instead. Using an eCD field with maximal amplitude equaling $g$ [$k=1$ in Eq. (18)], fidelities above $99\%$ are accessible starting from less than $400$ ns, while $\sim 99.9\%$ fidelity is accessible in $\sim 500$ ns. Allowing the maximal amplitude to reach $2g$ and $3g$ gives instead more interesting results: Fidelities above $99\%$ are attained in $\sim 65$ and $\sim 130$ ns, respectively, while fidelities above $99.9\%$ are achievable in $\sim 400$ ns. Let us recall that all these results do not require eCD oscillation frequencies larger than $7$ GHz. As in the BCM+eCD case of the previous subsection, one can see that the $7$ GHz ceiling is dominant for large timescales, while admitting larger $k$ in Eq. (18) determines how fast one can reach the “$7$-GHz regime.” Regimes relative to different maximal frequency ceilings, namely $\omega/2\pi=4,$ 5, 6, and $8$ GHz, can be seen from the black broken lines in Fig. 7: Increasing the maximal allowed frequency in general lowers the whole curve leading to higher fidelities for fixed duration.

This analysis shows that the LAD + eCD method is the most efficient protocol for the state transfer considered here in the timescales of practical interest, with respect to standard finite-time adiabatic driving, BCMs, and BCM + eCD. Comparing Figs. 6 and 7, we see that using the LAD + eCD strategy produces produces 99.9% fidelities more than twenty times faster than a simple LZ sweep.

IV.3 Choosing the system parameters

Let us now discuss how the results presented in the previous subsections depend on the specific experimental parameters used in the simulations (see the last paragraph of Sec. II.2). In general, for fixed detunings $\delta$ between qubits and resonator, raising the maximal values of the qubit-bus coupling $g$ would induce an improvement of all the methods discussed in this work in terms of the total duration needed for a desired infidelity, without changing much the comparison between different methods. This is so because, for the avoided crossing problem studied, raising the minimal gap naturally leads to smaller nonadiabatic errors and thus the global timescales for a given infidelity are shorter, provided $g$ does not become large enough to invalidate the RWA.

Changing the detunings of the qubits from the resonator can have a double effect. In fact, these have a role in (i) realizing the dispersive regime, which works when $g/\delta\ll 1$ and (ii) determining the strength of the resonator-mediated qubit-qubit coupling, which scales like $g^{2}/\delta_{\text{res}}$ where $\delta_{res}$ is the detuning of the two qubits from the resonator when in resonance with each other (see Appendix C). On the one hand, raising $\delta$ ulteriorly consolidates the dispersive regime, making the adiabatic elimination of states $\ket{0\uparrow\uparrow}$ and $\ket{2\downarrow\downarrow}$ a better approximation. On the other hand, it would reduce the mediated coupling and thus amplify the timescales needed for adiabaticity. Therefore, for achieving the best performance with adiabatic methods in our setup one should find the smallest qubit-resonator detunings which give a satisfactory dispersive-regime approximation, while maximizing $g$ for having a larger avoided crossing.

With the specific parameters used in this work, the RWA is a good approximation, in the sense that transitions outside the four-level subspace considered are negligible and conservation of the number of excitations does hold. This is also confirmed by numerical evidence. At the same time, the dispersive regime guarantees that the two levels involved in the avoided crossing exploited for the protocol are sufficiently decoupled from the other two. Nonetheless, counterrotating terms have still frequencies less than ten times greater than those of rotating ones [$|\omega_{k}-\omega_{r}|/2\pi\simeq 2$ GHz, $(\omega_{k}+\omega_{r})/2\pi\simeq 14$ GHz]. As a consequence, considering the full spectrum has an appreciable effect in renormalizing the width of the avoided crossing. More specifically, when $\omega_{r}<\omega_{k}$ the gap tends to reduce, while it becomes larger if $\omega_{r}>\omega_{k}$ as in our case. This is explained by a perturbative calculation in Appendix C. Importantly, taking into account the full spectrum can thus even cause an improvement of all the methods described in this work, thanks to the enlargement of the anticrossing. These features have also been seen in numerical simulations.

Let us now discuss the parameter $gf_{0}$, which quantifies how far the two qubits are from being in resonance at the beginning of the protocol. Smaller values are acceptable unless one starts very close to the avoided crossing, and as the value is decreased all the methods would tend to give the same results, as the whole driving gets closer to a local LZ linear regime. Larger values would make the LAD methods work even better with respect to BCMs. Indeed, by further accelerating in the beginning where nonadiabatic effects are almost negligible, the LADs would need a similar total time for achieving a given fidelity as with the parameters used here.

V Robustness

In this section, we discuss the robustness of the proposed quantum control method against possible sources of errors. In particular, we first study the effect of potential static shifts in the relevant parameters of the protocol and, following, we discuss damping and decoherence effects arising from the interaction of the qubit-resonator system with the environment. This analysis will be done for the case of the LAD + eCD sweep studied in Sec. IV, for $\omega/2\pi=7$ GHz, as represented in Fig. 7, since it has shown to give the best results in terms of fidelity for a given total protocol time.

We start by analyzing numerically how random static shifts in the frequency, amplitude and relative phase of the eCD control fields affect the final fidelity. With reference to the eCD Hamiltonian of Eqs. (16) and (17), this is done by adding a small perturbation to each one of the corresponding parameters, namely

[TABLE]

The perturbations $\bm{\epsilon}=\{\epsilon_{\omega}$, $\epsilon_{\phi}$, $\epsilon_{\mathcal{A}}\}$ are taken to be random variables uniformly distributed in an interval $[-\epsilon_{\text{max}},+\epsilon_{\text{max}}]$. We choose specifically $\epsilon_{\text{max}}=5\%$. The effect of the imperfections on the system evolution can be quantified by studying the average final infidelity $\langle\mathbb{I}_{f}(\bm{\epsilon})\rangle_{\epsilon}$, for a fixed total duration $t_{f}$, where the average $\langle\cdot\rangle_{\bm{\epsilon}}$ is taken over $N_{\bm{\epsilon}}$ realizations of the variables $\bm{\epsilon}$,

[TABLE]

and $\mathbb{I}_{f}(\bm{\epsilon})$ is the final infidelity in the presence of error $\epsilon$. The quantity we monitor is then the relative deviation between the average final infidelity $\langle\mathbb{I}_{f}(\bm{\epsilon})\rangle_{\epsilon}$ and the final infidelity $\mathbb{I}_{f}^{(0)}$ in absence of errors $[\bm{\epsilon}=\bm{0}]$:

[TABLE]

In Figure 8, this quantity is plotted for changing total duration and for errors present in all parameters, together with the maximal and minimal value and the standard deviation

[TABLE]

Both max/min value (asterisks and squares, respectively) and $\sigma$ (error bars) are represented so that one can estimate the worst-case scenario while figuring out what to expect on average.

From the top panel of Fig. 8, one can see that the average value is typically negative: As naively expected, unpredicted variations typically hinder the performance leading to higher infidelities for fixed duration. Nonetheless, the error bars, together with the maximal values, witness that there are some cases in which the fidelity increases in the presence of parameter biases. Considering the procedure described in Sec III.2 for developing the eCD field, this can happen when small errors lead to an accidental improvement of the approximate equality of the first Magnus terms Pet (2019). Moreover, since the eCD dynamics oscillates quickly around the desired one, small biases can bring it closer to the target in a given total time.

Worst-case values present infidelities up to $\sim$6 times greater than the unperturbed case, while the best case can improve up to 100%. The standard deviation is always less than $100\%$. Remarkably, this shows that, although the exact value of the infidelity is quite unstable, the order of magnitude of the final fidelity is rather robust: Errors within $5\%$ in system parameters never produce a jump of the order of magnitude.

To better inspect this feature, we study the average value of the logarithm of the final infidelity,

[TABLE]

This is shown in the lower panel of Fig. 8 for different values of the total duration, compared with the unperturbed values. The error bars indicate again one standard deviation, $\sigma(X)=\sqrt{\langle X^{2}\rangle_{\epsilon}-\langle X\rangle_{\epsilon}^{2}}$ with $X=\log_{10}[\mathbb{I}_{f}(\bm{\epsilon})]$. The average values are very close to the unperturbed curve, and $\sigma$ shows that, as previously observed, the difference between perturbed and unperturbed values is always smaller than one order of magnitude. This gives a good evidence of the robustness of the method: If one is interested in maintaining a certain fidelity threshold of the state transfer, small parameter biases (within $5\%$ here) need not be taken care of. For instance, from Fig. 8, we see that fidelities above $99.9\%$ are achievable in $>400$ ns even in the presence of imperfections.

Let us now turn to a discussion of dissipation and decoherence effects. To this end, we incorporate into the description of the evolution of the system the spontaneous emission and decoherence of the two qubits and photon leakage from the resonator. This is done by describing the evolution of the density matrix $\rho(t)$ of the system through the following master equation in Lindblad form Blais et al. (2007):

[TABLE]

where $\kappa$ is the damping rate for the resonator, $\gamma_{r}^{(k)}$ and $\gamma_{\phi}^{(k)}$ are the relaxation and dephasing rates, respectively, for qubit $k$ and $\mathcal{D}[X]=X\rho X^{\dagger}-\frac{1}{2}\{X^{\dagger}X,\rho\}$ is the Lindblad dissipator. The Hamiltonian $H(t)$ in Eq. (21) is the full Hamiltonian of the controlled system in the RWA, including the eCD correction. Figure 9 shows the final fidelity of the protocol for a representative value of the total duration $t_{f}=100$ ns and for realistic parametric ranges of the rates. This is obtained by numerically solving the master Eq. (21) with three photonic states of the resonator. For simplicity, equal rates for the two qubits are assumed, with $\gamma_{r}^{(1)}=\gamma_{r}^{(2)}\equiv\gamma$ and $\gamma_{\phi}^{(1)}=\gamma_{\phi}^{(2)}=\gamma/2$. For typical values $\kappa/2\pi=5$ kHz and $\gamma/2\pi=5$ kHz, the fidelity is well above $98\%$. Considering rates above $10-15$ kHz makes the fidelity drop below $98\%$, but even for large values it remains above $97\%$. These results show that the proposed control method is rather reliable even in the presence of damping and decoherence effects.

VI Conclusion

An assisted adiabatic passage has been theoretically proposed as a strategy for entangling two qubits in a cQED architecture, with the qubits dispersively coupled to a coplanar waveguide resonator. The method works by slowly driving the qubits into resonance with each other, so that the adiabatic path connects a separable initial eigenstate to a final Bell entangled state (eigenstate of the dispersive flip-flop interaction). The adiabatic passage is accelerated by adding fast oscillations in the qubit-resonator couplings which counteract nonadiabatic transitions, as recently proposed in Ref. Petiziol et al. (2018).

The method achieves both of the desired goals, high fidelities and a short gate time, while requiring only modest extra control resources. Furthermore, it inherits the robustness against unpredicted shifts in control parameters from traditional adiabatic methods, and proves to be stable with respect to damping and decoherence effects. As such, this is expected to be of great practical use in entangling qubits in a standard cQED infrastructure, or for other systems well described by a Jaynes-Cummings Hamiltonian.

It has been shown that best results are obtained by oscillating around a local adiabatic driving, that is, an adiabatic driving whose sweep rate is adapted so to satisfy a local adiabatic condition in each instant of time. Moreover, not only does the protocol work within the rotating-wave approximation, despite requiring fast oscillations, but the presence of the full spectrum can lead to better performances. This happens since, with the parameters used, including counter-rotating terms induces a Bloch-Siegert shift which can enlarge the anticrossing, thus globally reducing the timescales needed for adiabaticity.

Acknowledgements.

We gratefully acknowledge financial support within the Imperial College European Partners programme. This research benefits from the HPC (High Performance Computing) facility of the University of Parma, Italy.

Appendix A Matrices

In this appendix we introduce the following matrices, which are referred to in the main text in Sec. II.2 after Eq. (5), when the Jaynes-Cummings Hamiltonian of Eq. (2) is reduced to the four-by-four matrix of Eq. (3) by considering the subspace of two excitations:

[TABLE]

[TABLE]

Appendix B Local adiabatic driving

In this appendix, we recall the main ideas behind the local adiabatic driving (LAD) strategies, and we derive the sweep functions used in the main text, Eqs. (7) and (8). For relevance to our purpose, we treat directly the case of a two-state avoided-crossing problem, specifically with Hamiltonian ($\hbar=1$)

[TABLE]

with $g>0$, $f(t)\geq 0$ dimensionless, and $\partial_{t}f(t)\leq 0$. This writing of the two-level Hamiltonian is helpful for mapping its parameters to those used in the main text for the Hamiltonian of Eq. (3) and the sweep functions of Sec. II.3. Namely, $g$ and $f(t)$ can be directly mapped to $g$ and $f(t)$ of Eqs. (3) and (4), respectively, while for the parameter $g_{0}$ introduced in Sec. II.2 it holds $g_{0}=|x_{0}|$.

Let us suppose that the system state $\ket{\psi(t)}$ is initially in the instantaneous ground state $\ket{gs(0)}$ of the Hamiltonian $H(0)$, $\ket{\psi(0)}=\ket{gs(0)}$. The adiabatic theorem Messiah (1961) states that, if the adiabatic condition

[TABLE]

is satisfied, with $\ket{e(t)}$ being the excited instantaneous eigenvector at time $t$ and $E_{e}(t),E_{gs}(t)$ being the instantaneous eigenenergies, then the final fidelity is

[TABLE]

The idea which was put forward in Ref. Roland and Cerf (2002) is that the above reasoning can be applied to every arbitrarily small time interval $T$ of the total evolution. Then, one may adapt the speed of the evolution in such a way that the condition (23) is satisfied for the interval of width $T$. In the limit of vanishing $T$, this produces a local form of the adiabatic condition Roland and Cerf (2002); Amin (2009):

[TABLE]

Then, turning the inequality into an equality, one can solve the resulting differential equation to find $t(f)$ such that the condition is satisfied. This can then be inverted in some cases to find the sweep function $f(t)$. Here we discuss the use of two different adiabatic conditions to obtain the corresponding sweep functions $f(t)$.

B.1 Roland-Cerf function

With the Hamiltonian of Eq. (22), the instantaneous gap is

[TABLE]

and the nonadiabatic coupling gives

[TABLE]

The local adiabatic condition (24), assuming $f(t)$ to be monotonically non-increasing $df/dt\leq 0$ with $f(0)=f_{0}>0$, produces the relation

[TABLE]

Once integrated, taking into account that $f(t_{f})=0$, one obtains that the total time is

[TABLE]

The sweep function, introducing the rescaled time $s=t/t_{f}$, is

[TABLE]

B.2 Quantum adiabatic brachistochrone

In Ref. Rezakhani et al. (2009), it is pointed out that more accurate adiabatic conditions typically involve powers of the norm of the derivative of the Hamiltonian, rather than single matrix elements. For this reason and to be able to formulate the problem in a geometric context, Ref. Rezakhani et al. (2009) proposes the condition

[TABLE]

where $\lVert A\lVert=\sqrt{\text{tr}[AA^{\dagger}]}$ is the Hilbert-Schmidt norm. Acting as above to find $f(t)$, the latter choice, for $H(t)$ of Eq. (22), gives the differential equation

[TABLE]

The solution is

[TABLE]

for $s=t/t_{f}$ and $\alpha=\arctan(gf_{0}/x_{0})$. The total time is

[TABLE]

Let us recall that $f_{0}>0$ and $g>0$ so $t_{f}$ is actually positive for all values of $x_{0}$.

Appendix C Anticrossing formation

The purpose of this appendix is to discuss the formation of the anticrossing which is exploited for the adiabatic methods described in the main text. The qubit-qubit coupling mediated by the resonator, responsible for this phenomenon, is derived perturbatively in the RWA first, following Ref. Blais et al. (2004), and without the RWA afterward. Let us consider the full Rabi Hamiltonian

[TABLE]

Introducing the unitary

[TABLE]

with $\delta_{k}=\omega_{k}-\omega_{r}$, the Hamiltonian $U_{1}H_{\text{R}}U_{1}^{\dagger}$, to second order in $g/\delta_{k}$, becomes

[TABLE]

The first line in Eq. (26) would be the result in the RWA, while the second line constitutes the byproduct of keeping counter-rotating terms. The last term in the first line describes the production of the anticrossing, which is exploited in the adiabatic methods described in the main text. When the qubits are in resonance with each other, $\delta_{1}=\delta_{2}\equiv\delta$, we thus find a coupling $\frac{g^{2}}{\delta}$ (to second order), which for our parameters gives the value $\sim-7.1$ MHz, i.e., a minimal gap of $\sim 14.2$ MHz. In order to understand the contribution to the anticrossing of the counterrotating terms, let us make a second basis transformation determined by the unitary matrix

[TABLE]

where $\Delta_{k}=\omega_{k}+\omega_{r}$. This is constructed so to eliminate the remaining (counterrotating) $\propto g$ term in the second line of Eq. (26). Again to second order, one has:

[TABLE]

From the third term in the first line, we now see that the coupling is renormalized, due to the inclusion of counter-rotating terms, to $g^{2}\left(\frac{1}{\delta}-\frac{1}{\Delta}\right)$ (for $\Delta_{1}=\Delta_{2}\equiv\Delta$). Depending on whether the resonator has a higher or lower transition frequency with respect to the two qubits, the new shift can enhance or reduce the strength (absolute value) of the coupling. For our protocol, the first situation is the ideal one: with our parameters, the renormalized coupling is $\sim-8.2$ MHz. Keeping the same $\omega_{k}$ and $|\delta_{k}|$ for the second situation would give $\omega_{r}=3.8$ GHz and the coupling would reduce in strength to $\sim-5.5$ MHz.