Adiabatic Control of Parameters in Many-Body Localised Phases

This example shows how to code up the time-dependent disordered spin Hamiltonian:

\[\begin{split}H(t) &=& \sum_{j=0}^{L-2}\frac{J_{xy}}{2}\left(S^+_{j+1}S^-_{j} + \mathrm{h.c.}\right) + J_{zz}(t)S^z_{j+1}S^z_{j} + \sum_{j=0}^{L-1}h_jS^z_{j},\nonumber\\ J_{zz}(t) &=& (1/2 + vt)J_{zz}(0).\end{split}\]

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

UPDATE: The original version of this example discussed in the paper above has been deprecated as it is not compatible with python 3. Below we have implemented a new version of the example that uses more up-to-date features of QuSpin. We also leave a link down below to download the original script for comparison.

Script

download script (python 3)

download original script (python 2, depricated)

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
from __future__ import print_function, division
import sys,os
# line 4 and line 5 below are for development purposes and can be removed
qspin_path = os.path.join(os.getcwd(),"../../")
sys.path.insert(0,qspin_path)
##################################################################
#                            example 1                           #
#    In this script we demonstrate how to use QuSpin's           #
#    functionality to solve the time-dependent Schroedinger      #
#    equation with a time-dependent operator. We also showcase   #
#    a function which is used to calculate the entanglement      #
#    entropy of a pure state.                                    #
##################################################################
from quspin.operators import quantum_operator # Hamiltonians and operators
from quspin.basis import spin_basis_1d # Hilbert space spin basis
from quspin.tools.measurements import ent_entropy, diag_ensemble # entropies
from numpy.random import uniform,seed # pseudo random numbers
from joblib import delayed,Parallel # parallelisation
import numpy as np # generic math functions
from time import time # timing package
#
##### define simulation parameters #####
n_real=100 # number of disorder realisations
n_jobs=2 # number of spawned processes used for parallelisation
#
##### define model parameters #####
L=10 # system size
Jxy=1.0 # xy interaction
Jzz_0=1.0 # zz interaction at time t=0
h_MBL=3.9 # MBL disorder strength
h_ETH=0.1 # delocalised disorder strength
vs=np.logspace(-3.0,0.0,num=20,base=10) # log_2-spaced vector of ramp speeds
#
##### set up Heisenberg Hamiltonian with linearly varying zz-interaction #####
# define linear ramp function
def ramp(t,v):
	return (0.5 + v*t)
# compute basis in the 0-total magnetisation sector (requires L even)
basis = spin_basis_1d(L,m=0,pauli=False)
# define operators with OBC using site-coupling lists
J_zz = [[Jzz_0,i,i+1] for i in range(L-1)] # OBC
J_xy = [[Jxy/2.0,i,i+1] for i in range(L-1)] # OBC
# dictionary of operators for quantum_operator.
# define XXZ chain parameters
op_dict = dict(J_xy=[["+-",J_xy],["-+",J_xy]],J_zz=[["zz",J_zz]])
# define operators for local disordered field.
for i in range(L):
	op = [[1.0,i]]
	op_dict["hz"+str(i)] = [["z",op]]
# costruct the quantum_operator
H_XXZ = quantum_operator(op_dict,basis=basis,dtype=np.float64)
#
##### calculate diagonal and entanglement entropies #####
def realization(vs,H_XXZ,real):
	"""
	This function computes the entropies for a single disorder realisation.
	--- arguments ---
	vs: vector of ramp speeds
	H_XXZ: time-dep. Heisenberg Hamiltonian with driven zz-interactions
	basis: spin_basis_1d object containing the spin basis
	n_real: number of disorder realisations; used only for timing
	"""
	ti = time() # get start time
	#
	seed() # the random number needs to be seeded for each parallel process
	N = H_XXZ.basis.N
	basis = H_XXZ.basis
	hz = uniform(-1,1,size=N)
	# define parameters to pass into quantum_operator for 
	# hamiltonian at end of ramp
	pars_MBL = {"hz"+str(i):h_MBL*hz[i] for i in range(N)}
	pars_ETH = {"hz"+str(i):h_ETH*hz[i] for i in range(N)}
	# 
	pars_MBL["J_xy"] = 1.0
	pars_ETH["J_xy"] = 1.0
	# J_zz = 1 at end of the ramp for all velocities
	pars_MBL["J_zz"] = 1.0
	pars_ETH["J_zz"] = 1.0
	# diagonalize 
	E_MBL,V_MBL = H_XXZ.eigh(pars=pars_MBL)
	E_ETH,V_ETH = H_XXZ.eigh(pars=pars_ETH)
	# reset J_zz to be initial value:
	pars_MBL["J_zz"] = 0.5
	# get many-body bandwidth at t=0
	eigsh_args=dict(k=2,which="BE",maxiter=1E4,return_eigenvectors=False,pars=pars_MBL)
	Emin,Emax=H_XXZ.eigsh(**eigsh_args)
	# calculating middle of spectrum
	E_inf_temp=(Emax+Emin)/2.0
	# calculate nearest eigenstate to energy at infinite temperature
	E,psi_0=H_XXZ.eigsh(pars=pars_MBL,k=1,sigma=E_inf_temp,maxiter=1E4)
	psi_0=psi_0.reshape((-1,))

	run_MBL = []

	for v in vs:
		# update J_zz to be time-dependent operator
		pars_MBL["J_zz"] = (ramp,(v,))
		# get hamiltonian
		H = H_XXZ.tohamiltonian(pars=pars_MBL)
		# evolve state and calculate oberservables. 
		run_MBL.append(_do_ramp(psi_0,H,basis,v,E_MBL,V_MBL))

	run_MBL=np.vstack(run_MBL).T
	# reset J_zz to be initial value:
	pars_ETH["J_zz"] = 0.5
	# get many-body bandwidth at t=0
	eigsh_args=dict(k=2,which="BE",maxiter=1E4,return_eigenvectors=False,pars=pars_ETH)
	Emin,Emax=H_XXZ.eigsh(**eigsh_args)
	# calculating middle of spectrum
	E_inf_temp=(Emax+Emin)/2.0
	# calculate nearest eigenstate to energy at infinite temperature
	E,psi_0=H_XXZ.eigsh(pars=pars_ETH,k=1,sigma=E_inf_temp,maxiter=1E4)
	psi_0=psi_0.reshape((-1,))

	run_ETH = []

	for v in vs:
		# update J_zz to be time-dependent operator
		pars_ETH["J_zz"] = (ramp,(v,))
		# get hamiltonian
		H = H_XXZ.tohamiltonian(pars=pars_ETH)
		# evolve state and calculate oberservables. 
		run_ETH.append(_do_ramp(psi_0,H,basis,v,E_ETH,V_ETH))
	run_ETH=np.vstack(run_ETH).T
	# show time taken
	print("realization {0}/{1} took {2:.2f} sec".format(real+1,n_real,time()-ti))
	#
	return run_MBL,run_ETH
