Simulating Prethermalization Using Near-Term Quantum Computers

Quantum simulation is one of the most promising scientific applications of quantum computers. Due to decoherence and noise in current devices, it is however challenging to perform digital quantum simulation in a regime that is intractable with classical computers. In this work, we propose an experimental protocol for probing dynamics and equilibrium properties on near-term digital quantum computers. As a key ingredient of our work, we show that it is possible to study thermalization even with a relatively coarse Trotter decomposition of the Hamiltonian evolution of interest. Even though the step size is too large to permit a rigorous bound on the Trotter error, we observe that the system prethermalizes in accordance with previous results for Floquet systems. The dynamics closely resemble the thermalization of the model underlying the Trotterization up to long times. We make our approach resilient to noise by developing an error mitigation scheme based on measurement and rescaling of survival probabilities, which is applicable to time-evolution problems in general. We demonstrate the effectiveness of the entire protocol by applying it to the two-dimensional $XY$ model and we numerically verify its performance with realistic noise parameters for superconducting quantum devices. Our proposal thus provides a route to achieving quantum advantage for relevant problems in condensed-matter physics.

Popular Summary

Current quantum computers can already perform certain tasks more efficiently than classical computers. These tasks, however, are of limited practical use. It therefore remains an open challenge to identify an actually important problem for which near-term quantum computers outperform classical computers. A promising candidate in this context is the simulation of quantum dynamics. There are currently two main challenges for this task: performing the simulation in a resource-efficient manner and mitigating errors due to noise. In our work, we show that both of these issues can be overcome in the near future.

It is difficult to exactly simulate continuous dynamics on a digital quantum computer. When this is done approximately, the difference between the real model and the simulated dynamics can be interpreted as a source of heating. We demonstrate that this heating can be made sufficiently weak such that it does not overshadow the dynamics that we want to observe. Moreover, we propose an error-mitigation protocol that allows us to get rid of the noise up to a certain simulation time.

Our approach provides an avenue to demonstrate useful quantum advantage on noisy devices. It can be used to simulate quantum systems in two or even more dimensions, which can be challenging with classical techniques such as the quantum Monte Carlo algorithm. As a result, this approach can contribute to answering important open questions in condensed-matter physics.

Figure 1

Different regimes of the dynamics of local observables depending on the Trotter step $\tau$ and the evolution time $T$. The black lines separate three regimes: bounded Trotter error (bottom), Floquet prethermalization (middle), and chaotic dynamics (top). The lower line scales as $\mathcal{O}(max\{{\tau }^{-p/(d+1)},N^{-1}{\tau }^{-p}\})$, following error bounds for the $p\mathrm{th}$-order Trotter decomposition with system size $N$. The upper line is determined by the Floquet heating time and scales as $e^{\mathcal{O}(1/\tau )}$. The blue shaded area indicates constant maximum circuit depth, relevant for noisy quantum computers. The gray area is excluded due to the constraint that $T\ge \tau$. The red shading highlights where the total time $T$ exceeds a system-dependent (pre)thermalization time scale $T_{\mathrm{th}}$. The prethermalized expectation value problem is experimentally accessible in the purple intersection.

Figure 2

(a) Prethermal plateau of the two-dimensional (2D) $XY$ model for system size $N=4\times 4$. The initial state is $|X+⟩$. The colored lines show the time averages of the mean-squared in-plane magnetization for different Trotter step sizes $\tau$, corresponding to different driving frequencies $\omega =2\pi /\tau$, ranging from $4J$ to $10J$. The large circles stand for the starting and end points of the plateaus according to Definition 4 with tolerance $ϵ=0.05$ and a maximum value of $t_{2}J$ of ${10}^3$. The black dashed line represents the value in the diagonal ensemble of the initial Hamiltonian. (b) Comparison of the value at the prethermal plateau with the values in the microcanonical and diagonal ensemble values of the initial $XY$ Hamiltonian and in the diagonal ensemble of the first-order Magnus expansion. The energy on the $x$ axis is given by the expectation value of the $XY$ Hamiltonian with regard to initial states. The system size is $4\times 4$. The driving frequency is as $\omega =8J$ and the plateau value is taken from the time average at $t=20/J$, which is on the prethermal plateau for all computed initial states with tolerance $ϵ=0.05$. For the microcanonical ensemble, we average over an energy window of width $\delta =0.5J$ in the $m_z=0$ subspace (see Appendix pp1).

