Quantum tunneling through an arbitrary potential barrier

We are interested in solving the Schroedinger wave equation with scattering boundary conditions for an arbitrary potential function in one dimension. The boundary conditions (incident and reflected and transmitted wave) require the solution to be complex-valued. Our attitude will be to apply knowledge about general solution to a second-order ordinary differential equations (two independent solutions are required, and then matched to the physically imposed boundary conditions).

Textbook treatments begin with simple potential shapes (square barrier) for which the solution can be obtained exactly. For a symbolic solution along these lines we refer the reader to our book [Marko Horbatsch: Quantum Mechanics using Maple , chapter 3.2 (Springer 1995)]. Here we rely on numerical techniques included in Maple, namely to solve initial-value problems for systems of first-order differential equations. The tunneling problem requires us to solve a boundary-value problem. This will be achieved by generating two independent unphysical solutions based on inital values, and then to match a linear combination of these to the physical boundary conditions.

Maple6 has added capability to dsolve[numeric] to handle complex-valued differential equations. This permits to solve the Schroedinger equation for scattering boundary conditions directly. To remain compatible with Maple 5 we stay, however, with real-valued boundary conditions for the solution of the differential equation, and introduce the complex-valued nature of the physical solution in the step where linear combinations of two fundamental solutions are formed.

We could begin with simple barrier penetration. However, we assume that the reader is already familiar with the basics of tunneling. Therefore, we proceed with the more interesting double-barrier. An example of a double-barrier potential is:

> V:=exp(-5*x^2)+exp(-5*(x-2)^2);

> plot(V,x=-2..5);

We are interested in discovering so-called transmission resonances. We can safely assume that the interaction vanishes at x =-2 and x =4 for our example.

We work in units where hbar= m =1. This simplifies relationships between momentum and wavenumber (they become identical).

The wave number (momentum) can be obtained from the non-relativistic energy as:

> kE:=En->sqrt(2*En);

We need to propagate two independent solutions, and then to apply matching conditions. The Schroedinger wave equation is given as:

> SE:=En->-1/2*diff(u(x),x$2)+(V-En)*u(x)=0;

The boundary condition states that:

1) From the left we have an incident free plane wave with wavenumber k and of known amplitude, and a reflected free plane wave with unknown amplitude (whose magnitude squared will eventually become the reflection probability; (the choice of side from which the wave comes in is arbitrary).

2) To the right of the potential region we have an outgoing free plane wave with wavenumber k and unknown magnitude. The squared magnitude of the transmitted wave becomes the transmission probability (if the incident wave is normalized to unity).

We know from (2) that the transmitted flux goes like T exp(i kx ). Thus, the ratio of the derivative of the physical solution to the solution itself has to equal i k .

Therefore, the strategy is to generate two independent solutions, and to form the linear combination that satisfies the requested relationship.

The initial and final values for the independent varible range in the numerical integration of the differential equation are:

> x_i:=-2; x_f:=4;

We pick an energy value:

> En:=1/2;

An aside:

Experts can fix an oversight in Maple6.0, if they wish to use complex-valued boundary conditions:

For complex-valued solutions in Maple6 override: (Robert Israel contributed this fix to MUG).

`dsolve/numeric/init_y0`:=subs();numeric=complex(numeric),eval(`dsolve/numeric/init_y0`));

We pick two intial conditions that provide independent solutions:

> IC1:=u(x_i)=1,D(u)(x_i)=1;

> IC2:=u(x_i)=1,D(u)(x_i)=0;

> sol1:=dsolve({SE(En),IC1},u(x),numeric,output=listprocedure):

> u1:=subs(sol1,u(x)):

> u1p:=subs(sol1,diff(u(x),x)):

To demonstrate that it calculates something...

> u1(x_f);

> u1p(x_f);

Now the second solution:

> sol2:=dsolve({SE(En),IC2},u(x),numeric,output=listprocedure);

> # not needed as we picked real-valued initial conditions #`dsolve/numeric/init_y0`:=subs(numeric=complex(numeric),eval(`dsolve/numeric/init_y0`));

> #sol2:=dsolve({SE(En),IC2},u(x),numeric,output=listprocedure):

> u2:=subs(sol2,u(x)):

> u2p:=subs(sol2,diff(u(x),x)):

> u2(x_f);

> u2p(x_f);

We generated the two independent solutions using real boundary conditions.

