HW #10

1. Neutrino Oscillation

(a)  Time evolution operator

\(TraditionalForm\`e\^\(\(-i\)\ H\ t/\) = \(∑\_k\(\((\(-i\)\ t/)\)\^k\/\(k!\)\) 
              H\^k = ∑\_k\(\((\(-i\)\ t/)\)\^k\/\(k!\)\) \((U\ E\ U\^†)\)\^k\)\)

where E is the diagonal matrix of Hamiltonian eigenvalues (Subscript[E, n]) given in Eq. 1.  Expanding a bit further,

	\(TraditionalForm\`e\^\(\(-i\)\ H\ t/\) = \(\(∑\_k\(\((\(-i\)\ t/)\)\^k\/\(k!\)\) \((U\ E\ U\^†)\)·\((U\ E\ U\^†)\)· ... \)\(.\)\)·\((U\ E\ U\^†)\)·\((U\ E\ U\^†)\)\_k\).

We see that apart from the ends, every (U^†) has a U to the right of it.  A matrix is unitary if (U^†)=(U^(-1)) and ((((U^†)))^†) = U, so the expression simplifies to

	  \(TraditionalForm\`e\^\(\(-i\)\ H\ t/\) = \(\(U(∑\_k\(\((\(-i\)\ t/)\)\^k\/\(k!\)\) E\^n)\) U\^† = 
          U\ \(e\^\(\(-i\)\ E\ t/\)\) U\^†\)\)
	  
	  ⟹ \(e\^\(\(-i\)\ H\ t/\)\)=U (({{\(e\^\(\(-i\)\ \(E\_1\) t/\)\), 0, ⋯, 0}, {0, \(e\^\(\(-i\)\ \(E\_2\) t/\)\), ⋯, 0}, {⋮, ⋮, ⋱, ⋮}, {0, 0, ⋯, \(e\^\(\(-i\)\ \(E\_n\) 
                            t/\)\)}}))(U^†)
	  
by the properties of matrix multiplication.

(b)  Two-state transition probabilities

The probabilities are

	\(TraditionalForm\`P(
          1 → 
            2) = \(\(\(|\)\(⟨2 | 
                e\^\(\(-i\)\ H\ t/\) | 
                  1⟩\)\( | \^2\)\) = \(\(|\)\(⟨2 | 
                U\ \(e\^\(\(-i\)\ E\ t/\)\) U\^† | 
                  1⟩\)\( | \^2\)\)\)\)
	
	\(TraditionalForm\`P(
          2 → 
            1) = \(\(|\)\(⟨1 | 
              U\ \(e\^\(\(-i\)\ E\ t/\)\) U\^† | 
                2⟩\)\( | \^2\)\)\).
	
We will demonstrate that there is no T-violation by showing P(1 → 2) - P(2 → 1) = 
        0.   Let us have Mathematica do the work:

U2m=(E^(I θ))(({{Cos[θ] (E^(I ϕ)), -Sin[θ] 
                    (E^(I η))}, {Sin[θ] (E^(-I η)), Cos[θ] 
                    (E^(-I ϕ))}}));
E2m=(({{E1, 0}, {0, E2}}));
amp212=(({{0, 1}})).U2m.\(MatrixExp[\(-I\)\ E2m\ t/]\).Transpose[Conjugate[U2m]].(({{1}, {0}}));
p212 = 
        ComplexExpand[Conjugate[amp212]*amp212];
amp221=(({{1, 0}})).U2m.\(MatrixExp[\(-I\)\ E2m\ t/]\).Transpose[Conjugate[U2m]].(({{0}, {1}}));
p221 = 
        ComplexExpand[
          Conjugate[amp221]*amp221];
TrigExpand[
        p212 - p221][[1, 1]]

(c)  Three-state problem

Let us proceed in similar fashion as in the two-state problem above:

U3m=(({{1, 0, 0}, {0, Cos[θ23], Sin[θ23]}, {0, -Sin[θ23], Cos[θ23]}})).(({{Cos[θ13], 0, Sin[θ13] (E^(-I δ))}, {0, 1, 0}, {-Sin[θ13] (E^(I δ)), 0, Cos[θ13]}})).(({{Cos[θ12], Sin[θ12], 0}, {-Sin[θ12], Cos[θ12], 0}, {0, 0, 1}}));
E3m=(({{E1, 0, 0}, {0, E2, 0}, {0, 0, E3}}));
amp312=(({{0, 1, 0}})).U3m.\(MatrixExp[\(-I\)\ E3m\ t/]\).Transpose[Conjugate[U3m]].(({{1}, {0}, {0}}));
p312 = ComplexExpand[Conjugate[amp312]*amp312];
amp321=(({{1, 0, 0}})).U3m.\(MatrixExp[\(-I\)\ E3m\ t/]\).Transpose[Conjugate[U3m]].(({{0}, {1}, {0}}));
p321 = 
        ComplexExpand[Conjugate[amp321]*amp321];
Simplify[
        TrigExpand[p312 - p321]][[1, 1]]

Indeed, when δ ≠ 0, the two probabilities are different.

(d) [optional]  CPT conservation

Again as above, substituting \(TraditionalForm\`U -> \ \(U\^*\)\):

amp3c12=(({{0, 1, 0}})).Conjugate[U3m].\(MatrixExp[\(-I\)\ E3m\ t/]\).Transpose[U3m].(({{1}, {0}, {0}}));
p3c12 = 
        ComplexExpand[Conjugate[amp3c12]*amp3c12];
amp3c21=(({{1, 0, 0}})).Conjugate[U3m].\(MatrixExp[\(-I\)\ E3m\ t/]\).Transpose[U3m].(({{0}, {1}, {0}}));
p3c21 = 
        ComplexExpand[Conjugate[amp3c21]*amp3c21];

\(TraditionalForm\`P(1 -> 2) - P(\(1\&-\) -> \(2\&-\))\):

Simplify[TrigExpand[p312 - p3c12]][[1, 1]]

\(TraditionalForm\`P \((1 -> 2)\) - 
        P(\(2\&-\) -> \(1\&-\))\):

TrigExpand[p312 - p3c21][[1, 1]]

2. Periodic delta-function Potential

(a)  Matching conditions

Using the form of the wavefunction given in the problem,

	ψ(-ϵ) = A + B
ψ(+ϵ) = (e^(i k a))(
          A (e^(-i κ a)) + B (e^(i κ a)))
ψ' (-ϵ) = 
        i κ(A + B)
ψ' (+ϵ) = 
        i κ (e^(i k a))(
            A (e^(-i κ a)) - B (e^(i κ a)))

The wavefunction is continuous ψ(-ϵ) = ψ(+ϵ), but its derivative is discontinuous because of the delta-function potential, \(TraditionalForm\`ψ' \((\(+ϵ\))\) - ψ' \((\(-ϵ\))\) = \(\(2  m\ λ\)\/ℏ\^2\) \(ψ(0)\)\).  These give the conditions

	A + B = (e^(i k a))(
          A (e^(-i x a)) + B (e^(i x a)))
i κ (e^(i k a))(
              A (e^(-i κ a)) - B (e^(i k a))) - 
          i κ(A - B) = ((2  m λ)/(ℏ^2)) (A + 
              B)

which we now solve:

msol=\(Solve[{
A + 
                B == \(E\^\(I\ k\ a\)\) \((A\ E\^\(\(-I\)\ κ\ a\) + B\ E\^\(I\ κ\ a\))\), 
            I\ κ\ \(E\^\(I\ k\ a\)\) \((A\ E\^\(\(-I\)\ κ\ a\) - B\ E\^\(I\ κ\ a\))\) - 
                I\ κ \((A - 
                      B)\) == \(\(2\ m\ λ\)\/\^2\) \((A + 
                    B)\)
}, 
{B, 
            k}
]\);

Inserting the solution for k into the phase (e^(i k a)):

\(foo1 = 
      FullSimplify[E\^\(I\ k\ a\) /. msol, \ 
        Assumptions -> { > 0, 
            m > 0, λ > 0, \ κ > 0, a > 0}]\)

Making the suggested substitution \(TraditionalForm\`d = \^2/m\ λ\):

\(phase = 
      Expand[Simplify[foo1 /. λ -> \^2/\((m*d)\), 
          Assumptions -> { > 0, d > 0, \ κ > 0, 
              a > 0}]]\)

We can read off the expressions for the two roots:

	(e^(i k a)) = 
              cos κ a + (1/(κ d)) sin κ a±\(\@\(\(\(\(cos\)\(\ \)\)\^2\) κ\ a + \(1\/\((κ\ d)\)\^2\) \(sin\^2\) κ\ a - 
                1 + \(1\/\(κ\ d\)\) sin\ 2  κ\ a\)\).
		
Recognizing that sin 2 κ a = 
        2  sin κ a ·cos κ a, we may factor the radicand and pull out a (Sqrt[-1]) to give

	(e^(i k a)) = 
              cos κ a + (1/(κ d)) sin κ a±\(i \@\( 1 - \((cos\ κ\ a + \(1\/\(κ\ d\)\) sin\ κ\ a)\)\^2\)\)
	
as desired.

(b)  Zero potential limit

In the limit d -> ∞, which is nothing but a free particle without a potential, we have 

	(e^(i k a)) = cos κ a ± 
            i (Sqrt[ 1 - (cos^2) κ a]) = 
          (e^(±i κ a)), 
	
and so \(TraditionalForm\`κ = \(±\((k + \(2  π\ n\)\/a)\)\)\); equivalently, k is the momentum modulo ((2  π n)/a).  Therefore, κ and hence the energy grow continuously as a function of k.  This can be seen numerically with a d that is large enough:

\(\(kplot1 = 
        Plot[\(\(-I\)/a*Log[phase[\([1]\)]] /. {a -> 1, 
                d -> 100}\) /. κ -> π*x, {x, 0, 4}, 
          PlotRange -> {\(-Pi\), Pi}, 
          AxesLabel -> {\*""\!\(\(κ\\ a\)\/π\)"", "ka"}];\)\)

\(\(kplot2 = 
        Plot[\(\(-I\)/a*Log[phase[\([2]\)]] /. {a -> 1, 
                d -> 100}\) /. κ -> π*x, {x, 0, 4}, 
          PlotRange -> {\(-Pi\), Pi}, 
          AxesLabel -> {\*""\!\(\(κ\\ a\)\/π\)"", "ka"}];\)\)

Show[kplot1, kplot2];

As we can see, no band gaps -- every κ has a real k.

(c)  Weak potential

Looking at the equation 

	\(TraditionalForm\`e\^\(i\ k\ a\) = 
        cos\ κ\ a + \(\(1\/\(κ\ d\)\) 
              sin\ κ\ a ± 
            i \@\( 1 - \((cos\ κ\ a + \(1\/\(κ\ d\)\) sin\ κ\ a)\)\^2\)\)\), 

if the argument of the square root is negative, the l.h.s. becomes pure real and cannot satisfy the equation for a real k.  Therefore there is no solution when |cos κ a + (1/(κ d)) 
            sin κ a|>1.  When d is finite but large, the combination exceeds unity for κ a = n π + ϵ (ϵ > 0).  This can be seen by expanding it in terms of ϵ, 

	\(TraditionalForm\`cos\ \((n\ π + ϵ)\) = \(\((\(-1\))\)\^n\) \((1 - ϵ\^2\/2 + O(ϵ\^4))\)\), sin(
          n π + ϵ) = (((-1))^n) (ϵ + 
              O((ϵ^3))), 
	
	⟹ \(TraditionalForm\`cos\ κ\ a + \(1\/\(κ\ d\)\) 
            sin\ κ\ a = \(\((\(-1\))\)\^n\) \((1 + \(1\/\(κ\ d\)\) ϵ - ϵ\^2\/2 + O(ϵ\^3))\)\). 
 
The magnitude exceeds unity for 0 < ϵ < 2/κ d ≃ 
        2  a /n π d.  The gap must exist just above κ = n π/a for any n, while the gap becomes smaller for large n.  So there exists a band gap at κ = ±π/a.

(d)  Plots for weak and strong potentials

Let us plot for the given cases as we did in (b), first the weak d = 3  a:

\(\(kplotcw1 = 
        Plot[\(\(-I\)/a*Log[phase[\([\)\(1\)\(]\)]] /. {a -> 1, 
                d -> 3}\) /. κ -> π*x, {x, 0, 4}, 
          PlotRange -> {\(-Pi\), Pi}, 
          AxesLabel -> {\*""\!\(\(κ\\ a\)\/π\)"", "ka"}];\)\)

\(\(kplotcw2 = 
        Plot[\(\(-I\)/a*Log[phase[\([2]\)]] /. {a -> 1, 
                d -> 3}\) /. κ -> π*x, {x, 0, 4}, 
          PlotRange -> {\(-Pi\), Pi}, 
          AxesLabel -> {\*""\!\(\(κ\\ a\)\/π\)"", "ka"}];\)\)

Show[kplotcw1, kplotcw2];

We see a big gap at κ = 0, a smaller one at κ = π/a, a yet smaller one at κ = 2  π/a, and a gap you can barely see at κ = 3  π/a.  This is exactly what we predicted in part (c).

Now let's do the strong case d = a/3:

\(\(kplotcs1 = 
        Plot[\(\(-I\)/a*Log[phase[\([\)\(1\)\(]\)]] /. {a -> 1, 
                d -> 1/3}\) /. κ -> π*x, {x, 0, 4}, 
          PlotRange -> {\(-Pi\), Pi}, 
          AxesLabel -> {\*""\!\(\(κ\\ a\)\/π\)"", "ka"}];\)\)

\(\(kplotcs2 = 
        Plot[\(\(-I\)/a*Log[phase[\([2]\)]] /. {a -> 1, 
                d -> 1/3}\) /. κ -> π*x, {x, 0, 4}, 
          PlotRange -> {\(-Pi\), Pi}, 
          AxesLabel -> {\*""\!\(\(κ\\ a\)\/π\)"", "ka"}];\)\)

Show[kplotcs1, kplotcs2];

Obviously, there is significant distortion from the free case in part (b), with gaps at κ = n π/a much bigger than in the weak case above.

(e)  (Subscript[Z, 2]) symmetry in k

It is parity that changes the overall sign of k.  This can be seen from the explicit form of the wave function,

	ψ(x) = 
        A (e^(i κ x)) + B (e^(-i κ x)) 		for -a < x < 0
	
	ψ(x) = (e^(i k a))(
          A (e^(i κ(x - a))) + 
            B (e^(-i κ(x - a)))) 	for 0 < x < a .

The parity transformation gives

	ψ(x) = (e^(i k a))(
          A (e^(i κ(-x - a))) + 
            B (e^(-i κ(-x - a))))
=B (e^(i (k + κ) a)) 
            (e^(i κ x)) + 
          A (e^(i (k - κ) a)) 
            (e^(-i κ x))
=A' (e^(i κ x)) + 
          B' (e^(-i κ x))
	
and

	ψ(x) = 
        B (e^(i κ x)) + A (e^(-i κ x))
=(e^(-i k a))(
          B (e^(i (k + κ) a)) 
              (e^(i κ (x - a))) + 
            A (e^(i (k - κ) a)) 
              (e^(-i κ (x - a))))
=(e^(-i k a))(
          A' (e^(i κ(x - a))) + 
            B' (e^(-i κ(x - a))))	

respectively.  The two wave functions are related by the changes

	A -> A' = 
        B (e^(i (k + κ) a)), 
	
	B -> B' = 
        A (e^(i (k - κ) a)), 
	
	(e^(i k a)) -> (e^(-i k a)).