# Efficient Quantum Simulation of Electron-Phonon Systems by Variational Basis
State Encoder

Weitang Li  Tencent Quantum Lab, Tencent, Shenzhen, China  Jiajun Ren  College
of Chemistry, Beijing Normal Univerisity, Beijing, China  Sainan Huai  Tencent
Quantum Lab, Tencent, Shenzhen, China  Tianqi Cai  Tencent Quantum Lab,
Tencent, Shenzhen, China  Zhigang Shuai<EMAIL_ADDRESS>Department of
Chemistry, Tsinghua University, Beijing, China  School of Science and
Engineering, The Chinese University of Hong Kong, Shenzhen, China  Shengyu
Zhang<EMAIL_ADDRESS>Tencent Quantum Lab, Tencent, Shenzhen, China

###### Abstract

Digital quantum simulation of electron-phonon systems requires truncating
infinite phonon levels into $N$ basis states and then encoding them with qubit
computational basis. Unary encoding and the more compact binary/Gray encoding
are the two most representative encoding schemes, which demand $\order{N}$ and
$\order{\log{N}}$ qubits as well as $\order{N}$ and $\order{N\log{N}}$ quantum
gates respectively. In this work, we propose a variational basis state
encoding algorithm that reduces the scaling of the number of qubits and
quantum gates to both $\mathcal{O}(1)$. The cost for the scaling reduction is
a constant amount of additional measurement. The accuracy and efficiency of
the approach are verified by both numerical simulation and realistic quantum
hardware experiments. In particular, we find using 1 or 2 qubits for each
phonon mode is sufficient to produce quantitatively correct results across
weak and strong coupling regimes. Our approach paves the way for practical
quantum simulation of electron-phonon systems on both near-term hardware and
error-corrected quantum computers.

††preprint: APS/123-QED

Introduction Electron-phonon couplings are pervasive in quantum materials,
governing phenomena such as charge transport in semiconductors [1],
vibrational spectra [2], polaron formation [3], and superconductivity [4].
Classically, expensive numerical methods such as density matrix
renormalization group (DMRG) and quantum Monte-Carlo (QMC) are required to
accurately simulate electron-phonon systems due to the interior many-body
interaction [5, 6, 7, 8, 9]. Quantum computers hold promise for the simulation
of quantum systems with exponential speedup over classical computers [10]. In
the wake of the tremendous progress in the implementation of quantum computers
[11, 12] and the dawning of the noisy intermediate-scale quantum (NISQ) era
[13], how to solve electron-phonon coupling problems with quantum computers
has attracted a lot of research interest [14, 15, 16, 17].

A prominent problem for the digital quantum simulation of electron-phonon
systems is how to encode the infinite phonon states with finite quantum
computational basis states. The first step is usually truncating the infinite
phonon states into $N$ basis states $\\{\ket{m}\\}$ and then the second step
is encoding $\\{\ket{m}\\}$ into quantum computational basis $\\{\ket{n}\\}$.
The phonon basis states are usually the $N$ lowest harmonic oscillator
eigenstates or $N$ uniformly distributed grid basis states. There are two
established strategies to perform the encoding
$\\{\ket{m}\\}\mapsto\\{\ket{n}\\}$ [18, 19]. The first is unary encoding [20,
21], in which each $\ket{m}$ is encoded to $\ket{00\dots 1_{m}\dots 00}$, and
the total number of qubits required scales as $\order{N}$. The second is
binary encoding, in which each $\ket{m}$ is encoded to
$\prod_{i}\ket{\lfloor\frac{m}{2^{i}}\rfloor\mod 2}$ represented by
$\order{\log{N}}$ qubits [14, 15]. In terms of two-qubit gates required to
simulate quantum operators such as $\hat{b}^{\dagger}\pm\hat{b}$ and
$\hat{b}^{\dagger}\hat{b}$, unary encoding scales as $\order{N}$ and binary
encoding scales as $\order{N\log{N}}$ [19]. The features of unary encoding and
binary encoding are summarized in Table 1. Compared to the simulation of
electrons, the simulation of phonons consumes quantum resources in a much
faster manner, which becomes the bottleneck for efficient quantum simulation
of electron-phonon systems.

Table 1: Comparison of traditional encoding schemes and the proposed variational encoding in terms of encoding formula, the number of qubits $N_{\rm{qubit}}$ required and the number of quantum gates $N_{\rm{gate}}$ required to simulate common phonon operators such as $\hat{b}^{\dagger}\pm\hat{b}$ and $\hat{b}^{\dagger}\hat{b}$ .  Scheme | Formula | $N_{\rm{qubit}}$ | $N_{\rm{gate}}$
---|---|---|---
Unary | $\ket{m}\mapsto\ket{00\dots 1_{m}\dots 00}$ | $\order{N}$ | $\order{N}$
Binary | $\ket{m}\mapsto\prod_{i}\ket{\lfloor\frac{m}{2^{i}}\rfloor\mod 2}$ | $\order{\log{N}}$ | $\order{N\log{N}}$
Variational | $\sum_{m}C_{mn}\ket{m}\mapsto\ket{n}$ | $\order{1}$ | $\order{1}$

In this work, we propose a new basis encoding scheme called variational
encoding. Variational encoding maps linear combinations of $\ket{m}$ that are
most entangled to the simulated system into the computational basis, i.e.
$\sum_{m}C_{mn}\ket{m}\mapsto\ket{n}$, where $C_{mn}$ is determined by
variational principle. The advantage of our approach is that, by encoding only
the most entangled states and discarding the ones with little entanglement,
the size of $\\{\ket{n}\\}$ can be made irrelevant to the size of
$\\{\ket{m}\\}$. In other words, the number of qubits required scales as
$\order{1}$. Consequently, the scaling for the number of gates is also
$\order{1}$. Variational encoding is best suited to work in combination with
variational quantum algorithms such as variational quantum eigensolver (VQE)
[22, 23] and variational quantum dynamics (VQD) [24, 25]. Besides, the
variational encoding is also compatible with Trotterized time evolution and
quantum phase estimation (QPE) [10, 26, 27]. Numerical simulation and
experiments on realistic quantum hardware based on the Holstein model and
spin-boson model shows that using 1 or 2 qubits for each phonon mode is
typically sufficient for highly accurate results even in the strong coupling
regime.

Variational Basis State Encoder Encoder coefficients $C$ are determined by
variational principle for both static and dynamic cases. We start the
derivation using parameterized quantum circuit (PQC) and discuss how to
incorporate the variational encoder in Trotterized time evolution and QPE
later on. We use atomic units throughout the paper. More details for the
derivation can be found in the Appendix.

For each phonon mode $l$, encoded by $N_{l}$ qubits, define the variational
encoder $\hat{B}[l]$

$\hat{B}[l]=\sum_{m}\sum_{n=1}^{2^{N_{l}}}C[l]_{mn}\ket{n}_{l}\tensor{[}_{l}]{\bra{m}}{}\
,$ (1)

with orthonormal constraint $\hat{B}[l]\hat{B}[l]^{\dagger}=\hat{I}$. The
original Hamiltonian in $\ket{m}$ basis can then be encoded to $\ket{n}$ basis
$\hat{\tilde{H}}=\prod_{l}\hat{B}[l]\hat{H}\prod_{l}\hat{B}[l]^{\dagger}$.
Suppose the quantum circuit is parameterized by
$\ket{\phi}=\prod_{k}e^{i\theta_{k}\hat{R}_{k}}\ket{\phi_{0}}$, and then the
ground state Lagrangian with multipliers $\lambda_{lnn^{\prime}}$ is

$\mathcal{L}=\braket{\phi}{\hat{\tilde{H}}}{\phi}+\sum_{lnn^{\prime}}\lambda_{lnn^{\prime}}(\sum_{m}C[l]_{mn}C[l]_{mn^{\prime}}-\delta_{nn^{\prime}})\
.$ (2)

Taking the derivative with respect to $\theta_{k}$ leads to traditional VQE
with encoded Hamiltonian $\hat{\tilde{H}}$. Taking the derivative with respect
to $C[l]$ and setting it to 0 yields

