Digital quantum computation of fermion-boson interacting systems

We introduce a method for representing the low-energy subspace of a bosonic field theory on the qubit space of digital quantum computers. This discretization leads to an exponentially precise description of the subspace of the continuous theory thanks to the Nyquist-Shannon sampling theorem. The method makes the implementation of quantum algorithms for purely bosonic systems as well as fermion-boson interacting systems feasible. We present algorithmic circuits for computing the time evolution of these systems. The complexity of the algorithms scales polynomially with the system size. The algorithm is a natural extension of the existing quantum algorithms for simulating fermion systems in quantum chemistry and condensed-matter physics to systems involving bosons and fermion-boson interactions and has a broad variety of potential applications in particle physics, condensed matter, etc. Due to the relatively small amount of additional resources required by the inclusion of bosons in our algorithm, the simulation of electron-phonon and similar systems can be placed in the same near-future reach as the simulation of interacting electron systems. We benchmark our algorithm by implementing it for a two-site Holstein polaron problem on an Atos Quantum Learning Machine quantum simulator. The polaron quantum simulations are in excellent agreement with the results obtained by exact diagonalization.

Figure 1

(a) Energy spectrum ${\tilde{E}}_n$ of the discrete Hamiltonian ${\tilde{H}}_h$ [Eq. (47)] for $N_x=64$ and 128. The dashed line is the continuous harmonic oscillator spectrum [Eq. (24)]. (b) Overlap between the eigenvectors $|{\tilde{\varphi }}_n\rangle$ of ${\tilde{H}}_h$, and the HG functions projected on the $\left\{\right|x_i\rangle \}$ subspace $|{\chi }_n\rangle$ [Eq. (34)]. For $n<N_{ph}$ where $N_{ph}$ is a cutoff number increasing with increasing $N_x, {\tilde{E}}_n\approx n+\frac{1}{2}$ and $|{\tilde{\varphi }}_n\rangle \approx |{\chi }_n\rangle$.

Figure 2

$|\left([\tilde{X},\tilde{P}]-i\right)|{\tilde{\varphi }}_n\rangle |$ versus $n$ for different values of $N_x$. For $n<N_{ph}$, the commutation operator satisfies $|\left([\tilde{X},\tilde{P}]-i\right)|{\tilde{\varphi }}_n\rangle |<\epsilon$, with $\epsilon \lesssim 10exp[-\left(0.51N_x-0.765N_{ph}\right)]$. Up to an exponentially small error, the algebra generated by $\{\tilde{X},\tilde{P}\}$ on the low-energy subspace $\left\{\right|{\tilde{\varphi }}_n\rangle {\}}_{n<N_{ph}}$ of $\tilde{\mathcal{H}}$ is isomorphic with the algebra generated by $\{X,P\}$ on low-energy subspace of the harmonic oscillator $\left\{\right|{\varphi }_n\rangle {\}}_{n<N_{ph}}$.

Figure 3

(a) The size of the discrete space $N_x$ increases linearly with the size of the low-energy subspace $N_{ph}$. The full (open) symbols are extracted from Fig. 2 for $\epsilon ={10}^{-7}$ ($\epsilon ={10}^{-3}$). The dashed line is the limit given by the WKB turning points, Eq. (33). (b) (Logarithmic scale) $L$, the half-width of the HG functions ${\varphi }_n$, defined by $1-{\int }_{-L}^L{\left|{\varphi }_n\left(x\right)\right|}^{2}dx=\epsilon$ as a function of $n$. $L$ scales approximately as $\sqrt{n}$. The dashed line corresponds to the WKB equation for the turning points $L=\sqrt{2n+1}$.

Figure 4

The circuit for $e^{-i\theta {\tilde{X}}_n}|x_n\rangle$. On each qubit, a phase shift gate is applied. The angles of the phase shift gates are determined by writing $x_n={\sum }_{r=0}^{n_x-1}{x}_n^r{2}^r$, where ${\left\{x_n^r\right\}}_{r=\overline{0,n_x-1}}$ takes binary values. A phase factor $exp(-i{2}^{n_x-1}\theta )$ accumulates at every Trotter step.

Figure 5

The circuit for $e^{-i\theta {\tilde{X}}_n^2}|x_n\rangle$ requires $n_x$ phase shift gates (one on each qubit) and $n_x(n_x-1)/2$ controlled phase shift gates. The angles of the phase shift gates are determined by writing ${(x_n-2^{n_x-1})}^2={\sum }_{r=0}^{n_x-1}{x}_n^r\left(2^{2r}-2^{n_x+r}\right)+{\sum }_{r<s}{x}_n^r{x}_n^s{2}^{r+s+1}+2^{2n_x-2}$, where ${\left\{x_n^r\right\}}_{r=\overline{0,n_x-1}}$ is the binary representation of $x_n$, i.e., $x_n={\sum }_{r=0}^{n_x-1}{x}_n^r{2}^r$. A phase factor $exp\left(i{2}^{n_x-2}\theta \right)$ accumulates at every Trotter step.

Figure 6

The circuit for $e^{-i\theta {\tilde{X}}_n{\tilde{X}}_m}|x_n\rangle |x_m\rangle$ requires $2n_x$ phase shift gates (one on each qubit of the two boson registers) and $n_x^2$ controlled phase shift gates. The angles of the phase shift gates are determined by writing $(x_n-2^{n_x-1})(x_m-2^{n_x-1})={\sum }_{r,s=0}^{n_x-1}{x}_n^r{x}_m^s{2}^{r+s}-{\sum }_{r=0}^{n_x-1}(x_n^r+x_m^r)2^{r+n_x-1}+2^{2n_x-2}$, where ${\left\{x_n^r\right\}}_{r=\overline{0,n_x-1}}$ is the binary representation of $x_n$, i.e., $x_n={\sum }_{r=0}^{n_x-1}{x}_n^r{2}^r$. A phase factor $exp\left(i{2}^{n_x-2}\theta \right)$ accumulates at every Trotter step.

Figure 7

Circuit describing $e^{-i\theta c_i^{†}{c}_i{\tilde{X}}_n}|i\rangle \otimes |x_n\rangle$. The phase shift angle is proportional to the boson coordinate $x_n-2^{n_x-1}={\sum }_{r=0}^{n_x-1}{x}_n^r{2}^r-2^{n_x-1}$, where ${\left\{x_n^r\right\}}_{r=\overline{0,n_x-1}}$ take binary values.

Figure 8

The circuit implementing $e^{-i\theta \left(c_i^{†}{c}_j+c_j^{†}{c}_i\right){\tilde{X}}_n}$ is similar to the circuit for the fermion hopping term (see Fig. 9 of Ref. [7] or Table A1 of Ref. [6]). The angle of the $z$-rotation gates is $\theta x_n$.

Figure 9

The circuit implementing $e^{-i\theta H_{xp}}$. The unitary operator $U$ transforms the position vectors $|x\rangle$ into the eigenvectors $|\nu \rangle$ of $H_{xp}$. For every state $|\nu \rangle$ a phase shift with the angle $E_{\nu }\theta$ is applied to an ancilla qubit prepared in state $|1\rangle$.

Figure 10

The energy (a) and quasiparticle weight (b) for the two-site Holstein polaron versus electron-phonon coupling strength. (c) The phonon-number distribution in the polaron state for different values of the coupling strength. The open symbols are computed using the exact diagonalization technique while the full small circles are computed using the QPE algorithm on a quantum simulator. The state of each harmonic oscillator was stored in a $n_x=6$ qubit register.