Heating in Periodically Driven Spin Chains

This example shows how to code up the periodically-driven spin Hamiltonian:

\[\begin{split}H(t) &=& \left\{ \begin{array}{cl} J\sum_{j=0}^{L-1} \sigma^z_j\sigma^z_{j+1} + h\sum_{j=0}^{L-1}\sigma^z, & t\in[-T/4,\phantom{3}T/4] \\ \kern-8em g\sum_{j=0}^{L-1} \sigma^x_j, & t\in[\phantom{-} T/4,3T/4] \end{array} \right\} \mathrm{mod}\ T,\nonumber\\ &=& \sum_{j=0}^{L-1} \frac{1}{2}\left(J \sigma^z_j\sigma^z_{j+1} + h\sigma^z + g\sigma^x_j\right) + \frac{1}{2}\text{sgn}\left[\cos\Omega t\right]\left( J \sigma^z_j\sigma^z_{j+1} + h\sigma^z - g\sigma^x_j \right).\end{split}\]

Details about the code below can be found in SciPost Phys. 2, 003 (2017).

Script

download script

  1##################################################################################
  2#                            example 1                                           #
  3#     In this example we show how to use some of QuSpin's tools for studying     #
  4#     Floquet systems by analysing the heating in a periodically driven          #
  5#     spin chain. We also show how to construct more complicated multi-spin      #
  6#     interactions using QuSpin's interface.                                     #
  7##################################################################################
  8from quspin.operators import hamiltonian  # Hamiltonians and operators
  9from quspin.basis import spin_basis_1d  # Hilbert space spin basis
 10from quspin.tools.measurements import obs_vs_time, diag_ensemble  # t_dep measurements
 11from quspin.tools.Floquet import Floquet, Floquet_t_vec  # Floquet Hamiltonian
 12import numpy as np  # generic math functions
 13
 14#
 15##### define model parameters #####
 16L = 14  # system size
 17J = 1.0  # spin interaction
 18g = 0.809  # transverse field
 19h = 0.9045  # parallel field
 20Omega = 4.5  # drive frequency
 21
 22
 23#
 24##### set up alternating Hamiltonians #####
 25# define time-reversal symmetric periodic step drive
 26def drive(t, Omega):
 27    return np.sign(np.cos(Omega * t))
 28
 29
 30drive_args = [Omega]
 31# compute basis in the 0-total momentum and +1-parity sector
 32basis = spin_basis_1d(L=L, a=1, kblock=0, pblock=1)
 33# define PBC site-coupling lists for operators
 34x_field_pos = [[+g, i] for i in range(L)]
 35x_field_neg = [[-g, i] for i in range(L)]
 36z_field = [[h, i] for i in range(L)]
 37J_nn = [[J, i, (i + 1) % L] for i in range(L)]  # PBC
 38# static and dynamic lists
 39static = [["zz", J_nn], ["z", z_field], ["x", x_field_pos]]
 40dynamic = [
 41    ["zz", J_nn, drive, drive_args],
 42    ["z", z_field, drive, drive_args],
 43    ["x", x_field_neg, drive, drive_args],
 44]
 45# compute Hamiltonians
 46H = 0.5 * hamiltonian(static, dynamic, dtype=np.float64, basis=basis)
 47#
 48##### set up second-order van Vleck Floquet Hamiltonian #####
 49# zeroth-order term
 50Heff_0 = 0.5 * hamiltonian(static, [], dtype=np.float64, basis=basis)
 51# second-order term: site-coupling lists
 52Heff2_term_1 = [[+(J**2) * g, i, (i + 1) % L, (i + 2) % L] for i in range(L)]  # PBC
 53Heff2_term_2 = [[+J * g * h, i, (i + 1) % L] for i in range(L)]  # PBC
 54Heff2_term_3 = [[-J * g**2, i, (i + 1) % L] for i in range(L)]  # PBC
 55Heff2_term_4 = [[+(J**2) * g + 0.5 * h**2 * g, i] for i in range(L)]
 56Heff2_term_5 = [[0.5 * h * g**2, i] for i in range(L)]
 57# define static list
 58Heff_static = [
 59    ["zxz", Heff2_term_1],
 60    ["xz", Heff2_term_2],
 61    ["zx", Heff2_term_2],
 62    ["yy", Heff2_term_3],
 63    ["zz", Heff2_term_2],
 64    ["x", Heff2_term_4],
 65    ["z", Heff2_term_5],
 66]
 67# compute van Vleck Hamiltonian
 68Heff_2 = hamiltonian(Heff_static, [], dtype=np.float64, basis=basis)
 69Heff_2 *= -np.pi**2 / (12.0 * Omega**2)
 70# zeroth + second order van Vleck Floquet Hamiltonian
 71Heff_02 = Heff_0 + Heff_2
 72#
 73##### set up second-order van Vleck Kick operator #####
 74Keff2_term_1 = [[J * g, i, (i + 1) % L] for i in range(L)]  # PBC
 75Keff2_term_2 = [[h * g, i] for i in range(L)]
 76# define static list
 77Keff_static = [["zy", Keff2_term_1], ["yz", Keff2_term_1], ["y", Keff2_term_2]]
 78Keff_02 = hamiltonian(Keff_static, [], dtype=np.complex128, basis=basis)
 79Keff_02 *= np.pi**2 / (8.0 * Omega**2)
 80#
 81##### rotate Heff to stroboscopic basis #####
 82# e^{-1j*Keff_02} Heff_02 e^{+1j*Keff_02}
 83HF_02 = Heff_02.rotate_by(Keff_02, generator=True, a=1j)
 84#
 85##### define time vector of stroboscopic times with 100 cycles #####
 86t = Floquet_t_vec(Omega, 100, len_T=1)  # t.vals=times, t.i=init. time, t.T=drive period
 87#
 88##### calculate exact Floquet eigensystem #####
 89t_list = (
 90    np.array([0.0, t.T / 4.0, 3.0 * t.T / 4.0]) + np.finfo(float).eps
 91)  # times to evaluate H
 92dt_list = np.array(
 93    [t.T / 4.0, t.T / 2.0, t.T / 4.0]
 94)  # time step durations to apply H for
 95Floq = Floquet(
 96    {"H": H, "t_list": t_list, "dt_list": dt_list}, VF=True
 97)  # call Floquet class
 98VF = Floq.VF  # read off Floquet states
 99EF = Floq.EF  # read off quasienergies