$(1-\hat{P}[l])\braket{\phi}{\hat{\tilde{H}}^{\prime}[l]}{\phi}=0\ ,$ (3)

with projector $\hat{P}[l]=\hat{B}[l]^{\dagger}\hat{B}[l]$ and the half
encoded Hamiltonian $\hat{\tilde{H}}^{\prime}[l]=\prod_{k\neq
l}\hat{B}[k]\hat{H}\prod_{k}\hat{B}[k]^{\dagger}$. In practice, $\theta_{k}$
and $C[l]$ are solved iteratively until convergence. In the following, this
iteration is termed macro-iteration to avoid confusion with VQE iteration.

Next, we discuss the measurement required to solve $C[l]$ from Eq. 3. Suppose
the Hamiltonian can be written as $\hat{H}=\sum_{x}^{M}\hat{h}_{x}$ and
$\hat{h}_{x}=\prod_{k}\hat{h}[k]_{x}$, where $M$ is the total number of terms
in the Hamiltonian and $\hat{h}[k]_{x}$ acts on the $k$th degree of freedom.
The measurement of $\braket{\phi}{\hat{\tilde{H}}^{\prime}[l]}{\phi}$ boils
down to that of
$\braket{\phi}{n}_{l}\tensor{[}_{l}]{\bra{n^{\prime}}}{}\prod_{k\neq
l}\hat{\tilde{h}}[k]_{x}\ket{\phi}$, where
$\hat{\tilde{h}}[k]_{x}=\hat{B}[k]\hat{h}[k]_{x}\hat{B}[k]^{\dagger}$ is the
encoded local operator. For electron degree of freedom a dummy encoder
$\hat{B}[k]=\hat{I}$ is used for notational simplicity. The number of
additional measurements for the update of $C[l]$ is thus $\order{2^{N_{l}}M}$,
which is polynomial to the system size and does not increase with $N$. If the
number of phonon modes is assumed to be linear with $M$ and each $C[l]$ is
updated independently, then the total number of measurements for all $C[l]$ is
$\order{2^{N_{l}}M^{2}}$. The measurement overhead increases exponentially
with $N_{l}$. Due to arguments presented later, $N_{l}$ is usually small and
does not increase with system size. From numerical experiments, we find
$N_{l}\leq 2$ is sufficient to produce excellent results.

For time-dependent problems, it is convenient to define
$\ket{\psi}=\prod_{l}\hat{B}[l]^{\dagger}\ket{\phi}$ and use $\Theta_{K}$
denote both $\theta_{k}$ and $C[l]$. The Lagrangian with multipliers
$\lambda_{lnn^{\prime}}$ and $\gamma_{lnn^{\prime}}$ is then

$\displaystyle\mathcal{L}$
$\displaystyle=|i\sum_{K}\partialderivative{\ket{\psi}}{\Theta_{K}}\dot{\Theta}_{K}-\hat{H}\ket{\psi}|^{2}$
(4)
$\displaystyle\quad+\sum_{lnn^{\prime}}\lambda_{lnn^{\prime}}\Re{\sum_{m}C[l]^{*}_{mn}\dot{C}[l]_{mn^{\prime}}}$
$\displaystyle\quad+\sum_{lnn^{\prime}}\gamma_{lnn^{\prime}}\Im{\sum_{m}C[l]^{*}_{mn}\dot{C}[l]_{mn^{\prime}}}\
.$

The constraints ensure that $C[l]_{mn}$ remains orthonormal during time
evolution. Similar to the ground state problem, the equation of motion for
$\theta_{k}$ is the same as vanilla VQD with encoded Hamiltonian
$\hat{\tilde{H}}$

$\sum_{j}\Re{\partialderivative{\bra{\phi}}{\theta_{k}}\partialderivative{\ket{\phi}}{\theta_{j}}}\dot{\theta}_{j}=\Im{\partialderivative{\bra{\phi}}{\theta_{k}}\hat{\tilde{H}}\ket{\phi}}\
.$ (5)

The equation of motion for $C[l]$ reads

$i\rho[l]\dot{C}[l]^{*}=(1-\hat{P}[l])\braket{\phi}{\hat{\tilde{H}}^{\prime}[l]}{\phi}\
,$ (6)

where $\rho[l]_{nn^{\prime}}=\Tr{\braket{\phi}{n}\braket{n^{\prime}}{\phi}}$
is the reduced density matrix for the $N_{l}$ qubits of $\ket{\phi}$. Eq. 3
represents a $\dot{C}[l]=0$ stationary point during real and imaginary time
evolution. The measurement cost is the same as the ground state algorithm.

The VQD step described by Eq. 5 can be naturally replaced by a Suzuki-Trotter
time evolution step
$e^{-i\hat{\tilde{H}}\tau}\approx\prod_{x}^{M}e^{-i\hat{\tilde{h}}_{x}\tau}$
on a digital quantum simulator, so that Hamiltonian simulation is performed
via Trotterized time evolution instead of VQD. To update $C[l]$ based on Eq.
6, measurements on the circuit should be performed for every or every several
Trotter steps. The variationally encoded ground state can then be prepared by
adiabatic state preparation, whose energy is accessible by QPE using
$\hat{\tilde{H}}$.

It is instructive to observe that if the variational basis encoder is viewed
as a wavefunction ansatz $\ket{\psi}$, then the algorithm proposed in this
work can be viewed as a generalization for the local basis optimization method
for DMRG [28, 29], or a special case of the recently proposed quantum-
classical hybrid tensor network [30]. Thus, $\hat{B}[l]$ captures the
$2^{N_{l}}$ phonon states that are most entangled with the rest of the system.
For local Hamiltonian obeying the area law, the entanglement entropy between
one phonon mode and the rest of the system $S$ is a constant [31].
Consequently, $|\braket{\psi}{\Psi}|^{2}$, the fidelity between the
approximated encoded state and the target state has a lower bound of
$\frac{2^{N_{l}}}{e^{S}}$, which lays the theoretical foundation for the
effectiveness of the variational encoding approach to ground state and low-
lying excited state problems.

Simulations We first show numerical simulation results on a noiseless
simulator and verify the algorithm on a superconducting quantum computer at
the end of the section. The variational basis state encoder is first tested
for VQE simulation of the one-dimensional Holstein model [32, 33]

$\hat{H}=-\sum_{\langle
i,j\rangle}V\hat{a}^{\dagger}_{i}\hat{a}_{j}+\sum_{i}\omega\hat{b}^{\dagger}_{i}\hat{b}_{i}+\sum_{i}g\omega\hat{a}^{\dagger}_{i}\hat{a}_{i}(\hat{b}^{\dagger}_{i}+\hat{b}_{i})\
.$ (7)

where $V$ is the hopping coefficient, $\langle i,j\rangle$ denotes nearest
neighbour pairs with periodic boundary condition, $\omega$ is the vibration
frequency and $g$ is dimensionless coupling constant. In the following, we
assume $V=\omega=1$ and adjust $g$ for different coupling strengths. We
consider a 3-site system corresponding to $3(N_{l}+1)$ qubits. We use binary
encoding to represent traditional encoding approaches. Unary encoding is
expected to produce similar results with binary encoding only with different
quantum resource budgets. The ansatz used and more details of the simulation
are included in the Appendix.

Figure 1:  Numerical simulation results for the ground state of the Holstein
model. (a) Ground state energy by numerically exact DMRG, binary encoding, and
variational encoding with different coupling strength $g$; (b) Convergence of
ground state energy with respect to the macro-iteration for variational
encoding; (c) Ground state energy error for the variational encoding method at
different numbers of phonon basis states $N$; (d) The singular values for the
Schmidt decomposition between the last phonon mode and the rest of the system.

