introduction

Recently, quantum computation in the noisy intermediate scale quantum (NISQ) regime has been actively developed for the possibility of reaching practical quantum advantages without quantum error correction1,2. In this direction, quantum advantages with random gate sampling algorithms in superconducting systems have been demonstrated3,4. The possibility of surpassing conventional classical computation in solving useful and practical problems has been seriously explored5. In particular, a quantum-classical hybrid method that combines the advantages of classical and quantum computation, such as the variational quantum algorithm (VQA), has been theoretically and experimentally explored to reach higher performance than with only classical computers6,7,8,9. One of the critical issues in NISQ computation is to obtain results with high precision from quantum devices without quantum error correction10,11.

Various error-mitigation methods without using additional quantum resources for quantum error corrections have been proposed and applied to improve the computation results of quantum devices12,13,14,15. Representative error-mitigation methods include zero-noise extrapolation12,13,14,15,16,17, probabilistic error cancellation (PEC)12,13,14,15,18,19,20,21, measurement error mitigation22,23,24, symmetry verification25,26, virtual distillation27,28, subspace expansion29,30,31, N-representability methods32,33,34 and learning-based methods35,36. In particular, PEC is considered the general and systematic method for mitigating errors, which are not dependent on details of the physical implementation as long as errors are properly characterized. PEC consists of characterizing errors of each quantum gate, decomposing the inverse of the errors into combinations of basis operations, and sampling the operations with the characterized weights. In this way, we can estimate the expectation values of the desired observables without errors. PEC has been theoretically analyzed in detail and shown to be optimal for dephasing dominant noise situations12,13,20,21. However, it has been experimentally tested with only two-qubit systems of superconducting circuits and trapped ions18,19 and a multi-qubit system of the superconducting circuit via fitting a sparse error model motivated by the processor topology37.

Systematic error-mitigation schemes like PEC would enhance the performance of quantum computation in the NISQ regime drastically. Many important quantum algorithms have been developed based on the fact that quantum dynamics can be efficiently simulated in quantum computation, which also inspires VQA designs6,7,8,9. In digital quantum simulation, the time evolution of a certain Hamiltonian is typically implemented by the Trotter–Suzuki expansion38. Besides the errors in the expansion, imperfect primitive quantum gates also degrade the fidelity of the quantum simulation and eventually lead to faulty results. Due to this reason, it is an essential and urgent task to design and benchmark error-mitigation techniques that are both systematic and efficient in the NISQ era. As an instance for benchmark, the Fermi–Hubbard model39 is particularly interesting. This model is one of the fundamental models for interacting electrons, and it has been extensively studied in classical manners, either analytically40,41,42 or numerically43. As experimental quantum technologies develop, understanding and insights for strongly correlated quantum systems and high-temperature superconductivity44 can be explored through quantum simulation of the Fermi–Hubbard model eventually beyond the classical computational capability. This model has been utilized in experimental demonstration of various quantum platforms45,46,47,48. Recently, a combination of non-PEC methods has been applied to enhance simulation performance and observe the dynamics of the Fermi–Hubbard model in a superconducting system49.

In this paper, we experimentally study the possibility and limitations of PEC by applying the method to a four-qubit system with trapped ions50,51. In particular, we focus on simulating the dynamics of the 1D Fermi–Hubbard model by using multiple Trotter steps. We find the PEC method with a tomographically reconstructed error model significantly improves the fidelity of the simulation with up to three ion qubits. However, with four qubits, the improvement of fidelity by the PEC method alone is much worse than those for two- and three-qubit systems. This is because the tomographic characterization of the noisy gates is not precise enough to capture the time correlation and the cross-talk errors. With the application of total-spin and particle-number conservation15,25,26 as well as positivity constraint, however, the enhanced PEC method provides recovery for the ideal dynamics of the 1D spin-Fermi–Hubbard model, where the spin and the charge behave differently.

Results

Fermi–Hubbard model and experimental realization with trapped ions

Here we consider the extended Fermi–Hubbard (eFH) model, which describes a system of fermions moving in a one-dimensional lattice, as illustrated in Fig. 1a, b for one-component and two-component fermions. The one-component fermions can be considered spinless ones. The Hamiltonian reads

$$\begin{array}{l}{\hat{H}}_{{{{\rm{eFH}}}}}\,=\,-J\mathop{\sum}\limits_{l,\lambda }\left({\hat{c}}_{l,\lambda }^{{\dagger} }{\hat{c}}_{l+1,\lambda }+{{{\rm{h.c.}}}}\right)\\ \qquad\qquad +\,V\mathop{\sum}\limits_{l,\lambda ,{\lambda }^{{\prime} }}{\hat{n}}_{l,\lambda }{\hat{n}}_{l+1,{\lambda }^{{\prime} }}+U\mathop{\sum}\limits_{l}{\hat{n}}_{l,\uparrow }{\hat{n}}_{l,\downarrow },\end{array}$$
(1)

with \({\hat{n}}_{l,\lambda }\equiv {\hat{c}}_{l,\lambda }^{{\dagger} }{\hat{c}}_{l,\lambda }\), where \({\hat{c}}_{l,\lambda }\) (\({\hat{c}}_{l,\lambda }^{{\dagger} }\)) annihilates (creates) a fermion with spin λ (λ = , ) on the l-th site. Here J is the nearest-neighbor tunneling strength, and U and V respectively quantify the strength of the on-site and the nearest-neighbor interaction.

Fig. 1: Schematics of digital quantum simulation of interacting fermion models.
figure 1

(a) Spinless and (b) spin-\(\frac{1}{2}\) fermions moving in a one-dimensional lattice. The system dynamics can be unified by the extended Fermi–Hubbard model in Eq. (1), with the tunneling strength denoted by J and the on-site and the nearest-neighbor interaction strengths by U and V, respectively. c Trapped-ion system with four ion-qubits. The ions form a 1-D crystal and each ion represents a qubit. The ion qubits can be manipulated by Raman laser beams, consisting of a global beam and individual beams, with the former covering the whole ion crystal and each of the latter focusing on a single ion. Mølmer-Sørensen gate \(Y{Y}_{\varphi }(m,n)\equiv \exp \left(i\varphi {\hat{\sigma }}_{m}^{y}\otimes {\hat{\sigma }}_{n}^{y}\right)\) can be implemented between any qubit pair (Qm, Qn) by modulating the phase of the individual beams. d Quantum circuit for simulating the dynamics of a three-site single-component extended Fermi–Hubbard model. With the Trotter-Suzuki expansion, the evolution is approximately implemented by M identical steps. e Quantum circuit for a two-site double-component Fermi–Hubbard model. f Quantum circuit compilation. While the two-qubit entangling gate provided by the native gate set of the trapped-ion system is YYφ, XXφ and ZZφ can be constructed by YYφ and appropriate single-qubit rotations, such as \(\sqrt{Z}\equiv \exp \left(-i\frac{\pi }{4}{\hat{\sigma }}_{z}\right)\) and \(\sqrt{X}\equiv \exp \left(-i\frac{\pi }{4}{\hat{\sigma }}_{x}\right)\). The single-qubit gates (yellow squares) in (d, e) are the \({\hat{\sigma }}_{z}\)-rotations, i.e. \({Z}_{2\varphi }\equiv \exp \left(-i\varphi {\hat{\sigma }}_{z}\right)\).

