HW 6

1. WKB harmonic oscillator

a) Energy levels

(Hitoshi does this problem in his WKB notes, on page 8.)

The classical turning points a, b are where KE=0, i.e.
	E = V(x) = (1/2) m(ω^2) (x^2)
so
	a,  
      b  =  ∓(Sqrt[((2  E)/(m (ω^2)))]).
We apply or V(x) and the turning points to take the energy condition
	\(TraditionalForm\`∫\_a\%b\(\@\( 2 
                   m[E - V(x)]\)\) dx\  = \ \((n + 
              1\/2)\) π\ \)
and integrate the l.h.s.:

a0 = -(Sqrt[((2  E0)/(m (ω^2)))]); 
    b0 = -a0;
V0[x_] = ((m (ω^2) (x^2))/2);
Integrate[Sqrt[2  m (E0 - V0[x])], {x, a0, b0}, 
      Assumptions -> {m > 0, ω > 0, E0 > 0}]

We find
	\(TraditionalForm\`\(E\ π\)\/ω\  = \ \((n + 
              1\/2)\) π\ \)
such that

\(TraditionalForm\`E\  = \ \((n + 1\/2)\) \ ω\)
as desired.

b) SHO WKB wavefunctions

(Let us proceed in the units \(TraditionalForm\` = \(m = \(ω = 1\)\)\).)

The trick here is that we want to use the classically-allowed WKB wavefunction for SHO in between the turning points, and the "tunneling" WKB wavefunctions in the forbidden regions, which extend to ±∞.  This can be seen easily on a plot of the ground state, where the turning points are at ± 1:

Plot[{(π^(-1/4)) Exp[((-(x^2))/2)], 
          ((x^2)/2)}, {x, -2, 2}, 
        PlotStyle -> {Dashing[{}], Dashing[{ .05,  .05}]}];

Our SHO energies (from part (a)), potential function, and classical turning points are:

En[n_] = n + (1/2);
V[x_] = ((x^2)/2);
a[n_] = -(Sqrt[2 En[n]]);
b[n_] = -a[n];

We construct the wavefunction in allowed region using Eq. 30 in the notes.  We split up the integral around x = 0 and consider the two turning points a and b separately, so Mathematica doesn't give a trivial result:

ψa1[x_, 
          n_] = (((((((2  Abs[V'[a[n]]]))^(1/3))/(π (Sqrt[ 2 (En[n] - V[x])])))))^(1/2)) 
          Cos[Integrate[(Sqrt[2 (En[n] - V[t])]), {t, a[n],  x}, 
                Assumptions -> {n >= 0, a[n] < x, 
                    x < 0}] - (π/4)];

ψa2[x_, 
          n_] = (((((((2  Abs[V'[b[n]]]))^(1/3))/(π (Sqrt[ 2 (En[n] - V[x])])))))^(1/2)) 
          Cos[Integrate[(Sqrt[2 (En[n] - V[t])]), {t, x, b[n]}, 
                Assumptions -> {n >= 0, x > 0, 
                    x < b[n]}] - (π/4)];

Similarly, for the forbidden regions we use Eq. 31:

ψf1[x_, 
          n_] = (((((((2  Abs[V'[a[n]]]))^(1/3))/(π (Sqrt[ 2 (V[x] - En[n])])))))^(1/2)) (1/2) 
          Exp[-Integrate[(Sqrt[2 (V[t] - En[n])]), {t, x, a[n]}, 
                Assumptions -> {n >= 0, 
                    x < a[n]}]];
ψf2[x_, 
          n_] = (((((((2  Abs[V'[b[n]]]))^(1/3))/(π (Sqrt[ 2 (V[x] - En[n])])))))^(1/2)) (1/2) 
          Exp[-Integrate[(Sqrt[2 (V[t] - En[n])]), {t, b[n], x}, 
                Assumptions -> {n >= 0,  x > b[n]}]];

Ok, now what about the regions right around the classical turning points a and b?  We know that the wavefunctions above blow up near these points.  So, let's define a radius about the turning points inside which we treat specially:

ϵ = 0.5;

Next, we construct a piecewise function from the asymptotic wavefunctions above outside these regions (using 'UnitStep'):

ψwkbasy[x_, 
          n_] = (((-1))^n) ψf1[x, 
              n] (1 - 
                UnitStep[
                  x - (a[
                          n] - ϵ)]) + (((-1))^n) ψa1[x, n] 
            UnitStep[
              x - (a[n] + ϵ)] (1 - UnitStep[x]) + ψa2[
              x, n] UnitStep[
              x] (1 - UnitStep[x - (b[n] - ϵ)]) + ψf2[
              x, n] UnitStep[x - (b[n] + ϵ)];

Now, in the special regions, we use the Airy functions themselves, and again construct a piecewise function:

ψwkbtpa[x_, 
          n_] = (((-1))^n) 
            AiryAi[-(((2  Abs[V'[a[n]]]))^(1/3)) (x - a[n])] 
            UnitStep[
              x - (a[n] - ϵ)] (1 - 
                UnitStep[x - (a[n] + ϵ)]) + 
          AiryAi[(((2  Abs[V'[b[n]]]))^(1/3)) (x - b[n])] 
            UnitStep[
              x - (b[n] - ϵ)] (1 - 
                UnitStep[x - (b[n] + ϵ)]);

Finally, we add these two chunks to get the whole function:

ψwkb[x_, 
          n_] = ψwkbasy[x, n] + ψwkbtpa[x, n];

WAIT:  If we're using the asymptotic wavefunctions mated to the actual Airy functions in the transition regions, why not just use the Airy functions everywhere, expanding about the turning points?  A moment's reflection gives the answer:  this would work great coming in from ±∞ and passing through the turning points, but as you move into the classically-allowed region, the Airy function would not reach its cosine-shaped asymptote fast enough!  The reason we do it piecewise is that we want accurate results where E ≫ V, and live with nastiness in the transition regions.  (In the classically-forbidden regions, the Airy function should decay fast enough to match the asymptotic solution.)

To see this, let's construct this all-Airy function wavefunction by expanding around the two turning points:

ψair[x_, 
          n_] = (((-1))^n) 
            AiryAi[-(((2  Abs[V'[a[n]]]))^(1/3)) (x - 
                    a[n])] (1 - UnitStep[x]) + 
          AiryAi[(((2  Abs[V'[b[n]]]))^(1/3)) (x - b[n])] 
            UnitStep[x];

Finally, we wish to compare to the exact results:

ψexa[x_, 
          n_] = (1/(Sqrt[(2^n) Factorial[n]])) (π^(-1/4)) 
          Exp[((-(x^2))/2)] HermiteH[n, x];

Let's compare the exact and WKB solutions for n = 1:

Plot[{ψexa[x, 1], ψwkb[x, 1]}, {x, -4, 4}, 
        PlotStyle -> {GrayLevel[0], Hue[0]}];

Our WKB solutions looks great at |x| -> ∞ and x ≈ 0 where E ≪ V and E ≫ V respectively, as expected.  The Airy function matches nicely with the decaying wavefunctions in the forbidden regions, but does indeed mismatch the wavefunction in the allowed region.

Now, let's compare the exact and all-Airy solutions:

Plot[{ψexa[x, 1], ψair[x, 1]}, {x, -4, 4}, 
        PlotStyle -> {GrayLevel[0], Hue[0]}];

As expected, we get the opposite behavior.  Coming in from ±∞, the wavefunction matches the exact solution nicely.  However, as we move toward x = 0, the discrepancy builds up until you get a mismatch.

Now let's look at n = 10:

Plot[{ψexa[x, 10], ψwkb[x, 10]}, {x, -7, 7}, 
        PlotStyle -> {GrayLevel[0], Hue[0]}];

Fantastic!  Because the energy is much higher, the approximation is far more accurate (i.e, the state is much more classical), and on this plot the mismatch between the Airy function in the transition region and the wavefunction in the allowed region is not even visible.

Now let's look at the all-Airy solution for n = 10:

Plot[{ψexa[x, 10], ψair[x, 10]}, {x, -7, 7}, 
        PlotStyle -> {GrayLevel[0], Hue[0]}];

A bit better than n = 1.  The discrepancy builds up again coming toward x = 0 from the turning points; however, because n is even, there is no mismatch at x = 0, just a sharp point.

Finally, let us look at n = 20:

Plot[{ψexa[x, 20], ψwkb[x, 20]}, {x, -10, 10}, 
        PlotStyle -> {GrayLevel[0], Hue[0]}];

We see that the WKB solution has a bit of overshoot due to the higher energy; otherwise the agreement is quite good.

Finally, let us compare the exact and all-Airy solutions for n = 20:

Plot[{ψexa[x, 20], ψair[x, 20]}, {x, -10, 10}, 
        PlotStyle -> {GrayLevel[0], Hue[0]}];

The overshoot is still there, but the discrepancy is less approaching x = 0 as one would expect with the higher energy.

N.B.:  Of course, one could just use polynomial to handle the transition regions.  But, if we're going to use numerical methods, we might as well use the variational method instead of WKB.  :P

2. Classical limit of hydrogen atom

a) Correspondence between orbital motion and emitted photon

One can derive the Bohr model by taking the energy of an electron circling a proton (CGS units)
	E =  (1/2) m (v^2)  -  ((Z (e^2))/r)
and applying the quantization condition
	\(TraditionalForm\`L\_n = \(m\ v\ r = \ n\ \)\).
(This justification for this quantization is just units; we know that the ground state in fact has no angular momentum.  One can also derive the Bohr model by requiring whole wavelengths to fit on a circle of radius r.)
Substituting for v, we get
	\(TraditionalForm\`E = \(1\/2\) 
            m\ \((\(n\ \)\/\(m\ r\))\)\^2\  - \ \(Z\ e\^2\)\/r\).
We require one more condition to eliminate r.  Let us use the force balance
	F = ((m (v^2))/r) = ((Z (e^2))/(r^2))
which gives
	\(TraditionalForm\`r = \(\(Z\ e\^2\)\/\(m\ v\^2\) = \(\(Z\ e\^2\)\/m\) \((\(m\ r\)\/\(n\ \))\)\^2\)\)
and finally
	\(TraditionalForm\`r\_n = \ \(n\^2\ \^2\)\/\(Z\ \(e\^2\) m\)\).
Substituting this back in to the energy,
	\(TraditionalForm\`E\_n = \(\(1\/2\) \(\(n\^2\) \^2\)\/\(m\ r\_n\^2\) - \(Z\ e\^2\)\/r\_n = \ \(\(\(n\^2\) \^2\)\/\(2  
                    m\)\) \((\(Z\ \(e\^2\) m\)\/\(\(n\^2\) \^2\))\)\^2 - Z\ \(\(e\^2\)(\(Z\ \(e\^2\) m\)\/\(\(n\^2\) \^2\))\)\)\)
to obtain
	\(TraditionalForm\`E\_n = \(-\(\(\(Z\^2\) \(e\^4\) 
                m\)\/\(2 \( n\^2\) \^2\)\)\)\).
Similarly,
	\(TraditionalForm\`v\_n = \ \(\(n\ \)\/\(m\ r\_n\) = \ \(\(\(n\ \)\/m\) \(Z\ \(e\^2\) m\)\/\(\(n\^2\) \^2\) = \(Z\ e\^2\)\/\(n\ \)\)\)\).

Having derived these basic quantities, we can express the frequency of orbital motion in state n as
	\(TraditionalForm\`ν\_o\  = \ \(v\_n\/\(2  π\ r\_n\) = \(\(1\/\(2  π\)\) \(Z\ e\^2\)\/\(n\ \)\ \(Z\ \(e\^2\) m\)\/\(\(n\^2\) \^2\) = \ \(1\/\(2  π\)\) \(\(Z\^2\) \(e\^4\) m\)\/\(\(n\^3\) \^3\)\)\)\).
On the other hand, we know that the frequency of a photon due to a change from state n to n - k is
	\(TraditionalForm\`ν\_γ = \(\(E\_n - E\_\(n - k\)\)\/\(2  π\ \) = \ \(\(-\(1\/\(2  π\)\)\)\ \ \(\(\(Z\^2\) \(e\^4\) 
                    m\)\/\(2  \^3\)\) \((1\/n\^2 - 
                  1\/\((n - k)\)\^2)\) = 
            1\/\(2  π\)\ \ \(\(\(Z\^2\) \(e\^4\) 
                    m\)\/\(2 \( 
                        n\^2\) \^3\)\) \((1\/\((1 - k\/n)\)\^2 - 
                  1)\)\)\)\).
Approximating k ≪ n,
	\(TraditionalForm\`ν\_γ\  ≃ \ 
          1\/\(2  π\)\ \ \(\(\(Z\^2\) \(e\^4\) 
                  m\)\/\(2 \( n\^2\) \^3\)\) \((1 + \(2  k\)\/n - 
                1)\) = \ \((\(1\/\(2  π\)\) \(\(Z\^2\) \(e\^4\) m\)\/\(\(n\^3\) \^3\))\) \((k)\)\).

As indicated in the hint, we find that the frequency of the photon (Subscript[ν, γ]) is an integer multiple of the frequency of orbital motion (Subscript[ν, o]).

b) Correspondence between classical radiation and mean lifetime

Under classical dynamics the electron is being continually accelerated, so it should emit radiation.  To know the mean lifetime τ of a state n, we must know how fast it loses power.  Since the model is semiclassical, we can use a classical radiation formula.  We found above that the velocity goes like (1/n), so the motion at large n is non-relativistic and so we can just use the Larmor power formula (again in CGS)
	P = ((2  (e^2) (a^2))/(3  (c^3))).
From above,
	\(TraditionalForm\`a\_n = \ \(F\_n\/m = \ \(\(Z\ e\^2\)\/\(m\ r\_n\^2\) = \(\(\(Z\ e\^2\)\/m\) \((\(Z\ \(e\^2\) m\)\/\(\(n\^2\) \^2\))\)\^2 = \(\(Z\^3\) \(e\^6\) m\)\/\(\(n\^4\) \^4\)\)\)\)\).
Substituting,

	\(TraditionalForm\`P = \(2  Z\^6\ \(e\^14\) m\^2\)\/\(3 \( c\^3\) \(n\^8\) \^8\)\).
	
The mean lifetime should be

	\(TraditionalForm\`τ = \(\(E\_n - E\_\(n - 1\)\)\/P = \ \(\((\(\(Z\^2\) \(e\^4\) m\)\/\(\(n\^3\) \^2\))\) \((1)\)\ \(3 \( c\^3\) \(n\^8\) \^8\)\/\(2  Z\^6\ \(e\^14\) m\^2\) = \(3 \( c\^3\) \(n\^5\) \^6\)\/\(\(Z\^4\) \(e\^10\) m\)\)\)\).
	
Rearranging,

\(TraditionalForm\`1\/τ = \(2\/3\) \(e\^2\/c\) \(\((\(Z\ e\^2\)\/c)\)\^4\) \(\(m\ c\^2\)\/\) 1\/n\^5\)

as expected.

c) Comparison of classical and quantum lifetimes

We implement the correct expression above (with  in units of eV s):

lifetimec[
          n_] = ((((2/3) α (((Z α))^4) (mec2/hbar) (1/(n^5))))^(-1)) /. {α -> (137^(-1)), Z -> 1, 
            mec2 -> 511000, 
            hbar -> 6.582*(10^(-16))};

Computing for n = 2, 4, 6 we find respectively

{2, 4, 6} // lifetimec

These are accurate to within an order of magnitude of the given quantum values, and seem to improve with increasing n as we would expect.

3. Spin precession in magnetic field

a) Time evolution of eigenstates

The Schrödinger equation is

	\(\(\(\(i\  ∂\_t\  | ψ\)\(\ \)\(⟩\)\) = 
                          H | ψ\)\(\ \)\(⟩\)\)=(-g (e/(2  m c))S·B)|ψ ⟩=(-g ((e B)/(2  m c))(Subscript[S, z]))|ψ ⟩=-ω (Subscript[S, z]) | ψ ⟩ 

where we defined ω = g ((e B)/(2  m c)).
	
Then, for \(TraditionalForm\`\(\(\(\(\(|\)\(ψ\)\)\(\ \)\(⟩\)\)\  = \ \(\(|\)\(S\_z\)\); +\)\(\ \)\(⟩\)\)\) we have

	\(TraditionalForm\`\(\(\(\(i\  ∂\_t\  | 
              S\_z; +\)\(\ \)\(⟩\)\)\  = \ \(-\(ω\/2\)\) | 
            S\_z; +\)\(\ \)\(⟩\)\)\)
	
and for \(TraditionalForm\`\(\(\(\(\(|\)\(ψ\)\)\(\ \)\(⟩\)\)\  = \ \(\(|\)\(S\_z\)\); -\)\(\ \)\(⟩\)\)\) we have

	\(TraditionalForm\`\(\(\(\(i\  ∂\_t | 
              S\_z; -\)\(\ \)\(⟩\)\) = \ ω\/2 | S\_z; -\)\(\ \)\(⟩\)\)\).
	
So we have for the solutions

|(Subscript[S, z]); + , 
            t ⟩ =  
          (e^(+i (ω/2) t)) | 
            (Subscript[S, z]); +  ⟩
	|(Subscript[S, z]); - , 
            t ⟩ =  
          (e^(-i (ω/2) t)) | 
            (Subscript[S, z]); -  ⟩.

For a state aligned with the quantization axis, the "precession" is just a rotating phase.

b) Dynamics of \(TraditionalForm\`\(\(\(\(|\)\(S\_x\)\); +\)\(\ \)\(⟩\)\)\)

First, define the spin matrices:

\(\(S\_x = \(\/2\) {{0, 1}, {1, 0}};\)\)
\(\ \(S\_y = \(\/2\) {{0, \(-I\)}, {I, 
              0}};\)\)
\(\ \(S\_z = \(\/2\) {{1, 0}, {0, \(-1\)}};\)\)

Find the eigenstates of (Subscript[S, x]):

eigenx = Eigensystem[(Subscript[S, x])]

Our state \(TraditionalForm\`\(\(\(\(|\)\(S\_x\)\); +\)\(\ \)\(⟩\)\)\) is then

	\(TraditionalForm\`\(\(\(\(\(\(\(|\)\(S\_x\)\); \ +\)\(\ \)\(⟩\)\)\  = \ \(\(1\/\@2\ \((\(\(|\)\(S\_z\)\); +\ \ ⟩\)\)\(+\)\) | 
              S\_z; -\)\(\ \ \)\(⟩\)\)\(\ \)\()\)\)\)
	
and so

	\(TraditionalForm\`\(\(\(\(\(\(\(|\)\(S\_x\)\); \ +, 
              t\)\(\ \)\(⟩\)\)\  = \ \(\(1\/\@2\ \((\(\(|\)\(S\_z\)\); +, 
                    t\ \ ⟩\)\)\(+\)\) | 
              S\_z; -\ , 
          t\)\(\ \)\(⟩\)\)\(\ \)\()\)\)\).
	
In component form on the time-independent basis \(TraditionalForm\`\(\(\({\(\(|\)\(S\_z\)\); +\ \ ⟩, \(\(|\)\(S\_z\)\); -\)\(\ \ \)\(⟩\)\)\(}\)\)\),

|(Subscript[S, x]);  +, 
            t ⟩ ≐ (1/(Sqrt[2]))(({{(e^(+i (ω/2) t))}, {(e^(-i (ω/2) t))}})) .

c) Showing precession

Let us implement the above state in component form, as well as its conjugate-transpose:

ψsxp[
          t_] = (1/(Sqrt[2])) {{Exp[
                I (ω/2) t]}, {Exp[-I (ω/2) 
                  t]}};
ψsxpct[
          t_] = (1/(Sqrt[2])) {{Exp[-I (ω/2) t], 
              Exp[I (ω/2) t]}};

We calculate ⟨ (Subscript[S, x]);  +,  
        t  | (Overscript[S,  ->])   | 
          (Subscript[S, x]);  +, t ⟩ and put it into vector {x, y, z} form:

ExpToTrig[{(ψsxpct[t] . (Subscript[S, x]) . ψsxp[t])[[1, 
            1]], (ψsxpct[t] . (Subscript[S, y]) . ψsxp[t])[[1, 
            1]], (ψsxpct[t] . (Subscript[S, z]) . ψsxp[t])[[1, 
            1]]}]

This is the clockwise precession about the z-axis with angular frequency ω  =  g ((e B)/(2  m c)) of a vector of length \(TraditionalForm\`\/2\).  As expected, there is no z-component.  (The classical analogue is a gyroscope precessing clockwise in the plane; of course, here we can't simultaneously have components of x/y and z.)

4. Hamilton-Jacobi for light [optional]

a) Apply WKB approximation to equation of motion

Following the same steps we did for the Schrödinger equation, we first write the Maxwell's equation with (A^0) = (e^(i S/ℏ)):

	\(\((\(n\^2\/c\^2\) ∂\^2\/∂\ t\^2 - \(∇\& -> \)\^2)\) e\^\(i\ S/ℏ\)\)=(\(\(n\^2\/c\^2\) \((\(i\ S\&\(\(·\)\(·\)\)\)\/ℏ + \((\(i\ \(S\&·\)\)\/ℏ)\)\^2)\)\)-((i (((Overscript[∇,  ->]) )^2) S)/ℏ)-(((i (Overscript[∇,  ->]) S)/ℏ))2)(e^(i S/ℏ))=0 .
	
In the limit S ≫ ℏ, we can drop terms of O(S/ℏ) and keep those of ((O(S/ℏ))^2):

((n^2)/(c^2)) (((Overscript[S, ·]))^2)-((Overscript[∇,  ->]) S)2=0.
	
This is the "Hamilton-Jacobi" equation.

(Of course, in the Maxwell equation, there is no notion of ℏ.  What we are doing is a valid approximation when the variation of the phase is very fast compared to the variation of the index of refraction.  It is called "eikonal approximation" in optics.)

b) Equivalence to particle in given potential

The Hamiltonian of free particle in potential V((Overscript[x,  ->]) ) = -(1/(2  
                  m)) ((n((Overscript[x,  ->]) ))^2) is
	H = ((((Overscript[p,  ->]) )^2)/(2  m)) - 
          ((n^2)/(2  m)).
	
Applying this to the Hamilton-Jacobi equation, we obtain

	(1/(2  m)) ((((Overscript[∇,  ->]) S))^2) - 
          ((n^2)/(2  m)) + ∂(S/∂)t = 0.
	
Since the action has no explicit time dependence, we can make the Legendre transformation \(TraditionalForm\`S(
          t, \(x\& -> \)) = \(S\&~\)(E, \(x\& -> \)) - E\ t\) to get
	\(TraditionalForm\`\(1\/\(2  m\)\) \((\(∇\& -> \)S\&~)\)\^2 - 
          n\^2\/\(2  m\) = \(E = 0\)\)
	
and finally

\(TraditionalForm\`n\^2 - \((\(∇\& -> \)\(S\&~\))\)\^2 = 
        0\).

If n in part (a) has no time dependence, we can do the same Legendre transformation and obtain the equivalent "H-J" equation

\(TraditionalForm\`\((\(E\/c\) n)\)\^2 - \((\(∇\& -> \)\(S\&~\))\)\^2 = 0\).

c) Separate variables, integrate "action" variable

Assuming n((Overscript[x,  ->]) ) = n(x), the "Hamilton-Jacobi" equation does not have an explicit dependence on y or t.  Writing

	S(x, y, 
          t) = (Overscript[S,  ∼]) (x, (Subscript[p, y]), E) + (Subscript[p, y]) y - E t,
	
the equation becomes

	((n^2)/(c^2)) (E^2)-((d/(d x)) (Overscript[S,  ∼]) )2-((Subscript[p, y])^2)=0.
	
Therefore,

	\(TraditionalForm\`\(d\ \(S\& ∼ \)\)\/\(d\ x\) = \@\(\(n\^2\/c\^2\) E\^2 - p\_y\^2\)\)
	
and hence

\(TraditionalForm\`\(S\& ∼ \) = ∫\^x\(\@\(\(\( n(x')\)\^2\/c\^2\) E\^2 - p\_y\^2\)\) d\ x'\).
	
Here, the possible x-dependence of the index of refraction is emphasized.

d) Integral expressions for "angle" variables

Using \(TraditionalForm\`t = ∂\(\(\(\ \)\(S\)\)\& ∼ \)\/∂\ E\), \(TraditionalForm\`y = \(-\(∂\(\(\(\ \)\(S\)\)\& ∼ \)\/∂\ p\_y\)\)\), we find

t =  (∫^x)((((( n(x'))^2)/(c^2)) 
                  E)/(Sqrt[(((n(x'))^2)/(c^2)) (E^2) - ((Subscript[p, y])^2)])) 
            d x',
	\(TraditionalForm\`y = ∫\^x\ \(p\_y\/\@\(\(\(n(x')\)\^2\/c\^2\) E\^2 - p\_y\^2\)\) d\ x'\).

e) Work to Snell's law

Assuming that n(x) = (Subscript[n, 1]) for x < 0 and n(x) = (Subscript[n, 2]) for x > 0, and choosing the lower end of the integration at x = 0,

	\(TraditionalForm\`\(y\_ < \)(
          x) = \ \(p\_y\/\@\(\(n\_1\^2\/c\^2\) E\^2 - p\_y\^2\)\) x\) for x < 0
	
and

	\(TraditionalForm\`\(y\_ > \)(
          x) = \ \(p\_y\/\@\(\(n\_2\^2\/c\^2\) E\^2 - p\_y\^2\)\) x\) for x > 0.

Using the trigonometric relation sin α = 
        (1/(Sqrt[(cot^2) α + 1])), we find

	\(TraditionalForm\`sin\ \(α\_ < \) = \(1\/\@\(\(\(n\_1\^2\/c\^2\) E\^2 - p\_y\^2\)\/p\_y\^\(\(2\)\(\ \)\) + 1\) = \(c\ p\_y\)\/\(\(n\_1\) E\)\)\),
	\(TraditionalForm\`sin\ \(α\_ > \) = \(1\/\@\(\(\(n\_2\^2\/c\^2\) E\^2 - p\_y\^2\)\/p\_y\^\(\(2\)\(\ \)\) + 1\) = \(c\ p\_y\)\/\(\(n\_2\) E\)\)\),

and hence

\(TraditionalForm\`\(sin\ \(α\_ < \)\)\/\(sin\ \(α\_ > \)\) = n\_2\/n\_1\),
	
which is nothing but the Snell's law of refraction.