Load up the sympy python package to do symbolic mathematics
In [1]:
from IPython.display import display
from sympy.interactive import printing
printing.init_printing(use_latex=True)
from __future__ import division
import sympy as sym
from sympy import *
import matplotlib
matplotlib.rc('font',**{'family':'serif'})
We start by defining the function that calculates the truncated hamiltonian matrix. I call this Hamiltonian, but it is really the set of linear equations for the coefficients that describe the wave functions in terms of plane waves. We use the sympy package to define the matrix in terms of symbolic quantities. This allows directly getting the latex output for the matrix, which comes in handy when including it in a document.
In [2]:
def H0_sin( Npoints, q, v0 ):
def f(i,j):
if i==j:
return 4*(q+(i-Npoints))**2 + v0/2
elif i-j == 1:
return -v0/4
elif j-i == 1:
return -v0/4
else:
return 0.
mdim = (2*Npoints+1)
matrix = sym.Matrix(mdim,mdim,f)
return matrix
q, v0 = sym.symbols("q V_{0}")
H0_sin(2,q,v0)
Out[2]:
With this matrix in hand we can diagonalize it and get the eigenvalues, the eigenvalues of this matrix correspond to the energy for the given quasimomentum, lattice depth, and band index. We define functions to diagonalize the matrix. At the same time we define functions that calculate the energy of the harmonic oscillator states in each lattice site, when the siteis considered as an independent harmonic oscillator potential.
In [3]:
def eigenvalsH0( q, v0, NPoints=5):
r"""Calculate energy in 1D lattice
Given the quasimomentum and the lattice depth, this function
will calculate the set of energy eigenvalues for a number of
bands determined by NPoints. The calculation is for a 1D
lattice.
Parameters
----------
q : float
The value of quasimomentum. Units are 2\pi, so q can take
values in (-0.5,0.5]
v0 : float
The value of the lattice depth in recoils
NPoints : float, optional
Plane waves e^{imGx} will be considered in the expansion of
the Bloch states such that their momentum satisfies
|m|<NPoints. Here G is the unit reciprocal lattice vector
of the 1D lattice: G = 2*pi / a
A total of 2*NPoints + 1 plane waves is used.
Returns
-------
array
An array with the sorted energy eigenvalues. When they are
sorted the index of the array corresponds to the band index.
Other Parameters
----------------
Raises
------
See Also
--------
bandH0
Notes
-----
References
----------
For more information see [1]_.
.. [1] O. Morsch and M. Oberthaler, “Dynamics of Bose-Einstein
condensates in optical lattices,” Rev. Mod. Phys. 78, 179–215
(2006).
Examples
--------
>>> q=0.
>>> v0=7.0
>>> print eigenvalsH0(q, v0)
"""
q_Sym, v0_Sym = sym.symbols("q V_{0}")
h0 = np.array( np.array( H0_sin(NPoints,q_Sym,v0_Sym).subs( [(q_Sym,q), (v0_Sym,v0)] )), np.float )
return np.sort( np.linalg.eigvalsh(h0) )
def bandH0( NBand, qs, v0, NPoints=5):
r"""Version of eigenvalsH0 that works with an array input for q"""
Eq = np.array( [ eigenvalsH0( q, v0, NPoints=NPoints)[NBand] for q in qs] )
return Eq
def hosc( NBand, v0 ):
r"""Energy of the NBand^th h.o. state in a lattice with depth v0"""
return 2* np.sqrt( v0 ) * (NBand + 0.5)
We then make a plot of the band structure, energy as a function of quasimomentum
In [110]:
# Start matplotlib
from matplotlib import rc
rc('text', usetex=True)
plt.rcParams['text.latex.preamble'] = [
r'\usepackage{bm}', # for bold math
r'\usepackage{siunitx}', # i need upright \micro symbols, but you need...
r'\sisetup{detect-all}', # ...this to force siunitx to actually use your fonts
#r'\usepackage[textlf]{MyriadPro}',
#r'\usepackage{sansmath}', # load up the sansmath so that math -> helvet
#r'\sansmath', # <- tricky! -- gotta actually tell tex to use!
]
plt.rcParams['axes.linewidth'] = 0.6
plt.rcParams['patch.linewidth'] = 0.4
depths = [0., 4., 8., 12., 16.]
qs = np.linspace(-1,1,40)
fs = 12
fig = plt.figure(figsize=(2.16*len(depths),3.6))
gs = matplotlib.gridspec.GridSpec( 1,len(depths))
for i,d in enumerate(depths):
ax = fig.add_subplot( gs[0,i])
ax.set_xlabel("$q / (2\pi/a)$",fontsize=fs)
ax.set_ylabel("$E_{q} / E_{r} $",fontsize=fs)
ax.set_title('$%.1f E_{r}$'%d,fontsize=fs)
ax.set_xlim(-0.5,0.5)
ax.set_ylim(0.,16.)
ax.fill_between(qs, 0., d, facecolor='gray', alpha=0.3)
ax.plot( qs, hosc(0,d)*np.ones_like(qs), 'b:' )
qs = np.linspace(-0.5,0.5,101)
ax.plot( qs, bandH0( 0, qs, d, NPoints=5))
ax.plot( qs, bandH0( 1, qs, d, NPoints=5))
ax.plot( qs, bandH0( 2, qs, d, NPoints=5))
gs.tight_layout(fig, rect=[0.,0.,1.,1.0])
fig.savefig('BandStructure_figures/bands1d.png',dpi=300)
fig.savefig('BandStructure_figures/bands1d.pdf')
We define a way to calculate the effective mass of the system, to do this we fit the behavior of $E(q)$ near $q=0$
In [7]:
depth = 40.
qs = np.linspace(-0.2,0.2,40)
fig = plt.figure(figsize=(5,5))
gs = matplotlib.gridspec.GridSpec( 1,1)
ax = fig.add_subplot( gs[0,0])
ax.set_xlabel("$q / (2\pi/a)$",fontsize=16)
ax.set_ylabel("$E_{q} / E_{R} $",fontsize=16)
ax.set_title('$%.1f E_{R}$'%depth,fontsize=16)
#ax.set_xlim(-0.1,0.1)
#ax.set_ylim(1.,2.)
es = bandH0( 0, qs, depth, NPoints=5)
ax.plot( qs, es)
z = np.polyfit( qs, es, 2)
print "Effective mass at %.3g Er = %.3g m " %( depth, 4./z[0])
p = np.poly1d(z)
ax.plot( qs, p(qs))
gs.tight_layout(fig, rect=[0.,0.,1.,1.0])
#fig.savefig('BandStructure_figures/bands1d.png',dpi=240)
#fig.savefig('BandStructure_figures/bands1d.pdf')
#plt.close()
We save some interpolation data to easily access the effective mass as a function of lattice depth
In [5]:
def effective_mass( v0):
qs = np.linspace( -0.2, 0.2, 40 )
es = bandH0( 0, qs, v0, NPoints=5)
z = np.polyfit( qs, es, 2 )
return 4./z[0]
v0s = np.linspace(0.001, 20., 120)
ms = np.empty_like(v0s)
for ii,v0 in enumerate(v0s):
ms[ii] = effective_mass(v0)
np.savetxt('banddat/interpdat_effective_mass.dat', np.column_stack((v0s,ms)))
We define a function that allows us to calculate the band width for a certain lattice depth
In [4]:
def bands1d( v0, NBand=0 ):
r""" Returns an array of shape (2,) which has the bottom and top of the band energies
Parameters
----------
v0 : float
Depth of the lattice in recoils
NBand : optional
Band index (default = 0 )
"""
qbot = 0.
qtop = 0.5
Eq0 = eigenvalsH0( qbot, v0)[NBand]
Eq1 = eigenvalsH0( qtop, v0)[NBand]
if Eq0 > Eq1:
return np.array([Eq1,Eq0])
else:
return np.array([Eq0,Eq1])
def bands1dvec( v0, NBand=0):
r""" Returns an array of shape (len(v0,2) which has the bottom and top of the band energies
Parameters
----------
v0 : array-like
Array of lattice depths. Depths are in recoils
NBand : optional
Band index (default = 0 )
"""
band = numpy.empty( (len(v0), 2) )
for i,v in enumerate(v0):
band[i] = bands1d(v, NBand)
return band
And we proceed to make a plot of the lattice band structure as a function of lattice depth
In [108]:
# Start matplotlib
from matplotlib import rc
rc('text', usetex=True)
plt.rcParams['text.latex.preamble'] = [
r'\usepackage{bm}', # for bold math
r'\usepackage{siunitx}', # i need upright \micro symbols, but you need...
r'\sisetup{detect-all}', # ...this to force siunitx to actually use your fonts
#r'\usepackage[textlf]{MyriadPro}',
#r'\usepackage{sansmath}', # load up the sansmath so that math -> helvet
#r'\sansmath', # <- tricky! -- gotta actually tell tex to use!
]
plt.rcParams['axes.linewidth'] = 0.6
plt.rcParams['patch.linewidth'] = 0.4
fs = 12
fig = plt.figure(figsize=(4.4,2.6))
gs = matplotlib.gridspec.GridSpec( 1,1, left=0.11, right=0.74, bottom=0.16, top=0.96)
ax0 = fig.add_subplot( gs[0,0])
for ax in [ax0]:
ax.set_xlabel("$V_{0} (E_{R})$",fontsize=fs)
ax.grid(alpha=0.5)
ax.set_ylim(0.,15.)
ax.set_xlim(0.,20.)
ax0.set_ylabel("$E_{R}$",fontsize=fs,rotation=0, labelpad=10)
v0 = np.linspace(0.1,20,80.)
c = ['blue', 'red', 'green']
fc = ['lightblue','pink','limegreen']
for i in range(3):
Eb = bands1dvec(v0,i)
ax.fill_between( v0 , Eb[:,0], Eb[:,1], \
color=c[i], facecolor=c[i], alpha=0.7, zorder=4)
ax.plot( v0, Eb[:,1], color=c[i])
ax.plot( v0, Eb[:,0], color=c[i],\
label="Band = %d"%i)
if i==2: labelstr = 'h.osc.'
else: labelstr = None
ax.plot( v0, hosc(i,v0), color='black', label=labelstr, zorder=1 )
ax0.legend( bbox_to_anchor=(1.03,1.00), \
loc='upper left', numpoints=1, \
prop={'size':10}, handlelength=1.1, handletextpad=0.5 )
#gs.tight_layout(fig, rect=[0.,0.,0.8,1.0])
fig.savefig('BandStructure_figures/bands1d_V0.png',dpi=240)
fig.savefig('BandStructure_figures/bands1d_V0.pdf')
we then make plots to compare the band structure in the first quantized version and second quantized version of the Hamiltonian
In [36]:
from matplotlib import rc
rc('font', **{'family':'serif'})
rc('text', usetex=True)
fig = plt.figure(figsize=(3.24,2.7))
gs = matplotlib.gridspec.GridSpec( 1,1, left=0.16, right=0.95, bottom=0.17, top=0.94)
ax0 = fig.add_subplot( gs[0,0])
for ax in [ax0]:
ax.set_xlabel("$V_{0} (E_{R})$",fontsize=12)
ax.grid()
ax.set_ylim(0.,10.)
ax.set_xlim(0.,10.)
ax0.set_ylabel("$E_{R}$",fontsize=12,rotation=0)
v0 = np.linspace(0.1,10,80.)
c = ['blue', 'red', 'green']
fc = ['lightblue','pink','limegreen']
for i in range(2):
Eb = bands1dvec(v0,i)
ax.fill_between( v0 , Eb[:,0], Eb[:,1], \
color=c[i], facecolor=c[i], alpha=0.5, zorder=4)
ax.plot( v0, Eb[:,1], color=c[i], lw=2.)
ax.plot( v0, Eb[:,0], color=c[i], lw=2.,\
label="Band = %d"%i)
if i==1: labelstr = 'h.osc.'
else: labelstr = None
ax.plot( v0, hosc(i,v0), color='black', label=labelstr, zorder=1 )
ax0.legend( bbox_to_anchor=(0.,1.00), \
loc='upper left', numpoints=1, \
prop={'size':8}, handlelength=1.1, handletextpad=0.5 )
fig.savefig('BandStructure_figures/bands1d_V0_firstquant.png',dpi=300)
fig.savefig('BandStructure_figures/bands1d_V0_firstquant.pdf')
In [37]:
fig = plt.figure(figsize=(3.24,2.7))
gs = matplotlib.gridspec.GridSpec( 1,1, left=0.16, right=0.95, bottom=0.17, top=0.94)
ax0 = fig.add_subplot( gs[0,0])
for ax in [ax0]:
ax.set_xlabel("$V_{0} (E_{R})$",fontsize=12)
ax.grid()
ax.set_ylim(-1.,6.)
ax.set_xlim(0.,10.)
ax0.set_ylabel("$E_{R}$",fontsize=12,rotation=0)
v0 = np.linspace(0.1,10,80.)
c = ['blue', 'red', 'green']
fc = ['lightblue','pink','limegreen']
Eb = bands1dvec(v0,0)
E0 = (Eb[:,0] + Eb[:,1])/2
for i in range(2):
Eb = bands1dvec(v0,i)
ax.fill_between( v0 , Eb[:,0]-E0, Eb[:,1]-E0, \
color=c[i], facecolor=c[i], alpha=0.5, zorder=4)
ax.plot( v0, Eb[:,1]-E0, color=c[i], lw=2.)
ax.plot( v0, Eb[:,0]-E0, color=c[i], lw=2.,\
label="Band = %d"%i)
if i==1: labelstr = 'h.osc.'
else: labelstr = None
ax.plot( v0, hosc(i,v0)-E0, color='black', label=labelstr, zorder=1 )
#ax0.legend( bbox_to_anchor=(1.03,1.00), \
# loc='upper left', numpoints=1, \
# prop={'size':12}, handlelength=1.1, handletextpad=0.5 )
fig.savefig('BandStructure_figures/bands1d_V0_secondquant.png',dpi=240)
fig.savefig('BandStructure_figures/bands1d_V0_secondquant.pdf')
In [5]:
def eigenfunsH0( q, v0, x, NPoints=5 ):
r"""Calculate eigenfunctions in 1D lattice
Given the quasimomentum and the lattice depth, this function
will calculate the set of energy eigenfunctions for a number of
bands determined by NPoints. The calculation is for a 1D
lattice. The eigenfunctions will be returned as an array,
evaluated at points given in the parameter x.
Parameters
----------
q : float
The value of quasimomentum. Units are 2\pi, so q can take
values in (-0.5,0.5]
v0 : float
The value of the lattice depth in recoils
x : array-like of floats
The values of x at which the eigenfunctions will be evaluated
NPoints : float, optional
Plane waves e^{imGx} will be considered in the expansion of
the Bloch states such that their momentum satisfies
|m|<NPoints. Here G is the unit reciprocal lattice vector
of the 1D lattice: G = 2*pi / a
A total of 2*NPoints + 1 plane waves is used.
Returns
-------
psi_bands : list
A list with the eigenfuctions sorted by band index. Each
element of the list is an ndarray with length equal to len(x)
len(psi_bands) will be equal to 2*NPoints + 1
Notes
-----
The eigenfunctions are multiplied by an overall phase factor such
that they have matching phases at x0, where x0 should be the position
of a lattice site. This function has a value of x0=0.5 hard-coded
"""
q_Sym, v0_Sym = sym.symbols("q V_{0}")
h0 = np.array( np.array( H0_sin(NPoints,q_Sym,v0_Sym).subs( [(q_Sym,q), (v0_Sym,v0)] )), np.float )
w,v = numpy.linalg.eigh(h0)
a = 1.0 # We use x in units of the lattice spacing
k = np.pi/a
# m are the indices for the reciprocal lattice vectors
m = np.linspace(-NPoints,NPoints,2*NPoints+1)
# The sum of plane waves is carried out for each of the bands
psi_bands = []
for i in np.argsort(w):
coefs = v[:,i][:,None]
#if i==0:
# print coefs
planewaves = coefs * np.exp( 1j*(2*k)*np.outer(m,x) ) * np.exp(1j*q*2*k*x)
psi = planewaves.sum(axis=0)
# Adjust the wavefunctions so they have matching phases at x0
# where x0 is the position of a lattice site.
#
# This will turn out to be important when defining the Wannier
# wavefunctions.
#
x0 = 0.5
phase = -1*np.exp( 1j* np.angle( np.sum(coefs * np.exp( 1j*(2*k)*m*x0 ) * np.exp(1j*q*k*x0))))
psi = psi / phase
psi_bands.append((psi,w[i]))
return psi_bands
We can use the numeric version of the eigenfunctions to plot them for various lattice dephts
In [107]:
# Start matplotlib
from matplotlib import rc
rc('text', usetex=True)
plt.rcParams['text.latex.preamble'] = [
r'\usepackage{bm}', # for bold math
r'\usepackage{siunitx}', # i need upright \micro symbols, but you need...
r'\sisetup{detect-all}', # ...this to force siunitx to actually use your fonts
#r'\usepackage[textlf]{MyriadPro}',
#r'\usepackage{sansmath}', # load up the sansmath so that math -> helvet
#r'\sansmath', # <- tricky! -- gotta actually tell tex to use!
]
plt.rcParams['axes.linewidth'] = 0.7
plt.rcParams['patch.linewidth'] = 0.4
depths = [4., 8., 12., 16.]
c = ['blue', 'green','red','black','purple']
fs=12
fig = plt.figure(figsize=(7.0,3))
gs = matplotlib.gridspec.GridSpec( 1,2,\
left=0.10, right=0.86, bottom=0.14, top=0.93,
hspace=0., wspace=0.40)
ax0 = fig.add_subplot( gs[0,0])
ax1 = fig.add_subplot( gs[0,1])
for ax in [ax0,ax1]:
ax.set_xlabel("$x/a$",fontsize=fs)
ax.set_ylabel("$|\psi_{q}^{n}(x)|^{2}$",fontsize=fs, rotation=90, labelpad=8)
ax.set_xticks([-2, -1, 0, 1, 2])
ax.grid(alpha=0.5)
ax0.set_title('$q = 0 $',fontsize=fs)
ax1.set_title('$q = \pi/a$',fontsize=fs)
qbot = 0.
qtop = 0.5
x = np.linspace(-2,2,200)
for i,d in enumerate(depths):
psi0 = eigenfunsH0( qbot, d, x )[0]
psi1 = eigenfunsH0( qtop, d, x )[0]
ax0.plot( x, np.abs(psi0[0])**2, color = c[i] )
ax1.plot( x, np.abs(psi1[0])**2, color = c[i], label='$%.1d E_{R}$'%d)
ax0.fill_between( x, 0., np.sin(np.pi*x)**2, facecolor='gray', alpha=0.3)
ax1.fill_between( x, 0., np.sin(np.pi*x)**2, facecolor='gray', alpha=0.3)
ax1.legend( bbox_to_anchor=(1.05,1.00), \
loc='upper left', numpoints=1, \
prop={'size':10}, handlelength=1.1, handletextpad=0.5 )
#gs.tight_layout(fig, rect=[0.,0.,0.9,1.0])
fig.savefig('BandStructure_figures/eigenfuns1d.png',dpi=240)
fig.savefig('BandStructure_figures/eigenfuns1d.pdf')
We will consider the quasimomentun in the first Brilloun zone as $q \in [0,2\pi)\ $. Numerically, we will consider a lattice with $L\ $ sites and quasimomentum values than run from $0\ $ to $(L-1)/L\ $. In other words, we define the quasimomentum q without the factor of $2\pi\ $, so we need to be careful to write the factor down when it is necessary.
The expansion of the eigenstates of the problem in terms of plane waves is given by
$$ \psi_{q}^{n}(x) = \sum_{m \in \mathbb{Z} } c_{qm}^{n} e^{im2\pi x} $$and the Wannier functions are defined in terms of the eigenstates as
$$ w_{j}^{n} = \frac{1}{L} \sum_{q} e^{-iqx_{j}} \psi_{q}^{n}(x) $$Putting these two definitions together we arrive at
(see for instance http://alps.comp-phys.org/mediawiki/index.php/Documentation:Bosons_in_an_optical_lattice)
$$ w^{n} = \frac{1}{L} \sum_{m \in \mathbb{Z}} \sum_{q} c_{qm}^{n} e^{ix(mG+q)} $$Here $G=2\pi\ $ , and since the units of $q$ are also $2\pi\ $, as was stated above, we can write this as
$$ w^{n} = \frac{1}{L} \sum_{m \in \mathbb{Z}} \sum_{q} c_{qm}^{n} e^{ix2\pi(m+q)} $$The symmetry of the problem can be used to realize that
$$ c_{qm}^{n} = c_{q'm'}^{n} $$if
$$q+m = -(q'+m')\ \ \ $$which enables us to write the sum as $$ w^{n} = \frac{1}{L} \left( c_{00}^{n} + \sum_{m>0} \sum_{q>0} c_{qm}^{n} 2\cos(2\pi(m+q)x \right) $$
where we have to be careful to only sum over |(m+q)| once. This is the expression that we use to numerically determine the Wannier functions
In [6]:
def wannier( x, v0, band=0 ):
L = 200. # Lattice size
qset = np.arange( 0, L)/L
qsetSize = len(qset)
num_m = 20 # Number of plane waves
mset = np.arange(-num_m, num_m+1)
msetSize = len(mset)
# Diagonalize for all values of q
cCoefs = np.empty( (qsetSize, msetSize, msetSize) )
for qi,q in enumerate(qset):
h0 = np.array( np.array(H0_sin(num_m, q, v0)))
w,v = numpy.linalg.eigh(h0)
for j in range(msetSize):
if v[:,j][-1] > 0 :
cCoefs[qi][j] = v[:,j]
else:
cCoefs[qi][j] = -1.*v[:,j]
#Calculate Wannier
wan = np.zeros_like( x, dtype=np.float)
for qi,q in enumerate(qset):
for mj,m in enumerate(mset):
cqm = cCoefs[qi][band][mj]
if q + m == 0:
wan += cqm
if q + m > 0:
#print "added ",
if band%2 == 0:
wan += cqm * 2* np.cos( 2*np.pi* (q+m) * x )
else:
wan += cqm * 2* np.sin( 2*np.pi* (q+m) * x )
wan = wan / L
return wan
One can check that the Wannier functions come out (almost) properly normalized:
In [7]:
from scipy import integrate
x = linspace(-5, 5, 1000)
v0 = 5.
norm = integrate.simps(np.abs(wannier(x, v0))**2,x)
print norm
We then plot a few Wannier states for some lattice depths to show how they become increasingly localized
In [106]:
# Start matplotlib
from matplotlib import rc
rc('text', usetex=True)
plt.rcParams['text.latex.preamble'] = [
r'\usepackage{bm}', # for bold math
r'\usepackage{siunitx}', # i need upright \micro symbols, but you need...
r'\sisetup{detect-all}', # ...this to force siunitx to actually use your fonts
#r'\usepackage[textlf]{MyriadPro}',
#r'\usepackage{sansmath}', # load up the sansmath so that math -> helvet
#r'\sansmath', # <- tricky! -- gotta actually tell tex to use!
]
plt.rcParams['axes.linewidth'] = 0.7
plt.rcParams['patch.linewidth'] = 0.4
fs=12
fig = plt.figure(figsize=(4.2,2.68))
gs = matplotlib.gridspec.GridSpec( 1,1,
left=0.16, right=0.78, bottom=0.16, top=0.96)
ax0 = fig.add_subplot( gs[0,0])
for ax in [ax0]:
ax.set_xlabel("$x/a$",fontsize=fs)
ax.set_xticks([-3,-2, -1, 0, 1, 2, 3])
ax.grid(alpha=0.5)
ax0.set_ylabel("$w^{n}(x-x_{j})$",fontsize=fs,rotation=90, labelpad=4)
#ax0.set_title('Real part',fontsize=12)
#ax1.set_title('Imaginary part',fontsize=12)
x = np.linspace(-3,3,1001)
depths = [1., 4., 7., 14. ]
for i,d in enumerate(depths):
w = wannier(x, d)
labelstr = '$%.0fE_{R}$'%d
ax0.plot(x,np.real(w), label=labelstr,lw=1.8)
ax0.fill_between( x, 0., np.sin(np.pi*x)**2, facecolor='gray', alpha=0.3)
ax0.legend( bbox_to_anchor=(1.03,1.00), \
loc='upper left', numpoints=1, \
prop={'size':10}, handlelength=1.1, handletextpad=0.5 )
#gs.tight_layout(fig, rect=[0.,0.,0.85,1.0])
fig.savefig('BandStructure_figures/wannier1d_V0.png',dpi=300)
fig.savefig('BandStructure_figures/wannier1d_V0.pdf')
We can also plot the wannier states for different bands
In [105]:
# Start matplotlib
from matplotlib import rc
rc('text', usetex=True)
plt.rcParams['text.latex.preamble'] = [
r'\usepackage{bm}', # for bold math
r'\usepackage{siunitx}', # i need upright \micro symbols, but you need...
r'\sisetup{detect-all}', # ...this to force siunitx to actually use your fonts
#r'\usepackage[textlf]{MyriadPro}',
#r'\usepackage{sansmath}', # load up the sansmath so that math -> helvet
#r'\sansmath', # <- tricky! -- gotta actually tell tex to use!
]
plt.rcParams['axes.linewidth'] = 0.7
plt.rcParams['patch.linewidth'] = 0.4
fs=12
fig = plt.figure(figsize=(4.2,2.88))
gs = matplotlib.gridspec.GridSpec( 1,1,
left=0.16, right=0.78, bottom=0.16, top=0.91)
ax0 = fig.add_subplot( gs[0,0])
for ax in [ax0]:
ax.set_xlabel("$x/a$",fontsize=fs)
ax.set_xticks([-3,-2, -1, 0, 1, 2, 3])
ax.grid()
ax0.set_ylabel("$w^{n}(x-x_{j})$",fontsize=fs,rotation=90, labelpad=4)
ax0.set_title('$V_{0}=4 E_{R}$',fontsize=fs)
#ax1.set_title('Imaginary part',fontsize=12)
x = np.linspace(-3,3,1001)
v0 = 4.
bands = [0, 1, 2]
for i,b in enumerate(bands):
w = wannier(x, v0, band=b)
labelstr = 'Band %d'%b
ax0.plot(x,np.real(w), label=labelstr, alpha=(1-0.25*b),lw=1.8, zorder=5-b)
ax0.fill_between( x, 0., np.sin(np.pi*x)**2, facecolor='gray', alpha=0.3)
ax0.legend( bbox_to_anchor=(1.01,1.00), \
loc='upper left', numpoints=1, \
prop={'size':10}, handlelength=1.1, handletextpad=0.5 )
#gs.tight_layout(fig, rect=[0.,0.,0.82,1.0])
fig.savefig('BandStructure_figures/wannier1d_bands.png',dpi=300)
fig.savefig('BandStructure_figures/wannier1d_bands.pdf')
In [8]:
def bands3d( v0, NBand=0 ):
r"""Returns an array of shape (2,) which has the bottom and top of the band energies
Given the lattice depth it will calculate the bottom and the
top of the band with band index equal to NBand.
Parameters
----------
v0 : array-like
len(v0) must be equal to 3. Each entry corresponds to the
lattice depth in recoils along x, y, z respectively
NBand :
band index (default = 0 )
Returns
-------
array-like
array of shape (2,) which has the bottom of the band as
first entry and the top of the band as second entry.
"""
assert len(v0)==3
order = np.argsort( v0 )
l=0; m=0; n=0
while NBand > 0:
l+=1; NBand-=1
if NBand > 0:
m+=1; NBand-=1
if NBand > 0:
n+=1; NBand-=1
band1dindex = np.array([ l, m, n] )
qbot = 0.
qtop = 0.5
bandbot = 0.
bandtop = 0.
for i,o in enumerate(order):
in1d = band1dindex[i]
if in1d%2 ==0:
bandbot += eigenvalsH0(qbot, v0[o])[in1d]
bandtop += eigenvalsH0(qtop, v0[o])[in1d]
else:
bandbot += eigenvalsH0(qtop, v0[o])[in1d]
bandtop += eigenvalsH0(qbot, v0[o])[in1d]
return np.array([bandbot, bandtop])
def bands3dvec( v0, NBand=0):
r""" Returns an array of shape (len(v0,2) which has the bottom and
top of the band energies
Parameters
----------
v0 : array-like
The shape of v0 is (len(v0,3). Each row of v0 has 3
values. Each value corresponds to the lattice depth
in recoils along x, y, z respectively.
NBand : optional
Band index (default = 0 )
"""
band = numpy.empty( (len(v0), 2) )
for i,v in enumerate(v0):
band[i] = bands3d(v, NBand)
return band
def hosc3d( NBand, v0 ):
r"""Energy of the NBand^th excited state of the 3D harmonic oscillator"""
return 2* np.sqrt( v0 ) * (NBand + 1.5)
We make a plot with the band structure of the 3D lattice
In [104]:
# Start matplotlib
from matplotlib import rc
rc('text', usetex=True)
plt.rcParams['text.latex.preamble'] = [
r'\usepackage{bm}', # for bold math
r'\usepackage{siunitx}', # i need upright \micro symbols, but you need...
r'\sisetup{detect-all}', # ...this to force siunitx to actually use your fonts
#r'\usepackage[textlf]{MyriadPro}',
#r'\usepackage{sansmath}', # load up the sansmath so that math -> helvet
#r'\sansmath', # <- tricky! -- gotta actually tell tex to use!
]
plt.rcParams['axes.linewidth'] = 0.6
plt.rcParams['patch.linewidth'] = 0.4
fs = 12
fig = plt.figure(figsize=(4.4,2.6))
gs = matplotlib.gridspec.GridSpec( 1,1, left=0.11, right=0.74, bottom=0.16, top=0.96)
ax0 = fig.add_subplot( gs[0,0])
for ax in [ax0]:
ax.set_xlabel("$V_{0} (E_{R})$",fontsize=fs)
ax.grid()
ax.set_ylim(0.,15.)
ax.set_xlim(0.,20.)
ax0.set_ylabel("$E_{R}$",fontsize=fs,rotation=0, labelpad=10)
v0 = np.linspace(0.1,20,80.)
v03d = np.outer( v0, np.ones(3))
c = ['blue', 'red', 'green']
fc = ['lightblue','pink','limegreen']
for i in range(3):
Eb = bands3dvec(v03d,i)
ax.fill_between( v0 , Eb[:,0], Eb[:,1], \
color=c[i], facecolor=c[i], alpha=0.7, zorder=4)
ax.plot( v0, Eb[:,1], color=c[i])
ax.plot( v0, Eb[:,0], color=c[i],\
label="Band = %d"%i)
if i==2: labelstr = 'h.osc.'
else: labelstr = None
ax.plot( v0, hosc3d(i,v0), color='black', label=labelstr, zorder=1 )
ax0.legend( bbox_to_anchor=(1.03,1.00), \
loc='upper left', numpoints=1, \
prop={'size':10}, handlelength=1.1, handletextpad=0.5 )
#gs.tight_layout(fig, rect=[0.,0.,0.8,1.0])
fig.savefig('BandStructure_figures/bands3d_V0.png',dpi=240)
fig.savefig('BandStructure_figures/bands3d_V0.pdf')
The interpolation data is saved to file for easy retrieval without having to evaluate the entire thing
In [7]:
v0 = np.linspace(0.,120,320.)
NPoints = 5
evals0 = np.empty( (len(v0), 2*NPoints+1) )
evals1 = np.empty( (len(v0), 2*NPoints+1) )
qbot = 0.
qtop = 0.5
for i,v in enumerate(v0):
evals0[i] = eigenvalsH0(qbot, v, NPoints=NPoints)
evals1[i] = eigenvalsH0(qtop, v, NPoints=NPoints)
np.savetxt('banddat/interpdat_B1D_v0.dat', v0)
for n in range( 2*NPoints+1 ):
np.savetxt('banddat/interpdat_B1D_0_%d.dat'%n, evals0[:,n] )
np.savetxt('banddat/interpdat_B1D_1_%d.dat'%n, evals1[:,n] )
The interpolation data is retrieved from disk and we use it to define functions that calculate the band structure in the 3D lattice. The code that defines the interpolation functions is imported below
In [9]:
from Band3DInterp import *
In [123]:
%timeit eigenvalsH0(0., 10.)
%timeit interp0[0](10.)
In [124]:
v0 = np.linspace(0.1,20,80.)
v03d = np.outer( v0, np.ones(3))
band = 0
%timeit bands3dvec(v03d, NBand=band)
%timeit bands3dinterpvec(v03d, NBand=band)
In [125]:
print "Performance is increased by a factor of %d" % (4230 / 16.7)
We can use the interpolation version of the 3D bands calculation to plot the band structure
In [126]:
fig = plt.figure(figsize=(6.5,3.8))
gs = matplotlib.gridspec.GridSpec( 1,1)
ax0 = fig.add_subplot( gs[0,0])
for ax in [ax0]:
ax.set_xlabel("$V_{0} (E_{R})$",fontsize=16)
ax.grid()
ax.set_ylim(0.,15.)
ax.set_xlim(0.,20.)
ax0.set_ylabel("$E_{R}$",fontsize=16,rotation=0)
v0 = np.linspace(0.1,20,80.)
v03d = np.outer( v0, np.ones(3))
c = ['blue', 'red', 'green']
fc = ['lightblue','pink','limegreen']
for i in range(3):
Eb = bands3dinterpvec(v03d,i)
ax.fill_between( v0 , Eb[:,0], Eb[:,1], \
color=c[i], facecolor=c[i], alpha=0.7, zorder=4)
ax.plot( v0, Eb[:,1], color=c[i])
ax.plot( v0, Eb[:,0], color=c[i],\
label="Band = %d"%i)
if i==2: labelstr = 'h.osc.'
else: labelstr = None
ax.plot( v0, hosc3d(i,v0), color='black', label=labelstr, zorder=1 )
ax0.legend( bbox_to_anchor=(1.03,1.00), \
loc='upper left', numpoints=1, \
prop={'size':12}, handlelength=1.1, handletextpad=0.5 )
gs.tight_layout(fig, rect=[0.,0.,0.8,1.0])
fig.savefig('BandStructure_figures/bands3d_V0_interp.png',dpi=240)
fig.savefig('BandStructure_figures/bands3d_V0_interp.pdf')
In [4]:
def tCalc(v0, Xij=1., Nqs=40):
r"""Calculates the tunneling matrix element
Parameters
----------
v0 : float
Lattice depth in recoils
Xij : float, optional,
This is the separation between wannier states. For
nearest neighbors Xij = 1. The type is a float, but
really only integers should be fed into it.
default value = 1.0 (nearest neighbors)
Nqs : int, optional
Then number of sites on the lattice (along one of the
dimensions). This number defines the number of quasi-momentum
points that will be considered in the sum that calculates
the tunnling. The result should be independent of Nqs
default value = 40
Returns
-------
t: float
tunneling matrix element in recoils
"""
qset = np.arange( 0, Nqs)/Nqs
t = 0.
for q in qset:
t += eigenvalsH0(q, v0)[0] * np.exp( 1j*2*np.pi*q*Xij) / len(qset)
return -1*np.real(t)
Calculate the value of the tunneling in the limit of zero lattice depth:
In [8]:
v0 = 0.1
print tCalc(v0 )
v0 = 0.01
print tCalc(v0 )
v0 = 0.001
print tCalc(v0 )
v0 = 0.00001
print tCalc(v0 )
We would like to have an interpolation function for the nearest neighbor tunneling matrix element. Here we calculate the samples necessary to define the interpolation function.
In [60]:
v0 = np.linspace( 0.1, 80., 120.)
tv0 = np.array([tCalc(d) for d in v0])
np.savetxt('banddat/interpdat_t_v0.dat', v0)
np.savetxt('banddat/interpdat_t_tCalc.dat', tv0)
And below we load up the interpolation samples from disk and define the intepolartion function.
In [11]:
#Here the interpolation data is loaded from disk
from scipy.interpolate import interp1d
tInterp = interp1d( np.loadtxt('banddat/interpdat_t_v0.dat'), np.loadtxt('banddat/interpdat_t_tCalc.dat'))
We want to make a plot that illustrates the difference in tunneling matrix element between nearest-neighbor tunneling and beyond nearest-neighbor tunneling. We calculate the samples first:
In [13]:
v0 = np.linspace(0.1,8,18)
tCalc_2sites = np.array([ np.abs(tCalc(v0val, Xij=2.)) for v0val in v0 ])
tCalc_3sites = np.array([ np.abs(tCalc(v0val, Xij=3.)) for v0val in v0 ])
And then we make the plot
In [15]:
# Start matplotlib
from matplotlib import rc
rc('text', usetex=True)
plt.rcParams['text.latex.preamble'] = [
r'\usepackage{bm}', # for bold math
r'\usepackage{siunitx}', # i need upright \micro symbols, but you need...
r'\sisetup{detect-all}', # ...this to force siunitx to actually use your fonts
#r'\usepackage[textlf]{MyriadPro}',
#r'\usepackage{sansmath}', # load up the sansmath so that math -> helvet
#r'\sansmath', # <- tricky! -- gotta actually tell tex to use!
]
plt.rcParams['axes.linewidth'] = 0.6
plt.rcParams['patch.linewidth'] = 0.4
fs=12.
v0 = np.linspace(0.1,8,18)
fig = plt.figure(figsize=(4.4,2.66))
gs = matplotlib.gridspec.GridSpec( 1,1,
left=0.15, right=0.72, bottom=0.16, top=0.96)
ax0 = fig.add_subplot( gs[0,0])
for ax in [ax0]:
ax.set_xlabel("$V_{0} (E_{R})$",fontsize=fs)
ax.grid(alpha=0.5)
ax.set_ylim(0.,0.2)
#ax.set_xlim(0.,20.)
ax0.set_ylabel("$|t_{ij}|$ $(E_{R})$",fontsize=fs,rotation=90, labelpad=7)
c = ['blue', 'red', 'green']
fc = ['lightblue','pink','limegreen']
ax.plot( v0, tInterp(v0), label="$|\Delta_{ij}|/a=1$")
ax.plot( v0, tCalc_2sites, label="$|\Delta_{ij}|/a=2$" )
ax.plot( v0, tCalc_3sites, label="$|\Delta_{ij}|/a=3$ ")
ax0.legend( bbox_to_anchor=(1.02,1.00), \
loc='upper left', numpoints=1, \
prop={'size':10}, handlelength=1.1, handletextpad=0.5 )
#gs.tight_layout(fig, rect=[0.,0.,0.8,1.0])
fig.savefig('BandStructure_figures/tightbinding_V0_interp.png',dpi=300)
fig.savefig('BandStructure_figures/tightbinding_V0_interp.pdf')
We show the tunneling (nearest-neighbor) in a logarithmic scale alongside the Mathieu approximation to the tunneling. We use recoils and kHz
In [73]:
def tMathieu( v0 ):
return 4./np.sqrt(np.pi) * v0**(3./4.) * np.exp( -2*np.sqrt(v0))
In [102]:
# Start matplotlib
from matplotlib import rc
rc('text', usetex=True)
plt.rcParams['text.latex.preamble'] = [
r'\usepackage{bm}', # for bold math
r'\usepackage{siunitx}', # i need upright \micro symbols, but you need...
r'\sisetup{detect-all}', # ...this to force siunitx to actually use your fonts
#r'\usepackage[textlf]{MyriadPro}',
#r'\usepackage{sansmath}', # load up the sansmath so that math -> helvet
#r'\sansmath', # <- tricky! -- gotta actually tell tex to use!
]
plt.rcParams['axes.linewidth'] = 0.6
plt.rcParams['patch.linewidth'] = 0.4
fs=12.
fig = plt.figure(figsize=(7.7,2.94))
gs = matplotlib.gridspec.GridSpec( 1,2,\
left=0.09, right=0.84, bottom=0.15, top=0.96,
wspace=0.35, hspace=0.)
ax0 = fig.add_subplot( gs[0,0])
ax1 = fig.add_subplot( gs[0,1])
for ax in [ax0,ax1]:
ax.set_xlabel("$V_{0} (E_{R})$",fontsize=fs)
ax.grid(alpha=0.5)
#ax.set_xlim(0.,20.)
ax0.set_yscale("log")
ax0.set_ylim(0.002,0.6)
ax0.set_ylabel(r"$t/E_{r}$",fontsize=fs,rotation=90)
ax1.set_ylabel(r"$t/h\ \mathrm{(kHz)}$",fontsize=fs,rotation=90)
ax1.set_xlim(4., 10)
ax1.set_ylim(0.4, 2.6)
v0 = np.linspace(0.1,25,100.)
c = ['blue', 'red', 'green']
fc = ['lightblue','pink','limegreen']
ax0.plot( v0, tInterp(v0), label="Diag.")
ax0.plot( v0, tMathieu(v0), label="Mathieu")
ax1.plot( v0, tInterp(v0) * 29.0, label="Diag.") #29kHz = 1Er
ax1.plot( v0, tMathieu(v0) * 29.0, label="Mathieu")
ax1.legend( bbox_to_anchor=(1.03,1.00), \
loc='upper left', numpoints=1, \
prop={'size':10}, handlelength=1.1, handletextpad=0.5 )
#gs.tight_layout(fig, rect=[0.,0.,0.88,1.0])
fig.savefig('BandStructure_figures/tunneling_V0_Mathieu.png',dpi=300)
fig.savefig('BandStructure_figures/tunneling_V0_Mathieu.pdf')
We define the second derivative below
In [146]:
def wannierD2( x, v0, band=0 ):
L = 200. # Lattice size
qset = np.arange( 0, L)/L
qsetSize = len(qset)
num_m = 20 # Number of plane waves
mset = np.arange(-num_m, num_m+1)
msetSize = len(mset)
# Diagonalize for all values of q
cCoefs = np.empty( (qsetSize, msetSize, msetSize) )
for qi,q in enumerate(qset):
h0 = np.array( np.array(H0_sin(num_m, q, v0)))
w,v = numpy.linalg.eigh(h0)
for j in range(msetSize):
if v[:,j][-1] > 0 :
cCoefs[qi][j] = v[:,j]
else:
cCoefs[qi][j] = -1.*v[:,j]
#Calculate Wannier
wan = np.zeros_like( x, dtype=np.float)
for qi,q in enumerate(qset):
for mj,m in enumerate(mset):
cqm = cCoefs[qi][band][mj]
if q + m > 0:
wan += -1* cqm * 8*(( np.pi*(q+m) )**2)* np.cos( 2*np.pi* (q+m) * x )
wan = wan / L
return wan
The second derivative is useful in calculating the tunneling matrix element directely from the Wannier function overlap
$$ t = - \int w_{j}^{*}(x) \left( -\frac{\hbar^{2}}{2m} \frac{\partial^{2}}{\partial x^{2}} + V_{0}\sin^{2}(kx) \right) w_{j+1}(x) $$below we define a function that thas such a calculation and compare its results with the sum over energy eigenvalues method used above
In [147]:
from scipy import integrate
def tWannier( v0):
x0 = np.linspace(-5, 5, 1000)
x1 = np.linspace(-6, 4, 1000)
w0 = wannier( x0, v0 )
w1 = wannier( x1, v0 )
w1D2 = wannierD2( x1, v0 )
integrand = -1*( -1*w0*w1D2/(np.pi**2) + w0*v0*(np.sin(np.pi*x0)**2)*w1 )
tInt = integrate.simps( integrand, x0)
return tInt
In [150]:
v0list= [1.0, 3.0, 7.0, 10.]
print "v0\t"
print "----"
for v0 in v0list:
print "%4.2f\t%4.5f\t%4.5f\t%4.5f\t%4.5f" % (v0, tCalc(v0), tInterp(v0), tMathieu(v0), tWannier(v0) )
We see that for lower lattice depths the two methods differ. We look at the performance of the two methods and find that the one with the energy eigenvalues beats the other one by more than a factor of 10.
In [45]:
%timeit tWannier(7.0)
%timeit tCalc(7.0)
In [86]:
from scipy import integrate
def wF( v0 ):
x0 = np.linspace(-5, 5, 1000)
w0 = wannier( x0, v0 )
return 8/np.pi * (integrate.simps( w0**4, x0 ))**3
def wFho( v0 ):
return np.sqrt( 8*np.pi) * v0**(3./4.)
We calculate a table of values to be used in interpolation
In [87]:
v0 = np.hstack(( np.linspace(1e-3,0.1,4),np.linspace( 0.1, 80., 120.) ))
wFv0 = np.array([wF(d) for d in v0])
np.savetxt('banddat/interpdat_wF_v0.dat', v0)
np.savetxt('banddat/interpdat_wF_wF.dat', wFv0)
The interpolation data is loaded from disk and the function is defined
In [18]:
#Here the interpolation data is loaded from disk
from scipy.interpolate import interp1d
wFInterp = interp1d( np.loadtxt('banddat/interpdat_wF_v0.dat'), np.loadtxt('banddat/interpdat_wF_wF.dat'))
And we make a plot comparing the numerical result using Wannier functions and the harmonic approximation
In [101]:
# Start matplotlib
from matplotlib import rc
rc('text', usetex=True)
plt.rcParams['text.latex.preamble'] = [
r'\usepackage{bm}', # for bold math
r'\usepackage{siunitx}', # i need upright \micro symbols, but you need...
r'\sisetup{detect-all}', # ...this to force siunitx to actually use your fonts
#r'\usepackage[textlf]{MyriadPro}',
#r'\usepackage{sansmath}', # load up the sansmath so that math -> helvet
#r'\sansmath', # <- tricky! -- gotta actually tell tex to use!
]
plt.rcParams['axes.linewidth'] = 0.6
plt.rcParams['patch.linewidth'] = 0.4
fs=12.
fig = plt.figure(figsize=(7.7,2.94))
gs = matplotlib.gridspec.GridSpec( 1,2,\
left=0.09, right=0.84, bottom=0.15, top=0.96,
wspace=0.35, hspace=0.)
ax0 = fig.add_subplot( gs[0,0])
ax1 = fig.add_subplot( gs[0,1])
for ax in [ax0,ax1]:
ax.set_xlabel("$V_{0} (E_{R})$",fontsize=fs)
ax.grid()
#ax.set_xlim(0.,20.)
#ax.set_yscale("log")
#ax0.set_ylim(0.002,0.6)
ax0.set_ylabel("$U/a_{s}$ $(E_{R}/a)$",fontsize=fs,rotation=90)
ax1.set_ylabel("$U/a_{s}$ ($\mathrm{kHz}/a$)",fontsize=fs,rotation=90)
#ax1.set_xlim(5., 25)
#ax1.set_ylim(0.01, 4.)
v0 = np.linspace(0.1,25,100.)
c = ['blue', 'red', 'green']
fc = ['lightblue','pink','limegreen']
ax0.plot( v0, wFInterp(v0), label="Wannier")
ax0.plot( v0, wFho(v0), label="H.osc.")
ax1.plot( v0, wFInterp(v0) * 29.0, label="Wannier") #29kHz = 1Er
ax1.plot( v0, wFho(v0) * 29.0, label="H.osc")
ax1.legend( bbox_to_anchor=(1.03,1.00), \
loc='upper left', numpoints=1, \
prop={'size':10}, handlelength=1.1, handletextpad=0.5 )
#gs.tight_layout(fig, rect=[0.,0.,0.88,1.0])
fig.savefig('BandStructure_figures/wFactor_V0.png',dpi=240)
fig.savefig('BandStructure_figures/wFactor_V0.pdf')
In [100]:
# Start matplotlib
from matplotlib import rc
rc('text', usetex=True)
plt.rcParams['text.latex.preamble'] = [
r'\usepackage{bm}', # for bold math
r'\usepackage{siunitx}', # i need upright \micro symbols, but you need...
r'\sisetup{detect-all}', # ...this to force siunitx to actually use your fonts
#r'\usepackage[textlf]{MyriadPro}',
#r'\usepackage{sansmath}', # load up the sansmath so that math -> helvet
#r'\sansmath', # <- tricky! -- gotta actually tell tex to use!
]
plt.rcParams['axes.linewidth'] = 0.6
plt.rcParams['patch.linewidth'] = 0.4
fs=12.
fig = plt.figure(figsize=(7.7,2.94))
gs = matplotlib.gridspec.GridSpec( 1,2,\
left=0.09, right=0.84, bottom=0.15, top=0.96,
wspace=0.35, hspace=0.)
ax0 = fig.add_subplot( gs[0,0])
ax1 = fig.add_subplot( gs[0,1])
for ax in [ax0,ax1]:
ax.set_xlabel("$V_{0} (E_{R})$",fontsize=fs)
ax.grid()
#ax.set_xlim(0.,20.)
#ax.set_yscale("log")
#ax0.set_ylim(0.002,0.6)
ax0.set_ylabel(r'$\frac{U}{a_{s}}\ ( \frac{E_{r}}{100a_{0}})$', fontsize=fs, rotation=90 )
ax1.set_ylabel(r'$\frac{U/h}{a_{s}}\ ( \frac{\mathrm{kHz}}{100a_{0}})$', fontsize=fs, rotation=90 )
#ax0.set_ylabel(r'$U/a_{s}\ (E_{r}/100a_{0})$', fontsize=fs, rotation=90 )
#ax1.set_ylabel(r"$\frac{U}{a_{s}}$ $\left( \frac{\mathrm{kHz}}{100 a_{0}} \right)$",fontsize=fs,rotation=0)
#ax1.set_xlim(5., 25)
#ax1.set_ylim(0.01, 4.)
a0a = 100.* 5.29e-11 / (1064e-9/2)
v0 = np.linspace(0.1,12,100.)
c = ['blue', 'red', 'green']
fc = ['lightblue','pink','limegreen']
ax0.plot( v0, a0a*wFInterp(v0), label="Wannier")
ax0.plot( v0, a0a*wFho(v0), label="H.osc.")
ax1.plot( v0, a0a*wFInterp(v0) * 29.0, label="Wannier") #29kHz = 1Er
ax1.plot( v0, a0a*wFho(v0) * 29.0, label="H.osc")
ax1.legend( bbox_to_anchor=(1.03,1.00), \
loc='upper left', numpoints=1, \
prop={'size':10}, handlelength=1.1, handletextpad=0.5 )
#gs.tight_layout(fig, rect=[0.,0.,0.88,1.0])
fig.savefig('BandStructure_figures/wFactor_V0_a0.png',dpi=300)
fig.savefig('BandStructure_figures/wFactor_V0_a0.pdf')
In [162]:
# Calculation of the lattice depth at which next-nearest neighbor tunneling is 10 times slower than nearest neighbor tunneling
v0int = np.linspace( 2., 8., 40)
tun0 = tInterp(v0int)
tun1 = [ np.abs(tCalc(v0val, Xij=2.)) for v0val in v0int]
from scipy.interpolate import interp1d
tun1f = interp1d( v0int, tun1)
def g(x):
return tun1f(x)*20 - tInterp(x)
print brentq( g, 2.1, 6.)
In [19]:
def bandgap( v0 ):
v03d = np.outer( v0, np.ones(3))
Eb0 = bands3dinterpvec(v03d,0 )[:,1]
Eb1 = bands3dinterpvec(v03d,1 )[:,0]
return Eb1-Eb0
def tunnel(v0):
return tInterp(v0)
def Uonsite( v0, aS):
a0a = 5.29e-11 / (1064e-9/2)
return aS*a0a*wFInterp(v0)
# Find lines of constant U/t
def aSconstInt( v0, Ut):
return Ut * tInterp(v0) / Uonsite(v0,1.)
# Find the line where W = multiple*U , where W is the band gap.
def aSlimit( v0, multiple):
aS = bandgap(v0) / (multiple*Uonsite(v0,1.))
return aS
# Find the lattice depth where 12t = W/10
def f(x):
return bandgap(x)[0]/10. - 12.*tunnel(x)
from scipy.optimize import brentq
print brentq( f, 1., 10.)
def g(x):
return 2.*np.sqrt(x) - 12.*4/(np.sqrt(np.pi)) * (x**(3./4.)) * np.exp( -2.*np.sqrt(x))
print brentq( g, 1., 10.)
v0test = 3.
print type(tunnel(v0test))
print "v0 = {:0.2f}".format(v0test)
print " W = {:0.2f}".format(bandgap(v0test)[0])
print " U = {:0.2f}".format(Uonsite(v0test,200.))
print " t = {:0.2f}".format(float(tunnel(v0test)))
print " a_max = {:0.1f}".format( aSlimit(7.,10.)[0] )
# Start matplotlib
from matplotlib import rc
rc('text', usetex=True)
plt.rcParams['text.latex.preamble'] = [
r'\usepackage{bm}', # for bold math
r'\usepackage{siunitx}', # i need upright \micro symbols, but you need...
r'\sisetup{detect-all}', # ...this to force siunitx to actually use your fonts
#r'\usepackage[textlf]{MyriadPro}',
#r'\usepackage{sansmath}', # load up the sansmath so that math -> helvet
#r'\sansmath', # <- tricky! -- gotta actually tell tex to use!
]
plt.rcParams['axes.linewidth'] = 0.6
plt.rcParams['patch.linewidth'] = 0.4
fs=12.
fig = plt.figure(figsize=(5.1,2.94))
gs = matplotlib.gridspec.GridSpec( 1,1,\
left=0.12, right=0.66, bottom=0.15, top=0.96,
wspace=0.35, hspace=0.)
ax = fig.add_subplot( gs[0])
v0 = np.linspace(2.24,20,1000)
dash = [(4.8,2.8), (6.4,2,2,2),(2,0.001),(2,2.10)]
a10 = aSlimit(v0,10.)
a05 = aSlimit(v0,5.)
a03 = aSlimit(v0,3.)
ax.plot( v0, a10, label=r'$\Delta/U=10$', color='k', dashes=dash[0])
ax.plot( v0, a05, label=r'$\Delta/U=5$', color='k', dashes=dash[1])
ax.plot( v0, a03, label=r'$\Delta/U=3$', color='k', dashes=dash[3])
ax.fill_between( v0, np.ones_like(v0), a10, lw=0., color='gray', alpha=0.0)
ax.fill_between( v0, a10, a05, lw=0., color='gray', alpha=0.30)
ax.fill_between( v0, a05, a03, lw=0., color='gray', alpha=0.60)
ax.fill_between( v0, a03, 1e3*np.ones_like(v0), lw=0., color='gray', alpha=0.90)
v0fill = np.linspace(0.,2.242,100)
ax.fill_between( v0fill, np.ones_like(v0fill), 1e3*np.ones_like(v0fill), lw=0., color='gray', alpha=0.90 )
#plt.axvline( 3.025 ) # for 10 times slower nnn than nn tunneling
ax.axvline( 5.088, color='gray', label=r'$t=20\,t_{2}$' ) # for 20 times slower nnn than nn tunneling
from colorChooser import rgb_to_hex, cmapCycle
colors = [ cmapCycle(matplotlib.cm.hot, i, lbound=-1.5, ubound=4.5) for i in range(4) ]
v0U = np.linspace(0.2, 20., 180)
ax.plot( v0U, aSconstInt(v0U,1.), color=colors[3], label=r'$U/t=1$')
ax.plot( v0U, aSconstInt(v0U,5.), color=colors[2], label=r'$U/t=5$')
ax.plot( v0U, aSconstInt(v0U,10.), color=colors[1], label=r'$U/t=10$')
ax.plot( v0U, aSconstInt(v0U,20.), color=colors[0], label=r'$U/t=20$')
v0nnn=np.linspace(0., 5.088)
fill_between(v0nnn, np.ones_like(v0nnn), 1e3*np.ones_like(v0nnn), color="none", hatch="XX", edgecolor="gray", linewidth=0.0)
ax.axhline( 73., color='b', lw=0.5, label=r'$\mathrm{narrow\ Feshbach}$'+'\n'+r'$\mathrm{resonance}$')
ax.axhspan( 560., 1000., color='blue', lw=0.1, alpha=0.4, label=r'$\mathrm{collisional}$'+'\n'+'$\mathrm{losses}$')
ax.axvline( 832. )
l0 = ax.legend( bbox_to_anchor=(1.04,1.035), \
loc='upper left', numpoints=1, \
prop={'size':10}, handlelength=2., handletextpad=0.5 )
#ax.axvline( 7.0 )
#ell = matplotlib.patches.Ellipse(xy=(7., 320.), width=0.2, height=480., angle=0., lw=0. )
ell = matplotlib.patches.FancyBboxPatch(
(7.-0.2,316.-240.), 0.4, 480., lw=0.,
boxstyle=matplotlib.patches.BoxStyle("Round", pad=0.02))
ax.add_artist(ell)
ell.set_clip_box(ax.bbox)
ell.set_facecolor('green')
ell.set_alpha(0.5)
ax.set_xlim(0., 20.)
ax.set_ylim(10, 1000)
ax.set_yscale('log')
ax.set_xlabel(r'$V_{0}\ (E_{r})$', fontsize=fs)
ax.set_ylabel(r'$a_{s}\ (a_{0})$', fontsize=fs)
fig.savefig('BandStructure_figures/singleband.png', dpi=300)
fig.savefig('BandStructure_figures/singleband.pdf')
In [ ]:
class CustomXHatch(matplotlib.hatch.HorizontalHatch):
def __init__(self, hatch, density):
char_count = hatch.count('_')
if char_count > 0:
self.num_lines = int((1.0 / char_count) * density)
else:
self.num_lines = 0
self.num_vertices = self.num_lines * 2
In [215]:
matplotlib.hatch.HatchPatternBase
Out[215]: