Solution to 6.2

``Computational Physics'', J. M. Thijssen

(a)

Let us call the the region for (Subscript[A, n]),  (Subscript[B, n]) [ (n - 1) a + Δ/2  <  x <  
        n a  - Δ/2 ] region 1, (Subscript[C, n]), (Subscript[D, n]) [ n a - Δ/2  < x < 
        n a  + Δ/2] region 2, and (Subscript[A, n + 1]),  (Subscript[B, n + 1]) [ n a + Δ/2 < 
        x < (n  + 1) a - Δ/2] region 3.  We denote the wavefunctions in these three regions (Subscript[Ψ, 1]), (Subscript[Ψ, 2]), and (Subscript[Ψ, 3]). We have
(Subscript[Ψ, 1]) = 
        (Subscript[A, n]) (e^(i q(x - n a))) + (Subscript[B, n]) 
            (e^(-i q(x - n a))),
(Subscript[Ψ, 2]) = 
        (Subscript[D, n]) (e^(i k(x - n a))) + (Subscript[D, n]) 
            (e^(-i k(x - n a))),
(Subscript[Ψ, 3]) = 
        (Subscript[A, n + 1]) (e^(i q(x - (n + 1) a))) + (Subscript[B, n + 1]) 
            (e^(-i q(x - (n + 1) a))),
where q = (Sqrt[2  m E]), k = (Sqrt[2 m(E - (Subscript[V, 0]))]). The equations by matching at boundaries at n a  -  Δ/2 and n a + Δ/2 for continuous function and derivative continous are:
(Subscript[Ψ, 1 ])(
          n a - Δ/2)  =  (Subscript[Ψ, 2])(
          n a - Δ/2),  (Subscript[Ψ, 2 ])(
          n a + Δ/2)  =  (Subscript[Ψ, 3])(
          n a + Δ/2),
