Steady-state Peierls transition in nanotube quantum simulator

Quantum dots placed along a vibrating nanotube provide a quantum simulation platform that can directly address the electron-phonon interaction. This offers promising prospects for the search of new quantum materials and the study of strong correlation effects. As this platform is naturally operated by coupling the dots to an electronic reservoir, state preparation is straightforwardly achieved by driving into the steady state. Here we show that for intermediate electron-phonon coupling strength, the system with spin-polarized quantum dots undergoes a Peierls transition into an insulating regime which exhibits charge-density wave order in the steady state as a consequence of the competition between electronic Coulomb repulsive interactions and phonon-induced attractive interactions. The transport phenomena can be directly observed as fingerprints of electronic correlations. We also present powerful methods to numerically capture the physics of such an open electron-phonon system at large numbers of phonons. Our work paves the way to study and detect correlated electron-phonon physics in the nanotube quantum simulator with current experimentally accessible techniques.


INTRODUCTION
In the search for novel materials, the simulation of quantum matter is an extremely demanding computational task which is expected to profit substantially from the surge of quantum technologies. Quantum algorithms for programmable quantum computers offer the most flexible approaches 1-3 , but tailor-made quantum simulation devices are particularly well suited for large-scale simulations. For instance, tremendous efforts have been made to simulate graphene-like materials using synthetic quantum systems 4 , such as artificial semiconductor lattices 5 , cold atoms in optical lattices 6,7 , or photonic crystals 8 . However, in order to reach the complexity of real materials, quantum simulators need more controllable degrees of freedom. Specifically, apart from electrons, also phonons are fundamental ingredients of quantum materials, and the ubiquitous electron-phonon interaction is a central issue in condensed matter and material physics [9][10][11] . Many important and interesting phenomena are rooted in this interaction, dating back to the seminal work on Bardeen-Cooper-Schrieffer superconductivity 12 , or charge-density wave (CDW) order induced by the Peierls instability 13 .
An archetypal and simplified Hamiltonian with the electron-electron and electron-phonon interactions is the socalled Hubbard-Holstein model, which has been studied vastly [14][15][16][17][18][19][20][21] . Although many numerical and analytical approaches have been employed to solve this model, including quantum Monte Carlo 22 , density-matrix renormalization group (DMRG) 23 , variational ansatz 24,25 , dynamical meanfield theory (DMFT) 26 , and density-matrix embedding theory (DMET) 27 , it is still hard to study the strongly coupled electron-phonon systems for the large number of phonons. Approaching the electron-phonon models via quantum simulation and quantum computation techniques is challenging as well. With phonons being essentially absent in optical lattices, atomic quantum simulators rely on the explicit construction of dynamical lattices as recently proposed in refs. 28,29 or in the context of simulation of lattice gauge theory (cf. refs. 30,31 ), characterized by dynamical degrees of freedom on the bonds of the lattice. On the other hand, phonons occur as a natural ingredient in trapped ion quantum simulators. These systems are based on engineering the spin-phonon interactions, and hence appear to be quite natural candidates for the simulation of Holstein models 32,33 . However, mapping the trapped ions system onto electron-phonon problems via Jordan-Wigner transformation is not straightforward due to long-range couplings. Another approach could be the use of digital quantum computers. In this context, a scheme of mapping electron-phonon systems onto qubits has been proposed in ref. 34 .
On the other hand, there are also quantum simulation platforms which themselves consist of electrons and phonons and hence appear to be ideally suited for the study of quantum materials. One such platform that has recently been proposed in ref. 35 is quantum dots defined on a suspended carbon nanotube, where the phononic degrees of freedom are naturally provided by the flexural modes of the nanotube and the electrons localized in quantum dots interact with the mechanical modes electrostatically; see Fig. 1a for a sketch of the setup. Many parameters of the system can be tuned either at the fabrication stage or during the experiments (Fig. 1b). Especially, the short and controllable separation between the nanotube and gates allows one to reach a very large electronphonon coupling strength 36 . The nonlocal nature of phonon modes in this system also provides opportunities to explore physics beyond the Hubbard-Holstein model. Although the system is now limited in size, advances in nanofabrication make it promising to fabricate a carbon natotube with many quantum dots. On the other hand, even with the current small system, there are many interesting physical and technological phenomena associated, such as the phonon-induced pairing 35 and the nanomechanical qubit 37 .
As a simultaneous challenge and opportunity, one particular feature of this system is that the quantum dot setup can barely be viewed as an isolated quantum system, but rather  Figure 1. Setup and tuning range of parameters. a Schematic picture of the system. A suspended nanotube hosts four spin-polarized quantum dots. The back gates control the local chemical potential and hoppings between dots. The short separation between the nanotube and gates allows one to reach a very large electron-phonon coupling strength. The whole system couples to the left and right leads controlled by a bias voltage V bias . b The typical tuning range of different parameters. Here the fundamental frequency ω0 is divided by 2π, and other parameters are divided by 2π .
as an open quantum system which is coupled to a fermionic environment via leads. It would be important and experimentally relevant to be able to describe the transport and driving in this challenging system with many electronic and phononic degrees of freedom. While the open system properties of single and double dots coupled to a single phonon mode has been studied extensively [38][39][40][41][42][43][44][45][46][47][48][49][50][51][52][53][54] , the case of multiple dots has never been considered to our knowledge and may support interesting physics beyond single and double dots.
In this work, we study the correlated physics in a nanotube quantum simulator with four spin-polarized quantum dots coupled to the leads. We develop methods to attack the challenge of theoretically describing this electron-phonon open quantum system problem. Specifically, we generalize the shift method, developed in ref. 35 to facilitate the computational treatment of equilibrium systems with large phonon numbers, to the case of open quantum systems. This treatment then allows us to study transport through the device, a great opportunity to detect important system properties. In particular, the electric current provides an immediate smoking gun of a striking behavior of the device: Upon increasing the electron-phonon coupling strength, a Peierls instability turns the system into an insulator with alternating CDW pattern of empty and occupied quantum dots as a consequence of the competition between electronic Coulomb repulsive interactions and phonon-induced attractive interactions, which is beyond the standard Hubbard-Holstein model. Although such a CDW order is also present without coupling the system to leads, observing it as a steady state property in open quantum system provides a feasible route for the preparation as well as the detection of this intriguing phenomenon. Especially, despite the fact that the current always vanishes at low enough bias due to the finite charge gap in the considered finite-size system, the different behavior of critical bias voltage for nonvanishing current as we increase the electronphonon coupling strength directly reflects the different nature of insulating steady state, identifying the Peierls transition.