We first compare the accuracy of the variational encoding and the binary
encoding with $N_{l}=1$. It is clear from Fig. 1(a) that variational encoding
is significantly more accurate than binary encoding, especially at the strong
coupling regime. Within the setup, binary encoding uses only two phonon basis
states to describe each phonon mode, yet the variational encoding is allowed
to use up to 32 phonon basis states before combining them into the most
entangled states. We note that the quantum circuit used for variational
encoding and binary encoding is essentially the same. The number of macro-
iterations to determine $C[l]$ is found to be rather small, as shown in Fig.
1(b). Fully converged results are obtained within 5 iterations. In Fig. 1(c)
we show more details of the error for the variational approach. The simulation
error typically decreases exponentially with respect to the number of phonon
levels $N$ included in $C[l]$. It is worth noting that quantum computational
resources, including the number of qubits, the number of gates in the circuit,
and the number of measurements remain constant when $N$ is increased from 2 to
32. Furthermore, by using 2 qubits to encode each mode, it is possible to
further reduce the error at the $N\rightarrow\infty$ limit. When $g=3.0$, the
error is not sensitive to $N_{l}$, which implies that the error is dominated
by other sources such as limitations of the ansatz, instead of the small
$N_{l}$. Fig. 1(d) shows the singular values for the Schmidt decomposition
between the last phonon mode and the rest of the system by DMRG. The
exponential decay ensures the fast convergence of $N_{l}$. The von Neumann
entropy $S$ for the 3 systems is found to be 0.01, 0.25, and 0.65
respectively. We also note the $g=1.5$ case has the largest 3rd singular
value, which explains why setting $N_{l}=2$ significantly reduces the $g=1.5$
error in Fig. 1(c).

We now turn to the spin-relaxation dynamics of the spin-boson model [34],
described by the Hamiltonian

$\hat{H}=\frac{\epsilon}{2}\hat{\sigma}_{z}+\Delta\hat{\sigma}_{x}+\sum_{j}g_{j}\omega_{j}\hat{\sigma}_{z}(\hat{b}^{\dagger}_{j}+\hat{b}_{j})+\sum_{j}\omega_{j}\hat{b}^{\dagger}_{j}\hat{b}_{j}\
.$ (8)

where $\epsilon$ is the eigenfrequency and $\Delta$ is the tunnelling rate.
The coupling term has a similar form with Eq. 7 and is more commonly written
as $\sum_{j}c_{j}\hat{\sigma}_{z}\hat{x}_{j}$. For systems in the condensed
phase the coupling is usually characterized by the spectral density function
$\mathcal{J}(\omega)=\frac{\pi}{2}\sum_{j}\frac{c_{j}^{2}}{\omega_{j}}\delta(\omega-\omega_{k})$.
In the following we assume $\epsilon=0$ and $\Delta=1$. We first use VQD for
the simulation and discuss Trotterized time evolution at last. The variational
Hamiltonian ansatz [35] with 3 layers is used if not otherwise specified.

Figure 2:  Numerical simulation results for the spin-relaxation dynamics of
the spin-boson model. (a) Comparison between binary encoding with different
numbers of phonon basis states and variational encoding for a 1-mode spin-
boson model; (b) Variational encoding with different numbers of encoding
qubits $N_{l}$ for a 2-mode spin-boson model; (c) Comparison between binary
encoding and variational encoding for an 8-mode spin-boson model with sub-
Ohmic spectral density; (d) Trotterized time evolution with variational
encoding based on a 1-mode spin-boson model.

The performance of variational encoding and binary encoding is first compared
based on a 1-mode spin-boson model at the strong coupling ($\omega=1$ and
$g=3$) regime, shown in Fig. 2(a). Variational encoding with $N_{l}=1$
generates much more accurate dynamics than binary encoding with fewer qubits
and quantum gates. The simulation of binary encoding with $N_{l}>4$ is
prohibited by the deep circuit depth in the ansatz. The variational encoding
scheme is exceptionally efficient for this 1-mode model because Schmidt
decomposition guarantees that 2 variational bases for the phonon mode are
sufficient to exactly represent the system. In Fig. 2(b) a 2-mode model with
$\omega_{j}=\frac{1}{2},1$ and $g_{j}=\frac{1}{2},1$ is used. Variational
encoding with $N_{l}=1$ is accurate at $t<2$ but as the entanglement builds up
the dynamics deviate from the exact solution. Increasing $N_{l}$ to 2
effectively eliminates the error. Next, we move on to a more challenging model
with 8 modes, in which $\omega$ and $g$ are determined by discretizing a sub-
Ohmic spectral density
$\mathcal{J}(\omega)=\frac{\pi}{2}\alpha\omega^{s}\omega_{c}^{1-s}e^{-\omega/\omega_{c}}$
following the prescription in the literature [36]. The parameters are
$s=\frac{1}{4}$, $\omega_{c}=4$ and $\alpha=10$. As illustrated in Fig. 2(c)
variational encoding with $N_{l}=1$ captures the localization behavior yet
binary encoding with $N_{l}=1$ completely fails. The number of layers in the
variational Hamiltonian ansatz is 8 and 32 for variational and binary encoding
respectively. Fig. 2(d) demonstrates the possibility to incorporate
variational basis state encoder into Trotterized time evolution with
$\omega=g=1$ and $N_{l}=1$. The measurement and the evolution of $C[l]$ are
performed at each Trotter step.

Figure 3:  Quantum hardware experiments for the ground state energy of the
Holstein model with variational basis state encoder. (a) 3 qubits out of 9
qubits of a superconducting quantum computer and a one-parameter circuit are
used for the simulation; (b) Ground state energy by binary encoding and
variational encoding; (c) Convergence of ground state energy with respect to
the macro-iteration for variational encoding.

Lastly, we verify the accuracy and efficiency of the variational encoder
approach on a superconducting quantum computer. We consider the ground state
problem of a 2-site Holstein model described by Eq. 7 with $g=3$ and
$N_{l}=1$. The two electronic sites are represented by 1 qubit and the total
number of qubits for the system is thus 3. The quantum circuit for the
simulation is depicted in Fig. 3(a). The electronic degree of freedom is
mapped to the second qubit, and the two phonon modes are mapped to the first
and the third qubits respectively. There is one parameter to be determined by
VQE in the circuit and the same ansatz is used for both binary encoding and
variational encoding. More simulation details can be found in the Appendix. In
Fig. 3(b) we show the ground state energy by variational encoding from weak to
strong coupling, in analog to Fig. 1(a). The simulator result is based on the
parameterized quantum circuit described in Fig. 3(a) without considering gate
noise and measurement uncertainty. The results in Fig. 3(b) are consistent
with that in Fig. 1(a). The residual error is dominated by the intrinsic gate
noise in the quantum computer. In Fig. 3(c) we show the convergence with
respect to the macro-iteration for variational encoding. The algorithm is
resilient to the presence of quantum noise and measurement uncertainty. The
convergent energy is reached within 5 iterations.

Conclusion We proposed variational basis state encoder to encode phonon basis
states into quantum computational states for efficient quantum simulation of
electron-phonon systems. The proposed variational encoding approach requires
only $\order{1}$ qubits and $\order{1}$ quantum gates, which is significantly
better than traditional encoding schemes. The algorithm enables quantum
simulation of electron-phonon systems with smaller quantum processors and
shallower circuits. The additional measurement required to implement the
approach is found to be also $\order{1}$ with respect to the number of phonon
basis states and it scales quadratically with the number of Pauli strings in
the Hamiltonian. The accuracy of the approach is ensured by the finite
entanglement entropy between one phonon mode and the rest of the system in
common electron-phonon systems. Variational basis state encoder most naturally
works with variational quantum algorithms and is compatible with Trotterized
time evolution, adiabatic state preparation, and QPE. Numerical simulation and
quantum hardware experiments based on VQE of the Holstein model and dynamics
of the spin-boson model indicates that variational encoding is more accurate
and resource-efficient than traditional encoding methods. In particular, using
1 or 2 qubits to represent each phonon mode is sufficient for accurate
simulation even at the strong coupling regime where $N=64$ phonon basis states
are involved. The approach could also be extended to other quantum simulation
problems involving an infinite or large local Hilbert space.

We thank Jinzhao Sun and Shixin Zhang for helpful discussions. This work is
supported by the National Natural Science Foundation of China through grand
numbers 22273005 and 21788102. This work is also supported by Shenzhen Science
and Technology Program.

