Simple framework for systematic high-fidelity gate operations

Assuming the time dynamics of the targeted gate is known for each time $t\in[0,t_g]$, we can now equally define the erroneous dynamics based on equation (4) as

$$\begin{align} \mathcal{E}\left(t\right) = U^{\dagger}_\text{ideal}\left(t\right)U\left(t\right). \end{align} \tag{ 6 }$$

We can find the corresponding ideal and error Hamiltonian by plugging $\mathcal{E}(t)$ into equation (1) to arrive at [49, 55]

$$\begin{align} i\hbar\dot{\mathcal{E}}\left(t\right)& = \left[U_\text{ideal}^\dagger\left(t\right)HU_\text{ideal}\left(t\right)-i\hbar U_\text{ideal}^\dagger\left(t\right)\dot{U}_\text{ideal}\left(t\right)\right]\mathcal{E}\left(t\right). \end{align} \tag{ 7 }$$

Additionally, U_ideal and H_ideal are related via the time-ordered exponential

$$\begin{align} U_\text{ideal}\left(t\right) = \mathcal{T}\exp\left(-\frac{i}{\hbar}\int_0^{t} H_\text{ideal}\left(t^{\prime}\right) \mathrm{d}t^{\prime}\right). \end{align} \tag{ 8 }$$

Consequently, the associated ideal and error Hamiltonian are then given by [56]

$$\begin{align} H_\text{ideal}\left(t\right) & = i\hbar \dot{U}_\text{ideal}\left(t\right)U_\text{ideal}^\dagger\left(t\right), \end{align} \tag{ 9 }$$

$$\begin{align} H_\text{error}\left(t\right)& = U_\text{ideal}^\dagger\left(t\right)\left[H\left(t\right)-H_\text{ideal}\left(t\right)\right]U_\text{ideal}\left(t\right). \end{align} \tag{ 10 }$$

The formal solution to equation (7) at $t = t_g$ is again the time-ordered exponential given by

$$\begin{align} \mathcal{E} & = \mathcal{T}\exp\left(-\frac{i}{\hbar}\int_0^{t_g} H_\text{error}\left(t\right) \mathrm{d}t\right) \end{align} \tag{ 11 }$$

$$\begin{align} & = \exp\left(-\frac{i}{\hbar}\sum_{n = 1}^\infty \bar{H}_n\right) \end{align} \tag{ 12 }$$

$$\begin{align} &\approx \exp\left(-\frac{i}{\hbar}\int_0^{t_g} H_\text{error}\left(t\right) \mathrm{d}t\right). \end{align} \tag{ 13 }$$

From the first to the second line we applied the Magnus expansion with Magnus coefficients $\bar{H}_n$. We further assume small errors $\underset{t\in [0,t_g]}{\text{max}}||H_\text{error}(t)||_2\ll \hbar \pi/t_g$, where $||\cdot||_2$ denotes the spectral matrix norm, (otherwise we choose a closer desired gate) to ensure a fast converging Magnus series. In the last step, we truncated the Magnus expansion at lowest order. We show later that this order is sufficient to find parameters for quantum operations with gate fidelities in the order of $1-F\lt10^{-4}$ for three important applications. Since we assumed small errors, we can also expand the matrix exponential equation (13) up to linear order

$$\begin{align} \mathcal{E}\approx \unicode{x1D7D9}+\left(-\frac{i}{\hbar}\int_0^{t_g} H_\text{error}\left(t\right) \mathrm{d}t\right). \end{align} \tag{ 14 }$$

As a matter of fact, our error matrix corresponds to the Hamiltonian errors of the error generator [57] defined as $\text{log}(\chi)$, where χ is a superoperator that can be measured with quantum process tomography techniques. Furthermore, an analogous derivation can be performed if incoherent errors are included by replacing the system Hamiltonian with a stochastic Hamiltonian [58] or Liouvillian. We leave this to a future investigation and focus in this work on coherent errors.

To achieve a high-fidelity quantum gate we rewrite equation (14) into a minimization condition

$$\begin{align} \text{min}||\mathcal{E}-\unicode{x1D7D9}||_2 = \text{min}\bigg|\bigg|\frac{1}{\hbar}\int_0^{t_g} H_\text{error}\left(t^{\prime}\right) \mathrm{d}t^{\prime}\bigg|\bigg|_2. \end{align} \tag{ 15 }$$

There are multiple approaches to finding the minimum, such as brute force numerical Nelder-Mead [59], GRAPE [39] and CRAB algorithms [40] or a parameterization of the integration path [38, 45]. In general, there are n² free real parameters of a Hermitian matrix (or unitary matrix) that need to be optimized.

In this work, we use a different approach that makes use of two properties that are common in many quantum computing device architectures, sparse interactions and a priori knowledge of the system, that allow us to greatly increase the efficiency of finding optimized pulse shapes. Instead, we use the equivalence of matrix norms to rewrite equation (15) as

$$\begin{align} \text{min}\sum_{k = 1}^{n^2}|\text{tr}\left(\mathcal{E} O_k -O_k\right)|^2, \end{align} \tag{ 16 }$$

where the error channels O_k with $k = 1,\ldots, n^2$ describe a full basis set of n × n operators. A smart choice of the channels $\lbrace O_k\rbrace$ can often allow us to truncate the series after a few terms with minimal consequences. The error channels O_k can be constructed from the a priori knowledge of the system, either through fast tomographic methods [60] or through a trustworthy theoretical model. We focus in this paper on the latter.

We find the choice $O_k\equiv O_{k_1,k_2} = \vert k_1(0)\rangle \langle k_2(0)\vert$ with $k_{1,2} = 1,\ldots, n$ to work well, where $\vert k_{1,2}(t)\rangle$ are eigenstates of a dynamic invariant I(t) with respect to the ideal dynamics U_ideal. The dynamic invariant is a Hermitian operator I(t) that satisfies $i\hbar \dot{I}(t) = [H_\text{ideal}(t),I(t)]$ with $[H_\text{ideal}(0),I(0)] = 0 = [H_\text{ideal}(t_g),I(t_g)]$, where $[A,B] = AB-BA$ [61]. We note that $[H_\text{ideal}(t_g),I(t_g)] = 0$ is not required in general, but guarantees state transfers without final excitation. This allows us to conveniently decompose the ideal dynamics as [62]

$$\begin{align} U_\text{ideal}\left(t\right) = \sum_{m = 1}^n \mathrm{e}^{i\alpha_{m}\left(t\right)}\left\vert m\left(t\right)\right\rangle \left\langle m\left(0\right)\right\vert , \end{align} \tag{ 17 }$$

where the phase α_m is the Lewis–Riesenfeld phase

$$\begin{align} \alpha_m\left(t\right) = \frac{1}{\hbar}\int_0^t \left\langle m\left(t^{\prime}\right)\vert i\hbar\frac{\mathrm{d}}{\mathrm{d}t^{\prime}}-H_\text{ideal}\left(t^{\prime}\right)\vert m\left(t^{\prime}\right)\right\rangle \mathrm{d}t^{\prime}. \end{align} \tag{ 18 }$$

While finding the dynamic invariant I(t) in the general case is hard, there are two special cases of interest that allow for a great simplification. If U_ideal describes an adiabatic dynamic or if $[H_\text{ideal}(t_1),H_\text{ideal}(t_2)] = 0$ for all $t_1,t_2\in [0,t_g]$, one can find a common set of eigenstates for $H_\text{ideal}(t)$ and I(t). In the latter case, the eigenstates are time-independent, $\vert m(t)\rangle = \vert m(0)\rangle$, and also eigenstates of $U_\text{ideal}(t)$. Therefore, the phase can be simplified to $\alpha_{m}(t) = -\frac{1}{\hbar} \int_0^{t} \epsilon_m \mathrm{d}t^{\prime}$, where $\epsilon_m(t)$ is the mth eigenenergy of H_ideal. Note that all practical applications discussed in section 3 fall in the latter case. For the adiabatic case, one has also to add the geometric phase [56].

Using the dynamic invariant decomposition, we can now rewrite the error Hamiltonian into a complex parameter $g_{k_1,k_2}$, the phase $\alpha_{k_1,k_2}$ [62], and time-independent matrix elements $\vert k_1(0)\rangle \langle k_2(0)\vert$

$$\begin{align} H_{\text{error}}\left(t\right) = \sum_{k_1,k_2} g_{k_1,k_2}\left(t\right) \mathrm{e}^{i \alpha_{k_1,k_2}\left(t\right)}\left\vert k_1\left(0\right)\right\rangle \left\langle k_2\left(0\right)\right\vert . \end{align} \tag{ 19 }$$

Here the complex parameter

$$\begin{align} g_{k_1,k_2}\left(t\right) = \left\langle k_1\left(t\right)\right\vert H\left(t\right)-H_\text{ideal}\left(t\right)\left\vert k_2\left(t\right)\right\rangle \end{align} \tag{ 20 }$$

describes the transition matrix element caused by the erroneous dynamics. Similarly, we find $\alpha_{k_1,k_2}(t)$ to be the difference of the associated Lewis–Riesenfeld phases

$$\begin{align} \alpha_{k_1,k_2}\left(t\right) = \alpha_{k_2}\left(t\right)-\alpha_{k_1}\left(t\right). \end{align} \tag{ 21 }$$

In summary, the decomposition of the error Hamiltonian into operators $O_k\equiv O_{k_1,k_2} = \vert k_1(0)\rangle \langle k_2(0)\vert$ allows us to 'measure' the deviations from the ideal dynamics and quantifies the probability of making a coherent error generated by $\mathcal{E}$ [57]. The associated error rate is then given by

$$\begin{align} |\text{tr}\left(\mathcal{E} O_k -O_k\right)|^2 = \left|\frac{1}{\hbar} \int_0^{t_g} g_k\left(t\right) e^{i \alpha_k\left(t\right)} \mathrm{d}t \right|^2 \end{align} \tag{ 22 }$$

with $\alpha_k(t)\equiv \alpha_{k_1,k_2}(t)$ and $g_k(t)\equiv g_{k_1,k_2}(t)$. In the language of pulse optimization, $g_{k}(t)$ is the control signal. For later convenience, we also define a corresponding frequency

$$\begin{align} f_{k}\left(t\right) = \frac{\mathrm{d}}{\mathrm{d}t}\alpha_{k}\left(t\right). \end{align} \tag{ 23 }$$

The error rates can be related to the problem of transmitting a signal through a channel with finite frequency bandwidth. This can be shown by substituting [48, 63]

$$\begin{align} \frac{\mathrm{d}}{\mathrm{d}t}\alpha_{k}\left(t\right)\mathrm{d}t = f_k\left(t\right)\mathrm{d}t \equiv \nu_{f,k} t_g \mathrm{d}s \end{align} \tag{ 24 }$$

with a constant $\nu_{f,k}$ and assuming $f_k(t)\neq 0$. By choosing $s(t_g) = 1$, the static parameter $\nu_{f,k}$ can be seen as the averaged resonance frequency over the time interval $[0,t_g$] of the ideal system, $\nu_{f,k} t_g = \alpha_k(t_g)-\alpha_k(0)$. In certain scenarios it might be beneficial to fix $\nu_{f,k}$ equal to characteristic frequencies such as the idle resonance frequencies instead and allow $s(t_g)\neq 1$. In this manuscript, we use the upper convention. In both cases, the relation between real time t and dilated normal time s is given by integrating equation (24) arriving at [48]

$$\begin{align} t\left(s\right) & = \int_{s\left(0\right) = 0}^{s} \frac{\nu_{f,k} t_g }{f\left(s^{\prime}\right)} \mathrm{d}s^{\prime},\end{align} \tag{ 25 }$$

$$\begin{align} s\left(t\right) & = \int_{0}^{t} \frac{f\left(t^{\prime}\right)}{\nu_{f,k} t_g } \mathrm{d}t^{\prime}. \end{align} \tag{ 26 }$$

The inverted functions are best acquired using numerical interpolation [48], e.g. Mathematica directly provides $\tilde{g}(t)$ using the command $\text{Interpolation}[\text{Table}[\lbrace t[s], \tilde{g}[s]\rbrace, \lbrace s, 0, 1\rbrace]]$ with sufficient sampling.

For the following discussions, we focus on a single error rate, thus dropping the index k to increase readability. We can rewrite the error rate in equation (22) using the substitution (24) as

$$\begin{align} |\text{tr}\left(\mathcal{E} O -O\right)|^2 & = \left|\frac{1}{\hbar} \int_{s\left(0\right)}^{s\left(t_g\right)} \tilde{g}\left(t\left(s\right)\right) \mathrm{e}^{i \nu_f\, t_g s} \mathrm{d}s \right|^2, \end{align} \tag{ 27 }$$

$$\begin{align} & = S\left(\tilde{g}\left(\nu_f\, t_g\right)\right), \end{align} \tag{ 28 }$$

with $\tilde{g}(t(s)) = g(t(s)) \frac{\mathrm{d}t}{\mathrm{d}s}(s)$ and the associated energy spectral density S. From the first to second line, we replaced the integral which corresponds to a short-time Fourier transformation with the expression for an energy spectral power of the input signal $\tilde{g}(s)$ [48]. As a consequence, we have now shifted the task from minimizing the error rates to optimizing the energy spectral density, a task investigated in the field of signal processing, and which has been solved for many input signals. Below, we show a few examples of how signal processing can be used for finding optimized pulse shapes.

2.3.1. Synchronization

2.3.2. Static pulse shaping

We first discuss pulse shaping techniques in the case of a single experimental control parameter, e.g. baseband signals, corresponding to real transition matrix elements g(t) (constant phase factors can be factored out). In this case, fast operations with consistently low error rates can be achieved using window functions w(t) designed for optimized spectral concentration such as the discrete prolate spheroidal sequence (DPSS or Slepian), and the Dolph–Chebyshev window. Alternatively, if faster operations at the cost of larger coherent errors are desired, a Hamming window is the best choice. Unfortunately, the optimized window functions typically require a high computational cost, high bandwidth, and high time-resolution, thus restricting their practical use. Therefore, in practical applications in quantum computing, often approximations of the optimized windows are used, which reach almost equally small errors. An example is the Kaiser window that is often used to replace the DPSS window [65]

$$\begin{align} w\left(t\right) = \mathcal{N} \mathcal{I}_{0}\left(\frac{2\lambda}{t_g} \sqrt{t\left(t_g-t\right)} \right). \end{align} \tag{ 29 }$$

Here, $\mathcal{N}$ is a normalization constant defined via $\int_0^{t_g}w(t)\mathrm{d}t = t_g$ and $\mathcal{I}_{0}$ is the 0th order modified Bessel function. Alternatively, one can also use a Fourier series [48]

$$\begin{align} w\left(t\right) & = w_\text{even}\left(t\right) + w_\text{odd}\left(t\right) \end{align} \tag{ 30 }$$

with the even and odd decomposition

$$\begin{align} w_\text{even}\left(t\right) & = \sum_{n = 1}^{N}\lambda_{\text{even},n} \left[1-\cos\left(\frac{2\pi n t }{t_g}\right)\right], \end{align} \tag{ 31 }$$

$$\begin{align} w_\text{odd}\left(t\right) & = \sum_{n = 1}^{N}\lambda_{\text{odd},n} \left[1-\sin\left(\frac{2\pi n t }{t_g}\right)\right]. \end{align} \tag{ 32 }$$

The optimal Fourier coefficients $\left.\lambda_{\text{even}} = [1.0715,-0.0795,0.0043,0.0037]\right.$ and $\left.\lambda_{\text{odd}} = [0,0,0,0]\right.$ can be estimated from a Fourier expansion of the Slepian window or from direct numerical minimization of the error rate [48]. This approach has the advantage that by using only even components, a smooth pulse shape and $w(0) = w(t_g) = 0$ is guaranteed.

A simple and popular pulse shape is the Hann window (sometime also cosine window) ($\left.\lambda_{\text{even}} = [1,0,0,0]\right.$) or its generalization the Tukey window defined as [32]

$$\begin{align} w\left(t,\lambda\right) = \begin{cases} \frac{1}{2-\lambda}\left[1-\cos\left(\frac{2\pi t }{\lambda t_g}\right)\right] & 0 \leqslant t \leqslant \frac{\lambda t_g}{2} \\ \frac{2}{2-\lambda} & \frac{\lambda t_g}{2} < t< t_g -\frac{\lambda t_g}{2} \\ \frac{1}{2-\lambda}\left[1-\cos\left(\frac{2\pi \left(t_g-t\right) }{\lambda t_g}\right)\right] & t_g -\frac{\lambda t_g}{2}\leqslant t \leqslant t_g \end{cases}. \end{align} \tag{ 33 }$$

The Tukey window consists of two halves of a Hann window interleaved by a constant part and is formally the convolution of the Hann window and a rectangular window. For λ = 1 it reduces to the Hann window. Both the Hann and the Tukey window have the advantage that they can easily be synchronized due to their pronounced oscillatory spectral density. An interesting thought is also combining the Tukey window shape with the Fourier approximation, which we will leave for future investigations.

Once decided on a preferred pulse shape, the optimal pulse design to minimize the error rate is then given by setting

$$\begin{align} \tilde{g}\left(s\right) = \hbar A w\left(t_g s\right). \end{align} \tag{ 34 }$$

The amplitude A has to be estimated from the desired gate operation (9).

Figure 1(c) shows the resulting infidelity of a gate operation for different pulse shapes. As designed, the Kaiser window outperforms the Tukey and Hann windows in terms of performance, but the latter may be simpler to implement.

2.3.3. Dynamic pulse shaping

Next, we discuss the case of two orthogonal experimental control parameters, e.g. microwave amplitude and phase, corresponding to complex matrix transition elements, $g(t) = g_R(t)+i g_I(t)$. In this case, the steps for the upper mitigation strategies have to be simultaneously applied for the real and imaginary part. Again, the trivial case $g_R(t)\propto g_I(t)$ corresponds to a single control parameter and the constant phase factor can be factored out. In the general case, the upper methods may fail since the solutions for the real and imaginary components may be incompatible.

Such non-trivial complex matrix transition elements appear for example in the case of resonantly driven gates, i.e. single-qubit gates for spin qubits [12], controlled rotation gates [46], or simultaneous pulsing of both barrier gates for exchange-only qubits [12]. Within the rotating wave approximation, the phase of the MW signal translates into complex matrix transition elements $\tilde{g}(t(s)) = \tilde{g}_R(t(s))+i \tilde{g}_I(t(s))$. Here, $\tilde{g}_R(t(s))$ is the real and $\tilde{g}_I(t(s))$ the imaginary part corresponding to the I/Q quadrature of the MW signal.

Fortunately, such complex signals can actively be used to significantly reduce the error rates compared to window functions [44] by making full use of the additional degree of freedom. For example, the derivative removal by adiabatic gate (DRAG) [42, 66] and Wah-Wah [67–69] protocols both allow to suppress crosstalk from off-resonant drives beyond what conventional window functions can achieve. These protocols can be visualized by integrating equation (27) by parts with respect to the real part of the signal [70]

$$\begin{align} |\text{tr}\left(\mathcal{E}\,O-O\right)|^2 & = S\left(\frac{\frac{d}{ds}\tilde{g}_R\left(t\left(s\right)\right)}{i \nu_f\, t_g}+i \tilde{g}_I\left(t\left(s\right)\right)\right) \end{align} \tag{ 35 }$$

with boundary conditions $\tilde{g}_R(0) = \tilde{g}_R(t_g)$. We get a complete cancellation of the error rate with

$$\begin{align} \frac{\frac{d}{ds}\tilde{g}_R\left(t\left(s\right)\right)}{\nu_f\, t_g} = \tilde{g}_I\left(t\left(s\right)\right), \end{align} \tag{ 36 }$$

where the pulse shape of $\tilde{g}_R(t(s))$ can individually be optimized using window functions.

Figure 1(d) shows that using dynamic pulse shaping protocols (here DRAG) significantly reduces the error rate. For $\frac{dt}{ds} = \text{const}$ this exactly yields the DRAG condition $g_I\propto\dot{g}_R$. The advantage of this strategy is a strong suppression of the error rate at the cost of additional power consumption, which scales with the number of suppressed transitions [71]. Our framework allows generalizing this powerful method beyond microwave control to all systems with independent control over two orthogonal axes.

We first optimize single-qubit operations. Expected coherent errors for resonantly driven single-qubit gates are crosstalk and related spectator errors. For example, off-resonant driving due to microwave leakage or shared driving gates [73], non-linear driving and their impact such as higher harmonic generation [74], phase shifts and frequency shifts [75], and the impact of counter-rotating driving such as Bloch Siegert shifts [76]. To include non-linear driving effects into our theoretical description, we expand the modulated (effective) magnetic field $\boldsymbol{B}_j(t) = \boldsymbol{B}_{j,0}(t) + \sum_{k\neq 0} \boldsymbol{B}_{j,k}(t) e^{-2\pi i \nu_{D} k t},$ in terms of a Fourier series with respect to the drive frequency ν_D , where the Fourier components $\boldsymbol{B}_{j,0}(t)$ and $\boldsymbol{B}_{j,k}(t)$ are assumed to be slowly varying in the time interval $[0,\nu_{D}^{-1})$.

The Hamiltonian can be significantly simplified under the rotating wave approximation (RWA) where we keep stationary terms and disregard all terms which are modulated with frequency $k\nu_D$ with $|k| = 1,2,3\cdots$. Corrections from violations of the RWA scale with $\Omega_\text{Rabi}/\nu_D\sim 10^{-3}$ and are negligible for typical experimental conditions. However, higher-order corrections can become important for ultra-fast gate operations [77, 78] or driving at comparatively low frequencies [17].

In this section, we focus on the RWA case and leave the general investigation to the future. Corrections beyond the RWA can be included by using directly Hamiltonian (38) instead, or by applying the generalized RWA introduced in [79] and discussed in section 3.3 and appendix B.

We use the following Hamiltonian as a basis for our pulse optimization framework

$$\begin{align} H\left(t\right) = \frac{h}{2} \begin{pmatrix} 2\left(E_z -\nu_D\right)-\frac{\dot{\theta}}{\pi} & \left(B^x_{2,1}-i B^y_{2,1}\right)/2 & \left(B^x_{1,1}-i B^y_{1,1}\right)/2 & 0 \\ \left(B^x_{2,-1}+i B^y_{2,-1}\right)/2 & \Delta E_z & 0 & \left(B^x_{1,1}-i B^y_{1,1}\right)/2 \\ \left(B^x_{1,-1}+i B^y_{1,-1}\right)/2 & 0 & -\Delta E_z & \left(B^x_{2,1}-i B^y_{2,1}\right)/2 \\ 0 & \left(B^x_{1,-1}+i B^y_{1,-1}\right)/2 & \left(B^x_{2,-1}+i B^y_{2,-1}\right)/2 & -2\left(E_z-\nu_D\right)+\frac{\dot{\theta}}{\pi} \\ \end{pmatrix}, \end{align} \tag{ 39 }$$

where we kept time-dependent phases to account for shifts in resonance frequency.

Next we describe our targeted ideal operation [21] as

$$\begin{align} U_\text{ideal} = \mathrm{e}^{-i \pi \int_0^t dt^{\prime} \left[ B^x_{1,1}\left(t^{\prime}\right)S^x_1 + 2\left(\Delta f\left(t^{\prime}\right)-\Delta E_z\left(t^{\prime}\right)\right) S^z_2\right]}\mathrm{e}^{i\theta S^z_2} \end{align} \tag{ 40 }$$

with frequency detuning $\Delta f(t) = E_z(t)+\frac{\Delta E_z(t)}{2}-\nu_D$ and find H_error via equation (10). Our target gate operation describes a single-qubit Rabi oscillation on qubit Q1 and phase shifts on the non-driven qubit Q2. We include the largest expected error, a phase accumulation (S^z ₂) on the non-driven qubit known as ac Stark shift, into our target operation to keep the erroneous evolution small. Such phase-shift can be corrected easily via a virtual z gate on Q2 in an experimental realization.

We find (see appendix C.1) two pairs of dominant error channels described by the operators $O_1 = S_1^{\,z}\mp i S_1^{\,y}$ and $O_2 = S_2^x\pm i S_2^y$ with erroneous evolutions

$$\begin{align} g_1\left(t\right) & = \Delta f\left(t\right)-\frac{\dot{\theta}\left(t\right)}{2\pi}\mp i B^y_{1,1}\left(t\right), \end{align} \tag{ 41 }$$

$$\begin{align} g_2\left(t\right) & = B^x_{2,1}\left(t\right)\pm i B^y_{2,1}\left(t\right), \end{align} \tag{ 42 }$$

desired evolutions

$$\begin{align} f_1\left(t\right) & = \pm \pi B^x_{1,1}\left(t\right), \end{align} \tag{ 43 }$$

$$\begin{align} f_2\left(t\right) & = \pm \pi \left(\Delta f\left(t^{\prime}\right)-\Delta E_z\left(t\right) -\frac{\dot{\theta}}{2\pi}\right) \end{align} \tag{ 44 }$$

and error rates

$$\begin{align} |\text{tr}\left(\mathcal{E}O_1-O_1\right)|^2 & = S\left(\widetilde{g}_1\left[\nu_1 t_g\right]\right), \end{align} \tag{ 45 }$$

$$\begin{align} |\text{tr}\left(\mathcal{E}O_2-O_2\right)|^2 & = S\left(\widetilde{g}_2\left[\nu_2 t_g\right]\right). \end{align} \tag{ 46 }$$

The first error rate describes a shift in the rotation axis (x-direction) of qubit 1 giving rise to S^y ₁ and S^z ₁ errors. The best mitigation strategy is dynamic pulse shaping via a time-dependent phase $\theta(t)$, e.g. through chirping [80, 81]. The second error rate describes a spin-flip of the second qubit due to off-resonant driving, giving rise to S^x ₂ and S^y ₂ errors. The mitigation of the spin-flip errors requires either synchronization or dynamic pulse shaping.

The condition for the synchronization of a $R_{x,y}(\pi)$ gate on qubit j affecting qubit i (i ≠ j) with a rectangular pulse shape is (see also [47])

$$\begin{align} B^x_{i,1} = \frac{2 n+1}{2m}\sqrt{\frac{\left(B^x_{j,1}\right)^2}{4}+\Delta E_z^2} \end{align} \tag{ 47 }$$

with integer m and n. The condition for a synchronized $R_{x,y}(\frac{\pi}{2})$ gate is given by the substitution $n\rightarrow 2n$. For minimal gate time assuming an equally strong global drive $B^x_{i,1} = B^x_{j,1}$, the synchronization condition for a $R_{x}(\frac{\pi}{2})$ gate is simplified to

$$\begin{align} t_g = \frac{\sqrt{16 m^2-1}}{4\Delta E_z}, \end{align} \tag{ 48 }$$

which corresponds exactly to the minima in figure 1(b).

For static pulse shaping, $\theta(t) = 0$ and $B^y_{i,j} = 0$, constant drive frequency $\nu_D = E_z\pm\Delta E_z$, and negligible shift in resonance frequency $\Delta E_z(t) = \Delta E_z$, the first error rate equation (45) vanishes, and we only need to minimize equation (46) which is simplified to

$$\begin{align} |\text{tr}\left(\mathcal{E}O_2-O_2\right)|^2 & = S\left(B^x_{2,1}\left[\Delta E_z t_g\right]\right). \end{align} \tag{ 49 }$$

Optimal pulse shapes for a $R_{x}(\frac{\pi}{2})$ gate are then given by $B^x_{2,1}(t) = \frac{1}{4t_g} w(t)$, where we use the normalized window $\int_0^{t_g} w(t) dt = t_g$. Figure 1(c) displays the simulated infidelity for different pulse shapes as a function of gate time t_g .

Under the same assumptions such as constant drive frequency $\nu_D = E_z\pm\Delta E_z$ and negligible shift in resonance frequency $\Delta E_z(t) = \Delta E_z$, dynamic pulse shaping provides even faster gate times with small errors (see figures 1(c) and (d)). Applying the DRAG method [42] equations (45) and (46) combined with a Hann window, the optimized dynamic pulse shape is (see figure 1(d))

$$\begin{align} B^x_{2,1}\left(t\right) & = \frac{1}{4t_g} w\left(t\right) \left(1-\frac{5}{5+\left(4\Delta E_z t_g\right)^2}\right), \end{align} \tag{ 50 }$$

$$\begin{align} B^y_{2,1}\left(t\right) & = -\frac{\dot{B}^x_{2,1}\left(t\right)}{\Delta E_z}, \end{align} \tag{ 51 }$$

$$\begin{align} \dot{\theta}\left(t\right) & = - \frac{\left(\dot{B}^x_{2,1}\left(t\right)\right)^2}{\Delta E_z}. \end{align} \tag{ 52 }$$

The renormalization of the drive amplitude is due to the additional power in driving. The optimized pulse shape for driving qubit 2 is given by substituting $\Delta E_z\rightarrow -\Delta E_z$ in equations (50)–(52).

Next, we optimize the two-qubit gate. A crucial condition for high-fidelity two-qubit cz gates is an adiabatic turn on/off or pulse of the exchange interaction, which can give rise to substantial errors if violated [32]. While in principle an echo pulse sequence allows suppressing non-adiabatic errors [46] for a cz gate, the echo pulse is often inconvenient and introduces additional noise through the longer gate sequences.

Starting from Hamiltonian (38) we notice that the exchange interaction only affects the odd parity states $\lbrace\vert \uparrow\downarrow\rangle ,\vert \downarrow\uparrow\rangle \rbrace$ (see appendix A). Without loss of generality, this allows us to project the full dynamics on a two-level system, including the global phase of this subspace. Introducing a new set of Pauli operators $\sigma_x = \vert \uparrow\downarrow\rangle \langle \downarrow\uparrow\vert +\vert \downarrow\uparrow\rangle \langle \uparrow\downarrow\vert$, $\sigma_y = -i \vert \uparrow\downarrow\rangle \langle \downarrow\uparrow\vert +i\vert \downarrow\uparrow\rangle \langle \uparrow\downarrow\vert$, and $\sigma_z = \vert \uparrow\downarrow\rangle \langle \uparrow\downarrow\vert -\vert \downarrow\uparrow\rangle \langle \downarrow\uparrow\vert$, we find

$$\begin{align} H\left(t\right) = \frac{h}{2}\left( - J + \Delta E_z \,\sigma_z+J\, \sigma_x\right). \end{align} \tag{ 53 }$$

We define the targeted ideal gate as

$$\begin{align} U_\text{ideal} = \mathrm{e}^{i\pi\int_0^t J\left(t^{\prime}\right) dt^{\prime}} \mathrm{e}^{-i \pi\int_{0}^{t}dt^{\prime} \nu_\text{ST}\left(t^{\prime}\right)\widetilde{\sigma}_z} \end{align} \tag{ 54 }$$

which implements up to single-qubit phases a cz-gate at the time $\int_0^{t_g} J(t^{\prime}) dt^{\prime} = 1/2$. Here, the unitary $U_\text{ST} = e^{-\frac{i}{2} \tan^{-1}\left(\frac{J(t)}{\Delta E_z (t)}\right)\sigma_y}$ diagonalizes Hamiltonian (53) with resonance frequency $\nu_\text{ST}(t) = \sqrt{\Delta E_z(t)^2+J(t)^2}$ and $\widetilde{\sigma}_{x,z} = U_\text{ST}^\dagger\sigma_{x,z}U_\text{ST}$. Our target operation describes the adiabatic phase evolution due to the exchange interaction.

There is (see appendix C.2) a single dominant error channel causing swap-oscillations [28] that is described by the spin-flip operator $O = \vert \widetilde{\uparrow\downarrow}\rangle \langle \widetilde{\downarrow\uparrow}\vert$, where $\vert \widetilde{\uparrow\downarrow}\rangle$ and $\vert \widetilde{\downarrow\uparrow}\rangle$ are the eigenstates of Hamiltonian (53). The erroneous and targeted time evolutions are

$$\begin{align} g\left(t\right) & = -h\frac{\Delta E_z\, \dot{J}-\dot{\Delta E_z} J}{4\pi\nu_\text{ST}^2\left(t\right)}, \end{align} \tag{ 55 }$$

$$\begin{align} f\left(t\right) & = \pi \nu_\text{ST}\left(t\right). \end{align} \tag{ 56 }$$

For small exchange $J(t)\ll \Delta E_z$ and constant Zeeman splitting $\Delta E_z(t)\approx \Delta E_z$ we can simplify $g(t)\propto \dot{J} (t)$ and $\nu_\text{ST}(t)\approx \Delta E_z$ to find the error rate

$$\begin{align} P_{\vert \uparrow\downarrow\rangle \rightarrow\vert \downarrow\uparrow\rangle } \sim S\left(\dot{J}\left[\Delta E_z t_g\right]\right). \end{align} \tag{ 57 }$$

Remarkably, this optimization condition for an adiabatic cz gate is identical to the condition for minimizing single-qubit crosstalk in equation (49) under the replacement $B^x_{2,1}\rightarrow\dot{J}$ with the same invariant $t_g\times|\Delta E_z|$. The conditions for the synchronization of a cz gate with a rectangular pulse shape and minimal time t_g is [28]

$$\begin{align} t_g = \frac{\sqrt{4 m^2-1}}{2\Delta E_z\left(v_B\right)} \end{align} \tag{ 58 }$$

with integer m. Here $\Delta E_z(v_B)$ is the difference in resonance frequency during the pulse. Note that this is equivalent to the synchronization condition of a $R_{x,y}(\pi)$ gate (see equation (48)).

We show later that static pulse shaping is sufficient to get error rates for the cz-gate below $1-F\lt10^{-4}$. Due to the non-linear relation between barrier voltage and exchange interaction (see appendix A) the optimal pulse shape for the barrier voltage pulse shape $v_B(t)$ is then given by (see appendix D)

$$\begin{align} v_B\left(t\right)& = \frac{1}{\alpha}\text{log}\left(\frac{\sqrt{J\left(t\right)/J_\text{sat}}}{|1-J\left(t\right)/J_\text{sat}|}\right), \end{align} \tag{ 59 }$$

$$\begin{align} J\left(s\left(t\right)\right)& = \tan\left[2 A w\left(t_g s\right)+\tan^{-1}\left(\frac{J\left(0\right)}{\Delta E_z\left(0\right)}\right)\right]\Delta E_z, \end{align} \tag{ 60 }$$

using the relation between real-time t and time s. The amplitude A is given by the conditional phase condition $\int_0^{t_g} J(t^{\prime}) dt^{\prime} = 1/2$. For $J(t)\ll \Delta E_z$ and constant Zeeman splitting $\Delta E_z(t)\approx \Delta E_z$ we find

$$\begin{align} v_B\left(t\right)& = \frac{1}{2 \alpha}\text{log}\left(\frac{w\left(t\right)}{2 J_0 t_g}+1\right), \end{align} \tag{ 61 }$$

where we use $t = t_g s$ and the normalized window $\int_0^{t_g} w(t) dt = t_g$.

Another set of two-qubit gates can be accessed by driving the exchange interaction directly at the $\vert \widetilde{\uparrow\downarrow}\rangle \Longleftrightarrow\vert \widetilde{\downarrow\uparrow}\rangle$ resonance frequency $\nu_\text{ST} = \sqrt{\Delta E_z^2+J_0^2}$ [82, 83]. This gate exchanges the population (swapping) between the two states. For the resonant swap gate, dominant coherent errors are violations of the rotating wave approximation due to $1/t_g\sim \nu_\text{ST}\sim100$ MHz and the influence of higher harmonics due to the non-linear voltage-exchange relation.

In general, driving the barrier voltage yields $v_B(t) = v_{B,0}+v_{B,1}(t)\cos(2\pi \nu_\text{ST}t+\theta_j)$, which gives rise to

$$\begin{align} J\left(v_B\right) & = \mathcal{J}_0\left(v_B\right)+ \sum_{k> 0} 2\mathcal{J}_{k}\left(v_B\right) \cos\left(2\pi \nu_\text{ST} k t+ k \theta\right) , \end{align} \tag{ 62 }$$

$$\begin{align} \Delta E_z\left(v_B\right) & = \Delta \mathcal{E}_{z,0}\left(v_B\right) +\sum_{k> 0} 2\Delta \mathcal{E}_{z,k}\left(v_B\right) \cos\left(2\pi \nu_\text{ST} k t+ k \theta\right) , \end{align} \tag{ 63 }$$

where we have expressed J and $\Delta E_z$ in terms of Fourier coefficients with respect to the drive frequency $\nu_\text{ST}$. Due to the impact of the barrier voltage on the resonance frequency [32], we consider in our model $\Delta E_z \rightarrow \Delta E_z(v_{B}(t))$ with $\Delta E_z(v_{B,0}) = \Delta E_z$.

Without loss of generality, the dynamics is again projected on the odd-parity subspace spanned by $\lbrace\vert \uparrow\downarrow\rangle ,\vert \downarrow\uparrow\rangle \rbrace$ and described by Hamiltonian (53). To simplify the Hamiltonian, we perform a double basis transformation $U = \mathrm{e}^{-i (\pi\nu_\text{ST}t+\theta_j/2) \widetilde{\sigma}_z}\mathrm{e}^{-\frac{i}{2} \tan^{-1}\left(\frac{J_0}{\Delta E_z}\right)\sigma_y}$ before we apply our framework. The first transformation diagonalizes Hamiltonian (53) at t = 0, and the second moves us into the rotating frame with respect to the driving frequency and driving phase. The transformed and rotated Hamiltonian reads

$$\begin{align} H\left(t\right) & = -\frac{h }{2}J\left(t\right) +\frac{h}{2} \left( \frac{J_0\, J\left(t\right)+\Delta E_z\, \Delta E_z\left(t\right)}{\nu_\text{ST}}-\nu_\text{ST}-\frac{\dot{\theta}_j}{2\pi} \right)\widetilde{\sigma}_z \nonumber\\ & \quad + \frac{h}{2}\left(\frac{\Delta E_z J\left(t\right)-J_0 \Delta E_z\left(t\right) }{\nu_\text{ST}}\mathrm{e}^{2\pi i \nu_\text{ST} t +i \theta_j}\right) \widetilde{\sigma}_+ + h.c. \end{align} \tag{ 64 }$$

Within this frame, the target operation is then given by [82]

$$\begin{align} U_\text{ideal} = \mathrm{e}^{-i \pi \int_0^{t}dt^{\prime} \cos\left(2\pi \nu_\text{ST} t + \theta_j\right) \frac{\Delta E_z J\left(t\right)-J_0 \Delta E_z\left(t\right)}{\nu_\text{j}\left(t\right)}\widetilde{\sigma}_x } \end{align} \tag{ 65 }$$

which describes swap oscillations between the basis states.

We find (see appendix C.3) a single pair of error rates described by the operators $O = \widetilde{\sigma}_z\pm i \widetilde{\sigma}_y$. The erroneous transition matrix elements and the accumulated energy gap are then

$$\begin{align} g\left(t\right) & = \frac{J_0\,J\left(t\right)+\Delta E_z\, \Delta E_z\left(t\right)}{\nu_\text{ST}}-\nu_\text{ST}-\frac{\dot{\theta}_j}{2\pi} \mp i \sin\left(2\pi \nu_\text{ST} t + \theta_j\right) \frac{\Delta E_z \,J\left(t\right)-J_0 \Delta E_z\left(t\right)}{2\nu_\text{ST}}, \end{align} \tag{ 66 }$$

$$\begin{align} f\left(t\right) & = \pi\cos\left(2\pi \nu_\text{ST} t + \theta_j\right) \frac{\Delta E_z\, J\left(t\right)-J_0 \Delta E_z\left(t\right)}{\nu_\text{ST}}. \end{align} \tag{ 67 }$$

3.3.1. General case

In the general case, no closed-form analytical expressions can be derived for equation (62) and one has to rely on numerical techniques. Unfortunately, we cannot use the substitution (24) for the general Hamiltonian since f = 0 and proceed directly to the mitigation methods and omit the substitution. However, we note that all steps in our framework can be applied if we apply a (generalized) RWA [79] on Hamiltonian (64) to remove the oscillating components.

Our framework tells us for non-trivial complex matrix transition elements $g(t) = g_R(t)+i g_I(t)$ to use dynamical pulse shaping. We separate them into real g_R and imaginary g_I components and then perform integration by parts

$$\begin{align} &|\text{tr}\left(\mathcal{E}\,O-O\right)|^2 = \left|\int_0^{t_g}\left[g_R\left(t\right)+i g_I\left(t\right)\right]\mathrm{e}^{i \alpha\left(t\right)}dt \right|^2 \end{align} \tag{ 68 }$$

$$\begin{align} &\qquad\qquad\qquad\qquad\qquad\quad = \left|\int_0^{t_g}\left[g_R\left(t\right)+ f\left(t\right) \int_0^{t}g_I\left(t^{\prime}\right)dt^{\prime}\right] \mathrm{e}^{i \alpha\left(t\right)}dt \right|^2, \end{align} \tag{ 69 }$$

where $\alpha(t)$ is given by the antiderivative of f(t) (see equation (24)). Here, we restrict ourselves to solutions where the limits of integration vanish at $t = 0,t_g$. The optimized $\theta_j(t)$ is the solution of the integro-differential equation for the phase θ_j by plugging in equations (66) and (67).

3.3.2. Exponential exchange

We now consider the case of an exponential interaction $J(v_B)\propto \mathrm{e}^{2\alpha v_B}$ and linear frequency shifts, $\Delta E_z(v_B) = \Delta E_z + \beta\, (v_{B}-v_{B,0})$. This regime is experimentally accessed in [82, 83]. In this limit, exact analytical expressions can be derived. The Fourier coefficients are given by

$$\begin{align} \mathcal{J}_{k}\left(v_{B}\right) = J_0\mathcal{I}_k\left(2\alpha \left(v_{B}-v_{B,0}\right)\right) \end{align} \tag{ 70 }$$

for $k = 0,\pm 1$, where $\mathcal{I}_k$ denotes the modified Bessel function of order k. In this special case, the optimized dynamic pulse shape can be explicitly expressed

$$\begin{align} v_{B,1}\left(t\right)& = \frac{1}{2 \alpha}\text{log}\left(\frac{w\left(t\right)}{4 J_0 t_g}+1\right) \cos\left(2\pi \nu_\text{ST}+\theta_j\left(t\right)\right) \end{align} \tag{ 71 }$$

with the dynamic phase

$$\begin{align} \theta_{j}\left(t\right) & = 2\pi\int_0^t dt^{\prime} \frac{J_0\,\mathcal{J}_0\left(v_{B}\left(t^{\prime}\right)\right)+\Delta E_z^2-\nu^2_\text{ST}}{\nu_\text{ST}} + 2\pi\int_0^t dt^{\prime} \frac{\Delta E_z \mathcal{J}_1\left(t^{\prime}\right)- J_0 \beta\,v_{B,1}\left(t^{\prime}\right) }{\nu_\text{ST}^2}\nonumber\\[4pt] & \quad \times \bigg[\frac{\Delta E_z \mathcal{J}_0\left(t^{\prime}\right)}{2\pi\nu_{\text{ST}}^3}+\frac{\Delta E_z \mathcal{J}_1\left(t^{\prime}\right)}{4\pi\nu_{\text{ST}}^3} -\frac{J_0 \beta\,v_{B,1}\left(t^{\prime}\right)}{8\pi\nu_{\text{ST}}^3}-\sum_{k = 2}^\infty\frac{\Delta E_z \mathcal{J}_k\left(t^{\prime}\right)}{\left(k^2-1\right)\pi\nu_{\text{ST}}^3}\bigg]. \end{align} \tag{ 72 }$$

The first term originates from the conventional RWA and compensates the shift in resonance frequency $h\sqrt{\Delta E_z^2+J^2}$ due to non-linear exchange. The remaining term describes the driving-induced shift of the rotation angle and can be derived using a generalized RWA [79] (see appendix B).

Each large-scale qubit device needs consistently small single-qubit error rates to fulfill the requirements for error correction [86, 87]. A frequently reported number are infidelities of $1-F\lt0.5-1\%$ for the physical gate operations [86, 88], although based on many assumptions such as good initialization and readout and uncorrelated errors. Our simulations show that sufficiently small infidelities are within reach on state-of-the-art qubit devices using our framework. The simplest technique, synchronization of the Rabi-frequencies [46, 47], already shows small infidelities in our simulations. However, low-pass filters built into the electronic circuits to reduce high-frequency noise, are a detrimental error source for the synchronization technique. This is seen in figure 1(b) that shows the infidelity with (blue) and without (orange) a 150 MHz-Butterworth filter. The filter significantly reduces the dips as well as shifts the minima. On the other hand, our simulations (see figures 1(c) and (d)) show that high-fidelity operations can still be reached using static or dynamic pulse shaping to reach infidelities as small as 10⁻⁴ for gate times in the order of $t_g = 25$ ns with a frequency separation of $|\Delta E_z| = 100$ MHz. Note, that this, on the one hand, directly implies that gate times $t_g = 250$ ns are required if the frequency separation is reduced to $|\Delta E_z| = 10~\mathrm{MHz}$ as proposed in some architectures, e.g. [13]. On the other hand, a larger frequency separation allows for faster high-fidelity gate operations [35, 89] due to the infidelity being invariant for $t_g\times |\Delta E_z|$. In the former case, dynamic pulse shaping, such as the derivative removal by adiabatic gate (DRAG) protocol [42, 66], allows for an additional improvement.

Many experiments on single-spin qubits use the universal cz gate [28, 29] as their native high-fidelity two-qubit gate [19, 24, 32, 33, 35, 90] due to its simplicity and potential for scaling to larger arrays [26, 73]. The cz-class (cphase) gate

$$\begin{align} U_\text{cz} & = \mathrm{e}^{-i \Phi_{cz} \left(S_{1,z}S_{2,z}-\frac{1}{4}\right) } \end{align} \tag{ 74 }$$

$$\begin{align} \equiv \text{diag}\left(1,1,1,\mathrm{e}^{i \Phi_{cz}}\right) \end{align} \tag{ 75 }$$

with $\Phi_{cz} = 2\pi\int_0^{t_g}J(t)dt$ can be directly acquired using only single-qubit phase gates ¹ from the adiabatic phase evolution under the exchange interaction and can be transformed into a cnot gate using two single-qubit $R_{y}(\pm\frac{\pi}{2})$ gates. For the cz gate, the exchange interaction J(t) is pulsed, picking up a conditional phase $\Phi_{cz} = (2n+1)\pi$ with integer n. Since the cz gate is intended to be adiabatic innately, it is also directly (linearly) susceptible to low-frequency noise acting on the resonance frequencies. However, we claim that a substantial error in previous realizations [19, 91] is due to violations of the adiabaticity condition, which is more severely impacted due to the non-linear voltage-exchange relation and gives rise to bit-flip errors. Here, the adiabaticity condition is with respect to the frequency difference $|\Delta E_z|$ between the two spin qubits. Since the coherent errors of a cz two-qubit and resonant single-qubit gates are related and $t_g\times |\Delta E_z|$ is an invariant in the simulations, we know that the cz gate can, for example, be further improved using the Kaiser window or increasing the separation in qubit frequency. Figure 2(a) shows the simulated infidelity of a cz gate as a function of gate time t_g and the proportionality factor $\alpha\in [0,1]$ from the Tukey window, see equation (33). The oscillating infidelity as a function of exchange corresponds to an interference pattern of the diabatic contributions identical to the one observed in figures 1(b) and (c). On the other hand, a small frequency separation severely limits the performance, since long gate times $t_g = 300~\mathrm{ns}$ are required if the frequency separation is reduced to $|\Delta E_z| = 10~\mathrm{MHz}$ [92, 93]. In situations with a small frequency separation, a non-adiabatic cz gate can be realized as shown theoretically [28] and demonstrated experimentally [91] using the synchronization condition. Our simulations confirm such high-fidelity gates in figure 2(a) for λ = 0, where the Tukey pulse corresponds to a rectangular pulse. However, the diabatic implementation requires a precise timing, is sensitive to pulse imperfections such as filter effects, and is prone to dephasing due to low-frequency noise.

