Main

Quantum computation involves the execution of a large number of elementary operations that take a qubit register through the steps of a quantum algorithm6. A major challenge is to implement these operations with sufficient accuracy to arrive at a reliable outcome, even in the presence of decoherence and other error sources. The higher the accuracy, or fidelity, of the operations, the higher the likelihood that near-term applications for quantum computers come within reach7. Furthermore, for most presently known algorithms, the number of operations that must be concatenated will unavoidably lead to excessive accumulation of errors, and these errors must be removed using quantum error correction1. Correcting quantum errors faster than they occur is possible when the error probability per operation is below a certain threshold, known as the fault-tolerance threshold. For the widely considered surface code, for instance, the fault-tolerance threshold is between 0.6% and 1%, under certain assumptions, albeit at the cost of a large redundancy in the number of physical qubits2,3.

Among all the candidate platforms, electron spins in semiconductor quantum dots have advantages, such as their long coherence times8, small footprint9, the potential for scaling up10 and the compatibility with advanced semiconductor manufacturing technology4. Single-qubit operations of spin qubits in quantum dots achieve fidelities of 99.9% (refs. 11,12) but the two-qubit gate fidelities reported vary from 92% to 98% (refs. 13,14). This has limited the two-qubit Bell-state fidelities to 94% (ref. 15) and quantum algorithms implemented with spin qubits gave only coarsely accurate outcomes16,17. Pushing the two-qubit gate fidelity well beyond 99% requires not only low charge-noise levels and the elimination of nuclear spins by isotopic enrichment but also careful Hamiltonian engineering.

In this paper, using a precisely engineered two-qubit interaction Hamiltonian, we report the demonstration of single-qubit and two-qubit gates with fidelities above 99.5%. We use gate-set tomography (GST) not only to characterize the gates and to quantify the fidelity but also to improve the gate calibration. The high-fidelity gates allow us to compute the dissociation energy of molecular hydrogen with a variational quantum eigensolver (VQE) algorithm, reaching an accuracy for the dissociation energy of around 20 mHa, limited by readout errors.

We use a gate-defined double quantum dot in an isotopically enriched 28Si/SiGe heterostructure17 (Fig. 1a), with each dot occupied by a single electron (see Methods). The spin states of the electrons serve as qubits. The spin states are measured with the help of a sensing quantum dot (SQD), which is capacitively coupled to the qubit dots18. A micromagnet on top of the device provides a magnetic field gradient enabling electric-dipole spin resonance19 and separates the resonance frequencies of the qubits in the presence of an external magnetic field (~320 mT) to 11.993 GHz (Q1) and 11.890 GHz (Q2). Single-qubit X and Y gates are implemented by frequency-multiplexed microwave signals applied to gate MW and virtual Z gates are implemented by a phase update of the reference frame20. The plunger gates (LP and RP) control the chemical potentials of the quantum dots.

Fig. 1: Two-qubit device and symmetry operating point.
figure 1

a, Scanning electron microscopy images of a device similar to that used here, showing the quantum dot gate pattern and the micromagnet on top (the device used in the experiment has an additional screening gate above the fine gates17). The scale bar in the left panel denotes 500 nm. The scale bar in the right panel denotes 100 nm. b, Control paths for determining the symmetry operation point in the charge-stability diagram. (M, N) represent the number of electrons in the dots underneath the tip of LP and RP, respectively. a.u., arbitrary units. c, Pulse sequence schematic of a decoupled controlled-phase operation interleaved in a Ramsey interference sequence on Q1. d, Spin-up probability of Q1 after the Ramsey sequence in c, as a function of the detuning in the double-dot potential and the total duration of the barrier voltage pulses.

The native two-qubit gate for spin qubits uses the exchange interaction21,22, originating from the wave-function overlap of electrons in neighbouring dots. This selectively shifts the energy of the antiparallel spin states and, thus, enables an electrically pulsed adiabatic conditional Z (CZ) gate8,16,23. The barrier gate (B) controls the tunnel coupling between the dots, allowing the precise tuning of the exchange coupling from <100 kHz to 20 MHz. To minimize the sensitivity to charge noise, we activate the exchange coupling while avoiding a tilt in the double-dot potential24,25 (Fig. 1a). This symmetric condition can be determined accurately by decoupled adiabatic exchange pulses inside a Ramsey sequence (Fig. 1c, d). The tunnel barrier is controlled by simultaneously pulsing gate B and compensating LP and RP to avoid shifts in the electrochemical potentials24, constituting a virtual barrier gate. The detuning between quantum dots is controlled by additional offsets to the LP and RP pulses in opposite directions. As the decoupling pulses remove additional single-qubit phase accumulation from electron movement in the magnetic field gradient, the spin-up probability of Q1 results in a symmetric chevron pattern, with the symmetry point at the centre (Fig. 1d).

Among the various quantum benchmarking techniques, quantum process tomography (QPT) is designed to reconstruct all details in a target process6. Owing to the susceptibility of QPT to state preparation and measurement (SPAM) errors, self-consistent benchmarking techniques such as GST26 and alternative techniques such as randomized benchmarking27 have been developed. In contrast to randomized benchmarking, GST inherits the advantage of QPT in that it reports the detailed process, which allows us to isolate Hamiltonian errors from stochastic errors and to correct for such errors in the control signals (Extended Data Fig. 5). In addition, GST accounts for gate-dependent errors. We benchmark the fidelities of a universal gate set using GST26,28 (Fig. 2a). The gate set we choose contains an idle gate (I), sequentially operated single-qubit π/2 rotations about the \(\hat{x}\) and \(\hat{y}\) axes for each qubit (\({{\rm{X}}}_{{{\rm{Q}}}_{1}}\), \({{\rm{Y}}}_{{{\rm{Q}}}_{1}}\), \({{\rm{X}}}_{{{\rm{Q}}}_{2}}\) and \({{\rm{Y}}}_{{{\rm{Q}}}_{2}}\)) and a two-qubit controlled-phase (CZ) gate. A total of 36 fiducial sequences containing \(\{{\rm{n}}{\rm{u}}{\rm{l}}{\rm{l}},{{{\rm{X}}}_{{{\rm{Q}}}_{i}}}^{n=1,2,3},{{{\rm{Y}}}_{{{\rm{Q}}}_{j}}}^{n=1,3}\}\) on each qubit, where null (unlike the idle gate) has no waiting time, are used to tomographically measure the two-qubit state. These fiducials are interleaved by germ sequences and their powers up to a sequence depth of 16. Germs are short sequences of gates taken from the universal gate set (see Methods). They are repetitively executed to amplify different types of gate errors in the gate set, such that SPAM errors can be isolated. GST allows using a maximum-likelihood estimator to compute completely positive and trace-preserving process matrices for each element of the gate set6. The gate fidelity can be calculated by comparing the measured process using the Pauli transfer matrix (PTM), \({ {\mathcal M} }_{\exp }\), with the ideal PTM, \({ {\mathcal M} }_{{\rm{ideal}}}\), \({F}_{{\rm{gate}}}=({\rm{Tr}}({ {\mathcal M} }_{\exp }^{-1}{ {\mathcal M} }_{{\rm{ideal}}})+d)/\)[d(d+1)], where d is the dimension of the Hilbert space. These process matrices provide a detailed error diagnosis of the gate set, allowing for efficient feedback calibration29 (Fig. 2a). Analysing the error generator \( {\mathcal L} =\,\log ({ {\mathcal M} }_{\exp }{ {\mathcal M} }_{{\rm{ideal}}}^{-1})\) provides easy access to information. For example, coherent Hamiltonian errors can be isolated from incoherent stochastic errors and single-qubit errors can be isolated from each other and from two-qubit errors30.

Fig. 2: Gate-set tomography and single-qubit gate.
figure 2

a, Workflow of the GST experiment. Coloured blocks show the input and output fiducial sequences (Fidi and Fido, orange) and the germ sequences (green). A few examples of single-qubit germ sequences are listed. The outcome is used to adjust pulse parameters in the next run. b, c, PTMs of \({{\rm{X}}}_{{{\rm{Q}}}_{1}}\) and \({{\rm{Y}}}_{{{\rm{Q}}}_{1}}\) in the subspace of Q1. The red (blue) bars are theoretically +1 (−1) and are measured to be positive (negative). The brown (green) bars are theoretically 0 (0) but measured to be positive (negative). Pin and Pout are the input and output operators, respectively. d, Experimentally measured PTM of \({{\rm{Y}}}_{{{\rm{Q}}}_{1}}\otimes {{\rm{I}}}_{{{\rm{Q}}}_{2}}\) in the complete two-qubit space. The colour code is the same as in b, c.

