Operation and intrinsic error budget of a two-qubit cross-resonance gate

TL;DR

This paper provides a comprehensive analysis of the operation and intrinsic error sources of the two-qubit cross-resonance gate in superconducting qubits, offering methods to optimize gate duration and fidelity.

Contribution

It introduces a semi-analytical approach for accurate gate duration estimation and decomposes the intrinsic error budget, highlighting dominant error mechanisms at different drive amplitudes.

Findings

Semi-analytical method accurately predicts CNOT-equivalent gate duration.

Intrinsic fidelity can be less than 10^{-3} with simple pulse shapes.

Error sources vary with drive amplitude, from rotation imperfections to leakage.

Abstract

We analyze analytically, semi-analytically, and numerically the operation of Cross-Resonance (CR) gate for superconducting qubits (transmons). We find that a relatively simple semi-analytical method gives accurate results for the CNOT-equivalent gate duration and compensating single-qubit rotations. It also allows us to minimize the CNOT gate duration over the amplitude of the applied microwave drive and find dependence on the detuning between the qubits. However, full numerical simulations are needed to calculate intrinsic fidelity of the CR gate. We decompose numerical infidelity into contributions from various physical mechanisms, thus finding the intrinsic error budget. In particular, at small drive amplitudes the CR gate fidelity is limited by imperfections of the target-qubit rotations, while at large amplitudes it is limited by leakage. The gate duration and fidelity are analyzed numerically as functions of the detuning between qubits, their coupling, drive frequency, relative duration of pulse ramps, and microwave crosstalk. The effect of the echo sequence is also analyzed numerically. Our results show that the CR gate can provide intrinsic infidelity of less than $10^{-3}$ when a simple pulse shape is used.

Equations

$$H=H_{\mathrm{qb}}+H_g+H_{\epsilon },$$

$$H_{\mathrm{qb}}={\sum }_{n,m}(E_n^{(\mathrm{c})}+E_m^{(\mathrm{t})})|n,m⟩⟨n,m|,$$

$$E_n^{(\mathrm{c})}=E_n^{(\mathrm{c},\mathrm{lf})}-n{\omega }_{\mathrm{d}},E_m^{(\mathrm{t})}=E_m^{(\mathrm{t},\mathrm{lf})}-m{\omega }_{\mathrm{d}},$$

$$E_n^{(\mathrm{c})}=n(\Delta +\delta )-\frac{n(n-1)}{2}{\eta }_{\mathrm{c}},$$

$$E_m^{(\mathrm{t})}=m\delta -\frac{m(m-1)}{2}{\eta }_{\mathrm{t}},$$

$$\Delta \equiv {\omega }_{\mathrm{c}}-{\omega }_{\mathrm{t}},\delta \equiv {\omega }_{\mathrm{t}}-{\omega }_{\mathrm{d}}\approx 0,$$

$$H_g={\sum }_{n,m}g\sqrt{nm}|n,m-1⟩⟨n-1,m|+h.c.,$$

$$H_{\epsilon }={\sum }_{n,m}\epsilon (t)\sqrt{n}|n,m⟩⟨n-1,m|+h.c.,$$

$${\omega }_{\mathrm{t}}^{c0}=E_{\overline{|0,1⟩}}^{(\mathrm{lf})}-E_{\overline{|0,0⟩}}^{(\mathrm{lf})},{\omega }_{\mathrm{t}}^{c1}=E_{\overline{|1,1⟩}}^{(\mathrm{lf})}-E_{\overline{|1,0⟩}}^{(\mathrm{lf})},$$

$${\omega }_{zz}\equiv {\omega }_{\mathrm{t}}^{c1}-{\omega }_{\mathrm{t}}^{c0}=E_{\overline{|11⟩}}+E_{\overline{|00⟩}}-E_{\overline{|01⟩}}-E_{\overline{|10⟩}},$$

$${\omega }_{zz}\approx \frac{2g^2}{\Delta +{\eta }_{\mathrm{t}}}-\frac{2g^2}{\Delta -{\eta }_{\mathrm{c}}}.$$

$$\displaystyle H_{\rm CR}^{\rm ideal}-(H_{\rm qb}+H_{g})=\big{(}\tilde{\varepsilon}_{0}\overline{|0,1\rangle}\,\,\overline{\langle 0,0|}+\tilde{\varepsilon}_{1}\overline{|1,1\rangle}\,\,\overline{\langle 1,0|}$$

$$\displaystyle\hskip 96.73918pt+\tilde{\varepsilon}_{2}\overline{|2,1\rangle}\,\,\overline{\langle 2,0|}+...\big{)}+{\rm h.c.},$$

$$\frac{{\tilde{\epsilon }}_0-{\tilde{\epsilon }}_1}{2}{Z}_{\mathrm{c}}{X}_{\mathrm{t}}+\frac{{\tilde{\epsilon }}_0+{\tilde{\epsilon }}_1}{2}{I}_{\mathrm{c}}{X}_{\mathrm{t}},$$

$${\int }_0^{{\tau }_{\mathrm{p}}}[2{\tilde{\epsilon }}_1(t)-2{\tilde{\epsilon }}_0(t)]dt=\pi (\mathrm{mod}2\pi ),$$

$${\tilde{\epsilon }}_n=\overline{⟨n,1|}{H}_{\epsilon }\overline{|n,0⟩}.$$

$${\tilde{\epsilon }}_0=-\frac{g}{\Delta }\epsilon .$$

$${\tilde{\epsilon }}_1=\frac{g}{\Delta }\epsilon -\frac{\sqrt{2}g}{\Delta -{\eta }_{\mathrm{c}}}\sqrt{2}\epsilon =-\frac{g}{\Delta }\frac{\Delta +{\eta }_{\mathrm{c}}}{\Delta -{\eta }_{\mathrm{c}}}\epsilon .$$

$${\tilde{\epsilon }}_2=\frac{\sqrt{2}g}{\Delta -{\eta }_{\mathrm{c}}}\sqrt{2}\epsilon -\frac{\sqrt{3}g}{\Delta -2{\eta }_{\mathrm{c}}}\sqrt{3}\epsilon$$

$$=-\frac{g(\Delta +{\eta }_{\mathrm{c}})}{(\Delta -{\eta }_{\mathrm{c}})(\Delta -2{\eta }_{\mathrm{c}})}\epsilon ,$$

$${\tilde{\epsilon }}_n=\frac{ng\epsilon }{\Delta -(n-1){\eta }_{\mathrm{c}}}-\frac{(n+1)g\epsilon }{\Delta -n{\eta }_{\mathrm{c}}}$$

$$=-\frac{g(\Delta +{\eta }_{\mathrm{c}})}{[\Delta -(n-1){\eta }_{\mathrm{c}}](\Delta -n{\eta }_{\mathrm{c}})}\epsilon .$$

$${\tilde{\epsilon }}_1-{\tilde{\epsilon }}_0=\frac{2g{\eta }_{\mathrm{c}}}{\Delta ({\eta }_{\mathrm{c}}-\Delta )}\epsilon ,$$

$${\tilde{\epsilon }}_2-{\tilde{\epsilon }}_0=\frac{2g{\eta }_{\mathrm{c}}({\eta }_{\mathrm{c}}-2\Delta )}{\Delta ({\eta }_{\mathrm{c}}-\Delta )(2{\eta }_{\mathrm{c}}-\Delta )}\epsilon .$$

$${\tilde{\epsilon }}_0=\epsilon \frac{-1}{\Delta }g,$$

$${\tilde{\epsilon }}_1=\sqrt{2}\epsilon \frac{-1}{\Delta -{\eta }_{\mathrm{c}}}\sqrt{2}g+g\frac{-1}{-\Delta }\epsilon ,$$

$${\tilde{\epsilon }}_n=-\frac{\sqrt{n+1}\epsilon \sqrt{n+1}g}{\Delta -n{\eta }_{\mathrm{c}}}+\frac{\sqrt{n}g\sqrt{n}\epsilon }{\Delta -(n-1){\eta }_{\mathrm{c}}},$$

$${\tilde{\epsilon }}_n=\overline{{}_{\epsilon }⟨n,1|}{H}_g\overline{|n,0{⟩}_{\epsilon }},$$

$$\overline{|n,0{⟩}_{\epsilon }}=\overline{|n{⟩}_{\epsilon }}|0{⟩}_{\mathrm{t}},\overline{|n,1{⟩}_{\epsilon }}=\overline{|n{⟩}_{\epsilon }}|1{⟩}_{\mathrm{t}},$$

$$H_{\mathrm{qb}+\epsilon }^{(\mathrm{c})}\overline{|n{⟩}_{\epsilon }}=E_{\overline{|n{⟩}_{\epsilon }}}\overline{|n{⟩}_{\epsilon }}$$

$$H_{\mathrm{qb}+\epsilon }^{(\mathrm{c})}=\sum _n{E}_n^{(\mathrm{c})}|n⟩⟨n|+\sqrt{n}(\epsilon |n⟩⟨n-1|+{\epsilon }^{\ast }|n-1⟩⟨n|).$$

$$\overline{|n{⟩}_{\epsilon }}={\sum }_k{c}_k^{(n)}|k⟩,$$

$${\tilde{\epsilon }}_n=g{\sum }_k\sqrt{k}{c}_k^{(n)}{\left(c_{k-1}^{(n)}\right)}^{\ast }.$$

$$\overline{|0{⟩}_{\epsilon }}=\frac{|0⟩+\alpha |1⟩+\beta |2⟩+\gamma |3⟩+...}{\mathcal{N}},$$

Full text

I Introduction

Two decades have passed since the first superconducting qubit was created Nakamura1999 , and today superconducting quantum computing is a well-developed field with various types and uses of qubits Neill2018 ; Kandala2017 ; Chou2018 ; Hong2019 ; Hacohen-Gourgy2016 ; Satzinger2018 ; King2018 ; Kurpiers2018 ; Minev2018 ; Masuyama2018 . Currently, the most popular type of superconducting qubits is the transmon Koch2007 (including its Xmon and gmon modifications Barends2013 ; Chen2014 ), though other types of qubits (e.g., Chiorescu2003 ; Yan2016 ; Lin2018 ) are also of interest. Besides sufficiently good coherence of the qubits, quantum computing applications need high-fidelity gates forming a universal set Nielsen-Chuang-book . While single-qubit gates are already considered to be simple and accurate, current fidelity of two-qubit gates also exceeds 99% Barends2014 ; Sheldon2016 ; Hong2019 .

One of the high-fidelity two-qubit gates used for superconducting qubits is the Cross-Resonance (CR) gate Rigetti2010 ; Chow2011 . In this gate, two frequency-detuned qubits have a fixed coupling (usually via a resonator) and one of them (called control qubit) is driven by a microwave with frequency of the other one (target qubit). This induces Rabi oscillations of the target qubit, whose frequency depends on the state ($|0\rangle$ or $|1\rangle$) of the control qubit, thus entangling the two qubits and providing a natural way to realize CNOT operation. Since the CR gate uses only microwave control, it permits to use single-junction transmons, thus avoiding sensitivity to flux noise. However, the drawback is a relatively long gate duration compared with the gates based on tune-detune operation Barends2014 .

The idea of the CR gate was proposed in 2006 in Ref. Paraoanu2006 and then experimentally implemented for flux qubits in 2010 in Ref. deGroot2010 under the name of Selective Darkening (the difference compared with a simple CR gate is an additional active cancellation pulse applied to the target qubit). The CR terminology was introduced in 2010 in the theoretical paper Rigetti2010 and the first experiment under this name was realized with capacitively-shunted flux qubits in 2011 in Ref. Chow2011 with the fidelity of 81%. In 2012 the CR gate was applied to transmons Chow2012 , with resulting fidelity of 95%. Since that time the CR gate was used in numerous experiments by several groups (e.g., Corcoles2013 ; Corcoles2015 ; Takita2016 ; Sheldon2016 ; Kandala2017 ; Ware2015 ; Naik2019 ), with gradual increase of maximum fidelity. An important improvement of the CR operation was achieved by using the echo sequence Corcoles2013 ; Takita2016 , which not only increased the fidelity but also allowed protocols that avoid compensating one-qubit rotations in implementing the CNOT gate. The CR gate with duration of 160 ns and fidelity of 99.1% reported in Ref. Sheldon2016 was achieved by using both the echo sequence and active cancellation pulses applied to the target qubit.

In spite of extensive experimental use of the CR gate, its theoretical analysis has been rather limited. Besides the initial papers Paraoanu2006 ; Rigetti2010 outlying the main idea, the CR gate was analyzed in Ref. deGroot2012 with an account of the next level, briefly mentioned in Ref. Gambetta2013conf , and analyzed in detail in the recent paper Magesan2018 . There were also numerical studies Willsch2017 ; Kirchhoff2018 and related papers Economou2015 ; Allen2017 .

In this paper, we analyze the operation of the basic CR gate for transmons (using a simple pulse shape without the echo sequence, which is considered only in the Appendix) at three levels of complexity and accuracy: analytical, semi-analytical, and numerical. Some of the goals of our analysis are similar to those of Ref. Magesan2018 ; however, the approach is very different. After discussing the ideal theory of CR operation for transmons (using the Duffing oscillator model), we develop the next-order approximation somewhat similar to that of Ref. Magesan2018 (which still does not work well, as follows from comparison with full numerics), and then develop the semi-analytical approach, based on the numerical solution of a simple one-qubit time-independent Schrödinger equation.

The semi-analytical approach gives very accurate results (compared with full numerics) for the CNOT-equivalent gate duration and compensating single-qubit rotations. In particular, it can be used to find the shortest CNOT gate duration, corresponding optimal drive amplitude, and their dependence on detuning between the qubits. However, the semi-analytical approach cannot be used for finding intrinsic fidelity of the gate (neglecting decoherence), for which we use full numerical simulation. Our numerical simulation includes $7\times 5$ levels in the qubits (we replace qubit coupling via resonator with an equivalent direct coupling) and is based on Magnus expansion Magnus1954 for the evolution matrices. We use a simple pulse shape with cosine-shaped ramps and a flat middle part.

