Quantum logic with spin qubits crossing the surface code threshold – Nature

Measurement setup

The measurement setup and device are similar to those used in ref. ¹⁷. 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 v_LP/v_B = −0.081 and v_RP/v_B = 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. ¹⁷ 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 t_p is given by

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 t_XY1 and t_XY2, phase shifts caused by single-qubit gates on the addressed qubit itself ϕ₁₁ and ϕ₂₂, phase shifts caused by single-qubit gates on the unaddressed ‘victim qubit’ ϕ₁₂ and ϕ₂₁ (ϕ₁₂ is the phase shift on Q₁ induced by a gate on Q₂ and similar for ϕ₂₁), 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 θ₁ and θ₂ (in addition, we absorb into θ₁ and θ₂ 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 t_XY1 and t_XY2 are initially calibrated by an AllXY calibration protocol⁴³. The phases ϕ₁₁, ϕ₁₂, ϕ₂₁ and ϕ₂₂ are initially calibrated by measuring the phase shift of the victim qubit (Q₁ for ϕ₁₁ and ϕ₂₁; Q₂ for ϕ₂₂ and ϕ₁₂) in a Ramsey sequence interleaved by a pair of ([{{rm{X}}}_{{{rm{Q}}}_{i}},-{{rm{X}}}_{{{rm{Q}}}_{i}}]) gates on the addressed qubit (Q₁ for ϕ₁₁ and ϕ₁₂; Q₂ for ϕ₂₂ and ϕ₂₁) (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, Q₁) to prepare it into the |0⟩ or |1⟩ state, followed by a Ramsey sequence on the target qubit (Q₂) interleaved by an exchange pulse. The amplitude is precisely tuned to bring Q₂ completely out of phase (by 180°) between the two measurements (Extended Data Fig. 4d, e). The phase θ₂ is determined such that the phase of Q₂ changes by zero (π) when Q₁ 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 Q₂ as the control qubit and Q₁ as the target qubit to determine θ₁ (ref. ¹⁶).

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, t_XY1 and t_XY2 are adjusted according to the IX, IY, XI and YI errors. The phases ϕ₁₁, ϕ₁₂, ϕ₂₁ and ϕ₂₂ are adjusted according to the ZI and IZ errors. For the CZ gate, θ₁ and θ₂ 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 model²¹

with ({{bf{S}}}_{j}={({X}_{j},{Y}_{j},{Z}_{j})}^{{rm{T}}}/2,) where X_j, Y_j and Z_j 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 B_x,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 ΔE_z = gμ_B(B_z,2 − B_z,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 v_B. We note that v_B determines the compensation pulses applied to LP and RP for virtual barrier control. We model J as an exponential function^31,32

where J_res ≈ 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. ⁴⁴). We have

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 v_B → v_B + δv_B. 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

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⁻¹, β₁ = −2.91 ± 0.11 MHz V^−γ, β₂ = 67.2 ± 0.63 MHz V^−γ, γ = 1.20 ± 0.01 and J_res = 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: δv_B = 0.40 ± 0.01 mV, (delta {f}_{{Q}_{1}}=11pm 0.1{rm{kHz}}) and (delta {f}_{{Q}_{2}}=24pm 0.7{rm{kHz}}).

Numerical simulations

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

and iteratively compute the unitary propagator according to

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 (B_z,1 + B_z,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 v_B → v_B + δv_B. 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 v_B, 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 δv_B(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 = |01⟩⟨10| + |10⟩⟨01|, σ_y = −i|01⟩⟨10| + i|10⟩⟨01| and σ_z = |01⟩⟨01| − |10⟩⟨10|, we find

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

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 ΔE_z(t) ≈ ΔE_z and hJ(t) ≪ ΔE_z. The transition probability from state |↑↓⟩ to |↓↑⟩ using a pulse of length t_p is then given by³⁴

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 S_s 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(t_p) and

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

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 t_p = 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, following⁴⁷

Our numerical simulations predict an average gate infidelity 1 − F_gate < 10⁻⁶ without noise and 1 − F = 0.22 × 10⁻³ with the inclusion of noise through the fluctuations (delta {f}_{{Q}_{1}}), (delta {f}_{{Q}_{2}}) and δv_B, 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 δv_B 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 Q₁ and Q₂, 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 germs⁴⁸. 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. ³⁵).   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 germs³⁵ (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 package^49,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 L_max = 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 χ² hypothesis⁵⁰. This model violation is internally translated into a more accessible goodness ratio from 0 to 5, with 5 being the best⁴⁹, 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

It is related to the entanglement fidelity via (1-{F}_{{rm{ent}}}=frac{d+1}{d}(1-{F}_{{rm{gate}}})) (ref. ⁵¹), 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 process³⁰. 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. ³⁰. 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 by³⁰

where ({ {mathcal L} }_{{rm{sup }}}) is the error generator in Liouville superoperator form, P_m ∈ {I, X, Y, Z} are the extended Pauli matrices with m, n = 0, 1, 2, 3, 1_d 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 considered³⁰: the Jamiołkowski probability

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

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 isomorphism⁵². For small errors, the average gate infidelity can be approximated by³⁰

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

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. ³⁶ 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

where ({{bf{R}}}_{i}), M_i and Q_i 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

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

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

where ({{boldsymbol{sigma }}}_{i}=({{bf{r}}}_{i},{s}_{i})) is a multi-index describing the position ({{bf{r}}}_{i}) and the spin s_i of electron i. Such a second-quantized molecular Hamiltonian can be mapped onto qubits using the Jordan–Wigner (JW) or the BK transformation⁵. 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 (|O₁↑⟩, |O₁↓⟩) and the antibonding orbital state (|O₂↑⟩, |O₂↓⟩). The initial guess of the solution is the HF state in which both electrons occupy the |O₁⟩ 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 O_iS 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

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

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 H₂ 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 (h₀ through h₅). 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 weights³⁸.

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 computer⁵³. The UCC operator has a format

with

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 H₂ molecule, the UCC operator is transformed into a qubit operator as

where θ is a single parameter to variationally optimize.