Figure 2b, c shows the reduced PTMs of \({{\rm{X}}}_{{{\rm{Q}}}_{1}}\) and \({{\rm{Y}}}_{{{\rm{Q}}}_{1}}\) operations in the Q1 subspace and Fig. 2d shows the full PTM of \({{\rm{Y}}}_{{{\rm{Q}}}_{1}}\) in two-qubit space (\({{\rm{Y}}}_{{{\rm{Q}}}_{1}}\otimes {{\rm{I}}}_{{{\rm{Q}}}_{2}}\)) containing additional errors from decoherence and crosstalk on Q2 while operating Q1 (see Extended Data Figs. 1 and 2 for other PTMs) and from unintentional entanglement due to a residual exchange interaction. The average single-qubit gate fidelity is 99.72% in the single-qubit subspace (\({{\rm{X}}}_{{{\rm{Q}}}_{1}}\): 99.68%; \({{\rm{Y}}}_{{{\rm{Q}}}_{1}}\): 99.73%; \({{\rm{X}}}_{{{\rm{Q}}}_{2}}\): 99.61%; \({{\rm{Y}}}_{{{\rm{Q}}}_{2}}\): 99.87%; see Extended Data Fig. 2 for all error bars). A metric that is rarely reported is the single-qubit gate fidelity in the full two-qubit space, here 99.16% on average (see Methods and Extended Data Fig. 1). These results highlight that single-qubit benchmarking is not sufficient to identify all errors occurring during single-qubit operations. By analysing the error generators, we find that errors from uncorrelated dephasing of the idling qubit dominate the drop in single-qubit gate fidelity when characterized in the two-qubit space. Coherent, microwave-induced phase shifts—the main source of crosstalk errors—have been corrected by applying a compensating phase gate to the idling qubit (Extended Data Fig. 4). The elimination of idling errors and other crosstalk errors from the microwave drive, such as through heating effects, will be a crucial step to improve the quality of the single-qubit operations further.

For a high-fidelity adiabatic CZ gate, precise control of the exchange coupling, J, between the two qubits is required. Specifically, in order to avoid unintended state transitions due to non-adiabatic dynamics, we must be able to carefully shape the envelope of  J. We characterize J over a wide range using a Ramsey sequence interleaved by a virtual barrier pulse with incremental amplitude vB. Figure 3a shows the measured frequency shift of each qubit as functions of the barrier pulse amplitude and the state of the other qubit. The exchange interaction is modelled to be exponentially dependent on the barrier pulse amplitude \(J({v}_{{\rm{B}}})\propto {{\rm{e}}}^{2\alpha {v}_{{\rm{B}}}}\) (refs. 31,32). The micromagnet-induced single-qubit frequency shifts are approximated by linear functions within the voltage window of the CZ gate in the numerical simulations. By fitting the measured datasets simultaneously to theoretical models (see Methods), J can be extracted very precisely as the difference between the two conditional frequencies of each qubit16,33 (Fig. 3b). The barrier pulse \({v}_{{\rm{B}}}\propto \,\log ({A}_{{v}_{{\rm{B}}}}(1-\,\cos (2{\rm{\pi }}t/{t}_{{\rm{g}}{\rm{a}}{\rm{t}}{\rm{e}}}))/2)\) (Fig. 3d) compensates the exponential dependence such that  J (1 − cos(2πt/tgate)) follows a cosine window function, which ensures good adiabaticity34 (Fig. 3e). In addition, the virtual gates are calibrated such that the symmetric operation point is maintained for each barrier setting, minimizing the influence of charge noise via the double-dot detuning. The most relevant remaining noise sources include charge noise, affecting  J through fluctuations in the virtual barrier gate δvB, and fluctuating qubit frequencies \(\delta {f}_{{{\rm{Q}}}_{1}}\) and \(\delta {f}_{{{\rm{Q}}}_{2}}\) from charge noise entering through artificial spin–orbit coupling from the micromagnet and residual nuclear spin noise coupling through the hyperfine interaction. By analysing the decay of the Ramsey oscillations at each transition frequency, individual dephasing times \({T}_{2}^{\ast }\) can be extracted and, from there, also δvB, \(\delta {f}_{{{\rm{Q}}}_{1}}\) and \(\delta {f}_{{{\rm{Q}}}_{2}}\) (Fig. 3c).

Fig. 3: Hamiltonian engineering of exchange interaction.
figure 3

a, Frequency detuning of each qubit conditional on the state of the other qubit as a function of barrier pulse amplitude. The horizontal axis shows the real voltage applied to gate B. b, Exchange strength as a function of barrier pulse amplitude. The data are extracted directly from a. c, \({T}_{2}^{\ast }\) of each qubit conditional on the state of the other qubit as a function of barrier pulse amplitude (same colour code as in a). Each data point is averaged for about 8 min. By fitting the \({T}_{2}^{\ast }\) values to a quasistatic noise model (solid lines, see Methods), the low-frequency amplitudes of the fluctuations are estimated as \(\delta {f}_{{{\rm{Q}}}_{1}}=11\,{\rm{k}}{\rm{H}}{\rm{z}}\), \(\delta {f}_{{{\rm{Q}}}_{2}}=24\,{\rm{k}}{\rm{H}}{\rm{z}}\) and δvB = 0.4 mV. d, Shape of the barrier pulse, designed to achieve a high-fidelity CZ gate. e, The cosine-shaped J envelope seen by the qubits during the pulse shown in d.

Figure 4a shows an example GST pulse sequence that contains twice in a row the germ \([{\rm{C}}{\rm{Z}},{{\rm{X}}}_{{{\rm{Q}}}_{2}},{{\rm{Y}}}_{{{\rm{Q}}}_{1}},{\rm{C}}{\rm{Z}},{{\rm{Y}}}_{{{\rm{Q}}}_{2}},{{\rm{X}}}_{{{\rm{Q}}}_{1}}]\). The PTM of the CZ gate obtained from GST is shown in Fig. 4b. Using the detailed information from the error generator to fine-tune the calibration parameters, we can achieve a CZ fidelity of 99.65 ± 0.15% (Extended Data Figs. 4 and 5). Error bars included here and elsewhere are the 2σ ≈ 95% confidence intervals computed using the Hessian of the loglikelihood function35. The CZ error generator reveals that, at this point, incoherent errors dominate. The virtual barrier gate technique used here efficiently suppresses crosstalk errors during two-qubit gates. Therefore, we expect the CZ fidelity to be mostly affected by dephasing errors of idling qubits in a larger space, which can be corrected for using decoupling pulses. From the obtained PTMs, we can numerically estimate Bell-state fidelities by multiplications of the PTMs necessary to construct the corresponding state, giving an estimate of 97.75%–98.42%, neglecting SPAM errors, for the four Bell states (Fig. 4c and Extended Data Fig. 3).

Fig. 4: High-fidelity two-qubit gate.
figure 4

a, A sequence of pulses generated by the arbitrary waveform generators in an example GST sequence. The purple waveforms show the in-phase component of X/Y gates. The CZ gate is indicated by the orange pulse of gate B and the blue and red compensation pulses of gate LP and gate RP. b, Experimentally determined PTM of a CZ gate. The colour code is the same as in Fig. 2. c, Left, the quantum circuit used to reconstruct the Bell state \(|{\Psi }^{+}\rangle =(|01\rangle +|10\rangle )/\sqrt{2}\) based on the corresponding PTMs. Right, the real part of the reconstructed density matrix of the |Ψ+ state. The colour code is the same as in Fig. 2, except that red (blue) bars here are theoretically +0.5 (−0.5).

Next, we use the high-fidelity gate set in the context of an actual application, in order to provide a quantitative benchmark for future work under realistic conditions. Specifically, we implement a VQE algorithm to compute the ground-state energy of molecular hydrogen (H2) (Fig. 5a). In a VQE algorithm, a quantum processor is used to implement a classically inefficient subroutine (see Methods and Extended Data Fig. 6). The second quantized H2 Hamiltonian can be mapped onto two qubits under the Bravyi–Kitaev (BK) transformation H = h0II + h1ZI + h2IZ + h3ZZ + h4XX + h5YY. Here I, X, Y and Z are Pauli operators, for example, ZI is shorthand for Z  I, and the coefficients h0h5 are classically computable functions of the internuclear distance, R. Figure 5b shows the schematic of the VQE algorithm and its circuit implementation for a H2 molecule. The qubit is initialized in |01, which represents double occupation of the lowest molecular orbital, corresponding to the Hartree–Fock (HF) ground state. A parameterized ansatz state is then prepared by considering single and double excitation, which, after the BK transformation, yields \(|\psi (\theta )\rangle ={{\rm{e}}}^{-{\rm{i}}\theta {\rm{XY}}}|01\rangle \), with θ the parameter to variationally optimize. By performing partial tomography on the ansatz state with an initial guess θ0, the expectation value of the Hamiltonian for |ψ(θ0) can be calculated. A classical computer can efficiently compute the next guess θ1 as the new input for the quantum computer. This loop is iterated until the result converges. For a H2 molecule, there is only one parameter θ to optimize, thus, a scan of the entire parameter range of 2π with finite samples is sufficient to interpolate the smoothly changing measured expectation values. This emulates a real variational algorithm, where θ can be estimated to arbitrary precision by increasing the number of repetitions to suppress statistical fluctuations36. Figure 5c shows the partial tomography result after normalization of the visibility window. The data demonstrate high-quality phase control in the quantum circuits. The deviations in the odd-parity expectation values indicate correlations in the readout of the two qubits37. Figure 5d shows the energy curves of the H2 molecule from both theory38 and the VQE experiment. We observe a minimum energy at around 0.72 Å and an error of approximately 20 mHa at the theoretical bond length 0.7414 Å, mainly attributed to slow drift in the readout parameters. This accuracy matches the results obtained using superconducting and trapped ion qubits with comparable gate fidelities36,39.