After calculating the gate infidelity, we numerically decompose it into the contributions from various physical mechanisms, thus finding intrinsic error budget. We show that at small drive amplitudes, the error is dominated by imperfections of the unitary operation within the computational subspace. In contrast, at large drive amplitudes (which correspond to reasonably short gate durations) the infidelity is dominated by leakage. Particular leakage channels depend on detuning between the qubits. In this regime, the analytical estimate for the leakage probability agrees reasonably well with the results for the gate infidelity.

Using numerical results for the CNOT gate duration and fidelity, we analyze their dependence on various parameters, including detuning between qubits, their coupling, drive frequency, smoothness (relative duration of the pulse ramps), and microwave crosstalk. In the Appendix, we also analyze the effect of the echo sequence. Our results show that the CR gate can provide intrinsic infidelity of about $10^{-3}$ (and even less, comparable to $10^{-4}$) with a simple pulse shape.

The paper is organized as follows. In Sec. II we discuss the system and its Hamiltonian. In Sec. III we first consider the ideal operation of the CR gate, then derive the next-order analytics, and then develop the semi-analytical approach. The numerical method is discussed in Sec. IV. Numerical results for the CNOT-equivalent gate duration and compensating single-qubit rotations are discussed in Sec. V. Then in Sec. VI we analyze the error budget for the CNOT-gate intrinsic infidelity. In Sec. VII, we discuss the dependence of CNOT duration and infidelity on parameters. Finally, we conclude in Sec. VIII. In the Appendix, we analyze the echo-CR gate operation.

II System and Hamiltonian

In the CR gate, the control and target qubits (with frequencies $\omega_{\rm c}$ and $\omega_{\rm t}$, respectively) are usually detuned by 50–300 MHz and are permanently coupled via a resonator. However, in this paper, for simplicity we will consider a direct qubit-qubit coupling $g$ (Fig. 1) since the usual analysis of the CR gate Gambetta2013conf ; Magesan2018 also reduces the coupling via a resonator to an effective direct coupling. For the CR operation, the control qubit is rf-driven at the frequency of the target qubit, $\omega_{\rm d}\approx\omega_{\rm t}$. This produces an effective drive ($x$-rotation) of the target qubit, with the strength depending on the state of the control qubit. Such a process can be naturally used to realize the CNOT gate by calibrating the target-qubit rotation angle difference (between rotations for the control-qubit states $|0\rangle$ and $|1\rangle$) to be equal to $\pi$ and somehow compensating the target-qubit rotation for the control-qubit state $|0\rangle$. This compensation can be done, for example, by using the echo sequence Corcoles2013 ; Takita2016 ; Sheldon2016 or active cancellation deGroot2010 ; Sheldon2016 ; however, in this paper we will assume that the compensation is done afterward Rigetti2010 ; Chow2011 by applying single-qubit rotations (the echo sequence is considered only in the Appendix). We intentionally consider the simplest case in order to focus on developing a good understanding of the basic operation of the CR gate.

The operating principle of the CR gate can be understood classically, by replacing qubits with classical oscillators (Fig. 2). Since the drive is off-resonance with the control oscillator ($\omega_{\rm d}\not\approx\omega_{\rm c}$), it will produce very small forced oscillations at the drive frequency $\omega_{\rm d}$. However, since the target oscillator is on-resonance with this frequency ($\omega_{\rm d}\approx\omega_{\rm t}$), it will still get excited via the coupling $g$ with the control oscillator. Note that if the control oscillator is linear, then its own state (its oscillation with frequency $\omega_{\rm c}$) does not matter because of linearity. However, if the control oscillator is nonlinear, then its effective frequency depends on its own state (i.e., amplitude of $\omega_{\rm c}$-oscillations); therefore, the amplitude of the small forced oscillations of the control oscillator and consequently the excitation rate of the target oscillator will depend on the control-oscillator state. This simple classical picture explains the basic physical mechanism of the CR gate operation for transmons, which are slightly nonlinear oscillators. It also explains why the CR gate speed depends on nonlinearity of the control qubit and practically does not depend on the target-qubit nonlinearity.

For quantum analysis of the CR gate (Fig. 1), let us start with the rotating-frame Hamiltonian (the rotating frame is based on the drive frequency $\omega_{\rm d}$)

[TABLE]

where $H_{\rm qb}$ describes two uncoupled transmon qubits, $H_{g}$ describes their coupling, and $H_{\varepsilon}$ describes the microwave drive on the control qubit.

The uncoupled-qubit part can be written as

[TABLE]

where in the notation $|n,m\rangle$ the control-qubit state is at the left ($n=0,1,2,...$) and the target-qubit state is at the right ($m=0,1,2,...$), the control-qubit energies $E_{n}^{\rm(c)}$ in the rotating frame are related to the laboratory-frame energies $E_{n}^{\rm(c,\,lf)}$ via the drive frequency $\omega_{\rm d}$, and there is a similar relation for the target-qubit energies $E_{m}^{\rm(t)}$. We set $E_{0}^{\rm(c)}=E_{0}^{\rm(t)}=0$. For the energies $E_{n}^{\rm(c)}$ and $E_{m}^{\rm(t)}$, in this paper we use the Duffing (Kerr) oscillator approximation,

[TABLE]

where $\Delta$ is the detuning between the qubits, while $\eta_{\rm c}$ and $\eta_{\rm t}$ are anharmonicities of the control and target qubits, respectively (for transmons $\eta_{\rm c}>0$ and $\eta_{\rm t}>0$). A small mismatch $\delta$ between the drive frequency $\omega_{\rm d}$ and the bare frequency $\omega_{\rm t}$ of the target qubit can be used, e.g., to make the drive exactly resonant with the hybridized target qubit for the control-qubit states $|0\rangle$ or $|1\rangle$ (or in between). Note that $\Delta+\delta=\omega_{\rm c}-\omega_{\rm d}$. Instead of the approximation (4)–(6), it is possible to use numerical results for the transmon energies or at least the improved approximation Khezri2018 ; Sank2016 . However, we prefer the simple approximation for easier comparison with the previous theoretical analyses of the CR gate.

The qubit-qubit coupling Hamiltonian $H_{g}$ in general couples all pairs of the bare states $|n,m\rangle$ and $|n^{\prime},m^{\prime}\rangle$. However, in this paper, we use the simplest (traditional) approximation for transmons by keeping only the excitation-preserving terms, i.e., applying the Rotating Wave Approximation (RWA), and using the matrix elements for linear oscillators,

[TABLE]

additionally assuming (without loss of generality) that the coupling constant $g$ is real. Similarly, we use the RWA linear-oscillator matrix elements for the drive Hamiltonian (in the rotating frame),

[TABLE]

where the complex amplitude $\varepsilon$ of the drive depends on time, so that $\varepsilon(t)$ is the pulse shape of the CR gate, with $\varepsilon(t)=0$ before and after the gate. Instead of Hamiltonians (7) and (8), it is possible to use improved perturbative Hamiltonians Khezri2018 ; Sank2016 or numerical matrix elements for transmons, but in this paper, we use the simple traditional approximation. Here we do not consider the microwave crosstalk Chow2011 ; Chow2012 ; Sheldon2016 ; Magesan2018 ; however, it will be added in Sec. VII.

It is convenient to draw a diagram (Fig. 3) of bare levels $|n,m\rangle$, in which the left ladder of levels corresponds to the target-qubit state $|0\rangle$ ($m=0$), the next ladder corresponds to the target-qubit state $|1\rangle$, then $|2\rangle$, and so on. Note that for $\delta=0$, the left two ladders are at exactly equal energies. In Fig. 3 the coupling $H_{g}$ is represented by slanted blue arrows and the drive $H_{\varepsilon}$ corresponds to vertical orange arrows. For clarity, in Fig. 3 we show the case $\Delta>2\eta_{\rm c}$, while in experiments usually $0<\Delta<\eta_{\rm c}$. In such a case, all ladders turn down after the states $|1,m\rangle$ and the diagram becomes visually complicated, so for gaining intuition it is easier to use the case of Fig. 3.

Besides the bare states $|n,m\rangle$, we will also use the eigenstates of the Hamiltonian $H_{\rm qb}+H_{g}$ (without the drive), which we denote by an overline: $\overline{|n,m\rangle}$. The coupling $H_{g}$ affects the qubit frequencies, so instead of the bare frequency $\omega_{\rm t}$ of the target qubit, we have two eigenfrequencies: $\omega_{\rm t}^{\rm c0}$ and $\omega_{\rm t}^{\rm c1}$, depending on the control-qubit state ($|0\rangle$ and $|1\rangle$, respectively). They can be calculated as

[TABLE]

where $E_{\overline{|n,m\rangle}}^{\rm(lf)}$ is the laboratory-frame eigenenergy of the state $\overline{|n,m\rangle}$. We will call “$zz$-coupling” the difference between these frequencies,

[TABLE]

where this combination of eigenenergies is the same in the laboratory and rotating frames. The $zz$-coupling is mainly due to the repulsion of the energy level $|11\rangle$ from the levels $|02\rangle$ and $|20\rangle$, which gives the approximate value

[TABLE]

From Eq. (10) we see that the $zz$-coupling can be also defined as $\omega_{zz}=\omega_{\rm c}^{\rm t1}-\omega_{\rm c}^{\rm t0}$, where $\omega_{\rm c}^{\rm t0}$ and $\omega_{\rm c}^{\rm t1}$ are the eigenfrequencies of the control qubit for the target-qubit states $|0\rangle$ and $|1\rangle$, respectively. Nonzero $\omega_{zz}$ will be important for numerical results; however, it will be neglected for analytical and semi-analytical results in the next section; in particular, we will not distinguish between $\omega_{\rm t}^{\rm c0}$, $\omega_{\rm t}^{\rm c1}$, and $\omega_{\rm t}$.

III Analytical and semi-analytical analysis

III.1 Ideal CR gate operation

There is no drive, $\varepsilon=0$, before and after the CR gate operation. Therefore, the initial and final two-qubit states should be considered in the eigenbasis $\overline{|n,m\rangle}$ of the Hamiltonian $H_{\rm qb}+H_{g}$. The drive Hamiltonian $H_{\varepsilon}$ couples these eigenstates, providing an evolution used in the CR gate.

As follows from Fig. 3, in the rotating frame based on the drive frequency $\omega_{\rm d}$, there is a near-resonance condition between states $|n,0\rangle$ and $|n,1\rangle$, which leads to a near-resonance between eigenstates $\overline{|n,0\rangle}$ and $\overline{|n,1\rangle}$, while other pairs of states are off-resonance. Therefore, as long as the perturbation produced by $H_{\varepsilon}$ is small enough, it effectively couples only states $\overline{|n,0\rangle}$ and $\overline{|n,1\rangle}$, and for the ideal effective Hamiltonian $H_{\rm CR}^{\rm ideal}$ of the CR gate we can write

[TABLE]

where $\tilde{\varepsilon}_{0}$ is the amplitude of the effective drive on the target qubit when the control-qubit state is $|0\rangle$, $\tilde{\varepsilon}_{1}$ is the effective drive amplitude for the control-qubit state $|1\rangle$, etc. (for small $g$ there is almost no difference between the effective drive in the bare basis or eigenbasis). The effective drive amplitudes $\tilde{\varepsilon}_{n}$ depend on the actual drive amplitude $\varepsilon$ (in the linear approximation being proportional to $\varepsilon$).

Note that if we are interested only in the states $|0\rangle$ and $|1\rangle$ of the control qubit, then in the terminology of Refs. Chow2011 ; Chow2012 ; Magesan2018 ; Ware2015 the effective Hamiltonian (12) can be written as

[TABLE]

where the Pauli operators $Z_{\rm c}$ and $I_{\rm c}$ act on the control qubit and the operator $X_{\rm t}$ acts on the target qubit (here we assume real $\tilde{\varepsilon}_{0}$ and $\tilde{\varepsilon}_{1}$; otherwise, we also need $Y_{\rm t}$). The CNOT gate can be realized with this effective interaction by applying the drive pulse with duration $\tau_{\rm p}$, satisfying the condition

[TABLE]

complemented with two one-qubit rotations. The additional $x$-rotation of the target qubit by the angle $-\int_{0}^{\tau_{\rm p}}2\tilde{\varepsilon}_{0}(t)\,dt$ compensates the target-qubit rotation for the control-qubit state $|0\rangle$, also providing $x$-rotation by angle $\pi$ for the control-qubit state $|1\rangle$. Besides the $x$-rotation of the target qubit, the control qubit should be $z$-rotated by the angle $\pi/2$ (relative to the rotating frame of the control qubit). This is needed because the $x$-rotation of the target qubit by angle $\pi$ produces the operation $-iX$ instead of the desired (for CNOT) operation $X$, thus requiring additional phase factor $i$ for the control-qubit state $|1\rangle$ (the same factor exists in a one-qubit $X$ gate, but it is not important since it is an overall phase, in contrast to the phase difference in a controlled two-qubit operation). Note that the factors of 2 in Eq. (13) are needed because the Rabi frequency is twice larger than the drive matrix element in the Hamiltonian.

The effective drive amplitudes $\tilde{\varepsilon}_{n}$ in Eq. (12) can be easily found (in the ideal lowest-order case) by comparing Eqs. (1) and (12), which gives

[TABLE]

To calculate $\tilde{\varepsilon}_{n}$ to the lowest order, let us assume $\delta=0$, i.e., the drive resonant with the bare target qubit (the difference between bare and eigenfrequencies is not important for these approximate calculations). Then using $\overline{|0,0\rangle}=|0,0\rangle$ (see Fig. 3) and the first-order approximation $\overline{|0,1\rangle}=|0,1\rangle-(g/\Delta)\,|1,0\rangle$ (normalization correction is of the second order), we find the linear approximation

[TABLE]

Similarly, using the first-order approximations $\overline{|1,0\rangle}=|1,0\rangle+(g/\Delta)\,|0,1\rangle$ and $\overline{|1,1\rangle}=|1,1\rangle-[\sqrt{2}\,g/(\Delta-\eta_{\rm c})]\,|2,0\rangle+[\sqrt{3}\,g/(\Delta+\eta_{\rm t})]\,|0,2\rangle$ (see Fig. 3), we obtain approximation

[TABLE]

Also similarly, using approximations $\overline{|2,0\rangle}=|2,0\rangle+[\sqrt{2}\,g/(\Delta-\eta_{\rm c})]\,|1,1\rangle$ and $\overline{|2,1\rangle}=|2,1\rangle-[\sqrt{3}\,g/(\Delta-2\eta_{\rm c})]\,|3,0\rangle+[2g/(\Delta-\eta_{\rm c}+\eta_{\rm t})]\,|1,2\rangle$, we find

[TABLE]

and for an arbitrary control-qubit state $|n\rangle$, within the model (1)–(8) we obtain

[TABLE]

Since the target-qubit rotation for the control-qubit state $|0\rangle$ is usually compensated, most important are the differences of effective drive amplitudes from $\tilde{\varepsilon}_{0}$, e.g.,

[TABLE]

Note that these formulas depend on anharmonicity $\eta_{\rm c}$ of the control qubit but do not depend on the target-qubit anharmonicity $\eta_{\rm t}$. Also, for $\eta_{\rm c}=0$ we have $\tilde{\varepsilon}_{n}=\tilde{\varepsilon}_{0}=-(g/\Delta)\varepsilon$. These properties are in agreement with the classical description of the CR gate operation discussed in Sec. II.

Language of virtual-state transitions

Instead of using Eq. (14), we can find the effective drive amplitudes $\tilde{\varepsilon}_{n}$ (still to the first order) using the ideology of transitions via a virtual state. As seen in Fig. 3, it is possible to go from the state $|0,0\rangle$ to the resonant state $|0,1\rangle$ in two jumps: $|0,0\rangle\to|1,0\rangle\to|0,1\rangle$, which have transition amplitudes (matrix elements) $\varepsilon$ and $g$, with the intermediate state separated by the energy difference $\Delta$. Therefore, the amplitude of this transition (effective coupling between states $|0,0\rangle$ and $|0,1\rangle$) is

[TABLE]

which coincides with Eq. (15).

For the transition between states $|1,0\rangle$ and $|1,1\rangle$, there are two two-jump paths: via the state $|2,0\rangle$ (which is higher in energy by $\Delta-\eta_{\rm c}$) and via $|0,1\rangle$ (which is lower in energy by $\Delta$). Adding these two amplitudes, we obtain

[TABLE]

which coincides with Eq. (16).

Similarly, adding the amplitudes for the paths $|n,0\rangle\to|n+1,0\rangle\to|n,1\rangle$ and $|n,0\rangle\to|n-1,1\rangle\to|n,1\rangle$, we obtain

[TABLE]

which coincides with Eq. (18).

III.2 Next-order analytics

Numerical results for the effective drive amplitudes $\tilde{\varepsilon}_{n}$ (discussed later) show that $\tilde{\varepsilon}_{n}$ is proportional to the actual drive amplitude $\varepsilon$ [as expected from Eq. (18)] only in some range of $\varepsilon$-values. A minor deviation from the linearity at very small $\varepsilon$ (discussed later) is due to dependence of the target-qubit frequency on the control-qubit state – see Eq. (10). The deviation from linearity at large $\varepsilon$ is much more important for practice since it makes it impossible to shorten the CNOT gate duration beyond some value by simply increasing the drive amplitude.

In order to understand the reason for the deviation from linearity at large $\varepsilon$, in this section we develop the next-order analytics for $\tilde{\varepsilon}_{0}$ and $\tilde{\varepsilon}_{1}$, which gives corrections compared with Eqs. (15) and (16). Note that a similar next-order analytics has been developed in Ref. Magesan2018 , though in a very different way (after a misprint correction, the result of Ref. Magesan2018 coincides with our result).

The simple analytics (15)–(18) has been obtained from Eq. (14), which treats the drive Hamiltonian $H_{\varepsilon}$ as a small perturbation. However, for a large drive amplitude $\varepsilon$, the eigenbasis of $H_{\rm qb}+H_{g}$ is no longer the appropriate eigenbasis; instead, $H_{\rm qb}+H_{\varepsilon}$ is the main Hamiltonian, while $H_{g}$ is the perturbation. Note that in the linear approximation, the same $\tilde{\varepsilon}_{n}$ as in Eq. (14) can be obtained by exchanging the roles of $H_{g}$ and $H_{\varepsilon}$, i.e., by using

[TABLE]

where $\overline{|n,m\rangle_{\varepsilon}}$ denotes the eigenstate of $H_{\rm qb}+H_{\varepsilon}$. This equivalence is clear from the discussed above approach of virtual-state transitions, which treats $H_{g}$ and $H_{\varepsilon}$ on equal footing.

For a large $\varepsilon$, Eq. (24) is more appropriate than Eq. (14) to calculate $\tilde{\varepsilon}_{n}$. Even though the initial and final states should still be treated in the eigenbasis of $H_{\rm qb}+H_{g}$, during the front and rear ramps of the microwave pulse the appropriate eigenbases essentially transform into each other, leading to Eq. (24). While we do not have a rigorous justification of the approximation (24) (only a general understanding in the spirit of the adiabatic theorem), numerical results confirm its good accuracy.

Since for the eigenstates $\overline{|n,m\rangle_{\varepsilon}}$ used in Eq. (24) the ladders in Fig. 3 are uncoupled, we can write

[TABLE]

where $\overline{|n\rangle_{\varepsilon}}$ are the control-qubit eigenstates, which account for the drive. They satisfy the Schrödinger equation

[TABLE]

with the Hamiltonian for only the control qubit,

[TABLE]

Then solving this Schrödinger equation and finding the eigenstates,

[TABLE]

we find the effective drive amplitudes $\tilde{\varepsilon}_{n}$ from Eq. (24) as (see Fig. 3)

[TABLE]

Let us use this approach to find $\tilde{\varepsilon}_{0}$ up to the order $\varepsilon^{3}$ [instead of $\varepsilon^{1}$ in the linear approximation (15)], treating $H_{\varepsilon}$ as a perturbation of $H_{\rm qb}$. The eigenstate $\overline{|0\rangle_{\varepsilon}}$ of the control qubit can be written as

[TABLE]

where $\cal N$ is a normalization. Substituting this form into the Schrödinger equation (26) and equating the coefficients for the basis states $|0\rangle$, $|1\rangle$, and $|2\rangle$, we obtain

[TABLE]

where for brevity $E=E_{\overline{|0\rangle_{\varepsilon}}}$, $E_{n}=E_{n}^{\rm(c)}$, we used $E_{0}=0$ [as in Eq. (4)] and also assumed that $\varepsilon$ is real.

To the lowest order, assuming small $\varepsilon$ (therefore $\gamma\ll\beta\ll\alpha$ and $E\approx 0$) we crudely find

[TABLE]

Using these values for $\beta$ and $E$ in Eq. (32), we obtain a better approximation (up to $\varepsilon^{3}$) for $\alpha$:

[TABLE]

To find $\tilde{\varepsilon}_{0}$ with accuracy up to $\varepsilon^{3}$, we need $\alpha$ with accuracy up to $\varepsilon^{3}$, $\beta$ with accuracy up to $\varepsilon^{2}$ and $\cal N$ with accuracy up to $\varepsilon^{2}$, while $\gamma$ is not needed [see Eq. (29)]. Thus, we use Eq. (35) for $\alpha$, Eq. (34) for $\beta$, and ${\cal N}\approx 1+(\varepsilon/E_{1})^{2}/2$ to obtain

[TABLE]

The energy of state $\overline{|0\rangle_{\varepsilon}}$ (not needed for this derivation but needed later) is

[TABLE]

Finally, using Eqs. (29) and (36), we obtain

[TABLE]

with accuracy up to $\varepsilon^{3}$. Note that $E_{1}=E_{1}^{\rm(c)}=\Delta+\delta$, $E_{2}=E_{2}^{\rm(c)}=2(\Delta+\delta)-\eta_{\rm c}$, and we can neglect $\delta$ (i.e., use $\delta=0$). For a complex $\varepsilon$, we need to replace $\varepsilon^{2}$ in parentheses with $|\varepsilon|^{2}$.

Calculation of $\tilde{\varepsilon}_{1}$ up to the order $\varepsilon^{3}$ is similar and requires finding $\overline{|1\rangle_{\varepsilon}}$. Note that the calculations are easier if the energies are counted from $E_{1}$, because then the eigenenergy $E$ in equations similar to Eqs. (31)–(33) is small. The calculations give

[TABLE]

where $E_{nn^{\prime}}\equiv E_{n}-E_{n^{\prime}}=E_{n}^{\rm(c)}-E_{n^{\prime}}^{\rm(c)}$. The corresponding energy is

[TABLE]

Using Eqs. (29) and (39), we obtain

[TABLE]

with accuracy up to $\varepsilon^{3}$ (note the use of $E_{10}$ instead of $E_{01}$ in the preceding formulas). In this formula $E_{10}=\Delta+\delta$, $E_{21}=\Delta+\delta-\eta_{\rm c}$, $E_{31}=2(\Delta+\delta)-3\eta_{\rm c}$, and we can neglect $\delta$ (i.e., $\delta=0$). For a complex $\varepsilon$, we need to replace $\varepsilon^{2}$ in parentheses with $|\varepsilon|^{2}$.

Figure 4 shows the effective drive amplitudes $\tilde{\varepsilon}_{0}$ and $\tilde{\varepsilon}_{1}$ as functions of the actual drive amplitude $\varepsilon$ calculated in several ways for the following parameters (which are some typical experimental parameters): $g/2\pi=3$ MHz, $\eta_{\rm c}/2\pi=300$ MHz, $\Delta/2\pi=130$ MHz. The blue lines (which initially go down) show $\tilde{\varepsilon}_{0}$, the orange lines (which initially go up) show $\tilde{\varepsilon}_{1}$. The solid lines without symbols are calculated using Eqs. (38) and (41) (using $\delta=0$), while the straight dashed lines represent the simple linear approximation, Eqs. (15) and (16). The solid lines with symbols show the numerical results (the numerical procedure is described later in Sec. IV, for numerics we assume $\eta_{\rm t}=\eta_{\rm c}$ and $\omega_{\rm d}=\omega_{\rm t}^{\rm c0}$).

We see that the third-order approximation [Eqs. (38) and (41)] correctly describes the deviation of the dependences $\tilde{\varepsilon}_{0}(\varepsilon)$ and $\tilde{\varepsilon}_{1}(\varepsilon)$ from ideal straight lines at relatively small $\varepsilon$, but fails to fit well the case of relatively large $\varepsilon$. This is because higher-order terms become important even for moderate values of $\varepsilon$. The problem has a similarity with the poor performance of the perturbation approach in the analysis of the circuit QED measurement of transmons when the number of photons is comparable to the critical number.

Figure 5 shows the difference $\tilde{\varepsilon}_{1}-\tilde{\varepsilon}_{0}$ as a function of $\varepsilon$ (for the same parameters as in Fig. 4), also calculated in several ways. The dashed straight line corresponds to the simple formula (19) for the ideal operation. The numerical results are shown by the blue solid line with symbols. The green solid line is calculated using Eqs.  (38) and (41). We have checked (analytically and numerically) that this line coincides with the result given by Eq. (4.25) of Ref. Magesan2018 (after correction of a misprint in the initial version of Ref. Magesan2018 ; the translation of notations is $g=J_{\text{\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{Magesan2018}{\@@citephrase{(}}{\@@citephrase{)}}}}}$, $2\varepsilon=\Omega_{\text{\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{Magesan2018}{\@@citephrase{(}}{\@@citephrase{)}}}}}$, $\Delta=\Delta_{\text{\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{Magesan2018}{\@@citephrase{(}}{\@@citephrase{)}}}}}$, $\eta_{\rm c}=-\delta_{1\,\text{\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{Magesan2018}{\@@citephrase{(}}{\@@citephrase{)}}}}}$, $\tilde{\varepsilon}_{0}-\tilde{\varepsilon}_{1}=(ZX/2)_{{\rm coeff}\,\text{\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{Magesan2018}{\@@citephrase{(}}{\@@citephrase{)}}}}}$). Most importantly, from Fig. 5 we see that the third-order approximation correctly describes initial deviation of the CR gate speed from the linear-order result, but cannot be used for quantitative analysis in the practically interesting regime of large drive amplitudes.

III.3 Semi-analytical results

The method developed in the previous section can be naturally extended to arbitrary large drive amplitudes $\varepsilon$. For that the eigenstates $\overline{|n\rangle_{\varepsilon}}$ of the control-qubit Hamiltonian (27) can be found numerically, and then the effective drive amplitudes $\tilde{\varepsilon}_{n}$ can be calculated using Eq. (29). Since numerical diagonalization of a Hamiltonian for few levels is very easy (compared with full numerical simulation of the two-qubit evolution discussed in the next section), we call this method semi-analytical.

Figure 6 shows a comparison of the semi-analytical results (dashed lines) for $\tilde{\varepsilon}_{0}$ and $\tilde{\varepsilon}_{1}$ with the numerical results (solid lines with symbols, the numerical procedure is discussed in Sec. IV). The parameters are the same as in Figs. 4 and 5, except we use two values of the detuning: $\Delta/2\pi=130$ MHz and 190 MHz. In the semi-analytics, we use 7 levels of the control qubit. We see that the numerical results agree with semi-analytics very well for all values of the drive amplitude $\varepsilon$ (the lines are practically indistinguishable, except for the lowest lines at around 70 MHz, where a minor difference is caused by a resonance between levels $\overline{|0,1\rangle_{\varepsilon}}$ and $\overline{|1,2\rangle_{\varepsilon}}\,$). Similarly, we found a very good agreement for other values of the parameters as well. Therefore, the semi-analytical method based on Eqs. (26)–(29) seems to be a sufficiently simple and accurate way of analyzing the dependence of the CR gate speed on parameters.

Note that $\tilde{\varepsilon}_{n}$ in the semi-analytical method is proportional to the qubit-qubit coupling $g$ and also depends on two dimensionless ratios: $\Delta/\eta_{\rm c}$ and $\varepsilon/\eta_{\rm c}$ (assuming $\delta=0$). In the Duffing (Kerr) approximation (4), these two ratios fully define the eigenstates (28) (in a better approximation Khezri2018 ; Sank2016 the results would also depend on the dimensionless parameter $\eta_{\rm c}/\omega_{\rm c}$). Therefore, in our analysis the ratio $\tilde{\varepsilon}_{n}/g$ is a function of only two parameters: $\Delta/\eta_{\rm c}$ and $\varepsilon/\eta_{\rm c}$.

