Application Center - Maplesoft

# QUANTUM SCATTERING BY THE ONE-DIMENSIONAL POTENTIAL STEP IN MAPLE

You can switch back to the summary page by clicking here.

QUANTUM SCATTERING BY THE ONE-DIMENSIONAL
POTENTIAL STEP IN MAPLE

Alexei V. Tikhonenko

General Physics Department,
State Technical University of Nuclear Power Engineering,
Obninsk, Russia

Scattering problem of quantum particles by the one-dimensional potential step is solved in MAPLE.
Analytical formulas for reflection and transmission coefficients are obtained and visualized.

1. One-dimensional Schrodinger equation

 > restart:

Consider one-dimensional Schrodinger equation  [1]

 > SchrodingerEQN:=diff(psi(x),x\$2)+2*m/h^2*(E-U(x))*psi(x);

 (1.1)

where is wave function, m is mass and E is energy of particle; is Planck's constant .

This equation describes one-dimensional motion of quantum particles in the field of potential barrier :

 > U(x):=U0/(1+exp(-alpha*x));

 (1.2)

We will use notations:

 > EQ1:=k1=sqrt(2*m*E)/h; EQ2:=k2=sqrt(2*m*(E-U0))/h; E:=solve(EQ1,E); U0:=solve(EQ2,U0); U(x):=U0/(1+exp(-alpha*x));

 (1.3)

 (1.3)

 (1.3)

 (1.3)

 (1.3)

and obtain

 > SchrodingerEQN:=SchrodingerEQN;

 (1.4)

2. Solution of Schrodinger equation

Introduce new variable as

 > z(x)=-exp(-alpha*x); x=solve(z=-exp(-alpha*x),x);

 (2.1)

 (2.1)

Differentiation by new variable leads to

 > Diff(psi(x),x)=Diff(z(x),x)*Diff(psi(z),z); Diff(psi(x),x)=diff(-exp(-alpha*x),x)*diff(psi(z),z); Diff(psi(x),x) =-alpha*z*diff(psi(z),z); Diff(psi(x),x,x)=simplify(-alpha*z*diff(-alpha*z*diff(psi(z),z),z));

 (2.2)

 (2.2)

 (2.2)

 (2.2)

In this variable Schrodinger equation is

 > SchrodingerEQN:=alpha^2*z*(diff(psi(z),z)+z*diff(psi(z),`\$`(z,2)))+2*m/h^2*(1/2*k1^2*h^2/m-1/2*h^2*(-k2^2+k1^2)/m/(1-z))*psi(z);

 (2.3)

Next step is using substitution

 > psi(z):=F(z)*z^(-I*k2/alpha);

 (2.4)

and we have

 > SchrodingerEQN:=alpha^2*z*(diff(psi(z),z)+z*diff(psi(z),`\$`(z,2)))+2*m/h^2*(1/2*k1^2*h^2/m-1/2*h^2*(-k2^2+k1^2)/m/(1-z))*psi(z);

 (2.5)

 > (dsolve(SchrodingerEQN,F(z)));

 (2.6)

Thus solution of Schrodinger equation is

 > F(x):=subs(z=-exp(-alpha*x),_C1*hypergeom([-I*(k1+k2)/alpha,(-k2+k1)/alpha*I],[(alpha-2*I*k2)/alpha],z)+_C2*z^(2*I*k2/alpha)*hypergeom([-I/alpha*(-k2+k1), (k1+k2)/alpha*I],[1/alpha*(alpha+2*I*k2)],z)): psi(x):=expand(subs(z=-exp(-alpha*x),F(x)*z^(-I*k2/alpha)));

 (2.7)

or

 > psi(x):=(-1)^(-I*k2/alpha)*(exp(x*I*k2))*_C1*hypergeom([-I*(k1+k2)/alpha, (-k2+k1)/alpha*I],[(alpha-2*I*k2)/alpha],-exp(-alpha*x))+(-1)^(I*k2/alpha)*simplify((-exp(x*k2*I))*_C2*(-exp(-x*2*I*k2)))*hypergeom([-I/alpha*(-k2+k1), (k1+k2)/alpha*I],[1/alpha*(alpha+2*I*k2)],-exp(-alpha*x));

 (2.8)

