Out-of-Equilibrium Bose-Fermi Mixtures

This example shows how to code up the Bose-Fermi mixture Hamiltonian:

\[\begin{split}H(t) &=& H_\mathrm{b} + H_\mathrm{f}(t) + H_\mathrm{bf},\nonumber\\ H_\mathrm{b} &=& -J_\mathrm{b}\sum_{j}\left(b^\dagger_{j+1}b_j + \mathrm{h.c.}\right) - \frac{U_\mathrm{bb}}{2}\sum_j n^\mathrm{b}_j + \frac{U_\mathrm{bb}}{2}\sum_j n^\mathrm{b}_jn^\mathrm{b}_j,\nonumber\\ H_\mathrm{f}(t) &=& -J_\mathrm{f}\sum_{j}\left(c^\dagger_{j+1}c_j - c_{j+1}c^\dagger_j\right) + A\cos\Omega t\sum_j (-1)^j n^\mathrm{f}_j + U_\mathrm{ff}\sum_j n^\mathrm{f}_jn^\mathrm{f}_{j+1},\nonumber\\ H_\mathrm{bf} &=& U_\mathrm{bf}\sum_j n^\mathrm{b}_jn^\mathrm{f}_j\end{split}\]

Details about the code below can be found in SciPost Phys. 7, 020 (2019).

Script

download script

 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
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 10                                #
#   In this script we demonstrate how to use QuSpin's                  #	
#   tensor basis class to build Hamiltonians for mixtures of different #
#   species. We use this to study the non-equilibrium dynamics         #
#   in a Bose-Fermi mixture. The example also shows how to compute     #
#   the entanglement entropy shared between the species.               #
########################################################################
from quspin.operators import hamiltonian # Hamiltonians and operators
from quspin.basis import tensor_basis,spinless_fermion_basis_1d,boson_basis_1d # bases
from quspin.tools.measurements import obs_vs_time # calculating dynamics
from quspin.tools.Floquet import Floquet_t_vec # period-spaced time vector
import numpy as np # general math functions
import matplotlib.pyplot as plt # plotting library
#
##### setting up parameters for simulation
# physical parameters
L = 6 # system size
Nf, Nb = L//2, L # number of fermions, bosons
N = Nf + Nb # total number of particles
Jb, Jf = 1.0, 1.0 # boson, fermon hopping strength
Uff, Ubb, Ubf = -2.0, 0.5, 5.0  # bb, ff, bf interaction
# define time-dependent perturbation
A = 2.0
Omega = 1.0
def drive(t,Omega):
	return np.sin(Omega*t)
drive_args=[Omega]
#
###### create the basis
# build the two bases to tensor together to a bose-fermi mixture
basis_b=boson_basis_1d(L,Nb=Nb,sps=3) # boson basis
basis_f=spinless_fermion_basis_1d(L,Nf=Nf) # fermion basis
basis=tensor_basis(basis_b,basis_f) # BFM
#
##### create model
# define site-coupling lists
hop_b = [[-Jb,i,(i+1)%L] for i in range(L)] # b hopping
int_list_bb = [[Ubb/2.0,i,i] for i in range(L)] # bb onsite interaction
int_list_bb_lin = [[-Ubb/2.0,i] for i in range(L)] # bb interaction, linear term
# 
hop_f_right = [[-Jf,i,(i+1)%L] for i in range(L)] # f hopping right
hop_f_left = [[Jf,i,(i+1)%L] for i in range(L)] # f hopping left
int_list_ff = [[Uff,i,(i+1)%L] for i in range(L)] # ff nearest-neighbour interaction
drive_f = [[A*(-1.0)**i,i] for i in range(L)] # density staggered drive
#
int_list_bf = [[Ubf,i,i] for i in range(L)] # bf onsite interaction
# create static lists
static = [	
			["+-|", hop_b], # bosons hop left
			["-+|", hop_b], # bosons hop right
			["n|", int_list_bb_lin], # bb onsite interaction
			["nn|", int_list_bb], # bb onsite interaction
			#
			["|+-", hop_f_left], # fermions hop left
			["|-+", hop_f_right], # fermions hop right
			["|nn", int_list_ff], # ff nn interaction
			#
			["n|n", int_list_bf], # bf onsite interaction
			]
dynamic = [["|n",drive_f,drive,drive_args]] # drive couples to fermions only
#
###### set up Hamiltonian and initial states
no_checks = dict(check_pcon=False,check_symm=False,check_herm=False)
H_BFM = hamiltonian(static,dynamic,basis=basis,**no_checks)
# define initial Fock state through strings
s_f = "".join("1" for i in range(Nf)) + "".join("0" for i in range(L-Nf))
s_b = "".join("1" for i in range(Nb))
# basis.index accepts strings and returns the index which corresponds to that state in the basis list
i_0 = basis.index(s_b,s_f) # find index of product state in basis
psi_0 = np.zeros(basis.Ns) # allocate space for state
psi_0[i_0] = 1.0 # set MB state to be the given product state
print("H-space size: {:d}, initial state: |{:s}>|{:s}>".format(basis.Ns,s_b,s_f))
#
###### time evolve initial state and measure entanglement between species
t=Floquet_t_vec(Omega,10,len_T=10) # t.vals=times, t.i=initial time, t.T=drive period
psi_t = H_BFM.evolve(psi_0,t.i,t.vals,iterate=True)
# measure observable
Sent_args=dict(basis=basis,sub_sys_A="left")
meas = obs_vs_time(psi_t,t.vals,{},Sent_args=Sent_args)
# read off measurements
Entropy_t = meas["Sent_time"]["Sent_A"]
#
######
# configuring plots
plt.plot(t/t.T, Entropy_t)
plt.xlabel("$\\mathrm{driving\\ cycle}$",fontsize=18)
plt.ylabel('$S_\\mathrm{ent}(t)$',fontsize=18)
plt.grid(True)
plt.tick_params(labelsize=16)
plt.tight_layout()
plt.savefig('BFM.pdf', bbox_inches='tight')
#plt.show()
plt.close()