Transmission coefficient

Quantum Mechanics
                       K. Konishi, G. Paffuti

Study in the WKB approximation the problem of transmission and reflection coefficients, in particular consider the cases

(Subscript[V, 1])[x]= (V0/2)(1+Cos[2π (x/Δ)])Θ[((Δ/2)+x)]Θ[((Δ/2)-x)];

(Subscript[V, 2])[x]= V0 Exp[-  x^2/ (Δ^2)];

ClearAll["Global`*"];
Off[General::spell, General::spell1];

Introduction

Definitions

In this notebook we compute the transmission coefficient for a potential of the form  (Subscript[V, 0])(1+Cos[x]) , bounded to  the region of one peak, the result are compared with WKB approximation. The technique is easily generalized for generic potential barriers of the form shown below.

The asymptotic form of the wave functions for large |x| is

ψ = ({{Exp[i k x]+ ℛ Exp[-i k x], x→ - ∞}, {𝒯 Exp[i k x], x→ + ∞}})

The conditions describe an incident wave from the left of unit flux, a reflected wave and a transmitted wave.The phases of coefficients ℛ and 𝒯 depend on the definitions of the overall phase of ψ. The quantities

T = ((|𝒯|)^2) ;  R = ((|ℛ|)^2)

are respectively the trnsmission and reflection coefficients.

WKB results

In the WKB approximation is customary to choose the phases using the classical turning points as reference points,

(1/(Sqrt[p]))(Exp[i S[(Subscript[a, L]),x]] + ℛ Exp[-i S[(Subscript[a, L]),x]]);    x → -∞ 
𝒯(1/(Sqrt[p])) Exp[i S[(Subscript[a, R]),x]];    x → +∞

S[a,x] is the classical action from point a to point x in units of ℏ. Let us note that S[a,x] becomes large and negative as x → -∞, as + i k x in eq.(), this explains the identification of incident and reflected waves in eq.().

The general WKB result is (see notebooks [*])

𝒯 = (1/(Sqrt[1+(e^(2σ)) ]))(e^(-i δ))= ((e^(-σ))/(Sqrt[1+(e^(-2σ)) ]))(e^(-i δ));     ℛ = -i((e^σ)/(Sqrt[1+(e^(2σ)) ]))(e^(-i δ)) = -i(1/(Sqrt[1+(e^(-2σ)) ]))(e^(-i δ));

where (Subscript[x, 1])and (Subscript[x, 2]) are the classical inversion points and

σ = (1/( ℏ))(∫_(Subscript[x, 1])^(Subscript[x, 2]))|p[x]|𝕕x;            δ = (σ/π)Log[(σ/(π e))] + Arg[Γ((1/2)- i (σ/(π )))]

At lowest order only the result

T = ((|𝒯|)^2)≃ Exp[-2σ]

can be consistently deduced.

For scattering above the maximum (Subscript[V, M]) of the barrier the transmission coefficient goes rapidly to 1. For simplicity we can approximate the potential near the top with

V ∼ (Subscript[V, M]) - (1/2)(β^2)(x^2)

and we have. see notebook [],

T  = (1/(1 + (e^(-2π α)))) ∼ 1 - (e^(-2π α));       R =   ((e^(-2π α))/(1 + (e^(-2π α)))) ∼ (e^(-2π α)).

with

α=((|E - (Subscript[V, M])|)/ℏ)(Sqrt[(m/β)])

Numerical procedure

We will follow the following procedure. The ψ function is fixed with normalization 1 at large and positive x, i.e. its asymptotic form will be

left:  a exp(i k x) + b exp(-ikx);      right:   exp(ikx)

This will be taken as boundary condition for a sufficiently large and positive x = L.  With this choice of phases

𝒯 = (1/a);  ℛ = (b/a)

The coefficients 𝒯 and ℛ are computed from the ratio of transmitted and reflected waves with respect to incident wave. We assume for simplicity a symmetric potential, with (Subscript[a, L])= - (Subscript[a, R]). Comparing the values of ψ at x = ± l we get, see eq.():

((Exp[i k L])/(a Exp[-i k L]))=  ((𝒯 Exp[i S[(Subscript[a, R]),L]])/(Exp[i S[(Subscript[a, L]),-L]]))= ((𝒯 Exp[i S[(Subscript[a, R]),L]])/(Exp[-i S[(Subscript[a, R]),L]]))
((b Exp[i k(- L)])/(a Exp[-i k(- L)]))= ((ℛ Exp[-i S[(Subscript[a, L]),-L]])/(Exp[i S[(Subscript[a, L]),-L]]))= ℛ Exp[2 i S[(Subscript[a, R]),L]]

then

𝒯 = (1/a)Exp[-2i S[(Subscript[a, R]),L] + 2i kL]; ℛ = (b/a)Exp[-2 i S[(Subscript[a, R]),L]+ 2i kL].

The coefficient a and b are easily expressed in terms of the value of ψ[x] at large negative x. Let f[x] the solution of the Schrödinger equation, then it is easily verified that