We can map the fermion system to a qubit system by using the Jordan–Wigner transformation,

$$\begin{array}{l}{\hat{c}}_{l,\uparrow }\,=\,\mathop{\prod }\limits_{n=1}^{l-1}{\hat{\sigma }}_{n}^{z}{\hat{\sigma }}_{l}^{-},\quad {\hat{n}}_{l,\uparrow }=\frac{1}{2}\left(1-{\hat{\sigma }}_{l}^{z}\right),\\ {\hat{c}}_{l,\downarrow }\,=\,\mathop{\prod }\limits_{n=1}^{L+l-1}{\hat{\sigma }}_{n}^{z}{\hat{\sigma }}_{L+l}^{-},\quad {\hat{n}}_{l,\downarrow }=\frac{1}{2}\left(1-{\hat{\sigma }}_{L+l}^{z}\right),\end{array}$$
(2)

with \({\hat{\sigma }}_{n}^{\pm }=\left({\hat{\sigma }}_{n}^{x}\pm i{\hat{\sigma }}_{n}^{y}\right)/2\) and \({\hat{\sigma }}_{n}^{x,y,z}\) being the Pauli operators on the n-th qubit. Thus generally speaking, we need N = 2L qubits to simulate the dynamics of an L-site fermionic chain.

As shown in Fig. 1c, we use 171Yb+ ions for qubit systems. The 171Yb+ ions are confined in a linear Paul trap. The hyperfine levels in the 2S1/2 ground-state manifold, i.e. \(\left\vert F=0,{m}_{F}=0\right\rangle\) and \(\left\vert F=1,{m}_{F}=0\right\rangle\), are encoded as the qubit \(\left\{\left\vert 0\right\rangle ,\left\vert 1\right\rangle \right\}\), which can be initialized to \(\left\vert 0\right\rangle\) by optical pumping and distinctively detected by fluorescence measurement via multi-channel photo-multiplier tube52.

The qubit-system Hamiltonian becomes

$$\begin{array}{ll}{\hat{H}}_{{{{\rm{qb}}}}}\,=\,{\hat{H}}_{X}+{\hat{H}}_{Y}+{\hat{H}}_{Z},\\ {\hat{H}}_{X}\,=\,\frac{J}{2}\mathop{\sum }\limits_{l=1}^{L-1}\left({\hat{\sigma }}_{l}^{x}{\hat{\sigma }}_{l+1}^{x}+{\hat{\sigma }}_{L+l}^{x}{\hat{\sigma }}_{L+l+1}^{x}\right)\\ {\hat{H}}_{Y}\,=\,\frac{J}{2}\mathop{\sum }\limits_{l=1}^{L-1}\left({\hat{\sigma }}_{l}^{y}{\hat{\sigma }}_{l+1}^{y}+{\hat{\sigma }}_{L+l}^{y}{\hat{\sigma }}_{L+l+1}^{y}\right)\\ {\hat{H}}_{Z}\,=\,\frac{V}{4}\mathop{\sum }\limits_{l=1}^{L-1}\left[\left(1-{\hat{\sigma }}_{l}^{z}\right)\left(1-{\hat{\sigma }}_{l+1}^{z}\right)\right.\\ \qquad\quad \left.+\,\left(1-{\hat{\sigma }}_{L+l}^{z}\right)\left(1-{\hat{\sigma }}_{L+l+1}^{z}\right)\right]\\ \qquad\quad +\,\frac{U}{4}\mathop{\sum }\limits_{l=1}^{L}\left(1-{\hat{\sigma }}_{l}^{z}\right)\left(1-{\hat{\sigma }}_{L+l}^{z}\right).\end{array}$$
(3)

Note that according to the orientation of the Pauli operators, we divide the qubit-system Hamiltonian \({\hat{H}}_{{{{\rm{qb}}}}}\) into three parts, each of which contains mutually commutative terms.

To extract dynamical properties, the essential part is to simulate the evolution operator \(\hat{U}(t)=\exp \left(-i{\hat{H}}_{{{{\rm{qb}}}}}t\right)\) with a universal gate set available in a certain quantum computational platform. One of the well-known solutions is to use the Trotter-Suzuki decomposition. The first-order Trotter-Suzuki decomposition of \(\hat{U}(t)\) reads

$$\begin{array}{r}\hat{U}(t)={\left({e}^{-i{\hat{H}}_{Z}\delta t}{e}^{-i{\hat{H}}_{Y}\delta t}{e}^{-i{\hat{H}}_{X}\delta t}\right)}^{M}+{{{\mathcal{O}}}}\left(M\delta {t}^{2}\right),\end{array}$$
(4)

with the Trotter step size being δt = t/M. Figure 1d and e show examples of the circuits for simulating the dynamics of 3-site spinless and 2-site spin-1/2 fermionic systems, with the native gate set including the Mølmer-Sorensen gates, i.e. XXφ, YYφ and ZZφ, and single-qubit rotations. By applying the amplitude-shaped53 global and individual laser beams on the desired ions and driving the transverse motional modes54,55, we can realize the native Mølmer-Sørensen YY-gate (YYφ). The gate XXφ and ZZφ are composed with YYφ and corresponding single qubit rotations, as shown in Fig. 1f.

In this experiment, we consider two scenarios, one for spinless fermions, and the other for two-component fermions. In the former case, the on-site interaction strength naturally vanishes, i.e. U = 0, due to the fermionic anti-commutation relation. In this case, we use N = L ion-qubits to simulate the dynamics of a spinless fermionic chain with L sites. While for the latter case, we only consider the on-site interaction, since the interaction strength decreases as the distance increases and thus UV in most circumstances.