Fig. 5: Variational quantum eigensolver.
figure 5

a, Lowest two molecular orbitals of a H2 molecule, formed by the 1s orbitals of two hydrogen atoms. b, The quantum circuit to implement the VQE algorithm for a H2 molecule. The orange block prepares the HF initial state by flipping Q2. The circuit in green blocks creates the parameterized ansatz state. \(-{{\rm{X}}}_{{{\rm{Q}}}_{i}}\) and \(-{{\rm{Y}}}_{{{\rm{Q}}}_{j}}\) include virtual Z gates. CNOT gates are compiled as \([-{{\rm{Y}}}_{{{\rm{Q}}}_{2}},{\rm{C}}{\rm{Z}},{{\rm{Y}}}_{{{\rm{Q}}}_{2}}]\). To make use of the high-fidelity CZ gate, such compilation is preferred instead of using a single controlled-phase gate with incremental length for creating the parameterized ansatz state. c, Expectation values of the operators in the two-qubit Hamiltonian under BK transformation as a function of θ. Black solid lines show the predicted values. The coloured solid lines are sinusoidal fits to the data (and a constant fit for the case of ZZ). d, Potential energy of the H2 molecule at varying R. The VQE data are normalized to the theoretical energy at large R to directly compare the dissociation energy with the theoretical value. The inset shows the error in the normalized experimental data.

The two-qubit gate with fidelity above 99.5% and single-qubit gate fidelities in the two-qubit gate space above 99% on average place semiconductor spin qubit logic at the error threshold of the surface code. Recently, a two-qubit operation between nuclear spin qubits in silicon, mediated by an electron spin qubit, has been demonstrated to surpass 99% fidelity as well, further highlighting that semiconductor spin qubits offer precise two-qubit logic40. Independent studies have shown spin qubit readout with a fidelity above 98% in only a few μs (ref. 41), with further improvements underway42. Combining high-fidelity initialization, readout and control into a demonstration of fault tolerance poses several key challenges to be overcome. First, sufficiently large and reliable quantum dot arrays must be constructed, with good connectivity between the qubits. Second, the fidelities achieved in small-scale systems must be maintained across such larger systems, which will require reducing idling and crosstalk errors. The same advances will allow us to implement more sophisticated algorithms in the noisy intermediate-scale quantum era, such as solving energies involving excited states of more complex molecules.

Methods

Measurement setup

The measurement setup and device are similar to those used in ref. 17. We summarize a few key points and all the differences here. The gates LP, RP and B are connected to arbitrary waveform generators (AWGs, Tektronix 5014C) via coaxial cables. The position in the charge-stability diagram of the quantum dots is controlled by voltage pulses applied to LP and RP. Linear combinations of the voltage pulses applied to B, LP and RP are used to control the exchange coupling between the two qubits at the symmetry point. The compensation coefficients are vLP/vB = −0.081 and vRP/vB = 0.104. A vector signal generator (VSG, Keysight E8267D) is connected to gate MW and sends frequency-multiplexed microwave bursts (not necessarily time-multiplexed) to implement electric-dipole spin resonance (EDSR). The VSG has two I/Q input channels, receiving I/Q modulation pulses from two channels of an AWG. I/Q modulation is used to control the frequency, phase and length of the microwave bursts. The current signal of the sensing quantum dot is converted to a voltage signal and recorded by a digitizer card (Spectrum M4i.44), and then converted into 0 or 1 by comparing it to a threshold value.

Two differences between the present setup and that in ref. 17 are that (1) the programmable mechanical switch is configured such that gate MW is always connected to the VSG and not to the cryo-CMOS control chip and (2) a second AWG of the same model is connected to gate B, with its clock synchronized to the first AWG.

Gate calibration

In the gate set used in this work, \(\{{\rm{I}},{{\rm{X}}}_{{{\rm{Q}}}_{1}},{{\rm{Y}}}_{{{\rm{Q}}}_{1}},{{\rm{X}}}_{{{\rm{Q}}}_{2}},{{\rm{Y}}}_{{{\rm{Q}}}_{2}},{\rm{C}}{\rm{Z}}\}\), the duration of the I gate and the CZ gate are set to 100 ns, and we calibrate and keep the amplitudes of the single-qubit drives fixed and in the linear-response regime, where the Rabi frequency is linearly dependent on the driving amplitude. The envelopes of the single-qubit gates are shaped following a ‘Tukey’ window, as it allows adiabatic single-qubit gates with relatively small amplitudes, thus, avoiding the distortion caused by a nonlinear response. The general Tukey window of length tp is given by

$$W(t,r)=\{\begin{array}{cc}\frac{1}{2}\left[1-\,\cos \,\left(\frac{2{\rm{\pi }}t}{r{t}_{{\rm{p}}}}\right)\right] & 0\le t\le \frac{r{t}_{{\rm{p}}}}{2}\\ 1 & \frac{r{t}_{{\rm{p}}}}{2} < t < {t}_{{\rm{p}}}-\frac{r{t}_{{\rm{p}}}}{2}\\ \frac{1}{2}\left[1-\,\cos \,\left(\frac{2{\rm{\pi }}({t}_{{\rm{p}}}-t)}{r{t}_{{\rm{p}}}}\right)\right] & {t}_{{\rm{p}}}-\frac{r{t}_{{\rm{p}}}}{2}\le t\le {t}_{{\rm{p}}}\end{array},$$
(1)

where r = 0.5 for our pulses. Apart from these fixed parameters, there are 11 free parameters that must be calibrated: single-qubit frequencies \({f}_{{{\rm{Q}}}_{1}}\) and \({f}_{{{\rm{Q}}}_{2}}\), burst lengths for single-qubit gates tXY1 and tXY2, phase shifts caused by single-qubit gates on the addressed qubit itself ϕ11 and ϕ22, phase shifts caused by single-qubit gates on the unaddressed ‘victim qubit’ ϕ12 and ϕ21 (ϕ12 is the phase shift on Q1 induced by a gate on Q2 and similar for ϕ21), the peak amplitude of the CZ gate \({A}_{{v}_{{\rm{B}}}}\) and phase shifts caused by the gate voltage pulses used for the CZ gate on the qubits θ1 and θ2 (in addition, we absorb into θ1 and θ2 the 90° phase shifts needed to transform diag(1, i, i, 1) into diag(1, 1, 1, −1)).

For single-qubit gates, \({f}_{{{\rm{Q}}}_{1}}\) and \({f}_{{{\rm{Q}}}_{2}}\) are calibrated by standard Ramsey sequences, which are automatically executed every 2 h, at the beginning and in the middle (after 100 times the average of each sequence) of the GST experiment. The EDSR burst times tXY1 and tXY2 are initially calibrated by an AllXY calibration protocol43. The phases ϕ11, ϕ12, ϕ21 and ϕ22 are initially calibrated by measuring the phase shift of the victim qubit (Q1 for ϕ11 and ϕ21; Q2 for ϕ22 and ϕ12) in a Ramsey sequence interleaved by a pair of \([{{\rm{X}}}_{{{\rm{Q}}}_{i}},-{{\rm{X}}}_{{{\rm{Q}}}_{i}}]\) gates on the addressed qubit (Q1 for ϕ11 and ϕ12; Q2 for ϕ22 and ϕ21) (Extended Data Fig. 4).

The optimal pulse design presented in Fig. 3 gives a rough guidance of the pulse amplitude \({A}_{{v}_{{\rm{B}}}}\). In a more precise calibration of the CZ gate, an optional π-rotation is applied to the control qubit (for example, Q1) to prepare it into the |0 or |1 state, followed by a Ramsey sequence on the target qubit (Q2) interleaved by an exchange pulse. The amplitude is precisely tuned to bring Q2 completely out of phase (by 180°) between the two measurements (Extended Data Fig. 4d, e). The phase θ2 is determined such that the phase of Q2 changes by zero (π) when Q1 is in the state |0 (|1), corresponding to CZ = diag(1, 1, 1, −1) in the standard basis. The same measurement is then performed again with Q2 as the control qubit and Q1 as the target qubit to determine θ1 (ref. 16).

In such a ‘conventional’ calibration procedure of the CZ gate, we notice that the two qubits experience different conditional phases (Extended Data Fig. 4). We believe that this effect is caused by off-resonant driving from the optional π-rotation on the control qubit. Similar effects can also affect the calibration of the phase crosstalk from single-qubit gates.

This motivates us to use the results from GST as feedback to adjust the gate parameters. The error generators not only describe the total errors of the gates but also distinguish Hamiltonian errors (coherent errors) from stochastic errors (incoherent errors). We use the information on seven different Hamiltonian errors (IX, IY, XI, YI, ZI, IZ and ZZ) of each gate to correct all 11 gate parameters (Extended Data Fig. 5), except \({f}_{{{\rm{Q}}}_{1}}\) and \({f}_{{{\rm{Q}}}_{2}}\), for which calibrations using standard Ramsey sequences are sufficient. For single-qubit gates, tXY1 and tXY2 are adjusted according to the IX, IY, XI and YI errors. The phases ϕ11, ϕ12, ϕ21 and ϕ22 are adjusted according to the ZI and IZ errors. For the CZ gate, θ1 and θ2 are adjusted according to the ZI and IZ errors, and \({A}_{{v}_{{\rm{B}}}}\) is adjusted according to the ZZ error. The adjusted gates are then used in a new GST experiment.

Theoretical model

In this section, we describe the theoretical model used for the fitting, the pulse optimization and the numerical simulations. The dynamics of two electron spins in the (1,1) charge configuration can be accurately described by an extended Heisenberg model21

$$H=g{\mu }_{{\rm{B}}}{{\bf{B}}}_{1}\cdot {{\bf{S}}}_{1}+g{\mu }_{{\rm{B}}}{{\bf{B}}}_{2}\cdot {{\bf{S}}}_{2}+hJ({{\bf{S}}}_{1}\cdot {{\bf{S}}}_{2}-\frac{1}{4}),$$
(2)

with \({{\bf{S}}}_{j}={({X}_{j},{Y}_{j},{Z}_{j})}^{{\rm{T}}}/2,\) where Xj, Yj and Zj are the single-qubit Pauli matrices acting on spin j = 1, 2, μB the Bohr’s magneton, g ≈ 2 the g-factor in silicon and h is Planck’s constant. The first and second terms describe the interaction of the electron spin in dot 1 and dot 2 with the magnetic fields \({{\bf{B}}}_{j}={({B}_{x,j},0,{B}_{z,j})}^{{\rm{T}}}\) originating from the externally applied field and the micromagnet. The transverse components Bx,j induce spin-flips, thus, single-qubit gates if modulated resonantly via EDSR. For later convenience, we define the resonance frequencies by \(h{f}_{{Q}_{1}}=g{\mu }_{{\rm{B}}}{B}_{z,1}\) and \(h{f}_{{Q}_{2}}=g{\mu }_{{\rm{B}}}{B}_{z,2}\), and the energy difference between the qubits ΔEz = B(Bz,2 − Bz,1). The last term in the Hamiltonian of equation (2) describes the exchange interaction J between the spins in neighbouring dots. The exchange interaction originates from the overlap of the wave functions through virtual tunnelling events and is, in general, a nonlinear function of the applied barrier voltage vB. We note that vB determines the compensation pulses applied to LP and RP for virtual barrier control. We model J as an exponential function31,32

$$J({v}_{{\rm{B}}})={J}_{{\rm{res}}}{{\rm{e}}}^{2\alpha {v}_{{\rm{B}}}},$$
(3)

where Jres ≈ 20–100 kHz is the residual exchange interaction during idle and single-qubit operations and α is the lever arm. In general, the magnetic fields \({{\bf{B}}}_{j}\) depend on the exact position of the electron. We include this in our model \({B}_{z,j}\to {B}_{z,j}({v}_{{\rm{B}}})={B}_{z,j}(0)+{\beta }_{j}{v}_{{\rm{B}}}^{\gamma },\) where βj accounts for the impact of the barrier voltage on the resonance frequency of qubit j. The transition energies described in the main text are now given by diagonalizing the Hamiltonian from equation (2) and computing the energy difference between the eigenstates corresponding to the computational basis states {|00, |01, |10, |11} (ref. 44). We have

$$h{f}_{{{\rm{Q}}}_{1}}({{\rm{Q}}}_{2}=|0\rangle )= {\mathcal E} (|10\rangle )- {\mathcal E} (|00\rangle ),$$
(4)
$$h{f}_{{{\rm{Q}}}_{1}}({{\rm{Q}}}_{2}=|1\rangle )= {\mathcal E} (|11\rangle )- {\mathcal E} (|01\rangle ),$$
(5)
$$h{f}_{{{\rm{Q}}}_{2}}({{\rm{Q}}}_{1}=|0\rangle )= {\mathcal E} (|01\rangle )- {\mathcal E} (|00\rangle ),$$
(6)
$$h{f}_{{{\rm{Q}}}_{2}}({{\rm{Q}}}_{1}=|1\rangle )= {\mathcal E} (|11\rangle )- {\mathcal E} (|10\rangle ),$$
(7)

where \( {\mathcal E} (|\xi \rangle )\) denotes the eigenenergy of eigenstate |ξ and |0 = |↓ is defined by the magnetic field direction.

In the presence of noise, qubits start to lose information. In silicon, charge noise and nuclear noise are the dominating sources of noise. In the absence of two-qubit coupling and correlated charge noise, both qubits decohere largely independently of each other, giving rise to a decoherence time set by the interaction with the nuclear spins and charge noise coupling to the qubit via intrinsic and artificial (via the inhomogeneous magnetic field) spin–orbit interaction. We describe this effect by \({f}_{{Q}_{1}}\to {f}_{{Q}_{1}}+\delta {f}_{{Q}_{1}}\) and \({f}_{{Q}_{2}}\to {f}_{{Q}_{2}}+\delta {f}_{{Q}_{2}}\), where \(\delta {f}_{{Q}_{1}}\) and \(\delta {f}_{{Q}_{2}}\) are the  single-qubit  frequency fluctuations. Charge noise can additionally affect both qubits via correlated frequency shifts and the exchange interaction through the barrier voltage, which we model as vB → vB + δvB. In the presence of finite exchange coupling, one can define four distinct, pure dephasing times, each corresponding to the dephasing of a single qubit with the other qubit in a specific basis state. In a quasistatic approximation, the four dephasing times are then given by

$${T}_{2}^{\ast }({{\rm{Q}}}_{1}({{\rm{Q}}}_{2}=|0\rangle ))=\frac{1}{\sqrt{2}{\rm{\pi }}\sqrt{{\left[\frac{d(h{f}_{{{\rm{Q}}}_{1}}({{\rm{Q}}}_{2}=|0\rangle )))}{d{v}_{{\rm{B}}}}\right]}^{2}\delta {v}_{{\rm{B}}}^{2}+{\left[\frac{d(h{f}_{{{\rm{Q}}}_{1}}({{\rm{Q}}}_{2}=|0\rangle ))}{dh{f}_{{{\rm{Q}}}_{1}}}\right]}^{2}\delta {f}_{{{\rm{Q}}}_{1}}^{2}+{\left[\frac{d(h{f}_{{{\rm{Q}}}_{1}}({{\rm{Q}}}_{2}=|0\rangle ))}{dh{f}_{{{\rm{Q}}}_{2}}}\right]}^{2}\delta {f}_{{{\rm{Q}}}_{2}}^{2}}},$$
(8)
$${T}_{2}^{\ast }({{\rm{Q}}}_{1}({{\rm{Q}}}_{2}=|1\rangle ))=\frac{1}{\sqrt{2}{\rm{\pi }}\sqrt{{\left[\frac{d(h{f}_{{{\rm{Q}}}_{1}}({{\rm{Q}}}_{2}=|1\rangle )))}{d{v}_{{\rm{B}}}}\right]}^{2}\delta {v}_{{\rm{B}}}^{2}+{\left[\frac{d(h{f}_{{{\rm{Q}}}_{1}}({{\rm{Q}}}_{2}=|1\rangle ))}{dh{f}_{{{\rm{Q}}}_{1}}}\right]}^{2}\delta {f}_{{{\rm{Q}}}_{1}}^{2}+{\left[\frac{d(h{f}_{{{\rm{Q}}}_{1}}({{\rm{Q}}}_{2}=|1\rangle ))}{dh{f}_{{{\rm{Q}}}_{2}}}\right]}^{2}\delta {f}_{{{\rm{Q}}}_{2}}^{2}}},$$
(9)
$${T}_{2}^{\ast }({{\rm{Q}}}_{2}({{\rm{Q}}}_{1}=|0\rangle ))=\frac{1}{\sqrt{2}{\rm{\pi }}\sqrt{{\left[\frac{d(h{f}_{{\rm{Q2}}}({{\rm{Q}}}_{1}=|0\rangle )))}{d{v}_{{\rm{B}}}}\right]}^{2}\delta {v}_{{\rm{B}}}^{2}+{\left[\frac{d(h{f}_{{{\rm{Q}}}_{2}}({{\rm{Q}}}_{1}=|0\rangle ))}{dh{f}_{{{\rm{Q}}}_{1}}}\right]}^{2}\delta {f}_{{{\rm{Q}}}_{1}}^{2}+{\left[\frac{d(h{f}_{{{\rm{Q}}}_{2}}({{\rm{Q}}}_{1}=|0\rangle ))}{dh{f}_{{{\rm{Q}}}_{2}}}\right]}^{2}\delta {f}_{{{\rm{Q}}}_{2}}^{2}}},$$
(10)
$${T}_{2}^{\ast }({{\rm{Q}}}_{2}({{\rm{Q}}}_{1}=|1\rangle ))=\frac{1}{\sqrt{2}{\rm{\pi }}\sqrt{{\left[\frac{d(h{f}_{{{\rm{Q}}}_{2}}({{\rm{Q}}}_{1}=|1\rangle )))}{d{v}_{{\rm{B}}}}\right]}^{2}\delta {v}_{{\rm{B}}}^{2}+{\left[\frac{d(h{f}_{{{\rm{Q}}}_{2}}({{\rm{Q}}}_{1}=|1\rangle ))}{dh{f}_{{{\rm{Q}}}_{1}}}\right]}^{2}\delta {f}_{{{\rm{Q}}}_{1}}^{2}+{\left[\frac{d(h{f}_{{{\rm{Q}}}_{2}}({{\rm{Q}}}_{1}=|1\rangle ))}{dh{f}_{{{\rm{Q}}}_{2}}}\right]}^{2}\delta {f}_{{{\rm{Q}}}_{2}}^{2}}}.$$
(11)