RESULTS
System. The central component of our system, depicted in Fig. 1a, is a suspended carbon nanotube which hosts up to four electrons on equally spaced quantum dots spin-polarized by a large magnetic field, interacting among themselves through long-range Coulomb interactions and with the nanotube's flexural modes. Accordingly, the system Hamiltonian consists of a Hubbard-like electronic part a phononic part and the coupling between electrons and phonons Here d i is the spinless annihilation operator for electron on the i-th quantum dot, and n i = d † i d i is the electron number operator. Note that a single quantum dot cannot be occupied by two electrons due to the Pauli exclusion principle. The bosonic operator b α is the annihilation operator of the phonon mode α with frequency ω α = αω 0 being an integer multiple of the fundamental frequency ω 0 , and is the reduced Planck constant. The electrons couple to the phonons through the capacitance between the quantum dots and gates, which is modulated by the displacement of the nanotube. For the coupling strength g i,α , we assume the validity of the guitarstring model 35 , in which we expand the capacitance charge energy at small displacement in terms of different oscillating modes and integrate the coupling strength between the charge and mode displacement b α + b † α over the dot extension (assumed to be 1/4 of the nanotube length), yielding g i,α = g 0 (8/π)α −3/2 sin[(2i − 1)απ/8] sin(απ/8) with an overall coupling strength g 0 that can be tuned electrostatically 35 . We also assume that the electronic hopping t 0 is uniform and between nearest neighbors only. A local chemical potential ε is included to control the particle number. The parameters U j−i with j > i describe screened Coulomb interactions. An overview of realistic parameter ranges is illustrated in Fig. 1b.
For a quantum simulation of the transport behavior of the electron-phonon system, we need to take into account a fermionic environment which couples to the nanotube via two leads at the left (L) and right (R) end. The Hamiltonian of the environment reads for = L, R, where the electrons in leads are denoted as c and are labeled by the momentum k. The coupling between the system and leads is given by where d denotes a system electron in the dot at the left or right end. In the following, we consider the wide-band limit for the leads with constant density of states ν. We also assume the tunneling amplitudes to be energy-independent and symmetric, i.e., t k = t and t L = t R . Then the coupling strength between the system and leads can be captured by the tunneling rate Γ = 2πν|t L/R | 2 . Quantum master equation. The Hamiltonian H SE consists of four different tunneling processes between the system and leads: (i, ii) an electron entering the system from the left or right lead; (iii, iv) an electron leaving the system to the left or right lead. It is clear that these processes occur only in a particular region of the system and depend on the occupation of the environment level ε k = E n,in − E n+1,in+1 , where E n,in is the system energy for the i n -th eigenstate of H S with total electron number n. For systems with the tunneling rate Γ much smaller than the temperature k B T (k B is the Boltzmann constant), we can use the phenomenological position and energy resolving Lindblad master equation 55 to describe the reduced density matrix ρ of the system coupled to leads Here the Lindblad operators are given by for processes (i, ii) and Electron number Ne Figure 2. Role of electron-phonon coupling in the equilibrium system at atomic limit t0 = 0. a Electron number Ne in equilibrium as a function of local chemical potential ε/U1 and electron-phonon coupling strength g0/ √ ω0U1. b Energy of different electron configurations with fixed Ne = 2. Here we have used the single-mode approximation, which is verified by comparing the results obtained from single mode (solid lines) and multi (e.g., 100 considered here) modes (dots). The single-and multi-mode scenarios are found to agree with each other qualitatively well. The transition points between patterns for processes (iii, iv), where we have T ( ) n,in;n−1,in−1 = E n,in |d † |E n−1,in−1 after rewritting the tunneling Hamiltonian H SE into the eigenstate basis of H S . The Fermi-Dirac distribution f (x) = [e (x−µ )/kBT + 1] −1 with temperature T and chemical potential µ captures the occupation of environment levels in the leads. It is clear that the rate of adding an electron to the system should be proportional to the occupation in the leads, while it is hard to remove an electron from the system if the corresponding state has already been occupied in the leads. Here we use the bias voltage V bias to control the chemical potentials in leads and have µ L/R = ±eV bias /2 (e is the elementary electric charge), which is an important tuning knob in the open quantum system.
The Lindblad master equation describes the evolution into the steady state ρ ss , which is obtained in the infinite time limit, i.e. for t → ∞, or by diagonalizing the Liouvillian superoperator L. A challenge to solve this problem is how to treat the phononic degrees of freedom for its infinite-dimensional nature. Even after truncating the phononic Hilbert space, for realistic system parameters the phonon number remains beyond our numerical capability (limited to tens of phonons). For this, we have developed different methods which facilitate the description by shifting the phononic vacuum into a state with finite phonon number and only taking into account a small number of necessary (tilded) phononic states that are coupled to the electronic states effectively. In the following, we would first focus on the physical results and shall describe the details of shift method in the "Methods" section for interested readers.
We would like to mention that if the tunneling rate Γ is much smaller than the temperature k B T and the energy split- Electron number Ne eV bias /2πh 80 GHz 100 GHz 120 GHz 0 0.5 1 tings between states with the same number of electrons, a convenient simplification to the Lindblad master equation can be achieved by ignoring coherences 56,57 , and we only need to consider the diagonal terms of the density matrix, ρ n,in ≡ ρ n,in;n,in . This approximation leads to the Pauli master equation Although this approximation is not always valid, we can use the phonon number from Pauli master equation as a good starting point for the shift method; see Methods.
Role of electron-phonon coupling in equilibrium. The main application of our quantum simulation platform is to study the correlated electron-phonon physics. Hence we focus on the role played by the tunable electron-phonon coupling. Before turning our attention to the phonon-induced transport behavior in the steady state, we will first briefly show the role of electron-phonon coupling in the equilibrium system. The role of electron-phonon coupling in equilibrium can be most easily observed in the case of t 0 → 0. In this atomic limit, the electron-phonon problem can be solved analytically via the Lang-Firsov transformation U = e − i,α gi,αni(b † α −bα)/ ωα , which yields the effective system Hamiltoniañ Here the electron-phonon coupling is replaced by an effective phonon-induced long-range electron-electron interaction. The first role played by electron-phonon coupling is that the phonons mediate an attractive interaction, thus lowering the energy of states with more electrons; see Fig. 2a for the electron number N e = i n i in equilibrium. Second, even at a fixed electron number, the electron-phonon coupling has a strong effect on the distribution of electrons along the quantum dots, as shown in Fig. 2b for N  as a compromise between the two interactions for intermediate values of electron-phonon coupling.
A few remarks are in order. First, it is clear from the Lang-Firsov transformation that the phonon-induced electronelectron interaction is proportional to g 2 0 / ω 0 . For a large fundamental frequency, we also need a large electron-phonon coupling to obtain the same phonon-induced electron-electron interaction strength, but the underlying physics remains the same in different parameter regime. For this and to avoid a too large phononic Hilbert space for open quantum systems, we set ω 0 /2π = 3 GHz for most of our calculations, which is slightly larger than the typical value in experiments but is still reasonable; see Fig. 1b. Second, since the phonon-induced interaction decays in the scaling of O(α −4 ) for higher phonon modes, the higher modes would only affect the physics quantitatively but not qualitatively; see Fig. 2b for a comparison between the single-mode approximation and multi-mode scenario. Hence it is a good approximation to consider only the lowest phonon mode without losing the essential physics. In the following studies, we will always take this approximation and omit the subscript index of the lowest phonon mode (i.e., b = b 1 ) unless stated otherwise. Third, we note that the three equilibrium regimes discovered above are not limited to the atomic limit but still survive for a small hopping coefficient t 0 ; see Fig. 3 and Fig. 5. The finite hopping problem is solved via the equilibrium shift method proposed in ref. 35 ; see also Methods for a review. We fix the hopping coefficient as t 0 /2π = 5 GHz in the following studies of open quantum systems. For this value, the above discussed regimes are still clearly present.
Phonon-induced transport behavior. With the above insight for the equilibrium properties of the system, we now characterize the steady state. The results are presented in Fig. 3.
Since the local chemical potential ε also influences the electron population (cf. Fig. 2a) and the interesting physics mainly occurs in the regime with N e = 2, here we choose a local chemical potential such that the equilibrium state with two electrons starts from g 0 = 0. In Fig. 3a, we plot the lowest energies of different electron number sectors for the chosen local chemical potential ε/2π = −30 GHz. The equilibrium ground state changes from the two-electron sector to the four-electron sector upon increasing the electron-phonon coupling. The range of electron-phonon coupling strength with two electrons is also large for this local chemical potential, where different two-electron regimes can happen, although the range for four electrons is even larger and goes to infinity. We note that since the region for equilibrium state •••• is quite small compared to other two-electron regimes and is very close to the four-electron state for the local chemical potentials with two-electron states starting from g 0 = 0 in equilibrium (see Fig. 2a which characterizes the net rate of electron number flow through the system and vanishes in these regions; see Fig. 3c. While the insulating nature of four-electron state is a trivial consequence of Pauli blocking, it is a very intriguing behavior for N e = 2. In fact, this insulating regime coincides with a CDW pattern in the electronic configuration. To detect this pattern, we introduce the structure factor which takes the value C = 4 in the perfect CDW order. Indeed, the N e = 2 plateau in Fig. 3b coincides with CDW order; see Fig. 3d. As indicated by the black line, this regime also exhibits CDW order in equilibrium. Therefore, this pattern is not induced by the tunneling to leads but indeed by the competition between repulsive Coulomb interaction and phonon-induced attractive electron-electron interaction. We also introduce the order parameter  Fig. 3e, together with the equilibrium value of two-electron states. For small electron-phonon coupling, while the equilibrium system exhibits the pattern ••••, the order parameter O is highly reduced in the nonequilibrium steady state due to the enhanced mixing between patterns •••• and ••••/•••• and with the one-electron eigenstates. Hence it is hard to distinguish the weakly coupled regime from the CDW order via the order parameter O, unlike the CDW structure factor C. We note that the current is always small at low enough bias due to the finite charge gap in the considered finite-size system. To further identify the nature of insulating steady state at certain coupling strength, we define the critical bias voltage V bias,c , at which the corresponding current is 10 −3 GHz, a very small value for the used tunneling rate. Since the critical bias voltage is directly related to the charge gap and the changes of energy spectra as we increase the electron-phonon coupling are given by ( i g i,1 n i ) 2 / ω 0 in the atomic limit [cf. Eq. (11)], we would expect that the critical bias voltage is proportional to g 2 0 and depends on the electron distributions of the coupled states by leads. Indeed, this is also confirmed for a finite but small hopping t 0 in Fig. 4, where the critical bias voltages can be captured by a series of linear fittings. The different slope manifests the distinct nature of the corresponding insulating steady states.
Initially, the two-electron state in equilibrium mainly has the pattern •••• but is also perturbed by the CDW order due to the finite t 0 , which will be coupled to the one-electron states made of components ••••/•••• and ••••/•••• by leads (the coupling to three-electron states are highly suppressed by the energy gap; cf. Fig. 3a). Since the one-electron state dominated by pattern ••••/•••• has lower energy, the above two-electron state is mainly coupled to this state for small bias voltage, determining the charge gap. The slight g 0dependence of energies of these states gives a small slope for the critical bias voltage. As we increase the electron-phonon coupling, the Peierls instability from pattern •••• into the CDW order suddenly enhances the reduction of energy for two-electron state (cf. Fig. 2b), which is still mainly coupled to the one-electron state dominated by pattern ••••/••••, leading to the quick increase in critical bias voltage and a large slope. This feature provides an immediate smoking gun for the Peierls transition. As the coupling further increases, the energy of three-electron state would be lower than that of the one-electron state. The two-electron state now mainly couples to the three-electron state dominated by pattern ••••/•••• (the weight of pattern ••••/•••• will increase for larger g 0 due to the reduced energy), and the charge gap as well as the critical bias voltage start to decrease, although the insulating steady state still exhibit CDW order (Fig. 3d). The further reduced charge gap may even lead to a small increase in the current for a fixed large bias voltage; cf. Fig. 3c.
Finally, we consider the phonon number N p = Tr[ρ ss b † b] in the steady state using Pauli master equation; see Fig. 3f and discussions in Methods. An abrupt increase of N p is observed when the electron-phonon coupling becomes too strong to support the CDW structure, with little dependence on the applied bias voltage. Interestingly, for relatively large bias voltage, the phonon number also takes a large value for small electron-phonon coupling. This counter intuitive phenomenon (cf. Fig. 5a for phonon number in equilibrium) indeed is very different from the large phonon number at strong electron-phonon coupling. Actually the phonon number N p has two contributions, one coming from the polaronic displacement of the oscillator ∆x = Tr[ρ ss b] + Tr[ρ ss b † ], thus N p ≈ (∆x/2) 2 , and another from the distribution of phonons. The last part is related to the variance of phonon operator Fig. 3f), and measures the number of excited phonons in the system. The results indicate that for small coupling the system heats up. This is due to the current passing through the device, as one sees that in coincidence with the degeneracy point (g 0 / √ ω 0 U 1 ≈ 0.77) the variance has also a maximum and that the variance increases with the voltage in these regions. For strong electron-phonon coupling, the variance of phonon operator is quite small. This suggests that the large phonon number in this regime is only due to the polaronic displacement of the oscillator that is proportional to the electromechanical coupling constant g 0 . The system is actually in its ground state, since the current is blocked and cannot transfer heat to the oscillator.
We can also understand the large phonon number at small coupling from the perspective of the Pauli master equation. In the limit of g 0 → 0, the electronic and phononic degrees of freedom are indeed decoupled, which means that the leadinduced outgoing and incoming tunneling rates for an eigenstate of H S in Eq. (10) are uniform with respect to the phonon number in this limit. At small g 0 , the weak electron-phonon coupling only has limited influence on the eigenstates, hence the tunneling rates are still very close and do not have significant difference in the whole range of phonon numbers. This allows for the occupation of excited states with large phonon number at small electron-phonon coupling. Obviously, a larger bias voltage would allow more eigenstates to be coupled, leading to a larger phonon number.

