Error-Mitigated Quantum Simulation of Interacting Fermions with Trapped Ions

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 correction [1, 2]. In this direction, quantum advantages with random gate sampling algorithms in superconducting systems have been demonstrated [3, 4]. The possibility of surpassing conventional classical computation in solving useful and practical problems has been seriously explored [5]. 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 computers [6, 7, 8, 9]. One of the critical issues in NISQ computation is to obtain results with high precision from quantum devices without quantum error correction [10, 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 devices [12, 13, 14, 15]. Representative error-mitigation methods include zero-noise extrapolation [12, 13, 14, 15, 16, 17], probabilistic error cancellation (PEC) [12, 13, 14, 15, 18, 19, 20, 21], virtual distillation [22], etc. 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 situations [12, 13, 20, 21]. However, it has been experimentally tested with only two-qubit systems of superconducting circuits and trapped ions [18, 19] and a multi-qubit system of the superconducting circuit via fitting a sparse error model motivated by the processor topology [23].

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 designs [6, 7, 8, 9]. In digital quantum simulation, the time evolution of a certain Hamiltonian is typically implemented by the Trotter-Suzuki expansion [24]. 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 model [25] is particularly interesting. This model is one of the fundamental models for interacting electrons, and it has been extensively studied in classical manners, either analytically [26, 27, 28] or numerically [29]. As experimental quantum technologies develop, understanding and insights for strongly-correlated quantum systems and high-temperature superconductivity [30] 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 platforms [31, 32, 33, 34]. 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 system [35].

In this paper, we experimentally study the possibility and limitations of PEC by applying the method to a four-qubit system with trapped ions [36, 37]. 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 conservation [38, 39, 15] 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.

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 and Fig. 1b for one-component and two-component fermions. The one-component fermions can be considered spinless ones. The Hamiltonian reads

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 $\lambda$ ($\lambda=\uparrow,\downarrow$) 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.

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

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 ${}^{171}{\rm Yb}^{+}$ ions for qubit systems. The ${}^{171}{\rm Yb}^{+}$ ions are confined in a linear Paul trap. The hyperfine levels in the ${}^{2}S_{1/2}$ ground-state manifold, i.e. $\left|F=0,m_{F}=0\right>$ and $\left|F=1,m_{F}=0\right>$, are encoded as the qubit $\left\{\left|0\right>,\left|1\right>\right\}$, which can be initialized to $\left|0\right>$ by optical pumping and disgustingly detected by fluorescence measurement via multi-channel photo-multiplier tube [40].

The qubit-system Hamiltonian becomes

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

with the Trotter step being $\delta t=t/M$. Figs. 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}_{\varphi}$, ${YY}_{\varphi}$ and ${ZZ}_{\varphi}$, and single-qubit rotations. By applying the amplitude-shaped [41] global and individual laser beams on the desired ions and driving the transverse motional modes [42, 43], we can realize the native Mølmer-Sørensen $YY$-gate (${YY}_{\varphi}$). The gate ${XX}_{\varphi}$ and ${ZZ}_{\varphi}$ are composed with ${YY}_{\varphi}$ 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 $U\gg V$ in most circumstances.

Here, we introduce the PEC error-mitigation method, which systematically recovers ideal expectation values of target quantum circuits with the experimental noisy gates [14]. 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 the trapped ion system, the error of the two-qubit gate is larger than the single-qubit one [19], and we can 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, to obtain the Pauli transfer matrices (PTMs) of the Mølmer-Sørensen $YY$ gates on different qubit pairs, i.e. $R_{YY_{\varphi}^{\rm ns}(m,n)}$ for the qubit pair $(Q_{m},Q_{n})$. Together with the ideal PTM $R_{YY_{\varphi}}$, we can formally define an error operator for each of the noisy gates, with the PTM being

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 $P_{i}\in\{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,

with the inverse-error decomposition

Here the quasi-probabilities $q_{i}(m,n)$ are real and satisfy $\sum_{i}q_{i}(m,n)=1$. The fact that $q_{i}(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

where $p_{i}(m,n)=C_{YY_{\varphi}(m,n)}^{-1}\left|q_{i}(m,n)\right|$ are well-defined probabilities with the cost $C_{YY_{\varphi}(m,n)}=\sum_{i=0}^{15}\left|q_{i}(m,n)\right|$, 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 $N_{g}$ number of two-qubit gates $YY_{\varphi}(m,n)$, with the $g$-th gate labeled as $G_{g}$, 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 $q_{i,g}$ being the quasi-probability of the two-qubit Pauli operator $P_{i}$. The full set of the random circuits to be implemented in the experimental device denoted as $\left\{{\mathcal{S}}_{\mathbf{i}}:{\mathbf{i}}=\left(i_{1},\ldots,i_{N_{g}}\right)\right\}$, contains $16^{N_{g}}$ different circuits, each of which is obtained by adding $P_{i_{g}}$ right after $G_{g}$. To obtain the ideal expectation values, each random circuit ${\mathcal{S}}_{\mathbf{i}}$ is assigned with a probability $\propto\prod_{g}\left|q_{i_{g},g}\right|$ and a sign ${\rm sgn}_{\mathbf{i}}\equiv\prod_{g}{\rm sgn}\left(q_{i_{g},g}\right)$. 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 $N_{s}$ random circuits by sampling $P_{i_{g}}$ from the probability distribution $\propto\left|q_{i_{g},g}\right|$ and obtain the expectation value $\langle\mu\rangle_{\mathbf{i}}$ of an observable $\mu$ for the circuit ${\mathcal{S}}_{\mathbf{i}}$ by averaging 300 repetitions. The ideal expectation value $\langle\mu\rangle$ is then calculated by

where the cost is given by $C\equiv\prod_{g=1}^{N_{g}}C_{g}$ with the cost of the $g$-th gate being $C_{g}\equiv\sum_{i_{g}=0}^{15}\left|q_{i_{g},g}\right|$.

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

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_{\varphi}$ 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$.

Errors caused by gate imperfection accumulate as the number of Trotter steps increases and the fidelity of evolved dynamics would be degraded. We mitigate the errors and improve the results of dynamics. Since it is a two-qubit system, only the one native $YY_{\pi/4}$ gate needs to be mitigated. The ideal $YY_{\pi/4}$ gate can be decomposed as $R_{{YY}_{\pi/4}}=\sum_{i=0}^{15}q_{i}R_{P_{i}}R_{YY^{\rm ns}_{\pi/4}}$ under Pauli error assumption, which is visualized as Fig. 2b. The experimental $YY^{\rm ns}_{\pi/4}$ gate is characterized by QPT. The experimental results of the PTMs of the $YY^{\rm ns}_{\pi/4}$ gate and the deviation from the ideal one obtained by $R_{YY^{\rm ns}_{\pi/4}}-R_{{YY}_{\pi/4}}$ are shown in Figs. 2c and d, where the average gate fidelity of $R_{YY^{\rm ns}_{\pi/4}}$ is $0.9811\pm 0.0047$ [44, 45]. Here, we assume the initial state preparation is perfect and apply the detection error correction method [46] for the measurement results, which makes state preparation and measurement (SPAM) errors negligible. The quasi-probabilities, i.e. the coefficients $q_{i}$ of the two-qubit Pauli operator $P_{i}$ in the decomposition of $R_{\mathcal{E}}^{-1}$, are shown in Fig. 2e. With the values of the quasi-probabilities $q_{i}$, we obtain the cost $C_{YY_{\pi/4}}=1.083$ for $R_{YY^{\rm ns}_{\pi/4}}$.

Then, we generate $N_{s}=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|ab\right>}$ for two spinless fermionic modes, where the observables are $\mu=\left|ab\right>\langle ab|$ 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\leq P_{\left|ab\right>}\leq 1$ and $\sum_{a,b}P_{\left|ab\right>}=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 Figs. 2g and h for the cases without and with error mitigation, respectively. We use the initial two-qubits state of $(\left|11\right>+\left|10\right>)/\sqrt{2}$, which contains the superposition of two occupied fermions and a single occupied fermion. The initial states simultaneously reveal the time evolutions of both single and double fermions, since no interference occurs between both dynamics. As shown in both Figs. 2g and h, the state $\left|11\right>$ and $\left|00\right>$ do not change and there is an oscillation between the state $\left|01\right>$ and $\left|10\right>$, which is the hallmark of fermion dynamics [32]. 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\pm 0.0021$ and $0.9975\pm 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_{YY^{\rm ns}_{\pi/4}}$. The comparisons between the simulation with the noisy PTM and experimental data clearly validate the Pauli error assumption and accuracy of measured $R_{YY^{\rm ns}_{\pi/4}}$ (see Appendix B). 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 Appendix A).

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|101\right>+\left|110\right>)/\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|110\right>,\left|101\right>,\left|011\right>\}$ 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 Figs. 3a, b, d, and e, dashed-dot lines show the results of ideal numerical simulation. When $V=0$, the initial state $\left|110\right>$ mostly evolves to the state $\left|101\right>$ because of the Pauli exclusion principle of fermions. On the other hand, the initial state $\left|101\right>$ can evolve to both $\left|110\right>$ and $\left|011\right>$ states. Therefore, in the beginning, the overall population of the state $\left|110\right>$ reduces, and the population of the state $\left|011\right>$ increases as shown in Figs. 3a and b. When $V=2J$, we can see the population of $\left|110\right>$ is not rapidly decreasing and the state $\left|101\right>$ is not populated compared to the case of $V=0$. It is because the interaction is effectively attractive in our experiment.