Fitting qubit frequencies and dephasing times

The transition energies in equations (4)–(7) are fitted simultaneously to the measured results from the Ramsey experiment (see Fig. 3a). For the fitting, we use the NonLinearModelFit function from the software Mathematica with the least squares method. The best fits yield the following parameters: α = 12.1 ± 0.05 V−1, β1 = −2.91 ± 0.11 MHz Vγ, β2 = 67.2 ± 0.63 MHz Vγ, γ = 1.20 ± 0.01 and Jres = 58.8 ± 1.8 kHz.

The dephasing times in equations (8)–(11) are fitted simultaneously to the measured results from the Ramsey experiment (see Fig. 3c) using the same method. The best fits yield the following parameters: δvB = 0.40 ± 0.01 mV, \(\delta {f}_{{Q}_{1}}=11\pm 0.1{\rm{kHz}}\) and \(\delta {f}_{{Q}_{2}}=24\pm 0.7{\rm{kHz}}\).

Numerical simulations

For all numerical simulations, we solve the time-dependent Schrödinger equation

$${\rm{i}}\hbar \frac{{\rm{d}}}{{\rm{d}}t}|{\rm{\psi }}(t)\rangle =H|{\rm{\psi }}(t)\rangle $$
(12)

and iteratively compute the unitary propagator according to

$$U(t+\Delta t)={{\rm{e}}}^{-\frac{{\rm{i}}}{\hbar }H(t+\Delta t)}U(t),$$

where \(\hbar =h/(2{\rm{\pi }})\) is the reduced Planck’s constant. Here H(t + Δt) is discretized into N segments of length Δt such that H(t) is constant in the time interval [t, t + Δt]. All simulations are performed in the rotating frame of the external magnetic field (Bz,1 + Bz,2)/2 and neglecting the counter-rotating terms, making the so-called rotating-wave approximation. This allows us to choose Δt = 10 ps as a sufficiently small time step.

For the noise simulations, we included classical fluctuations of \({f}_{{Q}_{1}}\to {f}_{{Q}_{1}}+\delta {f}_{{Q}_{1}}\), \({f}_{{Q}_{2}}\to {f}_{{Q}_{2}}+\delta {f}_{{Q}_{2}}\) and vB → vB + δvB. We assume the noise coupling to the resonance frequencies \(\delta {f}_{{Q}_{1}}\) and \(\delta {f}_{{Q}_{2}}\) to be quasistatic and assume 1/f noise for vB, which we describe by its spectral density \({S}_{{v}_{{\rm{B}}}}(\omega )=\delta {v}_{{\rm{B}}}/\omega \), where ω is the angular frequency. To compute time traces of the fluctuation, we use the approach introduced in refs. 45,46 to generate time-correlated time traces. The fluctuations are discretized into N segments with time Δt such that δvB(t) is constant in the time interval [t, t + Δt), with the same Δt as above. Consequently, fluctuations that are faster than \({f}_{{\rm{\max }}}=\frac{1}{\Delta t}\) are truncated.

CZ gate

We realize a universal CZ = diag(1, 1, 1, −1) gate by adiabatically pulsing the exchange interaction using a carefully designed pulse shape. Starting from equation (2), the full dynamics can be projected on the odd-parity space spanned by |01 and |10. The entangling exchange gate is reduced in this subspace to a global phase shift, thus, the goal is to minimize any dynamics inside the subspace. Introducing a new set of Pauli operators in this subspace σx = |0110| + |1001|, σy = −i|0110| + i|1001| and σz = |0101| − |1010|, we find

$${H}_{{\rm{sub}}}(t)=\frac{1}{2}(-hJ({v}_{{\rm{B}}}(t))+\Delta {E}_{z}\,{\sigma }_{z}+hJ({v}_{{\rm{B}}}(t))\,{\sigma }_{x}).$$
(14)

In order to investigate the adiabatic behaviour, it is convenient to switch into the adiabatic frame defined by \({U}_{{\rm{ad}}}={{\rm{e}}}^{-\frac{{\rm{i}}}{2}{\tan }^{-1}\left(\frac{hJ({v}_{{\rm{B}}}(t))}{\Delta {E}_{z}}\right){\sigma }_{y}}.\)The Hamiltonian accordingly transforms as

$${H}_{{\rm{ad}}}={U}_{{\rm{ad}}}^{\dagger }(t){H}_{{\rm{sub}}}(t){U}_{{\rm{ad}}}(t)-{\rm{i}}\hbar {U}_{{\rm{ad}}}^{\dagger }(t){\dot{U}}_{{\rm{ad}}}(t)$$
(15)
$$\approx \frac{1}{2}(-hJ({v}_{{\rm{B}}}(t))+\varDelta {E}_{z}{\sigma }_{z}-\frac{{h}^{2}\dot{J}}{2{\rm{\pi }}\Delta {E}_{z}}{\sigma }_{y}),$$
(16)

where the first term is unaffected and describes the global phase accumulation due to the exchange interaction, the second term describes the single-qubit phase accumulations and the last term, \(f(t)={h}^{2}\dot{J}/(4{\rm{\pi }}\Delta {E}_{z})\), describes the diabatic deviation proportional to the derivative of the exchange pulse. From equation (15) and equation (16), we assumed a constant ΔEz(t) ≈ ΔEz and hJ(t)  ΔEz. The transition probability from state |↑↓ to |↓↑ using a pulse of length tp is then given by34

$${P}_{|\uparrow \downarrow \rangle \to |\downarrow \uparrow \rangle }\approx {|{\int }_{0}^{{t}_{{\rm{p}}}}f(t){{\rm{e}}}^{-\frac{{\rm{i}}}{\hbar }\Delta {E}_{z}t}{\rm{d}}t|}^{2}$$
(17)
$$\propto {S}_{{\rm{s}}}(f(t)).$$
(18)

From the first to the second line, we identify the integral by the (short-timescale) Fourier transform, allowing us to describe the spin-flip error probability by the energy spectral density Ss of the input signal f(t). Minimizing such errors is, therefore, identical to minimizing the energy spectral density of a pulse, a well-known and solved problem from classical signal processing and statistics. Optimal shapes are commonly referred to as window functions W(t) due to their property of restricting the spectral resolution of signals. A high-fidelity exchange pulse is consequently given by J(0) = J(tp) and

$${\int }_{0}^{{t}_{{\rm{p}}}}{\rm{d}}tJ({v}_{{\rm{B}}}(t))=1/4,$$
(19)

while setting \(J(t)={A}_{{v}_{{\rm{B}}}}W(t){J}_{{\rm{res}}}\) (ref. 34), with a scaling factor \({A}_{{v}_{{\rm{B}}}}\) that is to be determined. In this work, we have chosen the cosine window

$$W(t)=\frac{1}{2}\left[1-\,\cos \left(\frac{2{\rm{\pi }}t}{{t}_{{\rm{p}}}}\right)\right]$$
(20)

from signal processing, which has a high spectral resolution. The amplitude \({A}_{{v}_{{\rm{B}}}}\) follows from condition equation (19). For a pulse length of tp = 100 ns and a cosine pulse shape, we find \({A}_{{v}_{{\rm{B}}}}{J}_{{\rm{res}}}=10.06\,{\rm{MHz}}\). As explained in the main text, owing to the exponential voltage-exchange relation, the target pulse shape for J(t) must be converted to a barrier gate pulse, following47

$${v}_{{\rm{B}}}(t)=\frac{1}{2\alpha }\,\log ({A}_{{v}_{{\rm{B}}}}\,W(t)).$$
(21)

Our numerical simulations predict an average gate infidelity 1 − Fgate < 10−6 without noise and 1 − F = 0.22 × 10−3 with the inclusion of noise through the fluctuations \(\delta {f}_{{Q}_{1}}\), \(\delta {f}_{{Q}_{2}}\) and δvB, discussed in the previous section. The measured PTMs reveal much higher rates of incoherent errors, which we attribute to drifts in the barrier voltage on a timescale much longer than the timescale on which \(\delta {f}_{{Q}_{1}}\), \(\delta {f}_{{Q}_{2}}\) and δvB were determined.

Gate-set tomography analysis

We designed a GST experiment using the gate set \(\{{\rm{I}},{{\rm{X}}}_{{{\rm{Q}}}_{1}},{{\rm{Y}}}_{{{\rm{Q}}}_{1}},{{\rm{X}}}_{{{\rm{Q}}}_{2}},{{\rm{Y}}}_{{{\rm{Q}}}_{2}},{\rm{CZ}}\},\)where I is a 100-ns idle gate, \({{\rm{X}}}_{{{\rm{Q}}}_{1}}\) (\({{\rm{Y}}}_{{{\rm{Q}}}_{1}}\)) and \({{\rm{X}}}_{{{\rm{Q}}}_{2}}\) (\({{\rm{Y}}}_{{{\rm{Q}}}_{2}}\)) are single-qubit π/2 gates with rotation axis \(\hat{x}\) (\(\hat{y}\)) on Q1 and Q2, with durations of 150 ns and 200 ns, respectively, and CZ = diag(1, 1, 1, −1). A classic two-qubit GST experiment consists of a set of germs designed to amplify all types of error in the gate set when repeated and a set of 36 fiducials composed by the 11 elementary operations \(\{{\rm{null}},{{\rm{X}}}_{{{\rm{Q}}}_{1}},{{\rm{X}}}_{{{\rm{Q}}}_{1}}{{\rm{X}}}_{{{\rm{Q}}}_{1}},{{\rm{X}}}_{{{\rm{Q}}}_{1}}{{\rm{X}}}_{{{\rm{Q}}}_{1}}{{\rm{X}}}_{{{\rm{Q}}}_{1}},\) \({{\rm{Y}}}_{{{\rm{Q}}}_{1}},{{\rm{Y}}}_{{{\rm{Q}}}_{1}}{{\rm{Y}}}_{{{\rm{Q}}}_{1}}{{\rm{Y}}}_{{{\rm{Q}}}_{1}},{{\rm{X}}}_{{{\rm{Q}}}_{2}},{{\rm{X}}}_{{{\rm{Q}}}_{2}}{{\rm{X}}}_{{{\rm{Q}}}_{2}},{{\rm{X}}}_{{{\rm{Q}}}_{2}}{{\rm{X}}}_{{{\rm{Q}}}_{2}}{{\rm{X}}}_{{{\rm{Q}}}_{2}},{{\rm{Y}}}_{{{\rm{Q}}}_{2}},{{\rm{Y}}}_{{{\rm{Q}}}_{2}}{{\rm{Y}}}_{{{\rm{Q}}}_{2}}{{\rm{Y}}}_{{{\rm{Q}}}_{2}}\}\)required to carry out quantum process tomography of the germs48. We use a set of 16 germs \(\{{\rm{I}},{{\rm{X}}}_{{{\rm{Q}}}_{1}},{{\rm{Y}}}_{{{\rm{Q}}}_{1}},{{\rm{X}}}_{{{\rm{Q}}}_{2}},{{\rm{Y}}}_{{{\rm{Q}}}_{2}},{\rm{CZ}},{{\rm{X}}}_{{{\rm{Q}}}_{1}}{{\rm{Y}}}_{{{\rm{Q}}}_{1}},{{\rm{X}}}_{{{\rm{Q}}}_{2}}{{\rm{Y}}}_{{{\rm{Q}}}_{2}},{{\rm{X}}}_{{{\rm{Q}}}_{1}}{{\rm{X}}}_{{{\rm{Q}}}_{1}}{{\rm{Y}}}_{{{\rm{Q}}}_{1}},{{\rm{X}}}_{{{\rm{Q}}}_{2}}{{\rm{X}}}_{{{\rm{Q}}}_{2}}\)\({{\rm{Y}}}_{{{\rm{Q}}}_{2}},{{\rm{X}}}_{{{\rm{Q}}}_{2}}{{\rm{Y}}}_{{{\rm{Q}}}_{2}}{\rm{C}}{\rm{Z}},{{\rm{C}}{\rm{Z}}{\rm{X}}}_{{{\rm{Q}}}_{2}}{{\rm{X}}}_{{{\rm{Q}}}_{1}}{{\rm{X}}}_{{{\rm{Q}}}_{1}},{{\rm{X}}}_{{{\rm{Q}}}_{1}}{{\rm{X}}}_{{{\rm{Q}}}_{2}}{{\rm{Y}}}_{{{\rm{Q}}}_{2}}{{\rm{X}}}_{{{\rm{Q}}}_{1}}{{\rm{Y}}}_{{{\rm{Q}}}_{2}}{{\rm{Y}}}_{{{\rm{Q}}}_{1}},{{\rm{X}}}_{{{\rm{Q}}}_{1}}{{\rm{Y}}}_{{{\rm{Q}}}_{2}}{{\rm{X}}}_{{{\rm{Q}}}_{2}}{{\rm{Y}}}_{{{\rm{Q}}}_{1}}{{\rm{X}}}_{{{\rm{Q}}}_{2}}{{\rm{X}}}_{{{\rm{Q}}}_{2}},\) \({{\rm{C}}{\rm{Z}}{\rm{X}}}_{{{\rm{Q}}}_{2}}{{\rm{Y}}}_{{{\rm{Q}}}_{1}}{{\rm{C}}{\rm{Z}}{\rm{Y}}}_{{{\rm{Q}}}_{2}}{{\rm{X}}}_{{{\rm{Q}}}_{1}},{{\rm{Y}}}_{{{\rm{Q}}}_{1}}{{\rm{X}}}_{{{\rm{Q}}}_{1}}{{\rm{Y}}}_{{{\rm{Q}}}_{2}}{{\rm{X}}}_{{{\rm{Q}}}_{1}}{{\rm{X}}}_{{{\rm{Q}}}_{2}}{{\rm{X}}}_{{{\rm{Q}}}_{1}}{{\rm{Y}}}_{{{\rm{Q}}}_{1}}{{\rm{Y}}}_{{{\rm{Q}}}_{2}}\}\)(ref. 35).   Note  that  the null gate is the instruction for doing nothing in zero time, different from the idle gate. Simple errors such as errors in the rotation angle of a particular gate can be amplified by simply repeating the same gate. More complicated errors such as tilts in rotation axes can only be amplified by a combination of different gates. The germs and fiducials are then compiled into GST sequences, such that each sequence consists of two fiducials interleaved by a single germ or power of germs35 (as illustrated in Fig. 2a). The GST sequences are classified by their germ powers into lengths L = 1, 2, 4, 8, 16…, where a sequence of length n consists of n gates plus the fiducial gates. We note that the sequences used in GST are shorter than the sequences involved in other methods to self-consistently estimate the gate performance, such as randomized benchmarking. As a result, GST suffers less from drift in qubit frequencies and readout windows induced by long sequences of microwave bursts.

After the execution of all sequences, a maximum-likelihood estimation is performed to estimate the process matrices of each gate in the gate set and the SPAM probabilities. We use the open source pyGSTi Python package49,50 to perform the maximum-likelihood estimation, as well as to design an optimized GST experiment by eliminating redundant circuits and to provide statistical error bars by computing all involved Hessians. The circuit optimization allows us to perform GST with a maximum sequence length Lmax = 16 using 1,685 different sequences in total. The pyGSTi package quantifies the Markovian-model violation of the experimental data, counting the number of standard deviations exceeding their expectation values under the χ2 hypothesis50. This model violation is internally translated into a more accessible goodness ratio from 0 to 5, with 5 being the best49, where we obtain a 4 out of 5 rating, indicating remarkably small deviations from expected results. The total number of standard deviations exceeding the expected results for each L, as well as the contribution of each sequence to this number, can be found in the pyGSTi report, along with the supporting data.

From the GST experiment, we have extracted the PTM \({ {\mathcal M} }_{\exp }\) describing each gate in our gate set \(\{{\rm{I}},{{\rm{X}}}_{{{\rm{Q}}}_{1}},{{\rm{Y}}}_{{{\rm{Q}}}_{1}},{{\rm{X}}}_{{{\rm{Q}}}_{2}},{{\rm{Y}}}_{{{\rm{Q}}}_{2}},{\rm{CZ}}\}\). The PTM is isomorphically related to the conventionally used χ matrix describing a quantum process. A completely positive, trace-preserving, two-qubit PTM has 240 parameters describing the process. To obtain insight into the errors of the gates in the experiment, we first compute the error in the PTM given by \(E={ {\mathcal M} }_{\exp }{ {\mathcal M} }_{{\rm{ideal}}}^{-1}\), where we have adapted the convention to add the error after the ideal gate. The average gate fidelity is then conveniently given by