Figure 7 shows the dimensionless speed $(\tilde{\varepsilon}_{1}-\tilde{\varepsilon}_{0})/g$ of the CR gate as a function of the dimensionless drive amplitude $\varepsilon/\eta_{\rm c}$ for several values of the dimensionless detuning $\Delta/\eta_{\rm c}$ (the lines are calculated using the semi-analytical method). While the behavior at small $\varepsilon$ agrees with Eq. (19) (not shown), the behavior at large $\varepsilon$ mostly depends on whether the detuning $\Delta=\omega_{\rm c}-\omega_{\rm t}$ is negative or positive and on the integer part of the ratio $2\Delta/\eta_{\rm c}$ for positive $\Delta$. As seen in Fig. 7, at large $\varepsilon$ the lines group according to the interval to which $\Delta$ belongs: $(-\infty,0)$, $(0,\eta_{\rm c}/2)$, $(\eta_{\rm c}/2,\eta_{\rm c})$, $(\eta_{\rm c},3\eta_{\rm c}/2)$, $(3\eta_{\rm c}/2,2\eta_{\rm c})$, etc. (in Fig. 7 these groups of lines are labeled sequentially as I, II, III, etc.). We do not show the lines for $\Delta/\eta_{\rm c}=0$, 1/2, 1, 3/2, etc. because at these values there is a resonance between the levels, $E^{\rm(c)}_{n}=E^{\rm(c)}_{0}$ and $E^{\rm(c)}_{n-1}=E^{\rm(c)}_{1}$ for $n=2\Delta/\eta_{\rm c}+1$ [see Eq. (4) for $\delta=0$], and correspondingly the CR gate does not operate as intended (due to a very large leakage – see below), also leading to computational problems in the semi-analytical and numerical calculations.

It is simple to understand why the lines in Fig. 7 group into “bands” at large $\varepsilon$. In the solution of the Schrödinger equation (26) for the Hamiltonian (27) at large $\varepsilon$, the main effect is a strong level repulsion, which depends on the relative position (topology) of the bare energy levels $E_{n}^{\rm(c)}$ (i.e., which level is in between which levels; this topology does not change with $\varepsilon$ because of the adiabatic theorem). In contrast, the level repulsion does not depend much on a particular value of the initial bare level difference, since the effect of $\varepsilon$ dominates. Therefore, at very large $\varepsilon$ the eigenstates (28) do not depend on a particular value of $\Delta/\eta_{\rm c}$, but only on the integer part of $2\Delta/\eta_{\rm c}$ (for $\Delta>0$), which defines the topological structure of the ladder $E_{n}^{\rm(c)}$ (as mentioned above, the bare levels $E_{n}^{\rm(c)}$ intersect at integer values of $2\Delta/\eta_{\rm c}$). Consequently, at very large $\varepsilon$ the effective drive amplitudes (29) depend only on the integer part of $2\Delta/\eta_{\rm c}$. This is why there is a grouping of lines in Fig. 7.

We see that most of the lines in Fig. 7 (all the lines in the experimentally important groups II and III) have a maximum, with a relatively minor decrease after it (experimental results Chow2011 ; Ware2015 are somewhat similar). Experimentally, faster speed $(\tilde{\varepsilon}_{1}-\tilde{\varepsilon}_{0})$ means shorter CNOT gate duration and therefore there is no benefit to increase the drive amplitude beyond the maximum in Fig. 7.

The solid blue line in Fig. 8 shows the maximum value $(\tilde{\varepsilon}_{1}-\tilde{\varepsilon}_{0})_{\rm max}/g$ of the dimensionless speed (or the minimum value for the negative speed) as a function of the dimensionless detuning $\Delta/\eta_{\rm c}$. The dimensionless drive amplitude $\varepsilon/\eta_{\rm c}$ at which this maximum is reached is shown by the dashed orange line. We see that the maximum CR gate speed is reached for the detuning $\Delta$ between $\eta_{\rm c}/2$ and $\eta_{\rm c}$ (group III in Fig. 7). Note that our semi-analytical approach cannot be applied in close vicinities of the detunings $\Delta/\eta_{\rm c}=0$, 0.5, 1, 1.5, etc. The dependence of $(\tilde{\varepsilon}_{1}-\tilde{\varepsilon}_{0})_{\rm max}$ on $\Delta$ was measured experimentally Ware2015 and it showed a crudely similar behavior.

IV Numerical approach

Numerically, we simulate the quantum evolution due to the rotating-frame Hamiltonian (1)–(8), taking into account 7 levels in the control qubit (the same number as in the semi-analytics) and 5 levels in the target qubit, so that there are 35 levels in total. The simulation is based on matrix exponentiation for a time-dependent Hamiltonian, using the second-order Magnus expansion Magnus1954 . We also tried to use the fourth-order Runge-Kutta method, but found that for our typical parameters it is almost an order of magnitude slower to reach the same desired accuracy.

We start with diagonalization of the time-independent part $H_{\rm qb}+H_{g}$ of the Hamiltonian, and then the whole simulation is done in the eigenbasis $\overline{|n,m\rangle}$ of $H_{\rm qb}+H_{g}$, with the time-dependent drive Hamiltonian $H_{\varepsilon}(t)$ (expanded in the eigenbasis) causing the evolution. In this way, we obtain a $35\times 35$ unitary evolution matrix $V$ (in the eigenbasis $\overline{|n,m\rangle}\,$) for a given pulse of the drive amplitude $\varepsilon(t)$ with duration $\tau_{\rm p}$ (the pulse shape is discussed later). The matrix $V$ is then projected onto the computational two-qubit subspace, thus producing a $4\times 4$ matrix $M$, which is no longer unitary (here projection means the simple elimination of all other elements). Note that the reduced matrix $M$ is still defined in the eigenbasis, consisting of states $\overline{|0,0\rangle}=|0,0\rangle$, $\overline{|0,1\rangle}$, $\overline{|1,0\rangle}$ and $\overline{|1,1\rangle}$.

To find fidelity of an operation, we compare the reduced matrix $M$ with the desired $4\times 4$ unitary operation, which we denote $U$. The fidelity between $M$ and $U$ is defined as Zanardi2004 ; Pedersen2007

[TABLE]

where $d=4$ is the dimension of the two-qubit Hilbert space. This definition of the gate fidelity is equal to the final-state fidelity (squared overlap) averaged over all (pure) initial states in the two-qubit subspace; therefore, $F_{MU}$ corresponds to the fidelity in Randomized Benchmarking (assuming that the states leaked outside the computational subspace never come back).

Even though the final goal of the CR gate operation is to produce CNOT (after additional single-qubit rotations), in the numerical procedure the desired $U$ is obviously not the CNOT. Instead, for a given pulse $\varepsilon(t)$ (which produces some matrix $V$ and corresponding matrix $M$), we define $U$ as the closest two-qubit unitary (i.e., which maximizes the fidelity $F_{MU}$), restricted to the following class:

[TABLE]

where $|n\rangle\langle n|_{\text{c}}$ acts on the control qubit, while $X_{\text{t}}$ acts on the target qubit (as mentioned above, we use the eigenbasis of $H_{\rm qb}+H_{g}$ for both $M$ and $U$). The condition (43) means that the state of the control qubit does not change (in the eigenbasis). Also, for state $|0\rangle$ of the control qubit, the target qubit is rotated about $x$-axis by angle $\varphi_{0}$; similarly, for control-qubit state $|1\rangle$, the target qubit is rotated about $x$-axis by angle $\varphi_{1}$. Besides that, in Eq. (43) there are phases $\theta_{0}$ and $\theta_{1}$; disregarding the unimportant overall phase, this can be interpreted as $z$-rotation of the control qubit by angle $\theta_{1}-\theta_{0}$. Without loss of generality, we could assume $\theta_{0}=0$ (while keeping the same $\theta_{1}-\theta_{0}$) since this affects only the overall phase, and the definition (42) of the fidelity $F_{MU}$ is insensitive to the overall phase of $U$. Note that Eqs. (42) and (43) can be easily generalized to include the third state of the control qubit (to consider it as a qutrit); however, here we consider only the two-level subspace.

Thus, to find $U$ for a given pulse of $\varepsilon(t)$, we maximize $F_{MU}$ over parameters $\varphi_{0}$, $\varphi_{1}$, and $\theta_{1}-\theta_{0}$. Fortunately, these optimal angles are given by analytical formulas in terms of the matrix elements of $M$:

[TABLE]

where the rows (and columns) 1, 2, 3, and 4 of the matrix $M$ correspond to the states $\overline{|00\rangle}$, $\overline{|01\rangle}$, $\overline{|10\rangle}$, and $\overline{|11\rangle}$, respectively. In this way, for a given pulse $\varepsilon(t)$, the CR gate operation is characterized by 4 resulting parameters: the angles $\varphi_{0}$, $\varphi_{1}$, and $\theta_{0}-\theta_{1}$ of the unitary (43) and also infidelity $1-F_{MU}$, which is due to leakage outside of the computational two-qubit subspace and also due to the computational-space unitary not fitting well the class (43).

We consider the pulse shape $\varepsilon(t)$ of duration $\tau_{\rm p}$,

[TABLE]

which consists of the flat middle part with the real amplitude $\varepsilon_{\rm m}$ of the drive and two symmetric cosine-shaped ramps (so that there are no kinks), each with duration $\tau_{\rm r}$ – see Fig. 9. As discussed later, sufficiently long ramps are needed to reduce leakage outside of the computational subspace.

Effective drive amplitudes $\tilde{\varepsilon}_{0}$ and $\tilde{\varepsilon}_{1}$ used in Sec. III (Figs. 4–6) have been numerically calculated as the derivatives,

[TABLE]

while keeping the ramp duration $\tau_{\rm r}$, the middle-part amplitude $\varepsilon_{\rm m}$ (which replaces $\varepsilon$ in Sec. III), and other parameters fixed.

The small drive frequency detuning $\delta$ in numerical simulations is chosen in the following way. We first use the laboratory frame, i.e., $\delta=\omega_{\rm t}$ (for the Duffing oscillator model we can also use $\delta=\omega_{\rm t}=0$) and calculate the eigenfrequencies of the target qubit $\omega_{\rm t}^{\rm c0}$ and $\omega_{\rm t}^{\rm c1}$ – see Eq. (9). If we want the drive to be exactly on-resonance with the target qubit when the control qubit is $|0\rangle$, then we need to use $\omega_{\rm d}=\omega_{\rm t}^{\rm c0}$, which gives $\delta=\omega_{\rm t}-\omega_{\rm t}^{\rm c0}$. Similarly, if we want $\omega_{\rm d}=\omega_{\rm t}^{\rm c1}$ (the drive on-resonance with the target qubit when the control qubit is $|1\rangle$), then we use $\delta=\omega_{\rm t}-\omega_{\rm t}^{\rm c1}$. If we want the drive frequency exactly in between the two resonances, $\omega_{\rm d}=(\omega_{\rm t}^{\rm c0}+\omega_{\rm t}^{\rm c1})/2$, then we use $\delta=\omega_{\rm t}-(\omega_{\rm t}^{\rm c0}+\omega_{\rm t}^{\rm c1})/2$. Note that the frequency differences $\omega_{\rm t}-\omega_{\rm t}^{\rm c0}$ and $\omega_{\rm t}-\omega_{\rm t}^{\rm c1}$ do not depend on a choice of the rotating frame.

Since experimentally the CR gate is mainly used to realize CNOT, in the numerical simulations we are mostly interested in the operations equivalent to CNOT up to single-qubit rotations. In analyzing the CNOT-equivalent gates, we usually use the pulse shape in which the ramps occupy 30% of the whole pulse duration each, i.e., $\tau_{\rm r}=0.3\,\tau_{\rm p}$ in Eq. (48) (Fig. 21 is an exception). For a given middle-part amplitude $\varepsilon_{\rm m}$, we find the shortest pulse duration $\tau_{\rm p}$, for which

[TABLE]

This is what we call the CNOT gate duration $\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}(\varepsilon_{\rm m})$, neglecting durations of the additional single-qubit operations ($x$-rotation of the target qubit and $z$-rotation of the control qubit). We assume perfect fidelity of single-qubit operations; therefore, the CNOT infidelity is $1-F_{MU}$ for the pulse with duration $\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}(\varepsilon_{\rm m})$.

In the simulations we fully neglect decoherence. A crude estimate of the fidelity decrease $\Delta F$ due to energy relaxation and pure dephasing can be obtained by considering idle qubits, which decohere during time $\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}$. This gives the estimate

[TABLE]

where $T_{1}^{\rm(c)}$ and $T_{1}^{\rm(t)}$ are the energy relaxation times for the control and target qubits, and similarly $T_{2}^{\rm(c)}$ and $T_{2}^{\rm(t)}$ are the dephasing times (which include contributions due to the energy relaxation and pure dephasing, $1/T_{2}=1/2T_{1}+1/T_{\varphi}$). This estimate is obtained by summing the single-qubit Pauli error rates $t/2T_{1}$ and $t/2T_{\varphi}$ and then converting the result into the 2-qubit average gate fidelity by using the factor $4/5$. Note, however, that actual fidelity decrease $\Delta F$ can be significantly larger than the estimate (51) because the CR gate operation involves a significant population of the level $|2\rangle$ and even higher levels of the control qubit, which have poorer coherence than the level $|1\rangle$.

One run of the evolution simulation for a given pulse duration typically takes a few seconds on a desktop or a laptop computer. Finding $\tau_{\text{p}}^{\rm\scriptscriptstyle CNOT}$ requires several tens of runs, so a typical time to produce a line showing dependence of the CNOT gate operation on $\varepsilon_{\text{m}}$ is few hours. The simulation time significantly depends on the number of time steps in the pulse ramps; we have used 600 time steps for each ramp, which gives a quite good accuracy for the simulations. For quick (and much less accurate) simulations it is possible to use $\sim$100 time steps per ramp.

Let us list the main approximations used in our numerics: (i) neglected decoherence, (ii) Duffing-oscillator approximation for the transmon energy levels, (iii) linear-oscillator approximation for the transmon matrix elements, (iv) direct coupling of qubits instead of coupling via resonator, (v) RWA, (vi) using only $7\times 5$ levels (we checked that this is sufficient), (vii) no microwave crosstalk (except in Fig. 22), (viii) simple pulse shape without distortions, (ix) absence of neighboring qubits. In spite of a rather long list of approximations, we believe our simulation results give a reasonably accurate description of the intrinsic operation of the CR gate (neglecting decoherence, which in practice may give the biggest contribution to infidelity).

