HW #10

1. Variational Method

Plot of the potential

The potential in the problem is V = 50 ((((e^(-x)) - 1))^2).

V = 50 ((((E^(-x)) - 1))^2)

Plot[V, {x, -1, 1}]

It is basically quadratic around zero and is very steep.  One may guess that the harmonic oscillator is a pretty good approximation.

Initial guess

Series[V, {x, 0, 3}]

If we try to identify the first term with the harmonic oscillator potential (1/2) m ω (x^2), we find ω = 10 because m = 1. Then one may hope that it would give the ground-state energy (1/2) ℏ ω = 5.  However, this hope is not quite fulfilled.  Using the groundstate wave function of the harmonic oscillator,

ψ = ((((( m ω)/(π ℏ))))^(1/
              4)) (E^(-m ω (x^2)/(2 ℏ)))

ψ = ψ /. {m -> 1, ω -> 10, ℏ -> 
            1}

K = ((-(ℏ^2))/(2  m)) 
          Integrate[ψ D[ψ, {x, 
                  2}], {x, -∞, ∞}] /. {ℏ -> 1, m -> 1}

P = 
      Integrate[(ψ^2) 
          V, {x, -∞, ∞}]

Ebar = K + P

N[%]

The error is

((% - 39/8)/(39/8))

and is bigger than 5%, but is pretty good, accurate within 7.22%.
Therefore, we are motivated to try a Gaussian as a trial function.

Gaussian trial function

ψ = (1/((π^(1/4)) (Δ^(1/2)))) 
        (E^(-(x^2)/(2  (Δ^2))))

K = ((-(ℏ^2))/(2  m)) 
          Integrate[ψ D[ψ, {x, 
                  2}], {x, -∞, ∞}, 
            Assumptions -> Δ > 0] /. {ℏ -> 1, 
          m -> 1}

P = 
      Integrate[(ψ^2) V, {x, -∞, ∞}, 
        Assumptions -> Δ > 0]

Ebar = K + P

FindMinimum[Ebar, {Δ, (10^(-1/2))}]

((%[[1]] - 39/8)/(39/8))

It has improved, but not yet within 5%.

Linear times Gaussian

One way to improve it further is the following.  Note that the potential is not parity symmetric, and hence the groundstate wavefunction is not expected to be an even function.  Because the potential is lower on the right, we expect the wave function is skewed towards the right.  Therefore, we can try

ψ = (1 + k x) 
        (E^(-(x^2)/(2  (Δ^2))))

It is not normalized at this point, and we calculate its norm

norm2 = 
      Integrate[(ψ^2), {x, -∞, ∞}, 
        Assumptions -> Δ > 0]

K = ((-(ℏ^2))/(2  m)) 
          Integrate[ψ D[ψ, {x, 
                  2}], {x, -∞, ∞}, 
            Assumptions -> Δ > 0] /. {m -> 
            1, ℏ -> 1}

P = 
      Integrate[(ψ^2) V, {x, -∞, ∞}, 
        Assumptions -> Δ > 0]

Ebar = ((K + P)/norm2)

solution = 
      FindMinimum[
        Ebar, {Δ, (10^(-1/2))}, {k, 0}]

((solution[[1]] - 39/8)/(39/8))

This is correct at the 1% level!  The wave function is therefore

(ψ/(Sqrt[norm2])) /. solution[[2]]

Plot[(ψ/(Sqrt[norm2])) /. solution[[2]], {x, -1, 1}]

As expected, it is more or less a Gaussian, but skewed to the right.
Obviously, there are many ways to improve Gaussian.  I hope you found one successfully.

The analytic solution

This problem is actually a special case of the Morse potential

V = D ((((E^(-a u)) - 1))^2)

This potential is a form of the inter-atomic potential in diatomic molecules proposed by P.M. Morse, Phys. Rev. 34, 57 (1929).  The variable u is the distance between two atoms minus its equilibrium distance u = r - (Subscript[r, 0]).  The Schrödinger equation can be solved analytically.  Expanding it around the minimum,

Series[V, {u, 0, 4}]

Harmonic oscillator approximation gives ω = (Sqrt[((2  (a^2) D)/m)]).  The correct energy eigenvalues are known to be (Subscript[E, n]) = ℏ ω(
              n + (1/2)) - (((ℏ^2) (a^2))/(2  
                  m)) (((n + (1/2)))^2), where the second term is called the anharmonicanharmonic correction.  For large n, the second term will dominate and the energy appears to become negative.  Clearly, the bound state spectrum does not go forever, and ends at a certain value of n, quite different from the harmonic oscillator.  But this is expcted because the potential energy asymptotes to D for u -> +∞ and hence states for E > D must be unbound and hence have a continuous spectrum.