Now we assemble the complex-valued solution that satisfies the physical boundary condition of only a transmitted wave existing to the right of x_f. We need a single mixing coefficient, and determine it by imposing that the ratio of derivative over the solution equals i k .

> k:=kE(En);

> c2:=solve((u1p(x_f)+c2*u2p(x_f))/(u1(x_f)+c2*u2(x_f))=I*k,c2);

Now we can form the physical solution, which has arbitrary normalization at this point. We graph the real and imaginary parts:

> P1:=plot('Re(u1(x)+c2*u2(x))',x=-12..14,color=blue):

> P2:=plot('Im(u1(x)+c2*u2(x))',x=-12..14,color=red):

> plots[display](P1,P2);

This is a solution from which we should be able to extract the physical information. To the right of the potential we have obviously a travelling plane wave. It is travelling to the right with momentum k (we recognize the periodicity as about 2 Pi, corresponding to k =1, and we see that the imaginary part precedes the real part, corresponding to motion to the right) and we can read off the amplitude T .

On the left we must have the sum of an incident and reflected wave. This is characterized by unequal amplitudes for real and imaginary parts.

In the middle we have a complicated connecting piece. This is where dsolve[numeric] and forming the linear combo has done the work for us!

So far the incident wave has arbitrary normalization. Let us see whether we can figure the normalization out from the solution in region I (at x=x_i).

We match now at the x_i a solution which has independent amplitudes for the incoming and reflected waves, named N and R respectively:

> eq1:=N*exp(I*k*x_i)+R*exp(-I*k*x_i)=u1(x_i)+c2*u2(x_i);

> eq2:=I*(k*N*exp(I*k*x_i)-k*R*exp(-I*k*x_i))=u1p(x_i)+c2*u2p(x_i);

> sol:=solve({eq1,eq2},{N,R});

> assign(sol);

We observe that N and R are amplitudes, and turn out to be complex-valued.

The properly normalized reflection coefficient now becomes:

> r:=R/N;

and the transmission coefficient:

> t:=evalf((u1(x_f)+c2*u2(x_f))/exp(I*k*x_f)/N);

The check whether this all makes sense is to sum the magnitude squared for reflection and transmission:

> abs(r)^2+abs(t)^2;

Within numerical accuracy we preserve the norm; therefore we interpret as reflection and transmission probabilities:

> abs(r)^2,abs(t)^2;

Now we turn the calculation into a procedure: for given energy it calculates the reflection and transmission coefficients.

Our previously developed commands are used step-by-step as before:

> Tunnel:=proc(En) local IC1,IC2,sol1,sol2,u1,u1p,u2,u2p,k,c2,N,R,r,t,sol,eq1,eq2; global V,kE,x_i,x_f;

> IC1:=u(x_i)=1,D(u)(x_i)=1;

> IC2:=u(x_i)=1,D(u)(x_i)=0;

> sol1:=dsolve({SE(En),IC1},u(x),numeric,output=listprocedure);

> u1:=subs(sol1,u(x));

> u1p:=subs(sol1,diff(u(x),x));

> sol2:=dsolve({SE(En),IC2},u(x),numeric,output=listprocedure);

> u2:=subs(sol2,u(x));

> u2p:=subs(sol2,diff(u(x),x));

> k:=kE(En);

> c2:=solve((u1p(x_f)+c2*u2p(x_f))/(u1(x_f)+c2*u2(x_f))=I*k,c2);

> N:='N': R:='R':

> eq1:=N*exp(I*k*x_i)+R*exp(-I*k*x_i)=u1(x_i)+c2*u2(x_i);

> eq2:=I*(k*N*exp(I*k*x_i)-k*R*exp(-I*k*x_i))=u1p(x_i)+c2*u2p(x_i);

> sol:=solve({eq1,eq2},{N,R});

> assign(sol);

> r:=R/N;

> t:=evalf((u1(x_f)+c2*u2(x_f))/exp(I*k*x_f)/N);

> #print(" R= ",abs(r)^2," T= ",abs(t)^2," R+T= ",abs(r)^2+abs(t)^2," norm equal to 1.0 ?");

> [abs(r)^2,abs(t)^2]; end:

We check it for the previously tested energy value:

> Tunnel(1/2);

Now we are ready to scan the energy range from zero to above the peak height of the barrier. The results are accumulated in lists for subsequent plotting.

> Tr:=[]: Rf:=[]: for iE from 1 to 60 do: En:=iE/20; res:=Tunnel(En); Tr:=[op(Tr),[En,res[2]]]: Rf:=[op(Rf),[En,res[1]]]: od:

> P1:=plot(Rf,color=red): P2:=plot(Tr,color=blue):

> plots[display](P1,P2,labels=["E","{t,r}"]);

Observe the interesting features:

1) For sufficiently high energies only transmission occurs. Note, however, that classically this situation would be reached for E >1!

2) For E =0.66 perfect transmission and no reflection is observed. This is the phenomenon of a transmission resonance. The resonance structure has a width associated with it. Most physical phenomena that have a lifetime associated with them (unstable elementary particles, alpha decay, nuclear gamma lines, atomic optical transitions, excitations in semiconductors, etc.) can be understood as processes where the wave equation describing the matter (or radiation) field involves a tunneling process. Broad resonances are associated with short lifetimes, sharp resonances with long lifetimes, according to Heisenberg's energy-time uncertainty relation.

The transmission resonance is quite broad. We can make the walls of the two barriers thicker, and observe what the effect is on the transmission probability.

> V:=exp(-8*x^6)+exp(-8*(x-2)^6);

> plot(V,x=-2..5);

> Tr:=[]: Rf:=[]: for iE from 1 to 120 do: En:=iE/40; res:=Tunnel(En); Tr:=[op(Tr),[En,res[2]]]: Rf:=[op(Rf),[En,res[1]]]: od:

> P1:=plot(Rf,color=red): P2:=plot(Tr,color=blue):

> plots[display](P1,P2,labels=["E","{t,r}"]);

The transmission resonance is sharper now, and occurs just below the classical border value for the energy.

Exercise:

Observe what happens when the gap is widened between the potential humps. Be careful with the values of x_i, x_f, when changes are made to the potential.

We can include an option in Tunnel to graph the solution. Perhaps this will help to understand why reflection and transmission comes about.

> Tunnel:=proc(En,iplot) local IC1,IC2,sol1,sol2,u1,u1p,u2,u2p,k,c2,N,R,r,t,sol,eq1,eq2; global V,kE,x_i,x_f,PL;

> IC1:=u(x_i)=1,D(u)(x_i)=1;

> IC2:=u(x_i)=1,D(u)(x_i)=0;

> sol1:=dsolve({SE(En),IC1},u(x),numeric,output=listprocedure);

> u1:=subs(sol1,u(x));

> u1p:=subs(sol1,diff(u(x),x));

> sol2:=dsolve({SE(En),IC2},u(x),numeric,output=listprocedure);

> u2:=subs(sol2,u(x));

> u2p:=subs(sol2,diff(u(x),x));

> k:=kE(En);

> c2:=solve((u1p(x_f)+c2*u2p(x_f))/(u1(x_f)+c2*u2(x_f))=I*k,c2);

> N:='N': R:='R':

> eq1:=N*exp(I*k*x_i)+R*exp(-I*k*x_i)=u1(x_i)+c2*u2(x_i);

> eq2:=I*(k*N*exp(I*k*x_i)-k*R*exp(-I*k*x_i))=u1p(x_i)+c2*u2p(x_i);

> sol:=solve({eq1,eq2},{N,R});

> assign(sol);

> r:=R/N;

> t:=evalf((u1(x_f)+c2*u2(x_f))/exp(I*k*x_f)/N);

> if iplot=2 then print(" R= ",abs(r)^2," T= ",abs(t)^2," R+T= ",abs(r)^2+abs(t)^2," norm equal to 1.0 ?") fi;

> if iplot=1 then

> PL:=plot([V,'Re((u1(x)+c2*u2(x))/N)','Im((u1(x)+c2*u2(x))/N)'],x=x_i-10..x_f+10,color=[red,blue,green],title=cat("Energy= ",convert(evalf(En,3),string))):

> fi;

> [abs(r)^2,abs(t)^2]; end:

> Tunnel(0.8,1);

> plots[display](PL);

Can we run an animation as a function of energy? It takes a while to compute, but it can be done.

> PPL:='PPL': for iE from 1 to 120 do: En:=iE/40; res:=Tunnel(En,1); PPL[iE]:=PL: od:

> plots[display](seq(PPL[i],i=1..120),insequence=true);

The possibilities to explore scattering from different potentials with this worksheet appear to be unlimited! Transmission resonances appear to play a significant role in newly developed semiconductor devices, and may prove to provide opportunities for future generations of computing hardware.

>