V Numerical CNOT gate duration and single-qubit rotations

In numerical simulations, we use the qubit-qubit coupling $g/2\pi=3$ MHz (except in Fig. 20) and transmon anharmonicity $\eta_{\rm c}/2\pi=\eta_{\rm t}/2\pi=300$ MHz. Figure 10 shows the CNOT gate duration $\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}$ as a function of the drive amplitude $\varepsilon_{\rm m}$ in the flat middle part of a pulse (with $\tau_{\rm r}=0.3\,\tau_{\rm p}$), for several values of the qubit-qubit detuning: $\Delta/2\pi=-70$, 70, 130, and 190 MHz. For numerical results (solid lines) we choose $\omega_{\rm d}=\omega_{\rm t}^{\rm c0}$, i.e., $\delta=\omega_{\rm t}-\omega_{\rm t}^{\rm c0}$. The dashed lines (almost coinciding with the solid lines) show the result of the semi-analytical method, in which we use Eq. (13) and integrate over the pulse shape. We see that the semi-analytical method works very well; however, there are (barely) visible deviations at both small and large amplitudes $\varepsilon_{\rm m}$. We guess the slight deviation at large $\varepsilon_{\rm m}$ is because for a short pulse the non-adiabatic evolution during the ramps starts to play a noticeable role. The deviation at small $\varepsilon_{\rm m}$ is because here the $zz$-coupling (10) starts to play a relatively significant role.

For a more detailed analysis of the deviation at small $\varepsilon_{\rm m}$, solid lines in Fig. 11 show the numerical CNOT time $\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}$ for $\Delta/2\pi=130$ MHz and three values of the drive frequency: $\omega_{\rm d}=\omega_{\rm t}^{\rm c0}$ (blue line, drive on-resonance with the target qubit when the control qubit is $|0\rangle$, i.e., $\delta=\omega_{\rm t}-\omega_{\rm t}^{\rm c0}$), $\omega_{\rm d}=\omega_{\rm t}^{\rm c1}$ (orange line, on-resonance when the control qubit is $|1\rangle$, i.e., $\delta=\omega_{\rm t}-\omega_{\rm t}^{\rm c1}$), and exactly in between, $\omega_{\rm d}=(\omega_{\rm t}^{\rm c0}+\omega_{\rm t}^{\rm c1})/2$ (green line). The dashed line shows the semi-analytical result (actually, there are three dashed lines for the three values of $\delta$, but they are indistinguishable), and the dotted line shows the ideal result,

[TABLE]

which follows from Eqs. (13) and (19) after integration over the pulse shape (48) with $\tau_{\rm r}/\tau_{\rm p}=0.3$ (this integration gives the factor 0.7). We see that the solid lines noticeably deviate down from the dashed line for $\varepsilon_{\rm m}/2\pi$ less than $\sim$5 MHz, with the largest deviation for $\omega_{\rm d}=\omega_{\rm t}^{\rm c0}$ (blue line). Note that for $\varepsilon/2\pi=5$ MHz, Eq. (19) gives $(\tilde{\varepsilon}_{1}-\tilde{\varepsilon}_{0})/2\pi=0.41$ MHz, while $\omega_{zz}/2\pi=0.15$ MHz [see Eq. (11)], so the effect of $zz$-coupling is expected to be significant.

In more detail, for $\omega_{\rm d}=\omega_{\rm t}^{\rm c0}$ (blue line in Fig. 11) the drive is exactly on-resonance with the target qubit for the control-qubit state $|0\rangle$ and therefore the approximation $\varphi_{0}=0.7\,\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}\times 2\tilde{\varepsilon}_{0}(\varepsilon_{\rm m})$ using Eq. (15) should work well (as confirmed by numerics). In contrast, for the control-qubit state $|1\rangle$, the drive is detuned by $\omega_{zz}$ from the target-qubit frequency, which leads to Rabi oscillations with frequency $\sqrt{(2\tilde{\varepsilon}_{1})^{2}+\omega_{zz}^{2}}$ within the plane tilted by angle ${\rm atan}(\omega_{zz}/2\tilde{\varepsilon}_{1})$ from the $zy$ plane (the frequency and the plane are changing in time because of the pulse shape). The larger Rabi frequency (due to $\omega_{zz}$ contribution) leads to a slightly shorter CNOT gate duration than expected analytically, explaining the behavior of solid lines in Fig. 11 at small $\varepsilon_{\rm m}$. The same effect (rotation within a tilted plane) leads to a significant infidelity of the CNOT gate at small $\varepsilon_{\rm m}$ – see Sec. VI.

As discussed above, in order to realize the CNOT operation, the CR gate with the pulse duration $\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}$ should be complemented by single-qubit rotations. The target qubit should be rotated about $x$-axis by the angle $-\varphi_{0}$ to compensate the operator $e^{-i(\varphi_{0}/2)X_{\rm t}}$ in Eq. (43). Similarly, the control qubit should be rotated about $z$-axis to compensate the relative phase $\theta_{1}-\theta_{0}$ and the negative imaginary unit due to the relation $e^{-i(\pi/2)X_{\rm t}}=-iX_{\rm t}$. So, naively we would expect that $z$-rotation by the angle $\pi/2-(\theta_{1}-\theta_{0})$ is needed. However, the angle $\theta_{1}-\theta_{0}$ is numerically computed in the rotating frame of the drive, while experimental $z$-rotation should be in the rotating frame of the control qubit. For the latter frame based on frequency $\omega_{\rm c}^{\rm t0}=E_{\overline{|1,0\rangle}}^{\rm(lf)}-E_{\overline{|0,0\rangle}}^{\rm(lf)}$ (i.e., when the target qubit is $|0\rangle$, as usually done in experiments), we need to replace the phase difference $\theta_{1}-\theta_{0}$ with

[TABLE]

so in an experiment, the control qubit should be $z$-rotated by the angle $\theta_{0}^{\prime}-\theta_{1}^{\prime}+\pi/2$. Note that $\omega_{\rm c}^{\rm t0}-\omega_{\rm d}=\Delta+\delta+\omega_{\rm c}^{\rm t0}-\omega_{\rm c}$.

Figure 12 shows the angle $-\varphi_{0}$ (normalized by $\pi$) as a function of $\varepsilon_{\rm m}$ for two values of the detuning: $\Delta/2\pi=130$ and 190 MHz (we use $\omega_{\rm d}=\omega_{\rm t}^{\rm c0}$). Solid lines show the numerical results, dashed lines show the corresponding semi-analytical results (integrating $-2\tilde{\varepsilon}_{0}$ over the pulse shape), and horizontal dotted lines show the ideal result based on Eqs. (15) and (16): $-\varphi_{0}=\pi(\eta_{\rm c}-\Delta)/2\eta_{\rm c}$. We see that this ideal result for $-\varphi_{0}$ is never applicable. The deviation from it at large $\varepsilon_{\rm m}$ is described well by the semi-analytics and is due to the deviations from the ideal straight lines in Fig. 6. The deviation from the ideal result in Fig. 12 at small $\varepsilon_{\rm m}$ is due to $\omega_{zz}$ – the same effect as discussed above for Fig. 11.

Upper (blue) solid line in Fig. 13 shows the numerical result for the control-qubit rotation angle $\theta_{0}^{\prime}-\theta_{1}^{\prime}+\pi/2$ (normalized by $\pi$) as a function of $\varepsilon_{\rm m}$ for $\Delta/2\pi=130$ MHz and $\omega_{\rm d}=\omega_{\rm t}^{\rm c0}$. We see that it is quite different from the ideally expected value of $\pi/2$ (horizontal brown line). The difference is mainly due to two effects. First, strong drive $\varepsilon(t)$ causes the level repulsion (ac Stark shift) in the control qubit, which slightly changes the control-qubit frequency and produces the accumulated angle

[TABLE]

(Actually, because the state $|1\rangle$ is assumed to be at the bottom of the Bloch sphere, this produces the $z$-rotation of $-\theta_{\rm rep}$, so we need to apply the rotation of $+\theta_{\rm rep}$ to compensate it.) The angle $\theta_{\rm rep}$ calculated semi-analytically [see Eq. (26)] is shown by the dashed red line in Fig. 13; the same angle calculated using analytical results (37) and (40) is shown by the dotted red line (in both cases we use numerical $\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}$ to integrate over the pulse shape). The second contribution into $\theta_{0}^{\prime}-\theta_{1}^{\prime}$ is due to $zz$ coupling (10), which produces

[TABLE]

This is because the rotating frame here is based on $\omega_{\rm c}^{\rm t0}$, while $\omega_{\rm c}^{\rm t1}-\omega_{\rm c}^{\rm t0}=\omega_{zz}$ and both states of the target qubit participate equally (leading to the factor of $1/2$). The angle $\theta_{zz}$ is shown by the lower (green) solid line in Fig. 13. Adding the three contributions, $\theta_{\rm rep}+\theta_{zz}+\pi/2$, we obtain the dashed and dotted black lines (corresponding to the dashed and dotted red lines for $\theta_{\rm rep}$), which are quite close to the numerical result (solid blue line). This confirms the main physical mechanisms contributing to $\theta_{0}^{\prime}-\theta_{1}^{\prime}$ and also shows that the approximation based on Eqs.  (37) and (40) works quite well.

VI Error budget

In the previous section, we have discussed numerical results for the parameters of a CR-based CNOT gate: the duration $\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}$ as a function of the mid-pulse drive amplitude $\varepsilon_{\rm m}$ and the compensating single-qubit rotation angles $-\varphi_{0}$ and $\theta_{0}^{\prime}-\theta_{1}^{\prime}+\pi/2$. We have seen that these parameters can be obtained quite accurately by the semi-analytical method. In this section we discuss numerical results for the infidelity $1-F_{MU}$ of the CR-based CNOT gate (also as a function of $\varepsilon_{\rm m}$), neglecting decoherence and infidelity of single-qubit rotations (see Sec. IV for the definition of $F_{MU}$ and the calculation method). These results cannot be obtained semi-analytically and necessarily require full numerical simulations.

Figure 14 shows the numerical results for the CNOT gate infidelity $1-F_{MU}$ as a function of the mid-pulse drive amplitude $\varepsilon_{\rm m}$ for several values of the detuning: $\Delta/2\pi=-70$, 70, 130, and 190 MHz. As in the previous plots, we use $g/2\pi=3$ MHz, $\eta_{\rm c}=\eta_{\rm t}=300$ MHz, $\omega_{\rm d}=\omega_{\rm t}^{\rm c0}$, and $\tau_{\rm r}=0.3\,\tau_{\rm p}$. Most importantly, we see that the infidelity dependence on $\varepsilon_{\rm m}$ has a minimum, and at this minimum, the infidelity is crudely $10^{-3}$ for all lines. The second observation is that the minimum is not sharp and is reached for the values of the drive amplitude $\varepsilon_{\rm m}$ above which the CNOT duration $\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}$ does not become significantly shorter by a further increase of $\varepsilon_{\rm m}$ (see Figs. 7 and 10). Note that we do not take into account the effect of decoherence, which can be crudely (ideally) estimated by Eq. (51). If decoherence were added, then the minima in Fig. 14 would shift to higher values of $\varepsilon_{\rm m}$; however, since the corresponding decrease of $\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}$ is not significant, the benefit for fidelity is also not significant; moreover, using higher drive amplitudes could lead to other experimental problems. Therefore, we think the optimum values of $\varepsilon_{\rm m}$ in Fig. 14 should be somewhat close to experimental optima.

To understand the reason for the minima in Fig. 14 and to understand the physical mechanisms contributing to the infidelity (error budget), we have done the following more detailed calculations. We remind that we calculate fidelity $F_{MU}$ between the non-unitary $4\times 4$ matrix $M$ (which is the projection of the $35\times 35$ matrix of actual evolution onto the computational subspace) and the closest unitary matrix $U$, which belongs to the class (43) with $\varphi_{1}-\varphi_{0}=\pi$ ($\rm mod\,2\pi$) for a CNOT-equivalent gate. Now let us consider a bigger class of unitary matrices and define $\tilde{M}$ as the matrix, which is closest to $M$ out of unitaries satisfying the condition

[TABLE]

where $\tilde{U}_{0}^{\text{t}}$ and $\tilde{U}_{1}^{\text{t}}$ are any $2\times 2$ unitary matrices acting on the target qubit (here “closest” means that $\tilde{M}$ maximizes the fidelity $F_{M\tilde{M}}$). From this definition, $\tilde{M}$ does not change the control-qubit states $|0\rangle$ and $|1\rangle$ (i.e., does not allow any leakage from them), but its rotation of the target qubit is arbitrary. The matrix $U$ in Eq. (43) is more restrictive in the sense that it allows only $x$-rotations of the target qubit.

To clarify the error budget, we calculate additional infidelities $1-F_{M\tilde{M}}$ (between $M$ and $\tilde{M}$) and $1-F_{\tilde{M}U}$ (between $\tilde{M}$ and $U$) for the CNOT-equivalent CR operation. The idea is that the infidelity $1-F_{M\tilde{M}}$ is due to leakages (from the control-qubit states $|0\rangle$ and $|1\rangle$ to any state and also outside of the computational subspace for the target qubit). In contrast, the infidelity $1-F_{\tilde{M}U}$ is due to imperfect unitary rotations of the target qubit. From the physical approach of separation of the error into these different mechanisms, we would expect

[TABLE]

This is not an exact relation mathematically (the exact relation would require that certain elements of the quantum process tomography matrix for $M$ have exactly zero real and/or imaginary parts). However, numerical results (e.g., Fig. 15) show that this relation is very accurate.

In Fig. 15, the solid blue line, $1-F_{MU}$, is the same as the blue line in Fig. 14 ($\Delta/2\pi=130$ MHz). For the same case, the orange line shows $1-F_{M\tilde{M}}$, the green line shows $1-F_{\tilde{M}U}$, and the dashed blue line (practically indistinguishable from the solid blue line) shows the sum $(1-F_{M\tilde{M}})+(1-F_{\tilde{M}U})$. We see that the total infidelity $1-F_{MU}$ can really be decomposed into the infidelity $1-F_{M\tilde{M}}$ due to leakages and infidelity $1-F_{\tilde{M}U}$ due to imperfect unitaries acting on the target qubit. Similar decomposition into leakages and imperfect unitaries also works well for other lines in Fig. 14 (relative inaccuracy of the decomposition is typically about $10^{-3}$).

