Polaron dynamics of Bloch–Zener oscillations in an extended Holstein model

Recent developments in qubit engineering make circuit quantum electrodynamics devices promising candidates for the study of Bloch oscillations (BOs) and Landau–Zener (LZ) transitions. In this work, a hybrid circuit chain with alternating site energies under external electric fields is employed to study Bloch–Zener oscillations (BZOs), i.e. coherent superpositions of BOs and LZ transitions. We couple each of the tunable qubits in the chain to dispersionless optical phonons and build an extended Holstein polaron model with the purpose of investigating vibronic effects in the BZOs. We employ an extension of the Davydov ansatz in combination with the Dirac–Frenkel time-dependent variational principle to simulate the dynamics of the qubit chain under the influence of high-frequency quantum harmonic oscillators. Band gaps emerge due to energy differences in site energies at alternating qubit sites, and are shown to play key roles in tuning band structures and time periodic reconstructions of the wave patterns. In the absence of qubit–phonon interactions, the qubits undergo either standard BZOs or breathing modes, depending on whether the initial wave packet is formed by a broad or narrow Gaussian wave packet, respectively. The BZOs can get localized in space if the band gaps are sufficiently large. In the presence of qubit–phonon coupling, the periodic behavior of BZOs can be washed out and undergo dynamic localization. The influence of an ohmic bath on the dynamics of BZOs is investigated by means of a Markovian master equation approach. Finally, we calculate the von Neumann entropy as a measure of the entanglement between qubits and phonons.