DISCUSSION
In summary, we have investigated the steady state of a quantum simulation platform for correlated electron-phonon physics, where four spin-polarized quantum dots are equally placed along a suspended carbon nanotube and are coupled to a fermionic environment provided by two leads at the left and right end. Different regimes in the steady state as well as in equilibrium as a function of electron-phonon coupling strength have been explored. Particularly, for intermediate electron-phonon coupling, the Peierls instability drives the system into an insulating regime with CDW order as a consequence of the competition between electronic Coulomb repulsive interactions and phonon-induced attractive interactions. The different behavior of critical bias voltage for nonvanishing current further provides an immediate smoking gun for the different nature of insulating steady state, identifying the Peierls transition. To treat the electron-phonon open systems with large phonon numbers, we also developed a general-  Figure 5. Equilibrium two-electron states at finite hopping t0. a The phonon number Np increases as the electron-phonon coupling becomes stronger and is more smooth for a larger t0. b, c The order parameter O and CDW structure factor C get reduced, but the three equilibrium regimes discovered in the atomic limit still survive under small but finite t0. Here we set U1/2π = 200 GHz, U2/2π = 20 GHz, U3/2π = 2 GHz, and ω0/2π = 3 GHz. The local chemical potential ε is ignored for the two-electron case, and the bias voltage is set to zero for these plots.
ized shift method in the "Methods" section, which highly reduces the phononic degrees of freedom required to describe the problem. Our work shall largely stimulate the experimental advances in this field. Particularly, the proposed setup can be accessed via current techniques in nanofabrication. The high tunability of different parameters allows us to directly observe and detect the predicted phonon-induced correlated physics and transport behavior. Although our work focuses on the spinpolarized quantum dots, it would be interesting to further take into account the spin and/or valley degrees of freedom for the open electron-phonon systems in the future studies, which shall bring in new kind of orders and nonequilibrium behavior.

