Welcome to QuantumChaos’s documentation!

This project contains samples of the code used to write the following articles:

  • Guy Zisling, Lea F. Santos, and Yevgeny Bar Lev, How many particles make up a chaotic many-body quantum system? SciPost Phys. 10, 88 (2021).

  • Guy Zisling, Dante M. Kennes, and Yevegeny Bar Lev, Transport in Stark Many Body Localized Systems arXiv:2109.06196.

The code can be used for two main purposes:

  1. Build Hamiltonian’s (matrix representations) of spin-chain systems (mainly “XXZ” chain).

  2. Examine the statical and dynamical properties of those Hamiltonian’s.

Link to the github repo:

Written by Guy Zisling.

Indices and tables

Math Intro

The “XXZ” chain can be written in terms of spin matrices using

\[\hat{H}=\sum_{j=1}^{L-1}\frac{J_{xy}}{2}\left(\hat{S}_{j}^{+}\hat{S}_{j+1}^{-} + \hat{S}_{j}^{-}\hat{S}_{j+1}^{+}\right)+J_{z}\hat{S}_{j}^{z}\hat{S}_{j+1}^{z},\]
where \(\hat{S}^{\{x,y,z\}}\) are the Pauli matrices times a factor of \(\frac{1}{2}\).

Usually, this type of systems can be classified into two distinct classes:

  1. Many-body localized systems (MBL).

  2. Quantum Chaotic systems (thermalizing systems).

Some key differences are:

Feature

Quantum Chaotic

Many-body Localized

Memory of initial conditions

None

Some

Eigenvalues

Correlated

Uncorrelated

Off-diagonal elements distribution

Normal

Sharp distribution

Transport

Super-diffusive

Sub-diffusive

Entanglement entropy spreading

Power-law

Logarithmic

  • Static properties are usually examined by using the Hamiltonian eigen-basis, where \(E_\alpha = \left\langle {\phi_\alpha}\right .\left|{\hat{H}}\right|\left .{\phi_\alpha}\right \rangle\).

  • Dynamical properties are examined using the time evolution operator \(\hat{U}(t)=e^{-i\hat{H}t}\).

Building Hamiltonians

This module uses to build matrix representations of Hamiltonians.

The main logic behind the code is to represent the wavefunctions at an easy computational basis (and code it as binary numbers):

\[\left|\psi\right\rangle=\left|\uparrow\downarrow\downarrow\right\rangle=\texttt{[1,0,0]}.\]
We start building the Hamiltonian matrix by generating this basis (either in basis of 10 or binary numbers array). If we take for example three spins (\(L=3\)) and one excitation (\(N=1\)), we get the following basis:
\[\left\{\left|\downarrow\downarrow\uparrow\right\rangle , \left|\downarrow\uparrow\downarrow\right\rangle , \left|\uparrow\downarrow\downarrow\right\rangle\right\},\space \left\{\texttt{[1,0,0]}, \texttt{[0,1,0]}, \texttt{[0,0,1]} \right\} \]

Now we can see how a term in the Hamiltonian operates on this wavefunction

\[\hat{S}^-_0\hat{S}^+_1\left|\uparrow\downarrow\downarrow\right\rangle=\left|\downarrow\uparrow\downarrow\right\rangle\]
Then we could match that output with the origin \(\texttt{[1,0,0]}\leftrightarrow\texttt{[0,1,0]}\), so the matrix representation of this term in this basis is:

\[\begin{split}\hat{S}_{0}^{-}\hat{S}_{1}^{+}=\begin{array}{c|ccc} \psi & \left|\downarrow\downarrow\uparrow\right\rangle & \left|\downarrow\uparrow\downarrow\right\rangle & \left|\uparrow\downarrow\downarrow\right\rangle \\ \hline \left\langle \downarrow\downarrow\uparrow\right| & 0 & 0 & 0\\ \left\langle \downarrow\uparrow\downarrow\right| & 0 & 0 & 1\\ \left\langle \uparrow\downarrow\downarrow\right| & 0 & 0 & 0 \end{array}=\left(\begin{array}{ccc} 0 & 0 & 0\\ 0 & 0 & 1\\ 0 & 0 & 0 \end{array}\right)\end{split}\]

