Quantum algorithm for the simulation of open-system dynamics and thermalization

The quantum open-system simulation is an important category of quantum simulation. By simulating the thermalization process at zero temperature, we can solve the ground-state problem of quantum systems. To realize the open-system evolution on the quantum computer, we need to encode the environment using qubits. However, usually the environment is much larger than the system, i.e., numerous qubits are required if the environment is directly encoded. In this paper, we propose a way to simulate open-system dynamics by reproducing reservoir correlation functions using a minimized Hilbert space. In this way, we only need a small number of qubits to represent the environment. To simulate the $n\mathrm{th}-\mathrm{order}$ expansion of the time-convolutionless master equation by reproducing up to $n$-time correlation functions, the number of qubits representing the environment is $\sim \lfloor \frac{n}{2}\rfloor {log}_2\left(N_{\omega }{N}_{\beta }\right)$. Here, $N_{\omega }$ is the number of frequencies in the discretized environment spectrum, and $N_{\beta }$ is the number of terms in the system-environment interaction. By reproducing two-time correlation functions, i.e., taking $n=2$, we can simulate the Markovian quantum master equation. In our algorithm, the environment on the quantum computer could be even smaller than the system.

Figure 1

(a) The simulated dynamics of the system is determined by the Hamiltonian $H=H_{\mathrm{S}}+H_{\mathrm{E}}+\alpha H_{\mathrm{I}}$ and the state of the environment ${\rho }_{\mathrm{E}}$. On the quantum computer, instead of directly simulating the environment, we use a Hilbert space with a much lower dimension to represent the environment and simulate the dynamics driven by the Hamiltonian $\tilde{H}=H_{\mathrm{S}}+{\tilde{H}}_{\mathrm{E}}+\alpha {\tilde{H}}_{\mathrm{I}}$ and the environment state ${\tilde{\rho }}_{\mathrm{E}}$. Because the environment on the quantum computer has a finite size, we may need to introduce dissipation in order to relax the environment. (b) The evolution driven by $\tilde{H}$ and the dissipation is realized using a quantum circuit on the quantum computer. To obtain the final state of the system ${\rho }_{\mathrm{S}}\left(t\right)$, qubits representing the system are prepared in the initial state of the system ${\rho }_{\mathrm{S}}\left(0\right)$, qubits representing the environment are prepared in the initial state of the environment ${\tilde{\rho }}_{\mathrm{E}}$, and then the evolution is implemented using computation operations, i.e., state preparation, quantum gates, and measurement operation. Ancillary qubits may be needed in the simulation, e.g., for implementing the dissipation.

Figure 2

Level scheme of the environment simulation. Each energy level of the environment represents a frequency with respect to the initial state of the environment. These levels are degenerate. The initial state is encoded as a pure state on the level with the frequency zero. The system can release energy into the environment or absorb energy from the environment via transitions between these energy levels. Transitions are caused by the interaction ${\tilde{H}}_I={\sum }_{\beta }{A}_{\beta }\otimes {\tilde{B}}_{\beta }$.

Figure 3

An example of dissipation caused by the environment. (a) An atom coupled to the free space in the vacuum state. (b) A black material is placed at a distance from the atom and absorbs photons.

Figure 4

Simulation of the second-order equation. (a) The system is coupled to the environment via local interaction. Without dissipation, an excitation in the environment leaves the interaction region but never disappears. (b) The dissipation is switched on in the region without interaction. The excitation disappears when it leaves the interaction region.

Figure 5

Numerical results for the wave packet $N^{-1}\tilde{B}\left(t\right)|v\rangle$, where $N=\parallel \tilde{B}\left(0\right)|v\rangle \parallel$, i.e., the wave packet is normalized at $t=0$. The probability in the state $|x\rangle$ is plotted. We take $N_{\omega }=1001, \delta \omega =0.001, \mathrm{\Gamma }=0.004$, and $\overline{\gamma }=0.002$. The conditional dissipation is introduced between $x=500$ and 800. The wave packet initializes at $x=0$ when $t=0$. From left to right, the blue, green, yellow, and red solid curves correspond to time $t=(200,400,550,700)\times 2\pi /N_{\omega }\delta \omega$, respectively. The wave packet is computed using the quantum trajectory approach [39, 40], and each curve is obtained with 1000 instances. For comparison, dashed curves denote the wave packet when the dissipation is turned off. The black curve represents the correlation function $e^{-\overline{\gamma }\left|s\right|}$ with $s=x/c$, which vanishes at about $x=300$. Therefore, the wave packet centered at $x>300$ does not contribute to the correlation function. The wave packet vanishes (see the arrow) after it enters the dissipation zone (marked in gray).

Figure 6

Probability in the ground state, $p_{\mathrm{g}}$. We take $\mathrm{\Delta }=1, \alpha \sqrt{a}=0.01, \overline{\gamma }=10, \delta \omega =0.02, N_{\omega }=401$. The dissipation of the environment is introduced using the condition reinitialization protocol (see the end of Sec. 7c). The probabilities $p_{\mathrm{g}}$ for zero temperature (blue circles) and finite temperature ($\beta =1$, red circles) are computed using the quantum trajectory approach [39, 40]. We take 1000 instances for zero temperature and 5000 instances for finite temperature. The dissipation zone is $21\le x\le 380$ for zero temperature; further moving the dissipation zone towards $x=0$ (e.g., taking the dissipation zone $11\le x\le 390$) will change the correlation function. The dissipation zone is $3\le x\le 398$ for finite temperature, which is chosen to obtain the best fit to the Lindblad equation of the thermalization. Black curves represent the result of the corresponding Lindblad equation of the thermalization.