Integrability Breaking and Thermalising Dynamics in the 2D Transverse-Field Ising Model
This example shows how to code up the time-periodic 2D transverse-field Ising Hamiltonian:
\[\begin{split}H(t)=\Bigg\{ \begin{array}{cc}
H_{zz} +AH_x,& \qquad t\in[0,T/4), \\
H_{zz} -AH_x,& \qquad t\in[T/4,3T/4),\\
H_{zz} +AH_x,& \qquad t\in[3T/4,T)
\end{array}\end{split}\]
where
\[H_{zz} = -\sum_{\langle ij\rangle} S^z_iS^z_{j}, \qquad H_{x} = -\sum_{i}S^x_i.\]
Details about the code below can be found in SciPost Phys. 7, 020 (2019).
Script
1#####################################################################
2# example 9 #
3# In this script we demonstrate how to use QuSpin's #
4# general basis class to construct user-defined symmetry sectors. #
5# We study thermalisation in the 2D transverse-field Ising model #
6# with periodic boundary conditions. #
7#####################################################################
8from quspin.operators import hamiltonian, exp_op # operators
9from quspin.basis import spin_basis_1d, spin_basis_general # spin basis constructor
10from quspin.tools.measurements import obs_vs_time # calculating dynamics
11from quspin.tools.Floquet import Floquet_t_vec # period-spaced time vector
12import numpy as np # general math functions
13import matplotlib.pyplot as plt # plotting library
14
15#
16###### define model parameters ######
17L_1d = 16 # length of chain for spin 1/2
18Lx, Ly = 4, 4 # linear dimension of spin 1/2 2d lattice
19N_2d = Lx * Ly # number of sites for spin 1/2
20Omega = 2.0 # drive frequency
21A = 2.0 # drive amplitude
22#
23###### setting up user-defined symmetry transformations for 2d lattice ######
24s = np.arange(N_2d) # sites [0,1,2,....]
25x = s % Lx # x positions for sites
26y = s // Lx # y positions for sites
27T_x = (x + 1) % Lx + Lx * y # translation along x-direction
28T_y = x + Lx * ((y + 1) % Ly) # translation along y-direction
29P_x = x + Lx * (Ly - y - 1) # reflection about x-axis
30P_y = (Lx - x - 1) + Lx * y # reflection about y-axis
31Z = -(s + 1) # spin inversion
32#
33###### setting up bases ######
34basis_1d = spin_basis_1d(L_1d, kblock=0, pblock=1, zblock=1) # 1d - basis
35basis_2d = spin_basis_general(
36 N_2d,
37 kxblock=(T_x, 0),
38 kyblock=(T_y, 0),
39 pxblock=(P_x, 0),
40 pyblock=(P_y, 0),
41 zblock=(Z, 0),
42) # 2d - basis
43# print information about the basis
44print("Size of 1D H-space: {Ns:d}".format(Ns=basis_1d.Ns))
45print("Size of 2D H-space: {Ns:d}".format(Ns=basis_2d.Ns))
46#
47###### setting up operators in hamiltonian ######
48# setting up site-coupling lists
49Jzz_1d = [[-1.0, i, (i + 1) % L_1d] for i in range(L_1d)]
50hx_1d = [[-1.0, i] for i in range(L_1d)]
51#
52Jzz_2d = [[-1.0, i, T_x[i]] for i in range(N_2d)] + [
53 [-1.0, i, T_y[i]] for i in range(N_2d)
54]
55hx_2d = [[-1.0, i] for i in range(N_2d)]
56# setting up hamiltonians
57# 1d
58Hzz_1d = hamiltonian([["zz", Jzz_1d]], [], basis=basis_1d, dtype=np.float64)
59Hx_1d = hamiltonian([["x", hx_1d]], [], basis=basis_1d, dtype=np.float64)
60# 2d
61Hzz_2d = hamiltonian([["zz", Jzz_2d]], [], basis=basis_2d, dtype=np.float64)
62Hx_2d = hamiltonian([["x", hx_2d]], [], basis=basis_2d, dtype=np.float64)
63#
64###### calculate initial states ######
65# calculating bandwidth for non-driven hamiltonian
66[E_1d_min], psi_1d = Hzz_1d.eigsh(k=1, which="SA")
67[E_2d_min], psi_2d = Hzz_2d.eigsh(k=1, which="SA")
68# setting up initial states
69psi0_1d = psi_1d.ravel() # .astype(np.complex128)
70psi0_2d = psi_2d.ravel()
71#
72###### time evolution ######
73# stroboscopic time vector
74nT = 200 # number of periods to evolve to
75t = Floquet_t_vec(Omega, nT, len_T=1) # t.vals=t, t.i=initial time, t.T=drive period
76# creating generators of time evolution using exp_op class
77U1_1d = exp_op(Hzz_1d + A * Hx_1d, a=-1j * t.T / 4)
78U2_1d = exp_op(Hzz_1d - A * Hx_1d, a=-1j * t.T / 2)
79U1_2d = exp_op(Hzz_2d + A * Hx_2d, a=-1j * t.T / 4)
80U2_2d = exp_op(Hzz_2d - A * Hx_2d, a=-1j * t.T / 2)
81
82
83# user-defined generator for stroboscopic dynamics
84def evolve_gen(psi0, nT, *U_list):
85 yield psi0
86 for i in range(nT): # loop over number of periods
87 for U in U_list: # loop over unitaries
88 psi0 = U.dot(psi0)
89 yield psi0
90
91
92# get generator objects for time-evolved states
93psi_1d_t = evolve_gen(psi0_1d, nT, U1_1d, U2_1d, U1_1d)
94psi_2d_t = evolve_gen(psi0_2d, nT, U1_2d, U2_2d, U1_2d)
95#
96###### compute expectation values of observables ######
97# measure Hzz as a function of time
98Obs_1d_t = obs_vs_time(psi_1d_t, t.vals, dict(E=Hzz_1d), return_state=True)
99Obs_2d_t = obs_vs_time(psi_2d_t, t.vals, dict(E=Hzz_2d), return_state=True)
100# calculating the entanglement entropy density
101Sent_time_1d = basis_1d.ent_entropy(Obs_1d_t["psi_t"], sub_sys_A=range(L_1d // 2))[
102 "Sent_A"
103]
104Sent_time_2d = basis_2d.ent_entropy(Obs_2d_t["psi_t"], sub_sys_A=range(N_2d // 2))[
105 "Sent_A"
106]
107# calculate entanglement entropy density
108s_p_1d = np.log(2) - 2.0 ** (-L_1d // 2) / L_1d
109s_p_2d = np.log(2) - 2.0 ** (-N_2d // 2) / N_2d
110#
111###### plotting results ######
112plt.plot(
113 t.strobo.inds,
114 (Obs_1d_t["E"].real - E_1d_min) / (-E_1d_min),
115 marker=".",
116 markersize=5,
117 label="$S=1/2$",
118)
119plt.plot(
120 t.strobo.inds,
121 (Obs_2d_t["E"].real - E_2d_min) / (-E_2d_min),
122 marker=".",
123 markersize=5,
124 label="$S=1$",
125)
126plt.grid()
127plt.ylabel("$Q(t)$", fontsize=20)
128plt.xlabel("$t/T$", fontsize=20)
129plt.savefig("TFIM_Q.pdf")
130plt.figure()
131plt.plot(t.strobo.inds, Sent_time_1d / s_p_1d, marker=".", markersize=5, label="$1d$")
132plt.plot(t.strobo.inds, Sent_time_2d / s_p_2d, marker=".", markersize=5, label="$2d$")
133plt.grid()
134plt.ylabel("$s_{\\mathrm{ent}}(t)/s_\\mathrm{Page}$", fontsize=20)
135plt.xlabel("$t/T$", fontsize=20)
136plt.legend(loc=0, fontsize=16)
137plt.tight_layout()
138plt.savefig("TFIM_S.pdf")
139# plt.show()
140plt.close()