Perturbative Quantum Simulation

Approximation based on perturbation theory is the foundation for most of the quantitative predictions of quantum mechanics, whether in quantum many-body physics, chemistry, quantum field theory, or other domains. Quantum computing provides an alternative to the perturbation paradigm, yet state-of-the-art quantum processors with tens of noisy qubits are of limited practical utility. Here, we introduce perturbative quantum simulation, which combines the complementary strengths of the two approaches, enabling the solution of large practical quantum problems using limited noisy intermediate-scale quantum hardware. The use of a quantum processor eliminates the need to identify a solvable unperturbed Hamiltonian, while the introduction of perturbative coupling permits the quantum processor to simulate systems larger than the available number of physical qubits. We present an explicit perturbative expansion that mimics the Dyson series expansion and involves only local unitary operations, and show its optimality over other expansions under certain conditions. We numerically benchmark the method for interacting bosons, fermions, and quantum spins in different topologies, and study different physical phenomena, such as information propagation, charge-spin separation, and magnetism, on systems of up to 48 qubits only using an $8+1$ qubit quantum hardware. We demonstrate our scheme on the IBM quantum cloud, verifying its noise robustness and illustrating its potential for benchmarking large quantum processors with smaller ones.

Figure 1

Dynamics simulation of interacting (a) bosons, (b) fermions, and (c) quantum spin systems with different topologies. Gray dashed lines in site-edge diagrams manifest subsystem partitioning for PQS. We group 8 qubits as a subsystem to simulate $N=16$-qubit quantum systems using most $5\times {10}^5$ samples. Solid lines represent exact results from direct simulation. (a) Quantum walk (QW) of spinless bosons on a 1D array in the large on site repulsion limit (see Sec. VI in [44] for the Hamiltonian [64]). Two identical bosons are initially excited at the center. (a1) Density spreading $⟨{\widehat{n}}_j⟩=⟨{\widehat{b}}_j^{†}{\widehat{b}}_j⟩$ with bosonic operators $\widehat{b}$ under time evolution. (a2) The density distribution at site 9—13 (labeling from left to right). The nearest-neighbor Lieb-Robinson bounds (dashed) capture the density spreading [6, 65, 66]. The inset shows the errors for the average density and the average density-density correlator ${\widehat{\rho }}_{ij}=⟨{\widehat{b}}_i^{†}{\widehat{b}}_j^{†}{\widehat{b}}_i{\widehat{b}}_j⟩$ with respect to the exact results. (a3) Boson spatial antibunching in QW. The normalized correlator ${\widehat{\rho }}_{ij}/{\widehat{\rho }}_{ij}^{\max }$ at different $t$ [67, 68]. (b) Separation of charge and spin density (CSD) in a 1D Fermi-Hubbard model $H=-J{\sum }_{j,\sigma }({\widehat{c}}_{j,\sigma }^{†}{\widehat{c}}_{j+1,\sigma }+\mathrm{H}.\mathrm{c}.)+U{\sum }_j{\widehat{n}}_{j,\uparrow }{\widehat{n}}_{j,\downarrow }+{\sum }_{j,\sigma }{h}_{j,\sigma }{\widehat{n}}_{j,\sigma }$ (${\widehat{c}}_{j,\sigma }$, ${\widehat{c}}_{j,\sigma }^{†}$: fermionic operators with spin $\sigma$, $U=J=0.5$) [68]. Left: two partitioning strategies for small and large on-site potential $U$. The initial state is the ground state of a noninteracting Hamiltonian with quarter filling ($N_{\uparrow }=N_{\downarrow }=2$), in which the CSD are generated in the middle of the chain at $t=0$ [69, 70, 71]. (b1) The separation of charge (blue square) and spin (red diamond) densities. We characterize the separation speed from the middle as ${\kappa }_{\pm }={\sum }_{j=1}^N|j-(N+1)/2|(⟨{\widehat{n}}_{j,\uparrow }⟩\pm ⟨{\widehat{n}}_{j,\downarrow }⟩)$ for charge ($+$) and spin ($-$) degrees of freedom with $⟨{\widehat{n}}_j⟩=⟨{\widehat{c}}_j^{†}{\widehat{c}}_j⟩$ ($N=8$). The inset shows the errors under evolution. (b2) The difference of CSD under evolution. The relative separation is initially set as 0. (c) Information propagation of correlated Ising spin clusters with power law decay interactions $H_l^{\mathrm{loc}}={\sum }_{ij}{J}_{ij}{\widehat{\sigma }}_{l,i}^x{\widehat{\sigma }}_{l,j}^x+h{\sum }_j{\widehat{\sigma }}_{l,j}^z$ ($J_{ij}=|i-j{|}^{-1}$) in the subsystems and interaction $V^{\mathrm{int}}={\widehat{\sigma }}_{1,N}^x{\widehat{\sigma }}_{2,1}^x$ on the boundary. The initial state is prepared as $|{\psi }_0⟩={\widehat{\sigma }}_8^x{|0⟩}^{\otimes N}$. (c1) The signal of quasiparticle excitations at different sites, where the propagation is faster than the nearest-neighbor Lieb-Robinson velocity (dashed) [65, 72, 73]. (c2) The dynamics of the correlation function $C_d=⟨{\widehat{\sigma }}_8^z{\widehat{\sigma }}_{8+d}^z⟩-⟨{\widehat{\sigma }}_8^z⟩⟨{\widehat{\sigma }}_{8+d}^z⟩$. The inset shows the errors for the averaged quasiparticle excitations’ density and correlation functions.

