Transmission coefficient

Quantum Mechanics
K. Konishi, G. Paffuti

Study the reflection and trasmission coefficients in one dimension. In particular study the case of a quasi periodic potential.

Packages

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

Introduction

Let V(x) be a one dimensional potential and assume that V(x) → 0 at infinity. We want to study a scattering problem in which a plane wave is coming from the left and scatters. As a result we have a transmitted and a reflected wave, i.e. asymptotically the wave function ψ is described by

Exp[i k x] + (Subscript[a, R])Exp[-i k x] (Underscript[⟵, x→-∞])  ψ[x]  (Underscript[⟶, x→∞])   (Subscript[a, T]) Exp[i k x]

The normalization of the incoming wave has been arbitrarily set to 1. Reflection and transmission coefficients are defined as the probability to be reflected and transmitted respectively, computed as fractions between the transmitted and reflected fluxes:

ℛ = |(Subscript[a, R])((|)^2)  ;  𝒯 = |(Subscript[a, T])((|)^2);    ℛ + 𝒯 = 1.

Last equality follows from the unitarity (conservation of the total probability). Often in practical applications it is convenient to choose the normalization as

a Exp[i k x] + b Exp[-i k x] (Underscript[⟵, x→-∞])  ψ[x]  (Underscript[⟶, x→∞])    Exp[i k x]

Evidently

(Subscript[a, T])= (1/a);  (Subscript[a, R])= (b/a).

Eq.() for (Subscript[a, R]) is correct also for potentials having different limits for x → ±∞. In this case T can be directly obtained by unitarity, eq.().

Transmission matrix

In the text has been shown that coefficients (Subscript[a, R]) and (Subscript[a, T])  determine, using linearity, a general matrix Ω (transmission matrix) which connects left and right asymptotic behavior  of ψ:

Ω = (({{1/(Subscript[(SuperStar[a]), T]), -(Subscript[(SuperStar[a]), R])/(Subscript[(SuperStar[a]), T]) }, {-(Subscript[a, R])/(Subscript[a, T]), 1/(Subscript[a, T])}}));  (({{a'}, {b'}}))= Ω (({{a}, {b}}));
𝒯 = 1/((Abs[Ω[[1,1]]])^2);   ℛ = 1 - 𝒯;

where (a,b) are the coefficients of Exp[± i  k x] for x→-∞,  (a',b') play the same role for x→∞. The matrix Ω  is unimodular, i.e. det[Ω] = 1 as a consequence of unitarity. Its eigenvalues are in the relation (Subscript[λ, 1])(Subscript[λ, 2]) = 1.

The matrix formulation () is particularly useful for periodic potentials. Let Δ the period and Ω the transmission matrix for a single potential. The matrix for N potentials can be easily shown to be (see text)

(Ω^((N))) = (𝒰^(-N))(((𝒰 Ω))^N);  𝒰 =(({{Exp[i k Δ], 0 }, {0, Exp[-i k Δ]}})).

𝒰 is the translation operator needed to shift from one potential to the next one. For large values of N eq.() projects on the subspace of greatest eigenvalue of Ω, e.g. (Subscript[λ, 1]). It is easy to see that (Subscript[a, T])∼ (((1/(Subscript[λ, 1])))^N) so the transmission is killed. The transmission is possible only if 𝒰Ω has eigenvalues of modulus 1, i.e. Exp[i ϕ]. Expressing ϕ throw the parameters of the matrix the following bound is found

|((Cos[k Δ + δ])/(|(Subscript[a, T])|))| = |Cos[ϕ]| ≤ 1; with  (Subscript[a, T])= |(Subscript[a, T])|Exp[i δ].

or

1 - |Re[ ((Exp[i k Δ])/(Subscript[(SuperStar[a]), T]))]| ≥ 0

As a function of energy (or k) the condition () selects gaps in the possible transmission, modelling the presence of gaps in conduction bands of solids.

Parity

In the text is shown that for even potentials the following relation holds

(Subscript[(SuperStar[a]), R])(Subscript[a, T])+ (Subscript[(SuperStar[a]), T])(Subscript[a, R])= 0;  or    Arg[(Subscript[a, R])/(Subscript[a, T])] = ± π/2.

Numerical approach