3. Analysis of the solution

Note, that two linealy independent solutions of Schrodinger equation are complex conjugate each other. Now we choose solution corresponding to transmitted wave with dependence of x:

.

So we have:

 > psi[trans](x):=exp(k2*x*I)*hypergeom([-I*(k1+k2)/alpha,(-k2+k1)/alpha*I],[(alpha-2*I*k2)/alpha],-exp(-alpha*x));

 (3.1)

or

 > psi[trans](x):=exp(k2*x*I)*hypergeom([a,b],[c],-exp(-alpha*x));

 (3.2)

where used realtions:

 > a=-I*(k1+k2)/alpha; b=(-k2+k1)/alpha*I; c=expand((alpha-2*I*k2)/alpha);

 (3.3)

 (3.3)

 (3.3)

Obtained solution must be extend to x ---> -infinity. But Gauss hypergeometric function

is defined on interval . So we need to use Gauss relation [2]

or

 > C1=(GAMMA(c)*GAMMA(b-a))/((GAMMA(b)*GAMMA(c-a))); C2=(GAMMA(c)*GAMMA(a-b))/((GAMMA(a)*GAMMA(c-b))); hypergeom([a, b],[c],z) = C1*[(-z)^(-a)]*hypergeom([a, 1-c+a],[1-b+a],1/z)+C2*[(-z)^(-b)]*hypergeom([b, 1-c+b],[1-a+b],1/z);

 (3.4)

 (3.4)

 (3.4)

Now solution will have form:

 > psi(x):=exp(k2*x*I)*(C1*[(-z)^(-a)]*hypergeom([a, 1-c+a],[1-b+a],1/z)+C2*[(-z)^(-b)]*hypergeom([b, 1-c+b],[1-a+b],1/z));

 (3.5)

or

 > psi(x):=subs(z=-exp(-alpha*x),exp(k2*x*I)*(C1*subs(a = -I*(k1+k2)/alpha,[(-z)^(-a)])*hypergeom([a, 1-c+a],[1-b+a],1/z)+C2*subs(b = (-k2+k1)/alpha*I,[(-z)^(-b)])*hypergeom([b, 1-c+b],[1-a+b],1/z)));

 (3.6)

and

 > psi(x):=exp(-k1*x*I)*C1*hypergeom([a,1-c+a],[1-b+a],-1/exp(-alpha*x))+exp(k1*x*I)*C2*hypergeom([b,1-c+b],[1-a+b],-1/exp(-alpha*x));

 (3.7)

4. Probability flux densities

As

 > limit(-1/exp(-alpha*x), x=-infinity) assuming alpha>0;

 (4.1)

we obtained asymptotics solution at x = - infinity:

 > psi[as](x):=exp(-I*k1*x)*C1*hypergeom([a, 1-c+a],[1-b+a],0)+exp(k1*x*I)*C2*hypergeom([b, 1-c+b],[1-a+b],0);

 (4.2)

Physical  interpretation of obtained solutions (and its asymptotics) are

1) incoming wave:

 > psi[incom](x):=C2*exp(k1*x*I)*hypergeom([b, 1-c+b],[1-a+b],-1/exp(-alpha*x)); psi[incom_as](x):=C2*exp(k1*x*I)*hypergeom([b, 1-c+b],[1-a+b],0);

 (4.3)

 (4.3)

2) reflected wave:

 > psi[reflec](x):=C1*exp(-I*k1*x)*hypergeom([a, 1-c+a],[1-b+a],-1/exp(-alpha*x)); psi[reflec_as](x):=C1*exp(-I*k1*x)*hypergeom([a, 1-c+a],[1-b+a],0);

 (4.4)

 (4.4)

3) transmitted wave:

 > psi[trans](x):=exp(k2*x*I)*hypergeom([a, b],[c],-exp(-alpha*x)); psi[trans_as](x):=exp(k2*x*I)*hypergeom([a, b],[c],0);

 (4.5)

 (4.5)