The ground-state wave function is known to have the form

ψ = (E^(-d (E^(-a u)))) 
        (E^(-b a u/2))

where b = 2  d - 1, d = ((Sqrt[2  D m])/(a ℏ)).  Let us see that it satisfies the Schrödinger equation.

Simplify[(((-((ℏ^2)/(2  m)) D[ψ, {u, 2}] + D ((((E^(-a u)) - 1))^2) ψ))/ψ) /. {b -> 
              2  d - 1} /. {d -> ((Sqrt[2  D m])/(a ℏ))}]

The first term is \(TraditionalForm\`\(1\/2\) ℏ\ ω = ℏ\ a \@\( D\/\(2  m\)\)\) and the zero-point energy with the harmonic oscillator approximation.  The second term -(((ℏ^2) (a^2))/(8  m)) is the anharmonic correction.
The normalization is

norm2 = 
      Integrate[(ψ^2), {u, -∞, ∞}, 
        Assumptions -> a > 0 && d > 0]

Actually, this wave function could have been guessed if one pays a careful attention to the asymptotic behaviors.  For u -> ∞, the potential asymptotes to a constant D, and hence the wave function must damp exponentially (e^(-κ u)) with κ = (Sqrt[2  m | E | ])/ℏ.  For u -> -∞, the potential rises extremely steeply as D (e^(-2 a u)).  It suggests that the energy eigenvalue becomes quickly irrelevant, and the behavior of the wave function must be given purely by the rising behavior of the potential.  By dropping the energy eigenvalue and looking at the Schrödinger equation,
((-(ℏ^2))/(2  
                  m)) (((d^2) ψ)/(d (u^2))) + 
          D (e^(-2 a u)) ψ = 0,
and change the variable to y = (e^(-a u)), we find
\(TraditionalForm\`\(\(-ℏ\^2\)\/\(2  
                  m\)\) \((\(\(d\^2\) ψ\)\/\(d\ y\^2\) + \(1\/y\) \(d\ ψ\)\/\(d\ y\))\) + D\ ψ = 0\).
The second term in the parentheses is neglible for y -> ∞.  Therefore the wave function has the behavior ψ ∝ 
        (e^(-(Sqrt[2 m D]) y/ℏ a)).  Combining the behavior on both ends, the wave function has precisely the exact form given above.
This is the lesson: the one-dimensional potential problem is so simple that there are many ways to study the behavior of the wave function.  On the other hand, the real-world problem involves many more degrees of freedom.  In many cases, the Hamiltonian itself must be guessed.

Back to the Morse potential.  The case of the homework problem corresponds to D = 50, a = 1, m = 1, ℏ = 1.  Therefore the groundstate energy is \(TraditionalForm\`E\_0 = \(\(1\/2\) ℏ\ a \@\( D\/\(2  m\)\) - \(\(ℏ\^2\) a\^2\)\/\(8  m\) = \(5 - 1\/8 = 39\/8\)\)\).

Simplify[(1/(Sqrt[norm2])) ψ /. {b -> 
                2  d - 1} /. {d -> ((Sqrt[2  D m])/(a ℏ))} /. {D -> 50, a -> 1, m -> 1, ℏ -> 1}]

Plot[((800000000 (E^(-10 (E^(-u)) - ((19 u)/2))))/(567 (Sqrt[2431]))), {u, -1, 1}, 
      PlotPoints -> 50]

The variational linear times Gaussian wave function was

Plot[
      1.2847346074584307` (E^(-4.557919801537102` (x^2))) (1 + 0.7642931913579979` x), {x, -1, 1}, 
      PlotPoints -> 50]

Show[%, %%]

Quite close.

2. Neutrino Oscillation

(a)

The unitarity matrix is

{{1, 0, 0}, {0, Cos[(Subscript[θ, 23])], 
          Sin[(Subscript[θ, 23])]}, {0, -Sin[(Subscript[θ, 23])], 
          Cos[(Subscript[θ, 23])]}} // MatrixForm

{{Cos[(Subscript[θ, 13])], 0, Sin[(Subscript[θ, 13])] (E^(-I δ))}, {0, 
          1, 0}, {-Sin[(Subscript[θ, 13])] (E^(I δ)), 0, 
          Cos[(Subscript[θ, 13])]}} // MatrixForm

{{Cos[(Subscript[θ, 12])], Sin[(Subscript[θ, 12])], 0}, {-Sin[(Subscript[θ, 12])], 
          Cos[(Subscript[θ, 12])], 0}, {0, 0, 1}} // MatrixForm

U = {{1, 0, 0}, {0, Cos[(Subscript[θ, 23])], 
            Sin[(Subscript[θ, 23])]}, {0, -Sin[(Subscript[θ, 23])], 
            Cos[(Subscript[θ, 23])]}} . {{Cos[(Subscript[θ, 13])], 0, 
            Sin[(Subscript[θ, 13])] (E^(-I δ))}, {0, 1, 
            0}, {-Sin[(Subscript[θ, 13])] (E^(I δ)), 0, 
            Cos[(Subscript[θ, 13])]}} . {{Cos[(Subscript[θ, 12])], Sin[(Subscript[θ, 12])], 
            0}, {-Sin[(Subscript[θ, 12])], Cos[(Subscript[θ, 12])], 0}, {0, 0, 
            1}}

Udagger = Transpose[U /. {δ -> -δ}]

Simplify[Udagger . U]

For the two-by-two case, the unitarity matrix is

U23 = U /. {(Subscript[θ, 12]) -> 0, (Subscript[θ, 13]) -> 0}

% // MatrixForm

The main point in the calculation is that the Hamiltonian can be written as
H = U(c p + (((m^2) (c^3))/(2  p))) 
          (U^†),
and hence the time evolution operator is
(e^(-i H t/ℏ)) = 
          U (e^(-i(c p + (m^2) (c^3)/2  p) t/ℏ)) 
            (U^†)
=(e^(-i c p t/ℏ))U(({{(e^(-i ((Subscript[m, 1])^2) (c^3) t/2  
                              p ℏ)), 0, 0}, {0, (e^(-i ((Subscript[m, 2])^2) (c^3) t/2  
                              p ℏ)), 0}, {0, 0, (e^(-i ((Subscript[m, 3])^2) (c^3) t/2  
                              p ℏ))}}))(U^†)
The amplitude of our interest is

A = (E^(-I c p t/ℏ)) ({0, 1, 0} . U23 . 
              DiagonalMatrix[{(E^(-I ((Subscript[m, 1])^2) (c^3) t/(2  p ℏ))), (E^(-I ((Subscript[m, 2])^2) (c^3) t/(2  p ℏ))), 
                  (E^(-I ((Subscript[m, 3])^2) (c^3) t/(2  p ℏ)))}] . Transpose[U23] . {{0}, {1}, {0}})[[1]]

Astar = A /. {t -> -t}

Psurv = Simplify[ExpToTrig[Expand[A Astar]]]

This can be rewritten as
(Subscript[P, surv])
\(TraditionalForm\`\(\(=\)\(\(cos\^4\) θ\_23 + 
          2 \( cos\^2\) \(θ\_23\) \(sin\^2\) \(θ\_23\) 
            cos \(\((m\_2\^2 - m\_3\^2)\) \(c\^3\) t\)\/\(2  ℏ\ p\) + \(sin\^4\) θ\_23\)\)\)
\(TraditionalForm\`\(\(=\)\(\((\(cos\^2\) θ\_23 + \(sin\^2\) θ\_23)\)\^2 - 2 \( cos\^2\) \(θ\_23\) \(sin\^2\) θ\_23 + 
          2 \( cos\^2\) \(θ\_23\) \(sin\^2\) \(θ\_23\) 
            cos \(\((m\_2\^2 - m\_3\^2)\) \(c\^3\) t\)\/\(2  ℏ\ p\)\)\)\)
\(TraditionalForm\`\(\(=\)\(1 - 
          2 \( cos\^2\) \(θ\_23\) \(sin\^2\) \(\(θ\_23\)(
              1 - cos \(\((m\_2\^2 - m\_3\^2)\) \(c\^3\) t\)\/\(2  ℏ\ p\))\)\)\)\)
\(TraditionalForm\`\(\(=\)\(1 - 
          4 \( cos\^2\) \(θ\_23\) \(sin\^2\) \(θ\_23\) \(sin\^2\) \(\((m\_2\^2 - m\_3\^2)\) \(c\^3\) t\)\/\(4  ℏ\ p\)\)\)\)
=1 - (sin^2) 
            2  (Subscript[θ, 23]) (sin^2) (((((Subscript[m, 2])^2) - 
                        ((Subscript[m, 3])^2)) (c^3) 
                    t)/(4  ℏ p)).
This is nothing but Eq. (5) in the problem.

(b)

For p = 1 GeV, ((Subscript[m, 3])^2) - ((Subscript[m, 2])^2) = 
        2.5×(10^(-3)) (eV^2)/(c^2), and using ℏ = 6.582×(10^(-22)) MeV sec = 
          6.582×(10^(-16)) eV sec and c = 3.00×(10^10)  cm (sec^(-1)), the argument of (sin^2)is \(TraditionalForm\`\(2.5×10\^\(-3\) eV\^2/c\^2 \( c\^3\) t\)\/\(c\ 4\ \(10\^9\) eV\ 6.582\ \(10\^\(-16\)\) eV\ sec\) = \(\(950\ \ t\)\/sec = \(\(950\ c\ t\)\/\(c\ sec\) = \(\(317\ L\)\/\(\(10\^10\) cm\) = 3.17×10\^\(-3\) \(\(L\/km\)\(.\)\)\)\)\)\)
Taking (Subscript[θ, 23]) = (1/(Sqrt[2]))and hence (sin^2) 2  (Subscript[θ, 23]) = 1, we find
(Subscript[P, surv]) = 1 - (sin^2)(
              3.17×(10^(-3)) (L/km)) = (cos^2)(
            3.17×(10^(-3)) (L/km)).

Plot[Cos[3.17 (10^(-3)) L]^2, {L, 0, 10000}]

The oscillation occurs on the distances of thousands of kilometers.  This is what had been observed by the SuperKamiokande experiment, Phys. Rev. 93, 101801 (2004), even though the data shows nearly washed-out oscillation due to the so-so resolution in energy and distance.  A very macroscopic quantum phenomenon.

(c)

For the three-state case, we keep all angles, and we calculate the oscillation probability from (Subscript[ν, μ])to (Subscript[ν, e]),

Amue = (E^(-I c p t/ℏ)) ({1, 0, 0} . U . 
              DiagonalMatrix[{(E^(-I ((Subscript[m, 1])^2) (c^3) t/(2  p ℏ))), (E^(-I ((Subscript[m, 2])^2) (c^3) t/(2  p ℏ))), 
                  (E^(-I ((Subscript[m, 3])^2) (c^3) t/(2  p ℏ)))}] . Udagger . {{0}, {1}, {0}})[[1]]

Amuestar = % /. {t -> -t, δ -> -δ}

Pmue = ExpToTrig[Expand[Amue Amuestar]]

and from (Subscript[ν, e])to (Subscript[ν, μ]),

Aemu = (E^(-I c p t/ℏ)) ({0, 1, 0} . U . 
              DiagonalMatrix[{(E^(-I ((Subscript[m, 1])^2) (c^3) t/(2  p ℏ))), (E^(-I ((Subscript[m, 2])^2) (c^3) t/(2  p ℏ))), 
                  (E^(-I ((Subscript[m, 3])^2) (c^3) t/(2  p ℏ)))}] . Udagger . {{1}, {0}, {0}})[[1]]

Aemustar = % /. {t -> -t, δ -> -δ}

Pemu = ExpToTrig[Expand[Aemu Aemustar]]

Simplify[Pmue - Pemu]

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

(d)

As we did in the class, the time-reserval invariance states that \(TraditionalForm\`⟨a | e\^\(\(-i\)\ H\ t/ℏ\) | 
            b⟩ = ⟨\(b\&~\) | 
            e\^\(\(-i\)\ H\ t/ℏ\) | \(a\&~\)⟩\), and hece\(TraditionalForm\`P(a -> b) = P(b\&~ -> a\&~)\). Here, \(TraditionalForm\`a\&~\) is the time-reversed state of a, and \(TraditionalForm\`b\&~\) is the time-reversed state of b.  In our case, the time-reversed state has the opposite momentum.  (If you worry about the helicity, both the momentum and the spin flip, and remains the same, namely left-handed.)  However, the Hamiltonian depends only on the magnitude of the momentum and hence the oscillation probabilities also depend only on the magnitude of the momentum.  Therefore \(TraditionalForm\`P(b\&~ -> a\&~) = P(b -> a)\).  Choosing (Subscript[ν, e])and (Subscript[ν, μ])of the same momenta to be the initial and final states, we find that the time reversal invariance would predict P((Subscript[ν, μ]) -> (Subscript[ν, e])) = 
        P((Subscript[ν, e]) -> (Subscript[ν, μ])) .

The fact that the time-reversal violation requires δ ≠ 0makes sense because the time-reversal opeator involves the complex conjugation, and δ is the only complex parameter in the Hamiltonian.