## References

  * Bardeen and Shockley [1950] J. Bardeen and W. Shockley, Deformation potentials and mobilities in non-polar crystals, Phys. Rev. 80, 72 (1950).
  * Cardona and Thewalt [2005] M. Cardona and M. L. W. Thewalt, Isotope effects on the optical spectra of semiconductors, Rev. Mod. Phys. 77, 1173 (2005).
  * Salamon and Jaime [2001] M. B. Salamon and M. Jaime, The physics of manganites: Structure and transport, Rev. Mod. Phys. 73, 583 (2001).
  * Pickett [1989] W. E. Pickett, Electronic structure of the high-temperature oxide superconductors, Rev. Mod. Phys. 61, 433 (1989).
  * Jeckelmann and White [1998] E. Jeckelmann and S. R. White, Density-matrix renormalization-group study of the polaron problem in the Holstein model, Phys. Rev. B 57, 6376 (1998).
  * Nowadnick  _et al._ [2012] E. A. Nowadnick, S. Johnston, B. Moritz, R. T. Scalettar, and T. P. Devereaux, Competition between antiferromagnetic and charge-density-wave order in the half-filled Hubbard-Holstein model, Phys. Rev. Lett. 109, 246404 (2012).
  * Mishchenko  _et al._ [2015] A. S. Mishchenko, N. Nagaosa, G. De Filippis, A. de Candia, and V. Cataudella, Mobility of Holstein polaron at finite temperature: An unbiased approach, Phys. Rev. Lett. 114, 146401 (2015).
  * Cai  _et al._ [2021] X. Cai, Z.-X. Li, and H. Yao, Antiferromagnetism induced by bond Su-Schrieffer-Heeger electron-phonon coupling: A quantum Monte Carlo study, Phys. Rev. Lett. 127, 247203 (2021).
  * Li  _et al._ [2021] W. Li, J. Ren, and Z. Shuai, A general charge transport picture for organic semiconductors with nonlocal electron-phonon couplings, Nat. Commun. 12, 4260 (2021).
  * Lloyd [1996] S. Lloyd, Universal quantum simulators, Science 273, 1073 (1996).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell, _et al._ , Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019).
  * Xu  _et al._ [2022] S. Xu, Z.-Z. Sun, K. Wang, L. Xiang, Z. Bao, Z. Zhu, F. Shen, Z. Song, P. Zhang, W. Ren, _et al._ , Digital simulation of non-abelian anyons with 68 programmable superconducting qubits, arXiv preprint arXiv:2211.09802 (2022).
  * Preskill [2018] J. Preskill, Quantum Computing in the NISQ era and beyond, Quantum 2, 79 (2018).
  * Macridin  _et al._ [2018] A. Macridin, P. Spentzouris, J. Amundson, and R. Harnik, Electron-phonon systems on a universal quantum computer, Phys. Rev. Lett. 121, 110504 (2018).
  * Ollitrault  _et al._ [2020] P. J. Ollitrault, G. Mazzola, and I. Tavernelli, Nonadiabatic molecular quantum dynamics with quantum computers, Phys. Rev. Lett. 125, 260511 (2020).
  * Lee  _et al._ [2022] C.-K. Lee, C.-Y. Hsieh, S. Zhang, and L. Shi, Variational quantum simulation of chemical dynamics with quantum computers, J. Chem. Theory and Comput. 18, 2105 (2022).
  * Wang  _et al._ [2022] Y. Wang, J. Ren, W. Li, and Z. Shuai, Hybrid quantum-classical boson sampling algorithm for molecular vibrationally resolved electronic spectroscopy with Duschinsky rotation and anharmonicity, J. Phys. Chem. Lett. 13, 6391 (2022).
  * Sawaya  _et al._ [2020] N. P. Sawaya, T. Menke, T. H. Kyaw, S. Johri, A. Aspuru-Guzik, and G. G. Guerreschi, Resource-efficient digital quantum simulation of $d$-level systems for photonic, vibrational, and spin-$s$ hamiltonians, npj Quantum Inf. 6, 1 (2020).
  * Di Matteo  _et al._ [2021] O. Di Matteo, A. McCoy, P. Gysbers, T. Miyagi, R. Woloshyn, and P. Navrátil, Improving Hamiltonian encodings with the Gray code, Phys. Rev. A 103, 042405 (2021).
  * Somma  _et al._ [2003] R. Somma, G. Ortiz, E. Knill, and J. Gubernatis, Quantum simulations of physics problems, Int. J. Quantum Inf. 1, 189 (2003).
  * Miessen  _et al._ [2021] A. Miessen, P. J. Ollitrault, and I. Tavernelli, Quantum algorithms for quantum dynamics: a performance study on the spin-boson model, Phys. Rev. Res. 3, 043212 (2021).
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’brien, A variational eigenvalue solver on a photonic quantum processor, Nat. Commun. 5, 1 (2014).
  * McClean  _et al._ [2016] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, The theory of variational hybrid quantum-classical algorithms, New J. Physics 18, 023023 (2016).
  * Li and Benjamin [2017] Y. Li and S. C. Benjamin, Efficient variational quantum simulator incorporating active error minimization, Phys. Rev. X 7, 021050 (2017).
  * Yuan  _et al._ [2019] X. Yuan, S. Endo, Q. Zhao, Y. Li, and S. C. Benjamin, Theory of variational quantum simulation, Quantum 3, 191 (2019).
  * Abrams and Lloyd [1999] D. S. Abrams and S. Lloyd, Quantum algorithm providing exponential speed increase for finding eigenvalues and eigenvectors, Phys. Rev. Lett. 83, 5162 (1999).
  * Aspuru-Guzik  _et al._ [2005] A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon, Simulated quantum computation of molecular energies, Science 309, 1704 (2005).
  * Zhang  _et al._ [1998] C. Zhang, E. Jeckelmann, and S. R. White, Density matrix approach to local Hilbert space reduction, Phys. Rev. Lett. 80, 2661 (1998).
  * Guo  _et al._ [2012] C. Guo, A. Weichselbaum, J. von Delft, and M. Vojta, Critical and strong-coupling phases in one-and two-bath spin-boson models, Phys. Rev. Lett. 108, 160401 (2012).
  * Yuan  _et al._ [2021] X. Yuan, J. Sun, J. Liu, Q. Zhao, and Y. Zhou, Quantum simulation with hybrid tensor networks, Phys. Rev. Lett. 127, 040501 (2021).
  * Eisert  _et al._ [2010] J. Eisert, M. Cramer, and M. B. Plenio, Colloquium: Area laws for the entanglement entropy, Rev. Mod. Phys. 82, 277 (2010).
  * Holstein [1959a] T. Holstein, Studies of polaron motion: Part I. the molecular-crystal model, Ann. Phys. 8, 325 (1959a).
  * Holstein [1959b] T. Holstein, Studies of polaron motion: Part II. the “small” polaron, Ann. Phys. 8, 343 (1959b).
  * Leggett  _et al._ [1987] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. Fisher, A. Garg, and W. Zwerger, Dynamics of the dissipative two-state system, Rev. Mod. Phys. 59, 1 (1987).
  * Wecker  _et al._ [2015] D. Wecker, M. B. Hastings, and M. Troyer, Progress towards practical quantum variational algorithms, Phys. Rev. A 92, 042303 (2015).
  * Wang and Thoss [2010] H. Wang and M. Thoss, From coherent motion to localization: II. dynamics of the spin-boson model with sub-Ohmic spectral density at zero temperature, Chem. Phys. 370, 78 (2010).
  * Zhang  _et al._ [2022] S.-X. Zhang, J. Allcock, Z.-Q. Wan, S. Liu, J. Sun, H. Yu, X.-H. Yang, J. Qiu, Z. Ye, Y.-Q. Chen, _et al._ , Tensorcircuit: a quantum software framework for the nisq era, arXiv preprint arXiv:2205.10091 (2022).
  * Ren  _et al._ [2021] J. Ren, W. Li, T. Jiang, Y. Wang, and Z. Shuai, The Renormalizer package. https://github.com/shuaigroup/renormalizer (2021).
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nat. Methods 17, 261 (2020).

##  Appendix A I. Derivation of variational basis state encoder

###  A.1 Time-independent equation

We start with the Lagrangian defined in the main text, i.e. Eq. 2. Taking the
derivative with respect to $C[l]_{mn}$ and setting it to 0 yields

