Scattering from a spherical potential well
Our interest is in assembling the scattering amplitude for a spherical potential barrier or well. We calculate the phase shifts as a function of energy using the matching condition, and assemble the scattering amplitude. We use a simple unit system (h-bar=m=1), and measure length as a multiple of the scale set by the potential well.
> restart; with(orthopoly);
> R:=1; U0:=9/10; # (U0>0 = barrier, U0<0 = well)
> k:=1;
> kappa:=sqrt(k^2-U0);
Now we need the two types of Bessel functions.
> jn:=(n,x)->simplify(sqrt(Pi/2/x)*BesselJ(n+1/2,x));
> assume(r>0);
> jn(2,r);
> nn:=(n,x)->simplify((-1)^(n+1)*sqrt(Pi/2/x)*BesselJ(-n-1/2,x));
> nn(2,r);
> ul:=(l,x)->x*jn(l,x);
> vl:=(l,x)->x*nn(l,x);
The regular and irregular solutions to the radial equations in terms of the Ricatti-Bessel and Ricatti-Neumann functions.
We need their derivatives:
> ulp:=unapply('simplify(diff(ul(l,x),x))',l,x);
> ulp(2,r);
> vlp:=unapply('simplify(diff(vl(l,x),x))',l,x);
> vlp(2,r);
Now we are ready for the matching condition which determines the phaseshift.
The derivative functions defined above do not work with composite arguments yet:
> vlp(1,k*r);
We need to specify a simple argument, and then substitute. Our hope was that unapply would take care of that, but, of course, with the l -value unspecified nothing can be done. Therefore, we define the required function for specified l inside the loop. The trick is to force the evaluation of ulp and vlp before unapplying the answer (this does not work with the simple mapping construct).
> Lmax:=10;
> for l from 0 to Lmax do:
> rho:=kappa*R; Ulp:=unapply(simplify(ulp(l,r)),r): Vlp:=unapply(simplify(vlp(l,r)),r):
> taneta[l]:=evalf((Ulp(k*R)-Ulp(rho)*ul(l,k*R)/ul(l,rho))/(-Vlp(k*R)+Ulp(rho)*vl(l,k*R)/ul(l,rho))); od:
> print(taneta);
For very small phase shift (which happens for large l at finite wavenumber k ) the number can be undefined due to the lack of numerical precision.
For scattering from repulsive potential barriers the above matching condition works only when the particle can penetrate the barrier as the solution is assumed to be of a simple scattering type there (regular scattering solution). For scattering energies below the barrier height the particle would be tunneling into the potential barrier, and an exponentially dying function would be required.
> sigma:=4*Pi/k^2*add((2*l+1)*sin(arctan(taneta[l]))^2,l=0..3);
The partial-wave contributions:
> 4*Pi/k^2*[seq((2*l+1)*sin(arctan(taneta[l]))^2,l=0..Lmax)];
> for Lm from 0 to Lmax do: dsdO[Lm]:=1/4/k^2*evalc(abs(add((2*l+1)*(exp(2*I*arctan(taneta[l]))-1)*P(l,cos(theta)),l=0..Lm))^2); od:
> plot([seq(dsdO[Lm],Lm=0..4)],theta=0..Pi,color=[red,blue,green,brown,black]);
Exercise 1:
Explore the convergence of the partial-wave expansion of the differential cross section for different values of the wavenumber k (variation of the scattering energy).
Exercise 2:
Compare the differential cross section at fixed energy for repulsive and attractive potential well scattering.
Exercise 3:
Explore the dependence of the differential cross section at fixed energy on the size of the potential well (or barrier). Comment in particular on the amount of backscattering. Observe what happens in scattering from barriers of growing size when the scattering energy is just sufficient for barrier penetration without tunneling. What is the meaning of the structures in the differential cross section?
Energy dependence
Let us calculate the total scattering cross section as a function of energy in order to demonstrate the convergence of the partial wave expansion at low energies. We need to automate the calculation of the matching condition.
> with(plots):
Warning, the name changecoords has been redefined
> R:=2; U0:=-4; # (U0>0 = barrier, U0<0 = well)
> Lmax:=10;
> for ik from 1 to 10 do: k:=ik/10; kappa:=sqrt(k^2-U0); Ev[ik]:=k^2/2;
> for l from 0 to Lmax do:
> rho:=kappa*R; Ulp:=unapply(simplify(ulp(l,r)),r): Vlp:=unapply(simplify(vlp(l,r)),r):
> taneta[ik,l]:=evalf((Ulp(k*R)-Ulp(rho)*ul(l,k*R)/ul(l,rho))/(-Vlp(k*R)+Ulp(rho)*vl(l,k*R)/ul(l,rho)));
> sigma[ik,l]:=4*Pi/k^2*(2*l+1)*sin(arctan(taneta[ik,l]))^2;
> od:
> od:
> Lmax:=0;
> P1:=plot([seq([Ev[ik],log10(add(sigma[ik,l],l=0..Lmax))],ik=1..10)],style=point,color=red): display(P1);
> Lmax:=1;
> P2:=plot([seq([Ev[ik],log10(add(sigma[ik,l],l=0..Lmax))],ik=1..10)],style=point,color=blue): display(P1,P2);
> Lmax:=2;
> P3:=plot([seq([Ev[ik],log10(add(sigma[ik,l],l=0..Lmax))],ik=1..10)],style=point,color=green): display(P1,P2,P3);
> Lmax:=3;
> P4:=plot([seq([Ev[ik],log10(add(sigma[ik,l],l=0..Lmax))],ik=1..10)],style=point,color=black): display(P1,P2,P3,P4);
We see that at low scattering energies the cross section is dominated by isotropic (s-wave) scattering. As the energy increases the higher partial waves begin to contribute more to the total scattering cross section.
Exercise 4:
Explore the energy dependence of the total scattering cross section for attractive potential barriers of increasing depth and size. What is the interpretation of the minimum in the total cross section at low energies (which appears for potential wells of sufficient depth/size)?
>