All matrices generated are at the size of the basis.

build_hamiltonian.matt0(n, Jx, Jz, basis, bc)

Generates the XXZ Hamiltonian

Parameters

  • n (int) – number of spins

  • Jx (float) – xx interaction strength

  • Jz (float) – zz interaction strength

  • basis (np.array) – array of bits encoding the wave functions

  • bc (0 or 1) – boundary conditions 0 = open, 1 = close

Returns

xxz chain

Return type

sparse matrix

build_hamiltonian.matt0sz(i, basis)

Generates the operator \(\hat{S}^z_i\)

Parameters

  • i (int) – location of the operator

  • basis (np.array) – array of wave functions, encoded in bits formation

Returns

\(\hat{S}^z_i\)

Return type

sparse matrix

build_hamiltonian.matt_add_imp(H0, imp, h, basis, n=0)

Add a magentic impurity at a specific site of the chain \(h\hat{S}^z_{\textrm{imp}}\)

Parameters

  • H0 (sparse matrix) – existing Hamiltonian

  • imp (int) – impurity location

  • h (float) – impurity strength

  • basis (np.array) – array of wave functions, encoded in bits formation

Returns

The initial Hamiltonian plus the impurity term

Return type

sparse matrix

build_hamiltonian.matt_add_jz(H0, n, Jz, basis, bc)

Adds \(\hat{S}^z_i\hat{S}^z_{i+1}\) interaction for an existing Hamiltonian

Parameters

  • H0 (sparse matrix) – existing Hamiltonian

  • n (int) – number of spins in the chain

  • Jz (float) – interaction strength

  • basis (np.array) – array of numbers encoding the wave functions

  • bc (0 or 1) – boundary conditions 0 = open, 1 = close

Returns

The matrix \(\hat{H}_0\) with addition of \(zz\) interactions

Return type

sparse matrix

build_hamiltonian.matt_add_stark(H0, n, f, a, basis, bc)

Add a linear field of strength f to an existing Hamiltonian (stark potential).

\[\hat{H}_0\rightarrow\hat{H}_0 + \sum_{j=1}^{L}\left(\gamma j+\alpha j^{2}/L^{2}\right) \hat{S}_{j}^{z}\]

Parameters

  • H0 (sparse matrix) – existing Hamiltonian

  • n (int) – number of spins in the chain

  • f (float) – linear potential strength

  • a (float) – curvature strength

  • basis (np.array) – array of numbers encoding the wave functions

  • bc (0 or 1) – boundary conditions 0 = open, 1 = close

Returns

The initial Hamiltonian plus the Stark potential term

Return type

sparse matrix

Utils

For utils used in more than one module.

The transformation between wave-function bits and integer representation is defined as following:

\[\vec{x}\rightarrow\sum_{i=L}^1(1-x_i)\cdot 2^{L-i}=\boldsymbol{\mathrm{int}}.\]

See legacy_utils for more information about this transformation.

utils.bittoint(basis)

Basis in bits to integers (e.g \(\texttt{[0,1,1]}\rightarrow\boldsymbol{4}\))

Parameters

basis (np.array) – a list of bits that represents the spins (1 = spin pointing up, 0 = down)

Returns

an array of integers encoding the wave function

Return type

np.array

utils.blockex(n, j)

Generates a basis of n spins and j excitations (spins pointing up)

Parameters

  • n (int) – number of spins in the chain

  • j (int) – number of excitations in the basis

Returns

wave functions basis in bits

Return type

np.array

utils.ordtobit(ord, n)

Integers basis to bits basis (e.g \(\boldsymbol{4}\rightarrow\texttt{[0,1,1]}\))

Parameters

  • ord (np.array) – an array of integers that represents the spins

  • n – number of spins in the chain

Returns

a list of bits that represents the spins (1 = spin pointing up, 0 = down)

Return type

np.array

utils.ptime(tz)
Parameters

tz (float) – a time.time() output

Returns

a string with the time it took since the initial counter tz

Dynamical Tools

The cheapest way to calculate the time evolution of some operator \(\hat{O}\) is to use the ‘Dynamical Typicality’ concept (see Ref):