Besides introducing the error budget via the decomposition (57), we also tried a further decomposition by introducing another $4\times 4$ matrix $\tilde{M}^{\prime}$, which is a two-qubit unitary closest to $M$. In this case the infidelity $1-F_{M\tilde{M}^{\prime}}$ is due to leakages outside of the computational subspace, while the infidelity $1-F_{\tilde{M}^{\prime}\tilde{M}}$ is due to unitary transitions between states $|0\rangle$ and $|1\rangle$ of the control qubit (which in our terminology are also leakages for the CR gate). Then the infidelity consists of 3 components,

[TABLE]

We have checked that this relation is quite accurate numerically. However, for the cases we checked, the matrix $\tilde{M}^{\prime}$ was very close to either $M$ or $\tilde{M}$ (depending on the detuning $\Delta$, which determines the strongest leakage process). Therefore, usually only two terms in Eq. (58) are significant, and this is why we will continue using the simpler error decomposition (57) below.

As seen in Fig. 15, at small $\varepsilon_{\rm m}$ the CNOT gate infidelity is dominated by the imperfection of the unitaries $\tilde{U}_{0}^{\rm t}$ and $\tilde{U}_{1}^{\rm t}$ (contribution $1-F_{\tilde{M}U}$), while at large $\varepsilon_{\rm m}$ the leakage contribution $1-F_{M\tilde{M}}$ dominates (the same result for other detunings $\Delta$). This is because at small $\varepsilon_{\rm m}$ the effect of $zz$-coupling is important (as discussed in Sec. V), while at large $\varepsilon_{\rm m}$ the ramps of the pulse become short and high, making the process significantly non-adiabatic and causing leakages.

To clarify the dependence on $\varepsilon_{\rm m}$ of the imperfect-unitary contribution $1-F_{\tilde{M}U}$ (green line in Fig. 15), we draw this line again in Fig. 16 (now on semi-logarithmic scale). We also show the numerical results for the contributions $\Delta F_{\,\tilde{U},\,\rm c0}$ and $\Delta F_{\,\tilde{U},\,\rm c1}$ to this line from imperfections of the unitaries $\tilde{U}_{0}^{\rm t}$ (for the control-qubit state $|0\rangle$) and $\tilde{U}_{1}^{\rm t}$ (for the control-qubit state $|1\rangle$): orange and blue solid lines in Fig. 16, respectively. [Mathematically, the definitions are: $\Delta F_{\,\tilde{U},\,\rm c0}=4/5-(2/5)\,|{\rm Tr}(\tilde{U}_{0}^{\rm t}U_{0}^{\rm t\dagger})|$, $\Delta F_{\,\tilde{U},\,\rm c1}=4/5-(2/5)\,|{\rm Tr}(\tilde{U}_{1}^{\rm t}U_{1}^{\rm t\dagger})|$, where $U_{0}^{\rm t}$ and $U_{1}^{\rm t}$ are obtained from Eq. (43): $U_{0}^{\text{t}}=e^{-i(\varphi_{0}/2)X_{\rm t}}$ and $U_{1}^{\text{t}}=e^{-i(\varphi_{1}/2)X_{\rm t}}$.] At small $\varepsilon_{\rm m}$ the main contribution comes from imperfect $\tilde{U}_{1}^{\rm t}$. This is because we use the drive frequency resonant with the target qubit for the control-qubit state $|0\rangle$ ($\omega_{\rm d}=\omega_{\rm t}^{\rm c0}$), so for the the control-qubit state $|1\rangle$ the drive is off resonance by $\omega_{zz}$ [see Eq. (10)]. This detuning produces rotation of the target-qubit state about a tilted axis, instead of the desired $x$-axis. To check this explanation, we calculate analytically the corresponding infidelity,

[TABLE]

where the factor $25/40$ comes from the integration of $\varepsilon^{2}(t)$ over the pulse shape (48) with $\tau_{\rm r}/\tau_{\rm p}=0.3$ (such integration in only the denominator is an approximation). The dotted blue line in Fig. 16 shows this result with $\varphi_{1}$ and $\tilde{\varepsilon}_{1}$ calculated analytically: $\varphi_{1}=(\eta_{\rm c}+\Delta)/2\eta_{\rm c}$ and Eq. (16) for $\tilde{\varepsilon}_{1}$. The dashed blue line shows Eq. (59) with $\varphi_{1}$ and $\tilde{\varepsilon}_{1}$ calculated using the semi-analytical method. We see that both lines fit reasonably well the numerical result (solid blue line) at small $\varepsilon_{\rm m}$, with a better fit when using the semi-analytical values for $\varphi_{1}$ and $\tilde{\varepsilon}_{1}$. Even better fit (almost perfect, not shown) is when the Bloch-sphere evolution due to $\omega_{zz}$ and semi-analytical $\tilde{\varepsilon}(\varepsilon(t))$ is integrated over the pulse shape numerically, instead of using Eq. (59).

So, the contribution $\Delta F_{\,\tilde{U},\,\rm c1}$ to the gate infidelity due to imperfect unitary $\tilde{U}_{\rm t}^{\rm c1}$ is well explained quantitatively. In contrast, we do not have a simple analytical way to find the contribution $\Delta F_{\,\tilde{U},\,\rm c0}$ due to imperfect $\tilde{U}_{\rm t}^{\rm c0}$ (orange solid line in Fig. 16). Qualitatively, this contribution appears at large $\varepsilon_{\rm m}$ because the interplay between the level repulsions due to $g$ and $\varepsilon$ in Fig. 3 slightly changes the frequency $\omega_{\rm t}^{\rm c0}$, so that there is no longer exact resonance with the drive, and the Bloch-sphere evolution also becomes tilted, thus producing the infidelity. Note that small oscillations of the orange and blue solid lines in Fig. 16 at large $\varepsilon_{\rm m}$ are apparently related to the oscillations of the orange line in Fig. 15, which are discussed below.

The second contribution to the overall gate infidelity $1-F_{MU}$ (Fig. 15), which becomes dominating at large $\varepsilon_{\rm m}$, is the contribution $1-F_{M\tilde{M}}$ due to leakage produced by non-adiabaticity during the pulse ramps (orange line in Fig. 15). From the fidelity definition (42), we expect an approximate relation

[TABLE]

where $P_{\rm leak}^{\rm out}$ is the probability of leakage to outside the computational subspace, averaged over the initial two-qubit states, and $P_{\rm leak}^{\rm comp}$ is the averaged probability of leakage inside the computational subspace, which for our definition (56) means transitions between states $|0\rangle$ and $|1\rangle$ of the control qubit. Note that to average the leakage probability over all initial two-qubit states, it is sufficient to average it over the 4 basis states: $\overline{|0,0\rangle}$, $\overline{|0,1\rangle}$, $\overline{|1,0\rangle}$, and $\overline{|1,1\rangle}$. We have checked the relation (60) numerically and it works very well; however, usually either the first or second term in Eq. (60) strongly dominates over the other term (depending on the dominating leakage channel).

Figure 17 shows again the orange line $1-F_{M\tilde{M}}$ from Fig. 15 (now it is the thick orange line and the scale is semi-logarithmic) and also shows the numerical leakage probabilities (multiplied by the factor $1/4$ as in the averaging) for the processes $\overline{|0,0\rangle}\to\overline{|2,0\rangle}$, $\overline{|0,1\rangle}\to\overline{|2,1\rangle}$, $\overline{|0,0\rangle}\to\overline{|2,1\rangle}$, and $\overline{|0,1\rangle}\to\overline{|2,0\rangle}$ (thin solid lines). The sum of the thin solid lines is shown as the dotted brown line. It is very close to the thick orange line [as expected from Eq. (60)], indicating that these are the four dominating leakage channels. A minor difference between the thick orange line and dotted brown line at $\varepsilon_{\rm m}$ close to 80 MHz is due to the significance of the additional leakage channels $\overline{|0,1\rangle}\to\overline{|1,2\rangle}$ and $\overline{|0,0\rangle}\to\overline{|1,2\rangle}$ (at this frequency the states $\overline{|2,1\rangle}$ and $\overline{|1,2\rangle}$ become on-resonance); similarly, a visible difference at about 40 MHz is because the transitions $|0\rangle\leftrightarrow|1\rangle$ in the control qubit become relatively important.

The main leakage channels in Fig. 17 involve the transition $|0\rangle\to|2\rangle$ in the control qubit. This is expected because in the rotating frame $E_{0}^{\rm(c)}-E_{2}^{\rm(c)}=\eta-2\Delta$ is only 40 MHz for the parameters of Fig. 17. The matrix element of this transition via virtual state $|1\rangle$ is $-\sqrt{2}\,\varepsilon^{2}/\Delta$. Therefore, the non-adiabatic transition $\overline{|0\rangle_{\varepsilon}}\to\overline{|2\rangle_{\varepsilon}}$ during the front ramp with duration $\tau_{\rm r}$ of the pulse (48) has the amplitude proportional to the Fourier transform of $d(\varepsilon^{2})/dt$ at the frequency $\eta-2\Delta$. The standard calculations give an estimate of the leakage probability in the control qubit during the front ramp:

[TABLE]

This probability is shown in Fig. 17 by the black dashed line. We see that it gives a reasonable (crude) approximation of the infidelity $1-F_{M\tilde{M}}$ due to leakage (if also multiplied by the factor of $1/4$, it goes close to the top of oscillating blue and green solid lines). Such a good fit is somewhat surprising because in deriving Eq. (61) we assumed a fixed energy difference $E_{0}^{\rm(c)}-E_{2}^{\rm(c)}=\eta-2\Delta$, while for large $\varepsilon$ the difference of eigenenergies, $E_{\overline{|0\rangle_{\varepsilon}}}-E_{\overline{|2\rangle_{\varepsilon}}}$, becomes significantly larger, e.g., 60.7 MHz for $\varepsilon/2\pi=60$ MHz and 84.3 MHz for $\varepsilon/2\pi=80$ MHz. Therefore, we would expect that Eq. (61) should significantly overestimate the actual leakage (note the sixth power of $(\eta-2\Delta)$ in the denominator). As we checked, Eq. (61) still works well because the non-adiabatic transition mainly accumulates during the lower half of the ramp, when $\varepsilon(t)$ is not too large.

The oscillations of the solid lines in Fig. 17 (leakage probabilities multiplied by $1/4$) are easily understandable. The non-adiabatic leakage occurs during both front and rear ramps, and the transition amplitudes are added with a non-zero phase due to the energy difference $E_{\overline{|0\rangle_{\varepsilon}}}-E_{\overline{|2\rangle_{\varepsilon}}}$, accumulated between the ramps. So, the oscillations are due to constructive or destructive interference of the leakage contributions from the two ramps. We have checked that the $\varepsilon_{\rm m}$-difference between the peaks in Fig. 17 is consistent with estimates based on the numerical increase of $E_{\overline{|0\rangle_{\varepsilon}}}-E_{\overline{|2\rangle_{\varepsilon}}}$ with $\varepsilon$. The oscillations in the probabilities of individual leakage channels lead to the oscillations of their sum, thus explaining oscillations of $1-F_{M\tilde{M}}$ in Fig. 15 and oscillations of the overall CNOT gate infidelity $1-F_{MU}$ in Fig. 14 at large $\varepsilon_{\rm m}$.

Note that in Fig. 17 the leakage channels $\overline{|0,0\rangle}\to\overline{|2,0\rangle}$ and $\overline{|0,1\rangle}\to\overline{|2,1\rangle}$ (with non-changing state of the target qubit) have higher probabilities than for the leakage channels with changing state of the target qubit. This is because $|\varphi_{0}|\ll\pi$ in the interesting range of $\varepsilon_{\rm m}$ – see Fig. 12, so the target-qubit state does not change much during the pulse when the control-qubit state is $|0\rangle$. In the opposite limit, $|\pi-\varphi_{0}|\ll\pi$, we would expect all four leakage channels to have approximately the same strength, and also would not expect significant oscillations (because the two ramps would mainly contribute to different channels).

Besides the detailed analysis of leakage contribution $1-F_{M\tilde{M}}$ for the detuning $\Delta/2\pi=130$ MHz, we have also analyzed the leakage for other values of the detuning in Fig. 14. The case of $\Delta/2\pi=190$ MHz is similar (the same dominating leakage channels) and Eq. (61) still works well; however, $|\eta-2\Delta|$ is larger (80 MHz instead of 40 MHz), so the leakage becomes significant only at larger values of $\varepsilon_{\rm m}$. Also, a very large value of the green line in Fig. 14 at $\varepsilon_{\rm m}/2\pi\simeq 70$ MHz is due to an additional leakage channel $\overline{|0,1\rangle}\to\overline{|1,2\rangle}$, which is due to a resonance between these states. For the detuning $\Delta/2\pi$ of 70 MHz and $-70$ MHz (brown and red lines in Fig. 14), the main leakage is between states $\overline{|0\rangle_{\varepsilon}}$ and $\overline{|1\rangle_{\varepsilon}}$ of the control qubit (also due to non-adiabaticity during the ramps).

Thus, in this section we have shown that in our model the error budget of the CR gate consists of two main contributions: the imperfection of the unitary operation at small drive amplitudes and the leakage at large drive amplitudes. At the optimal drive amplitude, the infidelity is on the order of $10^{-3}$. However, we did not take into account contributions from decoherence and also from possible problems caused by a strong drive of the control qubit (e.g., due to a resonance between an impurity and $E_{\overline{|0\rangle_{\varepsilon}}}\,-E_{\overline{|n\rangle_{\varepsilon}}}\,$ or $E_{\overline{|1\rangle_{\varepsilon}}}\,-E_{\overline{|n\rangle_{\varepsilon}}}$). Note that a strong leakage makes the CR gate operation impractical for the detuning $\Delta$ close to [math], $\eta_{\rm c}/2$, $\eta_{\rm c}$, $3\eta_{\rm c}/2$, etc.