Probabilistic error cancellation in trapped ion system

Here, we introduce the PEC error-mitigation method, which systematically recovers ideal expectation values of target quantum circuits with the experimental noisy gates14. The essential idea is to first characterize experimental noisy gates, then identify the noise parts by comparing the ideal and the noisy gates, and then decompose the inverse of the noise parts into combinations of basis operations. Then, each of the ideal gates in target quantum circuits is replaced by a corresponding noisy one followed by basis operations sampled with probability distributions obtained from the inverse-noise decomposition. At last, the ideal expectation value can be obtained by adding up those of the sampled circuits with appropriate weights.

In circuit quantum computation, a computational task can be compiled into a quantum circuit consisting of single- and two-qubit gates. In our system, the error of the two-qubit gate is an order larger than that of the single-qubit gate19, and we mitigate only the errors in two-qubit gates.

To characterize the experimental noisy gates, we perform quantum process tomography (QPT) under the Pauli error assumption. Different from gate-set tomography, quantum process tomography does not account for state preparation and measurement (SPAM) errors. Here, we assume the initial state preparation is perfect and apply the detection error correction method for the measurement results56, which enables us to exclude SPAM errors. From the quantum process tomography, we obtain the Pauli transfer matrices (PTMs) of the Mølmer-Sørensen YY gates on different qubit pairs, i.e. \({R}_{Y{Y}_{\varphi }^{{{{\rm{ns}}}}}(m,n)}\) for the qubit pair (Qm, Qn). Together with the ideal PTM \({R}_{Y{Y}_{\varphi }}\), we can formally define an error operator for each of the noisy gates, with the PTM being

$$\begin{array}{r}{R}_{{{{{\mathcal{E}}}}}_{\varphi }(m,n)}={R}_{Y{Y}_{\varphi }^{{{{\rm{ns}}}}}(m,n)}{R}_{Y{Y}_{\varphi }}^{-1}.\end{array}$$
(5)

Note that under the Pauli-error assumption, both \({R}_{{{{\mathcal{E}}}}}\) and its inverse \({R}_{{{{\mathcal{E}}}}}^{-1}\) can be decomposed into linear superpositions of PTMs of two-qubit Pauli operators Pi {II, IX, IY, . . . , ZZ}. As a result, ideal PTMs can be rewritten in terms of the PTMs of the noisy gate and the corresponding error operators,

$$\begin{array}{r}{R}_{Y{Y}_{\varphi }}={R}_{{{{{\mathcal{E}}}}}_{\varphi }(m,n)}^{-1}{R}_{Y{Y}_{\varphi }^{{{{\rm{ns}}}}}(m,n)}\end{array}$$
(6)

with the inverse-error decomposition

$$\begin{array}{r}{R}_{{{{{\mathcal{E}}}}}_{\varphi }(m,n)}^{-1}=\mathop{\sum }\limits_{i=0}^{15}{q}_{i}(m,n){R}_{{P}_{i}}.\end{array}$$
(7)

Here the quasi-probabilities qi(m, n) are real and satisfy ∑iqi(m, n) = 1. The fact that qi(m, n) can be negative indicates that the inverse of the error is not physical and thus can not be implemented by a deterministic quantum operation. To cancel the effect of the error operator, we rewrite the decomposition as follows

$$\begin{array}{r}{R}_{{{{{\mathcal{E}}}}}_{\phi }(m,n)}^{-1}={C}_{Y{Y}_{\varphi }(m,n)}\mathop{\sum }\limits_{i=0}^{15}{p}_{i}(m,n){{{\rm{sgn}}}}\left[{q}_{i}(m,n)\right]{R}_{{P}_{i}},\end{array}$$
(8)

where \({p}_{i}(m,n)={C}_{Y{Y}_{\varphi }(m,n)}^{-1}\left\vert {q}_{i}(m,n)\right\vert\) are well-defined probabilities with the cost \({C}_{Y{Y}_{\varphi }(m,n)}=\mathop{\sum }\nolimits_{i = 0}^{15}\left\vert {q}_{i}(m,n)\right\vert\), and \({{{\rm{sgn}}}}(\cdot )\) is the sign function.

Then, we can use the Monte-Carlo sampling to compute the PEC error-mitigated results. Here we consider a target quantum circuit consisting of Ng number of two-qubit gates YYφ(m, n), with the g-th gate labeled as Gg, each of which is fully characterized and the inverse-error decomposition is written as \({R}_{{G}_{g}}={\sum }_{i}{q}_{i,g}{R}_{{P}_{i}}{R}_{{G}_{g}^{{{{\rm{ns}}}}}}\) with qi,g being the quasi-probability of the two-qubit Pauli operator Pi. The full set of the random circuits to be implemented in the experimental device denoted as \(\left\{{{{{\mathcal{S}}}}}_{{{{\bf{i}}}}}:{{{\bf{i}}}}=\left({i}_{1},\ldots ,{i}_{{N}_{g}}\right)\right\}\), contains \(1{6}^{{N}_{g}}\) different circuits, each of which is obtained by adding \({P}_{{i}_{g}}\) right after Gg. To obtain the ideal expectation values, each random circuit \({{{{\mathcal{S}}}}}_{{{{\bf{i}}}}}\) is assigned with a probability \(\propto {\prod }_{g}\vert {q}_{{i}_{g},g}\vert\) and a sign \({{{{\rm{sgn}}}}}_{{{{\bf{i}}}}}\equiv {\prod }_{g}{{{\rm{sgn}}}}({q}_{{i}_{g},g})\). As it is infeasible to implement all these circuits, we use the Monte Carlo sampling to obtain average values over the random circuits. Specifically, we generate Ns random circuits by sampling \({P}_{{i}_{g}}\) from the probability distribution \(\propto \vert {q}_{{i}_{g},g}\vert\) and obtain the expectation value 〈μi of an observable μ for the circuit \({{{{\mathcal{S}}}}}_{{{{\bf{i}}}}}\) by averaging 300 repetitions. The ideal expectation value 〈μ〉 is then calculated by

$$\begin{array}{r}\langle \mu \rangle =\frac{C}{{N}_{s}}\mathop{\sum }\limits_{s=1}^{{N}_{s}}{{{{\rm{sgn}}}}}_{{{{{\bf{i}}}}}_{s}}{\langle \mu \rangle }_{{{{{\bf{i}}}}}_{s}},\end{array}$$
(9)