So probability flux densities for asymptotics are [1]

 > with(tensor): J[incom](x):=simplify(evalc((I*h/(2*m))*(psi[incom_as](x)*conj(diff(psi[incom_as](x),x))-conj(psi[incom_as](x))*diff(psi[incom_as](x),x)))); J[reflec](x):=simplify(evalc((I*h/(2*m))*(psi[reflec_as](x)*conj(diff(psi[reflec_as](x),x))-conj(psi[reflec_as](x))*diff(psi[reflec_as](x),x)))); J[trans](x):=simplify(evalc((I*h/(2*m))*(psi[trans_as](x)*conj(diff(psi[trans_as](x),x))-conj(psi[trans_as](x))*diff(psi[trans_as](x),x))));

 (4.6)

 (4.6)

 (4.6)

5. Reflection and transmitted coefficients

Reflection coefficient is calculated by the formula [1]:

 > Reflection:=abs(J[reflec](x)/J[incom](x)); Transmitted:=abs(J[trans](x)/J[incom](x));

 (5.1)

 (5.1)

and using relations

 > a:=-I*(k1+k2)/alpha; b:=(-k2+k1)/alpha*I; c:=1-2*I/alpha*k2; C1:=(GAMMA(c)*GAMMA(b-a))/((GAMMA(b)*GAMMA(c-a))); C2:=(GAMMA(c)*GAMMA(a-b))/((GAMMA(a)*GAMMA(c-b)));

 (5.2)

 (5.2)

 (5.2)

 (5.2)

 (5.2)

we obtain

 > Reflection:=simplify((C1/C2)*conj(C1/C2)); Transmitted:=k2/k1*simplify((1/C2)*conj(1/C2));

 (5.3)

 (5.3)

and ascertain

 > T+R=simplify(convert(Reflection+Transmitted,exp));

 (5.4)

Now we can make up inverse transformation for constants and obtain reflection coefficient:

 > restart; Reflection:=(cosh(Pi/alpha*(-k2+k1))^2-1)/(cosh(Pi/alpha*(k1+k2))^2-1): Transmitted := sinh(2*Pi/alpha*k2)*sinh(2*Pi/alpha*k1)/(cosh(Pi/alpha*(k1+k2))^2-1): REFLACTION(E,alpha):=simplify(subs([k1=sqrt(2*m*E)/h,k2=sqrt(2*m*(E-U0))/h],Reflection)); TRASMITTED(E,alpha):=simplify(subs([k1=sqrt(2*m*E)/h,k2=sqrt(2*m*(E-U0))/h],Transmitted)); T+R=simplify(convert(REFLACTION(E,alpha)+TRASMITTED(E,alpha),exp));

 (5.5)

 (5.5)

 (5.5)

6. Visualization of reflection and transmitted coefficients

Specify systems parameter and visualized functions:

 > restart: U0:=10; m:=1; h:=1; alpha:=1; U1:=U0/(1+exp(-alpha*x)): REFLACTION1:=(cosh(Pi*2^(1/2)*((m*(E-U0))^(1/2)-(m*E)^(1/2))/alpha/h)^2-1)/(cosh(Pi*2^(1/2)*((m*E)^(1/2)+(m*(E-U0))^(1/2))/alpha/h)^2-1): TRASMITTED1:=sinh(2*Pi/alpha*2^(1/2)*(m*(E-U0))^(1/2)/h)*sinh(2*Pi/alpha*2^(1/2)*(m*E)^(1/2)/h)/(cosh(Pi*2^(1/2)*((m*E)^(1/2)+(m*(E-U0))^(1/2))/alpha/h)^2-1): alpha:=2; U2:=U0/(1+exp(-alpha*x)): REFLACTION2:=(cosh(Pi*2^(1/2)*((m*(E-U0))^(1/2)-(m*E)^(1/2))/alpha/h)^2-1)/(cosh(Pi*2^(1/2)*((m*E)^(1/2)+(m*(E-U0))^(1/2))/alpha/h)^2-1): TRASMITTED2:=sinh(2*Pi/alpha*2^(1/2)*(m*(E-U0))^(1/2)/h)*sinh(2*Pi/alpha*2^(1/2)*(m*E)^(1/2)/h)/(cosh(Pi*2^(1/2)*((m*E)^(1/2)+(m*(E-U0))^(1/2))/alpha/h)^2-1): alpha:=4; U3:=U0/(1+exp(-alpha*x)): REFLACTION3:=(cosh(Pi*2^(1/2)*((m*(E-U0))^(1/2)-(m*E)^(1/2))/alpha/h)^2-1)/(cosh(Pi*2^(1/2)*((m*E)^(1/2)+(m*(E-U0))^(1/2))/alpha/h)^2-1): TRASMITTED3:=sinh(2*Pi/alpha*2^(1/2)*(m*(E-U0))^(1/2)/h)*sinh(2*Pi/alpha*2^(1/2)*(m*E)^(1/2)/h)/(cosh(Pi*2^(1/2)*((m*E)^(1/2)+(m*(E-U0))^(1/2))/alpha/h)^2-1): alpha:=8; U4:=U0/(1+exp(-alpha*x)): REFLACTION4:=(cosh(Pi*2^(1/2)*((m*(E-U0))^(1/2)-(m*E)^(1/2))/alpha/h)^2-1)/(cosh(Pi*2^(1/2)*((m*E)^(1/2)+(m*(E-U0))^(1/2))/alpha/h)^2-1): TRASMITTED4:=sinh(2*Pi/alpha*2^(1/2)*(m*(E-U0))^(1/2)/h)*sinh(2*Pi/alpha*2^(1/2)*(m*E)^(1/2)/h)/(cosh(Pi*2^(1/2)*((m*E)^(1/2)+(m*(E-U0))^(1/2))/alpha/h)^2-1):

 (6.1)

 (6.1)

 (6.1)

 (6.1)

 (6.1)

 (6.1)

 (6.1)

