Optimal time-dependent lattice models for nonequilibrium dynamics

Lattice models are abundant in theoretical and condensed-matter physics. Generally, lattice models contain time-independent hopping and interaction parameters that are derived from the Wannier functions of the noninteracting problem. Here, we present a new concept based on time-dependent Wannier functions and the variational principle that leads to optimal time-dependent lattice models. As an application, we use the Bose-Hubbard model with time-dependent Wannier functions to study a quench scenario involving higher bands. We find a separation of times scales in the dynamics and show that under some circumstances the multi-band nonequilibrium dynamics of a quantum system can be obtained essentially at the cost of a single-band model.

Quantum phenomena of particles in periodic and quasi-periodic potentials are central themes in theoretical physics. Alone the question about the nature of the many-body ground state in such potentials has been investigated for decades [1]. Ultracold atoms in optical lattices are highly controllable realizations of such many-body systems and allow to directly monitor their nonequilibrium dynamics. Moreover, they promise to provide insight into the physics of real solids. Lattice models play a crucial role in the current physical understanding of such systems. At the heart of a lattice model is the idea of lattice site localized orbitals which are commonly referred to as Wannier functions [2,3]. In a lattice model Wannier functions are used as a single-particle basis for the many-body Hamiltonian. Among the most famous lattice models are certainly the single-band bosonic [4] and fermionic [5] versions of the Hubbard model, both of which predict phase transitions between insulating and superfluid phases. Recently there has been large interest in the multi-band physics of ultracold atoms by both theorists and experimentalists alike [6][7][8][9][10][11][12][13][14][15]. Interparticle interactions already couple the ground state to higher bands [6,9,11], and it was shown recently that even very weak interactions lead to multi-band physics in nonequilibrium problems [14].
Multi-band phenomena seem to be abundant and we would like to address the question of the optimal way to deal with them theoretically. This calls for optimal lattice models that incorporate multi-band physics effectively. Here, we present a completely new concept that leads to variationally optimal lattice models for multi-band nonequilibrium dynamics and statics problems. Our idea is to allow the Wannier functions of a given lattice model to depend on time and to determine the many-body dynamics from the variational principle. The concept is applicable to bosons and fermions as well as any number of lattice sites, particles and bands. As an application, we use the Bose-Hubbard model with timedependent Wannier functions to study a sudden quench of the interaction in a double-well problem that involves higher bands. We find a separation of time-scales in the dynamics.
On short time-scales higher-band excitations lead to rapid density oscillations, not captured on the single-band level. On longer time-scales we observe oscillations between condensed and fragmented states. The results are compared with numerically exact results of the time-dependent many-body Schrödinger equation. We thereby show that when the lattice is sufficiently deep and the interaction couples to higher bands the multi-band nonequilibrium dynamics of a quantum system can be obtained essentially at the cost of a single-band model. Already the physics of this example proves to be much richer than what the single-band Bose-Hubbard (BH) model can access.
Consider a one-dimensional (1D) lattice of M sites. We make the following ansatz for the time-dependent Wannier functions where w α k (x) denotes the conventional Wannier function of the band α at lattice site k and the d α k (t) its time-dependent amplitude. If ν is restricted to one, the ansatz (1) reduces to the conventional Wannier functions of the lowest band. Using (1) as a basis, the ansatz for the many-boson wave function becomes |Ψ(t) = n C n (t) | n; t , where the sum is over all permanents | n; t = |n 1 , n 2 , . . . , n M ; t of N bosons residing in M time-dependent Wannier functions. We writeb † k (t) (b k (t)) for the operator that creates (annihilates) a boson in the time-dependent Wannier function w k (x, t) andn k (t) =b † k (t)b k (t) for the corresponding occupation number operator. The bosonic commutation relations [b k (t), b † q (t)] = δ kq hold at all times. For simplicity we use dimensionless units in which = m = 1. Analogous to the derivation of the BH model, the Hamiltonian of the time-dependent BH (TDBH) model can then be derived from the full many-body is the one-body part of the Hamiltonian including the external potential V (x). The time-dependent parameters in Eq. (2) are defined as Contrary to the BH model, the hopping parameter J kk+1 (t), the onebody energy ǫ k (t) and the interaction parameter U kkkk (t) are now time-dependent and can vary from one lattice site to another. Note also that the hopping parameter J kk+1 (t) will generally not be real. It is then possible to derive coupled equations of motion for the timedependent Wannier functions and the coefficients of the many-body wave function from the variational principle [16,17]: where C(t) = {C n (t)}. The derivation and a discussion of Eqs. (3) are given in the appendix. The operatorsP k (t) = ν β=1 |w β k w β k | − |w k (t) w k (t)| appearing in Eq. (3) are projection operators, U kk (x, t) = λ 0 |w k (x, t)| 2 and H n n ′ (t) = n; t Ĥ T DBH n ′ ; t . Similarly, are matrix elements of the first-and second-order reduced density matrices. We stress that the optimal lattice model governed by Eqs. (1-3) preserves all conservation laws. In particular, particle-number is a good quantum number, energy is conserved when the original Hamiltonian is time-independent, and spatial symmetries of this Hamiltonian are preserved in time for initial conditions preserving them. Summarizing the above, we have obtained the TDBH model for which the wave function and the parameters of the Hamiltonian (2) evolve according to the variational principle. The dynamics is determined by the coupled equations of motion (3). Now we would like to demonstrate the above ideas in an apparently simple two-site example. The dynamics of bosons on two lattice sites has recently been studied intensively using the BH model [18][19][20][21][22][23][24][25][26][27][28]. A question that arises naturally in dynamical problems is the response of a quantum system to a perturbation. When the system is in the ground state and the perturbation is realized as a variation of one of the parameters in the Hamiltonian, the problem is known as a quantum quench. Quenches have recently received growing attention in the context of quantum gases in lattices [29][30][31][32]. Generally, the quench dynamics is studied using single-band Hubbard models or simplifications thereof. Here, we study a quench going beyond the BH model and would like to investigate what new physics appears.
Here, V (x) = V 0 cos 2 (x) is an external potential on an interval of length 2π that realizes a ring lattice of M = 2 sites (denoted L and R) with V 0 = 25E r , where E r = 1/2 is the recoil energy. The splitting 2J = 2.08 × 10 −3 of the ground and the first excited singleparticle state of V (x) determines the Rabi oscillation period t Rabi = π/J. To characterize the interaction strength we use the parameter λ = λ 0 (N − 1). The quench is then realized as a sudden change from λ = 0 to λ = 0.6 in a system of N = 20 bosons, initially prepared in the noninteracting ground state of V (x). Within the BH model this quench corresponds to a change from U/J = 0 to U/J = 25.8. We will now study and compare the dynamics of the TDBH and the BH model. For the TDBH model the number of bands ν used for the expansion of each time-dependent Wannier function in Eq. (1) is increased until convergence is reached. In the present example we found that the results for J(t), U(t), the one-particle density and the momentum distribution do not change visibly anymore for ν ≃ 5, and ν = 10 was used throughout for full convergence. Furthermore, the spatial symmetry of the problem implies that only bands with even Wannier functions, α = 1, 3, 5, 7, 9, contribute.
times. All TDBH computations were carried out on a standard desktop computer and none of the computations took longer than a few hours.
On physical grounds it can be anticipated that the sudden increase of the (repulsive) interaction strength will lead to breathing oscillations of the one-particle density ρ (1) (x|x; t) in each of the wells. However, as the final interaction strength is not very strong these density oscillations should have a small amplitude only. It is important to note that for the BH model the ratio of the interaction and the hopping parameter U/J is the only relevant parameter and also that it is constant in time after the quench. In the TDBH model on the other hand we expect U(t)/J(t) to vary, because the Wannier functions can change in time.
This is indeed the case. Fig. 1 shows the parameters U(t) and J(t) of the TDBH model and their ratio U(t)/J(t) as a function of time. As expected, right after the quench the density in each well expands and hence U(t) decreases, while J(t) increases. Both parameters then oscillate at a high frequency. Note the short time-scale. The interaction broadens the density and therefore U(t) (J(t)) is always smaller (larger) than its value at time t = 0. As J(t) We now turn to the question about the nature of the quantum state after the quench and to longer time-scales. We first focus on the eigenvalues of the first-order reduced density matrix, the natural occupation numbers n (1) i (t), which determine the degree to which the system is condensed or fragmented [33,34]. In Fig. 2  For details see the literature [17,35]. Note that the time-scale is now more than a hundred times larger compared to Fig. 1. It is clearly seen that all three results display oscillations between fragmented and partially condensed states. For the first few oscillations the three results essentially coincide, which implies that for this problem the natural occupation numbers are not very sensitive to the previously discussed multi-band excitations. However, after a few oscillations between condensed and fragmented states the BH result deviates substantially from the exact one. On the other hand, the TDBH result follows the exact one closely for many oscillations. In Fig. 2 (bottom) the accumulated error of the BH model, defined as the accumulated difference 1 1,exact (t ′ )| between the largest natural occupation number of the BH model and the numerically exact result are shown, together with the analogous quantity of the TDBH model. The accumulated error of the TDBH model is always smaller than that of the BH model and grows more slowly, which means that also on longer time-scales the TDBH model captures the true physics much better than the BH model.
We will now discuss the implications of the differences of the natural occupation numbers on observables. The rapid density oscillations discussed earlier have a very small amplitude, and consequently hardly any dynamics is visible in real space. However, a substantial dynamics occurs in momentum space. Fig. 3 shows the one-particle momentum distribution ρ (1) (k|k; t) of the exact, the TDBH and the BH result at t = 0 where they all coincide, and at later times where differences have developed. The momentum distribution of the TDBH model is always closer to the exact one and follows it for a long time, whereas that of the BH model deviates quickly. Here, the TDBH model allows to obtain the precise multi-band nonequilibrium dynamics essentially at the cost of a single-band model. Of course, in order to obtain the exact result the terms that were neglected so far in Eq. (2) have to be included into the model, e.g., the terms responsible for correlated tunneling. Thus, the advantages of the TDBH model over the BH model have been clearly proven.
In this work we have generalized the concept of lattice models by letting their Wannier functions become time-dependent. For such lattice models equations of motion can be derived from the variational principle, leading to variationally optimal lattice models that can efficiently incorporate multi-band physics. The concept is applicable to both nonequilibrium dynamics as well as statics. In the latter case the ground state can be computed using imaginary time-propagation. The variational principle ensures that for identical initial conditions the optimal lattice model improves on the original lattice model, if the latter contains the terms of the full many-body Hamiltonian relevant to the problem at hand. Since the numerical effort of the optimal and standard lattice models are comparable, the former can also be used to assess the quality of the latter: if the respective results differ noticeably, the standard lattice model is inapplicable. As an application, we have in disordered media or in the case of the fermionic Hubbard model. As a further direction we suggest to combine 1D optimal lattice models with the recently developed adaptive timedependent density-matrix renormalization group methods [36][37][38] to treat non-equilibrium dynamics of larger lattice systems, as is presently done utilizing standard lattice models.
Appendix: Derivation and discussion of the equations of motion of the TDBH model We will now derive Eqs. (3) and discuss their properties. Using the ansatz for the manybody wave function |Ψ(t) = n C n (t) | n; t given above, the action functional of the TDBH model reads where the Lagrange multipliers are introduced to ensure orthonormality between the timedependent Wannier functions. Note that Eq. (1) implies that time-dependent Wannier functions on different lattice sites are orthogonal by construction. In order to take functional derivatives it is useful to express the expectation value Ψ|Ĥ T DBH − i ∂ ∂t |Ψ as an explicit function of the amplitudes {d α k (t)} and expansion coefficients {C n (t)}. On the one hand Ψ|Ĥ T DBH − i ∂ ∂t |Ψ can be written as: Using Eq. (A.2) we now require stationarity of the action functional with respect to variations of the coefficients, 0 = δS δC * n (t) , and obtain: On the other hand Ψ|Ĥ T DBH − i ∂ ∂t |Ψ can also be written as: By including the additional phase factor e − t dt ′ w k (t)|ẇ k (t) into the definition of |w k (t) , we find that w k (t)|ẇ k (t) = 0, and that Eqs.  The TDBH result follows the exact one closely and in particular reproduces the frequencies of the oscillations well. The BH result deviates quickly from the exact one. Bottom: the accumulated error (see text) of the TDBH model is always smaller than that of the BH model and grows much more slowly. All quantities shown are dimensionless.