where the cost is given by \(C\equiv \mathop{\prod }\nolimits_{g = 1}^{{N}_{g}}{C}_{g}\) with the cost of the g-th gate being \({C}_{g}\equiv \mathop{\sum }\nolimits_{{i} = 0}^{15}\vert {q}_{{i},g}\vert\). Considering each expectation value 〈μi has the standard deviation \(\Delta\mu_{\mathbf i}\), the error bar of the ideal expectation value 〈μ〉 is obtained by error propagation as

$$\begin{array}{r}\Delta\mu=\frac{C}{N_s}\sqrt{{\sum \nolimits_{s = 1}^{{N}_{s}}}\Delta^{2}\mu_{{\mathbf{i}}_s}},\end{array}$$
(10)

which is approximately C times larger than the error bar of one general circuit with Ns × 300 repetitions12,13,14.

Error mitigation and simulation of two spinless fermionic modes

We now discuss the simulation of two spinless fermionic modes, where the qubit-system Hamiltonian Eq. (3) reduces to

$$\begin{array}{r}H=\frac{J}{2}\left({\hat{\sigma }}_{1}^{x}{\hat{\sigma }}_{2}^{x}+{\hat{\sigma }}_{1}^{y}{\hat{\sigma }}_{2}^{y}\right)+\frac{V}{4}\left({\hat{\sigma }}_{1}^{z}{\hat{\sigma }}_{2}^{z}+{\hat{\sigma }}_{1}^{z}+{\hat{\sigma }}_{2}^{z}\right).\end{array}$$
(11)

With the Trotter-Suzuki expansion of Eq. (4), the evolution of Hamiltonian for two spinless fermionic modes can be experimentally realized with the specific quantum circuit shown in Fig. 2a. For the singly occupied fermion in one of the sites, the fermion oscillates between two sites caused by the nearest-neighbor tunneling. For the doubly occupied fermions, there will be no dynamics due to the Pauli exclusion principle of fermions. We apply M = 8 Trotter steps and each step consists of three YYφ gates and ten single-qubit rotations. The whole quantum circuit consists of twenty-four two-qubit gates and eighty-one single-qubit gates, where the additional one is for the initial state preparation. We set the nearest-neighbor tunneling and interaction strength as V = 2J.

Fig. 2: Error-mitigated quantum simulation.
figure 2

a Quantum circuit of Trotter expansion for the simulation of dynamics of two spinless fermions. Two qubits are initialized to \((\left\vert 11\right\rangle +\left\vert 10\right\rangle )/\sqrt{2}\) after the state preparation, then we apply M times of Trotter steps consisting of three two-qubit gates YYns in each step, which is not perfect. Here, J and V represent the strength of nearest-neighbor tunneling and nearest-neighbor interaction. b Probabilistic error cancellation (PEC) scheme. The imperfections of the gate \({YY}_{\varphi }^{{{{\rm{ns}}}}}\) can be described by error channels acting on the ideal gate. The ideal gate YYφ can be realized by including the inverse of the error channels, which can be decomposed into combinations of basis operations with quasi-probabilities. Under the Pauli error assumption, there are 16 basis operations, where the list is shown in (ce). c Pauli transfer matrix and (d) Deviation of the experimental \({YY}_{\varphi }^{{{{\rm{ns}}}}}\) gate, which is characterized through QPT. e Quasi-probabilities in the decomposition of the inverse error operations of experimental \({YY}_{\varphi }^{{{{\rm{ns}}}}}\) gate. The error bars are around 0.002 on average. f Error-mitigated quantum circuits by the PEC. We implement the inverse error operations by sampling with the renormalized probabilities of quasi-probabilities as \({p}_{i}=| {q}_{i}| /\mathop{\sum }\nolimits_{i = 0}^{15}| {q}_{i}|\), where i indicates what number it corresponds to among basis operations. g Experimental data and (h) error-mitigated data for V = 2J. The solid points represent the state populations measured from the experiment and the error-mitigated data after applying PEC, MLE, and PS methods. The dash–dotted lines represent numerical simulation, where no Trotter errors exist. The colors of the lines correspond to the experimental data points. i Population fidelity of experiment data and error-mitigated data. The fidelity is calculated by \(| {\sum }_{k}\sqrt{{P}_{k,{{{\rm{ideal}}}}}{P}_{k}}{| }^{2}\), which k represents the state basis. The population fidelities after each error-mitigation method are presented. All shadows in (g)–(i) represent the standard deviation of 1000 samples generated from the raw experimental data, using the bootstrapping method.

Errors caused by gate imperfection accumulate as the number of Trotter steps increases and the fidelity of evolved dynamics is degraded. We mitigate the errors and improve the results of dynamics. Since it is a two-qubit system, only the one native YYπ/4 gate needs to be mitigated. The ideal YYπ/4 gate can be decomposed as \({R}_{{YY}_{\pi /4}}=\mathop{\sum }\nolimits_{i = 0}^{15}{q}_{i}{R}_{{P}_{i}}{R}_{Y{Y}_{\pi /4}^{{{{\rm{ns}}}}}}\) under Pauli error assumption, which is visualized as Fig. 2b. The experimental \(Y{Y}_{\pi /4}^{{{{\rm{ns}}}}}\) gate is characterized by QPT. The experimental results of the PTMs of the \(Y{Y}_{\pi /4}^{{{{\rm{ns}}}}}\) gate and the deviation from the ideal one obtained by \({R}_{Y{Y}_{\pi /4}^{{{{\rm{ns}}}}}}-{R}_{{YY}_{\pi /4}}\) are shown in Fig. 2c and d, where the average gate fidelity of \({R}_{Y{Y}_{\pi /4}^{{{{\rm{ns}}}}}}\) is 0.9811 ± 0.004757,58. The quasi-probabilities, i.e. the coefficients qi of the two-qubit Pauli operator Pi in the decomposition of \({R}_{{{{\mathcal{E}}}}}^{-1}\), are shown in Fig. 2e. With the values of the quasi-probabilities qi, we obtain the cost \({C}_{Y{Y}_{\pi /4}}=1.083\) for \({R}_{Y{Y}_{\pi /4}^{{{{\rm{ns}}}}}}\).