VII Dependence on parameters

A convenient way to present numerical results Kirchhoff2018 is to parametrically plot infidelity $1-F_{MU}$ versus CNOT gate duration $\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}$, with both quantities being functions of the drive amplitude $\varepsilon_{\rm m}$. Figure 18 presents such a plot for the results shown in Figs. 10 and 14 for the detuning values $\Delta/2\pi=-70$, 70, 130, and 190 MHz, while other parameters are $g/2\pi=3$ MHz, $\eta_{\rm c}/2\pi=\eta_{\rm t}/2\pi=300$ MHz, $\omega_{\rm d}=\omega_{\rm t}^{\rm c0}$, and $\tau_{\rm r}/\tau_{\rm p}=0.3$. An increase of $\varepsilon_{\rm m}$ corresponds to moving from right to left along the lines. The main observation is that for each $\Delta$ the line is very steep at the left, which naturally corresponds to a limit on decreasing the duration of the CNOT gate. Among the plotted lines, the line for the detuning of 190 MHz gives the shortest duration (consistent with our discussion in Sec. III.3); however, it does not give the best fidelity.

In Fig. 19 we show the results for three slightly different drive frequencies for the qubit detunings of 70 MHz (brown lines) and 190 MHz (green lines). Besides using the drive frequency resonant with the target qubit when the control-qubit state is $|0\rangle$ (solid lines, $\omega_{\rm d}=\omega_{\rm t}^{\rm c0}$, $\delta=\omega_{\rm t}-\omega_{\rm t}^{\rm c0}$), we also use the drive resonant with the target qubit when the control-qubit state is $|1\rangle$ (dotted lines, $\omega_{\rm d}=\omega_{\rm t}^{\rm c1}$), and also show the case of the drive frequency exactly in between (dashed lines, $\omega_{\rm d}=(\omega_{\rm t}^{\rm c0}+\omega_{\rm t}^{\rm c1})/2$). Note that $\omega_{\rm t}^{\rm c1}-\omega_{\rm t}^{\rm c0}=\omega_{zz}$ is 127 kHz and 200 kHz for the detunings of 70 MHz and 190 MHz, respectively [see Eq. (11)]. We see that this small change of the drive frequency practically does not affect the natural limit for the duration $\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}$; however, it may very significantly affect the infidelity. For example, for $\Delta/2\pi=70$ MHz the minimum infidelity is $1.7\times 10^{-4}$ for $\omega_{\rm d}=(\omega_{\rm t}^{\rm c0}+\omega_{\rm t}^{\rm c1})/2$, while it is $7.7\times 10^{-4}$ for $\omega_{\rm d}=\omega_{\rm t}^{\rm c0}$. We have checked that this small change of the drive frequency practically does not affect the leakage $1-F_{M\tilde{M}}$, but affects significantly the infidelity contribution $1-F_{\tilde{M}U}$ due to the imperfection of the unitary operation. This is exactly what is expected from the analysis in Sec. VI, since the unitary imperfection is due to imperfect resonance between the drive and the target qubit. As we see from the numerical results, the drive frequency $\omega_{\rm d}=(\omega_{\rm t}^{\rm c0}+\omega_{\rm t}^{\rm c1})/2$ typically gives a better optimized fidelity than for $\omega_{\rm d}=\omega_{\rm t}^{\rm c0}$ or $\omega_{\rm d}=\omega_{\rm t}^{\rm c1}$.

In all previous numerical plots, we assumed the qubit-qubit coupling $g/2\pi=3$ MHz. Figure 20 shows the CNOT gate infidelity for $g/2\pi=1.5$, 3, and 6 MHz, while $\Delta/2\pi=190$ MHz and other parameters are as in Fig. 18. As expected, for minima of the lines in Fig. 20, the CNOT gate duration decreases with increasing $g$ as $\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}\propto g^{-1}$; however, we see that the infidelity increases crudely as $g^{2}$ (consistent with the scaling $\omega_{zz}\propto g^{2}$, so that $\omega_{zz}\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}\propto g$). Increase of $g$ also increases undesired $zz$-interaction of idling qubits Galiautdinov2012 and therefore cannot be used as a simple way to reduce the CNOT gate duration in an experiment.

In Fig. 21 we numerically analyze dependence on the relative duration of the pulse ramp by changing the ratio $\tau_{\rm r}/\tau_{\rm p}$ in the pulse shape (48): $\tau_{\rm r}/\tau_{\rm p}=0.1$, 0.2, 0.3, 0.4, and 0.5 (other parameters are as in Fig. 18 with $\Delta/2\pi=190$ MHz.) For clarity we do not show all oscillations on the left, cutting the lines at some maxima of the oscillations. We see that at the right side of the graph (long $\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}$) it is better to have the longest possible ramp, $\tau_{\rm r}/\tau_{\rm p}=0.5$. However, in the optimal range of short $\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}$ with still small infidelity, the line with $\tau_{\rm r}/\tau_{\rm p}=0.3$ shows the best performance, which can be understood as the following trade-off. For longer ramps and the same $\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}$, we need to use larger $\varepsilon_{\rm m}$ that increases the leakage (even though the ramp is smoother), while for shorter ramps and the same $\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}$, the leakage is also increased because the ramp is too short and consequently non-adiabatic (even though $\varepsilon_{\rm m}$ is smaller). Thus, the ramps should be sufficiently smooth, but it is still beneficial to have a flat part of the pulse.

Finally, let us discuss the effect of the microwave crosstalk, which is always present in experiments Chow2011 ; Chow2012 ; Sheldon2016 . To include it, we need to add the crosstalk Hamiltonian

[TABLE]

which describes the microwave field applied directly to the target qubit. The crosstalk coefficient $c_{\rm ct}$ can in general be complex, and experimental results seem to indicate a complex $c_{\rm ct}$ Sheldon2016 ; Magesan2018 . However, for simplicity here we assume a real $c_{\rm ct}$. In the ideal and semi-analytical theory, the crosstalk with real $c_{\rm ct}$ does not affect the results, except adding the phase $\varphi_{\rm ct}=\int_{0}^{\tau{\rm{}_{p}}}2c_{\rm ct\,}\varepsilon(t)\,dt$ to both $\varphi_{0}$ and $\varphi_{1}$. However, it affects the numerical results because the crosstalk changes the Bloch-sphere angle of tilt caused by the effect of $\omega_{zz}$.

Figure 22 shows numerical results for several values of the (real) crosstalk coefficient $c_{\rm ct}$ (the parameters are as in Fig. 18 with $\Delta/2\pi=190$ MHz). We see that the crosstalk improves fidelity. This improvement can be easily understood. The crosstalk does not affect leakage, but it decreases the imperfect-unitary contribution $1-F_{M\tilde{M}}$ because now the drive detuning $\omega_{\rm d}-\omega_{\rm t}^{\rm c1}$ should be compared with a larger Rabi frequency $\tilde{\varepsilon}_{1}+c_{\rm ct}\varepsilon$ instead of $\tilde{\varepsilon}_{1}$, and therefore the resulting tilt of the Bloch-sphere evolution is smaller (the same effect for the control-qubit state $|0\rangle$ if the signs of $\tilde{\varepsilon}_{0}$ and $c_{\rm ct}\varepsilon$ coincide). Thus, somewhat unexpectedly, the microwave crosstalk can improve the gate fidelity. However, this improvement is significant only if the effect of $\omega_{zz}$ is significant. For example, if in an experiment the infidelity is mainly determined by leakage and decoherence, then the crosstalk will not improve fidelity. Also, we remind that here we considered only real crosstalk coefficients, and the results for a complex $c_{\rm ct}$ may be significantly different.

The effect of the echo sequence is analyzed numerically in the Appendix. The main observation is that the echo-CR gate is typically longer than the basic CR gate for the same level of infidelity (even assuming instantaneous single-qubit gates).

VIII Conclusion

In this paper, we have analyzed analytically, semi-analytically, and numerically the operation of the Cross-Resonance gate for superconducting qubits, focusing on using the CR gate to realize the CNOT operation Chow2011 ; Corcoles2015 ; Sheldon2016 . Our model has been based on the Hamiltonian (1)–(8) and simple pulse shape (no echo sequence).

The analytical theory gives Eq. (18) for the effective drive amplitude $\tilde{\varepsilon}_{n}$ of the target qubit, which depends on the control-qubit state $|n\rangle$. However, the analytics can be used only for sufficiently small (and not too small) physical drive amplitude $\varepsilon$. The next-order analytics (Sec. III.2) slightly widens the applicability range; however, it is still inapplicable in the practically interesting range of $\varepsilon$.

We have found that the speed of the CR gate and the compensating single-qubit rotations can be obtained very accurately by a sufficiently simple semi-analytical theory, discussed in Sec. III.3. This theory [Eqs. (26)–(29)] is based on solving a one-qubit time-independent Schrödinger equation. The semi-analytical theory depends on only two parameters: dimensionless qubit-qubit detuning $\Delta/\eta_{\rm c}$ (normalization is the control-qubit anharmonicity $\eta_{\rm c}$) and dimensionless drive amplitude $\varepsilon/\eta_{\rm c}$. The dimensionless speed of the CR gate as a function of these two parameters is shown in Fig. 7. The CR gate speed cannot be increased indefinitely by increasing the drive amplitude $\varepsilon$; the maximum speed (as a function of $\Delta/\eta_{\rm c}$) and the corresponding drive amplitude are shown in Fig. 8. As follows from Figs. 7 and 8, the best operation of the CR gate is expected for the detuning within the range $0.5\,\eta_{\rm c}<\Delta<\eta_{\rm c}$.

Full numerical approach (discussed in Sec. IV) is mainly needed to calculate intrinsic fidelity of the CNOT-equivalent CR gate. The numerical results depend on the pulse shape, for which we use the simple cosine-ramp model, Eq. (48). Most importantly, the infidelity $1-F_{MU}$ has a minimum as a function of the mid-pulse drive amplitude $\varepsilon_{\rm m}$. The minimum value depends on the detuning $\Delta$ (Figs. 14 and 18), and for typical parameters used in this paper the optimal infidelity is on the order of $10^{-3}$ (though it approaches $10^{-4}$ for some parameters, and in principle, the theoretical infidelity can be arbitrarily small for complicated pulse shapes Kirchhoff2018 ).

Our model does not include decoherence, so the error budget of the gate consists of two contributions: due to imperfect unitary operations and due to leakage [Eq. (57)]. The imperfect unitary dominates at small $\varepsilon_{\rm m}$. The mechanism of this imperfection is related to the $zz$-interaction of the qubits [Eq. (10)], which makes the target-qubit frequency dependent on the control-qubit state ($\omega_{\rm t}^{\rm c1}-\omega_{\rm t}^{\rm c0}=\omega_{zz}$), and therefore makes it impossible to use a drive frequency exactly on-resonance with the target qubit in both cases (for the control-qubit states $|0\rangle$ and $|1\rangle$). Because of this contribution into the error budget, at small $\varepsilon_{\rm m}$ the gate infidelity is very sensitive to small changes of the drive frequency (Fig. 19), and the microwave crosstalk can improve fidelity (Fig. 22).

The other contribution to the error budget, which dominates at large drive amplitudes $\varepsilon_{\rm m}$, is due to leakages. The main leakage is in the strongly-driven control qubit and is caused by non-adiabaticity during the ramps of the pulse. Depending on the ratio $\Delta/\eta_{\rm c}$ between the detuning and control-qubit anharmonicity, the main leakage mechanism can be either $|0\rangle\leftrightarrow|1\rangle$ (between the control-qubit states $|0\rangle$ and $|1\rangle$) or $|0\rangle\to|2\rangle$ or some other leakage channel. In particular, for $|\Delta|\ll\eta_{\rm c}$ there is a strong leakage $|0\rangle\leftrightarrow|1\rangle$ because the drive frequency is near resonance with the control qubit. Similarly, for $|\Delta-\eta_{\rm c}/2|\ll\eta_{\rm c}$ there is a strong leakage $|0\rangle\to|2\rangle$ because in this case the energy difference between states $|2\rangle$ and $|0\rangle$ of the control qubit is near resonance with doubled frequency of the drive. In some cases we also observed a significant leakage for the channel $\overline{|0,1\rangle}\to\overline{|1,2\rangle}$ (when it becomes near-resonance with the doubled drive frequency).

Strong leakage makes the CR gate operation impractical when detuning $\Delta$ is close to [math], $\eta_{\rm c}/2$, $\eta_{\rm c}$, $3\eta_{\rm c}/2$, etc. Therefore, we would expect the best operation of the CR gate for the detuning within the range $0.6\,\eta_{\rm c}<\Delta<0.8\,\eta_{\rm c}$ (see Figs. 7 and 8). Another reasonable range (which requires smaller drive amplitudes) is $0.2\,\eta_{\rm c}<\Delta<0.3\,\eta_{\rm c}$.

Crudely, the CNOT-equivalent gate duration $\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}$ (optimized over $\varepsilon_{\rm m}$) is comparable to $\pi/g$, with a coefficient somewhat larger or smaller than 1, depending on $\Delta/\eta_{\rm c}$ – see Fig. 18. The optimized intrinsic infidelity for a simple pulse shape (48) is crudely comparable to $(g/\eta_{\rm c})^{2}$ (as follows from scaling of $\omega_{zz}$ and $\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}$), with a coefficient typically about $10^{0}$–$10^{1.5}$, depending on the drive frequency, $\Delta/\eta_{\rm c}$, crosstalk, etc. We have found that to reduce the non-adiabatic leakage, the ramps of the pulse should occupy a significant fraction of the pulse duration; however, the flat part in the middle of the pulse should not be shortened to zero.

In this paper we have not focused on analyzing the echo sequence. However, some numerical results for the echo-CR gate are presented in the Appendix. In the echo sequence, there are four ramps instead of two; therefore, for the same total pulse duration, the ramps are shorter. This increases non-adiabatic leakages, so for the same infidelity, the echo-CR gate duration is typically longer than for the basic CR gate (Fig. 24).