Visualization of potential barrier at some values of :

 > plot([U1,U2,U3,U4],x=-4..4,color=[red,blue,gold,blue],thickness=[2,2,2,2]);

Dependences of reflection and transmitted coefficients by energy E:

 > plot([REFLACTION1,REFLACTION2,REFLACTION3,REFLACTION4],E=U0..U0+3,color=[red,blue,gold,blue],thickness=[2,2,2,2]); plot([TRASMITTED1,TRASMITTED2,TRASMITTED3,TRASMITTED4],E=U0..U0+3,color=[red,blue,gold,blue],thickness=[2,2,2,2]);

Dependences of reflection and transmitted coefficients by energy E and (3d surface plots):

 > restart: U0:=10; m:=1; h:=1; REFLACTION:=(cosh(Pi*2^(1/2)*((m*(E-U0))^(1/2)-(m*E)^(1/2))/alpha/h)^2-1)/(cosh(Pi*2^(1/2)*((m*E)^(1/2)+(m*(E-U0))^(1/2))/alpha/h)^2-1); TRASMITTED:=sinh(2*Pi/alpha*2^(1/2)*(m*(E-U0))^(1/2)/h)*sinh(2*Pi/alpha*2^(1/2)*(m*E)^(1/2)/h)/(cosh(Pi*2^(1/2)*((m*E)^(1/2)+(m*(E-U0))^(1/2))/alpha/h)^2-1);

 (6.2)

 (6.2)

 (6.2)

 (6.2)

 (6.2)

 > plot3d(REFLACTION,E=U0..U0+3,alpha=1..50, orientation=[-55,75],axes=BOXED,lightmodel=light4,style=PATCHCONTOUR); plot3d(TRASMITTED,E=U0..U0+3,alpha=1..50, orientation=[-55,75],axes=BOXED,lightmodel=light4,style=PATCHCONTOUR);

Dependences of reflection and transmitted coefficients by energy E and (Contour plots):

 > with(plots): contourplot(REFLACTION,E=U0..U0+3,alpha=1..50, contours=32,filled=true,coloring=[white,blue],grid=[35,35]); contourplot(TRASMITTED,E=U0..U0+3,alpha=1..50, contours=32,filled=true,coloring=[white,blue],grid=[35,35]);

References

1. Landau L.D., Lifshits E.V. Theoretical physics. V. 3. Quantum mechanics. - Moscow. Fizmatlit. 2003.

2. Handbook of Mathematical Functions. Edited by Milton Abramowitz and Irene A. Stegun. Issued June 1964. Tenth Printing, December. 1972.

Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the author for permission if you wish to use this application in for-profit activities.