Level statistics in many-body localized systems
Author Hannes Karlsson
The Hamiltonian for the Ising model with an applied magnetic field in \(z\)-direction $$ H = -J\sum_i \sigma^z_i \sigma^z_{i+1} - h\sum_i \sigma_i^z $$ commutes with all operators \(\sigma_i^z\), which constitutes an extensive set of conserved quantitues, and so the system is integrable. If one however adds a magnetic field in the \(x\)-direction $$ H = -J\sum_i \sigma^z_i \sigma^z_{i+1} - \sum_i (h\sigma_i^z + \gamma\sigma_i^x), $$ the model becomes ergodic. It is thus expected to follow the predictions of random matrix theory and the eigenstate thermalization hypothesis. Consider now the same Hamiltonian but with random interaction and field strengths. $$ H = \sum_i J_i \sigma^z_i \sigma^z_{i+1} - \sum_i (h_i\sigma_i^z + \gamma_i\sigma_i^x), $$ where, for the sake of concreteness, we take \(J_i\in [0.8,1.2]\), \(h_i\in [-W,W]\), \(\gamma_i=1\), and \(W\) is the disorder strength. This model is known to display many-body localisation for large enough \(W\), i.e, it hosts an extensive set of emergent (local) conserved quantities, often called l-bits. To study this, we compute the inverse participation ration, IPR\(=\sum_i |\psi_i|^4\), most often used to study Anderson localisation, for some of the low lying eigenstates. In a perfectly localised state, one would expect IPR\(=1\), while for a perfecly delocalised state, IPR\(=1/N\) with \(N\) being the number of basis states. IPR is however not a perfect measure of localisation for an interacting system, so we will also consider the level spacing ratio, \(r=min(\delta_n,\delta_{n+1})/max(\delta_n,\delta_{n+1})\), where \(\delta_n=E_{n+1}-E_n\). For an ergodic phase, we would observe a Gaussian (or Wigner-Dyson) distribution of the eigenvalues, with \(\langle r \rangle \approx 0.530\), while in the MBL phase (or regime), we would observe a Poisson distribution, with \(\langle r \rangle \approx 0.386\).
As we can see in the figure above, with \(L=10\) spins, the IPR is growing with the disorder \(W\), for all the observed eigenstates. We also see that the level spacing ratio goes from close to the value of a Poisson distribution for \(W=0\), then goes up to the value of the Gaussian distribution, where it reaches it maximal ergodicity, before going down again to a Gaussian distribution when the l-bits are formed.
#include <xdiag/all.hpp>
using namespace xdiag;
double level_spacing_ratio(arma::vec eigs) {
arma::vec diff = arma::vec(eigs.size() - 1);
for (int i = 0; i < eigs.size() - 1; i++) {
diff[i] = eigs[i + 1] - eigs[i];
}
arma::vec r = arma::vec(diff.size() - 1);
for (int i = 0; i < diff.size() - 1; i++) {
r[i] = std::min(diff[i], diff[i + 1]) / std::max(diff[i], diff[i + 1]);
}
return arma::sum(r) / r.size();
}
double ipr_func(arma::vec state) {
return arma::sum(arma::pow(arma::abs(state), 4));
}
int main() {
int L = 10;
int nres = 2;
int neigs = 4;
arma::vec Ws = arma::linspace(0, 5, 100);
arma::mat rs = arma::mat(nres, Ws.size());
arma::cube iprs = arma::cube(nres, neigs, Ws.size());
for (int i = 0; i < nres; i++) {
arma::vec J = arma::randu(L - 1, arma::distr_param(0.8, 1.2));
for (int j = 0; j < Ws.size(); j++) {
arma::vec h = arma::randu(L, arma::distr_param(-Ws(j), Ws(j)));
// construct Hamiltonian and Hilbert space
auto block = Spinhalf(L);
auto ops = OpSum();
for (int l = 0; l < L - 1; l++) {
ops += J(l) * Op("SzSz", {l, l + 1});
}
for (int l = 0; l < L; l++) {
ops += h(l) * Op("Sz", l);
ops += 0.5 * Op("S+", l);
ops += 0.5 * Op("S-", l);
}
arma::mat H = matrix(ops, block);
arma::mat eigvecs;
arma::vec eigvals;
arma::eig_sym(eigvals, eigvecs, H);
rs(i, j) = level_spacing_ratio(eigvals);
for (int k = 0; k < neigs; k++) {
arma::vec eigvec = eigvecs.col(k);
iprs(i, j, k) = ipr_func(eigvec);
}
}
}
arma::vec r = arma::sum(rs, 0);
arma::mat ipr = arma::sum(iprs, 0);
// save and plot
return 0;
}
using XDiag
using CairoMakie
using LinearAlgebra
using Statistics
function diagonalise(L,W,J,h)
# construct Hilbert space
hspace = Spinhalf(L)
# construct Hamiltonian
ops = OpSum()
for i in 1:L-1
ops += J[i] * Op("SzSz",[i,i+1])
end
for i in 1:L
ops += h[i] * Op("Sz",i)
ops += 0.5 * Op("S+",i)
ops += 0.5 * Op("S-",i)
end
# Lanczos
# res = eigs_lanczos(ops,hspace)
H = matrix(ops, hspace)
eig = eigen(Hermitian(H))
return eig.values, eig.vectors
end
function level_spacing_ratio(eigs)
diff = []
for i in 1:length(eigs)-1
push!(diff,eigs[i+1] - eigs[i])
end
r = []
for i in 1:length(diff)-1
push!(r,min(diff[i],diff[i+1])/max(diff[i],diff[i+1]))
end
return mean(r)
end
function ipr_func(state)
ipr = 0
for i in 1:length(state)
ipr += abs(state[i])^4
end
return ipr
end
function main()
# set model parameters
L = 10
Ws = range(0,step=0.05,length=100)
nres = 20
rs = Array{Float64}(undef,nres,length(Ws))
f = Figure()
neigs = 4 # number of low lying eigenstates to check IPR for
iprs = zeros(nres,neigs,length(Ws))
for ii in 1:nres
J = (rand(Float64,L-1) .* 0.4) .+ 0.8 # J between 0.8 and 1.2
for (jj,W) in enumerate(Ws)
println("Now computing W=$W")
h = (rand(Float64,L) .* 2*W) .- W # h between -W and W
evals, evecs = diagonalise(L,W,J,h)
rs[ii,jj] = level_spacing_ratio(evals)
for k in 1:neigs
evec = evecs[:,k]
iprs[ii,k,jj] = ipr_func(evec)
end
end
end
r = dropdims(mean(rs,dims=1),dims=1)
ipr = dropdims(mean(iprs,dims=1),dims=1)
ax1 = Axis(f[1,1],
xlabel=L"$W$",
ylabel="IPR",
)
for k in 1:neigs
scatter!(ax1,Ws,ipr[k,:])
end
ax2 = Axis(f[1,2],
xlabel=L"$W$",
ylabel=L"$\langle r \rangle$",
)
scatter!(ax2,Ws,r)
hlines!(ax2,0.386,label="Poisson",color=:red)
hlines!(ax2,0.530,label="Gaussian",color=:green)
axislegend(ax2)
save("../plots/mbl/L($L).png",f)
end
main()