Introduction
Bloch oscillations (BOs) and Landau-Zener (LZ) transitions are two fundamental phenomena in quantum mechanics associated with coherent transport of quantum particles and waves in periodic media under external driving forces [1,2]. In periodic potentials, there are continuous Bloch band structures and delocalized Bloch eigenfunctions that are uniformly distributed over the whole system. When an electric field is introduced in the form of a direct current (DC), the qubit chain exhibits discrete localized eigenstates that receive the name of Wannier-Stark ladders [3,4]. In the absence of scattering and dephasing effects, BOs refer to the periodic motion of charged carriers in reciprocal and real space [5]. When the variation of a linear potential due to a superimposed field is comparable with the smallest gaps to an adjacent band, the LZ transitions emerge and carriers may successively tunnel to higher-order bands [6,7]. Coherent superpositions of BOs and LZ transitions receive the name of Bloch-Zener oscillations (BZOs) [8][9][10]or Landau-Zener-Bloch oscillations [10], and have found applications in various fields, including quantum optics [11], solid-state physics [12], or atomic and molecular physics [13]. In fact, BZOs have already been directly observed in a variety of systems, such as light beams in two-dimensional, square photonic lattices [9], light waves in optical, binary superlattices [11], or rubidium Bose-Einstein condensates [14], to name a few.
In recent decades, a growing interest has been paid to tunable platforms that are based on circuit quantum electrodynamics (QED) because of their high-precision control, scalability, and long decoherence time [15]. For these reasons, the application of circuit QED setups as quantum simulators may provide relevant insights into the fundamental problems of quantum mechanics . These circuits can be implemented on various physical platforms, ranging from an individual artificial atom, to one-and two-dimensional lattices [34,35]. Versatile platforms can meet requirements of universal quantum computations and can be designed for specific quantum simulation tasks [36]. More importantly, rapid progresses in circuit QED devices have opened the door to the investigation of the interplay between BOs and LZ transitions, which is the focus of this work. Recently, several QED experiments have been dedicated to the analysis of BOs. In 2017, Ramasesh et al realized Bloch-oscillating quantum walks in a circuit QED protocol which is formed by a transmon qubit and a superconducting cavity mode [37]. Bahmani et al theoretically proposed a one-dimensional chain of coupled resonant circuits to probe the BOs and Wannier-Stark ladders [38]. Guo et al demonstrated the presence of BOs on a circuit processor of 5 superconducting transmon qubits [39].
Circuit QED experiments have also been performed to investigate LZ transitions. Izmalkov et al observed in 2004 the LZ transitions in a macroscopic circuit composed by an Al three-junction qubit and a Nb resonant tank [40]. Johansson et al measured the LZ transitions on an individual flux qubit within a superconducting chip of device qubits [41]. Hänggi and co-workers studied the the LZ transition probability at long times in circuit QED [42]. Neilinger et al demonstrated the LZ transitions in two experiments utilizing circuit QED arrangements with two types of superconducting qubits: the flux qubit on the basis of conventional Josephson junctions, and the phase-slip qubit based on niobium nitride nanowires [43]. Coherent transport in the BZOs can therefore be investigated by employing a circuit QED configuration.
Nevertheless, the quantum states of the qubits in the circuit QED platform will inevitably be affected by surrounding environments. The interaction between qubits and their environments has been probed in a variety of circuit QED setups [44][45][46]. Dissipative effects and the interaction between BOs, LZ transitions and phonons can be studied in hybrid circuit QED setups that incorporate for example mechanical resonators or microwave cavities [47,48]. O' Connell et al realized single-phonon control in a microwave-frequency mechanical resonator coupled to a flux qubit [49]. Riedinger et al demonstrated the correlated photon-phonon pairs based on a Duan-Lukin-Cirac-Zoller protocol [50]. Manenti et al measured average coherent phonon populations in a circuit QED architecture containing a superconducting qubit coupled to a surface acoustic wave cavity via piezoelectric effects [51]. In spite of these experimental efforts, the influence of quantized oscillators in the dynamics of BOs and LZ transitions is yet not fully understood.
A variational method that employs a multi-mode representation of a Davydov ansatz called multi-D2 has been developed for the simulation of open quantum systems with applications to fermions and bosons [52,53]. Here, the accuracy of the variational ansatz is increased by adding more modes until the value of the error, a measure based on a geometrical intepretation of the variational manifold, is maintained below a desired accuracy. In that respect we say that the method is numerically accurate within the theoretical model of the system. The method has been applied to a multitude of systems, including a Holstein model of polaron dynamics in the presence of external fields [54][55][56], phase transitions in the spin-boson model [57] or singlet fission model [58], a dissipative LZ model that is based on QED protocols [41,59] and other investigations in quantum optics [60,61].
In the present work, we employ the aforementioned multi-D2 ansatz to perform numerically accurate simulations of BZOs in a one-dimensional model of a qubit chain with alternating site energies under the influence of a constant electric field and proton. Our focus here will be the interaction between the qubits and the high frequency, optical phonons that have an intramolecular origin. The remainder of the paper is organized as follows. In section 2, we show the model Hamiltonian and explain the properties of the multi-D2 ansatz as the trial wave function. In section 3.1, the energy structures and transitions of a purely electronic qubit chain is studied. In section 3.2, the effects of different initial conditions and band gaps on the BZOs dynamics of a free excitation are investigated. In section 3.3, the effects of undamped, high-frequency phonon modes on the BZOs dynamics are studied. Finally, the influence of dephasing Ohmic bath on the dynamics of BZOs is examined in section 3.4. Conclusions are drawn in section 4.

Model
An schematic of the one-dimensional quantum circuit used in this work is shown in figure 1. The qubits are gap-tunable and the site energies of the model can thus be controlled [62]. The signs in qubit j − 1 and in qubit j illustrate that the magnetic fluxes at consecutive sites point at opposite directions, which also affects the sign of their site energies. The magnetic flux Φ e (t) represents an external driving force acting on the qubit array. The resultant biperiodic potential serves as an ideal platform for the investigation of BZOs [8]. The Hamiltonian of the one-dimensional qubit chain schematically shown in figure 1 is given bŷ where J is the nearest-neighbours hopping integral, δ is the difference in the site energies at alternating sites, F denotes an external electric field,â n (â † n ) is the annihilation (creation) operator of the electronic excitation at site n, and the distance between neighbouring sites is set to unity for simplicity. Periodic boundary conditions can be introduced by bending the linear array into a ring structure, leading to discontinuities in the applied electric potential at the boundaries due to the constant nature of the external field (DC-current) [4]. Thanks to the usage of a gauge transformed vector potential, A(t) = −Ft, this discontinuity problem can be avoided [63], and equation (1) can be recast as followŝ Taking into consideration the high-frequency optical phonons at each site, the Hamiltonian of the complete electronic-vibrational (vibronic) system can be written aŝ These undamped high-frequency phonon modes have been shown to strongly influence the excitation transfer and optical properties of organic semiconductors, photosynthetic molecular complexes [64], and qubit-resonator systems [65]. The phonon HamiltonianĤ ph and qubit-phonon coupling Hamiltonian H qu-ph terms take the following formĤ where g is the diagonal qubit-phonon coupling strength andb n (b † n ) denotes the annihilation (creation) operator of the phonons with site n. Next we transform the bosonic operators from real space to momentum space by expanding them in a Fourier serieŝ and rewrite the termsĤ qu ,Ĥ qu-ph of the Hamiltonian as follows, in which ω q is the phonon frequency of the mode with momentum q. Here we use dispersionless, Einstein phonons with ω q = ω 0 , and set the value of ω 0 to unity as the reference energy scale of this work. It is noted that the phonon related terms (equations (4) and (5)) in this work are the same with those in the standard Holstein model, which can be written asĤ Holstein = −J nâ † n (â n+1 +â n−1 ) + ω 0 nb † nb n − gω 0 nâ † nân b n +b † n and describes non-interacting electrons in a one-dimensional lattice coupled to dispersionless phonons [66,67]. The electronic terms in our work describe one electronic excitation in a double-band lattice under the external constant field, and is an extension of the electronic parts of the standard Holstein model. Therefore, our model can be treated as a type of extended Holstein model in terms of the Holstein language. In the Holstein system, a bare electron is dressed by a cloud of bosonic excitations, and form a quasiparticle called polaron. Polaron dynamics in the Holstein model have been continuingly studied using the following variational method [52][53][54][55][56].

