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
|