Load up the sympy python package to do symbolic mathematics 


```
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'})
```

# Band structure calculation for a 1D Lattice


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.  


```
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)
```

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. 

##Energy eigenvalues in 1D


```
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


```
# 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$


```
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


```
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


```
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


```
# 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


```
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')

```


```
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')
```

## Energy eigenstates for a 1D lattice

To calculate the eigenstates in 1D we diagonlize the coefficient matrix and use the set of coefficient eigenvectors to build up the eigentates as a sum of plane waves

### We can build up the eigenstates using arrays:


```
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


```
# 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')

```

# Calculation of Wannier states

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



## The wannier states are defined below


```
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:


```
from scipy import integrate
x = linspace(-5, 5, 1000)
v0 = 5. 
norm = integrate.simps(np.abs(wannier(x, v0))**2,x)
print norm
```

    0.999999999953


We then plot a few Wannier states for some lattice depths to show how they become increasingly localized


```
# 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


```
# 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')
```

# Band structure for a 3D lattice 


In a 3D optical lattice the hamiltonian is separable in the three cartesian coordinates.  So we can simply add up the energies and form products of the states from the solution of the 1D problem. 


```
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


```
# 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')
```

## Define an interpolation function to quickly calculate the band structure. The interpolation is for a 1D lattice and the 3D bands should be calculated by deciding which band index to use for each axis.


The interpolation data is saved to file for easy retrieval without having to evaluate the entire thing


```
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 


```
from Band3DInterp import *
```

## Compare the performance of the direct calculation and the interpolation


```
%timeit eigenvalsH0(0., 10.)
%timeit interp0[0](10.)
```

    100 loops, best of 3: 7 ms per loop
    10000 loops, best of 3: 28 us per loop



```
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)
```

    1 loops, best of 3: 4.23 s per loop
    100 loops, best of 3: 16.7 ms per loop



```
print "Performance is increased by a factor of %d" % (4230 / 16.7)
```

    Performance is increased by a factor of 253


We can use the interpolation version of the 3D bands calculation to plot the band structure


```
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')
```

# Calculation of tunneling coefficient


## Using the sum over energy eigenvalues


```
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: 


```
v0 = 0.1 
print tCalc(v0 )
v0 = 0.01
print tCalc(v0 )
v0 = 0.001
print tCalc(v0 )
v0 = 0.00001
print tCalc(v0 )
```

    0.20215147731
    0.202994190322
    0.203053269912
    0.203059485994


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.


```
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. 


```
#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:


```
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


```
# 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


```
def tMathieu( v0 ):
    return 4./np.sqrt(np.pi) * v0**(3./4.) * np.exp( -2*np.sqrt(v0))

```


```
# 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')
```

## Using the overlap between Wannier states

With the definition of the Wannier states that was defined above it is easy to compute derivatives of the Wannier functions:

$$ 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) $$ 



$$ \frac{\partial}{\partial x} w^{n} = -\frac{1}{L} \sum_{m>0} \sum_{q>0} c_{qm}^{n} 4\pi(m+q) \sin(2\pi(m+q)x) $$ 

$$ \frac{\partial^{2}}{\partial x^{2}} w^{n} = \frac{1}{L} \sum_{m>0} \sum_{q>0} c_{qm}^{n} 8\pi^{2}(m+q)^{2} \sin(2\pi(m+q)x) $$

We define the second derivative below


```
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


```
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
```


```
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) )
```

    v0	
    ----
    1.00	0.17812	0.17783	0.30542	0.17783
    3.00	0.11103	0.11134	0.16102	0.11103
    7.00	0.03942	0.03954	0.04889	0.03942
    10.00	0.01918	0.01923	0.02274	0.01918


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. 


```
%timeit tWannier(7.0)
%timeit tCalc(7.0)
```

    1 loops, best of 3: 10.5 s per loop
    1 loops, best of 3: 738 ms per loop


# Calculation of on-site interactions


```
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


```
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


```
#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


```
# 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')
```


```
# 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')
```

## Region of validity for the single band Hubbard model


```
# 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.)
```

    5.08836237817



```


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')
    
```


```
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
```


```
matplotlib.hatch.HatchPatternBase
```




    matplotlib.hatch.HatchPatternBase





fd2c18d65dbbd12d8e1236d38a3f5dca5070eb91
ipynb
Jupyter Notebook
BandStructure.ipynb
PedroMDuarte/thesis-lattice
e91d14e2a57af5a6f926511637b456dc6213d634
MIT
2019-05-02T15:52:11.000Z
2021-01-18T13:41:22.000Z
BandStructure.ipynb
PedroMDuarte/thesis-lattice
e91d14e2a57af5a6f926511637b456dc6213d634
MIT
BandStructure.ipynb
PedroMDuarte/thesis-lattice
e91d14e2a57af5a6f926511637b456dc6213d634
MIT
2020-04-22T18:37:35.000Z
2020-07-08T05:49:53.000Z
Qwen/Qwen-72B
1. YES
2. YES