Figure 2

Dynamics simulation of 1D 48-site spin chains. The subsystem and interaction Hamiltonians are $H_l^{\mathrm{loc}}={\sum }_i{\widehat{\sigma }}_{l,i}^x{\widehat{\sigma }}_{l,i+1}^x+{\sum }_i{\widehat{\sigma }}_{l,i}^z$ and $V_l^{\mathrm{int}}=f_l{\widehat{\sigma }}_{l,N}^x{\widehat{\sigma }}_{l+1,1}^x$, respectively, and the interactions on the boundary are randomly generated from $[0,J/2]$. (a) The average magnetization (in blue) $(1/N){\sum }_i⟨{\widehat{\sigma }}_i^z⟩$ and nearest-neighbor correlation function (in red) $[1/(N-1)]{\sum }_i⟨{\widehat{\sigma }}_i^z{\widehat{\sigma }}_{i+1}^z⟩$, compared with the TEBD method as a benchmark. The inset illustrates the geometry of the spin systems and the partitioning strategy where we group 8 adjacent qubits as subsystems. (b) The errors for the average magnetization and correlation using $5\times {10}^5$ samples.

Figure 3

Implementation and results of the dynamical phase transition of 8 interacting spins. The initial state $|{\psi }_0⟩={|0⟩}^{\otimes 8}$ is evolved under an 8-site Ising Hamiltonian $H={\sum }_j{\widehat{\sigma }}_j^z{\widehat{\sigma }}_{j+1}^z+0.5{\sum }_j{\widehat{\sigma }}_j^x$ with $T=1$. (a) Quantum circuit for 8-qubit simulation based on first-order Trotterization with four steps $t_0=1/4$. (b) The circuit block for a single-step evolution for time $t$ with parallelization. (c) An example for the implementation of PQS to simulate an 8-qubit system with operations on the $4+1$ qubit. The circuit blocks are similar as that in (b) with 4 qubits. When a generalized operation is inserted into a Trotter step, we divide the step into two evolutions and insert the operation between that. (d) The magnetization and Loschmidt amplitude in ferromagnetic and paramagnetic phases. Here, the Loschmidt amplitude $\mathcal{G}(t)={|⟨{\psi }_0|e^{-iHt}|{\psi }_0⟩|}^2$ characterizes the dynamical echo back to the initial state [76], as an indicator of dynamical phase transition when it decreases to 0. We compare the results of exact simulation (dashed line), PQS (numerics, circle), PQS using IBMQ [5 qubits in (c), upper triangle], and the direct simulation using IBMQ [8 qubits in (a), lower triangle]. We also show the results using measurement error mitigation for PQS (solid square) and direct simulation (solid diamond).