Periodic potential: Bloch states

The objective of this worksheet is to explore a periodic potential that can represent a one-dimensional lattice of potential wells wrapped around to form a circle.

In the Kronig-Penney model of conductivity Bloch's theorem leads to a derivation of energy bands. The eigenvalues are continuous as they depend on a real parameter (the wave propagation vector K). On the other hand it is possible to approximate the continuous energy eigenvalue regions using a finite number of potential wells (we could have Nw=10^24  as in a macroscopic solid). This situation is presented, e.g., in the QM text of Griffiths.

Here we use the cosine-potential familiar from the quantum pendulum. We impose the boundary condition that after Nw  potential wells the wavefunction (not just the probability density) repeats itself. From well to well we should have periodicity of the density (the wavefunctions are allowed to differ by a complex number of modulus 1 according to the Bloch theorem which utilizes the fact that the operator which represents displacement by one cell width commutes with the Hamiltonian.

We use the QuantumPendulum.mws  worksheet as a starting point. It was written in such a way that other potential energy assemblies than the cosine potential could be implemented for the eigenvalue calculations. We use, however, the specific cosine potential for the construction of the eigenfunctions, as this permits us to use the known solutions to the Mathieu differential equation.

>	restart; Digits:=15: with(LinearAlgebra):

>

a:=Pi;

>

V0:=2;

>	V:=x->V0*(1-cos(Pi*x/a));

>	Vmax:=V(a);

We pick an odd number of wells:

>

Nw:=3;

>	Vmw:=piecewise(seq(op([x>a*(2*j-1) and x<a*(2*j+1),V(x-2*a*j)]),j=-(Nw-1)/2..(Nw-1)/2));

>	PTS:=[seq((2*j-1)*a,j=-(Nw-1)/2..(Nw+1)/2)];

We want to have periodicity with 2a

>	Pmw:=plot(Vmw,x=PTS[1]..PTS[nops(PTS)],numpoints=500):

A periodic basis allows for left- and right-travelling states with periodicity conditions at a , and -a .

>	psi:=n->eval(exp(I*Pi*n*x/(Nw*a)));

>	int(conjugate(psi(0))*psi(2),x=PTS[1]..PTS[nops(PTS)]);

>	L:=int(conjugate(psi(1))*psi(1),x=PTS[1]..PTS[nops(PTS)]);

The kinetic energy for the plane-wave basis functions can be calculated in closed form. As an example we list the one for n=3 :

>	int(-diff(psi(3),x$2)*conjugate(psi(3)),x=PTS[1]..PTS[nops(PTS)])/L;

>	4*3^2*Pi^2/(L)^2;

Here we select the basis size. Note that the index labeling the plane-wave basis states ranges over positive and negative integers n  (corresponding to right-travelling and left-travelling waves).

>	n_max:=20; N:=2*n_max+1;

>	VME:=unapply(int(V(x)*exp(Pi*I*dn*x/(Nw*a)),x=x1..x2),x1,x2,dn);

>	#VME(PTS[1],PTS[2],21);

>	#limit(VME(PTS[1],PTS[2],s),s=0);

>	kh:=1; # for the cosine potential we put the strength into V(x).

>	H:=Matrix(N,N,shape=hermitian,datatype=complex[8]):

We should have a 1/sqrt(2*a) = 1/sqrt(L)  normalization factor in front of the plane-wave basis states! We can fix this by dividing the matrix by L=Nw*2a .

>	for i1 from 1 to N do: n1:=-n_max+i1-1: for i2 from 1 to i1 do: n2:=-n_max+i2-1: if i2=i1 or n1-n2=Nw then

>	H[i1,i2]:=evalf(1/L*kh*add(limit(VME(PTS[i0],PTS[i0+1],s),s=n1-n2),i0=1..Nw)): else H[i1,i2]:=evalf(1/L*kh*add(VME(PTS[i0],PTS[i0+1],n1-n2),i0=1..Nw)): fi: od: H[i1,i1]:=H[i1,i1]+evalf(4*n1^2*Pi^2/(L)^2): od:

>	Eigenvalues(H);

>	EV:=convert(%,list);

Compared to the single-well energy spectrum we observe that the lowest oscillator states are approximately replicated Nw  times. There is a unique ground state, and a pair of eigenvalues that appear to be degenerate.

We could construct the eigenfunctions from the Eigenvectors command. The calculations with these states are tedious (slow). We prefer to connect with the Mathieu differential equation. This requires us to rescale xi=2*x*Nw .

>	q:=2*V0*Nw^2;

>	MDE:=a->diff(u(x),x$2)+(a+2*q*cos(2*Nw*x))*u(x)=0;

Here a  is the eigenvalue scaled according to a=4*(EV-V0) , and the potential strength parameter q=2*V0 .

>	aM:=e->Nw^2*4*(e-V0);

>	sol:=dsolve(MDE(aM(EV[1])));

>	sol:=u(x) = _C1*MathieuC(-77049338729459/18000000000000,-4,3*x)+_C2*MathieuS(-77049338729459/18000000000000,-4,3*x);

The MathieuC  and MathieuS  functions are symmetric and antisymmetric respectively, i.e., behave like sine and cosine. We invoke a physics argument which states that the eigenstates of the problem are either symmetric or antisymmetric, as the Hamiltonian is symmetric. The ground state is symmetric.

>	plot([eval(rhs(sol),{_C1=1,_C2=0}),eval(rhs(sol),{_C2=1,_C1=0})],x=-Pi/2..Pi/2,color=[red,blue],view=[-Pi/2..Pi/2,-q/10..q/10]);

The antisymmetric solution is not acceptable, it is not periodic. The ground state is unique, it has a symmetric eigenfunction.

Let us plot a probability using the eigenenergy as a baseline, and superimpose it on the potential as a function of theta.

>	P0:=plot([EV[1],EV[1]+(eval(rhs(sol),{_C1=1,_C2=0,x=s/(2*Nw)}))^2],s=-Nw*Pi..Nw*Pi,color=[black,green]):

>	plots[display](Pmw,P0);

>	sol:=dsolve(MDE(aM(EV[2])));

>	sol:=u(x) = _C1*MathieuC(-191905256659097/45000000000000,-4,3*x)+_C2*MathieuS(-191905256659097/45000000000000,-4,3*x);

>

We might expect the first excited state is anti-symmetric. We find, however, since the 2nd and 3rd eigenenergies are so close that there have to be 2 eigenfunctions at this energy. They are the symmetric and anti-symmetric functions.

>	plot([eval(rhs(sol),{_C1=1,_C2=0}),eval(rhs(sol),{_C2=1E-2,_C1=0})],x=-Pi/2..Pi/2,color=[red,blue]);

Is it true that both the symmetric and the antisymmetric solutions work? Let us check the periodicity:

>	fS:=eval(rhs(sol),{_C1=1,_C2=0});

>	fA:=eval(rhs(sol),{_C1=0,_C2=1E-2});

>	evalf(limit(fS,x=-Pi/2));

>	evalf(limit(fS,x=Pi/2));

>	evalf(limit(diff(fS,x),x=-Pi/2));

>	evalf(limit(diff(fS,x),x=Pi/2));

The symmetric solution works! It is periodic to a high degree of precision.

>	evalf(limit(fA,x=-Pi/2));

>	evalf(limit(fA,x=Pi/2));

>	evalf(limit(diff(fA,x),x=-Pi/2));

>	evalf(limit(diff(fA,x),x=Pi/2));

The anti-symmetric solution also works! This is a feature that was observed for the rotator states in the single-well problem; in fact the lowest rotator states didn't display it.

At this point we should invoke physics as a guide. On the one hand we might argue that there is a symmetric and an antisymmetric state. The antisymmetric state, e.g., violates the expectation that the probability density is the same in all three wells, as the middle one comes out practically empty:

>	P1:=plot([EV[2],EV[2]+(eval(rhs(sol),{_C1=0,_C2=1E-2,x=s/(2*Nw)}))^2],s=-Nw*Pi..Nw*Pi,color=[black,magenta]):

>	plots[display](Pmw,P0,P1);

We can also argue that since the energy levels are degenerate, any linear combination of the two eigenfunctions found (and associated with this energy value) is also an acceptable eigenfunction.

In particular, we can construct a left-travelling and a right-travelling solution by using +I  and -I  in front of the anti-symmetric state, and a real constant in front of the symmetric one. We apply a trial-and-error strategy to ensure an equal peak height over the three wells for the probability density.

>	eval(fS,x=0);

>	f_Right:=fS+I*0.435*fA;

>	P1:=plot([EV[2],EV[2]+(eval(abs(f_Right)^2,{x=s/(2*Nw)}))^2],s=-Nw*Pi..Nw*Pi,color=[black,magenta]):

>	plots[display](Pmw,P0,P1);

To prove that these states represent travelling electrons (around the ring made up of three sites) we have to calculate the probability current density (we use units with mass m=1 ).

>	jPC:=psi->-I/2*(conjugate(psi)*diff(psi,x)-psi*diff(conjugate(psi),x));

>	#evalf(limit(jPC(f_Right)/f_Right,x=0.1));

Maple is having a hard time. Let's calculate the derivative numerically:

>	eta:=1E-3;

>	psi_d:=x0->(evalf(eval(f_Right,x=x0+eta)-eval(f_Right,x=x0-eta)))/2/eta;

Based on the current density we define the average local velocity in QM by dividing the probability current by the probability:

>	v_loc:=proc(x0) local valf,valf_d: global f_Right; valf:=evalf(eval(f_Right,x=x0));

>	valf_d:=psi_d(x0);

>	-I/2*(conjugate(valf)*valf_d-conjugate(valf_d)*valf)/(conjugate(valf)*valf); end:

>	v_loc(0.1);

>	plot(Re(v_loc(s/(2*Nw))),s=-Nw*Pi..Nw*Pi);

The local velocity is very high when the particle is tunneling through one of the barriers! What is it near the minimum in the potential?

>	plot(Re(v_loc(s/(2*Nw))),s=-Nw*Pi..Nw*Pi,view=[-Nw*Pi..Nw*Pi,0..0.1],numpoints=500);

This may seem strange: the quantum local velocity is lowest near the potential minima where a classical particle would be moving fastest.

Question:  does the high local velocity inside the potential barrier have anything to say about the tunneling time?

Concerning the interpretation of the local quantum velocity keep in mind that it is zero for oscillators (the eigenfunction is real-valued). It describes net flow.

In the atomic physics context: s-states and generally m=0  states are real-valued and have no local quantum velocity associated with them. The m>0  states have a net current flow about the z-axis in one direction, while the m<0  states have a net current in the opposite directions giving rise to magnetic moment interactions. In that sense these states involve current loops. Likewise, in the present problem there is net flow to the right (or to the left if the sign of the imaginary part of the wavefunction is reversed).

Exercise 1:

Compile information about the solutions with more wells (e.g., 5). Based on the case of Nw=3  make predictions concerning the nature of the states and then proceed to check out these predictions.

We now proceed with a demonstration how energy bands emerge when the potential chain is made up of more wells. We just compute the energy eigenvalues. The code is just a copy of the eigenvalue calculation provided above. When increasing the number of wells we get more replicas of the near-degenerate ground states. Thus, we have to increase the matrix size to obtain sufficiently many (and also accurate) eigenenergies.

>

>	restart; Digits:=15: with(LinearAlgebra):

>

a:=Pi;

We selected a standard interval for the dynamical variable, i.e., [-Pi,Pi) . The re-scaled potential energy is chosen such that in the oscillator limit the eigenenergies are given as the odd integers.

>

V0:=2;

>	V:=x->V0*(1-cos(Pi*x/a));

>	Vmax:=V(a);

We pick an odd number of wells:

>

Nw:=11;

>	Vmw:=piecewise(seq(op([x>a*(2*j-1) and x<a*(2*j+1),V(x-2*a*j)]),j=-(Nw-1)/2..(Nw-1)/2));

>	PTS:=[seq((2*j-1)*a,j=-(Nw-1)/2..(Nw+1)/2)];

>	Pmw:=plot(Vmw,x=PTS[1]..PTS[nops(PTS)],numpoints=500):

>	psi:=n->eval(exp(I*Pi*n*x/(Nw*a)));

>	L:=int(conjugate(psi(1))*psi(1),x=PTS[1]..PTS[nops(PTS)]);

Here we pick the matrix size:

>	n_max:=50; N:=2*n_max+1;

>	VME:=unapply(int(V(x)*exp(Pi*I*dn*x/(Nw*a)),x=x1..x2),x1,x2,dn);

>	H:=Matrix(N,N,shape=hermitian,datatype=complex[8]):

We should have a 1/sqrt(2*a) = 1/sqrt(L)  normalization factor! We can fix this by dividing the matrix by L=Nw*2a .

>	for i1 from 1 to N do: n1:=-n_max+i1-1: for i2 from 1 to i1 do: n2:=-n_max+i2-1: if i2=i1 or n1-n2=Nw then

>	H[i1,i2]:=evalf(1/L*add(limit(VME(PTS[i0],PTS[i0+1],s),s=n1-n2),i0=1..Nw)): else H[i1,i2]:=evalf(1/L*add(VME(PTS[i0],PTS[i0+1],n1-n2),i0=1..Nw)): fi: od: H[i1,i1]:=H[i1,i1]+evalf(4*n1^2*Pi^2/(L)^2): od:

>	Eigenvalues(H):

>	EV:=convert(%,list);

>	PEV:=plot([seq(EV[i],i=1..N)],x=PTS[1]..PTS[nops(PTS)],color=blue):

>	plots[display](Pmw,PEV, title="periodic stitching at the ends",view=[-L/2..L/2,0..1.2*Vmax]);

We observe that for the energy levels that correspond to excitations in the single-well problem the spread of eigenvalues is larger than for the ground state.

An investigation of the energy levels around E=2.6  shows that there are Nw=11  states which come as 5 degenerate pairs plus one extra state. The extra state is at the top of the 'band' in this case.

For each energy the wavefunctions are delocalized. One can say that tunneling plays an important role all the time.

Even at energies above the potential barriers the energy spectrum displays a band structure with Nw  states per band. For a (one-dimensional) macroscopic lattice the number Nw  is huge and the energy bands are practically continuous. The problem of conductivity (metals, insulators, semiconductors) can be understood in terms of this band structure.

Exercise 2:

Explore the band structure by increasing the number of wells. Will the gaps fill in as a result of increasing the number of wells Nw ? Make sure the matrix size is sufficiently large to yield enough eigenvalues and to ensure their accuracy.

>

>

>

Summary:

The Bloch theorem strictly speaking is applicable for an ifinite sequence of potential wells, in which case continuous energy bands emerge.

For a periodic chain of potential wells we can construct eigenstates that resemble the Bloch states in the sense that the degenerate energy levels allow eigenstates which have identical density over each well, and which correspond to net flow of particle density to the right or to the left.

The periodic boundary condition through which the lattice is stitched to form a ring conflicts with the Bloch condition. This is (probably) the reason for the possibility of eigenstates that don't agree with the Bloch theorem with densities that aren't necessarily periodic. For large Nw  the violation of the Bloch condition should not matter.