Universal logical operations in a silicon quantum processor

Device fabrication

The qubit device was fabricated on an isotopically purified ²⁸Si buffer layer (20 nm thick), epitaxially grown at 380 °C via ultrahigh-vacuum electron-beam deposition on a p-type natural silicon substrate (50–100 Ω cm). This layer maintains a ²⁹Si concentration below 130 ppm, effectively isolating qubits from the nuclear spin bath of the natural silicon substrate and enhancing spin coherence times. Using STM-assisted hydrogen depassivation lithography at 77 K (Scienta Omicron INFINITY SPM Lab), we deterministically patterned the hydrogen-passivated silicon surface with atomic precision. The exposed regions were then dosed with phosphine (PH₃) at 5 × 10⁻⁷ mbar for 5 min (room temperature), followed by thermal incorporation at 330 °C for 1 min to activate phosphorus donors. This process creates three tunnel-coupled quantum dot structures (P-atom clusters) positioned 17.0 nm away from a single-electron transistor charge sensor for individual qubit addressability and read-out. Subsequently, the doped nanostructure was encapsulated with a 20 nm epitaxial ²⁸Si capping layer, grown at 250 °C (0.7 nm h⁻¹ growth rate). This low-temperature epitaxial process ensures minimal dopant diffusion while maintaining crystalline perfection. After removing the device from the STM system, a polymethyl methacrylate resist and electron-beam lithography were used to define a series of small holes (200 nm in diameter) aligned with the phosphorus patches. This step was followed by reactive ion etching and aluminium lead deposition. Through these holes, the leads penetrated the dopant layer, establishing direct electrical contact between the metal and the phosphorus patches. The ohmic contacts between the buried P-dopant device and aluminium electrodes were formed through silicon vias. Then, an aluminium microwave antenna was fabricated atop the device, separated by a 10 nm atomic-layer-deposited Al₂O₃ dielectric layer to prevent current leakage and enable high-fidelity spin manipulation.

Estimation of uncertainties

All reported uncertainties correspond to 1 s.d. from the mean, estimated through 2,000 Monte Carlo bootstrap resampling trials. Specifically, we model the measured counts for each spin state (for example, \(\left|\Downarrow \Downarrow \Downarrow \Downarrow \right\rangle\) to \(\left|\Uparrow \Uparrow \Uparrow \Uparrow \right\rangle\)) across N measurement repetitions as following a multinomial distribution. To clarify the representation of the uncertainties, we provide the following illustrative example: the value 96.5(20)% corresponds to 96.5 ± 2.0%, where the number in parentheses denotes the s.d.

[[4, 2, 2]] detection code

The logical code words \(\left|{{\rm{L}}}_{1}{{\rm{L}}}_{{\rm{2}}}\right\rangle\) in the computational basis are

The stabilizers performing the syndrome extraction of the physical phase-flip and bit-flip errors are \({\widehat{S}}^{X}={X}_{1}{X}_{2}{X}_{3}{X}_{4}\) and \({\widehat{S}}^{Z}={Z}_{1}{Z}_{2}{Z}_{3}{Z}_{4}\), respectively. All logical basis states in equation (1) have even parity, while any single physical qubit errors lead to a state with odd parity.

The correspondence between logical operators and physical operators is as follows: the logical operation X_LI_L corresponds to applying X gates on physical qubits 1 and 3 (X₁I₂X₃I₄). Similarly, I_LX_L maps to X₁X₂I₃I₄. By contrast, the Z_LI_L operation corresponds to Z₁Z₂I₃I₄, while I_LZ_L translates to Z₁I₂Z₃I₄. The logical Hadamard operation H_LH_L requires H gates on all four physical qubits (H₁H₂H₃H₄), followed by swapping of physical qubits 2 and 3 (SWAP₂₃). Finally, the logical CNOT_L is implemented as a SWAP₁₂ operation, where the logical L₁ is the control qubit.

Compared with the logical operations above, the implementation of the S_LI_L gate is less intuitive. The S_LI_L gate is implemented by applying physical S gates to the first two physical qubits, followed by a CZ gate between them (Fig. 3a). The physical S gates accumulate the target phase on the logical qubit L₁, while the CZ gate cancels the extra phase accumulated on physical qubits 1 and 2. For example, when the two physical S gates are applied for the first two physical qubits, with the initial state being \({\left|{\psi }_{0}\right\rangle }_{{\rm{L}}}\,=\,{\left|00\right\rangle }_{{\rm{L}}}+{\left|10\right\rangle }_{{\rm{L}}}\), the corresponding physical qubit state becomes \(\left|{\psi }_{1}\right\rangle =\left|0000\right\rangle +{{\rm{e}}}^{i{{\uppi }}}\left|1111\right\rangle +\,{{\rm{e}}}^{i{{\uppi}}/2}\left|0101\right\rangle +{{\rm{e}}}^{i{{\uppi }}/2}\left|1010\right\rangle\). Then, the CZ gate is applied to correct the phase accumulated on the \(\left|1111\right\rangle\) state. As a result, the target logical qubit state is obtained: \({\left|{\psi }_{2}\right\rangle }_{{\rm{L}}}={\left|00\right\rangle }_{{\rm{L}}}+{{\rm{e}}}^{i{{\uppi }}/2}{\left|10\right\rangle }_{{\rm{L}}}\).

