

\section{Device}


\begin{figure*}[htbp] 
\center{\includegraphics[width=1\linewidth]{figures/Figure_1.eps}}
\caption{\textbf{a.} Scanning electron microscope images of the silicon quantum processor showing the quantum dot gate pattern and the micromagnet on top. \textbf{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. \textbf{c.} Pulse sequence schematic of a decoupled controlled-phase operation interleaved in a Ramsey interference sequence on Q1. \textbf{d.} Spin-up probability of Q1 after the Ramsey sequence in \textbf{c}, as a function of the detuning in the double dot potential and the total duration of the barrier voltage pulses.}
\label{fig:device}
\end{figure*}



We use a gate-defined double quantum dot in an isotopically enriched $^{28}$Si/SiGe heterostructure (Fig.~\ref{fig:device}a)~\cite{xue2020cmos}, with each dot occupied by one 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 dots~\cite{elzerman2004single}. A micromagnet on top of the device provides a magnetic field gradient enabling electric-dipole spin resonance (EDSR)~\cite{pioro2008electrically}, and separates the resonance frequencies of the qubits in the presence of an external magnetic field ($\sim$\unit[320]{mT}) to \unit[11.993]{GHz} (Q1) and \unit[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 frame~\cite{vandersypen2005nmr}. The plunger gates (LP and RP) control the chemical potentials of the quantum dots.\par

The native two-qubit gate for spin qubits utilizes the exchange interaction~\cite{loss1998quantum, petta2005coherent}, originating from the wave-function overlap of electrons in neighbouring dots. This selectively shifts the energy of the anti-parallel spin states and thus allows for an electrically pulsed adiabatic \textsc{cphase} gate~\cite{meunier2011efficient,veldhorst2015two,watson2018programmable}. The barrier gate (B) controls the tunnel coupling between the dots, allowing to precisely tune the exchange coupling from < \unit[100]{kHz} to \unit[20]{MHz}. In order to minimize the sensitivity to charge noise, we activate the exchange coupling while avoiding a tilt in the double dot potential~\cite{martins2016noise,reed2016reduced,zhang_2017_tilt} (Fig.~\ref{fig:device}a). 
This symmetric condition can be determined accurately by decoupled adiabatic exchange pulses inside a Ramsey sequence (Fig.~\ref{fig:device}c-d). 
The tunnel barrier is controlled by simultaneously pulsing gate B and compensating on LP and RP to avoid shifts of the electrochemical potentials~\cite{martins2016noise}. The detuning between quantum dots is controlled by additional offsets on 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 center (Fig.~\ref{fig:device}e).\par





\section{Conclusion}


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 logic~\cite{madzik2021precision}. Independent studies have shown spin qubit readout with a fidelity above 98\% in only a few $\mu$s~\cite{zheng2019rapid}, with further improvements underway~\cite{schaal2020fast}. With a modest effort in reducing crosstalk errors and in extending the device designs, we are optimistic that the individually demonstrated advantages of semiconductor spin qubits can be combined into a fault-tolerant and highly-integrated quantum computer. The same advances will allow us to implement more sophisticated algorithms in the NISQ era, such as solving energies involving excited states of more complex molecules.

\textbf{Acknowledgements}
We acknowledge fruitful discussions with P. Cerfontaine, C. Bureau-Oxton, M. Madzik, A. Morello, J. Helsen, B. Terhal, M. Veldhorst, and all the members of the spin qubit team, and technical assistance by O. Benningshof, M. Sarsby, R. Schouten and R. Vermeulen. This research was funded by the Dutch Ministry for Economic Affairs through the allowance for Topconsortia for Knowledge and Innovation (TKI) and the Army Research Office (ARO) under grant numbers W911NF-17-1-0274. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the ARO or the US Government. The US Government is authorized to reproduce and distribute reprints for government purposes notwithstanding any copyright notation herein.

\textbf{Author contributions} X.X. performed the experiment with help from N.S. and B.U., M.R developed the theory model and analyzed the data with X.X., N.S. fabricated the quantum dot device, A.S. and G.S. designed and grew the Si/SiGe heterostructure, X.X. and L.M.K.V. conceived the project, L.M.K.V. supervised the project, X.X., M.R. and L.M.K.V. wrote the manuscript with input from all authors.

\textbf{Competing interests}
The authors declare no competing interests.

\textbf{Data availability}
Data supporting this work are available at zenodo, https://doi.org/10.5281/zenodo.5044450.

\textbf{Code availability}
The codes used for data acquisition and processing are from the open source python packages QCoDeS, which is available at https://github.com/QCoDeS/Qcodes, QTT, which is available at https://github.com/QuTech-Delft/qtt, and PycQED, which is available at https://github.com/DiCarloLab-Delft/PycQED{\textunderscore}py3. The codes used for the design and analysis of the gate set tomography experiment are from pyGSTi, which is available at https://github.com/pyGSTio/pyGSTi. And the codes used for the design and analysis of the variational quantum eigensolver experiment are from OpenFermion, which is available at https://github.com/quantumlib/OpenFermion.

\section{Fitting qubit frequencies and dephasing times}
The transition energies Eqs.~\eqref{eq:resonance1}-\eqref{eq:resonance4} are fitted simultaneously to the measured results from the Ramsey experiment (see Fig.~\ref{fig:pulse_engineer}a). For the fitting we use the \texttt{NonLinearModelFit} function from the software Mathematica with the least square method. The best fits yield the following parameters $\alpha=\unit[12.1\pm 0.05]{V^{-1}}$, $\beta_1=\unit[-2.91\pm 0.11]{MHz^\gamma/V^\gamma}$, $\beta_2=\unit[67.2\pm 0.63]{MHz^\gamma/V^\gamma}$, and $\gamma=\unit[1.20\pm 0.01]{}$, and $J_\text{res}=\unit[58.8\pm 1.8]{kHz}$.\par

The dephasing times Eqs.~\eqref{eq:T2star_1}-\eqref{eq:T2star_4} are fitted simultaneously to the measured results from the Ramsey experiment (see Fig.~\ref{fig:pulse_engineer}c) using the same method. The best fits yield the following parameters $\delta v_B=\unit[0.40\pm 0.01]{m V}$, $\delta f_{Q1}=\unit[11\pm 0.1]{kHz}$, and $\delta f_{Q2}=\unit[24\pm 0.7]{kHz}$. 
\section{Gate calibration}

In the gate set used in this work, $\{I$, $X_{Q1}$, $Y_{Q1}$, $X_{Q2}$, $Y_{Q2}$, $\textsc{cphase} \}$, the duration of the $I$ gate and the \textsc{cphase} are set to \SI{100}{\nano\second}, and we calibrate and keep the amplitudes of the single-qubit drives fixed and in the linear response regime. The envelope of the single-qubit gates are shaped following a ``tukey'' window, as it allows adiabatic single-qubit gates with relatively small amplitudes, thus avoids distortion caused by nonlinear response.
The general tukey window of length $t_p$ is given by
\begin{align}
    W(t,r) = \begin{cases} 
      \frac{1}{2}\left[1-\cos\left(\frac{2\pi t }{r t_p}\right)\right] & 0 \leq t \leq \frac{r t_p}{2} \\
      1 & \frac{r t_p}{2} < t\leq t_p -\frac{r t_p}{2} \\
       \frac{1}{2}\left[1-\cos\left(\frac{2\pi (t_p-t) }{r t_p}\right)\right] & t_p -\frac{r t_p}{2}\leq t \leq t_p,
   \end{cases}
   \label{eq:tukey_window}
\end{align}
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_{Q1}$ and $f_{Q2}$, burst lengths for single-qubit gates $t_{XY1}$ and $t_{XY2}$, phase shifts caused by single-qubit gates on the addressed qubit itself $\phi_{11}$ and $\phi_{22}$, phase shifts caused by single-qubit gates on the unaddressed ``victim qubit'' $\phi_{12}$ and $\phi_{21}$ ($\phi_{12}$ is the phase shift on Q1 induced by a gate on Q2 and similar for $\phi_{21}$), the peak amplitude of the \textsc{cphase} gate $A_{v_B}$, and phase shifts caused by the gate voltage pulses used for \textsc{cphase} gate on the qubits $\theta_1$ and $\theta_2$ (in addition, we absorb into $\theta_1$ and $\theta_2$ the 90 degree phase shifts needed to transform $\text{diag}(1,i,i,1)$ into $\text{diag}(1,1,1,-1)$).\par
For single-qubit gates, $f_{Q1}$ and $f_{Q2}$ are calibrated by standard Ramsey sequences, which are automatically executed every two hours. The EDSR burst times $t_{XY1}$ and $t_{XY2}$ are initially calibrated by an AllXY calibration protocol~\cite{reed2013entanglement}. The phases $\phi_{11}$, $\phi_{12}$, $\phi_{21}$, and $\phi_{22}$ are initially calibrated by measuring the phase shift of the victim qubit (Q1 for $\phi_{11}$ and $\phi_{21}$; Q2 for $\phi_{22}$ and $\phi_{12}$) in a Ramsey sequence interleaved by a pair of [$X_{Qi}$, $-X_{Qi}$] gates on the addressed qubit (Q1 for $\phi_{11}$ and $\phi_{12}$; Q2 for $\phi_{22}$ and $\phi_{21}$) (Extended Data Fig.~\ref{fig:Manual_calibration}). \par
The optimal pulse design presented in Fig.~\ref{fig:pulse_engineer} gives a rough guidance of the pulse amplitude $A_{v_B}$. In a more precise calibration of the \textsc{cphase} gate, an optional $\pi$-rotation is applied to the control qubit (e.g. Q1) to prepare it into the $\ket{0}$ or $\ket{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 degree) between the two measurements (Extended Data Fig.~\ref{fig:Manual_calibration} d-e). The phase $\theta_2$ is determined such that the phase of Q2 changes by zero ($\pi$) when Q1 is in the state $\ket{0}$ ($\ket{1}$),  corresponding to the $\textsc{cphase}=\text{diag}(1,1,1,-1)$ in the standard basis. The same measurement is then performed again with the Q2 as the control qubit and Q1 as the target qubit to determine $\theta_1$~\cite{watson2018programmable}.\par
In such a ``conventional'' calibration procedure of the \textsc{cphase} gate, we notice that the two qubits experience different conditional phases (Extended Data Fig.~\ref{fig:Manual_calibration}). We believe that this effect is caused by off-resonant driving from the optional $\pi$-rotation on the control qubit. Similar effects can also affect the calibration of the phase crosstalk from single-qubit gates.\par

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 7 different Hamiltonian errors ($IX$, $IY$, $XI$, $YI$, $ZI$, $IZ$ and $ZZ$) of each gate, to correct all 11 gate parameters, except $f_{Q1}$ and $f_{Q2}$, 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 $\phi_{11}$, $\phi_{12}$, $\phi_{21}$, and $\phi_{22}$ are adjusted according to the $ZI$ and $IZ$ errors. For the \textsc{cphase} gate, $\theta_1$ and $\theta_2$ are adjusted according to the $ZI$ and $IZ$ errors, and $A_{v_B}$ is adjusted according to the $ZZ$ error. The adjusted gates are then used in a new GST loop.
\section{Measurement setup}

The measurement setup and device are similar to the one used in Ref.~\cite{xue2020cmos}. 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 on LP and RP. Linear combinations of the voltage pulses on 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, 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.\par

Two differences between the present setup and the one in Ref.~\cite{xue2020cmos} 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; 2) a second AWG of the same model is connected to gate B with its clock synchronized with the first AWG.\par

\section{Numerical simulations}
\label{sec:num_sim}

For all numerical simulations,  we solve the time-dependent Schr{\"o}dinger equation 
\begin{align}
    i\hbar \ket{\psi(t)} = H \ket{\psi(t)}
\end{align}
and iteratively compute the unitary propagator according to
\begin{align}
    U(t+\Delta t) = e^{-\frac{i}{\hbar} H(t+\Delta t)}U(t).
\end{align}
Here, $H(t+\Delta t$ is discretized into $N$ segments of length $\Delta t$ such that $H(t)$ is constant in the time-interval $\left[t,t+\Delta t\right]$.
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 (RWA). This allows us to chose $\Delta t=\unit[10]{ps}$ as a sufficiently small time step. 

For the noise simulations, we included classical fluctuations of $f_{Q1}\rightarrow f_{Q1} + \delta f_{Q1}$, $f_{Q2}\rightarrow f_{Q2} + \delta f_{Q2}$, and $v_B\rightarrow v_B + \delta v_B$. We assume the noise coupling to the resonance frequencies $\delta f_{Q1}$ and $\delta f_{Q2}$ to be quasistatic and assume $1/f$ noise for $v_B$ which we describe by its spectral density $S_{\delta v_B}(\omega)= \delta v_B/\omega$. To compute time traces of the fluctuation we use the approach introduced in Refs.~\cite{yang_achieving_2019,koski_strong_2020} to generate time-correlated time traces. The fluctuations are discretized into $N$ segments with time $\Delta t$ such that $\delta v_B(t)$ is constant in the time interval $\left[t,t+\Delta t\right)$, with the same $\Delta t$ as above. Consequently, fluctuations which are faster than $f_\text{max}=\frac{1}{\Delta t}$ are truncated. 
\section{\textsc{cphase} gate}
We realize a universal $\textsc{cphase}=\text{diag}(1,1,1,-1)$ gate by adiabatically pulsing the exchange interaction using a carefully designed pulse shape. Starting from Eq.~\eqref{eq:ham_matrix}, the full dynamics can be projected on the odd-parity space spanned by $\ket{01}$ and $\ket{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
$\sigma_x=\ket{01}\bra{10}+\ket{10}\bra{01}$, $\sigma_y=-i \ket{01}\bra{10}+i\ket{10}\bra{01}$, and $\sigma_z=\ket{01}\bra{01}-\ket{10}\bra{10}$, we find 
\begin{align}
H_{\text{sub}}(t)=\frac{1}{2}\big( - J\left(v_B(t)\right) + \Delta E_z \,\sigma_z+J\left(v_B(t)\right)\, \sigma_x\big).
\label{eq:HheisSub}
\end{align}
In order to investigate the adiabatic behaviour, it is convenient to switch into the adiabatic frame defined by $U_\text{ad}=e^{-\frac{i}{2} \tan^{-1}\left(\frac{J\left(v_B(t)\right)}{\Delta E_z}\right)\sigma_y}$. The Hamiltonian accordingly transforms as
\begin{align}
H_{\text{ad}}&=U_\text{ad}^\dagger(t) H_{\text{sub}}(t) U_\text{ad}(t) -i \hbar U_\text{ad}^\dagger(t)\dot{U}_\text{ad}(t)\label{eq:Ham_ramp_part1}\\
&\approx \frac{1}{2}\bigg(- J\left(v_B(t)\right) +  \Delta E_z\,\sigma_z
-\frac{\hbar\dot{J}}{\Delta E_z}\sigma_y \bigg),
\label{eq:Ham_ramp_part2}
\end{align}
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)=\hbar\dot{J}/(2\Delta E_z)$ describes the diabatic deviation proportional to the derivative of the exchange pulse. From Eq.~\eqref{eq:Ham_ramp_part1} to Eq.~\eqref{eq:Ham_ramp_part2} we assumed a constant $\Delta E_z(t)\approx \Delta E_z$, and $J(t)\ll \Delta E_z$. The transition probability from state $\ket{\uparrow\downarrow}$ to $\ket{\downarrow\uparrow}$ using a pulse of length $t_p$ is then given by~\cite{martinis_fast_2014}
\begin{align}
P_{\ket{\uparrow\downarrow}\rightarrow\ket{\downarrow\uparrow}} &\approx \left|\int_{0}^{t_p} f(t) e^{-\frac{i}{\hbar} \Delta E_z t }\, dt \right|^2 \label{eq:ESD_part1} \\
&= S\big(f(t)\big).
\label{eq:ESD}
\end{align}
From the first to the second line we replaced the integral by the (short time-scale) Fourier transform, allowing us to describe the spin-flip error probability by the energy spectral density (ESD) of the input signal $f(t)$. Minimizing such errors is therefore identical to minimizing the ESD 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 to restrict the spectral resolution of signals. 
A high-fidelity exchange pulse is consequently given by $J(0)=J(t_p)$ and 
\begin{align}
    \int_0^{t_p}dt J\left(v_B(t)\right)/\hbar=\pi/2,
    \label{eq:cphase_condition}
\end{align}
while setting $J(t)= A_{v_B} W(t) J_\text{res}$~\cite{martinis_fast_2014}, with a scaling factor $A_{v_B}$ that is to be determined. In this work, we have chosen the cosine window
\begin{align}
    W(t)=\frac{1}{2}\left[1-\cos\left(\frac{2\pi t }{t_p}\right)\right],
\end{align}
from signal processing which has a high spectral resolution. The amplitude $A_{v_B}$ follows from condition Eq.~\eqref{eq:cphase_condition}. For a pulse length of $t_p=\unit[100]{ns}$ and a cosine pulse shape we find $A_{v_B} J_\text{res}=\unit[10.06]{MHz}$. Our numerical simulations predict an average gate infidelity $1-F_\text{gate}<10^{-6}$ without noise and $1-F=0.22\times 10^{-3}$ with the inclusion of noise through the fluctuations $\delta f_{Q1}$, $\delta f_{Q2}$, and $\delta v_B$ discussed in section~\ref{sec:num_sim}. 
As explained in the main text, due to the exponential voltage-exchange relation the target pulse shape for J(t) must be converted to a barrier gate pulse, following~\cite{russ_soon_to_be_appear}
\begin{align}
    v_B(t)=\frac{1}{2\alpha}\text{log}\left(A_{v_B}\, W(t)\right).
\end{align}
The numerical simulations with the fitted noise parameters in the simplified model from Section~\ref{sec:num_sim} predict a gate fidelity above 99.97\%. The measured PTMs reveal significantly 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_{Q1}$, $\delta f_{Q2}$, and $\delta v_B$ were determined. 
 
\section{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 well-described by an extended Heisenberg model~\cite{loss1998quantum}
\begin{align}
H &=J\,(\vv{S}_1\cdot \vv{S}_2-\frac{1}{4}) +  g\mu_B\vv{B}_1\cdot\vv{S}_1+  g\mu_B\vv{B}_2\cdot\vv{S}_2,
\label{eq:ham_matrix}
\end{align}
with $\vv{S}_j = \hbar(X_j,Y_j,Z_j)^T/2$, where $X_j,Y_j,Z_j$ are the single-qubit Pauli-matrices acting on spin $j=1,2$, $\mu_B$ the Bohr's magneton, $g\approx 2$ the g-factor in silicon, and $\hbar=h/(2\pi)$ the reduced Planck constant. The first and second term describe the interaction of the electron spin in dot 1 and dot 2 with the magnetic fields $\vv{B}_j=(B_{x,j},0,B_{z,j})^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 $hf_{Q1}=g\mu_B B_{z,1}$ and $hf_{Q2}=g\mu_B B_{z,2}$, and the energy difference between the qubits $\Delta E_z= g\mu_B (B_{z,2}- B_{z,1})$. The last term in the Hamiltonian of Eq.~\eqref{eq:ham_matrix} describes the exchange interaction $J$ between the spins in neighboring dots. The exchange interaction originates from the overlap of the wave-functions through virtual tunneling events and is in general a non-linear function of the applied barrier voltage $v_B$. 
We note that $v_B$ determines the compensation pulses on LP and RP for virtual barrier control.
We model $J$ as an exponential function~\cite{cerfontaine2020high,pan_resonant_2020}
\begin{align}
J(v_B) = J_\text{res} e^{2\alpha v_B},
\label{eq:J_sim_expression}
\end{align}
where $J_\text{res}\approx\unit[20-100]{kHz}$ is the residual exchange interaction during idle and single-qubit operations and $\alpha$ the lever arm.
In general the magnetic fields $\vv{B}_j$ depend on the exact position of the electron. We include this in our model $B_{z,j}\rightarrow B_{z,j}(v_B)= B_{z,j}(0) + \beta_j v_B^\gamma$, where $\beta_j$ amounts 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 Hamiltonian from Eq.~\eqref{eq:ham_matrix} and computing the energy difference between the eigenstates corresponding to the computational basis states $\lbrace\ket{00},\ket{01},\ket{10},\ket{11}\rbrace$~\cite{russ_high-fidelity_2018}. We have
\begin{align}
    hf_\text{Q1 (Q2$=\ket{0}$)} &= \mathcal{E}(\ket{10}) - \mathcal{E}(\ket{00}),\label{eq:resonance1} \\
    hf_\text{Q1 (Q2$=\ket{1}$)} &= \mathcal{E}(\ket{11}) - \mathcal{E}(\ket{01}),\label{eq:resonance2} \\
    hf_\text{Q2 (Q1$=\ket{0}$)} &= \mathcal{E}(\ket{01}) - \mathcal{E}(\ket{00}),\label{eq:resonance3} \\
    hf_\text{Q2 (Q1$=\ket{1}$)} &= \mathcal{E}(\ket{11}) - \mathcal{E}(\ket{10}),\label{eq:resonance4}
\end{align}
where $\mathcal{E}(\ket{\xi})$ denotes the eigenenergy of eigenstate $\ket{\xi}$ and $\ket{0}=\ket{\downarrow}$ is defined by the magnetic field direction.\par

In the presence of noise, qubits start to loose information. In silicon, charge noise and nuclear noise are the dominating noise sources. 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_{Q1}\rightarrow f_{Q1} + \delta f_{Q1}$ and $f_{Q2}\rightarrow f_{Q2} + \delta f_{Q2}$, where $\delta f_{Q1}$ and $\delta f_{Q1}$ are the single-qubit frequency fluctuations. Charge noise additionally can affect both qubits via correlated frequency shifts and the exchange interaction through the barrier voltage, which we model as $v_B\rightarrow v_B + \delta 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
\begin{widetext}

\begin{align}
    T_2^\star(\text{Q1 (Q2$=\ket{0}$)}) &= \frac{1}{\sqrt{2}\pi\sqrt{\left[\frac{d\left(hf_\text{Q1 (Q2$=\ket{0}$)})\right)}{d v_B}\right]^2\delta v_B^2
    + \left[\frac{d\left(hf_\text{Q1 (Q2$=\ket{0}$)}\right)}{d hf_{Q1}}\right]^2 \delta f_{Q1}^2
    +\left[\frac{d\left(hf_\text{Q1 (Q2$=\ket{0}$)}\right)}{d hf_{Q2}}\right]^2\delta f_{Q2}^2 }},\label{eq:T2star_1} \\
    T_2^\star(\text{Q1 (Q2$=\ket{1}$)}) &= \frac{1}{\sqrt{2}\pi\sqrt{\left[\frac{d\left(hf_\text{Q1 (Q2$=\ket{1}$)})\right)}{d v_B}\right]^2\delta v_B^2
    + \left[\frac{d\left(hf_\text{Q1 (Q2$=\ket{1}$)}\right)}{d hf_{Q1}}\right]^2 \delta f_{Q1}^2
    +\left[\frac{d\left(hf_\text{Q1 (Q2$=\ket{1}$)}\right)}{d hf_{Q2}}\right]^2\delta f_{Q2}^2 }}, \label{eq:T2star_2} \\
    T_2^\star(\text{Q2 (Q1$=\ket{0}$)}) &= \frac{1}{\sqrt{2}\pi\sqrt{\left[\frac{d\left(hf_\text{Q2 (Q1$=\ket{0}$)})\right)}{d v_B}\right]^2\delta v_B^2
    + \left[\frac{d\left(hf_\text{Q2 (Q1$=\ket{0}$)}\right)}{d hf_{Q1}}\right]^2 \delta f_{Q1}^2
    +\left[\frac{d\left(hf_\text{Q2 (Q1$=\ket{0}$)}\right)}{d hf_{Q2}}\right]^2\delta f_{Q2}^2 }}, \label{eq:T2star_3} \\
    T_2^\star(\text{Q2 (Q1$=\ket{1}$)}) &= \frac{1}{\sqrt{2}\pi\sqrt{\left[\frac{d\left(hf_\text{Q2 (Q1$=\ket{1}$)})\right)}{d v_B}\right]^2\delta v_B^2
    + \left[\frac{d\left(hf_\text{Q2 (Q1$=\ket{1}$)}\right)}{d hf_{Q1}}\right]^2 \delta f_{Q1}^2
    +\left[\frac{d\left(hf_\text{Q2 (Q1$=\ket{1}$)}\right)}{d hf_{Q2}}\right]^2\delta f_{Q2}^2 }}.
    \label{eq:T2star_4} 