In experiment, we set $\varphi=\frac{J\delta t}{2}=\frac{\pi}{8}$ of the $YY_{\varphi}$ gates and apply $M=4$ Trotter steps for both scenarios. With nearest-neighbor interactions, each Trotter step consists of six $YY_{\varphi}$ gates and twenty single-qubit rotations. The whole quantum circuit consists of twenty-five two-qubit gates, where the additional one is for the initial state preparation, and eighty single-qubit gates. 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_{\pi/8}(1,2)$ and $YY_{\pi/8}(2,3)$, both of which require characterization and inverse-error decomposition. In the PEC scheme, both ${YY_{\pi/8}(1,2)}$ and ${YY_{\pi/8}(2,3)}$ gates are individually characterized by QPT and the inverse errors are decomposed independently. The average gate fidelities are $0.9779\pm 0.0088$ for $R_{YY^{\rm ns}_{\pi/8}(1,2)}$ and $0.9748\pm 0.0085$ for $R_{YY^{\rm ns}_{\pi/8}(2,3)}$ (see Appendix B).

For the case of $V=0$, the experimental results without and with error mitigation are shown in Figs. 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\pm 0.0161$ and $0.9942\pm 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 Appendix C). Similar conclusions are obtained for the scenario with the nearest-neighbor interaction of $V=2J$, as shown in Figs. 3d, e and f. The fitted population fidelities per gate without and with error mitigation are $0.9827\pm 0.0142$ and $0.9889\pm 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_{YY_{\pi/8}(1,2)}=1.157$ and $C_{YY_{\pi/8}(2,3)}=1.171$. Therefore, we sample more random circuits as $N_{s}=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 Figs. 3c and 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 Figs. 3c and 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.