Then, we generate Ns = 1000 random circuits by Monte Carlo sampling illustrated in Fig. 2f. With the method described in Eq. (9), we obtain the dynamical evolutions of state populations \({P}_{\left\vert ab\right\rangle }\) for two spinless fermionic modes, where the observables are \(\mu =\left\vert ab\right\rangle \left\langle ab\right\vert\) with occupied fermion a(b) = 1 or no occupied fermion a(b) = 0 for each mode. The error-mitigated state populations after PEC can be negative values, which are unphysical. We address this problem by imposing constraints based on general properties of probabilities, i.e. \(0\le {P}_{\left\vert ab\right\rangle }\le 1\) and \({\sum }_{a,b}{P}_{\left\vert ab\right\rangle }=1\), and obtain physical state populations by the maximum-likelihood estimation (MLE). Specifically, we minimize the norm distance between the physical populations and the PEC error-mitigated ones. On top of PEC and MLE, we also perform post-selection (PS) with respect to the symmetry constraints of the model Hamiltonian, i.e. the conservation law of the number of fermions. To demonstrate the effect of each of the error-mitigation techniques, i.e. PEC, MLE, and PS, we apply them step by step and obtain a set of error-mitigated state populations after each step.

The experimental results of fermionic dynamics are shown in Fig. 2g and h for the cases without and with error mitigation, respectively. We use the initial two-qubits state of \((\left\vert 11\right\rangle +\left\vert 10\right\rangle )/\sqrt{2}\), which contains the superposition of two occupied fermions and a single occupied fermion. This initial state simultaneously reveals the time evolutions of both single and double fermions, since no interference occurs between both dynamics. As shown in both Fig. 2g and h, the state \(\left\vert 11\right\rangle\) and \(\left\vert 00\right\rangle\) do not change and there is an oscillation between the states \(\left\vert 01\right\rangle\) and \(\left\vert 10\right\rangle\), which is the hallmark of fermion dynamics46. We can see that without error mitigation, the contrast between the experimental populations (solid points) clearly decays compared with the numerical simulation under the Trotter step (hollow points) as shown in Fig. 2g, and the error-mitigated results are closer to the ideal numerical simulations as shown in Fig. 2h. The population fidelities depending on Trotter steps without and with the three error-mitigation schemes are shown in Fig. 2i. The fidelity of the gate without and with error mitigation are 0.9895 ± 0.0021 and 0.9975 ± 0.0009, which are obtained by exponential fitting of data.

We note that the main improvement in fidelities comes from the PEC method. It is because our characterization of the gate is a good representation of the actual gate. We confirm it by simulating the dynamics of two spinless modes with the experimentally constructed PTM of \({R}_{Y{Y}_{\pi /4}^{{{{\rm{ns}}}}}}\). The comparisons between the simulation with the noisy PTM and experimental data clearly validate the Pauli error assumption and accuracy of measured \({R}_{Y{Y}_{\pi /4}^{{{{\rm{ns}}}}}}\) (see Methods). We also note that MLE reduces the error bar of the results after PEC, which may compensate for the drawback of PEC in the aspect of variance. All the error bars are estimated by bootstrapping method (see Methods).

Simulation of three spinless fermionic modes

We implement the time evolution of three spinless fermionic modes as shown in Fig. 1a. The quantum circuit in Fig. 1d, obtained by the Trotter-Suzuki expansion, simulates the dynamics of the system. We initialize the three-qubit state to \((\left\vert 101\right\rangle +\left\vert 110\right\rangle )/\sqrt{2}\), which includes two fermions in three sites. The fermions propagate in the system caused by the nearest-neighbor tunneling to other sites. The dynamical evolution would be restricted in the subspace spanned by \(\{\left\vert 110\right\rangle ,\left\vert 101\right\rangle ,\left\vert 011\right\rangle \}\) because of the conservation law of the number of fermions.

The dynamics can be different depending on the existence of nearest-neighbor interaction different from the previous case with two sites. In the experiment, we choose V = 0 and V = 2J for the scenarios without and with nearest-neighbor interactions, respectively. In Fig. 3a, b, d, and e, dash–dotted lines show the results of ideal numerical simulation. When V = 0, the initial state \(\left\vert 110\right\rangle\) mostly evolves to the state \(\left\vert 101\right\rangle\) because of the Pauli exclusion principle of fermions. On the other hand, the initial state \(\left\vert 101\right\rangle\) can evolve to both \(\left\vert 110\right\rangle\) and \(\left\vert 011\right\rangle\) states. Therefore, in the beginning, the overall population of the state \(\left\vert 110\right\rangle\) reduces, and the population of the state \(\left\vert 011\right\rangle\) increases as shown in Fig. 3a and b. When V = 2J, we can see the population of \(\left\vert 110\right\rangle\) is not rapidly decreasing and the state \(\left\vert 101\right\rangle\) is not populated compared to the case of V = 0. It is because the interaction is effectively attractive in our experiment.

Fig. 3: Dynamics of three spinless fermionic modes.
figure 3

Output state population with the initial state of \((\left\vert 101\right\rangle +\left\vert 110\right\rangle )/\sqrt{2}\) ac for no interaction V = 0 and df for V = 2J. a, d Experimentally measured populations without error mitigation and (b), (e) with error mitigation. The solid points represent the experimentally measured state populations and the dashed–dotted lines represent the numerical simulation, which contains the average Trotter errors of 3.2% in (a) and (b) and 10.4% in (d) and (e). c, f The population fidelities after each error-mitigation method as PEC, MLE, and PS. All shadows in (a)–(f) represent the standard deviation of 1000 samples generated from the raw experimental data, using the bootstrapping method.

In experiment, we set \(\varphi =\frac{J\delta t}{2}=\frac{\pi }{8}\) of the YYφ gates and apply M = 4 Trotter steps for both scenarios. With nearest-neighbor interactions, each Trotter step consists of six YYφ gates and twenty single-qubit rotations. The whole quantum circuit consists of twenty-five two-qubit gates and eighty-one single-qubit gates, where the additional single- and two-qubit gates are for the initial state preparation. Then, we apply the error-mitigation methods to the dynamic evolution of both scenarios. There are two types of YY gates on different pairs of qubits, i.e., YYπ/8(1, 2) and YYπ/8(2, 3), both of which require characterization and inverse-error decomposition. In the PEC scheme, both YYπ/8(1, 2) and YYπ/8(2, 3) gates are individually characterized by QPT and the inverse errors are decomposed independently. The average gate fidelities are 0.9779 ± 0.0088 for \({R}_{Y{Y}_{\pi /8}^{{{{\rm{ns}}}}}(1,2)}\) and 0.9748 ± 0.0085 for \({R}_{Y{Y}_{\pi /8}^{{{{\rm{ns}}}}}(2,3)}\) (see Methods).

