dynamite can be used for not just the obvious spin chain problems, but anything that can be mapped onto a set of spins. Here we will build a model of interacting Majoranas.
There are multiple ways to define a Majorana creation/annihilation operator in a spin basis. In particular, we want to satisfy the anticommutation relation
$$\{ \chi_i, \chi_j \} = 2 \delta_{ij}$$for $i \neq j$. It turns out we can do so with the following mapping:
$$\chi_i = \frac{1}{2} \sigma_{\lfloor i/2 \rfloor}^{x/y} \prod_{k}^{\lfloor i/2 \rfloor - 1} \sigma^z_k$$where that first Pauli matrix is $\sigma^x$ if $i$ is even, and $\sigma^y$ if $i$ is odd.
This basis can be shown fairly easily to satisfy the anticommutation relation we desired. Now let's implement it in dynamite!
In [8]:
from dynamite.operators import sigmax, sigmay, sigmaz, index_product
In [9]:
# product of sigmaz along the spin chain up to index k
k = 4
index_product(sigmaz(), size=k)
Out[9]:
In [10]:
# with that, we can easily build our operator
def majorana(i):
k = i//2
edge_op = sigmay(k) if (i%2) else sigmax(k)
bulk = index_product(sigmaz(), size=k)
return edge_op*bulk
In [11]:
# let's check it out!
majorana(8)
Out[11]:
Looks like exactly what we wanted! We can even check that the anticommutation relation holds:
In [13]:
from dynamite.operators import zero, identity
def anticommutator(a, b):
return a*b + b*a
def check_anticom():
print('i', 'j', 'correct', sep='\t')
print('=======================')
for i in range(3):
for j in range(3):
if i == j:
correct_val = 2*identity()
else:
correct_val = zero()
print(i, j, anticommutator(majorana(i), majorana(j)) == correct_val, sep='\t')
check_anticom()
It was instructive to build it ourselves, but dynamite actually has a Majorana operator built-in, for ease of use. It is the same as ours:
In [14]:
# rename our function, so that we can set majorana to be the dynamite one
my_majorana = majorana
from dynamite.extras import majorana
In [15]:
majorana(8)
Out[15]:
In [16]:
majorana(8) == my_majorana(8)
Out[16]:
We want to build the model
$$H_{\text{SYK}} = \sum_{i,j,k,l} J_{ijkl} \cdot \chi_i \chi_j \chi_k \chi_l$$where the $\chi_i$ represent a Majorana creation/annihilation operator for particle index $i$, and the $J_{ijkl}$ are some random coefficients.
First we must import the things we need:
In [17]:
from dynamite.operators import op_sum, op_product, index_sum
We need to generate all combinations of indices for i,j,k,l, without repeats. Sounds like a task for Python's itertools:
In [18]:
from itertools import combinations
def get_all_indices(n):
'''
Get all combinations of indices i,j,k,l for a system of n Majoranas.
'''
return combinations(range(n), 4)
In [19]:
# does it do what we expect?
for n,idxs in enumerate(get_all_indices(6)):
print(idxs)
if n > 5:
break
print('...')
Looks good! Now let's use that to build the Hamiltonian:
In [20]:
import numpy as np
from numpy.random import seed, normal
# abbreviate
maj = majorana
def syk_hamiltonian(n, random_seed=0):
'''
Build the SYK Hamiltonian for a system of n Majoranas.
'''
# so the norm scales correctly
factor = np.sqrt(6/(n**3))/4
# it's very important to have the same seed on each process if we run in parallel!
# if we don't set the seed, each process will have a different operator!!
seed(random_seed)
return op_sum(factor*normal(-1,1)*maj(i)*maj(j)*maj(k)*maj(l) for i,j,k,l in get_all_indices(n))
Let's try it for a (very) small system!
In [21]:
syk_hamiltonian(5)
Out[21]:
Neat, looks good! Why don't we build it for a bigger system, say 16 Majoranas? (which lives on 8 spins)
In [22]:
H = syk_hamiltonian(16)
Yikes, that was awfully slow for such a small system size. The problem is that the individual Majorana operators are being rebuilt for every term of the sum, and there are a lot of terms. Maybe we can do better by precomputing the Majorana operators. We also use op_product and operator.scale to avoid making unnecessary copies.
In [23]:
def syk_hamiltonian_fast(n, random_seed=0):
'''
Build the SYK Hamiltonian for a system of n Majoranas.
'''
factor = np.sqrt(6/(n**3))/4
seed(random_seed)
majs = [maj(i) for i in range(n)]
return op_sum(op_product(majs[i] for i in idxs).scale(factor*normal(-1,1)) for idxs in get_all_indices(n))
In [24]:
# make sure they agree
assert(syk_hamiltonian(10) == syk_hamiltonian_fast(10))
In [25]:
# check which one is faster!
from timeit import timeit
orig = timeit('syk_hamiltonian(16)', number=1, globals=globals())
fast = timeit('syk_hamiltonian_fast(16)', number=1, globals=globals())
print('syk_hamiltonian: ', orig, 's')
print('syk_hamiltonian_fast:', fast, 's')
That's a huge speedup!
One last thing to note. It may seem odd that we've never actually specified a spin chain length during this whole process. Don't we need to tell dynamite how many spins we need, and thus how big to make our matrices? If the spin chain length is not specified, dynamite just assumes it to extend to the position of the last non-trivial Pauli operator:
In [28]:
m8 = majorana(8)
print('spin chain length:', m8.get_length())
We can use operator.table() to take a look at it:
In [29]:
print(m8.table())
The last non-identity operator is on spin index 4, so a 5-spin chain makes sense.