\end{align}
\end{widetext}
\section{VQE}

We follow the approach of Ref.~\cite{hempel2018quantum} to using the VQE algorithm to compute the ground state 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 ($\hbar = 1$) reads 
\begin{align}
    H = &-\sum_{i} \frac{\nabla_{R_i}^2}{2M_i} -\sum_{j} \frac{\nabla_{r_j}^2}{2} - \sum_{i,j} \frac{Q_i}{\abs{R_i - r_j}}\nonumber\\
    &+\sum_{i,j>i} \frac{Q_iQ_j}{\abs{R_i - R_j}} + \sum_{i,j>i} \frac{1}{\abs{r_i - r_j}},
\end{align}
where $R_i$, $M_i$ and $Q_i$ are the position, mass and charge of the $i$-th nuclei, and $r_j$ is the position of the $j$-th 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. Consequentially, the electronic Hamiltonian can be simplified to
\begin{align}
    H_e = &-\sum_{i} \frac{\nabla_{r_i}^2}{2} - \sum_{i,j} \frac{Q_i}{\abs{R_i - r_j}} + \sum_{i,j>i} \frac{1}{\abs{r_i - r_j}}.
\end{align}
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
\begin{align}
    H_e = \sum_{pq}h_{pq}a_p^{\dagger}a_q + \sum_{pqrs}h_{pqrs}a_p^{\dagger}a_q^{\dagger}a_ra_s.