Quantum state tomography

The density matrix of a four-qubit state can be expressed as follows:

where c_i,j,k,l denotes the expansion coefficients of the density matrix and \({M}_{i},{M}_{j},{M}_{k},{M}_{l}\in \left\{I,X,Y,Z\right\}\) are the Pauli operators acting on the qubits i, j, k and l. There are 4⁴ = 256 linear independent operators in total, from I ⊗ I ⊗ I ⊗ I, X ⊗ I ⊗ I ⊗ I, …to Z ⊗ Z ⊗ Z ⊗ Z.

To reconstruct the density matrix, 255 expectation values are to be determined in the experiment (except for c_1,1,1,1 = 1). For qubit j in the I base measurement, no prerotation is applied after the target state \(\left|\psi \right\rangle\) is prepared; for qubit k in the X base measurement, a \({R}_{y}^{k}(-{{\uppi }}/2)\) prerotation is applied; a \({R}_{x}^{l}({{\uppi }}/2)\) prerotation is applied for qubit l in the Y base measurement; and we use I base projection outcomes for the Z base calculation rather than performing R_x(π) pretotation to any qubit, to prevent rotation error⁵². All the above prerotations are realized by NMR pulses, which are conditional on the electron spin being in the \(\left|\downarrow \right\rangle\) state.

After that, the maximum likelihood estimation method is used to restrict the measured density matrix to be physical (Hermitian, positive semidefinite and trace preserving) while providing the closest measurement probabilities. According to the Cholesky decomposition, a physical density matrix can be written as \(\rho ={T}^{\dagger }T/{\rm{Tr}}({T}^{\dagger }T)\). The matrix T is a complex lower triangle matrix, which has 2^2D (D = 4 for four qubits) parameters t₁, t₂, …t₂₅₆. The loss function is defined as follows:

where P_i are the measured probabilities projected at basis \(\left|{\psi }_{i}\right\rangle\). We minimize this loss function to obtain the final and physical density matrix. All the fidelity in this paper is calculated by \(F={(\mathrm{Tr}\sqrt{\sqrt{{\rho }^{\mathrm{ideal}}}\rho \sqrt{{\rho }^{\mathrm{ideal}}}})}^{2}\), where ρ^ideal is the corresponding ideal density matrix.

Logical density matrix reconstruction

The complete quantum state tomography reconstructs the physical density matrix, which generally contains components outside the logical subspace defined by the stabilizers \({\widehat{S}}^{X}\) and \({\widehat{S}}^{Z}\). The parity projection is fundamentally enabled by the projection operator (Π + I^⊗4)/2, where Π represents the stabilizers of the [[4, 2, 2]] code, and I^⊗4 denotes the four-qubit identity operator. We can project the physical density matrix into three subspaces defined by the +1 eigenspace of \({\widehat{S}}^{X}\), the +1 eigenspace of \({\widehat{S}}^{Z}\) and their simultaneous eigenspace. The logical density matrix is then obtained after applying the projection.

Quantum process tomography

A quantum operation can be described by the operator-sum representation, expressed as

where the operation elements {α_i} satisfy the completeness condition:

To experimentally determine the operators α_i, we adopt an equivalent representation using a fixed basis \(\{{\widetilde{\alpha }}_{i}\}\) spanning the space of d² complex matrices (for N-qubit systems, d = 2^N). This gives the χ-matrix representation:

where the d² × d² matrix χ fully characterizes the quantum process. For the single-qubit system, the standard basis operators are chosen as follows:

Experimentally determining an arbitrary single-qubit operation requires measuring 12 parameters through input states \(\left|0\right\rangle\), \(\left|1\right\rangle\), \(\left|+\right\rangle =(\left|0\right\rangle +\left|1\right\rangle )/\sqrt{2}\) and \(\left|+i\right\rangle =(\left|0\right\rangle +i\left|1\right\rangle )/\sqrt{2}\). Quantum state tomography on the output states yields four matrices:

These correspond to \({\rho }_{j}^{{\prime} }={\mathcal{E}}({\rho }_{j})\) with \({\rho }_{1}=\left|0\right\rangle \left\langle 0\right|\), ρ₂ = ρ₁X, ρ₃ = Xρ₁ and ρ₄ = Xρ₁X. The linear decomposition of \({\mathcal{E}}({\rho }_{j})\) in the basis states ρ_k leads to the following relation:

where coefficients \({\beta }_{jk}^{mn}\) are determined by the basis states and operators. For the single-qubit case, these relations are encapsulated in the transformation matrix:

The χ-matrix is then constructed through block matrix operations:

The process fidelity definition we used here is

where χ^ideal is the ideal process matrix.

Classical optimizer

For the classical optimizer, we employed the Nelder–Mead algorithm⁵³ for optimization, which is a simplex-based direct-search method suitable for unconstrained minimization of objective functions. Compared with other direct-search methods such as gradient descent that may exhibit better convergence on smooth functions, Nelder–Mead demonstrates robustness for noisy functions and can explore nearby valleys, thereby avoiding local minima trapping and overcoming non-smoothness of objective functions.

We improved the Nelder–Mead algorithm to suit our specific experimental situation. To enforce the algorithm’s angular domain constraint of [−π, π), we applied periodic boundary conditions by computationally wrapping the optimization parameter θ (Fig. 5a) into this interval during each iterative optimization step. Furthermore, we optimized the Nelder–Mead algorithm’s reflection, expansion, contraction and shrink coefficients to achieve convergence to the global minimum with fewer iterations.

Clifford fitting

To reduce the impact of noise in near-term quantum computations, error mitigation techniques process noisy results to more accurately estimate ideal expectation values⁵⁴. Clifford fitting is also referred to as Clifford data regression^50,55. The core idea is to learn a correction function, f^CF, using training data generated from efficiently simulatable Clifford circuits. For these circuits, both the noisy expectation values (X^noisy) measured on hardware and the exact expectation values (X^exact) computed classically are obtained. The function f^CF is determined by fitting these data pairs. Subsequently, this learned function is applied to the noisy result measured for the actual target circuit, for which an exact value may be hard to compute. Specifically, given the noisy measurement \({X}_{\psi }^{\,\mathrm{noisy}}\) for a target state \(\left|\psi \right\rangle\) (prepared by circuit U_ψ) and observable X, Clifford data regression aims to find f^CF to approximate the ideal value \({X}_{\psi }^{\,\mathrm{exact}}=\langle \psi | X| \psi \rangle\) via the following relation:

where \(\vec{a}\) represents a set of parameters characterizing the function f^CF.

In the VQE experiment, we implement Clifford fitting in a quantum parameterized circuit. First, according to the method in ref. ⁵⁰, we distinguish between observables that vary with parameters and those that remain constant. Constant observables cannot be analysed using the Clifford fitting method. Let M denote the number of parameterized gates. Then, by randomly selecting K gates from M, we independently and uniformly sample their K parameters from the range [−π, π), while the remaining M − K gates are converted into identity operations. This process generates m random circuits \({\{{S}_{i}^{(K)}\}}_{i=1}^{m}\). The number m should be sufficient (typically over 10) to ensure an adequate fitting performance of \(\vec{a}\) in the linear regression. The dataset consists of pairs of noisy and exact expectation values, \({{\mathcal{T}}}_{\psi }=\{({X}_{\psi }^{\,\mathrm{noisy}},{X}_{\psi }^{\mathrm{exact}})\}\), generated from a set of quantum circuits \({\{{S}_{i}^{(K)}\}}_{i=1}^{m}\). Thus, the parameters \(\vec{a}\) of the correction function f^CF can be fitted by the the linear regression function

The parameters a₁ and a₂ are typically determined by minimizing the least-squares cost function

There is theoretical justification for the linear model, as it can perfectly correct for certain noise models, such as global depolarizing noise.

Finally, after determining the optimal parameters \({\vec{a}}^{* }=({a}_{1}^{* },{a}_{2}^{* })\), the error-mitigated expectation value for the original target circuit is predicted using equation (13):

Symmetry verification

Molecular Hamiltonians in quantum simulations inherently possess symmetry constraints such as particle number conservation. Symmetry verification through subspace expansion provides an effective error mitigation strategy. Following the methodology in ref. ⁵⁶, we implement this approach for the [[4, 2, 2]] code. For a general [[n, k, d]] code, the code space is defined by the projection operator:

where l = n − k; \({P}_{i}=\frac{1}{2}(I+{u}_{i})\) projects onto the +1 eigenspace of stabilizer generators u_i ∈ G = 〈u₁, …, u_l〉; and S denotes the stabilizer group. For a code-space Hamiltonian H_c = ∑_iγ_iΓ_i, where γ_i denotes the coefficient of the ith term, and Γ_i represents the corresponding Pauli string operator. the symmetry-verified expectation value becomes the following:

where \(c={\rm{Tr}}[\widehat{P}\rho ]\) normalizes the density matrix ρ projected onto the symmetry-preserving subspace.