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

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 $\mathcal{O}(N)$ and $\mathcal{O}(\log{N})$ qubits as well as $\mathcal{O}(N)$ and $\mathcal{O}(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)$ for systems obeying the area law of entanglement entropy. 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 one or two 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.

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 {|m } and then the second step is encoding {|m } into quantum computational basis {|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 {|m } → {|n } [22,23]. The first is unary encoding [24,25], in which each |m is encoded to |00 . . . 1 m . . . 00 , and the total number of qubits required scales as O(N ). The second is binary encoding, in which each |m is encoded to * zgshuai@tsinghua.edu.cn † shengyzhang@tencent.com i | m 2 i mod 2 represented by O(log N ) qubits [16][17][18]. In terms of two-qubit gates required to simulate quantum operators such asb † ±b andb †b , unary encoding scales as O(N ) and binary encoding scales as O(N log N ) [23]. 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 electronphonon systems. In this work, we propose a new basis encoding scheme called variational encoding. Variational encoding maps linear combinations of |m that are most entangled to the simulated system into the computational basis, i.e. m C mn |m → |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 {|n } can be made irrelevant to the size of {|m }. In other words, the number of qubits required scales as O(1). Consequently, the scaling for the number of gates is also O(1). The premise of the scaling reduction is the area law of entanglement entropy.
Variational encoding is best suited to work in combination with variational quantum algorithms such as variational quantum eigensolver (VQE) [26,27] and variational quantum dynamics (VQD) [28,29]. Besides, the variational encoding is also compatible with Trotterized time evolution and quantum phase estimation (QPE) [12,30,31]. 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.

II. VARIATIONAL BASIS STATE ENCODER
In this section, we present a more rigorous formulation of the variational basis state encoder. To encode each phonon mode l, encoded by N l qubits, we define the variational basis state encoderB[l] as followŝ with orthonormal constraint or equivalently In this paper, we use |m to represent phonon states and |n to represent qubit states. Eq. 1 can be rewritten aŝ and it is clear thatB performs m C mn |m → |n , The original Hamiltonian in |m basisĤ can then be encoded to |n basis using the following expressioñ For both static and dynamic cases, encoder coefficients C are determined by variational principle. In the remainder of the section, we will derive the equation for C. We use atomic units throughout the paper.

A. Time-independent equation
Suppose the quantum circuit is parameterized by |φ = k e iθ kRk |φ 0 , and then the ground state Lagrangian with multipliers λ lnn is Taking the derivative with respect to θ k immediately leads to traditional VQE with encoded HamiltonianĤ ∂ φ|Ĥ|φ ∂θ k = 0 (7) Taking the derivative with respect to C[l] mn and setting it to 0 yields whereĤ [l] is the the half-encoded Hamiltoniañ Multiply Eq. 8 with C[l] mn and use the C[l] orthonormal condition Eq. 3 to get λ Define projector (11) Substitute λ (Eq. 10) into Eq. 8 then yields (1 −P [l]) φ|Ĥ [l]|φ = 0 . Here To summarize, circuit parameters θ k are solved by VQE according to Eq. 7, and variational parameters C[l] are determined by solving Eq. 13 classically. Because Eq. 7 contains C[l] and Eq. 13 contains θ k , θ k and C[l] are solved iteratively until convergence. In the following, this iteration is termed macro-iteration to avoid confusion with VQE iteration.

B. Quantum circuit measurement
In this section, we discuss the quantum circuit measurement required to solve C[l] from Eq. 13. The key quantity to be computed is matrix G[l] mn , defined as Suppose the Hamiltonian can be written as a sum of direct productsĤ where M is the total number of terms in the Hamiltonian andĥ [k] x acts on the kth degree of freedom. Similar to the encoded Hamiltonian, the encoded local operator is denoted asĥ For electron degree of freedom a dummy encoderB[k] = I is used for notational simplicity. G[l] mn can then be written as G[l] mn then becomes Thus to evaluate G[l] mn it is sufficient to measure J[l] xnn . |n l n | l in general is not Hermitian, but the real and imaginary parts can be measured separately with (|n l n | l + |n l n| l ) and i(|n l n| l − |n l n | l ). Assuming the number of measurement shots for each Pauli string is N shots , the number of measurements to determine J[l] is thus O 2 N l M N shots , which is polynomial to the system size and does not increase with N . After J[l] is measured, evaluating G[J] and the left-hand side of Eq. 13 scales as O 2 N l N 2 M by matrix multiplication on a classical computer. Considering the measurement of a parameterized quantum circuit takes much longer time than a float-point number operation on classical computers, the classical workload is negligible compared to the additional measurements for reasonable values of N shots and N , such as N shots = 4096 and N = 64. Thus the reduction in quantum resources is not achieved by increasing classical resources [32].
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 O 2 N l M 2 N shots . 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 ≤ 2 is sufficient to produce excellent results.

C. Time-dependent equation
For time-dependent problems, it is convenient to define and use Θ K denote both θ k and C[l]. The Lagrangian with multipliers λ lnn and γ lnn is then The constraints ensure that C[l] mn remains orthonormal during time evolution. Taking the derivative with respect toΘ K yields The subsequent derivation involves a more intricate process similar to that of Eq. 13. Elaborate details are documented in Appendix A. Here, we provide an outline of the crucial steps. First, consider the case where Θ K = θ k and we find that the equation of motion for θ k is the same as vanilla VQD with encoded HamiltonianĤ Next, consider the case where Θ K = C[l], which ultimately leads to the following equation forĊ[l] where ρ[l] nn = Tr{ φ|n l n | l φ } is the reduced density matrix for the N l qubits of |φ . Eq. 13 represents aĊ[l] = 0 stationary point during real and imaginary time evolution. The measurement cost is the same as the ground state algorithm.
While we have relied on parameterized quantum circuits (PQC) for our derivation thus far, it is worth noting that incorporating the variational encoder into Trotterized time evolution and QPE is a straightforward extension. The VQD step described by Eq. 23 can be naturally replaced by a Suzuki-Trotter time evolution step e −iĤτ ≈ M x e −iĥxτ 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. 24, 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 us-ingĤ.

D. Variational basis state encoder as an ansatz
It is instructive to observe that if the variational basis encoder is viewed as a wavefunction ansatz |ψ , then the algorithm proposed in this work can be viewed as a generalization for the local basis optimization method for DMRG [33,34], or a special case of the recently proposed quantum-classical hybrid tensor network [35]. Thus,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 [36]. Consequently, | ψ|Ψ | 2 , the fidelity between the approximated encoded state and the target state has a lower bound of 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.

A. Numerical simulation on a noiseless simulator
The variational basis state encoder is first tested for VQE simulation of the one-dimensional Holstein model [37,38] whereâ andb are annihilation operators for electron and phonon respectively, V is the hopping coefficient, i, j denotes nearest neighbour pairs with periodic boundary condition, ω is the vibration frequency and g is dimensionless coupling constant. In the following, we assume V = ω = 1 and adjust g for different coupling strengths. We consider a three-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. We devise the following ansatz where L is the number of layers and L = 3 is adopted. More details of the simulation are included in the Appendix B. 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 → ∞ 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 three 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 [39], described by the Hamiltonian where is the eigenfrequency and ∆ is the tunneling rate. The coupling term has a similar form with Eq. 25 and is more commonly written as j c jσzxj . For systems in the condensed phase the coupling is usually characterized by the spectral density function J (ω) = π 2 j c 2 j ωj δ(ω − ω k ). In the following we assume = 0 and ∆ = 1. We first use VQD for the simulation and discuss Trotterized time evolution at last. The variational Hamiltonian ansatz [40] with three layers is used if not otherwise specified.
The performance of variational encoding and binary encoding is first compared based on a one-mode spinboson model at the strong coupling (ω = 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 onemode 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 two-mode model with ω j = 1 2 , 1 and g j = 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 ω and g are determined by discretizing a sub-Ohmic spectral density J (ω) = π 2 αω s ω 1−s c e −ω/ωc following the prescription in the literature [41]. The parameters are s = 1 4 , ω c = 4 and α = 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 ω = g = 1 and N l = 1. The measurement and the evolution of C[l] are performed at each Trotter step.

B. Verification on a superconducting quantum processor
In this section, we verify the accuracy and efficiency of the variational encoder approach on a superconducting quantum processor [42,43]. We consider the ground state problem of a 2-site Holstein model described by Eq. 25 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 Supplemental Material. 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.

IV. CONCLUSION
We proposed a variational basis state encoder to encode phonon basis states into quantum computational states for efficient quantum simulation of electronphonon systems. The proposed variational encoding approach requires only O(1) qubits and O(1) quantum gates for systems obeying the area law of entanglement entropy, which is significantly better than traditional encoding schemes and 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 O(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. The 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 spinboson model indicate that variational encoding is more accurate and resource-efficient than traditional encoding methods. In particular, using one or two qubits to represent each phonon mode is sufficient for accurate simulation even at the strong coupling regime where N = 32 phonon basis states are involved. The approach could also be extended to other quantum simulation problems involving an infinite or large local Hilbert space.
In this section, we derive the time-dependent equation for C[l]. For time-dependent problems, C[l] in general is complex where both D[l] and E[l] are real. The minus sign is for convenience expressingB † |φ . From the definition we have The starting point is Eq. 22 We first consider the case of Θ K = θ k , and then which means at the ∂L ∂θ k = 0 minimum, we have Using Eq. A2 the last two terms become which is zero because where the constraint m C[l] mnĊ [l] * mn = 0 is used. Thus the simplified equation of motion reads or equivalently In short, the equation of motion for θ k is the same as vanilla VQD with encoded HamiltonianĤ .
Here the orthonormal condition is again used. Substitute the equation back into Eq. A10.
Use this equation to eliminate λ and γ in Eq. A12, we get the equation of motion for C[l] which can be simplified to 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. 13.

Appendix B: Numerical simulation details
All numerical quantum circuit simulation is performed using the TensorCircuit [44] package and the Ten-CirChem [45] package without considering noise. Classical DMRG simulation is performed using the Renormalizer package [46]. 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 → ∞ limit. We use Gray code for binary encoding as an improvement to the standard approach [22]. For both ground state simulation and dynamics simulation, For the VQE simulation of the Holstein model, the circuit parameters θ are optimized by the L-BFGS-G method implemented in SciPy package [47]. 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. 13 is solved by the DF-SANE method implemented in SciPy [47]. Since this is a non-linear equation, we provide three 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,B[l]ĥ[l] xB [l] † spans the whole Hermitian matrix space. Thus forĥ[l] x the whole Pauli matrix set {X, Y, Z, I} ⊗N l is added to the ansatz. To obtain the quantities required to calculate θ k , the Jacobian of the wavefunction φ( θ) is firstly calculated by auto-differentiation, and then the r.h.s and l.h.s of Eq. 5 in the main text are calculated by matrix multiplication. How to measure the quantities in realistic quantum circuits is well described in the literature [19]. To calculateĊ[l] it is necessary to take the inverse of ρ[l] which is sometimes ill-conditioned. We add 1 × 10 −5 to the diagonal elements of ρ[l] for regularization. The time evolution of θ k and C[l] is carried out using the RK45 method implemented in SciPy [47]. We observe that the gradient of θ 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. II. Single qubit gate parameters. ωr is the resonant frequency of the readout cavity for each qubit. ωj,max (j = 1 ∼ 9) are the maximum resonant frequencies when qubits are biased at the sweet spot. ω j,idle (j = 1 ∼ 9) are the idle frequencies for implementing the single-qubit operations. αj (j = 1 ∼ 9) are the qubits' anharmonicities. T1, 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 |g (|e ) when it is prepared in |g (|e ), labeled by F0,j and F1,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 esq are measured with randomized benchmarking (RB).
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.b †b = m m |m m| is then encoded tô It is then possible to express the encoded operator aŝ where c 1i = (F 00 + F 11 )/2 , c 1x = F 01 = F 10 , c 1z = (F 00 − F 11 )/2 .
The ansatz is compiled into the following quantum circuit with 4 CNOT gates.
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(θ)/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 θ = 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.