\end{align}
The anti-symmetry under exchange is retained through the anti-commutation relation of the operators.
The weights of the two sums are given by the integrals 
\begin{align}
    &h_{pq} = \int d\sigma\psi_p^*(\sigma)(\frac{\nabla_{r_i}^2}{2} - \sum_{i} \frac{Q_i}{\abs{R_i - r}})\psi_q^*(\sigma), \\
    &h_{pqrs} = \int d\sigma_1d\sigma_2\frac{\psi_p^*(\sigma_1)\psi_q^*(\sigma_2)\psi_s(\sigma_1)\psi_r(\sigma_2)}{\abs{r_1 - r_2}}.
\end{align}
Such a second-quantized molecular Hamiltonian can be mapped onto qubits using the Jordan-Wigner (JW) or the Bravyi-Kitaev (BK) transformation~\cite{mcardle2020quantum}. The JW transformation directly encodes the occupation number (0 or 1) of the $i$-th spin-orbital into the state ($\ket{0}$ or $\ket{1}$) of the $i$-th 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.\par

Taking molecular hydrogen in the Hartree-Fock basis as an example, we are interested in investigating the bonding ($\ket{O_1\uparrow}$, $\ket{O_1\downarrow}$) and the anti-bonding orbital state ($\ket{O_2\uparrow}$, $\ket{O_2\downarrow}$). The initial guess of the solution is the Hartree-Fock (HF) state in which both electrons occupy the $\ket{O_1}$ orbital. The JW transformation encodes the HF initial state as $\ket{0011}$, representing $\ket{N_{O_2\downarrow} N_{O_2\uparrow} N_{O_1\downarrow} N_{O_1\uparrow}}$ from left to right, where $N_{O_iS}$ is the occupation of the $O_iS$ spin-orbital with $S=\uparrow,\downarrow$. The BK transformation encodes the HF initial state as $\ket{0001}$, where the first and the third qubit (counting from the right) encode the occupation number of the first and third spin-orbital ($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{align}
    H_{JW} = &~g_0I + g_1Z_1 + g_2Z_2 + g_3Z_3 + g_4Z_4\nonumber\\
             &~+ g_5Z_1Z_2 + g_6Z_1Z_3 + g_7Z_1Z_4\nonumber\\ 
             &~+ g_8Z_2Z_3 + g_9Z_2Z_4 + g_{10}Z_3Z_4\nonumber\\
             &~+ g_{11}Y_1X_2X_3X_4 + g_{12}Y_1Y_2X_3X_4\nonumber\\
             &~+ g_{13}X_1X_2Y_3Y_4 + g_{14}X_1Y_2Y_3Y_4,\\
    H_{BK} = &~g_0I + g_1Z_1 + g_2Z_2 + g_3Z_3\nonumber\\
             &~+ g_4Z_1Z_2 + g_5Z_1Z_3 + g_6Z_2Z_4\nonumber\\
             &~+ g_7Z_1Z_2Z_3 + g_8Z1Z_3Z_4 + g_{9}Z_2Z_3Z_4\nonumber\\
             &~+ g_{10}Z_1Z_2Z_3Z_4 + g_{11}X_1Z_2X_3\nonumber\\ 
             &~+ g_{12}Y_1Z_2Y_3 + g_{13}Y_1Z_2Y_3Z_4.
\end{align}
We see that due 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
\begin{align}
    H_{BK}^\text{reduced} = &~h_0I + h_1Z_1 + h_2Z_2 + h_3Z_1Z_2 \nonumber\\
                        &~+ h_4X_0X_1 + h_5Y_0Y_1.
\end{align}
This reduced representation requires only two qubits to simulate the hydrogen molecule. The HF initial state thus becomes $\ket{01}$. We emphasize that such a reduction of the BK Hamiltonian is not a special case for $H_2$ molecule but is connected to symmetry considerations to reduce the complexity of systems, in a scalable way.\par


The variational quantum eigensolver (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_0$ through $h_5$). The weights are only functions of the internuclear separation ($R$) and can be computed efficiently by a classical computer. Here we use the \texttt{OpenFermion} python package to compute these weights~\cite{mcclean2020openfermion}.\par
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 widely believed to be a powerful approach and cannot be efficiently executed on a classical computer~\cite{taube2006new}. The UCC operator has a format
\begin{align}
    U_{UCC}(\vv{\theta}) = e^{\sum_n (T_n(\vv{\theta}) - T_n^{\dagger}(\vv{\theta}))},
\end{align}
with 
\begin{align}
    &T_1(\vv{\theta}) = \sum_{m,i} \vv{\theta}_i^m a_m^{\dagger}a_i,\\
    &T_2(\vv{\theta}) = \sum_{m,n,i,j} \vv{\theta}_{i,j}^{m,n} a_m^{\dagger}a_n^{\dagger}a_ia_j
\end{align}
representing single-excitation and double-excitation of the electrons. The indices $i, j$ label the occupied spin-orbitals and $m, n$ are the labels of the unoccupied spin-orbitals. The vector $\vv{\theta}$ is the set of all parameters to optimize. In the case of a $H_2$ molecule, the UCC operator is transformed into a qubit operator as
\begin{align}
    U_{UCC}^{BK}(\vv{\theta}) = e^{-i\theta XY},
\end{align}
where $\theta$ is a single parameter to variationally optimize.
\section{Gate set tomography analysis}
We designed a customized gate set tomography (GST) experiment using the gate set $\{ I$, $X_{Q1}$, $Y_{Q1}$, $X_{Q2}$, $Y_{Q2}$, $\textsc{cphase}\}$, where $I$ is a $\unit[100]{ns}$ idle gate, $X_{Q1}$ ($Y_{Q1}$) and $X_{Q2}$ ($Y_{Q2}$) are single qubit $\pi/2$ gates on qubit 1 and 2 with rotation axis $\hat{x}$ ($\hat{y}$), and $\textsc{cphase}=\text{diag}(1,1,1,-1)$. A classic two-qubit GST experiment consists of a set of germs designed to amplify any error rate in the sequence when repeated,  and a set of 36 fiducials composed by the 11 elementary operations $\{\textit{null}$, $X_{Q1}$, $X_{Q1}X_{Q1}$, $X_{Q1}X_{Q1}X_{Q1}$, $Y_{Q1}$, $Y_{Q1}Y_{Q1}Y_{Q1}$, $X_{Q2}$, $X_{Q2}X_{Q2}$, $X_{Q2}X_{Q2}X_{Q2}$, $Y_{Q2}$, $Y_{Q2}Y_{Q2}Y_{Q2}\}$ required to do quantum process tomography of the germs~\cite{greenbaum2015introduction}. Note, that the $\textit{null}$ gate is the instruction for doing nothing in zero time, different from the idle gate. 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.~\ref{fig:1QGST}a of the main text)~\cite{nielsen2020gate}. The GST sequences are classified by their germ powers into lengths $L=1,2,4,8,16\cdots$, where a sequence of length $n$ consists of $n$ gates plus the fiducial gates. After the execution of all sequences a maximum-likelihood estimation (MLE) is performed to estimate the process matrices of each gate in the gate set and the SPAM probabilities. We use the open source \texttt{pyGSTi} python package~\cite{nielsen_python_2019,nielsen_probing_2020} to perform the MLE, as well as to design a reduced 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_\text{max}=16$ using 1685 different sequences in total.
The \texttt{pyGSTi} package quantifies the Markovian-model violation of the experimental data counting the number of standard deviations exceeding their expectation values under the $\chi^2$ hypothesis~\cite{nielsen_probing_2020}. This model violation is internally translated into a more accessible goodness ratio from $0-5$ with 5 being the best~\cite{nielsen_python_2019}, where we get a 4 out of 5 rating indicating remarkably small deviations from expected results.