#
##### evolve state and evaluate entropies #####
def _do_ramp(psi_0,H,basis,v,E_final,V_final):
	"""
	Auxiliary function to evolve the state and calculate the entropies after the
	ramp.
	--- arguments ---
	psi_0: initial state
	H: time-dependent Hamiltonian
	basis: spin_basis_1d object containing the spin basis (required for Sent)
	E_final, V_final: eigensystem of H(t_f) at the end of the ramp t_f=1/(2v)
	"""
	# determine total ramp time
	t_f = 0.5/v 
	# time-evolve state from time 0.0 to time t_f
	psi = H.evolve(psi_0,0.0,t_f)
	# calculate entanglement entropy
	subsys = range(basis.L//2) # define subsystem
	Sent = basis.ent_entropy(psi,sub_sys_A=subsys)["Sent_A"]
	# calculate diagonal entropy in the basis of H(t_f)
	S_d = diag_ensemble(basis.L,psi,E_final,V_final,Sd_Renyi=True)["Sd_pure"]
	#
	return np.asarray([S_d,Sent])
#
##### produce data for n_real disorder realisations #####
# __name__ == '__main__' required to use joblib in Windows.
if __name__ == '__main__':

	# alternative way without parallelisation
	data = np.asarray([realization(vs,H_XXZ,i) for i in range(n_real)])
	"""
	data = np.asarray(Parallel(n_jobs=n_jobs)(delayed(realization)(vs,H_XXZ,basis,i) for i in range(n_real)))
	"""
	#
	run_MBL,run_ETH = zip(*data) # extract MBL and data
	# average over disorder
	mean_MBL = np.mean(run_MBL,axis=0)
	mean_ETH = np.mean(run_ETH,axis=0)
	#
	##### plot results #####
	import matplotlib.pyplot as plt
	### MBL plot ###
	fig, pltarr1 = plt.subplots(2,sharex=True) # define subplot panel
	# subplot 1: diag enetropy vs ramp speed
	pltarr1[0].plot(vs,mean_MBL[0],label="MBL",marker=".",color="blue") # plot data
	pltarr1[0].set_ylabel("$s_d(t_f)$",fontsize=22) # label y-axis
	pltarr1[0].set_xlabel("$v/J_{zz}(0)$",fontsize=22) # label x-axis
	pltarr1[0].set_xscale("log") # set log scale on x-axis
	pltarr1[0].grid(True,which='both') # plot grid
	pltarr1[0].tick_params(labelsize=16)
	# subplot 2: entanglement entropy vs ramp speed
	pltarr1[1].plot(vs,mean_MBL[1],marker=".",color="blue") # plot data
	pltarr1[1].set_ylabel("$s_\\mathrm{ent}(t_f)$",fontsize=22) # label y-axis
	pltarr1[1].set_xlabel("$v/J_{zz}(0)$",fontsize=22) # label x-axis
	pltarr1[1].set_xscale("log") # set log scale on x-axis
	pltarr1[1].grid(True,which='both') # plot grid
	pltarr1[1].tick_params(labelsize=16)
	# save figure
	fig.savefig('example1_MBL.pdf', bbox_inches='tight')
	#
	### ETH plot ###
	fig, pltarr2 = plt.subplots(2,sharex=True) # define subplot panel
	# subplot 1: diag enetropy vs ramp speed
	pltarr2[0].plot(vs,mean_ETH[0],marker=".",color="green") # plot data
	pltarr2[0].set_ylabel("$s_d(t_f)$",fontsize=22) # label y-axis
	pltarr2[0].set_xlabel("$v/J_{zz}(0)$",fontsize=22) # label x-axis
	pltarr2[0].set_xscale("log") # set log scale on x-axis
	pltarr2[0].grid(True,which='both') # plot grid
	pltarr2[0].tick_params(labelsize=16)
	# subplot 2: entanglement entropy vs ramp speed
	pltarr2[1].plot(vs,mean_ETH[1],marker=".",color="green") # plot data
	pltarr2[1].set_ylabel("$s_\\mathrm{ent}(t_f)$",fontsize=22) # label y-axis
	pltarr2[1].set_xlabel("$v/J_{zz}(0)$",fontsize=22) # label x-axis
	pltarr2[1].set_xscale("log") # set log scale on x-axis
	pltarr2[1].grid(True,which='both') # plot grid
	pltarr2[1].tick_params(labelsize=16)
	# save figure
	fig.savefig('example1_ETH.pdf', bbox_inches='tight')
	#
	#plt.show() # show plots