Multi-D2 ansatz
The variational ansatz that we employ in this work is the multi-D2 Davydov ansatz [52,55,56], a superposition of electronic excitations and coherent states that takes the following form where M is the multiplicity of the ansatz and measures the number of different coherent states per site and i labels the ith coherent state. ψ in ∈ C denotes the time-dependent excitation amplitude of the nth qubit and λ iq ∈ C is the phonon displacement of the qth mode. The time evolution of the time-dependent variational where the Lagrangian L is to be found in appendix A. The integration of equations of motion for the variational variables u i can result in numerical instabilities [57]. An apoptosis procedure has been implemented to circumvent the singularity problem of the numerical solver that makes the multi-D2 ansatz a rather stable tool in the study of polaron dynamics [68]. Using this ansatz, we can provide numerically stable and accurate dynamics for both electronic and bosonic degrees of freedom. The validity of the variational approach has been extensively investigated and confirmed in previous works [41,52,55]. The accuracy of the multi-D2 ansatz is again certified in this work for sufficiently high multiplicity M. More details on the accuracy test can be found in appendix B. Observables of interest include the time evolution of the excitation probability of the qubit state, P qu (t, n), and the time evolution of phonon displacement, X ph (t, n), as follows In the presence of the external electric field, the current of the composite system can be described as Furthermore, the expected value c(t) of the excitation's position is utilized to characterize the centroid of the wave packet, and the standard deviation of the excitation wave packet σ(t) is used to describe how far it moves away from the origin, The shape of the initial standard deviation σ 0 has a significant impact on the mean value c(t) of the excitations' position [55]. Based on the multi-D2 trial state, the qubit energy, phonon energy, qubit-phonon coupling energy, and total energy of the system can be respectively written by In order to better characterize the effects of phonons on the qubit system, the bipartite quantum entanglement is investigated. In particular, we calculate the von Neumann entropy as a measure of the amount of entanglement between the qubit system and the optical phonons [69]. The time evolution of the von Neumann entropy can be written as where the reduced single-excitation density matrix,ρ qu (t) = Tr phρ (t), is the partial trace of the full density matrix,ρ (t) = D M 2 (t) D M 2 (t) , by tracing over the environmental phonon degrees of freedom at time t.