$\braket{\phi}{n}_{l}\tensor{[}_{l}]{\bra{m}}{}\hat{\tilde{H}}^{\prime}[l]\ket{\phi}+\sum_{n^{\prime}}\lambda_{lnn^{\prime}}C[l]_{mn^{\prime}}=0\
.$ (9)

Multiply with $C[l]_{mn^{\prime\prime}}$

$\sum_{m}C[l]_{mn^{\prime\prime}}\braket{\phi}{n}_{l}\tensor{[}_{l}]{\bra{m}}{}\hat{\tilde{H}}^{\prime}[l]\ket{\phi}+\sum_{n^{\prime}}\lambda_{lnn^{\prime}}\sum_{m}C[l]_{mn^{\prime}}C[l]_{mn^{\prime\prime}}=0\
,$ (10)

and use the $C[l]$ orthonormal condition
$\sum_{m}C[l]_{mn}C[l]_{mn^{\prime}}=\delta_{nn^{\prime}}$ to get $\lambda$

$\lambda_{lnn^{\prime}}=-\sum_{m}C[l]_{mn^{\prime}}\braket{\phi}{n}_{l}\tensor{[}_{l}]{\bra{m}}{}\hat{\tilde{H}}^{\prime}[l]\ket{\phi}\
.$ (11)

Substitute $\lambda$ into Eq. 9 yields

$\braket{\phi}{n}_{l}\tensor{[}_{l}]{\bra{m}}{}\hat{\tilde{H}}^{\prime}[l]\ket{\phi}-\sum_{n^{\prime}m}C[l]_{mn^{\prime}}\braket{\phi}{n}_{l}\tensor{[}_{l}]{\bra{m}}{}\hat{\tilde{H}}^{\prime}[l]\ket{\phi}C[l]_{m^{\prime}n^{\prime}}=0\
.$ (12)

Using

$\hat{P}=\hat{B}[l]^{\dagger}B[l]=\sum_{mm^{\prime}}\sum_{n}\ket{m}_{l}C[l]_{mn}C[l]_{m^{\prime}n}\tensor{[}_{l}]{\bra{m^{\prime}}}{}\
,$ (13)

to simplify the notation of the second term

$\braket{\phi}{n}_{l}\tensor{[}_{l}]{\bra{m}}{}\hat{\tilde{H}}^{\prime}[l]\ket{\phi}-\braket{\phi}{n}_{l}\tensor{[}_{l}]{\bra{m}}{}\hat{P}\hat{\tilde{H}}^{\prime}[l]\ket{\phi}=0\
.$ (14)

Rearranging and rewriting in matrix form, we get the equation for $C[l]$

$(1-\hat{P}[l])\braket{\phi}{\hat{\tilde{H}}^{\prime}[l]}{\phi}=0\ .$ (15)

###  A.2 Quantum circuit measurement

In this section, we discuss the quantum circuit measurement required to solve
$C[l]$ from Eq. 3. The key quantity to be computed is matrix $G[l]_{mn}$,
defined as

$G[l]_{mn}=\braket{\phi}{n}_{l}\tensor{[}_{l}]{\bra{m}}{}\hat{\tilde{H}}^{\prime}[l]\ket{\phi}\
.$ (16)

Express $\hat{H}$ in sum-of-product form
$\hat{H}=\sum_{x}^{M}\prod_{k}\hat{h}[k]_{x}$, using notations in the main
text, and we get

$\displaystyle G[l]_{mn}$
$\displaystyle=\sum_{x}^{M}\braket{\phi}{n}_{l}\tensor{[}_{l}]{\bra{m}}{}\prod_{k\neq
l}\hat{B}[k]\prod_{k}\hat{h}[k]_{x}\prod_{k}\hat{B}[k]^{\dagger}\ket{\phi}$
(17)
$\displaystyle=\sum_{x}^{M}\braket{\phi}{n}_{l}\tensor{[}_{l}]{\bra{m}}{}\hat{h}[l]_{x}\hat{B}[l]^{\dagger}\prod_{k\neq
l}\hat{\tilde{h}}[k]_{x}\ket{\phi}\ .$

Here we assume $\prod_{k\neq l}\hat{\tilde{h}}[k]_{x}$ can be expressed by a
constant amount of Pauli strings. Represent $\hat{h}[l]_{x}$ in operator form

$\hat{h}[l]_{x}=\sum_{mm^{\prime}}h[l]_{xm^{\prime}m}\ket{m^{\prime}}_{l}\tensor{[}_{l}]{\bra{m}}{}\
.$ (18)

$G[l]_{mn}$ then becomes

$G[l]_{mn}=\sum_{x}^{M}\sum_{m^{\prime}n^{\prime}}h[l]_{xmm^{\prime}}C[l]_{m^{\prime}n^{\prime}}\braket{\phi}{n}_{l}\tensor{[}_{l}]{\bra{n^{\prime}}}{}\prod_{k\neq
l}\hat{\tilde{h}}[k]_{x}\ket{\phi}\ .$ (19)

Thus to evaluate $G[l]_{mn}$ it is sufficient to measure
$\braket{\phi}{n}_{l}\tensor{[}_{l}]{\bra{n^{\prime}}}{}\prod_{k\neq
l}\hat{\tilde{h}}[k]_{x}\ket{\phi}$.
$\ket{n}_{l}\tensor{[}_{l}]{\bra{n^{\prime}}}{}$ in general is not Hermitian,
and the real and imaginary parts can be measured by

$\displaystyle\Re{\braket{\phi}{n}_{l}\tensor{[}_{l}]{\bra{n^{\prime}}}{}\prod_{k\neq
l}\hat{\tilde{h}}[k]_{x}\ket{\phi}}$
$\displaystyle=\frac{1}{2}\bra{\phi}(\ket{n}_{l}\tensor{[}_{l}]{\bra{n^{\prime}}}{}+\ket{n^{\prime}}_{l}\tensor{[}_{l}]{\bra{n}}{})\prod_{k\neq
l}\hat{\tilde{h}}[k]_{x}\ket{\phi}\ ,$ (20)
$\displaystyle\Im{\braket{\phi}{n}_{l}\tensor{[}_{l}]{\bra{n^{\prime}}}{}\prod_{k\neq
l}\hat{\tilde{h}}[k]_{x}\ket{\phi}}$
$\displaystyle=\frac{1}{2}\bra{\phi}i(\ket{n^{\prime}}_{l}\tensor{[}_{l}]{\bra{n}}{}-\ket{n}_{l}\tensor{[}_{l}]{\bra{n^{\prime}}}{})\prod_{k\neq
l}\hat{\tilde{h}}[k]_{x}\ket{\phi}\ .$

To evaluate all matrix elements in $G[l]$, the total number of measurements
required scales as $\order{2^{N_{l}}M}$.

###  A.3 Time-dependent equation

In this section, we derive the time-dependent equation for $C[l]$. For time-
dependent problems, $C[l]$ in general is complex

$C[l]=D[l]-iE[l]\ ,$ (21)

where both $D[l]$ and $E[l]$ are real. The minus sign is for convenience
expressing $\hat{B}^{\dagger}\ket{\phi}$. From the definition we have

$\partialderivative{\ket{\psi}}{E[l]_{mn}}=i\partialderivative{\ket{\psi}}{D[l]_{mn}}\
.$ (22)

The starting point is the Lagrangian Eq. 4 defined in the main text. Taking
the derivative with respect to $\dot{\Theta}_{K}$ yields

$\displaystyle\partialderivative{\mathcal{L}}{\dot{\Theta}_{K}}$
$\displaystyle=\sum_{J}\partialderivative{\bra{\psi}}{\Theta_{J}}\partialderivative{\ket{\psi}}{\Theta_{K}}\dot{\Theta}_{J}+\sum_{J}\partialderivative{\bra{\psi}}{\Theta_{K}}\partialderivative{\ket{\psi}}{\Theta_{J}}\dot{\Theta}_{J}$
(23)
$\displaystyle\quad+i\partialderivative{\bra{\psi}}{\Theta_{K}}\hat{H}\ket{\psi}-i\bra{\psi}\hat{H}\partialderivative{\ket{\psi}}{\Theta_{K}}$
$\displaystyle\quad+\sum_{lnn^{\prime}}\lambda_{lnn^{\prime}}\Re{\sum_{m}C[l]^{*}_{mn}\partialderivative{\dot{C}[l]_{mn^{\prime}}}{\dot{\Theta}_{K}}}$
$\displaystyle\quad+\sum_{lnn^{\prime}}\gamma_{lnn^{\prime}}\Im{\sum_{m}C[l]^{*}_{mn}\partialderivative{\dot{C}[l]_{mn^{\prime}}}{\dot{\Theta}_{K}}}\
.$

We first consider the case of $\Theta_{K}=\theta_{k}$, and then

$\displaystyle\partialderivative{\mathcal{L}}{\dot{\theta}_{k}}$
$\displaystyle=\sum_{J}\partialderivative{\bra{\psi}}{\Theta_{J}}\partialderivative{\ket{\psi}}{\theta_{k}}\dot{\Theta}_{J}+\sum_{J}\partialderivative{\bra{\psi}}{\theta_{k}}\partialderivative{\ket{\psi}}{\Theta_{J}}\dot{\Theta}_{J}$
(24)
$\displaystyle\quad+i\partialderivative{\bra{\psi}}{\theta_{k}}\hat{H}\ket{\psi}-i\bra{\psi}\hat{H}\partialderivative{\ket{\psi}}{\theta_{k}}$
$\displaystyle=2\sum_{J}\Re{\partialderivative{\bra{\psi}}{\theta_{k}}\partialderivative{\ket{\psi}}{\Theta_{J}}}\dot{\Theta}_{J}-2\Im{\partialderivative{\bra{\psi}}{\theta_{k}}\hat{H}\ket{\psi}}\
,$

which means at the $\partialderivative{\mathcal{L}}{\dot{\theta}_{k}}=0$
minimum, we have

$\sum_{J}\Re{\partialderivative{\bra{\psi}}{\theta_{k}}\partialderivative{\ket{\psi}}{\Theta_{J}}}\dot{\Theta}_{J}=\Im{\partialderivative{\bra{\psi}}{\theta_{k}}\hat{H}\ket{\psi}}\
.$ (25)

Substitute $\Theta_{J}$ with $\theta_{k}$, $D[l]_{mn}$ and $E[l]_{mn}$

$\displaystyle\sum_{J}\Re{\partialderivative{\bra{\psi}}{\theta_{k}}\partialderivative{\ket{\psi}}{\Theta_{J}}}\dot{\Theta}_{J}$
$\displaystyle=\sum_{j}\Re{\partialderivative{\bra{\psi}}{\theta_{k}}\partialderivative{\ket{\psi}}{\theta_{j}}}\dot{\theta}_{j}$
(26)
$\displaystyle\quad+\sum_{lmn}\Re{\partialderivative{\bra{\psi}}{\theta_{k}}\partialderivative{\ket{\psi}}{D[l]_{mn}}}\dot{D}[l]_{mn}$
$\displaystyle\quad+\sum_{lmn}\Re{\partialderivative{\bra{\psi}}{\theta_{k}}\partialderivative{\ket{\psi}}{E[l]_{mn}}}\dot{E}[l]_{mn}\
.$

Using Eq. 22 the last two terms become

$\sum_{lmn}\Re{\partialderivative{\bra{\psi}}{\theta_{k}}\partialderivative{\ket{\psi}}{D[l]_{mn}}}\dot{D}[l]_{mn}+\sum_{lmn}\Re{\partialderivative{\bra{\psi}}{\theta_{k}}\partialderivative{\ket{\psi}}{E[l]_{mn}}}\dot{E}[l]_{mn}=\sum_{lmn}\Re{\partialderivative{\bra{\psi}}{\theta_{k}}\partialderivative{\ket{\psi}}{D[l]_{mn}}\dot{C}[l]_{mn}^{*}}\
,$ (27)

which is zero because

$\sum_{mn}\partialderivative{\bra{\psi}}{\theta_{k}}\partialderivative{\ket{\psi}}{D[l]_{mn}}\dot{C}[l]_{mn}^{*}=\sum_{mn}\partialderivative{\bra{\phi}}{\theta_{k}}\hat{B}[l]\ket{m}_{l}\tensor{[}_{l}]{\bra{n}}{}\dot{C}[l]_{mn}^{*}\ket{\phi}=0\
,$ (28)

where the constraint $\sum_{m}C[l]_{mn}\dot{C}[l]_{mn^{\prime}}^{*}=0$ is
used. Thus the simplified equation of motion reads

$\sum_{j}\Re{\partialderivative{\bra{\psi}}{\theta_{k}}\partialderivative{\ket{\psi}}{\theta_{j}}}\dot{\theta}_{j}=\Im{\partialderivative{\bra{\psi}}{\theta_{k}}\hat{H}\ket{\psi}}\
,$ (29)

or equivalently

$\sum_{j}\Re{\partialderivative{\bra{\phi}}{\theta_{k}}\partialderivative{\ket{\phi}}{\theta_{j}}}\dot{\theta}_{j}=\Im{\partialderivative{\bra{\phi}}{\theta_{k}}\hat{\tilde{H}}\ket{\phi}}\
.$ (30)

In short, the equation of motion for $\theta_{k}$ is the same as vanilla VQD
with encoded Hamiltonian $\hat{\tilde{H}}$ .

Next we consider the case of $\Theta_{K}=D[l]$ and $\Theta_{K}=E[l]$. After
some complex algebra, we have

$i\sum_{J}\partialderivative{\bra{\psi}}{D[l]_{mn}}\partialderivative{\ket{\psi}}{\Theta_{J}}\dot{\Theta}_{J}+i\frac{1}{2}\sum_{n^{\prime}}\lambda_{ln^{\prime}n}C[l]_{mn^{\prime}}^{*}-\frac{1}{2}\sum_{n^{\prime}}\gamma_{ln^{\prime}n}C[l]_{mn^{\prime}}^{*}=\partialderivative{\bra{\psi}}{D[l]_{mn}}\hat{H}\ket{\psi}\
.$ (31)

Similar to the case of $\Theta_{K}=\theta_{k}$, substitute $\Theta_{J}$ with
$\theta_{k}$, $D[l]_{mn}$ and $E[l]_{mn}$

$\displaystyle\sum_{J}\partialderivative{\bra{\psi}}{D[l]_{mn}}\partialderivative{\ket{\psi}}{\Theta_{J}}\dot{\Theta}_{J}$
$\displaystyle=\sum_{k}\partialderivative{\bra{\psi}}{D[l]_{mn}}\partialderivative{\ket{\psi}}{\theta_{k}}\dot{\theta}_{k}+\sum_{km^{\prime}n^{\prime}}\partialderivative{\bra{\psi}}{D[l]_{mn}}\partialderivative{\ket{\psi}}{D[k]_{m^{\prime}n^{\prime}}}\dot{C}[k]_{m^{\prime}n^{\prime}}^{*}$
(32)
$\displaystyle=\sum_{k}\partialderivative{\bra{\psi}}{D[l]_{mn}}\partialderivative{\ket{\psi}}{\theta_{k}}\dot{\theta}_{k}+\sum_{n^{\prime}}\partialderivative{\bra{\psi}}{D[l]_{mn}}\partialderivative{\ket{\psi}}{D[l]_{mn^{\prime}}}\dot{C}[l]_{mn^{\prime}}^{*}\
.$

Here the orthonormal condition is again used. Substitute the equation back
into Eq. 31.

$i\sum_{k}\partialderivative{\bra{\psi}}{D[l]_{mn}}\partialderivative{\ket{\psi}}{\theta_{k}}\dot{\theta}_{k}+i\sum_{n^{\prime}}\partialderivative{\bra{\psi}}{D[l]_{mn}}\partialderivative{\ket{\psi}}{D[l]_{mn^{\prime}}}\dot{C}[l]_{mn^{\prime}}^{*}+\frac{1}{2}\sum_{n^{\prime}}(i\lambda_{ln^{\prime}n}-\gamma_{ln^{\prime}n})C[l]_{mn^{\prime}}^{*}=\partialderivative{\bra{\psi}}{D[l]_{mn}}\hat{H}\ket{\psi}\
,$ (33)