For the case of V = 0, the experimental results without and with error mitigation are shown in Fig. 3a and b, respectively. We can see the error-mitigated results are more consistent with ideal numerical simulations. The fitted population fidelities per gate without and with error mitigation are 0.9792 ± 0.0161 and 0.9942 ± 0.0067, as shown in Fig. 3c. We note that the main improvement comes from the PEC method. This is consistent with the comparison results between the simulation with the noisy PTMs of both YY gates and experimental data, which indicates the measured PTMs properly characterize the two-qubit gate errors (see Methods). Similar conclusions are obtained for the scenario with the nearest-neighbor interaction of V = 2J, as shown in Fig. 3d–f. The fitted population fidelities per gate without and with error mitigation are 0.9827 ± 0.0142 and 0.9889 ± 0.0080. The two-qubit gate errors in the three-qubit system are larger than that in the two-qubit system, which leads to the larger cost as \({C}_{Y{Y}_{\pi /8}(1,2)}=1.157\) and \({C}_{Y{Y}_{\pi /8}(2,3)}=1.171\). Therefore, we sample more random circuits as Ns = 1500 to reduce the error bars of the final error-mitigated results.

We note that the fitted population fidelities after PEC could be unphysical because of the possible negative quasi-probabilities and the cost being over 1. As shown in Fig. 3c, f, fidelities greater than 1 after applying PEC exist at the beginning of the three-fermion experiments. With the help of the MLE method based on the general properties of probabilities, we can make the population fidelities physical and reasonable, which is, below 1 as shown in green dots of Fig. 3c, f. In addition, MLE reduces the error bars of the results after PEC, which may compensate for the drawback of PEC in the aspect of variance. After applying PEC, the standard deviation of the fourth Trotter step become unexpectedly large. However, the MLE method suppresses the standard deviation by a few factors in this case.

Simulation of four fermionic modes

We implement the simulation of four fermionic modes with spins by encoding qubits Q1Q3 and Q2Q4 as two fermionic sites, where each site can contain two fermionic modes with different spins, as shown in Fig. 1b. We initialize the four qubits state to \((\left\vert 1001\right\rangle +\left\vert 1010\right\rangle )/\sqrt{2}\), which is corresponding to the superposition state between \({\left\vert \uparrow \right\rangle }_{1}{\left\vert \downarrow \right\rangle }_{2}\) and \({\left\vert \uparrow \right\rangle }_{1}{\left\vert \downarrow \right\rangle }_{1}\) as shown in Fig. 1b. The two fermions interact and exchange on the four modes caused by the nearest-neighbor tunneling and on-site interaction on subspace states \(\{\left\vert 0101\right\rangle ,\left\vert 0110\right\rangle ,\left\vert 1001\right\rangle ,\left\vert 1010\right\rangle \}\) with conserved fermions and spins. The corresponding quantum circuit is illustrated in Fig. 1e. Here, we expect the different dynamic behaviors of spin and charge depending on the on-site interaction strength. Spin and charge of the fermions are defined as spin = nj, − nj, and charge = nj, + nj, with nj,() being the fermion number of spin-up (down) on the j-th site.

To demonstrate the different behavior of spin and charge, we consider two scenarios U = 0 and U = 2J, which are related without and with on-site interactions, respectively. For the case of no on-site interaction, the spin and the charge of fermions should show the same dynamics. But if there is on-site interaction, the dynamics for spin and charge should show a different behavior, which is related to spin-charge separation for low energy excitation49,59,60,61. For the case of no on-site interaction, U = 0, the dynamics of fermions with spins are similar to those without spins, which are discussed in Fig. 2g and h. The state \({\left\vert \uparrow \right\rangle }_{1}{\left\vert \downarrow \right\rangle }_{2}\)(\(\left\vert 1001\right\rangle\)) is going back and forth with the other state \({\left\vert \uparrow \right\rangle }_{2}{\left\vert \downarrow \right\rangle }_{1}\)(\(\left\vert 0110\right\rangle\)). Likewise, the state \({\left\vert \uparrow \right\rangle }_{1}{\left\vert \downarrow \right\rangle }_{1}\)(\(\vert 1010\rangle\)) is oscillating with the other state \({\left\vert \uparrow \right\rangle }_{2}{\left\vert \downarrow \right\rangle }_{2}\)(\(\left\vert 0101\right\rangle\)). This leads to equal amounts of charge and spin change at each site. However, when there exists on-site interaction, the dynamics of those fermionic states are no longer simple oscillations, and the state \({\left\vert \uparrow \right\rangle }_{1}{\left\vert \downarrow \right\rangle }_{2}\), which is not affected by the on-site interaction, and the state \({\left\vert \uparrow \right\rangle }_{1}{\left\vert \downarrow \right\rangle }_{1}\), which is strongly influenced by the on-site interaction, show completely different dynamics.

In experiment, we set \(\varphi =\frac{J\delta t}{2}=\frac{\pi }{8}\) and apply M = 4 Trotter steps for both scenarios. With on-site interactions, each step consists of six YYφ gates and twenty single-qubit rotations. The whole quantum circuit consists of twenty-five two-qubit gates and eighty-two single-qubit gates, where the additional single- and two-qubit gates are for the initial state preparation. In order to apply error mitigations, in particular, the PEC scheme, we characterize four types of YY gates on different pairs of qubits, i.e., YYπ/8(1, 2), YYπ/8(3, 4), YYπ/8(1, 3), and YYπ/8(2, 4), individually and decompose their inverse errors independently. The average gate fidelities are 0.9755 ± 0.0096 for \({R}_{Y{Y}_{\pi /8}^{{{{\rm{ns}}}}}(1,2)},0.9706\pm 0.0085\) for \({R}_{Y{Y}_{\pi /8}^{{{{\rm{ns}}}}}(3,4)},0.9720\pm 0.0092\) for \({R}_{Y{Y}_{\pi /8}^{{{{\rm{ns}}}}}(1,3)}\) and 0.9744 ± 0.0091 for \({R}_{Y{Y}_{\pi /8}^{{{{\rm{ns}}}}}(2,4)}\) (see Methods).

For the case of U = 0, the experimental results without and with error mitigation are shown in Fig. 4a and b, respectively. We can see the error-mitigated results are more consistent with the numerical simulation. The fitted population fidelities per gate without and with error mitigation are 0.9658 ± 0.0426 and 0.9993 ± 0.0053, respectively, as shown in Fig. 4c. In Fig. 4d, the same behaviors of spin and charge, where the amounts of charge and spin changes are the same in the dynamics, are clearly shown with error mitigation, where the net charge and the net spin at site 1 are displayed. We can observe similar improvements in the dynamics with on-site (U = 2J) interaction, as shown in Fig. 4e, and f. In Fig. 4g, the fitted population fidelities per gate without and with error mitigation are 0.9551 ± 0.0285 and 0.9925 ± 0.0081, respectively. With on-site interaction, the different behaviors of spin and charge are more clearly shown after error mitigations as shown in Fig. 4h.

