import numpy as np
import pandas as pd
import ipywidgets as ipw

from scipy.interpolate import interp1d
from scipy.linalg import expm
from bokeh.palettes import RdYlBu11
from bokeh.transform import linear_cmap
from bokeh.io import push_notebook, output_notebook, show
from bokeh.layouts import row, column
from bokeh.plotting import figure, curdoc
from bokeh.models import ColumnDataSource, ColorBar, BasicTickFormatter, BasicTicker, FixedTicker, FuncTickFormatter
from bokeh.client import push_session

from ipywidgets import interact

import matplotlib.pyplot as plt
import matplotlib.animation as anim
output_notebook()

Loading BokehJS ...

Wavepacket construction and propagation along the 1D chain

Introduction

In this notebook, we demonstrate the wavepacket construction and propagation, which is at the heart of the semiclassical approximation. For simplicty, we will consider the familiar 1D chain described using the tight-binding Hamiltonian. We will have to take a sufficiently large number of sites $N$ and impose periodic boundary conditions. For simplicity, in the code below we assume lattice constant $a = 1$.

N = 50; # Number of atomic sites
M = int((N)/2); # Number of atomic sites on each side of the origin
a = 1.0 # Real space distance between each site
lattx = np.arange(-a*M,a*M,a); # 
latty = np.zeros(lattx.size);
lattice = figure(height=100, width=700, title="1D atomic chain",x_range=[-1.1*a*M,1.1*a*M]);
source1 = ColumnDataSource(data = {'xVal': lattx, 'yVal': latty});
lattice.scatter('xVal','yVal',source = source1, size = 6, line_color = 'black', fill_color = 'black');
lattice.xaxis.axis_label = 'r (a)';
show(lattice)

Eigenvalues

The solution for the 1D chain with one atom per unit cell is known since the previous lectures:

$$ E(k) = -2\gamma \cos (ka) $$

The plot below reproduces the corresponsing band structure. Each of the $N$ points on the band corresponds to an eigenvalue of the $N$-site chain with periodic boundary (Born-von karman) condition $\psi(r) = \psi(r + Na)$. We will also assume $\gamma = -1$ in the code for simplicity.

k0 = 2*np.pi/(N*a); # Calculate the smallest allowed k-vector
k = np.linspace(-M*k0,(M-1)*k0,N);
energy = -2*np.cos(k);

# Calcualte resulting dispersion relation using periodic boundary condition (PC)
band = figure(height = 400, width = 700, title='Band structure of the 1D chain and the eigenvalues of the N-site chain', x_range = [-1.1*k0*M/np.pi,1.1*k0*M/np.pi], y_range= [1.2*energy.min(),1.2*energy.max()]);
source2 = ColumnDataSource(data = {'xVal': k/np.pi, 'yVal' : energy});
band.scatter('xVal','yVal',source = source2,size = 6);
band.line('xVal','yVal',source = source2);
band.xaxis.axis_label='k ( \u03c0 / a )';
band.yaxis.axis_label='Energy/|gamma|';
show(band)

# define the initial wavefunction (with an artificial Gaussian dist. in k-space)
sigm = 3.0*k0; # Distribution width
kk0 = 12.0*k0; # Distribution center
kk = np.linspace(-M*k0,(M-1)*k0, 200); # Create a smoother grid for illustration purpose
phi_k0 = 1/(np.sqrt(2.0*np.pi)*sigm)*np.exp(-1/2*( (kk-kk0)/sigm )**2);

# Initialize an electron at the center of the 1D chain with a gaussian distribution to indicates its probability density
dist0 = figure(height = 400, width = 700, title='Initial wavefunction in real space', x_range = [-1.1*k0*M/np.pi,1.1*k0*M/np.pi], y_range= [1.2*energy.min(),1.2*energy.max()]);
source3 = ColumnDataSource(data = {'xVal': kk/np.pi, 'yVal': phi_k0/phi_k0.max()});
dist0.line('xVal','yVal', source = source3, line_color = 'red');
dist0.scatter('xVal','yVal',source = source2, size = 6, line_color = 'black', fill_color = 'black');
dist0.xaxis.axis_label='k-space distance (arb. unit)';
dist0.yaxis.axis_label='Probability density';
show(dist0)

# Construct hamiltonian using occupation number basis and set initial wavefunction
asd = np.zeros((N-1,N-1),complex);
np.fill_diagonal(asd,-1);
b = np.append(asd,np.zeros((N-1,1),complex), axis = 1);
c = np.append(np.zeros((1,N),complex),b,axis = 0);
hamil = c + np.transpose(c);
np.fill_diagonal(hamil,0);
psi0 = np.zeros((N,1),complex);
rho0 = np.zeros((N,1), float);

phi_k0 = 1/(np.sqrt(2*np.pi)*sigm)*np.exp(-1/2*( (k-kk0)/sigm )**2); # Calculate the initial wavefunction using really k points ('k' instead of 'kk').
for nn in range(0,psi0.size):
    psi_temp = 0.0;
    for kn in range(k.size):
        psi_temp = psi_temp + phi_k0[kn]*np.exp(-1j*(k[kn]*(nn-M)*a));
    psi0[nn] = psi_temp;
    rho0[nn] = np.abs(psi0[nn]);

hamil[N-1][0] = -1; #Apply periodic boundry condition
hamil[0][N-1] = -1;