Figure 3

(a) Quantum circuit to map the expectation value of a (unitary) observable onto a survival probability. The initial state is prepared with $V$, $U=U_1,U_2,U_3$, or $U_4$ is a single step in the Trotter decomposition, and $\mathcal{N}$ denotes a local noise channel. (b),(c) Dependence of $L_{\mathbb{1}}^{{\mathcal{N}}_p}(t)$ and $L_A^{{\mathcal{N}}_p}(t)$ on the circuit depth $D$ and system size $N$ in the presence of depolarizing noise with error probability $p=0.3\mathrm{\%}$. The initial state is $|\psi ⟩=|X+⟩$. The observable $A=4S_i^x{S}_{i+1}^x$ is a correlator in the center of the lattice. The black dashed lines represent the scaling predicted by Eq. (11).

Figure 4

(a) The mitigation error $s_A^{{\mathcal{N}}_p}$ for the range of data in Fig. 3 where $L_{\mathbb{1}}^{{\mathcal{N}}_p}(t)>0.01$. We choose this cutoff due to the limited number of trajectories in the simulation, which limits the significant digits. (b) The root-mean-square of $s$, $\sqrt{\sum _{\mathrm{data}}{s}^2/\sum _{\mathrm{data}}}$, evaluated over the window of circuit depth $[D-16,D+16]$ in (a).

Figure 5

The time average of $L_{{\sigma }_i^x{\sigma }_{i+1}^x,\psi }(t)$ at $Jt\approx 7.85$ with $\omega =8J$, or $\tau =\pi /4J$, corresponding to ten Trotter steps. The black crosses show the noiseless result. The red points were obtained by applying our error-mitigation strategy to noisy simulations with a single-qubit depolarization rate $p=0.3\mathrm{\%}$. Error bars indicate the statistic errors due to fluctuations of different Monte Carlo trajectories, propagated from the standard deviations of $L_A^{{\mathcal{N}}_p}(t)$ and $L_{\mathbb{1}}^{{\mathcal{N}}_p}(t)$. The system size $N=4\times 4$.

Figure 6

(a) Trotterization of the 2D $XY$ Hamiltonian. The circles represent qubits and the rectangles on the bounds represent two-qubit evolution gates on neighboring qubits. (b) Expectation value of the initial Hamiltonian during the Floquet evolution for the initial state $|X+⟩$ evolved. The parameters are the same as Fig. 2. (c),(d) Simulation of evolving Floquet $XY$ model with the Magnus expansion truncated to $n\mathrm{th}$ order. Their differences from Floquet evolution are plotted. The initial state is $|X+⟩$ and $A=m_x^2+m_y^2$. The system size is $N=4\times 3$. (c) $\omega =4J$, (d) $\omega =8J$. Note the difference in scale on the $y$ axis.

Figure 7

Simulation results for phase damping (a1)–(a4) and amplitude damping (b1)–(b4) noise and $p=0.3\mathrm{\%}$. The plotted quantities are the same as shown in Figs. 3 and 4.

Figure 8

The relative estimation error of $L_{\mathbb{1}}^{{\mathcal{N}}_p}(t)$ as a function of the number $n$ of Monte Carlo samples for circuit depth $D=40$ and $80$. The estimation error is defined as the standard deviation of the ensemble of expectation values from the trajectories divided by $\sqrt{n}$. The figures show the relative estimation error, i.e., the ratio of the estimation error to the estimated value (mean) of $L_{\mathbb{1}}^{{\mathcal{N}}_p}(t)$. The system size is $4\times 4$ and the initial state is $|X+⟩$.