$${F}_{{\rm{gate}}}=\frac{{\rm{Tr}}({ {\mathcal M} }_{\exp }^{-1}{ {\mathcal M} }_{{\rm{ideal}}})+d}{d(d+1)}.$$
(22)

It is related to the entanglement fidelity via \(1-{F}_{{\rm{ent}}}=\frac{d+1}{d}(1-{F}_{{\rm{gate}}})\) (ref. 51), where d is the dimension of the two-qubit Hilbert space. Although the PTM \( {\mathcal M} \) perfectly describes the errors, it is more intuitive to analyse the corresponding error generator \( {\mathcal L} =\,\log (E)\) of the process30. The error generator \( {\mathcal L} \) relates to the error PTM E in a similar way as a Hamiltonian H relates to a unitary operation U = e−iH. The error generator can be separated into several blocks. A full discussion about the error generator can be found in ref. 30. In this work, we have used the error generator to distinguish the dynamics originating from coherent Hamiltonian errors, which can be corrected by adjusting gate parameters (see Extended Data Fig. 5), and from noisy/stochastic dynamics, which cannot be corrected easily. The coherent errors can be extracted by projecting \( {\mathcal L} \) onto the 4 × 4-dimensional Hamiltonian space H. In the Hilbert–Schmidt space, the Hamiltonian projection is given by30

$${H}_{mn}=-\frac{{\rm{i}}}{{d}^{2}}{\rm{Tr}}[({P}_{m}^{{\rm{T}}}\otimes {P}_{n}^{{\rm{T}}}\otimes {{\bf{1}}}_{d}-{{\bf{1}}}_{d}\otimes {P}_{m}\otimes {P}_{n}){ {\mathcal L} }_{{\rm{\sup }}}],$$
(23)

where \({ {\mathcal L} }_{{\rm{\sup }}}\) is the error generator in Liouville superoperator form, Pm {I, X, Y, Z} are the extended Pauli matrices with m, n = 0, 1, 2, 3, 1d is the d-dimensional identity matrix and d = 4 is the dimension of the two-qubit Hilbert space. To improve the calibration of our gate set, we use the information of seven different Hamiltonian errors (IX, IY, XI, YI, ZI, IZ and ZZ). To estimate coherent Hamiltonian errors and incoherent stochastic errors, two new metrics are considered30: the Jamiołkowski probability

$${\epsilon }_{{\rm{J}}}( {\mathcal L} )=-{\rm{Tr}}({\rho }_{{\rm{J}}}( {\mathcal L} )|\Psi \rangle \langle \Psi |)),$$
(24)

which describes the amount of incoherent error in the process, and the Jamiołkowski amplitude

$${\theta }_{{\rm{J}}}( {\mathcal L} )={\Vert (1-|\Psi \rangle \langle \Psi |){\rho }_{{\rm{J}}}( {\mathcal L} )|\Psi \rangle \Vert }_{2},$$
(25)

which approximately describes the amount of coherent Hamiltonian errors (Extended Data Table 1). Here \({\rho }_{{\rm{J}}}( {\mathcal L} )=( {\mathcal L} \otimes {1}_{{d}^{2}})[|\varPsi \rangle \langle \Psi |]\) is the Jamiołkowski state and |Ψ is a maximally entangling four-qubit state that originates from the relation of quantum processes to states in a Hilbert space twice the dimension via the Choi–Jamiołkowski isomorphism52. For small errors, the average gate infidelity can be approximated by30

$$1-{F}_{{\rm{gate}}}=\frac{d}{d+1}[{\epsilon }_{{\rm{J}}}( {\mathcal L} )+{\theta }_{{\rm{J}}}{( {\mathcal L} )}^{2}].$$
(26)

For a comparison of the performance of the single-qubit gates with previous experiments reporting single-qubit gate fidelities, we compute the fidelities projected to the single-qubit space from the PTMs or the error generators. In Fig. 2 and Extended Data Fig. 2, single-qubit gate fidelities are estimated by projecting the PTMs onto the corresponding subspace. Let \({{\mathscr{P}}}_{j}\) be the projector on the subspace of qubit j, then the fidelity is given by

$${F}_{{\rm{sub}}}=\frac{{\rm{Tr}}({{\mathscr{P}}}_{j}{ {\mathcal M} }_{\exp }^{-1}{{\mathscr{P}}}_{j}{ {\mathcal M} }_{{\rm{ideal}}})+(d/2)}{(d/2)((d/2)+1)}.$$
(27)

Error bars for the fidelity projected to the subspace are computed using standard error propagation of the confidence intervals of \({ {\mathcal M} }_{\exp }\) provided by the pyGSTi package. A more optimistic estimation for the fidelities in the single-qubit subspace is given by projecting the error generators instead of the PTMs.

Variational quantum eigensolver

We follow the approach of ref. 36 to using the VQE algorithm to compute the ground-state energy of molecular hydrogen, after mapping this state onto the state of two qubits. We include this information here for completeness. The Hamiltonian of a molecular system in atomic units is

$$\begin{array}{c}H=-\sum _{i}\frac{{\nabla }_{{{\bf{R}}}_{i}}^{2}}{2{M}_{i}}-\sum _{j}\frac{{\nabla }_{{{\bf{r}}}_{j}}^{2}}{2}-\sum _{i,j}\frac{{Q}_{i}}{|{{\bf{R}}}_{i}-{{\bf{r}}}_{j}|}\\ +\sum _{i,j > i}\frac{{Q}_{i}{Q}_{j}}{|{{\bf{R}}}_{i}-{{\bf{R}}}_{j}|}+\sum _{i,j > i}\frac{1}{|{{\bf{r}}}_{i}-{{\bf{r}}}_{j}|},\end{array}$$
(28)

where \({{\bf{R}}}_{i}\), Mi and Qi are the position, mass and charge, respectively, of the ith nuclei and \({{\bf{r}}}_{j}\) is the position of the jth electron. The first two sums describe the kinetic energies of the nuclei and electrons, respectively. The last three sums describe the Coulomb repulsion between nuclei and electrons, nuclei and nuclei, and electrons and electrons, respectively. As we are primarily interested in the electronic structure of the molecule, and nuclear masses are a few orders of magnitude larger than the electron masses, the nuclei are treated as static point charges under the Born–Oppenheimer approximation. Consequently, the electronic Hamiltonian can be simplified to

$${H}_{{\rm{e}}}=-\sum _{i}\frac{{\nabla }_{{{\bf{r}}}_{i}}^{2}}{2}-\sum _{i,j}\frac{{Q}_{i}}{|{{\bf{R}}}_{i}-{{\bf{r}}}_{j}|}+\sum _{i,j > i}\frac{1}{|{{\bf{r}}}_{i}-{{\bf{r}}}_{j}|}.$$
(29)

Switching into the second-quantization representation, described by fermionic creation and annihilation operators, \({a}_{p}^{\dagger }\) and aq, acting on a finite basis, the Hamiltonian becomes

$${H}_{{\rm{e}}}=\sum _{pq}{h}_{pq}{a}_{p}^{\dagger }{a}_{q}+\sum _{pqrs}{h}_{pqrs}{a}_{p}^{\dagger }{a}_{q}^{\dagger }{a}_{r}{a}_{s},$$
(30)

where p, q, r and s label the corresponding basis states. The antisymmetry under exchange is retained through the anticommutation relation of the operators. The weights of the two sums are given by the integrals

$${h}_{pq}=\int {\rm{d}}{\boldsymbol{\sigma }}{{\rm{\psi }}}_{p}^{\ast }({\boldsymbol{\sigma }})(\frac{{\nabla }_{{{\bf{r}}}_{i}}^{2}}{2}-\sum _{i}\frac{{Q}_{i}}{|{{\bf{R}}}_{i}-{\bf{r}}|}){\psi }_{q}({\boldsymbol{\sigma }}),$$
(31)
$${h}_{pqrs}=\int {\rm{d}}{{\boldsymbol{\sigma }}}_{1}{\rm{d}}{{\boldsymbol{\sigma }}}_{2}\frac{{{\rm{\psi }}}_{p}^{\ast }({{\boldsymbol{\sigma }}}_{1}){{\rm{\psi }}}_{q}^{\ast }({{\boldsymbol{\sigma }}}_{2}){{\rm{\psi }}}_{s}({{\boldsymbol{\sigma }}}_{1}){{\rm{\psi }}}_{r}({{\boldsymbol{\sigma }}}_{2})}{|{{\bf{r}}}_{1}-{{\bf{r}}}_{2}|},$$
(32)