Following the same strategy with the derivation of the time-independent
equation, multiply Eq. 33 with $C[l]_{mn}$

$i\sum_{k}\braket{\phi}{n}_{l}\tensor{[}_{l}]{\bra{n^{\prime}}}{}\partialderivative{\ket{\phi}}{\theta_{k}}\dot{\theta}_{k}+\frac{1}{2}(i\lambda_{ln^{\prime}n}-\gamma_{ln^{\prime}n})=\sum_{m}C[l]_{mn^{\prime}}\partialderivative{\bra{\psi}}{D[l]_{mn}}\hat{H}\ket{\psi}\
,$ (34)

where $\sum_{m}C[l]^{*}_{mn^{\prime}}C[l]_{mn}=\delta_{n^{\prime}n}$ and
$\sum_{m}\dot{C}[l]^{*}_{mn^{\prime}}C[l]_{mn}=0$ are used. Then, multiply
again with $C[l]_{mn}^{*}$

$i\sum_{k}\partialderivative{\bra{\psi}}{D[l]_{mn}}\partialderivative{\ket{\psi}}{\theta_{k}}\dot{\theta}_{k}+\frac{1}{2}\sum_{n^{\prime}}(i\lambda_{ln^{\prime}n}-\gamma_{ln^{\prime}n})C[l]_{mn^{\prime}}^{*}=\hat{P}[l]\partialderivative{\bra{\psi}}{D[l]_{mn}}\hat{H}\ket{\psi}\
.$ (35)

Use this equation to eliminate $\lambda$ and $\gamma$ in Eq. 33, we get the
equation of motion for $C[l]$

$i\sum_{n^{\prime}}\partialderivative{\bra{\psi}}{D[l]_{mn}}\partialderivative{\ket{\psi}}{D[l]_{mn^{\prime}}}\dot{C}[l]_{mn^{\prime}}^{*}=(1-\hat{P}[l])\partialderivative{\bra{\psi}}{D[l]_{mn}}\hat{H}\ket{\psi}\
,$ (36)

which can be simplified to Eq. 6. The measurement required for time evolution
is in the same order as the static VQE algorithm.

In the end, we note that imaginary time evolution might be another approach to
finding the ground state, in addition to the iterative method described in the
main text. Imaginary time evolution might also be a feasible approach to
determine $C[l]$ as an alternative to solving Eq. 3.

##  Appendix B II. Numerical noiseless simulations

All numerical quantum circuit simulation is performed using the TensorCircuit
[37] package without considering noise. Classical DMRG simulation is performed
using the Renormalizer package [38]. We use harmonic oscillator eigenstates
for phonon basis states. Using positional states might affect the performance
of traditional encodings because of the truncation, however, we expect
variational encoding to be insensitive to the choice of phonon basis states at
the $N\rightarrow\infty$ limit. We use Gray code for binary encoding as an
improvement to the standard approach [18]. For both ground state simulation
and dynamics simulation, $C[l]$ is initialized as $C[l]_{mn}=\delta_{mn}$.

For the VQE simulation of the Holstein model, the following ansatz is used

$\ket{\phi}=\prod_{l}^{L}\left\\{\prod_{\langle
j,k\rangle}e^{\theta_{ljk}(\hat{a}^{\dagger}_{j}\hat{a}_{k}-\hat{a}^{\dagger}_{k}\hat{a}_{j})}\prod_{j}e^{\theta_{lj}\hat{a}^{\dagger}_{j}\hat{a}_{j}(\hat{b}^{\dagger}_{j}-\hat{b}_{j})}\right\\}\ket{\phi_{0}}\
.$ (37)

where $L$ is the number of layers and $L=3$ is adopted. The advantage of Eq.
37 is enforcing real-valued wavefunction. The circuit parameters
$\vec{\theta}$ are optimized by the L-BFGS-G method implemented in SciPy
package [39]. The parameter gradient is calculated by auto-differentiation.
The initial values for the parameters are set to zero at the first round of
the macro-iteration. In subsequent macro-iterations, the previously optimized
parameters are used as the initial value for faster convergence. Eq. 3 is
solved by the DF-SANE method implemented in SciPy [39]. Since this is a non-
linear equation, we provide 3 initial guesses and adopt the one with the
lowest energy. The solved $C[l]$ sometimes does not satisfy the orthonormal
condition due to numerical imprecision and the orthonormal condition is
enforced by QR decomposition in each macro-iteration.

For the VQD simulation of the spin-boson model, the variational Hamiltonian
ansatz used is more complex than the VQE simulation. Because $C[l]$ is
complex, $\hat{B}[l]\hat{h}[l]_{x}\hat{B}[l]^{\dagger}$ spans the whole
Hermitian matrix space. Thus for $\hat{h}[l]_{x}$ the whole Pauli matrix set
$\\{X,Y,Z,I\\}^{\otimes N_{l}}$ is added to the ansatz. To obtain the
quantities required to calculate $\theta_{k}$, the Jacobian of the
wavefunction $\phi(\vec{\theta})$ is firstly calculated by auto-
differentiation, and then the r.h.s and l.h.s of Eq. 5 is calculated by matrix
multiplication. How to measure the quantities in realistic quantum circuits is
well described in the literature [16]. To calculate $\dot{C}[l]$ it is
necessary to take the inverse of $\rho[l]$ which is sometimes ill-conditioned.
We add $1\times 10^{-5}$ to the diagonal elements of $\rho[l]$ for
regularization. The time evolution of $\theta_{k}$ and $C[l]$ is carried out
using the RK45 method implemented in SciPy [39]. We observe that the gradient
of $\theta_{k}$ is usually much larger than $C[l]$. Thus it is possible to
evolve the two sets of parameters separately, which deserves further
investigation. For Trotterized time evolution, $N=16$ and a time step of 0.01
are used.

##  Appendix C III. Experiments on a superconducting quantum processor

###  C.1 Device parameters

The superconducting quantum processor, as shown in Fig. 3(a), is composed of
nine computational transmon qubits with each pair of neighboring qubits
mediated via a tunable coupler, forming a cross-shaped architecture. Each
computational qubit has an independent readout cavity for state measurement
and $XY$/$Z$ control lines for state operation. High-fidelity simultaneous
single-shot readout for all qubits are achieved with the help of the
multistage amplification with the Josephson parametric amplifier (JPA)
functioning as the first stage of the amplification. The fundamental device
parameters including qubit parameters and gate parameters are outlined in
Table. 2 and Table. 3, where the parasitic $ZZ$ interaction between qubits is
suppressed by the coupler.