From the gate set tomography experiment, we have extracted the Pauli transfer matrix (PTM) $\mathcal{M}_\text{exp}$ describing each gate in our gate set $\{I$,  $X_{Q1}$, $Y_{Q1}$, $X_{Q2}$, $Y_{Q2}$, $\textsc{cphase}\}$. The PTM is isomorphically related to the conventionally used $\chi$-matrix describing a quantum process. A completely positive trace-preserving (CPTP) two-qubit PTM has 240 parameters describing the process. To get insight in the errors of the gates in the experiment, we first compute the error in the PTM given by $E=\mathcal{M}_\text{exp}\mathcal{M}_\text{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
\begin{align}
F_\text{gate} = \frac{\text{Tr}(\mathcal{M}_\text{exp}^{-1}\mathcal{M}_\text{ideal})+d}{d(d+1)}.
\label{eq:averageGate}
\end{align}
It is related to the entanglement fidelity via $F_\text{ent}=\frac{d+1}{d}F_\text{gate}$~\cite{white_measuring_2007}, where $d$ is the dimension of the two-qubit Hilbert space. While the PTM $\mathcal{M}$ perfectly describes the errors, it is more intuitive to analyze the corresponding error generator $\mathcal{L}=\log(E)$ of the process~\cite{blume2021taxonomy}. The error generator relates $\mathcal{L}$ to the error PTM $E$ in a similar way as a Hamiltonian relates $H$ to a unitary operation $U=e^{-i H}$. The error generator can be separated into several blocks. A full discussion about the error generator can be found in Ref.~\cite{blume2021taxonomy}. 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 Figure~\ref{fig:GST_optimization}, from noisy/stochastic dynamics which cannot be corrected easily. The coherent errors can be extracted by projecting $\mathcal{L}$ onto the $4\times 4$-dimensional Hamiltonian space $H$. In the Hilbert-Schmidt space, the Hamiltonian projection is given by~\cite{blume2021taxonomy}
\begin{align}
    H_{mn} = -\frac{i}{d^2}\text{Tr}\left[\left(P_m^T\otimes P_n^T\otimes \boldsymbol{1}_d - \boldsymbol{1}_d\otimes P_m\otimes P_n \right)\mathcal{L}_\text{sup}\right],
\end{align}
where $\mathcal{L}_\text{sup}$ is the error generator in Liouville superoperator form, $P_m\in{I,X,Y,Z}$ are the extended Pauli matrices with $m,n=0,1,2,3$, $\boldsymbol{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 7 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~\cite{blume2021taxonomy}; the Jamiolkowski probability 
\begin{align}
    \epsilon_J(\mathcal{L}) = -\text{Tr}(\rho_J(\mathcal{L})\ket{\Psi}\bra{\Psi})),
    \label{eq:Jam_prob}
\end{align}
which describes the amount of incoherent error in the process, and the Jamiolkowski amplitude
\begin{align}
    \theta_J(\mathcal{L}) = ||(1-\ket{\Psi}\bra{\Psi})\rho_J(\mathcal{L})\ket{\Psi}||_2,
    \label{eq:Jam_amp}
\end{align} which approximately describes the amount of coherent Hamiltonian errors (Extended Data Table.~\ref{tab:GST_results}). Here, $\rho_J(\mathcal{L}) = (\mathcal{L}\otimes\boldsymbol{1}_{d^2})[\ket{\Psi}\bra{\Psi}]$ is the Jamiolkowski state and $\ket{\Psi}$ is a maximally entangling four-qubit state which originates from the relation of quantum processes to states in a Hilbert space twice the dimension via the Choi-Jamiolkowski isomorphism~\cite{jamiolkowski_linear_1972}. For small errors, the average gate infidelity can be approximated by~\cite{blume2021taxonomy}
\begin{align}
    1-F_\text{gate}=\frac{d}{d+1}\left[\epsilon_J(\mathcal{L})+\theta_J(\mathcal{L})^2\right].
\end{align}
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.~\ref{fig:1QGST} and Extended Data Fig.~\ref{fig:Subspace_PTM}, single-qubit gate fidelities are estimated by projecting the PTMs onto corresponding subspace. Let $\mathcal{P}_{j}$ be the projector on the subspace of qubit $j$ then the fidelity is given by
\begin{align}
F_\text{sub} = \frac{\text{Tr}(\mathcal{P}_{j}\mathcal{M}_\text{exp}^{-1}\mathcal{P}_{j}\mathcal{M}_\text{ideal})+(d/2)}{(d/2)((d/2)+1)}.
\label{eq:averageGate_sub}
\end{align}
Error bars for the fidelity projected to the subspace are computed using standard error propagation of the confidence intervals of $\mathcal{M}_\text{exp}$ provided by the \texttt{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.



\section*{Methods}
\input{Measurement_setup}
\input{Gate_calibration}
\input{Theoretical_model}
\input{Fitting}
\input{Numerical_simulations}
\input{Optimal_pulse_shapes}
\input{gate_set_tomography}
\input{VQE}
\vfill

\input{Extended_figures}
\clearpage