Fig. 4: Dynamics of four fermionic modes mapped by two fermion sites and spins.
figure 4

Ouput population with the initial state of \((\left\vert 1001\right\rangle +\left\vert 1010\right\rangle )/\sqrt{2}\) ad for no interaction U = 0 and eh for U = 2J. a, e Experimentally measured population without error mitigations and (b, f) with error mitigation. The solid points represent the experimentally measured state populations and the dashed–dotted lines represent the numerical simulation, where no Trotter errors in (a) and (b) and the average Trotter errors of 10.4% in (e) and (f). c, g The population fidelities after each error-mitigation method as PEC, MLE, and PS. d, h Fermionic spin and charge population. The hollow and filled symbols represent the experimental data without and with error mitigations, respectively. The dashed–dotted lines represent the ideal numerical simulation with Trotter errors. All shadows represent the standard deviation of 1000 samples generated from the raw experimental data, using the bootstrapping method.

We note that the PEC scheme alone does not recover the dynamics for the four fermionic modes as shown in Fig. 4c and g different from the simulations with two and three fermionic modes. The main improvements in fidelities come from the post-selection methods based on the fermion-number-conservation assumption. One reason is that the experimental PTMs of two-qubit gates in the four-qubit system do not properly characterize the imperfections of those gates. It can be seen from the inconsistency of the numerically simulated data with the PTMs of two-qubit gates and experimental results (see Methods). The PTMs in our experiments cannot accommodate infidelities from cross-talks between qubits and time-correlated errors, that is, the performance of the gates during characterization can be different from those of actual gates in the quantum simulation. We note that coherent errors beyond the Pauli-error model assumption also are not captured by our experimental PTMs, which partially characterize errors related to Pauli errors. However, errors from cross-talks and time-correlations would be more important factors because coherent errors have no obvious reason to be seriously worse with system size. We also note that the number of sampling is greatly increased for applying the PEC method in the four-qubit system. With the measured fidelities, the costs increase to \({C}_{Y{Y}_{\pi /8}(1,2)}=1.211,{C}_{Y{Y}_{\pi /8}(3,4)}=1.228,{C}_{Y{Y}_{\pi /8}(1,3)}=1.228\) and \({C}_{Y{Y}_{\pi /8}(2,4)}=1.223\). The experimental Monte-Carlo sampling as Ns = 2000 is not sufficient, which is limited by the time-consuming running of massive random circuits in physical hardware. After applying PEC, the standard deviation of the fourth Trotter step becomes unexpectedly large, which is related to the insufficient samples as the Trotter step increases.

Discussion

We have applied the PEC method to improve the simulation of interacting fermions on up to four trapped ion qubits. We have observed the improvement of the simulation close to an order of magnitude in fidelities. With four qubits, we are able to observe the different dynamics of spin and charges for the interacting fermions.

Although PEC is the fundamental method in the experiment, several other error migration methods have been developed and applied. One of them is using the positive probability constraint. One intrinsic problem in the PEC method is that the probabilities of quantum states can be negative and not necessarily normalized. We impose constraints such that the probabilities of fermion states in the dynamics are positive and normalized, and determine the probabilities through the maximum-likelihood method. This method guarantees that observables are physically reasonable and can suppress statistical errors.

We also apply the symmetry constraint such that the total number of fermions and spins are conserved. In the simulation with four qubits, the majority of infidelity comes from the populations outside the number-conserving states. The constraints of total fermions conservation resolve the problem of population leakages. However, it is questionable if the number-conservation constraints would be scalable since the leakage can increase with the system and circuit depth.

From our experimental demonstration, although PEC incorporating constraints recoveries the dynamics of four trapped-ion qubits, it is clear that we need further improvements on the method for large-scale implementation. In particular, time and spatial correlations cause remaining errors after PEC as pointed out in Ref. 35. Here, time correlation refers to the difference between gates when characterizing and when implementing. This can be caused by changes in laser parameters such as amplitude, phase, and frequency over time, which determine the performance of the gate. The spatial correlation is basically the crosstalk of the gates to the neighboring ions, which can come from the spillover of focused laser beams. Recently, several proposals and experiments have been reported to tackle these correlations such as using sparse Pauli-Lindblad models37, learning-based error-mitigation35, matrix product operator representation62, etc. With all these techniques combined together, NISQ computers can potentially achieve the accuracy required for gaining the quantum advantage in practical applications, such as model inference for nuclear magnetic resonance (NMR) spectroscopy63 and large-depth QAOA64.

Methods

Estimating error bars with the bootstrapping method

We statistically estimated the uncertainty of experimental data with and without error mitigations by using the bootstrapping method. In order to estimate the variance of experimental data, we repeatedly and randomly take out samples from our original data and apply PEC, MLE, and PS methods step by step to mitigate the population errors of each sample. In our work, the size of the samples is the same as the original data sets, and the total number of samples for each data point is usually chosen as 1000. Then we calculate the standard deviation of the state population and the population fidelities after each step with the sets of samples. In order to avoid the appearance of unphysical fidelities (larger than 100% or less than 0% after MLE and PS methods), we separated the sets of samples into two parts according to the fidelities of the original data set and calculated the standard deviation of each side, which resulted in unsymmetrical error bars. Figure 5 shows the fidelity distribution (step 4 for the simulation of three spinless fermions, V=0) after each error-mitigation step, which indicates an increasing deviation from the normal distribution after more error-mitigation steps are applied.

Fig. 5: Fidelity distribution at the fourth Trotter step in the simulation of three spinless fermions.
figure 5

All the samples are generated from the bootstrap method after a PEC method, b MLE method, and c PS method. The black lines in all the figures represent the original calculated results. Here each figure contains 1000 samples.

The PTMs and quasi-probabilities for three and four qubit system

The quantum circuits for the simulation of the three spinless fermionic modes (three-site single-component) and four fermionic modes (two-site double-component Fermi–Hubbard model) are shown in Fig. 1d, e. Here up to six two-qubit Mølmer-Sørensen gates YYφ are applied in each Trotter step. We classify the gates depending on related qubit pairs and benchmark each type of gate with the QPT method. Under Pauli error assumption, we measure 15 different components in \({R}_{YY^{\rm ns}_{\varphi }}\) to get the PTMs of each type of gate as shown in Fig. 6a–f. The quasi-probabilities for YYφ gates in the simulation of three spinless fermions and two fermions with spins are shown in Fig. 6g–l, which decides the probabilities of choosing Pauli operations in the Monte Carlo method.

