The purpose of this tutorial is to teach the interested user to construct bosonic Hamiltonians using QuSpin. To this end, below we focus on the Bose-Hubbard model (BHM) of a 1d chain. The Hamiltonian is $$ H = -J\sum_{j=0}^{L-1}(b^\dagger_{j+1}b_j + \mathrm{h.c.})-\mu\sum_{j=0}^{L-1} n_j + \frac{U}{2}\sum_{j=0}^{L-1}n_j(n_j-1)$$ where $J$ is the hopping matrix element, $\mu$ -- the chemical potential, and $U$ -- the interaction strength. We label the lattice sites by $j=0,\dots,L-1$, and use periodic boundary conditions.
First, we load the required packages:
In [ ]:
from quspin.operators import hamiltonian # Hamiltonians and operators
from quspin.basis import boson_basis_1d # Hilbert space boson basis
import numpy as np # generic math functions
Next, we define the model parameters:
In [2]:
##### define model parameters #####
L=6 # system size
J=1.0 # hopping
U=np.sqrt(2.0) # interaction
mu=2.71 # chemical potential
In order to construct the Hamiltonian of the BHM, we need to construct the bosonic basis. This is done with the help of the constructor boson_basis_1d. The first required argument is the chain length L. As an optional argument one can also specify the number of bosons in the chain Nb. We print the basis using the print() function.
In [3]:
##### construct Bose-Hubbard Hamiltonian #####
# define boson basis with 3 states per site L bosons in the lattice
basis = boson_basis_1d(L,Nb=L) # full boson basis
print(basis)
If needed, we can specify the on-site bosonic Hilbert space dimension, i.e. the number of states per site, using the flag sps=int. This can help study larger systems of they are dilute.
In [4]:
basis = boson_basis_1d(L,Nb=L,sps=3) # particle-conserving basis, 3 states per site
print(basis)
Often times, the model under consideration has underlying symmetries. For instance, translation invariance, parity (reflection symmetry), etc. QuSpin allows the user to construct Hamiltonians in symmetry-reduced subspaces. This is done using optional arguments (flags) passed to the basis constructor.
For instance, if we want to construct the basis in the $k=0$ many-body momentum sector, we do this using the flag kblock=int. This specifies the many-body momentum of the state via $k=2\pi/L\times\texttt{kblock}$.
Whenever symmetries are present, the print() function returns one representative from which one can obtain all 'missing' states by applying the corresponding symmetry operator. It is important to note that, physically, this representative state stands for the linear combination of vectors in the class, not the state that is displayed by print(basis).
In [5]:
basis = boson_basis_1d(L,Nb=L,sps=3,kblock=1) # ... and zero momentum sector
print(basis)
Additionally, the BHM features reflection symmetry around the middle of the chain. This symmetry block-diagonalises the Hamiltonian into two blocks, corresponding to the negative and positive eigenvalue of the parity operator. The corresponding flag is pblock=+1,-1.
In [6]:
basis = boson_basis_1d(L,Nb=L,sps=3,kblock=0,pblock=1) # ... + zero momentum and positive parity
print(basis)
Now that we have constructed the basis in the symmetry-reduced Hilbert space, we can construct the Hamiltonian. It will be hepful to cast it in the fllowing form:
$$H= -J\sum_{j=0}^{L-1}(b^\dagger_{j+1}b_j + \mathrm{h.c.})-\left(\mu+\frac{U}{2}\right)\sum_{j=0}^{L-1} n_j + \frac{U}{2}\sum_{j=0}^{L-1}n_jn_j $$We start by defining the site-coupling lists. Suppose we would like to define the operator $\sum_j \mu_j n_j$. To this, end, we can focus on a single summand first, e.g. $2.71 n_{j=3}$. The information encoded in this operator can be summarised as follows:
In QuSpin, the first two points are grouped together, defininging a list [mu_j,j]=[2.71,3], while the type of operator we specify a bit later (see parantheses). We call this a site-couling list. Summing over multiple sites then results in a nested list of lists:
In [7]:
# define site-coupling lists
hop=[[-J,i,(i+1)%L] for i in range(L)] #PBC
interact=[[0.5*U,i,i] for i in range(L)] # U/2 \sum_j n_j n_j
pot=[[-mu-0.5*U,i] for i in range(L)] # -(\mu + U/2) \sum_j j_n
print(hop)
#print(interact)
#print(pot)
The site coupling lists specify the sites on which the operators act, yet we need to tell QuSpin which operators are to act on these (pairs of) sites. Thus, we need the following operator strings which enter the static and dynamic lists used to define the Hamiltonian. Since the BHM is time-independent, we use an empty dynamic list
In [8]:
# define static and dynamic lists
static=[['+-',hop],['-+',hop],['n',pot],['nn',interact]]
dynamic=[]
print(static)
Building the Hamiltonian with QuSpin is now a one-liner using the hamiltonian constructor
In [9]:
# build Hamiltonian
H=hamiltonian(static,dynamic,basis=basis,dtype=np.float64)
print(H.todense())
when the Hamiltonian is constructed, we see three messages saying that it passes three type of symmetries. QuSpin does checks under the hood on the static and dynamic lists to determine if they satisfy the requested symmetries in the basis. They can be disabled by parsing the following flags to the hamiltonian constructor: check_pcon=False, check_symm=False and check_herm=False.
We can now diagonalise H, and e.g. calculate the entanglement entropy of the ground state.
In [10]:
# calculate eigensystem
E,V=H.eigh()
E_GS,V_GS=H.eigsh(k=2,which='SA',maxiter=1E10) # only GS
print("eigenenergies:", E)
#print("GS energy is %0.3f" %(E_GS[0]))
# calculate entanglement entropy per site of GS
subsystem=[i for i in range(L//2)] # sites contained in subsystem
Sent=basis.ent_entropy(V[:,0],sub_sys_A=subsystem,density=True)['Sent_A']
print("GS entanglement per site is %0.3f" %(Sent))
In [11]:
psi_k=V[:,0]
psi_Fock=basis.get_vec(psi_k)
print(psi_k.shape, psi_Fock.shape)
In [ ]: