HW #1 Prob 1: Stern Gerlach deflection

Remove["Global`*"]

The following cell is simply a collection of the (SI) values for many commonly used quantities

ℏ=QuantityMagnitude[UnitConvert[Quantity["ReducedPlanckConstant"]]];
c=QuantityMagnitude[UnitConvert[Quantity["SpeedOfLight"]]];
e=QuantityMagnitude[UnitConvert[Quantity["ElementaryCharge"]]];
k=QuantityMagnitude[UnitConvert[Quantity["BoltzmannConstant"]]];
u=QuantityMagnitude[UnitConvert[Quantity["AtomicMassUnit"]]];
μB=QuantityMagnitude[UnitConvert[Quantity["BohrMagneton"]]];

Now we work on the specific problem and get some specific numbers

mAg=u QuantityMagnitude[ElementData["Silver","AtomicWeight"]];

dBdz=10;
temp=1227+273
length=1;
dist=1;

Find the time the sliver atom spends in the inhomogeneous magnetic field

energy=(3/2) k temp;
v=Sqrt[2 energy/mAg];
time=length/v;

Find the z- velocity from the force

force=μB dBdz;
accel=force/mAg;
vz=accel time;

Finally, find the deflection angle and the separation between the two beams (in mm) at the screen

θ=vz/v;
sep=1000 2 θ dist

HW #6 Prob 1: Plotting the amplitude of Rabi’s formula

Remove["Global`*"]

Write down the amplitude

amp=(ω1^2/4)/((ω0-ω)^2+(ω1^2/4));

Substitute for scaled version of the variable frequency

ampx=amp/.ω->x ω0;

Substitute for scaled version of the variable field ratio

ampr=Simplify[ampx/.ω1->r ω0]

Make the plot

Plot[{ampr/.r->0.01,ampr/.r->0.001,ampr/.r->0.0001},{x,0.975,1.025},
PlotRange->{0,1},PlotStyle->{Black,Red,Blue},
PlotLegends->{0.01,0.001,0.0001}]

HW #7 Prob 2: Linear combination of square well states

Remove["Global`*"]
$Assumptions={x∈Reals,t∈Reals,m>0,ℏ>0,L>0};

Set up the energies and eigenfunctions

En=ℏ^2Pi^2n^2/(2m L^2);
E1=En/.n->1;
E2=En/.n->2;

ψn=Sqrt[2/L]Sin[Pi x n/L];
ψ1=ψn/.n->1;
ψ2=ψn/.n->2;

Write the wave function at t=0

a1=(1+I)/2;
a2=1/Sqrt[2];
ψt0=a1 ψ1+a2 ψ2;
ψt0Star=Conjugate[a1] ψ1+a2 ψ2;

Check that it is properly normalized

Integrate[ψt0Star ψt0,{x,0,L}]

Put in the time dependence

ψt=a1 Exp[-I E1 t/ℏ]ψ1+a2 Exp[-I E2 t/ℏ] ψ2

ψtStar=Conjugate[a1] Exp[I E1 t/ℏ]ψ1+a2 Exp[I E2 t/ℏ]ψ2;

Now get the expectation value of the energy

Integrate[ψtStar(-ℏ^2/(2m))D[ψt,{x,2}],{x,0,L}]

Calculating the expectation value of position is straightforward

Integrate[ψtStar x ψt,{x,0,L}]

TeXForm[%]

The probability of measuring E1 is the square of the projection onto ψ1

Refine[Abs[Integrate[ψ1 ψt,{x,0,L}]]]^2

HW #8 Prob 1: Uncertainty principle for SHO

Remove["Global`*"]
$Assumptions={x∈Reals,t∈Reals,m>0,ℏ>0,ω>0};

Calculate the ground state quantities

gswf=(m ω/(Pi ℏ))^(1/4)Exp[-m ω x^2/(2ℏ)];
Integrate[gswf x gswf,{x,-∞,∞}]
Integrate[gswf x^2 gswf,{x,-∞,∞}]
(ℏ/I)Integrate[gswf D[gswf,x],{x,-∞,∞}]
-ℏ^2Integrate[gswf D[gswf,{x,2}],{x,-∞,∞}]

HW #8 Prob 2: Plot of time dependent oscillator

Remove["Global`*"]
$Assumptions={x∈Reals,t∈Reals,m>0,ℏ>0,ω>0};

Write down the n=0 and n=1 wave functions

u0=(m ω/(Pi ℏ))^(1/4)Exp[-m ω x^2/(2ℏ)];
u1=((4/Pi)(m ω/ ℏ)^3)^(1/4)x Exp[-m ω x^2/(2ℏ)];

Write down the probability density function

ρ=(1/2)u0^2+(1/2)u1^2+u0 u1 Cos[ω t];

Check that it is normalized

Integrate[ρ,{x,-∞,∞}]

Redefine with time scaled by the period T=2π/ω

ρScaled=Simplify[(ρ/.{t->2 Pi τ/ω,x->y/Sqrt[m ω/ℏ]})/Sqrt[m ω/ℏ]];
Integrate[ρScaled,{y,-∞,∞}]

Plot as in Figure 7.10, for t=0, T/4, and T/2 (i.e. τ=0, 1/4, and 1/2)

Plot[
{ρScaled/.τ->0,
ρScaled/.τ->1/4,
ρScaled/.τ->1/2},{y,-3,3},
PlotStyle->{Black,Blue,Red},
PlotLegends->{t=0,t=T/4,t=T/2}]

Now make the animation

Animate[
Plot[ρScaled/.τ->time,{y,-3,3},PlotRange->{0,0.7}],
{time,0,3},AnimationRunning->False]

HW #10 Prob 1: Deuteron well depth

Remove["Global`*"]

Write down some values needed for the calculation and plot. Units are MeV and fm

mN=940;
a=1.7;
hbarc=200;
eBind=2.2;

Solving the transcendental equation. Make a plot, then solve.

eVal=mN a^2 eBind/hbarc^2

Plot[
{Tan[Sqrt[eVal(v-1)]],-Sqrt[v-1]},
{v,2,100}]

vSol=v/.FindRoot[Tan[Sqrt[eVal(v-1)]]==-Sqrt[v-1],{v,20}]

V0=eBind vSol

To get the wave function, define k0 and q, then solve for the coefficient ratio and  normalize

k0=Sqrt[eVal(vSol-1)]/a;
q=Sqrt[eVal]/a;

rCoef=Exp[q a]Sin[k0 a];
normInt=NIntegrate[Sin[k0 r]^2,{r,0,a}]+rCoef^2 NIntegrate[Exp[-q r]^2,{r,a,∞}];
CoefA=1/Sqrt[normInt];
CoefC=rCoef CoefA;

Now make the plot

rMax=5a;

wfPlot=Show[
Plot[CoefA Sin[k0 r],{r,0,a},
PlotRange->{{0,rMax},{0,0.6}},
AxesLabel->{"r (fm)"}],
Plot[CoefC Exp[-q r],{r,a,rMax}],
ListPlot[{{a,0},{a,1}},Joined->True,PlotStyle->{Red,Dashed}]
]

nrgPlot=Show[
ListPlot[{{0,-V0},{a,-V0},{a,0},{rMax,0}},
PlotRange->{{0,rMax},{-1.05V0,0}},
AxesLabel->{"r (fm)"},
Joined->True,PlotStyle->{Black,Solid}],
ListPlot[{{0,-eBind},{rMax,-eBind}},
Joined->True,PlotStyle->{Green,Solid}]
]

GraphicsGrid[{{wfPlot},{nrgPlot}}]

HW #11 Prob 1: Deuteron well depth, again

Remove["Global`*"]

Start the same way as HW#10 Prob 1, and get the second solution. This sets the value of V0.

mN=940;
μ=mN/2;
a=1.7;
hbarc=200;
eBind=2.2;

eVal=2μ a^2 eBind/hbarc^2;

Plot[
{Tan[Sqrt[eVal(v-1)]],-Sqrt[v-1]},
{v,2,180}]

vSol=v/.FindRoot[Tan[Sqrt[eVal(v-1)]]==-Sqrt[v-1],{v,{20,150}}]

v1=vSol[[1]];
v2=vSol[[2]];

V0=v2 eBind

Now solve the new transcendental equation for the energy eigenvalues. Plot it first.

eVal=2μ a^2 V0/hbarc^2;
Plot[
{Tan[Sqrt[eVal(1-1/v)]],-Sqrt[v-1]},
{v,0,180}]

Plot[
{Tan[Sqrt[eVal(1-1/v)]],-Sqrt[v-1]},
{v,0,2}]

Calculate the energy eigenvalues

vSol=v/.FindRoot[Tan[Sqrt[eVal(1-1/v)]]==-Sqrt[v-1],{v,{1.5,150}}];
e1=V0/vSol[[1]]
e2=V0/vSol[[2]]

Form the wave functions and plot them

k0=Sqrt[2μ(V0-e)/hbarc^2];
q=Sqrt[2μ e/hbarc^2];
rCoef=Exp[q a]Sin[k0 a];

rCoef1=rCoef/.e->e1;
wf1A=Sin[k0 r]/.e->e1;
wf1B=Exp[-q r]/.e->e1;
normInt=NIntegrate[wf1A^2,{r,0,a}]+rCoef1^2 NIntegrate[wf1B^2,{r,a,∞}];
Coef1A=1/Sqrt[normInt];
Coef1B=rCoef1 Coef1A;

rCoef2=rCoef/.e->e2;
wf2A=Sin[k0 r]/.e->e2;
wf2B=Exp[-q r]/.e->e2;
normInt=NIntegrate[wf2A^2,{r,0,a}]+rCoef2^2 NIntegrate[wf2B^2,{r,a,∞}];
Coef2A=1/Sqrt[normInt];
Coef2B=rCoef2 Coef2A;

rMax=5a;
Show[
Plot[Coef1B wf1B,{r,a,rMax},PlotStyle->Black,
PlotRange->{{0,rMax},{-0.6,1}},AxesLabel->{"r (fm)"}],
Plot[Coef1A wf1A,{r,0,a},PlotStyle->Black],
Plot[Coef2A wf2A,{r,0,a},PlotStyle->Red],
Plot[Coef2B wf2B,{r,a,rMax},PlotStyle->Red],
ListPlot[{{a,-0.6},{a,1}},Joined->True,PlotStyle->{Blue,Dashed}]
]

HW #12 Prob 2: H-atom perturbation from finite proton size

Remove["Global`*"]

$Assumptions={a0>0,e>0,R>0};

Write down the wave functions from (10.43) and (10.44)

Z=1

R10=2(Z/a0)^{3/2}Exp[-Z r/a0];
R20=2(Z/(2a0))^{3/2}(1-Z r/(2a0))Exp[- r/(2a0)];
R21=(1/Sqrt[3])(Z/(2a0))^{3/2}(Z r/a0)Exp[- r/(2a0)];

Write down the perturbation

pert=e^2/r-(3/2)(e^2/R^3)(R^2-r^2/3);

Do the integrals

Δ10=Integrate[r^2 R10 pert R10,{r,0,R}];
Δ20=Integrate[r^2 R20 pert R20,{r,0,R}];
Δ21=Integrate[r^2 R21 pert R21,{r,0,R}];

Expand the results in powers of R/a0

Δx=Normal[Series[Δ10/.R->x a0,{x,0,5}]];
Δx[[1]]/.x->R/a0

Δx=Normal[Series[Δ20/.R->x a0,{x,0,5}]];
Δx[[1]]/.x->R/a0

Δx=Normal[Series[Δ21/.R->x a0,{x,0,5}]];
Δx[[1]]/.x->R/a0

HW #14 Prob 2: Variational Principle for 1D linear potential

Remove["Global`*"]

$Assumptions={a>0,m>0,g>0,h>0,ℏ>0};

Write down the two trial wave functions

ψ1=x Exp[-x/(2a)];
ψ2=x Exp[-x^2/(2a^2)];
ψs={ψ1,ψ2};

Calculate the expectation values of the Hamiltonian

expH=Integrate[
ψs(-ℏ^2/(2m)D[ψs,{x,2}]+λ x ψs),
{x,0,∞}]/
Integrate[ψs^2,{x,0,∞}]

Minimize these with respect to a

expHderv=D[expH,a];

a1=a/.Solve[expHderv[[1]]==0,a][[2]]

a2=a/.Solve[expHderv[[2]]==0,a][[2]]

Extract the ground state energy approximations and compare them

Eg1=expH[[1]]/.a->a1

Eg2=expH[[2]]/.a->a2

Simplify[Eg2/Eg1]

N[%]

As expected, the second trial function gives a smaller (and therefore better) answer

For a “bouncing ball” λ=mg, where x is the height to which the ball bounces. Get the height for the ground state.

h=Simplify[(Eg2/.λ->m g)/( m g)]

Calculate the height (in μm) for a neutron bouncing on the Earth’s surface

g=9.8;
ℏ=QuantityMagnitude[UnitConvert[Quantity["ReducedPlanckConstant"]]];
m=QuantityMagnitude[UnitConvert[Quantity["NeutronMass"]]];

10^6N[h]