Spectral properties of BZOs phenomenon
We start by discussing the spectral properties, band structure and miniband transitions of the Hamiltonian of equation (1). We consider the Hamiltonian (1) in the Bloch representation to study energy properties in the reciprocal space. In this representation, the electronic ground state of the periodic system is usually described in terms of extended Bloch orbitals, simultaneous eigenstates of the periodic Hamiltonian and of the direct lattice translations. In the absence of an external electric field (F = 0), the eigenenergies are given by where μ = 0, 1 is the miniband index, k is the wave vector index, and η = sgn(δ). The electronic band structure in a reduced Brillouin zone of the qubit Hamiltonian of equation (1)  In the presence of an external field (F = 0), the transition matrix elements ofĤ qu between minibands μ = 0, 1 with momentum k, k can be expressed as φ 1,k Ĥ qu φ 0,k = T k δ π k − k , where |φ μ,k represents a Bloch miniband μ of momentum k and the delta function δ π k − k ensures that only direct interband transitions take place. The integrated transition matrix element T k is Transitions occur mainly around the edge of the reduced Brillouin zone and, due to the large energy difference between minibands, almost no transitions take place at the center of the zone. This effect has been verified in experiments [9].
As an example, we take the value J = 0.2 for the electronic coupling strength between qubits in an external electric field F = 0.1 to study the transition amplitudes T k , whose absolute value is shown in figure 2(b). The case with the smallest energy difference between consecutive sites (δ = 0.2287) exhibits the largest values, while the transition amplitudes vanish for a gap-free case. In the Wannier representation, the solution to the eigenvalue problemĤ qu |ψ ν,n = E ν,n |ψ ν,n in a two-band model of the qubit chain gives rise to Wannier-Stark states |ψ ν,n that are indexed by what is known as the mini-ladder index ν = 0, 1 and are 2n periodic in the lattice. The two energy ladders are equally spaced and have an energy offset of 2E 0 (δ), The spectrum of E ν,n (δ) is shown in figure 2(c) and an enlarged plot at zero energy is also shown in figure 2(d). The orange curve in figure 2(d) corresponds to the results of a numerical simulation with  [70]. When δ = 0, the excitation undergoes BZOs as a combined phenomenon involving both BOs and LZ transitions [16]. The time periodicity of BZOs is dominated by the LZ transition, and the time period is usually expressed as T BZ = sT LZ = rT re , with r, s ∈ N and T re is known as the reconstruction period, and reference [16] describes these time periodicities in great detail. The energy spacing between the Wannier-Stark ladders is reflected in the period of the LZ transition T LZ = 2π /(F − 2E 0 (δ)). The reconstruction period is T re = π /F = T B /2 because the size of the reduced Brillouin zone for δ = 0 is half of that for δ = 0. When E 0 (δ) = 0, we can obtain a relation of T BZ = 2T re = T B .

BZOs dynamics of a free excitation
Next we investigate the dynamics of BZOs under a wide range of parameters. We remark how the sensitive dependence of electronic dynamics on the energy difference between consecutive sites δ can be exploited to manipulate the evolution of the wave packet. For simplicity, we assume that the initial wave packet in the qubit chain is shaped as a Gaussian function of the form in which C represents a normalization coefficient and σ 0 is the initial standard deviation of the Gaussian wave function [55]. We shall assume thereon that the initial state is centered at the site n = 0, with zero momentum. We start by studying the purely electronic case in the absence of phonons. Figure 3 shows the excitation dynamics for various values of δ and a broad initial state of spatial width σ 0 = 2. When the qubit chain is energetically homogeneous (δ = 0) we obtain the distinctive pattern of BOs with oscillation of the center-of-mass and essentially constant wave shape, as shown in figure 3(a). When an energetic bias is introduced δ = 0, interband tunneling induces coherent superpositions of BOs and LZ transitions. As the LZ transitions take place mainly when the wave packet is close to the boundary of the energy band, transitions zones can be characterized at times t = T BZ /4 and 3T BZ /4 in one T BZ period. Before t = T BZ /4, the wave packet mainly propagates in a lower miniband and the travelling direction is orthogonal to the band curve. At T BZ /4, part of the wave tunnels to a higher miniband and travels along different spatial paths. Around 3T BZ /4, the LZ transitions occur again and the two branches are combined to reconstruct the original wave shape. The BZOs pattern is most obvious in the case of δ = 0.2287 where the interband transitions dominate, as shown in figures 2(b) and 3(b). As depicted in figures 3(c) and (d) for δ = 0.5772 and 0.8301, the wave packet starts to get localized for higher values of δ, although the characteristic structure of BZOs can still be inferred. Further localization of the wave packet ensues as δ/J ∼ 1, as shown in more detail in appendix C In addition to the initially broad Gaussian wave packet, we also consider a narrow initial state in figure 4 with standard deviation σ 0 = 0.1. Here, the wave packet is initially localized at site n = 0. When the chain is energetically homogeneous (δ = 0), the characteristic motion of a breathing mode is shown in figure 4(a). The wave packets widen and shrink in space around the origin in a symmetric manner with precise time periodicity. The dynamics of the qubit chain in the presence of an energetic bias δ = 0 is shown in figures 4(b)-(d). Interestingly, the spatial symmetry is destroyed, while the temporal periodicity is preserved. The localization effects are also clearly shown in the plots with large values of δ. In addition, we show in figure 5 the time evolution of various observables of interest, such as the mean value of the position of the excitation c(t), the standard deviation of the initial state σ(t), and the currents j(t). For the broad, Gaussian initial condition, the values of c(t) in the case of BOs is larger than those of BZOs as shown in figure 5(a). The case of δ = 0.2287 has the largest displacement σ(t) due to its stronger tunneling, as shown in figure 5(c). We show in figure 5(e) that the amplitude of the current j(t) gets reduced and oscillations become more prominent as the energetic bias δ is increased. When the particle is initially localized at the origin as in the case of a narrow Gaussian wave packet, the center of the wave packet oscillates for cases with nonzero δ, as shown in figure 5(b), the amplitude of σ(t) decreases as δ increases, and the current j(t) is no longer equal to zero for cases with δ = 0, as shown in figure 5(e).