A crude estimate of the CR gate infidelity contribution due to decoherence (within the computational subspace only) is given by Eq. (51). However, for large drive amplitudes, the bare state $|2\rangle$ and higher states of the control qubit are significantly occupied; their decoherence is typically much faster than in the computational subspace and therefore can significantly increase the gate infidelity. Another potential mechanism for experimental CR gate infidelity is due to two-level systems (TLSs) produced by impurities. A strong drive changes the effective energy levels in the control qubit by about one hundred MHz, and therefore the TLSs can become on-resonance with the control qubit during the ramp of the pulse. Moreover, large drive amplitude makes multi-photon processes easily possible, and therefore TLSs can become resonant with combinational frequencies for many channels (similar to the situation for fast measurement of a qubit).

Note that in this paper we did not explicitly take into account the resonator, providing coupling between the qubits; instead we replaced it with an equivalent direct coupling. It would be interesting to repeat our numerical simulations, taking into account the resonator levels explicitly. However, we do not expect a significant modification of the results because the additional non-adiabatic effects should be suppressed by typically large detuning between the resonator and qubits. It would also be interesting to include decoherence into the simulations; however, it is not obvious what a proper model is for decoherence of higher levels, relevant to actual experimental situations.

We hope that some experimental group will carry out detailed measurements of the CNOT gate duration and error budget of the CR gate as functions of the drive amplitude, drive frequency, detuning, and pulse shape. It will be interesting (and important for the CR gate application in quantum computing) to compare experimental results with our theoretical findings.

Acknowledgements.

The authors thank Juan Atalaya, Leonid Pryadko, Vinay Ramasesh, Machiel Blok, Ravi Naik, Irfan Siddiqi, Matthew Ware, and Britton Plourde for useful discussions. The work was supported by ARO grant No. W911NF-18-1-0178.

Appendix A Effect of the echo sequence

In the main text we focused on analysis of the basic CR gate. Here we analyze the echo-CR gate Sheldon2016 ; Corcoles2013 ; Takita2016 .

The echo-CR gate Corcoles2013 is essentially a sequence of two same-shape basic CR gates with halved rotation angles, $|\varphi_{1}-\varphi_{0}|=\pi/2$, with the control qubit state flipped ($|0\rangle\leftrightarrow|1\rangle$) in between the two halves of the procedure and with the flipped phase ($\varepsilon\to-\varepsilon$) of the applied microwave drive for the second half of the procedure. In this case the $x$-rotation angles for the target qubit become $\varphi_{1}=\varphi_{1}^{(1)}+\varphi_{0}^{(2)}=\varphi_{1}^{(1)}-\varphi_{0}^{(1)}$ and $\varphi_{0}=\varphi_{0}^{(1)}+\varphi_{1}^{(2)}=\varphi_{0}^{(1)}-\varphi_{1}^{(1)}$, where the superscripts refer to the first or second half of the procedure. Consequently, $\varphi_{1}=-\varphi_{0}$ and also $\varphi_{1}-\varphi_{0}=\pi$ when $\varphi_{1}^{(1)}-\varphi_{0}^{(1)}=\pi/2$. This is what is often called $ZX_{\pi/2}$ gate Chow2011 ; Corcoles2013 ; Sheldon2016 . Note that the relations $\varphi_{0}^{(2)}=-\varphi_{0}^{(1)}$ and $\varphi_{1}^{(2)}=-\varphi_{1}^{(1)}$ are because of the symmetry of the procedure and the phase shift by $\pi$ for the drive, $\varepsilon(t+\tau_{\rm p}/2)=-\varepsilon(t)$, where $\tau_{\rm p}/2$ is the time difference between the two halves of the procedure. Also note that the control qubit should be flipped back after the second pulse (though this flip can sometimes be compiled into the overall sequence of an algorithm).

Because of the symmetry, the echo sequence eliminates the need to apply the compensating $x$-rotation of the target qubit (by $-\varphi_{0}$, as assumed in the main text). It also eliminates the need to apply compensating $z$-rotation of the control qubit (by $\theta_{0}-\theta_{1}+\pi/2$, as assumed in the main text). This significantly reduces experimental complexity. Nevertheless, if we want to produce CNOT gate from the echo-CR gate $ZX_{\pi/2}$, we still need to apply additional $x$-rotation by $\pi/2$ for the target qubit and $z$-rotation by $\pi/2$ for the control qubit.

The analytical (Sec. III.1) and semi-analytical theory (Sec. III.3) for the echo-CR gate does not change compared with the basic CR gate because of the symmetry: we can just use the total gate duration with the same drive amplitude. In particular, Figs. 7 and 8 are still applicable without any change. However, numerical results are different because the ramps for the echo-CR gate of the same duration are shorter and correspondingly the leakage is typically bigger.

In the numerical simulations, for the first half of the procedure ($0\leq t\leq\tau_{\rm p}/2$) we use the pulse shape (48) with substitutions $\tau_{\rm p}\to\tau_{\rm p}/2$ and $\tau_{\rm r}\to\tau_{\rm r}/2$, while for the second half ($\tau_{\rm p}/2\leq t\leq\tau_{\rm p}$) we use the inverted shape, $\varepsilon(t)=-\varepsilon(t-\tau_{\rm p}/2)$. Therefore, $\tau_{\rm p}$ is the total pulse duration and $\tau_{\rm r}$ is the total duration of the two front ramps (or two rear ramps). We assume ideal instantaneous $\pi$-rotations of the control qubit (about $x$ axis) at time moments $\tau_{\rm p}/2$ and $\tau_{\rm p}$, and also ideal instantaneous rotations converting the $ZX_{\pi/2}$ gate into CNOT: $x$-rotation by $\pi/2$ for the target qubit and $z$-rotation by $\pi/2$ for the control qubit. The fidelity is calculated using Eq. (42), which compares the actual gate with the ideal $ZX_{\pi/2}$ gate, i.e., Eq. (43) with $\theta_{0}-\theta_{1}=0$ and $|\varphi_{1}|=\pi/2$ (the relation $|\varphi_{1}-\varphi_{0}|=\pi$ is achieved by varying $\tau_{\rm p}$).

Figure 23 shows a comparison between the CR gate performances with and without echo sequence. The dashed lines (without echo) are the same as lines in Fig. 18, except now we use the drive frequency $\omega_{\rm d}=(\omega_{\rm t}^{\rm c0}+\omega_{\rm t}^{\rm c1})/2$. The solid lines show the results for the echo-CR gate with the same parameters (colors correspond to particular detunings $\Delta$). Most importantly, we see that the steep increase at the left for the solid lines occurs at larger $\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}$. This means that for the same infidelity $1-F_{MU}$, the echo-CR gate is longer than the basic CR gate (even not including the durations of the additional $\pi$-pulses). As mentioned above, this is because the leakage is a more severe problem for the echo-CR gate: four ramps instead of two make their durations shorter, and this significantly increases nonadiabaticity during ramps.

Figure 24 summarizes our semi-analytical and numerical results for the basic CR gate and the echo-CR gate. It shows the duration of the CNOT-equivalent gate $\tau_{\rm p}^{\rm\scriptscriptstyle CNOT}$ versus the detuning $\Delta$ between the control and target qubits, using both the dimensionless and dimensional scales. We use the coupling $g/2\pi=3$ MHz, anharmonicity $\eta_{\rm c}=\eta_{\rm t}=300$ MHz, relative ramp duration $\tau_{\rm r}/\tau_{\rm p}=0.3$, and the drive frequency $\omega_{\rm d}=(\omega_{\rm t}^{\rm c0}+\omega_{\rm t}^{\rm c1})/2$. For different values of $g$, $\eta_{\rm c}$ and $\eta_{\rm t}$, the numerical results on the dimensionless scale are not expected to change significantly (semi-analytical results would remain exactly the same).

The dashed line in Fig.  24 shows the CNOT gate duration (minimized over the drive amplitude $\varepsilon$) as follows from the semi-analytical theory with a rectangular pulse shape. This line is the same as the solid line in Fig. 8 on the inverted scale (duration instead of speed). The solid line in Fig.  24 shows the optimized CNOT gate duration for the smooth pulse shape (48) with $\tau_{\rm r}/\tau_{\rm p}=0.3$, also obtained semi-analytically. We see that this duration is longer than for the rectangular pulse, but not by the naive factor $1/(1-\tau_{\rm r}/\tau_{\rm p})=1.43$ (the ratio is significantly less than this factor because of the pulse-shape integration over the lines in Fig. 7, for which the region near maximum is most important). The semi-analytical solid and dashed lines are the same for the echo-CR and basic CR gates.

The symbols in Fig.  24 show numerical durations of the CNOT-equivalent gates for three levels of the infidelity $1-F_{MU}$: 0.3% (crosses), 1% (triangles), and 3% (circles), using the pulse shape (48) with $\tau_{\rm r}/\tau_{\rm p}=0.3$. Green (thicker) symbols are for the echo-CR gate, while orange (thinner) symbols are for the basic CR gate. The symbols are presented for the detunings $\Delta/2\pi$ of $-200$, $-170$, $-130$, $-70$, $-40$, $40$, $70$, $100$, $130$, $170$, $190$, $210$, $230$, and $250$ MHz. As expected, all symbols are above the solid line. For the detunings of $-130$ and 100 MHz, the symbols are quite close to the solid line. This means that at these detunings, the leakage is not yet too strong for the near-optimal values of the drive amplitudes (e.g., compare lines in Figs. 10 and 14 for $-70$ MHz); the reason for a weak leakage is discussed later. In contrast, for some detunings, the symbols in Fig.  24 are much above the solid line. For example, for the detuning of 170 MHz, the semi-analytics predicts the duration of 70 ns, while for 1% infidelity, the numerical results give 115 ns for the basic CR gate and 141 ns for the echo-CR gate. Such a big difference indicates a very significant leakage for the near-optimum drive amplitudes. As expected, in this case the echo-CR gate requires a significantly longer duration than the basic CR gate for the same level of infidelity.

Figure 24 can be used to estimate the range of best detunings, which provide close-to-shortest CNOT gate durations. For the basic CR gate with 1% infidelity, the best detuning range is crudely $0.6<\Delta/\eta_{\rm c}<0.65$, which for our parameters provides the fastest CNOT gate duration of around 90 ns (excluding single-qubit pulses). For the echo-CR gate with 1% infidelity, the best range is the same, and it provides the fastest duration of about 110 ns. Another reasonably well performing range of detunings is around $\Delta/\eta_{\rm c}\simeq 0.25$; for our parameters it gives CNOT gate duration of about 120 ns for the basic CR gate and 130 ns for the echo-CR gate (for 1% infidelity). One more reasonable range is around $\Delta/\eta_{\rm c}\simeq-0.25$; the corresponding CNOT gate durations are 130 ns (without echo) and 150 ns (with echo).

The CNOT gate duration is typically longer if we require a smaller infidelity. Correspondingly, the crosses in Fig. 24 typically are significantly higher than triangles or circles. Moreover, in some cases (e.g., no echo, $250$ MHz and $-200$ MHz) there are no crosses because the infidelity level of 0.3% is never reached. However, for some detunings (e.g., $100$ MHz and $-130$ MHz) the duration is almost the same for the three considered levels of infidelity; as mentioned above, this indicates low leakage for relatively large drive amplitudes. There is even one weird case (no echo, $70$ MHz), where the order of the symbols is reversed (this is because of the unusual behavior of the dashed brown line in Fig. 23 at the left, which relates to the increasing behavior of the orange line in Fig. 10).

As mentioned above, the difference between the symbols and the solid line in Fig. 24 is mainly determined by the leakage. We have checked that the main leakage channel for the detuning within the range $1/3\lesssim\Delta/\eta_{\rm c}\lesssim 2/3$ is the transition $|0\rangle\to|2\rangle$ in the control qubit. This leakage channel is very strong when $\Delta/\eta_{\rm c}$ is close to $0.5$ (because of the resonance in the rotating frame), thus making impossible the practical operation of the CR gate at $\Delta/\eta_{\rm c}\simeq 0.5$. The leakage $|0\rangle\to|2\rangle$ becomes weaker for detunings farther away from this resonance point. For $\Delta/\eta_{\rm c}\gtrsim 2/3$, the leakage channel $|1\rangle\to|2\rangle$ in the control qubit becomes more important (exact resonance at $\Delta/\eta_{\rm c}=1$). Similarly, for $\Delta/\eta_{\rm c}\lesssim 1/3$, the leakage channel $|0\rangle\leftrightarrow|1\rangle$ in the control qubit becomes more important (exact resonance at $\Delta/\eta_{\rm c}=0$). Thus, the leakage is relatively low for the detuning $\Delta$ near $(1/3)\eta_{\rm c}$ or $(2/3)\eta_{\rm c}$. The trade-off between the lower leakage and shorter semi-analytical durations determines the best detuning ranges in Fig. 24.

Note that besides the leakage channels $|0\rangle\to|2\rangle$, $|1\rangle\to|2\rangle$, and $|0\rangle\leftrightarrow|1\rangle$ for the control qubit, there are also important leakage channels, which involve level $|2\rangle$ of the target qubit. For example, for detunings of 210 and 230 MHz, the leakage is dominated by a near-resonance between levels $|01\rangle$ and $|12\rangle$. Also, for detunings of $-200$, $-170$, and $-130$ MHz, the main leakage is due to a near-resonance between levels $|11\rangle$ and $|02\rangle$. Overall, the interplay between different leakage channels is rather complicated, leading to a rather complicated behavior of numerical results in Fig. 24.

Since the echo-CR gate is more affected by the leakage than the basic CR gate, a natural hypothesis is that its operation can be improved by using smoother ramps, in particular, by changing the relative duration of the front (and rear) ramps from the value $\tau_{\rm r}/\tau_{\rm p}=0.3$ (used in Fig. 24) to the maximum possible value $\tau_{\rm r}/\tau_{\rm p}=0.5$. To check this hypothesis, we have also simulated the echo-CR gate operation with $\tau_{\rm r}/\tau_{\rm p}=0.5$, but the results were inconclusive: sometimes this makes CNOT time slightly shorter, sometimes slightly longer. For example, for the detuning of 190 MHz (our shortest-duration point), changing $\tau_{\rm r}/\tau_{\rm p}$ from 0.3 to 0.5 increases the CNOT time by 4 ns for 0.3% infidelity and by 1 ns for 1% infidelity, but decreases it by 12 ns for 3% infidelity. So, crudely, we think that the pulse shape ratio $\tau_{\rm r}/\tau_{\rm p}=0.3$ is still reasonable for the echo-CR gate.