We implement the simulation of four fermionic modes with spins by encoding qubits $Q_{1}Q_{3}$ and $Q_{2}Q_{4}$ 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|1001\right>+\left|1010\right>)/\sqrt{2}$, which is corresponding to the superposition state between $\left|\uparrow\right>_{1}\left|\downarrow\right>_{2}$ and $\left|\uparrow\right>_{1}\left|\downarrow\right>_{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|0101\right>,\left|0110\right>,\left|1001\right>,\left|1010\right>\}$ 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 ${\rm spin}=n_{j,\uparrow}-n_{j,\downarrow}$ and ${\rm charge}=n_{j,\uparrow}+n_{j,\downarrow}$ with $n_{j,\uparrow(\downarrow)}$ 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 excitation [47, 48, 35, 49]. 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 Figs. 2g and h. The state $\left|\uparrow\right>_{1}\left|\downarrow\right>_{2}$($\left|1001\right>$) is going back and forth with the other state $\left|\uparrow\right>_{2}\left|\downarrow\right>_{1}$($\left|0110\right>$). Likewise, the state $\left|\uparrow\right>_{1}\left|\downarrow\right>_{1}$($\left|1010\right>$) is oscillating with the other state $\left|\uparrow\right>_{2}\left|\downarrow\right>_{2}$($\left|0101\right>$). 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|\uparrow\right>_{1}\left|\downarrow\right>_{2}$, which is not affected by the on-site interaction, and the state $\left|\uparrow\right>_{1}\left|\downarrow\right>_{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_{\varphi}$ gates and twenty single-qubit rotations. The whole quantum circuit consists of twenty-five two-qubit gates, where the additional one is for the initial state preparation, and eighty single-qubit gates. 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_{\pi/8}(1,2)$, $YY_{\pi/8}(3,4)$, $YY_{\pi/8}(1,3)$, and $YY_{\pi/8}(2,4)$, individually and decompose their inverse errors independently. The average gate fidelities are $0.9755\pm 0.0096$ for $R_{YY^{\rm ns}_{\pi/8}(1,2)}$, $0.9706\pm 0.0085$ for $R_{YY^{\rm ns}_{\pi/8}(3,4)}$, $0.9720\pm 0.0092$ for $R_{YY^{\rm ns}_{\pi/8}(1,3)}$ and $0.9744\pm 0.0091$ for $R_{YY^{\rm ns}_{\pi/8}(2,4)}$ (see Appendix B).

For the case of $U=0$, the experimental results without and with error mitigation are shown in Figs. 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\pm 0.0426$ and $0.9993\pm 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 Figs. 4e, and f. In Fig. 4g, the fitted population fidelities per gate without and with error mitigation are $0.9551\pm 0.0285$ and $0.9925\pm 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.

We note that the PEC scheme alone does not recover the dynamics for the four fermionic modes as shown in Figs. 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 Appendix C). 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_{YY_{\pi/8}(1,2)}=1.211$, $C_{YY_{\pi/8}(3,4)}=1.228$, $C_{YY_{\pi/8}(1,3)}=1.228$ and $C_{YY_{\pi/8}(2,4)}=1.223$. The experimental Monte-Carlo sampling as $N_{s}=2000$ is not sufficient, which is limited by the time-consuming running of massive random circuits in physical hardware.

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, the drift of laser parameters and spectral crowding causes correlations resulting in remaining errors after PEC. As demonstrated in this work, combining with constraint methods can ease this issue. Recently, several proposals and experiments have been reported to tackle correlations such as using sparse Pauli-Lindblad models [23], learning-based error-mitigation [50], matrix product operator representation [51], 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) spectroscopy [52] and large-depth QAOA [53].

This work was supported by the Innovation Program for Quantum Science and Technology under Grants No. 2021ZD0301602, and the National Natural Science Foundation of China under Grants No. 92065205, No. 11974200, No. 12225507, No. 12088101, and No. 12204535.

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. Fig. 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.

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. 1 (d) and (e). Here up to six two-qubit Mølmer-Sørensen gates $YY_{\varphi}$ 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. 6 and Fig. 7. The quasi-probabilities for ${YY_{\varphi}}$ gates in the simulation of three spinless fermions and two fermions with spins are shown in Fig. 8, which decides the probabilities of choosing Pauli operations in the Monte Carlo method.

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 Figs. 9- 11 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. 2(i). 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. 10 and  11 have obvious deviations from the experimental data, which can be considered as a result of the unexpected errors of gates such as crosstalk error.