Effects of the optical phonon
In this subsection, we proceed to investigate the BZOs dynamics of the qubit chain in the presence of dispersionless optical phonons, using the linear, vibronic coupling model described by the extended Holstein Hamiltonian of equation (3). In all following simulations, we take δ = 0.2287 as the energy difference between consecutive sites as we did in the purely electronic case for an initially broad or narrow wave packet in figures 3(b) and 4(b), respectively. This value of δ is chosen so that T BZ = T B . In addition we fix the number of sites to N = 30, with electronic coupling (J = 0.2) and the strength of the external field (F = 0.1).
In analogy with the previous section, we first study the time evolution in figure 6 of a broad Gaussian packet with standard deviation σ 0 = 2. The characteristic patterns of BZOs are washed out as the qubit-phonon coupling is increased from g = 0.1 (figure 6(a)) to g = 1.4 ( figure 6(d)). When the qubit-phonon coupling is sufficiently weak, the contour profile of the BZOs can still be recognized, as shown in figure 6(a). For intermediate coupling, the structure of BZOs start to smear out after time 3t BZ /4,  as displayed in figures 6(b) and (c). In the strong coupling regime, the center of the wave packet becomes localized, as visualized in figure 6(d). The time evolution of the phonon displacement X ph (t, n) at each site is shown in figure 7 for the initially broad Gaussian state. Initially the phonons are in the vacuum state. As the electronic excitation travels along the chain, phonons are emitted across the chain following the dynamics of the electronic wave packet. Subsequently, the emitted phonons are scattered by the electronic excitation, leading to the aforementioned vague patterns of BZOs. The BZOs frequency ω BZ is equal to the BOs frequency ω B since the value δ = 0.2287 was chosen to satisfy T BZ = T B . Since the external field  F = 0.1, ω BZ (ω BZ = F) is ten times smaller than the phonon's frequency ω 0 = 1, ten peaks can be found in one period of the BZOs.
Complementary to the evolution of an initially broad wave packet, we show in figure 8 the time evolution of an initially narrow wave packet in the presence of phonons for various values of the qubit-phonon coupling strength. As the qubit-phonon coupling strength is increased in figures 8(a)-(c), the spatial symmetry of the left and right branches around the origin is broken, and irregular excitation motion ensues. The localization effect by the strong coupling on the evolution of the excitation populations can be clearly seen in figure 8(d). We notice how the range of the colorbars increases with g in figure 9, indicating a higher phonon emission rate for an increasing qubit-phonon coupling strength. Additionally, figure 10 illustrates the effects of phonons on the expected value of the particle's position c(t), the standard deviation σ(t) of the wave packet and the currents j(t). Our results show that the excitation of the qubit states can be manipulated by the external field and qubit-phonon coupling.  In order to better understand the interaction between the qubit chain and the surrounding phonons, we have analyzed the time evolution of various energy variables, such as the kinetic energy of the excitation E qu (t), the total energy of the high-frequency phonons E ph (t) and the interaction energy between the qubits and the phonons E qu-ph (t). The total energy E tot (t) of the whole system is conserved if no external field is added (F = 0). As shown in figure 11, the total energy varies with time when an external field is applied to the qubit chain. In the presence of an external field F = 0.1, the qubit chain undergoes either BZOs or a breathing mode when the qubit-phonon interaction is turned off (g = 0), as shown in figures 11(a) and (b) for the initially broad and narrow cases, respectively. When g > 0 the energy of the high-frequency phonons quickly grows from zero and starts to oscillate with high-frequency, as evidenced by the high value of the ratio ω 0 /J = 5, as shown in figures 11(c) and (d) for the broad and narrow case, respectively.  Next we proceed to evaluate the quantum correlations between qubits and phonons by analyzing the time evolution of the von Neumann entropy, given by equation (17). First, we vary the qubit-phonon coupling strength and calculate the von Neumann entropy S vN in figure 12 for different values of the qubit-phonon coupling and band gaps. It can be found that large band gaps δ restrict the entanglement between the qubit and phonons. For an initially broad Gaussian wave packet and fixed δ, increasing the qubit-phonon coupling strength g the large amplitude of the oscillatory dynamics of the von Neumann entanglement entropy are dramatically reduced as g ∼ 1. This can be interpreted as a signature of the dynamic localization of the particle. In contrast, when the excitation starts from a narrow Gaussian wave packet and g J, the entanglement entropy increase steadily S vN (t) ∝ t κ for some κ < 1. In all the simulations so far we have taken J = 0.2ω 0 as the electronic coupling strength between neighbouring qubits (we remind the reader that ω 0 = 1 is taken as energy unit throughout this work). That is to say, the studied cases above focus on the anti-adiabatic regime (ω 0 J). We now vary the adiabaticity ratio ω 0 /J for fixed qubit-phonon coupling strength g = 0.2 to clarify the different energy scales involved in the dynamics. As shown in figures 13(a) and (d), the effects of phonons on the qubits are accentuated in the anti-adiabatic limit when compared with the adiabatic case. In the adiabatic regime ω 0 /J = 1, energy resonance between the electronic and vibrational subsystems are established. In the adiabatic regime (ω 0 /J = 1), we show the dynamics of S vN (t) for weak (figures 12(b) and (e)) and strong (figures 12(c) and (f)) qubit-phonon coupling strength g. The effect of qubit-phonon interactions in the entanglement entropy is overall less pronounced as in the anti-adiabatic case as a result of the higher electronic coupling strength. In comparison with the anti-adiabatic case of figure 12, the entanglement entropy reaches higher values in the adiabatic case for strong values of the qubit-phonon coupling strength, as shown in figures 13(c) and (f).