\[\frac{1}{\mathcal{D}}\text{Tr}\hat{O}\approx\left\langle{\psi}\right .\left|{\hat{O}}\right| \left .{\psi}\right \rangle,\]
and to use the Krylov time-space evolution, explained in the following Ref and built in scipy.sparse.linalg.expm_multiply.

briefly, instead of diagonalizing the full matrix \(\hat{H}\) in order to calculate the time evolution operator \(\hat{U}(t)=e^{-i \hat{H} t}\), one asses the time evolution by averaging over random states \(\left| \psi (t) \right\rangle = \hat{U}(t)\left|\psi(0)\right\rangle\), without explicitly calculating \(\hat{U}(t)\).

dynamical_tools.entropy(v, n, basis)

Calculates the entanglement entropy of a vector v

Parameters

  • v (np.array) – vector

  • n (int) – number of spins in the chain

  • basis (np.array) – array of numbers encoding the wave functions

Returns

entanglment entropy

Return type

(float)

dynamical_tools.entropy_vs_time(H, n, basis, t, k, dt)

Calculate the entanglement entropy

Parameters

  • H (sparse matrix) – Hamiltonian

  • n (int) – number of spins in the chain

  • basis (np.array) – array of numbers encoding the wave functions

  • t (float) – ending time of the calculation

  • k (int) – number of point to sample in time

  • dt (float) – time intervals for the Krylov time evolution

Returns

[time array, entangelment entropy array]

Return type

(np.array)

dynamical_tools.evol_kryl(H, v, dt, T)

Calculates the time evolution of a vector using Krylov sub-space

\[\vec{v}(t) = e^{-i\hat{H}t}\vec{v}\]

Parameters

  • H (sparse matrix) – Hamiltonian

  • v (np.array) – starting vector

  • dt (float) – time interval

  • T (int) – number of steps

Returns

vector

Return type

(np.array)

dynamical_tools.msd(H, n, basis, t, k, dt, seed=False)

Calculates the spins mean square displacement after a spin flip in the middle of the chain

\[G_{n}\left(t\right)=\frac{1}{\mathcal{D}}\textrm{Tr}\left[\hat{S}_{n}^{z}\left(t\right)\hat{S}_{L/2}^{z}\right]\]
\[x^{2}\left(t\right)=\sum_{n}n^{2}\left(G_{n}\left(t\right)-G_{n}\left(0\right)\right)\]

Parameters

  • H (sparse matrix) – Hamiltonian

  • n (int) – number of spins in the chain

  • basis (np.array) – array of numbers encoding the wave functions

  • t (float) – ending time of the calculation

  • k (int) – number of point to sample in time

  • dt (float) – time intervals for the Krylov time evolution

  • seed (bool) – fix the seed

Returns

time array, msd array, Gr profile array

Return type

(np.array, np.array, np.array)

Static Tools

Static properties are usually examined by using the Hamiltonian eigen-basis, where

static_tools.diag_elements(H0, n, basis)

ETH diagonal elements test, used to plot \(\left\langle {\phi_\alpha}\right .\left|{\hat{O}}\right|\left .{\phi_\alpha}\right \rangle\) as function of the energy \(E_\alpha\).

Parameters

  • H0 (sparse matrix) – the Hamiltonian

  • n (int) – number of spins

  • basis (np.array) – array of numbers encoding the wave functions

Returns

normalized energies (from 0 to 1), corresponding diagonal element

Return type

(np.array, np.array)

static_tools.offdiag(H0, n, basis, dw=0.05, de=0.05, is_unfold=False)

ETH for off-diagonal elements test, used to identify the chaoticity of an Hamiltonian as a function of energy. Used to plot how close to gaussian dist are the off-diagonal elements \(\hat{O}_{\alpha\beta}={\left.\left\langle {\phi_\alpha}\right .\left|{\hat{O}}\right|\left .{\phi_\beta}\right \rangle\right.}_{\alpha\neq\beta}\) as a function of the frequency \(\omega\equiv E_{\beta} - E_{\alpha}\).

The gaussianity test is defined as \(\Gamma_{\hat{O}}\left(\omega\right)=\frac{\overline{\left|O_{ \alpha\beta}\right|^{2}}}{\overline{\left|O_{\alpha\beta}\right|}^{2}}\) which is \(\pi/2\) for Gaussian distribution.