Table 2: Single qubit gate parameters. $\omega_{r}$ is the resonant frequency of the readout cavity for each qubit. $\omega_{j,max}\,(j=1\sim 9)$ are the maximum resonant frequencies when qubits are biased at the sweet spot. $\omega_{j,idle}\,(j=1\sim 9)$ are the idle frequencies for implementing the single-qubit operations. $\alpha_{j}\,(j=1\sim 9)$ are the qubits’ anharmonicities. $T_{1}$, $T_{2,idle}$ and $T_{2E,idle}$ are the corresponding energy relaxation time, Ramsey dephasing time and echoed dephasing time for the qubits measured at the idle frequency. The readout fidelities are typically characterized by detecting each qubit in $\ket{g}$ ($\ket{e}$) when it is prepared in $\ket{g}$ ($\ket{e}$), labeled by $F_{0,j}$ and $F_{1,j}$. To mitigate the error coming from the readout infidelity, the outcomes are reconstructed with the calibration matrix through the Bayes’ rule. Single-qubit errors $e_{sq}$ are measured with randomized benchmarking (RB). |  $Q_{0}$ |  $Q_{1}$ |  $Q_{2}$ |  $Q_{3}$ |  $Q_{4}$ |  $Q_{5}$ |  $Q_{6}$ |  $Q_{7}$ |  $Q_{8}$
---|---|---|---|---|---|---|---|---|---
$\omega_{r}$ (GHz) |  $6.874$ |  $6.825$ |  $6.931$ |  $6.901$ |  $6.845$ |  $6.786$ |  $6.991$ |  $6.961$ |  $6.806$
$\omega_{j,max}$ (GHz) |  $4.003$ |  $4.215$ |  $4.479$ |  $4.689$ |  $4.470$ |  $4.479$ |  $4.657$ |  $4.512$ |  $4.362$
$\omega_{j,idle}$ (GHz) |  $3.988$ |  $4.187$ |  $4.464$ |  $4.668$ |  $4.404$ |  $4.359$ |  $4.641$ |  $4.498$ |  $4.223$
$\alpha_{j}/2\pi$ (MHz) |  $-260$ |  $-258$ |  $-255$ |  $-250$ |  $-254$ |  $-258$ |  $-253$ |  $-257$ |  $-264$
$T_{1}$ ($\mu$s) |  $35.3$ |  $31.6$ |  $29.5$ |  $27.7$ |  $33.9$ |  $34.3$ |  $33.3$ |  $22.1$ |  $31.8$
$T_{2,idle}$ ($\mu$s) |  $11.0$ |  $10.2$ |  $32.6$ |  $38.2$ |  $9.1$ |  $5.6$ |  $43.1$ |  $24.1$ |  $4.3$
$T_{2E,idle}$ ($\mu$s) |  $48.2$ |  $38.4$ |  $47.8$ |  $44.2$ |  $31.6$ |  $21.8$ |  $56.8$ |  $32.9$ |  $18.6$
$F_{0,j}$ (%) |  $96.9$ |  $97.4$ |  $98.6$ |  $98.9$ |  $98.7$ |  $98.4$ |  $96.3$ |  $97.2$ |  $94.1$
$F_{1,j}$ (%) |  $93.7$ |  $94.3$ |  $92.5$ |  $94.3$ |  $94.5$ |  $94.6$ |  $92.7$ |  $92.4$ |  $90.9$
$e_{sq}$ (%) |  $0.07$ |  $0.32$ |  $0.06$ |  $0.07$ |  $0.08$ |  $0.05$ |  $0.06$ |  $0.15$ |  $0.08$
Table 3: Two qubits gate parameters. $\omega_{c,idle}$ are the idle frequencies for each coupler where the $ZZ$ interaction between neighboring computational qubits are maximally suppressed. $\xi_{ZZ}$ is the residual $ZZ$ interaction between each qubit pairs. Two-qubit gates are implemented with the controlled-Z (CZ) and the corresponding gate errors $e_{tq,CZ}$ are characterized with RB. |  $Q_{0}-Q_{1}$ |  $Q_{0}-Q_{2}$ |  $Q_{0}-Q_{3}$ |  $Q_{0}-Q_{4}$ |  $Q_{1}-Q_{5}$ |  $Q_{2}-Q_{6}$ |  $Q_{3}-Q_{7}$ |  $Q_{4}-Q_{8}$
---|---|---|---|---|---|---|---|---
$\omega_{c,idle}$ (GHz) |  $5.020$ |  $5.445$ |  $5.570$ |  $5.335$ |  $5.325$ |  $5.595$ |  $5.695$ |  $5.355$
$|\xi_{ZZ}|$ (kHz) |  $18.0$ |  $10.0$ |  $5.0$ |  $8.0$ |  $2.0$ |  $3.0$ |  $5.0$ |  $2.0$
$e_{tq,CZ}$ (%) |  $1.57$ |  $2.22$ |  $1.99$ |  $2.47$ |  $0.91$ |  $1.04$ |  $1.2$ |  $0.96$

###  C.2 Experimental details

We use 3 qubits out of the 9-qubit computer for the 2-site Holstein model

$\hat{H}=-V(a^{\dagger}_{1}a_{2}+a^{\dagger}_{2}a_{1})+\omega
b^{\dagger}_{1}b_{1}+\omega b^{\dagger}_{2}b_{2}+g\omega
a^{\dagger}_{1}a_{1}(b^{\dagger}_{1}+b_{1})+g\omega
a^{\dagger}_{2}a_{2}(b^{\dagger}_{2}+b_{2})\ .$ (38)

The electronic degree of freedom is mapped to the second qubit. Thus,
$a^{\dagger}_{1}a_{1}$ is mapped to $\frac{1}{2}(1+Z_{1})$ and
$a^{\dagger}_{2}a_{2}$ is mapped to $\frac{1}{2}(1-Z_{1})$. The phonon modes
are mapped to the first and the third qubit. With binary encoding and
$N_{l}=1$, the Hamiltonian in the Pauli string form reads

$\hat{H}=-VX_{1}+\frac{1}{2}\omega(1-Z_{0})+\frac{1}{2}\omega(1-Z_{2})+\frac{1}{2}g\omega(1+Z_{1})X_{0}+\frac{1}{2}g\omega(1-Z_{1})X_{2}\
.$ (39)

For variational encoding, we assume $C[l]=C$. That is, the two modes share the
same variational encoder. This is a reasonable assumption for translational
invariant systems. Supposing $\hat{b}^{\dagger}\hat{b}$ and
$\hat{b}^{\dagger}+\hat{b}$ are mapped to the following form

$\displaystyle\hat{B}(\hat{b}^{\dagger}\hat{b})\hat{B}^{\dagger}$
$\displaystyle=c_{1i}I+c_{1x}X+c_{1z}Z$ (40)
$\displaystyle\hat{B}(\hat{b}^{\dagger}+\hat{b})\hat{B}^{\dagger}$
$\displaystyle=c_{2i}I+c_{2x}X+c_{2z}Z\ ,$

the encoded Hamiltonian is then

$\displaystyle\hat{H}$
$\displaystyle=-VX_{1}+\omega(c_{1i}I_{0}+c_{1x}X_{0}+c_{1z}Z_{0})+\omega(c_{1i}I_{2}+c_{1x}X_{2}+c_{1z}Z_{2})$
(41)
$\displaystyle\quad+\frac{1}{2}g\omega(1+Z_{1})(c_{2i}I_{0}+c_{2x}X_{0}+c_{2z}Z_{0})+\frac{1}{2}g\omega(1-Z_{1})(c_{2i}I_{2}+c_{2x}X_{2}+c_{2z}Z_{2})\
.$

We use the following ansatz for the parameterized quantum circuit

$\ket{\phi}=\prod_{j=1}^{2}e^{\theta_{j}a^{\dagger}_{j}a_{j}(b^{\dagger}_{j}-b_{j})}\frac{1}{\sqrt{2}}\left(\ket{000}+\ket{100}\right)\
,$ (42)

Because $C[1]=C[2]$, the parameter space can be further simplified by setting
$\theta_{1}=\theta_{2}$. With binary encoding, the ansatz transforms to

$\ket{\phi}=e^{i\theta Y_{2}}e^{-i\theta Z_{1}Y_{2}}e^{i\theta
Y_{0}}e^{i\theta Z_{1}Y_{0}}H_{1}\ket{0}\ .$ (43)

The ansatz is compiled into the following quantum circuit with 4 CNOT gates.

@C=1.0em @R=0.2em @!R
q_0 : & q_0 : R_Z(-π2) H R_Z(-θ) R_Z(-θ) H R_Z(π2)
q_1 : q_1 : H -1 -1 1 1
q_2 : q_2 : R_Z(-π2) H R_Z(θ) R_Z(-θ) H R_Z(π2)

Each energy term is measured by 8192 shots, and the uncertainty is obtained by
repeating the measurement 5 times and taking the standard deviation. For the
update of $C[l]$, 4096 shots are performed for each term. Local readout error
mitigation is applied for all results presented unless otherwise stated.

In Fig. 4 we plot the energy landscape $E(\theta)/V$ in VQE with binary
encoding. Both raw data and data with local readout error mitigation (EM) are
presented for the energy expectation from quantum hardware. The mitigated
landscape is in decent agreement with the perfect simulator. A minimum at
around $\theta=0.6$ is clearly visible. We note that the perfect simulator is
also based on the $N_{l}=1$ ansatz and $N$ is far smaller than what is
physically demanded. Thus the minimum presented by the perfect simulator can
not be recognized as the ground truth.

Figure 4:  VQE energy landscape for the 2-site Holstein model with binary
encoding. For the data from quantum hardware, both raw data and data with
readout error mitigation are presented. The error bar indicates the
measurement uncertainty.