Zoom In Zoom Out Reset image size

We now compare our results to recent experiments demonstrating high-fidelity cz-gate operations with infidelities $1-F = 4\times 10^{-3}$ [32] and $1-F = 2\times 10^{-3}$ [35]. Both experiments are performed in isotopically enriched silicon quantum dots, which are prone to (low-frequency) charge noise. Low-frequency noise couples to the cz-gate through the diagonal matrix elements $S^z_{1,(2)}$ and $S^z_{1}S^z_{2}$ via the qubit frequencies $\Delta E_z$ and via the exchange interaction J. In figure 2 we highlighted the pulse shape, expected noise level, and pulse duration (star) of the extracted parameters in [32]. Our simulation predicts coherent errors as low as $1-F = 10^{-5}$, and total errors of $1-F = 2\times 10^{-4}$ are achievable. We speculate that the discrepancy of simulations and experiment is related to heating effects [75]. While we have not performed simulations using the parameters in [35] we can nevertheless predict the coherent infidelity using the invariant $t_g\times |\Delta E_z|$. The conversion results in a rectangular pulse with an effective gate duration $\widetilde{t}_g = 158~\mathrm{ns}$ due to the frequency difference $|\Delta E_z| = 396~\mathrm{MHz}$ (see appendix G for details). Our simulations predict a coherent infidelity of $1-F = 5\times 10^{-4}$ without and $1-F = 8\times 10^{-4}$ with incoherent noise sources.

Figure 2(b) shows the cz gate infidelity in general as a function of charge noise amplitude A and gate time t_g , showing the importance of the interplay between coherent and incoherent errors. Too fast gates suffer from coherent errors while too slow gates are prone to incoherent errors with an optimum depending on the charge noise amplitude. Considering realistic values from [32] (see dashed line in figure 2(b)) one can clearly observe that coherent errors from non-adiabaticity are dominating for $t_g\lt40~\mathrm{ns}$ and infidelities as low as $1-F = 10^{-4}$ are possible.

While arbitrary single-qubit gates combined with cz form a universal gate set for quantum circuits, it is often more efficient to include additional gates into the gate set. A frequently necessary gate in qubit architectures with nearest-neighbor couplings only is the swap gate, as it enables long-range qubit–qubit communication as well as read-out [13, 73, 82, 83]. The swap-class gate

$$\begin{align} U_\text{swap}\left(\phi\right) = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & \mathrm{e}^{i\Phi_\text{swap}} & 0 \\ 0 & \mathrm{e}^{i\Phi_\text{swap}} & 0 & 0 \\ 0 & 0 & 0 & 1 \\ \end{pmatrix}. \end{align} \tag{ 76 }$$

can be either accessed directly through a diabatic exchange pulse for small frequency differences $|\Delta E_z|\ll J$ [29, 91] or in general by driving the exchange interaction at the frequency difference $|\Delta E_z|$ between the corresponding qubits. Here we defined the swap-class gate via the additional phase $\Phi_\text{swap} = \pi \int_0^{t_g} J(t)dt$. For $\Phi_\text{swap} = n\pi$, we recover the classical swap-gate while the iswap, $\Phi_\text{swap} = (2n+1)\pi/2$, maximally entangles the qubits during the swap. Figure 3(a) shows the simulated infidelity of a swap-class gate for different static and dynamic pulse shapes, assuming perfect phase compensation. Due to the non-linear exchange interaction and short target gate times $|\Delta E_z|\sim 1/t_g$, dynamic pulse shaping greatly enhances the performance of the gates. For situations when the explicit phase $\Phi_\text{swap}$ matters, e.g. compiling Clifford gates using iswap gate, we provide two methods. The simplest method to obtain, for example, an iswap gate is to append a cz gate such that $\Phi_\text{swap}+\Phi_\text{cz}/2 = (2n+1)\pi/2$ [82]. Remarkably, one can also perform the compensation cz gate simultaneously with the swap gate by combining an ac and dc control signal. This comes with the advantage of a faster gate time at the cost of additional calibrations.

Zoom In Zoom Out Reset image size

The swap-class is directly (linear) susceptible to low-frequency noise coupling in via the exchange interaction J. Additionally, both discussed implementations have in common that they require careful calibration to compensate for the adiabatic phase acquisition from exchange, making them (at least) equally susceptible to low-frequency charge noise as the conventional cz gate. Figure 4 compares the infidelity of the different two-qubit gate implementations discussed in this paper. It is clearly visible that the cz gate always outperforms the swap-class gate. The lower fidelity of the swap gate is due to the overall larger conditional phase picked up, $\Phi_\text{swap}\approx 2\Phi_\text{cz}$.

Zoom In Zoom Out Reset image size

Plugging in $\nu_D = E_z+\frac{\Delta E_z}{2} - \Delta f(t)$ into equation (39) we get

$$\begin{align} H_\text{total}& = \frac{h}{2} \begin{pmatrix} -\Delta E_z+2\Delta f\left(t\right)-\frac{\dot{\theta}}{\pi} & \left(B^x_{2,1}-i B^y_{2,1}\right)/2 & \left(B^x_{1,1}-i B^y_{1,1}\right)/2 & 0 \\ \left(B^x_{2,-1}+i B^y_{2,-1}\right)/2 & \Delta E_z & 0 & \left(B^x_{1,1}-i B^y_{1,1}\right)/2 \\ \left(B^x_{1,-1}+i B^y_{1,-1}\right)/2 & 0 & -\Delta E_z & \left(B^x_{2,1}-i B^y_{2,1}\right)/2 \\ 0 & \left(B^x_{1,-1}+i B^y_{1,-1}\right)/2 & \left(B^x_{2,-1}+i B^y_{2,-1}\right)/2 & \Delta E_z-2\Delta f\left(t\right)+\frac{\dot{\theta}}{\pi} \\ \end{pmatrix}, \end{align} \tag{ C.1 }$$

$$\begin{align} H_\text{ideal}& = \frac{h}{2} \begin{pmatrix} -\Delta E_z+\Delta f\left(t\right)-\frac{\dot{\theta}}{2\pi} & 0 & B^x_{1,1}/2 & 0 \\ 0 & \Delta E_z-\Delta f\left(t\right)+\frac{\dot{\theta}}{2\pi} & 0 & B^x_{1,1}/2 \\ B^x_{1,-1}/2 & 0 & -\Delta E_z+\Delta f\left(t\right)-\frac{\dot{\theta}}{2\pi} & 0 \\ 0 & B^x_{1,-1}/2 & 0 & \Delta E_z-\Delta f\left(t\right)+\frac{\dot{\theta}}{2\pi} \\ \end{pmatrix}, \end{align} \tag{ C.2 }$$

