TEBD for time-dependent Hamiltonians

[abhiroop]

Hello,
I am trying to implement a program which calculates the time-dependent Hubbard Hamiltonian with hopping t*exp(-i*phi(t)) if the system is driven by some periodic electric field.

In tenpy/models/fermions_hubbard.py we put t as well as t2 (defined later)

Code: Select all

      for u1, u2, dx in self.lat.nearest_neighbors:
            self.add_coupling(t, u1, 'Cdu', u2, 'Cu', dx, 'JW', True)
            self.add_coupling(t2, u1, 'Cdu', u2, 'Cu', -dx, 'JW', True)  # t2  has been put in place of t
            self.add_coupling(t, u1, 'Cdd', u2, 'Cd', dx, 'JW', True)
            self.add_coupling(t2, u1, 'Cdd', u2, 'Cd', -dx, 'JW', True)  # h.c.
            self.add_coupling(V, u1, 'Ntot', u2, 'Ntot', dx)

Then, in tenpy/algorithms/tebd.py, in the function tebd.run(), we set Nsteps = 1
Finally, in the main program, we can have a loop for time steps

Code: Select all

    time_steps = 10000   # or any desired number of steps
    for tsteps in range(time_steps):
         t =cmath.exp(-1.0*math.sin(omega*deltat)*1j)    # here phi(t) = math.sin(omega*deltat) related to my problem, anything can be assigned
         t2 =cmath.exp(1.0*math.sin(omega*deltat)*1j)
         model_params = dict(L=L, t=t, t2=t2, U=U, V=V, mu=0., bc_MPS='finite', cons_N='N', cons_Sz='Sz')
         M = Hubbard_efield_Chain(model_params)
         eng = tebd.Engine(psi, M, tebd_params)
         eng.run()
         psi = eng.psi

I guess this will be useful for simulating time-driven (Floquet type) systems.

[Johannes]

Looks like a good thing.
Actually, if you update your local repository, you'll find that the code for the FermiHubbardModel has been moved to tenpy.models.hubbard.FermiHubbardModel. There it says

Python: Select all

        for u1, u2, dx in self.lat.nearest_neighbors:
            self.add_coupling(-t, u1, 'Cdu', u2, 'Cu', dx)
            self.add_coupling(-np.conj(t), u2, 'Cdu', u1, 'Cu', -dx)  # h.c.
            self.add_coupling(-t, u1, 'Cdd', u2, 'Cd', dx)
            self.add_coupling(-np.conj(t), u2, 'Cdd', u1, 'Cd', -dx)  # h.c.
            self.add_coupling(V, u1, 'Ntot', u2, 'Ntot', dx)

First of all, note that the sign of the hopping changed, so positive t means the usual \(-t c^\dagger_i c_j + h.c.\).
Second, note that it uses -np.conj(t), which allows to choose t complex, as in your case. In that way, we don't need to define a second t2 at all, which gives less room for errors.

In general, I planed to provide a way to do such a time dependent simulation for a given, time dependent model paramater (whichever it is), but didn't find time to do it yet. I'd like to base it on a general simulation class, see Issue #46.

Thank you for providing the example how to do that with TEBD, it looks perfectly right.
I'd like to mention that in the tebd_params, you should choose N_steps=1 and order=1 to update the time-dependent Hamiltonian at each small time step. I've never really thought about the correct way to do TEBD with time-dependent Hamiltonians second-order or even fourth-order in the time step. If you have any insight in that direction, please let me know

[GRomeroH]

Hi there, I hope that you are doing well. I found this discussion very helpful. Indeed, I implemented a time-dependent version of the Bose-Hubbard model following the first suggestion

Python: Select all

for x in range(len(tspan)):
        
        t1 = t * np.cos(omegaD * tspan[x])
        t2 = t * np.cos(omegaD * tspan[x])
       
        "Full modulation"
        t_array = np.array([(t1 if i % 2 == 0 else t2) for i in range(L-1)])
        model_params = dict( filling =filling,  t=t_array, U=U, V=V, mu=mu, bc_MPS=bc_MPS, L=L, verbose=False)
        M = BoseHubbardChain(model_params)
        eng = tebd.Engine(psi, M, tebd_params)
        start_time = time.time()
        eng.run()
        psi = eng.psi

It works pretty well indeed. I am using it for studying Floquet prethermalization in the Bose-Hubbard chain using the second-order Suzuki-Trotter expansion. Though I am computing the full dynamics, I was wondering whether it is possible to compute the evolution operator at one period of the drive so as to obtain the stroboscopic dynamics, do you think is it possible to code it?