# Define time evolution operator
dt = 0.1; # Define time step (has to be small if using asymptotic exponential for evolution operator)
endt = 5000;

# Option 1: Use first order Taylor expansion to approximate the evolution operator
diaganol = np.ones(N);
evol = np.diag(diaganol)+1.j*hamil*dt; # first order Taylor expansion

## Option 2: Use Pade approximation in scipy package to approximate the evolution operator
# evol = expm(1.j*hamil*dt);

ts = np.arange(0,endt*dt,dt);
psis = np.empty(ts.size,dtype = object); # Initiate a container for the wavefunctions at different time slice
phis = np.empty(ts.size,dtype = object); # Initiate a container for the phase at different time slice
rhos = np.empty(ts.size,dtype = object); # Initiate a container for the probability density at different time slice
psis[0]=psi0;
for ii in range(0,ts.size-1):
    psis[ii+1] = np.dot(evol,psis[ii]);
for ii in range(0,ts.size-1):
    phis[ii] = np.angle((psis[ii]));
    rhos[ii] = np.abs(psis[ii]);

# Plot the phase of the wavefunction at time = t*dt
t = 0;
phi_t = np.squeeze(phis[t]);
rho_t = interp1d(lattx,np.squeeze(rhos[t]), kind = 'cubic');
wavefunc_t = figure(height=300, width=700, title="Electron wavefunction",x_range=[-1.1*a*M,1.1*a*M]);
source4 = ColumnDataSource(data = {'xVal': lattx, 'yVal' : latty, 'cVal': phi_t});
lattx_smooth = np.linspace(-a*M,a*(M-1),num = 200, endpoint = True);
source5 = ColumnDataSource(data = {'xVal_s': lattx_smooth, 'yVal_s': rho_t(lattx_smooth)/rho_t(lattx_smooth).max()});

#phis_smooth = np.empty(ts.size,dtype = object); # Initiate a container for the phase at different time slice
#for ii in range(0,ts.size-1): # Option 1: use smoothed curve to show phase variation
#    phis_smooth[ii] = interp1d(lattx,np.squeeze(phis[ii]), kind = 'cubic');
#wavefunc_t.line('xVal','yVal', source = source4); 

mapper = linear_cmap(field_name = 'cVal', palette = RdYlBu11, low = min(phi_t), high = max(phi_t));
wavefunc_t.scatter('xVal','yVal',source = source4, size = 6, line_color = mapper, fill_color = mapper);
wavefunc_t.line('xVal_s','yVal_s',source = source5, line_color = 'red');
color_bar = ColorBar(color_mapper = mapper['transform'], width = 300, location = (150,0), orientation = "horizontal");
wavefunc_t.add_layout(color_bar, 'below');
show(wavefunc_t)

def update(time):
    time = int(time/dt);
    rho_t = interp1d(lattx,np.squeeze(rhos[time]), kind = 'cubic');
    source4.data['cVal'] = np.squeeze(phis[time]);
    source5.data['yVal_s'] = rho_t(lattx_smooth)/rho_t(lattx_smooth).max();
    
    push_notebook(handle = handle)

handle = show(wavefunc_t, notebook_handle = True);

interact(update, time = ipw.FloatSlider(min = 0, max = (endt-2)*dt, step = dt, value = 0));

interactive(children=(FloatSlider(value=0.0, description='time', max=499.8), Output()), _dom_classes=('widget-…

Questions

When answering the following questions, please think of the answer first and then modify and run the code

• What happens upon changing the sign of $k_0$?
• Explain the behavior of wavepacket when $k_0 = 0$.
• Explain the changes of the wavepacket at sufficient large $t$.
• How does the evolution of the wavepacket change if periodic boundary condition are removed (i.e. break the bond between sites 1 and $N$)?
• What happens if a different on-siter energy $E$ is set for one of the sites?

# # This animation is not functioning due to the current incompatibility between bokeh and tornado, should be resolved in furture updates.

# def auto_update():
#     global time;
#     source_tt.data['cVal'] = np.squeeze(phis[t])
#     source_tt.data['yVal'] = interp1d(lattice,np.squeeze(rhos[t]), kind = 'cubic')
#     source_tt.trigger('data', source_tt.data, source_tt.data)
#     if time < endt.size - 1:
#         time = time + 1
#     else:
#         time = 0
#         source_tt.data['cVal'] = np.squeeze(phis[0])
#         source_tt.data['yVal'] = interp1d(lattice,np.squeeze(rhos[0]), kind = 'cubic')

# # This animation is not functioning due to the current incompatibility between bokeh and tornado, should be resolved in furture updates.

# session = push_session(curdoc())
# time = 0; # set initial frame
# source_tt = ColumnDataSource(data = {'xVal': lattx, 'yVal' : rhos[0], 'cVal': np.zeros(lattx.size)})
# wavefunc_tt = figure(plot_height=300, plot_width=700, title="Evolution of electron wavefunction",x_range=[-1.1*a*M,1.1*a*M]);

# mapper = linear_cmap(field_name = 'cVal', palette = Spectral6, low = 0, high = 1);
# wavefunc_tt.line(x = 'xVal', y = 'yVal', source = source_tt, line_color = 'red');
# wavefunc_tt.circle('xPos','yPos',source = source4, size = 10, line_color = mapper, fill_color = mapper);

# curdoc().add_periodic_callback(auto_update, 67)

# session.show()
# session.loop_until_closed()