Legal Notice, Disclaimer, and Book Ordering Information:

(* Legal Notice, Disclaimer, and Book Ordering Information:

Last Modified: FIO June 8, 2002

Mathematica For Physics: 2nd Edition

Robert L. Zimmerman and Fredrick I. Olness

WebSite:  
  www.physics.smu.edu/~olness
  darkwing.uoregon.edu/~phys600/

MathSource Number: 0206-862

ISBN 0-201-53796-6
For ordering information, call 1-800-282-0693
Addison-Wesley Publishing Company

Communication with the authors:
Fredrick I. Olness:  olness@mail.physics.smu.edu
Robert L. Zimmerman: bob@zim.uoregon.edu

Copyright 2002, Addison-Wesley Publishing Company, Inc. The material in
this file may be distributed freely  so long as the content remains
unchanged, and the copyright and  reference notices are included. The
publisher grants permission for the noncommercial use of these programs
and program segments.  All other uses require the prior written consent
of the publisher.    

All rights are reserved with regard to the  material in the
"Mathematica For Physics" text, and may not be reproduced, stored in a
retrieval system, or transmitted,  in any form or by any means,
electronic, mechanical, photocopying,  recording, or otherwise, without
the prior written permission  of the publisher. 

The programs and applications presented in this book have been included
for  their instructional value. They have been tested with care but are
not  guaranteed for any particular purpose. The publisher does not
offer any  warranties or representations, nor does it accept any
liabilities with  respect to the programs or applications. 

 *)

Contains correction on:  
    2nd Edition:Chpt.9,Problem 9.1.5 p.547: The One-Dimensional Harmonic Oscillator  (FIO, 19 Apr 2004)  
      Part b. Definition of Table was previoiusly missing. 
      
2nd Edition:Chpt.9,Problem 9.3.4 p.569: Particle in a Finite, Spherical Well  (FIO, 19 Apr 2004)  
     Part b. Mathematica Version 5 has difficult with the "ℓ" character. Replace with "ell"

9. Quantum Mechanics

9.1  Introduction

Off[General::"spell"]; 
Off[General::"spell1"];

9.2 One-Dimensional Schroedinger's Equation

Problem 1: Particle bound in an Infinite Potential Well

Required Packages

Needs["Graphics`Legend`"]

Solution:

Clear["Global`*"];

Part a

eq1= 0==(V - En) ψ[x]-(((ℏ^2) ψ''[x])/(2 m));

eq2 = eq1 /. V -> 0

krule =  En -> (((ℏ^2) (k^2))/(2 m));

schrod = Solve[eq2 /. krule, ψ''[x] ][[1, 1]] /. 
        Rule -> Equal

dsol= DSolve[schrod,ψ,x][[1,1]]

Part b

sol1 = 
      Reduce[{ψ[0] == 0, ψ[a] == 0} /. dsol]

En==(((k^2) (ℏ^2))/(2 m))/.k -> ((n π)/a)

ψ[x_, n_] =ψ[x ]  /. dsol/.{C[1] -> 0, k -> ((n π)/a)}

normEq =  1  ==  Integrate[ψ[x, n]^2, {x, 0, a}] // 
          Simplify[#, Element[n, Integers]] & // Solve[#, C[2]][[2, 
              1 ]] &

ψ[x_, n_]  = ψ[x, n] /. normEq // Simplify

Part c

table= Table[Integrate[ψ[x,i]ψ[x,j],{x,0,a}],{i,4},{j,4}];

TableForm[table, 
      TableHeadings -> {Table[ψ[i], {i, 4}], 
          Table[ψ[j], {j, 4}]}]

Plot[ 
        Table[ψ[n, x], {n, 4}] /. {a -> 1} // 
          Evaluate, 
{x, 0, 1}, 
        PlotStyle -> {{}, Dashing[{ .05}], Dashing[{ .03}], 
            Dashing[{ .01}]}, 
PlotLabel -> 
          StyleForm[ψ, FontSize -> 14, 
            FontWeight -> 
              Bold], 
PlotLegend -> {"n=1", "n=2\
", "n=3", "n=4" }, 
        LegendPosition -> {1, -1/2}];

Problem 2: Particle bound in a Finite Potential Well

Solution

Clear["Global`*"];

Part a

schrodeq=0==(V - En) ψ[x]-(((ℏ^2) ψ''[x])/(2 m));

EWrule={En -> (((kW^2) (ℏ^2))/(2 m)) + V};

ψ[x] /. DSolve[schrodeq, ψ[x], x][[1, 1]] /. 
            EWrule // ExpToTrig  // Simplify  // PowerExpand

kWrule= kW -> (((Sqrt[2]) (Sqrt[m (En - V)]))/ℏ)/.V -> -V0/.En -> -Wn

ψW[x_] = cSym Cos[kW x] + cAsym Sin[kW x];

ELRrule={En -> (((kLR^2) (ℏ^2))/(2 m))};

ψ[x] /. DSolve[schrodeq, ψ[x], x][[1, 1]] /. 
            V -> 0 /. ELRrule  // 
        Simplify// PowerExpand

kLRrule= kLR -> (((Sqrt[2]) (Sqrt[En m]))/ℏ) /.En -> -Wn //PowerExpand

qLRrule= qLR -> (((Sqrt[2]) (Sqrt[m Wn]))/ℏ);

ψ[
                  x] /. DSolve[schrodeq, ψ[x], x][[1, 1]] /. 
              V -> 0 /. ELRrule  /. kLR -> I qLR // 
        Simplify// PowerExpand

ψR[x_] =  
        cR (E^(-qLR x))   ;

ψL[x_] =  
        cL (E^(+qLR x)) ;

Part b

eq1= {  (ψL[ x]-ψW[ x]==0)  /. {x-> -a}, 
        (ψW[ x]-ψR[ x]==0)  /. {x-> +a},
        (ψL'[x]-ψW'[x]==0)  /. {x-> -a},
        (ψW'[x]-ψR'[x]==0)  /. {x-> +a}  }

(eq2 = eq1 /. cAsym -> 0 ) // ColumnForm

eq3 = 
      Reduce[{ eq2 ,  cL ≠ 0, cR ≠ 0, 
            kW ≠ 0, qLR  ≠ 0,  cSym ≠ 0, 
            Cos[ a kW] ≠ 0 } // Flatten  , {cL, 
          cR} ] ; eq3 //ColumnForm

symSol =  Solve[eq3  , {cL, cR} ]  // 
          Simplify//Flatten

symEn = Tan[a kW] ==   qLR/kW;

eq4 = 
        Reduce[{ eq1 /. cSym -> 0 ,  cL ≠ 0, 
              cR ≠ 0, kW ≠ 0, qLR ≠ 0, 
              cAsym ≠ 0, Cos[ a kW] ≠ 0 } // 
            Flatten  , {cL, cR} ] ;
eq4 // ColumnForm

asymSol =  Solve[eq4 , {cL, cR} ]  // 
          Simplify//Flatten

asymEn = Tan[a kW] ==  - kW/qLR;

Part c

values = {a -> 1, m -> 1, ℏ -> 1, 
          V0  -> { 100,  200,  500, ∞}} ;

nRule=Wn -> 
            V0 - (((n^2) (π^2) (ℏ^2))/(8 (a^2) m));

kRules = {kW -> (( (Sqrt[2 m (V0 - Wn )]))/ℏ), qLR -> (( (Sqrt[2 (Wn) m]))/ℏ)};

eq5 =symEn/.kRules /.nRule//PowerExpand

\(\(pt1 = 
        Plot[{ArcTan[eq5[\([1]\)]], ArcTan[eq5[\([2]\)]]} /. values // 
            Evaluate, 
{n, 0, 
            9\ }, 
Frame -> True, 
          FrameTicks -> {{{1, "n=1"}, {3, "n=3"}, {5, "n=5"}, {7, "n=7"}, {9, "n=9"}}, {0, 0.5, 1, 
                1.5}}, 
PlotRange -> { .4, 
              1.6}, 
DisplayFunction -> 
            Identity];\)\)

text={
Text[ "V0=\!\(
\*InterpretationBox[\(∞\),
DirectedInfinity[ 1]]\)",  {6, 1.5}],
Text[ "V0=500",  {6, 1.3}],
Text[ "V0=200",  {6, 1.1}],
Text[ "V0=100",  {6, 0.8}]};

Show[pt1, Graphics[text] , GridLines -> Automatic, 
        DisplayFunction -> $DisplayFunction];

eq6 =asymEn/.kRules /.nRule//PowerExpand

text={
Text[ "V0=\!\(
\*InterpretationBox[\(∞\),
DirectedInfinity[ 1]]\)",  {7, 0}],
Text[ "V0=500",  {7, - .4}],
\(Text[\ "V0=200", \ {7, \(- .7\)}]\),Text[ "V0=100",  {7, - .9 }]};

\(\(Plot[{ArcTan[eq6[\([1]\)]], ArcTan[\ eq6[\([2]\)]]}\  /. values // 
          Evaluate, 
{n, 0.0\ , 
          9}, 
Frame -> True, \ 
        FrameTicks -> {{{2, "n=2"}, {4, "n=4"}, {6, "n=6"}, {8, "n=8"}}, {\(-1.5\), \(-1\), 0, 0.5, 1, 
              1.5}}, 
PlotRange -> {\(-1\)\ , \  .1}, \ 
GridLines -> 
          Automatic, 
Epilog -> \ 
          text\ \ ];\)\)

Part d

nValues[eq_, potential_, 
        guess_] := (n /. FindRoot[
              eq /. {m -> 1, a -> 1, ℏ -> 1, 
                    V0 -> potential} // Evaluate, {n, 
                guess}][[1]])

symGuess= {0.9, 2.9, 4.9, 6.9, 8.9};

symValues= {nValues[eq5,100,#]
           ,nValues[eq5,200,#]
           ,nValues[eq5,500,#]}&  /@  symGuess;

asymGuess = {1.9 , 3.9 , 5.9, 7.9 };

asymValues={nValues[eq6,100,#]
           ,nValues[eq6,200,#]
           ,nValues[eq6,500,#]}& /@  asymGuess;

\(\((Partition[Sort[{symValues, asymValues} // Flatten], 3] // 
        TableForm[#, TableSpacing -> {0, 2}, 
            TableHeadings -> {{"n=1", "n=2", "n=3", "n=4", "n=5", "n=6", "n=7", "n=8", "n=9"}, {"V0=100", "V0=200", "V0=500"}}] &)\)\)

\(Partition[Sort[{symValues, asymValues} // Flatten], 3] // \ 
      TableForm[#, 
TableSpacing -> {0, 2}, 
          TableHeadings -> {{"n=1", "n=2", "n=3", "n=4", "n=5", "n=6", "n=7", "n=8", "n=9"}, {"V0=100", "V0=200", "V0=500"}}
] &\)

Part e

\(\(\(symRules = {kRules, nRule, symSol, a -> 1, ℏ -> 1, 
          m -> 1, V0 -> 100, cSym -> 100, cAsym -> 0} // 
        Flatten\)\(
\)
    \)\)

Clear[sψ,aψ];
(* Symmetric*) 
sψ[x_ /;x<-1,n0_]:= (ψL[x]//.symRules//.{n->n0});
sψ[x_ /;-1<=x<1,n0_]:= (ψW[x]//.symRules//.{n->n0});
sψ[x_ /;x>=1,n0_]:= (ψR[x]//.symRules//.{n->n0});

symEnergy= nValues[eq5, 100,#]&/@symGuess;

plotsym= 
Plot[ sψ[x,#] //.symRules //Evaluate
        ,{x,-2,2}
        ,PlotLabel->  NumberForm[#,2]
        ,Axes->None
        ,DisplayFunction->Identity ]&  /@  symEnergy;

Show[GraphicsArray[{plotsym[[{1,2,3}]],plotsym[[{4,5}]]}]];

asymRules = {kRules, nRule, asymSol, a -> 1, ℏ -> 1, 
          m -> 1, V0 -> 100, cAsym -> 100, cSym -> 0} // 
        Flatten

(* Asymmetric*) 
aψ[x_ /;x<-1,n0_]:= (ψL[x]//.asymRules//.{n->n0});
aψ[x_ /; -1<=x<1,n0_]:= (ψW[x]//.asymRules//.{n->n0});
aψ[x_ /;x>=1,n0_]:= (ψR[x]//.asymRules//.{n->n0});

asymEnergy = nValues[eq6, 100, #] & /@ asymGuess ;

plotasym= 
  Plot[ aψ[x,#] //.asymRules //Evaluate
          ,{x,-2,2}
          ,PlotLabel->NumberForm[#,2]
          ,Axes->None
          ,DisplayFunction->Identity]& /@ asymEnergy;

Show[GraphicsArray[{plotasym[[{1,2}]],plotasym[[{3,4}]]}]];

Problem 3: Particle Hitting a Finite Step Potential

Required Packages

Needs["Graphics`Legend`"]

Solution

Clear["Global`*"];

Part a

ψL[
          x_] = (E^(I kL x)) + (E^(-I kL x)) R;

kLRule=kL -> (Sqrt[2]) ((Sqrt[En m])/ℏ);

ψR[
          x_] = (E^(I kR x)) T;

kRRule=kR->(((Sqrt[2]) (Sqrt[m (En - V0)]))/ℏ);

boundary = {ψL[x] == ψR[x], ψL'[
              x] == ψR'[x]} //. x -> 0

RTrule = Solve[boundary, {R, T}] // Flatten // Simplify

Part b

\(\(\(\ \ \)\(\(exp_\^*\) := \(\(exp\)\(\ \)\(/.\)\({Complex[a_, 
              b_] :> Complex[a, \(-b\)]}\)\(\ \)\)\)\)\)

\(\(\(\((a - I\ b + E\^\(\(-I\)\ θ\))\)\^*\)\(\ \)\)\)

Conjugate[
      a - I b + (E^(-I θ))]

∇f_ := D[f, # ] &   /@   {x, y, 
          z}

∇ψ[x  ]

\(flux[ψ_] := \(ℏ\/\(2  I\ m\)\) \((\(ψ\^*\)\ ∇ψ\  - \ ∇\(ψ\^*\)\ ψ)\)\  // Simplify\)

ψL[x] // flux//Simplify

incFlux = flux[ψL[x]] /.{R -> 
            0}

Rflux = (incFlux - flux[ψL[x]]) // Simplify

Tflux = flux[ψR[x]]

{RR,TT}= { Rflux[[1]] /incFlux[[1]],
           Tflux[[1]] /incFlux[[1]] }

{RR,TT}={RR,TT}/. RTrule

RR + TT  // Simplify

values = {m -> 1, ℏ -> 1, V0 -> 1};

\(\(Plot[\(\(\({RR, TT, TT + RR} /. RTrule\) /. kLRule\) /. kRRule\) /. 
            values // Evaluate, 
{En, 1, 
          2}, 
AxesLabel -> {"E", "Coefficient"}, 
PlotStyle -> {Dashing[{0.06, 
                0.02}], 
            Dashing[{0.02}], {}}, 
PlotLegend -> {"RR", "TT", "TT+RR"}, 
LegendPosition -> {1, \(-1\)/2}, 
PlotRange -> {0, 1}];\)\)

Part c

kRule = {kLRule, kRRule};

wave[x_, 
          En_] = {ψL[x], ψR[x]} /. RTrule /. kRule /. 
            values  // Simplify ;

ψ[x_, En_] := 
      Which[ x ≤ 0, wave[x, En][[1]], 
        x > 0, wave[x, En][[2]]]

waveplot2D[Energy_] := 
        Plot[{ψ[x, Energy] // Re, ψ[x, Energy] // 
              Im}, 
{x, -10, 
            10}, 
PlotStyle -> {{}, Dashing[{0.02}]}, 
          AxesLabel -> {"x", "y"}];

waveplot2D[1.1];

waveplot2D[0.9];

Part d

Plot3D[ψ[x, Energy] // Re, {x, -10, 10}, {Energy,  .5, 1.5}, 
        PlotPoints -> 50, 
        AxesLabel -> {"x", "Energy", "ψ"}];

Plot3D[ψ[x, Energy] // Im, {x, -10, 10}, {Energy,  .5, 1.5}, 
        PlotPoints -> 50, 
        AxesLabel -> {"x", "Energy", "ψ"}];

Part e

timePlot[time_, En_] := 
        Plot[{Exp[-I En time] ψ[x, En] // Re, 
            Exp[-I En time] ψ[x, En] // 
              Im}, 
{x, -10, 
            10}, 
PlotStyle -> {{}, {Thickness[
                  0.008]}}, 
PlotRange -> {-2, 
              2}];

timePlot[1. , 0.9];

frames = 1; 
Energy = 0.9;
timestep = ((2  π)/((frames) Energy));
\(Do[
\(timePlot[n\ timestep\ , 
          Energy];\)
, {n, 1, 
        frames}
]\)

frames = 1;
Energy = 1.1;
timestep = ((2  π)/((frames) Energy));
\(Do[
\(timePlot[n\ timestep\ , 
          Energy];\)
, {n, 1, 
        frames}
]\)

Problem 4: Particle Propagating Towards a Rectangular Potential

Solution

Clear["Global`*"];

Part a

ψL[
          x_] = (E^(I kLR x)) + (E^(-I kLR x)) R;

ψR[
          x_] = (E^(I kLR x)) T;

ψW[x_] = 
        B (E^(-I kW x)) + 
          A (E^(I kW x));

kRules={kW -> (Sqrt[2]) ((Sqrt[m (En - V0)])/ℏ),kLR->(Sqrt[2]) (((Sqrt[En m]))/ℏ)};

Part b

eq1 = {(ψL[x] - ψW[x] == 
                  0) /. {x -> -a}, 
     (ψW[x] - ψR[x] == 
                  0) /. {x -> +a}, 
     (ψL'[x] - ψW'[x] == 
                  0) /. {x -> -a}, 
     (ψW'[x] - ψR'[x] == 0) /. {x -> +a}} // 
          Simplify ;

ABRTrule = Solve[eq1, {A, B, R, T}] // FullSimplify // 
          Flatten;
ABRTrule // ColumnForm

Part c

\(\(\(\ \ \)\(\(exp_\^*\) := \(\(exp\)\(\ \)\(/.\)\({Complex[a_, 
              b_] :> Complex[a, \(-b\)]}\)\(\ \)\)\)\)\)

∇f_ := D[f, # ] &  /@  {x, y, z}

\(flux[ψ_] := \(ℏ\/\(2  I\ m\)\) \((\(ψ\^*\)\ ∇ψ\  - \ ∇\(ψ\^*\)\ ψ)\)\  // Simplify\)

incFlux=flux[ψL[x] /. R -> 0 /. 
          ABRTrule ]

Rflux=  (incFlux-flux[ψL[x]/.ABRTrule]) //ExpToTrig //Simplify

Tflux = flux[ψR[x] /. ABRTrule] // 
          ExpToTrig //Simplify

{RR, TT} = {Rflux[[1]]/incFlux[[1]], 
          Tflux[[1]]/incFlux[[1]]};
{RR, TT} // ColumnForm

TT + RR == 1   // Simplify

Part d

values = {ℏ -> 1, m -> 1, 
          a -> 1 };

Plot[{RR, TT} /. kRules /. {V0 -> -100} /. values // 
          Evaluate,  {En, 1, 
          200},  
AxesLabel -> {"En", "RR, \
TT"}, PlotStyle -> {{}, Dashing[{0.02, 0.02}]}];

Plot[{RR, TT} /. kRules /. {V0 -> +1} /. values // 
          Evaluate,  {En, 0, 
          5},  
AxesLabel -> {"En", "RR, \
TT"}, PlotStyle -> {{}, 
            Dashing[{0.02, 0.02}]}, 
PlotPoints -> 
          50
];

Part e

\(\(Plot3D[\(\(TT /. kRules\)\  /. values\) /. a -> 1 // 
          Evaluate, {En, 0, 10}, {V0, \(-20\), 
          0}, \ 
PlotPoints -> 
          40, 
\ 
        Ticks -> {Range[0, 10, 5], 
            Range[\(-20\), 0, 10], {0, 1/2, 
              1}\ \ }, 
AxesLabel -> {"Energy", "Potential", "TT"}];\)\)

\(\(Plot3D[\(\(\(RR /. kRules\) /. ABRTrule\) /. values\) /. a -> 1 // 
          Evaluate, {En, 0, 10}, {V0, \(-20\), 
          0}, 
PlotPoints -> 40, 
\ 
        Ticks -> {Range[0, 10, 5], 
            Range[\(-20\), 0, 10], {0, 1/2, 
              1}}, 
AxesLabel -> {"Energy", "Potential", "RR"}];\)\)

Part f

{ψL[x_, aa_], ψW[x_, aa_], ψR[x_, 
            aa_]} =  {ψL[x], ψW[x], ψR[x]}  /. 
                ABRTrule /. kRules //. {ℏ -> 1,  m -> 1, 
              a -> aa}  // Simplify;

ψ[x_, a_]:=Which[
        x  ≤  -a,  ψL[x, a], -a < x ≤  
          a,  ψW[x, a], a < x,  ψR[x, a]]

Plot[  ψ[x, 2]  /. { En -> 1,  
              V0 -> -12} // Re, {x, -15, 15}, 
        Epilog -> {Thickness[0.02], 
            Line[{{-2, 0}, {2, 0}}]}];

\(\(\(\ \)\(Plot[\ ψ[x, 2] \(ψ[x, 2]\^*\)\  //. {\ 
            En -> 1, \ V0 -> \(-12\)}\ , {x, \(-15\), 15}, 
        Epilog -> {Thickness[0.02], 
            Line[{{\(-2\), 0}, {2, 0}}]}]\ ;\)\)\)

Part g

timePlot[t_, En_, 
          V0_] := 
Plot[ {(E^(-I En t/ℏ)) ψ[x, 2 ]  // 
                  Re, (E^(-I En t/ℏ)) ψ[x, 2] // Im} //. { ℏ -> 1}  // 
            Evaluate, 
{x, -15, 15}, 
          PlotStyle -> {{}, {Dashing[{0.02}], 
                Thickness[
                  0.008]}}, 
PlotRange -> {-2, 
              2}, 
Epilog -> {Thickness[0.02], 
              Line[{{-2, 0}, {2, 0}}]}];

frames=1;
 V0=1.1;
 En=1;
 timestep=2 π/(frames En);
 Do[timePlot[n timestep ,En,V0],{n,1,frames}]

Problem 5: The One-Dimensional Harmonic Oscillator

Solution

Clear["Global`*"];

Part a

eq1=0==(V - En) ψ(x)-(((ℏ^2) ψ''(x))/(2 m)) /.V->((k (x^2))/2)

eq2 = DSolve[eq1, ψ[x], x][[1]]  // PowerExpand  // 
        Simplify

notation = {k -> (ω^2) m,  
          En   ->  ((1/2) + 
                  n) ℏ ω,  ω -> (ℏ/(m (α^2)))};

dsol  = eq2[[1]]  //. notation // PowerExpand // 
        Simplify

ψ0[x_, n_] = ψ[x] /. dsol  /. C[2] -> 0 /. 
        C[1] -> 1

energylevels=En /. notation/.{ω->(Sqrt[(k/m)])}

energylevels /. n -> Range[0, 3, 1]

Part b

norm[n_]   := 
      Integrate[ψ0[x , n] ψ0[x , 
              n] , {x, -∞, ∞}, 
          Assumptions -> { α > 0 , Element[n, Integers]} ] // 
        Sqrt

Table[
      norm[n] ==     Sqrt[  2^
              n   n!  α Sqrt[π]  ] , {n, 0, 
        5}]

ψ[x_, n_]=  ((ψ0[x, n])/((π^(1/4)) (Sqrt[(2^n) α n!])))

Table[
        Integrate[ψ[x, n]^2 , {x, -∞, ∞},  
          Assumptions -> { α > 0}], {n, 0, 3}]  // 
      PowerExpand

table = ψ[x, #] & /@ Range[0, 3] // Simplify

\(TableForm[
      table, 
TableHeadings -> {{\ "ψ[x,0]", "ψ[x,1]", "ψ[x,2]", "ψ[x,3]"\ \ \ }}]\)

Plot[  Table[ψ[x, n]/2  +   n, {n, 0, 
                7}]   /. {α -> 1}  // 
          Evaluate , 
{x, -6, 
          6} , 
Epilog -> Plot[
              0.6  z^2, {z, -6, 
                6}, 
PlotStyle -> 
                Thickness[ .01], 
DisplayFunction -> 
                Identity
][[1]]  
];

Plot[  Table[ψ[x, n]^2 +  n, {n, 0, 
                7}]   /. {α -> 1}  // Evaluate , {x, -5, 
          5}, 
Epilog -> Plot[
              0.6  z^2, {z, -6, 
                6}, 
PlotStyle -> 
                Thickness[ .01], 
DisplayFunction -> 
                Identity
][[1]] 
];

9.3 Three-Dimensional Schroedinger's Equation

Problem 1:  Three-Dimensional Harmonic Oscillator in Cartesian Coordinates

Required Packages

Needs["Calculus`VectorAnalysis`"]

Solution

Clear["Global`*"];

Part a

SetCoordinates[Cartesian[x, y, z]];

ψ[x_,y_,z_] = X[x] Y[y] Z[z];

eq1 =  ((-(ℏ^2))/(2  m)) Laplacian[ψ[x, y, 
              z]] +  (Vx[x] + Vy[y] + Vz[z]  -  En) ψ[x, y, 
            z]

eq2 = eq1/ψ[x, y, z] // ExpandAll

Enrule = {En -> Ex + Ey + Ez} ;

select[eq_,Ei_,var_]:= Select[eq,(!FreeQ[#,var]|| !FreeQ[#,Ei ])&]

fieldEq = 
0 == { 
            select[eq2 /. Enrule, Ex, x],   select[eq2 /. Enrule, Ey, y], 
            select[eq2 /. Enrule, Ez, z]} // Thread; 
    fieldEq  // ColumnForm

Part b

potential={Vx(x) -> ((kx (x^2))/2), 
            Vy(y) -> ((ky (y^2))/2), Vz(z) -> ((kz (z^2))/2)};

{xEq, yEq, zEq} = fieldEq /. potential;
{xEq, yEq, zEq} // ColumnForm

notation={kx -> m ω(x^2), 
            Ex -> (nx + 
                    (1/2)) ωx ℏ, ωx -> (ℏ/(m α(x^2)))};

xSol =  DSolve[xEq //. notation, X[x], x][[1]] // 
        Simplify

X[x_, nx_] =X[x]   /. xSol//.{(Subscript[c, 2])->0,(Subscript[c, 1])->(1/((π^(1/4)) (Sqrt[(2^nx) αx nx!])))}

Table[ 
        Integrate[X[x, nx]^2  , {x, -∞, ∞}, 
          GenerateConditions -> False] , {nx, 0, 
          5}]//PowerExpand

Y[y_, ny_] = X[y, ny]  /. {αx -> αy}

Z[z_, nz_] = X[z, nz] /. {αx -> αz}

Part c

αRule =  {αx -> 1,  αy -> 
            1};

wavePlot = 
        Plot3D[X[x, 3] Y[y, 
                2] /. αRule, 
{x, -4, 
            4}, {y, -4, 4}, 
PlotPoints -> 
            35, 
DisplayFunction -> 
            Identity];

vPlot=Plot3D[0.15 ((x^2) + (y^2))-2,
{x, -4, 4},{y, -4, 4},
PlotPoints -> 35,
DisplayFunction -> Identity];

Show[{vPlot, wavePlot}, 
        DisplayFunction -> $DisplayFunction ];

Problem 2:  Schroedinger's Equation for Spherically Symmetric Potentials

Required Packages

Needs["Calculus`VectorAnalysis`"]

Solution

Clear["Global`*"];

Part a

SetCoordinates[Spherical[r, θ, ϕ]];

ψ[r,θ,ϕ] =R[r] Θ[θ] Φ[ϕ];

eq1= -((ℏ^2)/(2 m))  Laplacian[ψ[r, θ, ϕ]]+ (V  -  En) ψ[r, θ, ϕ]

eq2=eq1 (((r^2) (sin^2)(θ))/ψ)[
                  r, θ, ϕ]//ExpandAll

select[eq_,var_]:= Select[eq,(!FreeQ[#,var])&]

eq3=select[eq2, ϕ]==(( (ℏ^2) m(ϕ^2))/(2 m))

ϕEq=Solve[eq3,Φ''[ϕ] ][[1,1]] /.Rule->Equal

eq4=DSolve[ϕEq,Φ[ϕ],ϕ][[1]]

eq4 /. {C[1] -> 1,  C[2] -> {I, -I}} // 
          MapAll[TrigToExp, #] &  // Map[Thread, #, 2] &  // 
      Flatten

ϕRule=Φ[ϕ]->(E^(I mϕ ϕ));

Part b

eq5 =  (eq2 /. (ϕEq /. Equal -> Rule ))/
          Sin[θ]^2//Expand

eq6=select[eq5, r]==-((ℓ (ℓ + 
                      1) (ℏ^2))/(2 m))

rEq= Solve[eq6,R''[r] ][[1,1]] /. Rule->Equal // ExpandAll

eq7= select[eq5, θ]== ((ℓ (ℓ + 
                    1) (ℏ^2))/(2 m))

θEq =  Solve[
              eq7, Θ''[θ]][[1, 1]]  /. 
          Rule -> Equal//ExpandAll

Part c

vChange[ Eq_, ψ_, x_, z_, f_] :=  
      Eq /. {D[ψ[x], {x, 
                n_}] :> 
Nest[( (1/
                        D[f, z])  D[#, z]) &, ψ[z], n], ψ[
              x] :> ψ[z], x :> f}

eq8 =  vChange[θEq, Θ, θ, z, 
            ArcCos[z]]  // 
          ExpandAll// Simplify

eq9 = DSolve[eq8 , Θ[z] , z] // Flatten

LegendreQ[2 , 1, z]

θRule = 
      eq9 /. {C[1] -> 1,  C[2] -> 0}  // Simplify

Part d

norm[ℓ, m]  Θ[
                z] Φ[ϕ]  /. θRule /. ϕRule /.  z -> 
          Cos[θ] /. mϕ->m

norm[ℓ_, m_]=(1/(2 (Sqrt[π]))) (Sqrt[(((2 ℓ + 1) (ℓ - m)!)/((m + ℓ)!))]);

Y[ℓ_, 
        m_, θ_, ϕ_] = norm[ℓ, m] Θ[
                  z] Φ[ϕ]  /. θRule /. ϕRule /.  z -> 
            Cos[θ] /. mϕ->m

Table[ 
          SphericalHarmonicY[ℓ, m, θ, ϕ] == 
            Y[ℓ, m, θ, ϕ], {ℓ, 0, 
            4}, {m, -ℓ, ℓ}] // 
        Simplify//PowerExpand

Problem 3: Particle in an Infinite, Spherical Well

Required Packages

Needs["NumericalMath`BesselZeros`"]
 Needs["Graphics`Legend`"]

Solution

Clear["Global`*"];

Part a

rEq=0==((2 m R(r) (En - V(r)) (r^2))/(ℏ^2))+R''(r) (r^2)+2 R'(r) r-ℓ (ℓ + 1) R(r) /.{V(r) -> 0, 
            En -> (((k^2) (ℏ^2))/(2 m))}

dsol =  DSolve[rEq , R[r], r] // Flatten // Simplify

R[r_, ℓ_, 
        k_:  1]  =  R[r]  /. 
          dsol  /. {C[1] -> (Sqrt[1 + 2  ℓ]),  
          C[2] -> 0}

Integrate[ R[r, ℓ, k]^2 , {r, 0, ∞} , 
      Assumptions -> { k > 0, ℓ > -1/2}]

\(\(TableForm[\ \(R[r, #, k] &\) /@ 
            Range[0, 
              3], 
TableHeadings -> {\ {\ "ℓ=0", "ℓ=1", "ℓ=2", "ℓ=3"\ \ }}\ ]\  // Simplify\) // PowerExpand\)

Part b

?BesselJZeros

zeros = BesselJZeros[1/2 + #, 5] &  /@   Range[0, 
          4]

\(TableForm[zeros, \ 
      TableSpacing -> \ 
        1\ , \ \ \ TableHeadings -> {{"S-Wave", "P-Wave", "D-Wave", "F-Wave", "G-Wave"}, 
{"First", "Second", "Third", "Fourth", "Fifth"}}]\)

levels = 
        Table[
{{i - 0.4, zeros[[i, j]]}, {i + 0.4, 
              zeros[[i, j]]}}, 
{i, 1, 
            Length[zeros]}, {j, 1, Length[zeros[[1]]]}] ;

Show[
{Graphics[{Thickness[0.01], 
              Map[Line, levels, {2}]}], 
Graphics[
            Text["S", {1, 1}]], Graphics[Text["P", {2, 1}]], 
          Graphics[Text["D", {3, 1}]], 
Graphics[
            Text["F", {4, 1}]], 
          Graphics[
            Text["G", {5, 
                1}]]
}, 
Frame -> 
          True, 
PlotRange -> {{0, 5.9}, {-1, 
              25}}];

Part c

\(\(\(\ \)\(Plot[
        R[r, 0, k] /. {k -> \((zeros[\([0 + 
                        1]\)])\)}
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  // Evaluate, 
{r, \ 0\ , 
          1}, \ \ PlotStyle -> {{}, Dashing[{0.01}], Dashing[{0.02}], 
            Dashing[{ .03}], 
            Dashing[{ .05}]}\ , 
PlotLabel -> 
          StyleForm["ψ", 
            FontSize -> 
              14], 
PlotLegend -> \ \ {k1, k2, k3, k4, 
            k5}\ \ , 
LegendPosition -> {1, \(-1\)/
              2}, 
Epilog -> 
          Text[StyleForm["S(ℓ=0)-State", 
              FontSize -> 14], {0.6, 
              2}]\ 
];\)\)\)

rPlots = 
Table[
Plot[ 
            R[r, ℓ, zeros[[+1 + ℓ, n]]] // 
              Evaluate, {r, 0.01, 1},  Ticks -> None,  
            PlotRange -> 
              All, 
PlotLabel -> {n, ℓ}, 
DisplayFunction -> 
              Identity
],  {n, 1, 5}, {ℓ, 0, 
            4}];

rPlots // GraphicsArray  // Show;

Problem 4: Particle in a Finite, Spherical Well

Required Packages

Needs["Graphics`Legend`"]

Solution

Clear["Global`*"];

Part a

rEq=0==(k^2) R(r) (r^2)+R''(r) (r^2)+2 R'(r) r-ℓ (ℓ + 1) R(r) /.En -> -Wn/.V -> -V0

rRule =  DSolve[rEq , R[r], r]  // Flatten // 
        ExpandAll

kRule={kIn -> (((Sqrt[2]) (Sqrt[m]) (Sqrt[V0 - Wn]))/ℏ), 
            kOut -> ((I (Sqrt[2]) (Sqrt[m]) (Sqrt[Wn]))/ℏ)};

Series[ R[r] /. rRule /. ℓ -> # , {r, 0, 
              1}] &   /@   Range[0, 3] // ColumnForm

waveIn[kIn_, ℓ_, 
        r_] = R[r] /. rRule /. {C[1] -> cIn, C[2] -> 0, 
              k -> kIn} // Simplify // PowerExpand

\(\(TableForm[\(waveIn[kIn\ , #, r\ ] &\) /@ 
            Range[0, 
              3], 
TableHeadings -> {\ {\ "ℓ=0", "ℓ=1", "ℓ=2", "ℓ=3"\ \ }}\ ]\  // Simplify\) // PowerExpand\)

eq1 = R[r] /. rRule /. k  -> I qOut // 
        PowerExpand

eq1 /. ℓ -> 0 // TrigToExp // Simplify // 
        PowerExpand//Collect[#, { 
            Exp[a_]  }] &

eq2 = eq1 /. C[2] -> 1 C[1]

waveOut[qOut_,ℓ_,r_]= eq1 /.rRule /.{C[1]-> cOut,C[2]-> +I cOut}

\(\(\(TableForm[\(\((\(\(waveOut[qOut\ , #, r\ ] // TrigToExp\)\  // 
                Simplify\) // PowerExpand\ \ )\) &\) /@ 
        Range[0, 
          3], 
TableHeadings -> {\ {\ "ℓ=0", "ℓ=1", "ℓ=2", "ℓ=3"\ \ }}\ ]\)\(\ \ \)\)\)

Part b

bc1 = (waveIn[kIn, ℓ, r] == 
            waveOut[qOut, ℓ, r])//.{r -> 
            a}

cOutRule =  Solve[bc1, cOut] /. cIn -> 1  // Flatten // 
        Simplify

bc2= (D[waveIn[kIn,ℓ,r],r] == D[waveOut[qOut,ℓ,r],r]) //.{r->a}//Simplify;

allRules = {kIn -> (((Sqrt[2]) (Sqrt[m (V0 - Wn)]))/ℏ), 
            qOut -> (Sqrt[2]) (Sqrt[((m Wn)/(ℏ^2))]),   cOutRule,  ℏ -> 1, a -> 1, m -> 1, V0 -> 110, cIn -> 1} // 
          Flatten ;

bc2Sol[ℓ_] =  bc2 //. cOutRule /. cIn -> 1 //. 
              allRules  //. allRules  // Simplify;

eq3 = bc2Sol[0] // Simplify

Plot[eq3[[1]]  // Evaluate, {Wn, 0, 110}];

sValues =  FindRoot[ eq3 , {Wn, #}] &   /@   {10, 40, 70, 
          90}

eq4  =  bc2Sol[1]  // Simplify  ;

eq5 =  bc2Sol[2]  // Simplify ;

Plot[{eq4[[1]], eq5[[1]]} // 
          Evaluate, 
{Wn, 0, 
          110}, 
PlotStyle -> {{}, 
            Dashing[{0.01}]} ];

pValues =   FindRoot[ eq4 , {Wn, #}] &   /@   {20, 60, 80, 
          100}

dValues =   FindRoot[ eq5, {Wn, #}] &   /@   {10, 40, 70, 
          90}

\(TableForm[{sValues, pValues, dValues}, \ 
      TableSpacing -> \ 
        1\ , \ \ \ TableHeadings -> {{"S-Wave", "P-Wave", "D-Wave"\ \ }, {"First", "Second", "Third", "Fourth"\ }}]\)

Part c

{in[r_, ℓ_, Wn_], 
          out[r_, ℓ_, Wn_]} = {waveIn[kIn, ℓ, r], 
            waveOut[qOut, ℓ, r]} //. allRules;

wave[r_, ℓ_, Wn_] = 
        Which[r ≤ 1 , in[r, ℓ, Wn], r > 1 , 
          out[r, ℓ, Wn]];

Plot[ wave[r, 0, Wn]  /. sValues  // Re // 
          Evaluate , 
{r, 0, 1.3},  
        PlotStyle -> {{}, Dashing[{0.01}], Dashing[{0.02}], 
            Dashing[{ .03}]},  
PlotLegend ->  
          Evaluate[ Wn /. sValues[[#]] & /@ 
              Range[4]] , 
LegendPosition -> {1, -1/2}, 
Epilog -> 
          Line[{{1, -3}, {1, 3}}]] ;

Plot[ wave[r, 1, Wn]  /. pValues // Re // 
          Evaluate    , 
{r, 0, 1.3},  
        PlotStyle -> {{}, Dashing[{0.01}], Dashing[{0.02}], 
            Dashing[{ .03}]}, 
PlotLegend ->  
          Evaluate[ Wn /. pValues[[#]] & /@ 
              Range[4]] , 
LegendPosition -> {1, -1/2}, 
Epilog -> 
          Line[{{1, -3}, {1, 3}}]
];

Plot[ wave[r, 2, Wn]  /. dValues // Re // 
          Evaluate    , 
{r, 0, 1.3},  
        PlotStyle -> {{}, Dashing[{0.01}], Dashing[{0.02}], 
            Dashing[{0.03}]}, 
PlotLegend ->  
          Evaluate[ Wn /. dValues[[#]] & /@ 
              Range[3]] , 
LegendPosition -> {1, -1/2}, 
Epilog -> 
          Line[{{1, -3}, {1, 3}}]
];

Problem 5: The Hydrogen Atom in Spherical Coordinates

Required Packages

Needs["Graphics`ParametricPlot3D`"]

Solution

Clear["Global`*"];

Part a

eq1=0==((2 m R(
                        r) (En - V(r)) (r^2))/(ℏ^2))+R''(r) (r^2)+2 R'(r) r-ℓ (ℓ + 1) R(r)/.V[r]->-((e^2)/(4 π r ϵ0))/.En -> -Wn

eq2=DSolve[eq1,{R[r]},r][[1]]

substitutions={m -> ((4 π ϵ0 (ℏ^2))/(a0 (e^2))), Wn -> (((e^4) m)/(32 (n^2) (π^2) ϵ(0^2) (ℏ^2)))};

eq3 = eq2 //. substitutions   // PowerExpand // 
        Simplify

eq4 = eq3 /. {C[1] -> 0 }

rNorm=(2^(ℓ + 1)) ((((1/(a0 n))))^ℓ) (Sqrt[(((n - ℓ - 1)!)/((a0^3) (n^4) (n + ℓ)!))]);

rwave[n_, ℓ_, 
        r_]  = R[r] /. eq4[[1]] /. C[2] -> rNorm // 
          PowerExpand // Simplify

Integrate[ (r^2) {  rwave[1, 0, r ], rwave[2, 0, r ], rwave[2, 1, r ]}^2, {r, 0, ∞}, 
      Assumptions -> {Re[a0] > 0}]

Part b

eq5 = Table[
          rwave[n, ℓ, r],  {n, 1, 3}, {ℓ, 0, 
            n - 1} ]// TableForm[#,  
            TableHeadings -> {{"n=1", "n=2", "n=3"}, {"\
ℓ=0", "ℓ=1", "ℓ=2"}}] &

plot[n_, ℓ_, rmax_] := 
Plot[
          (r^2) rwave[n, ℓ, r]^2 /. a0 -> 1 // 
            Evaluate, 
{r, 0.01, 
            rmax}, 
Ticks -> {{1, 4, 8, 14, 20}, {0.1, 
                0.3, 0.5}}, 
DisplayFunction -> 
            Identity];

Show[GraphicsArray[ Table[plot[n,ℓ,n*7],{n,1,3},{ℓ,0,n-1}]]];

Part c

plot[ℓ_, 
          m_] := 
SphericalPlot3D[
SphericalHarmonicY[ℓ, m, θ, ϕ]  // Abs // 
            Evaluate, 
{θ, 0, π}, {ϕ, 0, 
            2 π}, 
Boxed -> 
            False, 
Axes -> 
            False, 
DisplayFunction -> 
            Identity];

pt1=Table[plot[ℓ,m] ,{ℓ,0,2} ,{m,-ℓ,ℓ}];
Show[GraphicsArray[pt1]];

9.4 Exercises

Exercise 1:  Infinite Potential Well with Rectangular Perturbation

Exercise 2:  Tilted Square Well

Exercise 3: The Wentzel-Kramers-Brillouin Approximation

Exercise 4: Plot of the Electron Probability Density

Exercise 5: Perturbed Harmonic Oscillator

Exercise 6: Perturbation Theory

Exercise 7:  Separation of variables in Cylindrical Coordinates