\(TraditionalForm\`\(\(\(\(Ψ'\)\_\(\(1\)\(\ \)\)\) \((n\ a - Δ/
                2)\)\  = \ \(\(Ψ'\)\_2\) \((n\ a - Δ/2)\)\)\(,\)\)\) \(TraditionalForm\`\(\(Ψ'\)\_\(\(2\)\(\ \)\)\) \((n\ a + Δ/
                2)\)\  = \ \(\(Ψ'\)\_3\) \(\((n\ a + Δ/2)\)\(.\)\)\)   
This gives us the following 4 equations.  Solving for (Subscript[A, n + 1]), (Subscript[B, n + 1]), eliminating (Subscript[C, n]) and (Subscript[D, n]), we get  (A and B is (Subscript[A, n]) and (Subscript[B, n]), A1 and B1 is (Subscript[A, n + 1]), (Subscript[B, n + 1])):

eq1 = 
        A*Exp[-I*q*Δ/2] + 
            B*Exp[I*q*Δ/2] == 
          Cn*Exp[-I*k*Δ/2]  + 
            Dn*Exp[I*k*Δ/2];
eq2 = 
        q*A*Exp[-I*q*Δ/2] - 
            q*B*Exp[I*q*Δ/2] == 
          k*Cn*Exp[-I*k*Δ/2]  - 
            k*Dn*Exp[I*k*Δ/2];
eq3 = 
        Cn*Exp[I*k*Δ/2] + 
            Dn*Exp[-I*k*Δ/2] == 
          A1*Exp[I*q*(Δ/2 - a)]  + 
            B1*Exp[-I*
                  q*(Δ/2 - 
                      a)];
eq4 = 
        k*Cn*Exp[I*k*Δ/2] - 
            k*Dn*Exp[-I*k*Δ/2] == 
          q*A1*Exp[I*q*(Δ/2 - a)]  - 
            q*B1*Exp[-I*q*(Δ/2 - a)];

res  = 
      Solve[{eq1, eq2, eq3, eq4}, {A1, B1}, {Cn, Dn}]  //  
        Simplify

Now (Subscript[T, ij]) for part (a) is computed as following:

T11  =  
      Simplify[(4*k/q)*A1  /.  res[[1]]  /.  {A -> 1,  
            B ->  0}]

T22  = 
      Simplify[ (4*k/q)*B1  /.  
            res[[1]]  /.  {A ->  0,  
            B ->  1}]

T12  =  
      Simplify[(4*k/q)*A1  /.  res[[1]]  /.  {A -> 0,  
            B ->  1}]

T21  =  
      Simplify[(4*k/q)*B1  /.  res[[1]]  /.  {A -> 1,  
            B ->  0}]

Check (Subscript[T, ij]) against with textbook result on page 142.  (checked)

Simplify[
      T11  -  Exp[
            I*q*(a - Δ)]*(Exp[
                  I*k*Δ]*(1 + k/q)^2  -  
              Exp[-I*k*Δ]*(1 - k/q)^2)]

Simplify[
      T12  -  (-2*I*Exp[I*q*a] (1 - k^2/q^2)*
            Sin[k*Δ])]

(b)

The determinant and trace of the matrix T has a simple result.

T = (q/(4  k))*{{T11, T12}, {T21, T22}};

Det[T]  //  Simplify

tr  =  FullSimplify[T[[1, 1]] + T[[2, 2]]]

So the eigenvalues of T have the property:
(Subscript[λ, 1]) (Subscript[λ, 2]) = 1  =  
          Det(T),  (Subscript[λ, 1]) + (Subscript[λ, 2])= Tr(T)= tr.
We note that (Subscript[Ψ, 1])and (Subscript[Ψ, 3]) are the wavefuctions at different period of the potential, thus by Bloch theorem, they most be related as (Subscript[Ψ, 1])(
          x - a) = (e^(-i p a)) (Subscript[Ψ, 3])(x), in term of the coefficients A and B,this means that (Subscript[A, n + 1]) =  (e^(i p a)) (Subscript[A, n]) and (Subscript[B, n + 1]) = (e^(i p a)) (Subscript[B, n]).  Since
T (({{(Subscript[A, n])}, {(Subscript[B, n])}}))=(({{(Subscript[A, n + 1])}, {(Subscript[B, n + 1])}})), we conclude that the eigenvalues of T must be (e^(i p a)), where p is the lattice momentum in Bloch theorem.  For each energy E the two eigenvalues correspond to two momenta with equal magnitude and opposite in sign.  Let (Subscript[λ, 1]) = (e^(i p a)), (Subscript[λ, 2]) = (e^(-i p a)), (Subscript[λ, 1])*(Subscript[λ, 2]) = 1 is automatically satisfied, we have (Subscript[λ, 1]) + (Subscript[λ, 2]) =  
        2 cos(p a) = tr.  A plot of the spetrum for p a vs E is given here. [Choosing the following constants]

nvalues = {V0  -> 1.0, a -> 1.0, Δ -> 0.2,  
          m -> 1.0};
spectrap[e_]  =  
      Simplify[ArcCos[tr/2]  /.  {q  ->  Sqrt[2*m*e],  
              k ->  Sqrt[2*m*(e - V0)]}  /.  
          nvalues]

ParametricPlot[{{spectrap[e], e}, {-spectrap[e], e}}, {e, 0.196,  
        30.0}]

Note that the plot is in the first Brillouin zone from -π to π.  The formula for p(E) appears also valid for  E < (Subscript[V, 0]), but I don't have a prove.

(c)

In APW approach, we take a plane wave in the region V(x) = 0, and in the barrier region, we solve the Schroedinger equation to match the boundary condition.  This is done in one unit cell from -a/2 to -a/2 (corresponding to n = 0).    Let χ(x)  =  (e^(i q x)) for |x|>Δ/
            2 and χ(x)  =  (Subscript[Ψ, 2])(x) for |x|<Δ/
            2, Bloch theorem requires that 
χ(-a/2)  =  χ(
            a/2) (e^(-i p a)),  where p is lattice momentum.
This condition requires that  q  -  p  =  2  π m/a, m is integer {...,-2,-1,0,1,2}.   Each of the q values will thus indexed by m.  [m has nothing to do with n]   The set of basis functions are (Subscript[χ, m])(x)  =  (e^(i (Subscript[q, m]) x))for |x|>Δ/
            2.  For |x|<Δ/
            2, we determine the coefficient C and D by continuity:
(e^(-i (Subscript[q, m]) Δ/2))  =  
        C (e^(-i k Δ/2)) + 
          D (e^(i k Δ/2)),
(e^(i (Subscript[q, m]) Δ/2))  =  
        C (e^(i k Δ/2)) + 
          D (e^(-i k Δ/2)).   
Sove for C and D we find:

res  =  
      Solve[{Exp[-I*qm*Δ/2]  ==  
            Cm*Exp[-I*k*Δ/2]  +  
              Dm*Exp[I*k*Δ/2], 
Exp[
              I*qm*Δ/2]  ==  
            Cm*Exp[I*k*Δ/2]  +  
              Dm*Exp[-I*k*Δ/2]}, {Cm, Dm}]

c  =  Cm  /.  res[[1, 1]]

d = Dm  /.  res[[1, 2]]

{Cm, Dm} = Simplify[ExpToTrig[{c, d}]]

This is in agreement with the book on page 143.   Now we need to compute the matrix elements S and H.  \(TraditionalForm\`S\_\(m\ l\) = ∫\_\(\(-a\)/2\)\%\(a/2\)\(\(χ\_m\^*\)(x)\) \(\(χ\_l\)(x)\) dx\).  We divide the integral in two regions, the external one for |x|>Δ/
            2 and internal one for \(TraditionalForm\`\(\(|\)\(x\)\(|\)\(\ \)\(\(<\)\(\ \)\(Δ/2. \)\)\)\)  Similarly for (Subscript[H, m l]).  We have 
\(TraditionalForm\`\(S\^ext\)\_\(m\ l\) = ∫\_\(\(\(a/2\)\(>\)\) | x | \(\(>\)\(Δ/2\)\)\)\(e\^\(\(-iq\_m\) 
                  x\)\) \(e\^\(i\ \(q\_l\) x\)\) dx\).   qsub[m] stands for (Subscript[q, m.]) \(TraditionalForm\`\(S\^int\)\_\(m\ l\) = ∫\_\(\(|\)\(x\)\(|\)\(\(<\)\(Δ/2\)\)\)\((\(C\_m\) 
                  e\^\(\(-i\)\ k\ x\) + \(D\_m\) 
                  e\^\(i\ k\ x\))\) \((\(C\_l\) e\^\(i\ k\ x\) + \(D\_l\) 
                  e\^\(\(-i\)\ k\ x\))\) dx\).
Note that there will be no m index on k.
We'll define (Subscript[C, m])=Csub[m] and (Subscript[D, m])=Dsub[m] and (Subscript[q, m])=qsub[m] later.
p is lattice momentum, but k = (Sqrt[2  m(E - (Subscript[V, 0]))]).  They are different.

outres  =  
      Simplify[ExpToTrig[
          Integrate[
              E^(I*(qsub[l]  -  qsub[m])*x), {x, -a/
                  2, -Δ/2}]  +  
            Integrate[
              E^(I*(qsub[l]  -  qsub[m])*x), {x, Δ/2, 
                a/2}] ]]

So[m_, 
        l_]  := ((a (Sin[(l - m) π] - Sin[(((l - m) π Δ)/a)]))/((l - m) π))  /;  
        m ≠ l
So[m_, m_]  :=  (a - Δ)

inres  =  
      Simplify[ExpToTrig[
          Integrate[(Csub[m]*E^(-I*k*x) + 
                  Dsub[m]*E^(I*k*x))*
(Csub[l]*
                    E^(I*k*x) + 
                  Dsub[l]*E^(-I*k*x)  ), {x, -Δ/
                2, Δ/2}]]]

Si[m_, l_]  =  inres

Clearly, the expressions for (S^ext) and (S^int) are wrong in part (c) on page 144 (even before substitute the values of (Subscript[C, m]), (Subscript[D, m]), and (Subscript[q, m])).

(d)

Now we compute \(TraditionalForm\`H\_ml\  = \ ∫\_\(\(-a\)/2\)\%\(a/2\)\(\(χ\_m\^*\)(
              x)\)\ \((\ \(-\ \((1/2)\)\)\ \ d\^2/dx\^2\  + \ 
                V(x))\)\ \(\(χ\_l\)(x)\)\ dx\).     We divide the integral into three terms, the internal contribution |x|<Δ/
            2, external contribution |x|>Δ/
            2, and a term due to second derivative to a discontinuity delta-function contribution.  
((Subscript[H, ml])^o)=\(TraditionalForm\`∫\_\(\(\(a/2\)\(>\)\) | x | \(\(>\)\(Δ/2\)\)\)\ \(\(χ\_m\^*\)(
              x)\) \((\((1/2)\) q\_l\^2)\)\ \(\(χ\_l\)(
              x)\)\ dx\  = \ \((1/2)\) 
          q\_l\^2\ \((S\^o)\)\_ml\).   [If we do an integration by part, we can make it symmetric in m and l.

Ho[m_, l_] :=  (1/2) qsub[l]^2 * So[m, l];

\(TraditionalForm\`H\_ml\^int\  = \(∫\_\(\(|\)\(x\)\(|\)\(\(<\)\(Δ/2\)\)\)\ \(\(χ\_m\^*\)(
                x)\) \((\((1/2)\) k\^2 + V\_0)\)\ \(\(χ\_l\)(
                x)\)\ dx = \ \((\((1/2)\) k\^2 + V\_0)\)*
            S\_ml\^int\)\), again the result is different from the textbook.

Hi[m_, l_] :=  ((1/2)*k^2  +  V0)*Si[m, l];

We find expression for ((Subscript[H, ml])^int):

Hi[m, l]

For the discontinuity, we need to quote the following general results:
define θ(x)  =  0 for x< 0, and 1 for x > 0.   Then θ(
          x)  =  (∫_(-∞)^x) δ(
              y) dy, thus θ' (x)  =  δ(x).    We need to compute the jump of the first derivative of (Subscript[χ, l])(x) at x = -Δ/2 and Δ/2.  If the discontinuous jump in d(Subscript[χ, l])(x)/dx is dΨ, then (d^2) (Subscript[χ, l])(x)/d (x^2) has dΨ δ(x - (Subscript[x, 0])), its integral is just dΨ.  Since
(Subscript[χ, 1])(x) = (e^(i (Subscript[q, l])x)), |x| > Δ/2, and \(TraditionalForm\`\(χ\_2\)(
          x)\  = \ \(C\_\(\(l\)\(\ \)\)\) e\^\(i\ k\ x\) + \(D\_l\) 
            e\^\(\(-i\)\ k\ x\), \ \(\(|\)\(x\)\(|\)\(\ \)\(\(<\)\(\ \)\(Δ/2\)\)\)\).  From this we find the discontinuity at -Δ/2 to be
dΨ(-Δ/
            2)  =  (Subscript[χ, 2])(-Δ/
                2)  -  (Subscript[χ, 1])(-Δ/2)  =  
          i k((Subscript[C, l]) (e^(-i k Δ/2)) - (Subscript[D, l]) 
                    (e^(i k Δ/2)))  -  
            i (Subscript[q, l]) (e^(-i (Subscript[q, l]) Δ/2)),
dΨ(Δ/
            2)  =  (Subscript[χ, 1])(Δ/
                2)  -  (Subscript[χ, 2])(Δ/2)  =  
          i (Subscript[q, l]) (e^(i (Subscript[q, l]) Δ/2)) -   i k((Subscript[C, l]) 
                    (e^(i k Δ/2)) - (Subscript[D, l]) 
                    (e^(-i k Δ/2))).
Don't forget we also have a continuous piece from \(TraditionalForm\`\(χ\_m\^*\)(x)\), the total contribution to ∂(Subscript[H, ml]) is \(TraditionalForm\`\(-\((1  \^2/2  
              m)\)\) \((\(\(χ\_m\^*\)(\(-Δ\)/
                  2)\)*\ \(dΨ(\(-Δ\)/
                  2)\)\  + \(\(χ\_m\^*\)(Δ/
                  2)\)*\ \(dΨ(Δ/2)\)\ )\)\):

pH[m_, 
        l_] :=  -(1/
              2) (Exp[
                I*qsub[m]*Δ/2]*(I*
                    k*(Csub[l]*Exp[-I*k*Δ/2]  -  
                        Dsub[l]*Exp[I*k*Δ/2])  -  
                  I*qsub[l]*
                    Exp[-I*
                        qsub[l]*Δ/
                          2])  +  
Exp[-I*
                  qsub[m]*Δ/2]*(-I*
                    k*(Csub[l]*Exp[I*k*Δ/2]  -  
                        Dsub[l]*Exp[-I*k*Δ/2])  +  
                  I*qsub[l]*
                    Exp[I*qsub[l]*Δ/2]) )

There is an expression for ∂(Subscript[H, m, l]):

ExpToTrig[pH[m, l]]  //  Simplify

The total Hamiltonian and total S matrix are

Ham[m_, l_] :=  
        Ho[m, l]  +  Hi[m, l]  +  pH[m, l];
Sam[m_, l_]  :=  So[m, l]  + Si[m, l]

Before any numerical computation, we need to define these:

Csub[m_]  :=  
        Csc[k Δ] Sin[
            (1/2) (k + 
                  qsub[m]) Δ];
Dsub[m_] := 
        Csc[k Δ] Sin[
            (1/2) (k - 
                  qsub[m]) Δ];
qsub[m_]  :=  p  +  2*Pi*m/a;

Now we have to solve the "eigenvalue" problem. This is not really an eigenvalue problem, since Hand S depends on E through k = (Sqrt[2 (E - (Subscript[V, 0]))]). Following the textbook suggestion, we fix a value p, try to find E such that det(H - E S)  =  0.   For a test, we use a 5 by 5 matrix, at (Subscript[V, 0]) = a = 1, Δ=0.2, and k = 2.0.

HmES  =  
        Table[ (Ham[m, l]  -  e*Sam[m, l])  /. {k -> 
                  Sqrt[2.0*(e - V0)], 
                p -> 2.0}  /.  {V0  -> 1.0, 
              a -> 1.0, Δ -> 0.2},  {m, -2, 
            2}, {l, -2, 2}];

This is a nonlinear function, too complicated to show:

poly = Chop[Det[HmES]];

Plot[{poly, 0}, {e, 0, 14},  
      PlotStyle -> {RGBColor[0, 0, 1], RGBColor[1, 0, 0]}]

We see for p=2.0, the solution for E is about 2.3 and 9.5.  Looking at the first plot, we see the it is about right.  We can be more qualitative by generation many points for different values of p.  But I am quite sure the formulas are correct.   I'll stop here.