$$\begin{align} H_\text{error}& = U^\dagger_\text{ideal} \frac{h}{2} \begin{pmatrix} \Delta f\left(t\right)-\frac{\dot{\theta}}{2\pi} & \left(B^x_{2,1}-i B^y_{2,1}\right)/2 & -i B^y_{1,1}/2 & 0 \\ \left(B^x_{2,-1}+i B^y_{2,-1}\right)/2 & \Delta f\left(t\right)-\frac{\dot{\theta}}{2\pi} & 0 & -i B^y_{1,1}/2 \\ i B^y_{1,-1}/2 & 0 & -\Delta f\left(t\right)+\frac{\dot{\theta}}{2\pi} & \left(B^x_{2,1}-i B^y_{2,1}\right)/2 \\ 0 & i B^y_{1,-1}/2 & \left(B^x_{2,-1}+i B^y_{2,-1}\right)/2 & -\Delta f\left(t\right)+\frac{\dot{\theta}}{2\pi} \\ \end{pmatrix}U_\text{ideal}. \end{align} \tag{ C.3 }$$

The corresponding error channels follow from $O_k = \vert k_1\rangle \langle k_2\vert$ with $k_1\neq k_2$, where $\vert k_{1,2}\rangle$ are eigenstates of equation (C.2). As Hamiltonian (C.1) describes an uncoupled two-qubit system, the dynamics of each qubit can be treated independently, reducing the number of non-zero channels to two per qubit. The dynamics of qubit 1 and qubit 2 are described by a x-rotation and z-rotation that leads to the following error channels

$$\begin{align} O_1 = S_1^{\,z}\mp i S_1^{\,y}, \end{align} \tag{ C.4 }$$

$$\begin{align} O_2 = S_2^x\pm i S_2^y. \end{align} \tag{ C.5 }$$

Plugging in $B_{\perp,1} = B_{\perp,2} = \nu_D = \dot{\theta} = 0$ into equation (38) and diagonalizing we get the Hamiltonians

$$\begin{align} H_\text{tot}& = \frac{h}{2} \begin{pmatrix} 2E_{z} & 0 & 0 & 0\\ 0 & \sqrt{\Delta E_z^2+J^2} - J & -\frac{\Delta E_z\, \dot{J}-\dot{\Delta E_z} J}{2\pi\nu_\text{ST}^2} & 0 \\ 0 & -\frac{\Delta E_z\, \dot{J}-\dot{\Delta E_z} J}{2\pi\nu_\text{ST}^2} & -\sqrt{\Delta E_z^2+J^2} -J & 0 \\ 0 & 0 & 0 & -2E_{z} \\ \end{pmatrix}, \end{align} \tag{ C.6 }$$

$$\begin{align} H_\text{ideal}& = \frac{h}{2} \begin{pmatrix} 2E_{z} & 0 & 0 & 0\\ 0 & \sqrt{\Delta E_z^2+J^2} - J & 0 & 0 \\ 0 & 0 & -\sqrt{\Delta E_z^2+J^2} -J & 0 \\ 0 & 0 & 0 & -2E_{z} \\ \end{pmatrix}, \end{align} \tag{ C.7 }$$

$$\begin{align} H_\text{error}& = \frac{h}{2} \begin{pmatrix} 0 & 0 & 0 & 0\\ 0 & 0 & -\frac{\Delta E_z\, \dot{J}-\dot{\Delta E_z} J}{2\pi\nu_\text{j}^2}\mathrm{e}^{-i\pi \int_0^t \sqrt{\Delta E_z^2+J^2} dt} & 0 \\ 0 & -\frac{\Delta E_z\, \dot{J}-\dot{\Delta E_z} J}{2\pi\nu_\text{j}^2}\mathrm{e}^{i\pi \int_0^t \sqrt{\Delta E_z^2+J^2} dt} & 0 & 0 \\ 0 & 0 & 0 & 0 \\ \end{pmatrix}. \end{align} \tag{ C.8 }$$

The corresponding error channels follows from $O_k = \vert k_1\rangle \langle k_2\vert$ with $k_1\neq k_2$, where $\vert k_{1,2}\rangle$ are eigenstates of equation (C.7). Hamiltonian (C.6) can be effectively reduced to a single-qubit system in the space spanned by the odd parity states $\lbrace\vert \uparrow\downarrow\rangle ,\vert \downarrow\uparrow\rangle \rbrace$ (see appendix A). In this subspace, equation (C.7) describes a z-rotation, giving rise to

$$\begin{align} O = \left(\widetilde{\sigma}_x\pm i \widetilde{\sigma}_y\right)/2 = \left\vert \widetilde{\uparrow\downarrow}\right\rangle \left\langle \widetilde{\downarrow\uparrow}\right\vert , \end{align} \tag{ C.9 }$$

where the tilde labels the instantaneous eigenstates.

For the resonant swap gate the Hamiltonians are

$$\begin{align} H_\text{tot}& = \frac{h}{2} \begin{pmatrix} 2E_{z} & 0 & 0 & 0\\ 0 & \frac{J_0J\left(t\right)+\Delta E_z \Delta E_z\left(t\right)-\nu_\text{ST}^2}{\nu_\text{ST}} -\frac{\dot{\theta_j}}{2\pi} - J & \frac{\Delta E_z\, J\left(t\right)-J_0 \Delta E_z}{\nu_\text{ST}}\mathrm{e}^{2\pi i \nu_\text{ST} t +i \theta_j} & 0 \\ 0 & \frac{\Delta E_z\, J\left(t\right)-J_0 \Delta E_z}{\nu_\text{ST}}\mathrm{e}^{-2\pi i \nu_\text{ST} t -i \theta_j} & -\frac{J_0J\left(t\right)+\Delta E_z \Delta E_z\left(t\right)-\nu_\text{ST}^2}{\nu_\text{ST}} +\frac{\dot{\theta_j}}{2\pi} - J & 0 \\ 0 & 0 & 0 & -2E_{z} \\ \end{pmatrix}, \end{align} \tag{ C.10 }$$

$$\begin{align} H_\text{ideal}& = \frac{h}{2} \begin{pmatrix} 2E_{z} & 0 & 0 & 0\\ 0 & 0 & \frac{\Delta E_z\, J\left(t\right)-J_0 \Delta E_z}{\nu_\text{ST}}\cos\left(2\pi \nu_\text{ST} t +\theta_j\right) & 0 \\ 0 & \frac{\Delta E_z\, J\left(t\right)-J_0 \Delta E_z}{\nu_\text{ST}}\cos\left(2\pi \nu_\text{ST} t +\theta_j\right) & 0 & 0 \\ 0 & 0 & 0 & -2E_{z} \\ \end{pmatrix}, \end{align} \tag{ C.11 }$$

$$\begin{align} H_\text{error}& = U^\dagger_\text{ideal}\frac{h}{2} \begin{pmatrix} 0 & 0 & 0 & 0\\ 0 & \frac{J_0J\left(t\right)+\Delta E_z \Delta E_z\left(t\right)-\nu_\text{ST}^2}{\nu_\text{ST}} -\frac{\dot{\theta_j}}{2\pi} - J & i\frac{\Delta E_z\, J\left(t\right)-J_0 \Delta E_z}{\nu_\text{ST}}\sin\left(2\pi \nu_\text{ST} t +\theta_j\right) & 0 \\ 0 & -i\frac{\Delta E_z\, J\left(t\right)-J_0 \Delta E_z}{\nu_\text{ST}}\sin\left(2\pi \nu_\text{ST} t +\theta_j\right) & -\frac{J_0J\left(t\right)+\Delta E_z \Delta E_z\left(t\right)-\nu_\text{ST}^2}{\nu_\text{ST}} +\frac{\dot{\theta_j}}{2\pi} - J & 0 \\ 0 & 0 & 0 & 0 \\ \end{pmatrix}U_\text{ideal}. \end{align} \tag{ C.12 }$$

The corresponding error channels follows from $O_k = \vert k_1\rangle \langle k_2\vert$ with $k_1\neq k_2$, where $\vert k_{1,2}\rangle$ are eigenstates of equation (C.11). Hamiltonian (C.10) can be effectively reduced to a single-qubit system in the space spanned by the odd parity states $\lbrace\vert \uparrow\downarrow\rangle ,\vert \downarrow\uparrow\rangle \rbrace$ (see appendix A). In this subspace, equation (C.11) describes a x-rotation, giving rise to

$$\begin{align} O = \widetilde{\sigma}_z\pm i \widetilde{\sigma}_y, \end{align}$$

where the tilde labels the instantaneous eigenstates.