METHODS
Shift method for equilibrium systems at finite hopping. In Eq. (11), we have utilized the Lang-Firsov transformation to transform the equilibrium electron-phonon Hamiltonian in the atomic limit into a purely electronic problem. However, the electrons and phonons cannot be decoupled at finite hopping, since the hopping coefficient transforms as t 0 → t 0 e −(gi,1−gi+1,1)(b † −b)/ ω0 and will gain a phonon-dependent phase. One common approach is to truncate the infinitedimensional phononic Hilbert space to a finite-dimensional Hilbert space and then perform the numerical diagonalization. However, this method is limited to small phonon numbers. To overcome this problem, a shift method has been developed in ref. 35 .
The basic idea of this method is to shift the phononic vacuum to a state with finite number of phonons by the transformation b →b + S and only take into account a small number of effective phononic states that are coupled to the electrons (e.g., we consider a Hilbert space with maximal 30 tilded phonons). The shift parameter S will be updated iteratively according to the replacement S ← −[ (b † + S * )(b + S) ] 1/2 until convergence, where the expectation · is performed on the ground state of the shifted Hamiltonian H S [d i ,b, S]. In the practical calculation, we can utilize the phonon number  Figure 6. Comparison between different methods for solving the steady state. Here we choose a large fundamental frequency ω0/2π = 80 GHz such that the phonon number in the system is small. We truncate the phononic Hilbert space to a finite-dimensional Hilbert space with maximal 30 (or 100) phonons for the direct method (or Pauli master equation approach). On the other hand, the shift methods with and without updating the shift parameter are performed within a Hilbert space with maximal 5 tilded phonons. Here we set eV bias /2π = 100 GHz, and other parameters are the same as in Fig. 3. N p = ( i g i,1 n i / ω 0 ) 2 in the atomic limit to obtain a good initial guess S 0 = − N p for the shift parameter.
As an application of this method and a complement to Fig. 2b, we show the two-electron equilibrium states at finite hopping in Fig. 5, with a particular interest in the fate of CDW order at large values of hopping coefficient. First, unlike the phonon number in steady state (cf. Fig. 3f), the equilibrium phonon number N p = (b † + S * )(b + S) always increases as the electron-phonon coupling increases; see Fig. 5a. There is also a clear discontinuity for small t 0 at the transition point from CDW order to the •••• pattern due to the first order nature of this transition, which is smoothed by large hopping coefficient. For the electronic properties, the three equilibrium regimes discovered in the atomic limit survive under small but finite t 0 and distinguish from each other clearly, as indicated in the order parameter O (Fig. 5b) and CDW structure factor C (Fig. 5c), although the corresponding values are reduced. On the other hand, at large hopping coefficient, the boundary between the pattern •••• and CDW configuration is strongly blurred. This is attributed to the fact that the energy gap between these two states is small and changes slowly as we increase the electron-phonon coupling strength; cf. Fig. 2b. Therefore, these two states will strongly affect each other for a wide range of g 0 , blurring the distinctions. Generalized shift methods for steady states. Now we generalize the shift method to open electron-phonon systems. We note that the steady state ρ ss can be considered as the ground state of an effective Hamiltonian L † L. As in the shift method for equilibrium ground states, we first make the transformation b →b + S to shift the phononic vacuum to a state with phonon number The above procedures are repeated until convergence of the shift parameter. It is also a good initial guess for the shift parameter to use the phonon number from Pauli master equation.
We would like to mention that the dimensionality of the shifted phononic Hilbert space in general depends on the chemical potential and temperature in the leads and the phonon frequency, which determine how many phonons will be coupled by the environment. For a large bias voltage and temperature compared with the phonon frequency, the number of tilded phonons may be still large. In this case, although the effective phononic degrees of freedom have been highly reduced compared to the bare phonons, it is still not easy to solve the steady state. Therefore, the above method is limited to relatively small bias voltage and temperature.
On the other hand, we note that the fermionic environment directly couples to the electronic degrees of freedom in the system. To solve the steady state for a large tilded phononic Hilbert space, we further make an approximation for the shift method, i.e., we solve the phonon number from the Pauli master equation and use it as a shift parameter for the shift method to develop the possible coherence between electronic degrees of freedom using the time evolution or iterative method until the electronic properties converge, which is much faster than the convergence of phonon number N p = Tr{ρ ss [S](b † + S * )(b + S)}. In this method, the shift parameter will not be updated, and the phonon number from Pauli master equation is assumed to be close to the real one. For the regime with small electron-phonon coupling strength or ignorable coherence, this approximation should work well, since the corresponding coherence between phononic degrees of freedom is weak. We note that even for the regime beyond this limitation, the electronic properties obtained from this method are still more reasonable than the results obtained from the Pauli master equation, as part of the coherence has been captured.
In Fig. 6, we benchmark our generalized shift methods by observing the steady state phonon number N p , electron number N e , and current I. Here we choose a very large fundamental frequency ω 0 such that the phonon number in the system is small and it is still possible to solve the steady state directly (i.e., without using the shift methods). Compared with the direct method by solving the equation L[ρ ss ] = 0, the results obtained from Pauli master equation matches qualitatively. However, the difference may be large near g 0 / √ ω 0 U 1 ≈ 0.77, which corresponds to the equilibrium transition point from two-electron states to the four-electron states (cf. Fig. 3a). Here the energy levels are close to each other, and the coherence cannot be ignored. Due to the lack of coherence, the current from Pauli master equation is in general larger than the true value even if the phonon number and electron number are nearly the same. We note that the phonon number from Pauli master equation only differs from the real one significantly in the region where the electron-phonon coupling is strong and the coherence is relevant.
On the other hand, the results from shift method matches quantitatively with the exact results, and the effective dimension of phononic Hilbert space has been highly reduced compared with the direct method (e.g., from 30 to 5 in our calculation). Hence the shift method works quite well. We also plot the electron number and current obtained from shift method without updating the shift parameter, where the shift parameter is provided by the phonon number from Pauli master equation and we evolve the shifted Lindblad master equation until the electronic properties converge. The matching with the full shift method as well as the exact results suggests that this method is also a good approximation.

DATA AVAILABILITY
The data supporting the results in this work are available from the corresponding author upon reasonable request.

CODE AVAILABILITY
The code supporting the results in this work are available from the corresponding author upon reasonable request.