Parameters

  • H0 – Hamiltonian or eigenvalues and eigenvectors

  • n (int) – number of spins

  • basis (np.array) – array of wave functions, encoded in bits formation

  • dw (float) – size of frequency bins

  • de (float) – partial part of the spectrum that is being examined (from 0 to 1)

  • is_unfold (bool) – take unfolded eigenvalues

Returns

\(\omega\), \(\Gamma_{\hat{O}}\left(\omega\right)\)

Return type

(np.array, np.array)

static_tools.offdiag_dist(H0, n, basis, e_number=200, bin_num=200, normed=False)

ETH for off-diagonal elements test, used to plot the distribution of the off-diagonal elements of some local observable. return the histogram of \(\hat{O}_{\alpha\beta}={\left.\left\langle {\phi_\alpha}\right .\left|{\hat{ O}}\right|\left .{\phi_\beta}\right \rangle\right.}_{\alpha\neq\beta}\) built from e_number of eigenstates near the middle of the spectrum.

Parameters

  • H0 (sparse matrix) – the Hamiltonian

  • n (int) – number of spins

  • basis (np.array) – array of wave functions, encoded in bits formation

  • e_number (int) – number of eigenvalues used to calculate the

  • bin_num (int) – number of bins for the histogram

  • normed (bool) – normalize by the system size and Hilbert dimension

Returns

an array with the histogram x,y data (same number of points at each axis)

Return type

(np.array)

static_tools.offdiag_stats(H0, n, basis, e_number=200)

ETH for off-diagonal elements test, used to plot the distribution of the off-diagonal elements of some local observable. return the variance and kurtosis of \(\hat{O}_{\alpha\beta}={\left.\left\langle {\phi_\alpha}\right .\left|{\hat{ O}}\right|\left .{\phi_\beta}\right \rangle\right.}_{\alpha\neq\beta}\) built from e_number of eigenstates near the middle of the spectrum.

Parameters

  • H0 (sparse matrix) – the Hamiltonian

  • n (int) – number of spins

  • basis (np.array) – array of numbers encoding the wave functions

  • e_number (int) – number of eigenvalues used to calculate the

Returns

variance, kurtosis

Return type

(float, float)

static_tools.rfold(H)

A chaos matrix the arises from the eigenvalues spacing statstics

Defining the level spacing \(s_{\alpha}=E_{\alpha+1}-E_{\alpha}\) one can define the ‘r-metric’ by

\[r_{\alpha}=\min\left(\frac{s_{\alpha}}{s_{\alpha-1}},\frac{s_{\alpha-1}}{s_{\alpha}}\right),\]
where
\[\begin{split}\left\langle r\right\rangle \approx\begin{cases} 0.39 & \text{Poisson dist.}\\ 0.536 & \text{Wigner Dyson dist.} \end{cases}\end{split}\]

Parameters

H (np.array) – matrix or eigenvalues

Returns

\(\langle r \rangle\)

Return type

float

static_tools.unfold(H, cut=2000, normed=False)

Returns the unfolded eigenvalues (used to remove the bais of the denisity of states in the level spacing statistics).

Parameters

  • H (np.array) – a matrix or the eigenvalues of it

  • cut (int) – how many eigenvales to cut

  • normed (bool) – normalize the density of staes

Returns

unfolded eigenvalues

Return type

np.array

Legacy Utils

Old versions of some utils, that have been replaced with a more efficient src_code. This module is used for benchmarking, and educational purposes.

The main idea behind this module is that matrices are generating using Kronecker product (denoted by \(\otimes\)), this is much less efficient but a bit easier to grasp than the methods used in build_hamiltonian module.

For example, a term can be built as follow:

\[\hat{S}^z_i = \mathbb{I}_{2^i}\otimes \hat{S}^z \otimes \mathbb{I}_{2^{L-i-1}},\]
where \(L\) is the chain length, and \(\mathbb{I}_N\) is the identity matrix of dim \(N\).

Using this method the output is the matrix representation on the ‘computational basis’, for example for two spins (\(L=2\)):