((f')/(f ))= ik ((a-b)/(a+b)) ;      S = ((f')/(i k f)) ⟹       (b/a)= ((1-S)/(1+S))exp(2 i k x) ⟹ 
            a  = (e^(- ik x))(1/2)(f + ((f')/ik)) =(e^(- ik x))(f/2)(1 + S)

Which allows an immediate evaluation of a in terms of the solution.

Numerical example

We write a couple of potentials which can be used by the reader as examples

V0=50.;Δ = 2.;

Vexample1[x_]:= (V0/2)(1+Cos[2π (x/Δ)])UnitStep[((Δ/2)+x)((Δ/2)-x)];

Vexample2[x_]:= V0 Exp[-  x^2/ (Δ^2)];

Here there are their shapes, quite similar:

Block[{xmax = 2.5 Δ},
GraphicsRow[{ Plot[Vexample1[x],{x,-xmax,xmax},AxesLabel-> {"x","V[x]"},
PlotRange-> All,Ticks-> None],
 Plot[Vexample2[x],{x,-xmax,xmax},AxesLabel-> {"x","V[x]"},PlotRange-> All,Ticks-> None]
},ImageSize-> 72 6]]

In our example we will take the second potential.

V[x_]:= Vexample1[x]

We first choose a list of energies at which the computation must be done

dV1 = V0/20;  dV2 = V0/50; dV3 = V0/300;
eneList = Join[Range[10*dV1,V0/2-dV1,dV1],Range[V0/2,V0 -  dV2 -dV3,dV2],Range[V0-dV2,V0-dV3,dV3]];

then a limiting point to the right where the potential is negligible

L = 5 Δ;

The function below compute ℛ and 𝒯 using the above procedure. The boundary conditions are given at x = L and the function is recunstructed by a backward integration to x = -L.

getRT[ene_,L_]:= Module[{k,φLeft,φRight,f,df,Sf,t,r},
k = (Sqrt[2.ene]);
{φLeft,φRight} ={ Exp[-I k L],  Exp[I k L]};
solution=NDSolve[{y''[x] + k^2y[x] -2 V[x]y[x] ==0,
y[L]== φRight,y'[L]==I k φRight },y,{x,L,-L}];
f = y[-L]/.solution[[1]]; df =  y'[-L]/.solution[[1]];
Sf = (df/(I k f));  
{t,r} ={(2/(f(1+Sf)))φLeft, ((1-Sf)/((1+Sf)))((φLeft)^2)};
{t,r}]

The list numTR contains, for each energy, the coefficients 𝒯 and ℛ

numTR = getRT[#1,L]&/@eneList;
{𝒯num,ℛnum} = Transpose[numTR];

WKB computation

We have to compute (Subscript[a, R]) , the action S[(Subscript[a, R]),L] and the action S[(Subscript[a, L]),(Subscript[a, R])] (in the forbidden region)

fzero[x_?NumberQ,ene_,V_]:=  ene - V[x];

getParametersWKB[ene_,V_,L_]:= Module[{aR,S,σ},
aR = Re[x/.FindRoot[fzero[x,ene,V] == 0, {x,(Δ/4.),0.,5Δ}]];
S =Re[ NIntegrate[(Sqrt[2(ene - V[x])]),{x,aR,L},AccuracyGoal-> 7]];
σ =2* Re[ NIntegrate[(Sqrt[2(-ene + V[x])]),{x,0,aR},AccuracyGoal-> 7]];
{aR,S,σ}];

The list paramWKB contains, for each energy, {aR, S in allowed region, σ in forbidden region}.

paramWKB=getParametersWKB[#1,V,L]&/@eneList;
{aR,S,σ} = Transpose[paramWKB];

Now we compute |𝒯| and δ in the WKB approximation using equations ()

modTwkb= (Sqrt[((Exp[-  2#1])/(1+  Exp[- 2 #1]))])&/@σ;

δwkb= Mod[ ((#1)/π)Log[((Abs[#1])/(π E))] +
 Arg[Gamma[(1/2) - I ((#1)/π)]],-2Pi]&/@ σ;

The same phase can be extracted from numerical computation, using ()

k = (Sqrt[2 eneList]);

δ = - Arg[𝒯num Exp[-2I S +2 I k L]];

Here is a comparison of the results

NumberForm[TableForm[Transpose[{eneList/V0,Abs[𝒯num],modTwkb,δ,δwkb}],TableHeadings->{{}, {"E/V0","\[LeftBracketingBar]𝒯\[RightBracketingBar]","\[LeftBracketingBar](Subscript[𝒯, WKB])\[RightBracketingBar]","δ","(Subscript[δ, WKB])"}},
TableSpacing-> {1,3}],4,ExponentFunction->(If[-10<#≤ -3,#,Null]&)]

In graphical form, continuous lines are WKB predictions:

ListPlot[{Transpose[{eneList,Abs[𝒯num]}],Transpose[{eneList,modTwkb}]},
Joined->{False, True},AxesLabel-> {"E","\[LeftBracketingBar] 𝒯 \[RightBracketingBar]"}]

ListPlot[{Transpose[{eneList,δ}],Transpose[{eneList,δwkb}]},
Joined->{False, True},AxesLabel-> {"E","δ"}]