100#
101##### calculate initial state (GS of HF_02) and its energy
102EF_02, psi_i = HF_02.eigsh(k=1, which="SA", maxiter=1e4)
103psi_i = psi_i.reshape((-1,))
104#
105##### time-dependent measurements
106# calculate measurements
107Sent_args = {"basis": basis, "chain_subsys": [j for j in range(L // 2)]}
108# meas = obs_vs_time((psi_i,EF,VF),t.vals,{"E_time":HF_02/L},Sent_args=Sent_args)
109# """
110# alternative way by solving Schroedinger's eqn
111psi_t = H.evolve(psi_i, t.i, t.vals, iterate=True, rtol=1e-9, atol=1e-9)
112meas = obs_vs_time(psi_t, t.vals, {"E_time": HF_02 / L}, Sent_args=Sent_args)
113# """
114# read off measurements
115Energy_t = meas["E_time"]
116Entropy_t = meas["Sent_time"]["Sent"]
117#
118##### calculate diagonal ensemble measurements
119DE_args = {"Obs": HF_02, "Sd_Renyi": True, "Srdm_Renyi": True, "Srdm_args": Sent_args}
120DE = diag_ensemble(L, psi_i, EF, VF, **DE_args)
121Ed = DE["Obs_pure"]
122Sd = DE["Sd_pure"]
123Srdm = DE["Srdm_pure"]
124#
125##### plot results #####
126import matplotlib.pyplot as plt
127import pylab
128
129# define legend labels
130str_E_t = "$\\mathcal{E}(lT)$"
131str_Sent_t = "$s_\mathrm{ent}(lT)$"
132str_Ed = "$\\overline{\mathcal{E}}$"
133str_Srdm = "$\\overline{s}_\mathrm{rdm}$"
134str_Sd = "$s_d^F$"
135# plot infinite-time data
136fig = plt.figure()
137plt.plot(t.vals / t.T, Ed * np.ones(t.vals.shape), "b--", linewidth=1, label=str_Ed)
138plt.plot(t.vals / t.T, Srdm * np.ones(t.vals.shape), "r--", linewidth=1, label=str_Srdm)
139plt.plot(t.vals / t.T, Sd * np.ones(t.vals.shape), "g--", linewidth=1, label=str_Sd)
140# plot time-dependent data
141plt.plot(t.vals / t.T, Energy_t, "b-o", linewidth=1, label=str_E_t, markersize=3.0)
142plt.plot(t.vals / t.T, Entropy_t, "r-s", linewidth=1, label=str_Sent_t, markersize=3.0)
143# label axes
144plt.xlabel("$\\#\ \\mathrm{periods}\\ l$", fontsize=18)
145# set y axis limits
146plt.ylim([-0.6, 0.7])
147# display legend
148plt.legend(loc="lower right", ncol=2, fontsize=18)
149# update axis font size
150plt.tick_params(labelsize=16)
151# turn on grid
152plt.grid(True)
153# save figure
154plt.tight_layout()
155plt.savefig("example2.pdf", bbox_inches="tight")
156# show plot
157# plt.show()
158plt.close()