\[\{\left|\uparrow\uparrow\right\rangle, \left|\uparrow\downarrow\right\rangle, \left|\downarrow\uparrow\right\rangle, \left|\downarrow\downarrow\right\rangle\}.\]

If we want to use only a sub-sector of the \(N=2^L\) matrix we need to code this basis into ordered numbers, to do so we use the following transformation

\[\left|\downarrow\uparrow\right\rangle=\texttt{[0,1]}=\sum_{i=L}^1(1-x_i)\cdot 2^{L-i}=0\cdot2^0+1\cdot2^1=\boldsymbol{2}.\]

For example if we want to take only the zero magnetization sector (where the same number of spins point up and down) out of the following matrix

\[\begin{split}\hat{S}_{1}^{z}=\mathbb{I}_{2}\otimes\hat{S}^{z}=\left(\begin{array}{cccc} 0.5 & 0 & 0 & 0\\ 0 & -0.5 & 0 & 0\\ 0 & 0 & 0.5 & 0\\ 0 & 0 & 0 & -0.5 \end{array}\right)\end{split}\]

We need to cut the following entries

\[\{\left|\uparrow\downarrow\right\rangle, \left|\downarrow\uparrow\right\rangle\}=\{\boldsymbol{1}, \boldsymbol{2}\},\]

so the matrix representation in this basis would be

\[\begin{split}\hat{S}_{1,\mathrm{basis}}^{z}=\left(\begin{array}{cc} -0.5 & 0\\ 0 & 0.5 \end{array}\right)\end{split}\]

legacy_utils.outerr(a, i, n)

Generate an \(2^n\) dim matrix representing matrix a in place i in the full hilbert space

Parameters

  • a (np.array) – matrix

  • i (int) – location in the spin chain

  • n (int) – spin chain length

Returns

\(\hat{a}_i\) \(2^n imes2^n\) matrix repreesntation

Return type

np.array

legacy_utils.xxz_block_add_imp(H, n, imp, h, basis)

Add impurity to existing Hamiltonian \(\hat{H}\rightarrow \hat{H} + h\cdot\hat{S}^z_{\mathrm{imp}}\)

Parameters

  • H (np.array) – initial Hamiltonian

  • n (int) – number of spins

  • imp (int) – impurity location

  • h (float) – impurity strength

  • basis (np.array) – an array encoding the wave functions

Returns

The Hamiltonian matrix H plus the impurity term.

Return type

np.array

legacy_utils.xxz_block_add_stark(H, n, f, a, basis)

Add Stark potential to existing Hamiltonian

Parameters

  • H (np.array) – initial Hamiltonian

  • n (int) – number of spins

  • f (float) – linear potential strength

  • a (float) – curvature strength

  • basis (np.array) – an array encoding the wave functions

Returns

The Hamiltonian matrix H plus the Stark potential

Return type

np.array

legacy_utils.xxz_block_stark(n, Jx, Jz, f, a, basis, bc)

Generates the Stark Hamiltonian (xxz chain with linear potential)

Parameters

  • n (int) – number of spins

  • Jx (float) – xx interaction strength

  • Jz (float) – zz interaction strength

  • f (float) – linear potential strength

  • a (float) – curvature strength

  • basis (np.array) – an array encoding the wave functions

  • bc (int) – boundary conditions (0 - open, 1 - closed)

Returns

the Stark Hamiltonian matrix representation on the basis given

Return type

np.array

Unit Test

Test module, to use it run python -m unittest -v src_code.unit_test.Test

class unit_test.Test(methodName='runTest')
classmethod setUpClass()

Initialize variables, creates XXZ matrices for benchmarking, diagonalize the matrices.

test_diag_elements()

Make sure diag_elements runs properly and returns a tuple

test_impurity()

Benchmark the addition of impurity

test_msd()

Make sure msd functions runs properly and return a tuple

test_offdiag()

Make sure offdiag runs properly and returns a tuple

test_offdiag_dist()

Make sure offdiag_dist runs properly and returns an np.array

test_rfold()

Make sure rfold is runs properly and returns a float

test_stark()

Benchmark the Stark potentials

test_xxz()

Benchmark the basic XXZ matrices