Dissipative dynamics of BZOs in an ohmic environment
Recent developments in circuit QED have shown that qubits can be experimentally coupled to an ohmic phonon bath [41]. We proceed to study the effects of an Ohmic vibrational environment on the dynamics of BZOs. We introduce the Hamiltonian of the phonon reservoirĤ R (equation (4)) with a linear vibronic interaction termĤ qu-R (equation (5) in whichâ n (â † n ) is the annihilation (creation) operator of the excitation at site n,b nξ (b † nξ ) is the annihilation (creation) operator of the ξth phonon mode with frequency ω nξ , and g nξ denotes the coupling strength of the ξth phonon mode at site n. We notice that the interaction term of equation (23) has the structurê H qu-R = nÂ n ⊗B n withÂ n =â † nâ n andB n = ξ g n b + nξ +b nξ operating on the electronic and vibrational subsystems, respectively. After bringing the number of oscillators in the environment to infinity, the coupling between the qubits and the bosonic reservoir can be characterized by the corresponding spectral density function, J R (ω) = π ξ g 2 nξ δ ω − ω ξ = 2αωe −ω/ω c , where α is the dimensionless coupling strength and ω c is the cut-off frequency of the ohmic environment. In order to simulate the dissipative dynamics of BZOs, we employ the following Markovian master equation [71], where the reduced density matrix,ρ qu (t) = Tr Rρ (t), is the partial trace of the density matrix of the total system over the reservoir's degrees of freedom, and the dissipator D takes the form andΠ( ) = k | k k | is a projector onto the subspace of qubit eigenstates with energy k . The parameter γ nm (ω) sets the dissipation rate at each frequency and is given by in which the reservoir correlation functions B † n (s)B m (0) = Tr ph ρ ph e iH ph sB n e −iH ph sB m in the interaction picture. We assume that each qubit is subject to independent (uncorrelated) environments, i.e. γ nm (ω) = γ mn (ω) = δ n,m C(ω), equation (25) simplifies significantly and where n(ω) = exp ω/k B T − 1 −1 represents the mean phonon number of a phonon mode with a frequency ω at temperature T, and J(ω) is given by the Ohmic spectral density. The polaron dynamics of BZOs in the presence of an Ohmic environment are shown in figure 14, where the excitation probabilities P qu (t, n) =ρ nn (t)(n = 1, 2, . . . , N) are calculated. In the current paper, T = 0 is adopted as QED devices work at extremely low temperatures and energies of electronic excitation and phonons are high in comparison with the thermal energy k B T in a wide range [44]. When the inital state is

Conclusion
In this work, we have built an extended Holstein model to describe a one-dimensional qubit chain with alternating site energies under the influence of a constant, electric field and high-frequency, dispersionless phonons. To investigate the transient polaron dynamics in the extended Holstein model, we have employed a time-dependent variational method based on an extension of the Davydov ansatz that becomes numerically accurate for sufficiently high values of the multiplicity parameter. This quantity is a measure of the complexity of the trial wave function and is related to the number of variational degrees of freedom in the ansatz. We have first studied the spectral properties and the band structure of the electronic Hamiltonian under various conditions. A tight-binding model (equation (1)) with near-neighbours interactions in the reciprocal and Wannier-Stark representations has been employed. The energy difference between consecutive site energies equals to the gap between the two minibands at the edges of the Brillouin zone. We have discussed the dependence of the energy offset between the Wannier-Stark energy ladders on the band gaps, and calculated the values for which the commensurate condition is satisfied, so that the time period of BZOs is equal to the period of BOs. We observe in our simulations that the band gap has a strong influence on the dynamics of BZOs.
After considering the coupling to the dispersionless optical phonon, we have carefully checked the validity of the multi-D2 ansatz by performing convergence tests and relative error calculations in order to demonstrate the accuracy of the method employed. Overall, the BZOs patterns and the corresponding breathing modes are smeared out by phonons. For weak qubit-phonon coupling strengths, the time periodicity of electronic motion is maintained and the spatial periodicity starts to break down when the coupling strength becomes comparable to the electronic interaction between qubits. Finally, we have numerically simulated the dissipative dynamics of BZOs employing a Markovian master equation and a vibrational ohmic environment. , and its value is smaller than that of P qu (t, n) by two orders of magnitude. The comparison here presents the nearly converged results calculated from D M=32 2 ansatz, showing the superior accuracy of our variational approach. Second, validity of the variational approach is carefully checked by testing how faithfully the calculations obey the Schrödinger equation. At the moment t, |Ψ(t) serves as the real wave function, and we assume the trial state |D M 2 (t) = |Ψ(t) . Then a deviation vector δ(t) is introduced to measure the accuracy of dynamics obtained by the multi-D2 trial state, in which the vectors γ(t) and χ(t) follow the Dirac-Frenkel time-dependent variational dynamics γ(t) = ∂|D M 2 /∂t in equation (10) and the Schrödinger equation χ(t) = ∂|Ψ(t) /∂t = 1 i Ĥ |Ψ(t) , respectively. Based on the Schrödinger equation as well as the relationship |Ψ(t) = |D M 1,2 (t) at the time t, the deviation vector δ(t) can be written as Thus, we can use the amplitude of the deviation vector Δ(t) = δ(t) to estimate the accuracy of the results calculated by the multi-D2 ansatz. To better verify the deviation in the parameter space (δ, J, g), we define a dimensionless relative deviation Σ as follows Σ = max{Δ(t)} mean{N err (t)} , t ∈ [0, t max ]. (B3) in which N err (t) = χ(t) denotes the amplitude of the time derivative of the wave function, As illustrated in figure 17, the relative deviation Σ drops when the multiplicity M of the multi-D2 ansatz grows. As supported by the relationship of Σ ∼ M −β with an exponent of β = 0.77 in the inset of figure 17,