One of the simplest possibility is to solve the Schrödinger equation imposing the asymptotic behavior given in eq.(). The solution for x→-∞ gives a and b. Given ψ[x] = a Exp[i k x] + b Exp[-i k x] and ψ'[x] = i k ( a Exp[i k x] - b Exp[-i k x]) we can compute easily

(b/a)= ((1- ξ)/(1+ξ))Exp[2i k x];    a = Exp[-i k x]((ψ[x])/2)(1+ξ);  where:  ξ = ((ψ'[x])/(i k ψ[x])).

From eq.() coefficients (Subscript[a, R]) and (Subscript[a, T]) follow.

Example 1: reflection above a barrier

Let us consider the potential

V0=1; Δ=1.; m=1;
V[x_]:= (V0/(1+(E^(-x/Δ))));

We choose natural units, ℏ = m=1 and for definiteness we have selected the parameters V0 and Δ as 1.

Plot[V[x],{x,-10 Δ,10 Δ},AxesLabel->{"x","V[x]"}]

For energy E > V0 we aspect a transmitted and a reflected wave (reflection above the barrier). The wave number of the incoming and transmitted waves are:

ℏ kL = (Sqrt[2 m E]);   ℏ kR = (Sqrt[2m(E - V0)]).

We have selected this example as a check for the method, the known analytical solution for the reflection coefficient is

ℛ =(((((Sinh[π Δ (kL - kR)])/(Sinh[π Δ (kL + kR)]))))^2) .

ene = Table[V0+0.01j,{j,1,30}];
kL = Sqrt[2 m ene]; kR = Sqrt[2m(ene-V0)];

To solve numerically the equation we choose two distant points, one on the left and one on the right, xL and xR. xR is used to put boundary conditions, at xL the wave function ψ and its derivative are computed, then eq.() and () are used to find (Subscript[a, R]) and (Subscript[a, T]).

xR = 15 Δ; xL = -xR;

ξ = Table[sol = NDSolve[{y''[x] + kL[[j]]^2y[x] -2m V[x]y[x] ==0,
y[xR]== (E^(I kR[[j]] xR)),y'[xR]==I kR[[j]] (E^(I kR[[j]] xR))},y,{x,xR,xL}];
((y'[xL])/(I kL[[j]] y[xL]))/.sol[[1]],{j,1,Length[ene]}];

Let us compare the numerical

ℛ = ((Abs[((1- ξ)/(1+ξ))])^2);

and the analytical reflection coefficient

ℛTh =(((((Sinh[π Δ(kL-kR)] )/(Sinh[π Δ(kL+kR)] ))))^2);

Show[ListPlot[Transpose[{ene,ℛ}],PlotStyle->PointSize[0.015`],AxesLabel->{"E/V0","ℛ"}],ListPlot[Transpose[{ene,ℛTh}],Joined->True,PlotStyle->RGBColor[1,0,0]]]

Example 2: Periodic potential

Let us now consider the following periodic potential

V0=1; Δ = 2; m=1;
V[x_,nP_]:= (V0/2)(1+Cos[2π (x/Δ)])UnitStep[((nP/2)Δ + x)((nP/2)Δ - x)]

Plot[V[x,11],{x,-11 Δ,11 Δ},PlotRange->All,AxesOrigin->{-20,0},Epilog->Text["V[x]",{18,0.5`}]]

With the  procedure used in the previous example let us define a function which compute aR and aT for a list of energies (ene) and for an arbitrary number of copies of the potential (numpeaks)

calcCoefficient[ene_,numPeaks_]:=Module[{k,ψ,ξ,xR,xL,aR,aT},
xR = numPeaks *Δ/2; xL = -xR;
k = (Sqrt[2 ene]);
 ris = Table[sol = NDSolve[{y''[x] + k[[j]]^2y[x] -2m V[x,numPeaks]y[x] ==0,
y[xR]== (E^(I k[[j]] xR)),y'[xR]==I k[[j]] (E^(I k[[j]] xR))},y,{x,xR,xL}];
{y[xL],((y'[xL])/(I k[[j]] y[xL]))}/.sol[[1]],{j,1,Length[ene]}];
ψ = ris[[All,1]]; ξ = ris[[All,2]];
aR = ((1-ξ)/(1+ξ))Exp[2I k xL]; aT = (2/(ψ(1+ξ)))Exp[I k xL];
{aR,aT}]

Here is the list of energies and corresponding coefficients

numPeaks=11;
ene =Join[ Range[0.4V0,1.1V0,0.005 V0],Range[1.105 V0,4 V0,0.008 V0]];
{aR,aT}=calcCoefficient[ene,numPeaks];

ListPlot[Transpose[{ene,((Abs[aT])^2)}],PlotRange->All,Joined->True,AxesLabel->{"E","T"}]

It is already clear, for numpeaks = 11, the presence of a gap between E=1.5 and E=2. Let us stress that this energy range is above the maximum of the potential (V0 = 1). Below we plot the first part of the graph to show the oscillating behavior of the transmission coefficient.

pl2=ListPlot[Transpose[{ene,((Abs[aT])^2)}],PlotRange->{{0.4`,1.1`},Automatic},AxesLabel->{"E","T"}]

Parity

The potential is even, so relation () must holds: (Subscript[(SuperStar[a]), R])(Subscript[a, T])+ (Subscript[(SuperStar[a]), T])(Subscript[a, R])= 0

Max[Abs[Conjugate[aR]*aT + Conjugate[aT]*aR]]

Transmission matrix

Let us compute the coefficients aR and aT for a single copy of the potential.

{aR1,aT1}=calcCoefficient[ene,1];

We can now compute the transmission matrix using eq.() and ():

calcOmega[ene_,{aR_,aT_},n_]:= Module[{k,aRT,𝒰,𝒰N},
aRT = Transpose[{aR,aT}];
k = Sqrt[2 ene];
Ω = Apply[{Conjugate[{1/#2,-#1/#2}],{-#1/#2,1/#2}}&,Transpose[{aR,aT}],{1}];
If[n==1,Ω,
𝒰= Apply[DiagonalMatrix[{#1,#2}]&,Transpose[{Exp[I  k Δ],Exp[- I  k Δ]}],{1}];
𝒰N = Apply[DiagonalMatrix[{#1,#2}]&,Transpose[{Exp[-I n k Δ],Exp[ I n k Δ]}],{1}];
 Table[𝒰N[[j]].MatrixPower[𝒰[[j]].Ω[[j]],n],{j,1,Length[ene]}]]]

We used the Mathematica function Apply. Used at level {1} of a list, as in our case, its behavior is
     Apply[ F, {{ (Subscript[x, 1]),(Subscript[y, 1])},{ (Subscript[x, 2]),(Subscript[y, 2])},...}, {1}] = {F[(Subscript[x, 1]),(Subscript[y, 1])], F[(Subscript[x, 2]),(Subscript[y, 2])],... }
In our case the function F is the (symbolic) form of the matrix, #1 and #2 are the entries which must be aR and aT. We used Transpose[{aR,aT}] in order to produce a list in which elements are the couple of coefficients for each energy. An equivalent procedure would be:
Ω=Table[{{1/Conjugate[aT1[[j]]],-Conjugate[aR1[[j]]/aT1[[j]]]},{-aR1[[j]]/aT1[[j]],1/aT1[[j]]}},{j,1,Length[ene]}];
The reader can check the equivalence of the results. The same procedure has been used for 𝒰 and 𝒰N.

We compute the transmission matrix for 11 potentials

ΩN = calcOmega[ene,{aR1,aT1},numPeaks];

and compare the extracted transmission coefficient (red line) with the previous computation finding perfect agreement

Show[pl2,ListPlot[Transpose[{ene,(1/((Abs[ΩN[[All,1,1]]])^2))}],Joined->True,PlotStyle->RGBColor[1,0,0]]]

Many barriers ( 30 )

To  sharpen the presence of the gap let us consider the case of 30 potentials. In this case the usefulness of Ω is obvious.

ΩN = calcOmega[ene,{aR1,aT1},30];

Eq.() describes the conditions to be satisfied to permit transmission: the quantity bound defined below must be positive:

bound = 1 - Abs[Re[ Exp[I  Δ(Sqrt[2 ene])]/Conjugate[aT1]]];

This is in perfect coincidence with the computed transmission coefficient (obtained from ΩN )

Show[ListPlot[Transpose[{ene,(1/((Abs[ΩN[[All,1,1]]])^2))}],Joined->True,PlotRange->{{0.5`,3},Automatic},AxesLabel->{"E","T"}],ListPlot[Transpose[{ene,bound}],Joined->True,PlotStyle->{Dashing[{0.02`}],Red}]]