Hello, I have been trying to find a way to replicate the quasienergy spectrum of a periodically driven dimers chain, as depicted in Fig 3 of https://doi.org/10.1103/PhysRevLett.110.200403, using tkwant. I started with a basic time-independent model import kwant from matplotlib import pyplot import numpy as np a = 1.0 b = 0.2 * a lat = kwant.lattice.general([(0, a), (a, 0)], [(1, 0), (b, 0)]) lat_a, lat_b = lat.sublattices sym = kwant.TranslationalSymmetry((-a, 0)) dimers = kwant.Builder(sym) dimers[lat_a(-1, 0)] = 5 dimers[lat_b(-1, 0)] = 3 dimers[[kwant.builder.HoppingKind((0, 0), lat_a, lat_b)]] = -1.5 dimers[[kwant.builder.HoppingKind((0, 1), lat_a, lat_b)]] = -1.7 kwant.plot(dimers) dimers = dimers.finalized() bands = kwant.physics.Bands(dimers) momenta = np.linspace(-np.pi, np.pi, 101) energies = [bands(k) for k in momenta] pyplot.plot(momenta, energies) pyplot.grid("on") Now I want to introduce the time-dependent vector potential def vec_potential(time, freq, e_field): q = -1 A_0 = q * e_field / freq return A_0 * np.sin(freq * time) however, I have been stuck over how exactly to couple the two together and correctly modify the hoppings. I am also uncertain whether the existing kwant/tkwant API will allow me to extract quasi-energy spectrum directly. Any hints will be greatly appreciated. Best wishes, Valentin