Fig. 6: PTMs and quasi-probabilities measured in three and four qubits system.
figure 6

a, b PTMs in the simulation of three spinless fermions. cf PTMs in the simulation of two fermions with spins. g, h Quasi-probabilities in the simulation of three spinless fermions. il quasi-probabilities in the simulation of two fermions with spins.

The PTM simulation to understand experimental error models

To characterize our noisy two-qubit gates, a partial QPT with Pauli error assumption is applied with 15 measurements out of 144. To verify the assumption that the errors of our gates are consistent with the increasing of Trotter steps, we perform the numerical simulation of the extended Fermi–Hubbard model with the PTMs of the noisy gates. The simulation results are shown in Fig. 7 for two spinless fermions, three spinless fermions, and two fermions with spins, respectively. For the simulation of two spinless fermions, with no additional leakage to other states out of the selected subspace, we can get a consistent simulation result with the experimental data, which further convinced the significant improvement of fidelity with the PEC method as shown in the introduction of Trotter errors Fig. 2i. For the simulation results for three spinless fermions and two fermions with spins, the population in other states out of the ideal subspace is introduced because of the Trotter errors. The leakage will be further increased by the noisy gates in the experiment, which needs to be mitigated by MLE and PS methods. As a result, the simulation results shown in Fig. 8c–j have obvious deviations from the experimental data, which can be considered as a result of the unexpected errors of gates such as crosstalk error.

Fig. 7: Experimental raw data compared with error-included simulation.
figure 7

We simulated the dynamics with noisy gates measured from QPT. Populations and fidelities of experimental raw data and simulation in (a, b) two spinless fermions, (c, d) three spinless fermions V = 0, (e, f) three spinless fermions V = 2J, (g, h) three spinless fermions U = 0 and (i, j) three spinless fermions U = 2J. The solid points represent the state populations measured from the experiment, and the dashed–dotted lines represent numerical simulation data with noisy gates. All the shadows in (a)–(j) represent the standard deviation of 1000 samples generated from the raw experimental data, using the bootstrapping method. k Trotter error rates from numerical simulation.

Fig. 8: Other state population of raw data and after applying PEC+MLE method.
figure 8

Populations in a two-qubit, b three-qubit V = 0, c three-qubit V = 2J, d four-qubit U = 0, e four-qubit U = 2J case.

The infidelities of two-qubit gates from experimental imperfections

We study the infidelities of the gates using numerical calculations on experimental imperfections. The infidelities are expected to mainly originate from the slow drift of laser frequency with respect to qubit energy level, intensity fluctuation, and slow frequency drifts of vibrational modes. On the other hand, infidelities from the spontaneous emission, solution imperfections, off-resonant coupling to the carrier, and crosstalk are expected to be smaller than 10−3. We note that SPAM errors should be negligible in our experiment, since we assume the initial state preparation is perfect and apply the scheme of detection-error correction56.

As shown in the main text, the average fidelities of two-qubit gates from quantum process tomography in two-, three-, and four-ion systems are 0.981(5), 0.976(6), and 0.973(5), respectively. We note that the infidelities of the two-qubit gates in our experiments are comparable to those in refs. 65,66,67,68. The typical coherence time of qubits is around 6 ms, and the dephasing time of the center of mass vibrational mode for a single ion is around 4 ms. We estimate around 1% infidelity comes from the slow drift of laser frequency close to 2 kHz; 0.8% from the frequency drifts of vibrational modes around 1.5 kHz; and 0.1% from 5% of Rabi-frequency fluctuations for two-qubit gates of the two-ion chain. With three and four ions, the durations of two-qubit gates are increased as shown in Table 1. Accordingly, the drifts of laser and mode frequencies become larger, which leads to larger infidelities for three- and four-ion chains. By stabilizing these slow drifts, the fidelities of the gates in multiple ions can be increased to over 0.99.

Table 1 Experimental parameters for trapped-ion systems with N = 2, 3 and 4 ions and corresponding two-qubit gates.

The anaysis of crosstalk errors

The crosstalk errors come from the spillover of an individually focused laser beam on the other ion. In the experiment, the inter-ion distance is around 5 μm and the beam waist of the individual beam is 1.5 μm. The measured ratios of Rabi frequencies from the left and the right ions to the center one in a three-ion chain are 0.024 and 0.011, respectively. We analyze the crosstalk error with the model that the spillover laser beams are applied to the neighboring ions with the same frequency and phase on the target ion with the different intensities quantified by the Rabi ratios.

In our model, the laser-ion interaction Hamiltonian with spillover laser beams can be written as

$$\begin{array}{l}\hat{H}_{{{{\rm{cross}}}}}\,=\,2\mathop{\sum}\limits_{j,m}{\eta }_{j,m}\cos {\omega }_{0}t\left({\hat{a}}_{m}^{{\dagger} }+{\hat{a}}_{m}\right){\hat{\sigma }}_{x}^{j}\\ \qquad\qquad \times \,\left[{{{\Omega }}}_{j}(t)\cos \left[\mu t+{\phi }_{j}(t)\right]+\underbrace{\gamma {\sum }_{i\in \langle i,j\rangle }{{{\Omega }}_{i}(t)\cos \left[\mu t+{\phi }_{i}(t)\right]}}_{\begin{array}{c}{{{\rm{laser}}}}\,{{{\rm{crosstalk}}}}\end{array}}\right],\end{array}$$
(12)

with γ 1 being the crosstalk ratio. Here all the symbols have the same meaning as those used in Ref. 51. The linearity of the constraints on the displacement of αj,m, j-th ion on m-th mode, guarantees that αj,m = 0 after considering laser crosstalk. The effect of the laser crosstalk manifests in the deviation of the geometric phase for the pair of j-th and \({j}^{{\prime}}\)-th ions, \({\theta }_{j,{j}^{{\prime} }}\) to the target values. Table 2 shows the average gate infidelities induced by the laser crosstalk. Finally, we note that the effect of cross-talk errors can be canceled by adjusting the subsequent gate phases or adding an additional gate to reduce these coherent errors69.

Table 2 Estimated average gate infidelities induced by laser intensity spillover for the two-qubit gates YYπ/8.