where \({{\boldsymbol{\sigma }}}_{i}=({{\bf{r}}}_{i},{s}_{i})\) is a multi-index describing the position \({{\bf{r}}}_{i}\) and the spin si of electron i. Such a second-quantized molecular Hamiltonian can be mapped onto qubits using the Jordan–Wigner (JW) or the BK transformation5. The JW transformation directly encodes the occupation number (0 or 1) of the ith spin orbital into the state (|0 or |1) of the ith qubit. The number of qubits required after JW transformation is, thus, the same as the number of spin orbitals that are of interest. The BK transformation, on the other hand, encodes the information in both the occupation number and parities, whether there is an even or odd occupation in a subset of spin orbitals.

Taking molecular hydrogen in the HF basis as an example, we are interested in investigating the bonding (|O1, |O1) and the antibonding orbital state (|O2, |O2). The initial guess of the solution is the HF state in which both electrons occupy the |O1 orbital. The JW transformation encodes the HF initial state as |1100, representing \(|{N}_{{O}_{1}\downarrow }{N}_{{O}_{1}\uparrow }{N}_{{O}_{2}\downarrow }{N}_{{O}_{2}\uparrow }\rangle \) from left to right, where \({N}_{{O}_{i}S}\) is the occupation of the OiS spin orbital with S = ↑, ↓. The BK transformation encodes the HF initial state as |1000, where the first and the third qubits (counting from the right) encode the occupation number of the first and third spin orbitals (\({N}_{{O}_{1}\uparrow }=1\) and \({N}_{{O}_{2}\uparrow }=0\)), the second qubit encodes the parity of the first two spin orbitals (\(({N}_{{O}_{1}\uparrow }+{N}_{{O}_{1}\downarrow })\) mod 2 = 0) and the fourth qubit encodes the parity of all four spin orbitals (\(({N}_{{O}_{1}\uparrow }+{N}_{{O}_{1}\downarrow }+{N}_{{O}_{2}\uparrow }+{N}_{{O}_{2}\downarrow })\) mod 2 = 0). With the standard transformation rules for fermionic creation and annihilation operators, the system Hamiltonian becomes a four-qubit Hamiltonian

$$\begin{array}{c}{H}_{{\rm{JW}}}=\,{g}_{0}{\rm{I}}+{g}_{1}{{\rm{Z}}}_{1}+{g}_{2}{{\rm{Z}}}_{2}+{g}_{3}{{\rm{Z}}}_{3}+{g}_{4}{{\rm{Z}}}_{4}\\ \,+{g}_{5}{{\rm{Z}}}_{1}{{\rm{Z}}}_{2}+{g}_{6}{{\rm{Z}}}_{1}{{\rm{Z}}}_{3}+{g}_{7}{{\rm{Z}}}_{1}{{\rm{Z}}}_{4}\\ \,\,+{g}_{8}{{\rm{Z}}}_{2}{{\rm{Z}}}_{3}+{g}_{9}{{\rm{Z}}}_{2}{{\rm{Z}}}_{4}+{g}_{10}{{\rm{Z}}}_{3}{{\rm{Z}}}_{4}\\ \,+{g}_{11}{{\rm{Y}}}_{1}{{\rm{X}}}_{2}{{\rm{X}}}_{3}{{\rm{Y}}}_{4}+{g}_{12}{{\rm{Y}}}_{1}{{\rm{Y}}}_{2}{{\rm{X}}}_{3}{{\rm{X}}}_{4}\\ \,\,+{g}_{13}{{\rm{X}}}_{1}{{\rm{X}}}_{2}{{\rm{Y}}}_{3}{{\rm{Y}}}_{4}+{g}_{14}{{\rm{X}}}_{1}{{\rm{Y}}}_{2}{{\rm{Y}}}_{3}{{\rm{X}}}_{4},\end{array}$$
(33)
$$\begin{array}{l}{H}_{{\rm{BK}}}=\,{g}_{0}{\rm{I}}+{g}_{1}{{\rm{Z}}}_{1}+{g}_{2}{{\rm{Z}}}_{2}+{g}_{3}{{\rm{Z}}}_{3}\\ \,+{g}_{4}{{\rm{Z}}}_{1}{{\rm{Z}}}_{2}+{g}_{5}{{\rm{Z}}}_{1}{{\rm{Z}}}_{3}+{g}_{6}{{\rm{Z}}}_{2}{{\rm{Z}}}_{4}\\ \,+{g}_{7}{{\rm{Z}}}_{1}{{\rm{Z}}}_{2}{{\rm{Z}}}_{3}+{g}_{8}{{\rm{Z}}}_{1}{{\rm{Z}}}_{3}{{\rm{Z}}}_{4}+{g}_{9}{{\rm{Z}}}_{2}{{\rm{Z}}}_{3}{{\rm{Z}}}_{4}\\ \,+{g}_{10}{{\rm{Z}}}_{1}{{\rm{Z}}}_{2}{{\rm{Z}}}_{3}{{\rm{Z}}}_{4}+{g}_{11}{{\rm{X}}}_{1}{{\rm{Z}}}_{2}{{\rm{X}}}_{3}\\ \,+{g}_{12}{{\rm{Y}}}_{1}{{\rm{Z}}}_{2}{{\rm{Y}}}_{3}+{g}_{13}{{\rm{X}}}_{1}{{\rm{Z}}}_{2}{{\rm{X}}}_{3}{{\rm{Z}}}_{4}+{g}_{14}{{\rm{Y}}}_{1}{{\rm{Z}}}_{2}{{\rm{Y}}}_{3}{{\rm{Z}}}_{4}.\end{array}$$
(34)

The subscripts are used to label the qubits. We see that, owing to the symmetry of the represented system in HBK, qubit 2 and qubit 4 are never flipped, allowing us to reduce the dimension of the Hamiltonian to

$$\begin{array}{c}{H}_{{\rm{BK}}}^{{\rm{reduced}}}=\,{h}_{0}{\rm{I}}+{h}_{1}{{\rm{Z}}}_{1}+{h}_{2}{{\rm{Z}}}_{2}+{h}_{3}{{\rm{Z}}}_{1}{{\rm{Z}}}_{2}\\ \,+{h}_{4}{{\rm{X}}}_{1}{{\rm{X}}}_{2}+{h}_{5}{{\rm{Y}}}_{1}{{\rm{Y}}}_{2}{h}_{0}\\ \,+{h}_{1}{\rm{ZI}}+{h}_{2}{\rm{IZ}}+{h}_{3}{\rm{ZZ}}\\ +{h}_{4}{\rm{XX}}+{h}_{5}{\rm{YY}},\end{array}$$
(35)

where qubit 1 has been relabelled as qubit 2 and qubit 3 has been relabelled as qubit 1. The HF initial state is, therefore, reduced to |01 and the Hamiltonian is rephrased to be consistent with the partial tomography expression in Fig. 5. This reduced representation requires only two qubits to simulate the hydrogen molecule. We emphasize that such a reduction of the BK Hamiltonian is not a special case for the H2 molecule but is connected to symmetry considerations to reduce the complexity of systems, in a scalable way.

VQE is a method to compute the ground-state energy of the Hamiltonian. The total energy can be directly calculated by measuring the expectation value of each Hamiltonian term. This can be done easily by partial quantum state tomography. All the expectation values are then added up with a set of weights (h0 through h5). The weights are only functions of the internuclear separation (R) and can be computed efficiently by a classical computer. Here we use the OpenFermion Python package to compute these weights38.

The main task of the quantum processor is, then, to encode the molecular spin-orbital state into the qubits. The starting point is the HF initial state, which is believed to largely overlap with the actual ground state. In order to find the actual ground state, the initial state needs to be ‘parameterized’ into an ansatz to explore a subspace of all possible states. We apply the unitary coupled cluster (UCC) theory to the parameterized ansatz state, which is used to describe many-body systems and cannot be efficiently executed on a classical computer53. The UCC operator has a format

$${U}_{{\rm{UCC}}}({\boldsymbol{\theta }})={{\rm{e}}}^{\sum _{n}({T}_{n}({\boldsymbol{\theta }})-{T}_{n}^{\dagger }({\boldsymbol{\theta }}))},$$
(36)

with

$${T}_{1}({\boldsymbol{\theta }})=\sum _{m,i}{{\boldsymbol{\theta }}}_{i}^{m}{a}_{m}^{\dagger }{a}_{i},$$
(37)
$${T}_{2}({\boldsymbol{\theta }})=\sum _{m,n,i,j}{{\boldsymbol{\theta }}}_{i,j}^{m,n}{a}_{m}^{\dagger }{a}_{n}^{\dagger }{a}_{i}{a}_{j}$$
(38)

representing single and double excitation of the electrons. The indices i and j label the occupied spin orbitals and m and n are the labels of the unoccupied spin orbitals. The vector θ is the set of all parameters to optimize. In the case of a H2 molecule, the UCC operator is transformed into a qubit operator as

$${U}_{{\rm{UCC}}}^{{\rm{BK}}}({\boldsymbol{\theta }})={{\rm{e}}}^{-{\rm{i}}\theta {\rm{XY}}},$$

where θ is a single parameter to variationally optimize.