130B Quantum Physics

Part 1. Path Integral Quantization

From Classical to Quantum

Historical Review

History: What is the Nature of Light?

The corpuscular (particle) theory: light is composed of steady stream of particles carrying the energy and travelling along rays in the speed of light.

The wave theory: light is wave-like, propagating in the space and time.

The Wave-Particle Wars in History

Concluding Remarks

Light: originally considered as wave, also behaves like particle.

Electrons, α particles (( ^4)He nucleus): originally considered as particles, also behave like wave.

Quantization of Light

Geometric Optics

Fermat’s Principle: Light always travels along the path of extremal optical path length.

The optical path length is defined by

where n is the refractive index of the medium and ds is an infinitesimal displacement along the ray.

The optical path length is simply related to the light travelling time T by L=c T, where c is the speed of light in vacuum. So extremization of either of them will be equivalent.

Eikonal equation (Newton’s law of light)

n(d/(dt))((n^2)((dx)/(dt)))=(c^2)∇n.

Derive Eq. () from Fermat’s principle Eq. ().

L=(∫_A^B)n(x)ds
=∫n(x)|(OverDot[x])|dt.

n(x) is the refractive index at position x,

ds=|(OverDot[x])|dt is the infinitesimal arc length along the path.

ℒ(x,(OverDot[x]))=n(x)|(OverDot[x])|

(d/(dt))(((∂L)/(∂(OverDot[x]))))-((∂L)/(∂x))=0.

((∂L)/(∂(OverDot[x])))=n(x)((OverDot[x])/(|(OverDot[x])|)),
((∂L)/(∂x))=∇n(x)|(OverDot[x])|,

(d/(dt))(n((OverDot[x])/(|(OverDot[x])|)))=∇n|(OverDot[x])|.

n(d/(dt))((n^2)((OverDot[x])/(n|(OverDot[x])|)))=n|(OverDot[x])|∇n
=>n(d/(dt))((n^2)((OverDot[x])/c))=c∇n
=>n(d/(dt))((n^2)((dx)/(dt)))=(c^2)∇n.

Examples of light rays in the medium by solving Eq. ():

Refraction (Snell’s law)

Module[{x0=-1,y0=-1,θ=30Degree,n,tend=5,sol,data,xmin,xmax,ymin,ymax,colbg=ColorData["Legacy"]["DodgerBlue"],colray=ColorData["Legacy"]["Firebrick"]},n[x_,y_]:=1.366+0.366Tanh[25y];sol=NDSolve[{n[x[t],y[t]]D[n[x[t],y[t]]^2D[x[t],t],t]==D[n[x[t],y[t]],x[t]],n[x[t],y[t]]D[n[x[t],y[t]]^2D[y[t],t],t]==D[n[x[t],y[t]],y[t]],x[0]==x0,y[0]==y0,x'[0]==Cos[θ]/n[x0,y0],y'[0]==Sin[θ]/n[x0,y0]},{x,y},{t,0,tend}];data=First@Cases[ParametricPlot[Evaluate[{x[t],y[t]}/.sol],{t,0,tend},PlotPoints->10,MaxRecursion->3],Line[pts_]:>pts,∞];{{xmin,xmax},{ymin,ymax}}={Min[#],Max[#]}&/@Thread[data];Show[ContourPlot[n[x,y],{x,xmin,xmax},{y,ymin,ymax},Contours->16,ContourStyle->None,ColorFunction->(Blend[{White,colbg},Rescale[#,{0.9,2.0}]]&),ColorFunctionScaling->False],Graphics[{AbsoluteThickness[1.2],colray,Arrowheads[{{0.001,0.2,Arrowhead[]},{0.001,0.8,Arrowhead[]}}],Arrow[data]}],AspectRatio->Automatic,PlotRangePadding->0.05,FrameTicks->None,ImageSize->120]]

Total reflection

Module[{x0=-1,y0=1,θ=-45Degree,n,tend=6.6,sol,data,xmin,xmax,ymin,ymax,colbg=ColorData["Legacy"]["DodgerBlue"],colray=ColorData["Legacy"]["Firebrick"]},n[x_,y_]:=1.3+0.3Tanh[3y];sol=NDSolve[{n[x[t],y[t]]D[n[x[t],y[t]]^2D[x[t],t],t]==D[n[x[t],y[t]],x[t]],n[x[t],y[t]]D[n[x[t],y[t]]^2D[y[t],t],t]==D[n[x[t],y[t]],y[t]],x[0]==x0,y[0]==y0,x'[0]==Cos[θ]/n[x0,y0],y'[0]==Sin[θ]/n[x0,y0]},{x,y},{t,0,tend}];data=First@Cases[ParametricPlot[Evaluate[{x[t],y[t]}/.sol],{t,0,tend},PlotPoints->10,MaxRecursion->3],Line[pts_]:>pts,∞];{{xmin,xmax},{ymin,ymax}}={Min[#],Max[#]}&/@Thread[data];Show[ContourPlot[n[x,y],{x,xmin,xmax},{y,ymin-0.5,ymax},Contours->16,ContourStyle->None,ColorFunction->(Blend[{White,colbg},Rescale[#,{0.9,2.0}]]&),ColorFunctionScaling->False],Graphics[{AbsoluteThickness[1.2],colray,Arrowheads[{{0.001,0.2,Arrowhead[]},{0.001,0.8,Arrowhead[]}}],Arrow[data]}],AspectRatio->Automatic,PlotRangePadding->0.05,FrameTicks->None,ImageSize->140]]

Gradient-index (GRIN) optics

Module[{x0=-5,θ=0Degree,n,tend=13,sols,datas,colbg=ColorData["Legacy"]["DodgerBlue"],colray=ColorData["Legacy"]["Firebrick"]},n[x_,y_]:=1+0.8Exp[-((x/2)^10+y^2/14)];sols=NDSolve[{n[x[t],y[t]]D[n[x[t],y[t]]^2D[x[t],t],t]==D[n[x[t],y[t]],x[t]],n[x[t],y[t]]D[n[x[t],y[t]]^2D[y[t],t],t]==D[n[x[t],y[t]],y[t]],x[0]==x0,y[0]==#,x'[0]==Cos[θ]/n[x0,#],y'[0]==Sin[θ]/n[x0,#]},{x,y},{t,0,tend}]&/@{-2,-1,0,1,2};datas=Cases[ParametricPlot[Evaluate[{x[t],y[t]}/.sols],{t,0,tend},PlotPoints->10,MaxRecursion->3],Line[pts_]:>pts,∞];
Show[ContourPlot[n[x,y],{x,-5,5},{y,-5,5},Contours->16,ContourStyle->None,PlotRange->All,ColorFunction->(Blend[{White,colbg},Rescale[#,{0.9,2.0}]]&),ColorFunctionScaling->False],Graphics[{AbsoluteThickness[1.2],colray,Arrowheads[{{0.001,0.2,Arrowhead[]},{0.001,0.6,Arrowhead[]}}],Arrow/@datas}],AspectRatio->Automatic,PlotRangePadding->0.2,FrameTicks->None,ImageSize->120]]

Physical Optics

Huygens’ Principle: Every point on the wavefront acts as a secondary source emitting spherical wavelet. The new wavefront is formed by the coherent superposition of these wavelets.

Block[{y0=-10,r0=10,dr=0.1,dθ=0.005π,dθ2=0.05π,θmax=0.2π,θb=0.16π,wv,wc},
wv[r_,θ_,dk_]:=Block[{col},col=Lighter[Yellow,1-dk];{EdgeForm[col],col,Polygon[ReIm[I Flatten[-y0+{r  Exp[I (θ+dθ{-1,1}/2)],(r+dr)Exp[I (θ+dθ{1,-1}/2)]}]]]}];
wc[r_,θ_,dk_]:=Table[{Opacity[Mod[a,1] dk],Yellow,Annulus[ReIm[I(-y0+r Exp[I θ])],{a-dr,a}]},{a,2-dr,2dr,-dr}];
Graphics[{Table[wv[r,θ,Mod[r,1]LogisticSigmoid[20(θb^2-θ^2)]LogisticSigmoid[0.5r Sin[θ]]],{θ,-θmax,θmax,dθ},{r,r0+3,r0+5-dr,dr}],
Block[{r=r0+3},Table[wc[r,θ,LogisticSigmoid[-0.5r Sin[θ]]],{θ,dθ2,-θmax+dθ2,-dθ2}]],
Table[wv[r,θ,Mod[r,1]LogisticSigmoid[20(θb^2-θ^2)]],{θ,-θmax,θmax,dθ},{r,r0,r0+3-dr,dr}],
Translate[{Point[{0,0}],Text["\!\( A\)",{0,0},{0,1.2}]},ReIm[I(-y0+(r0+3) Exp[-I dθ2])]],Translate[{Point[{0,0}],Text["\!\( B\)",{0,0},{0,-1}]},ReIm[I(-y0+(r0+4) Exp[I (-dθ2+0.01π)])]]},ImageSize->220]]

Wave propagation: The wave amplitude (Subscript[ψ, B]) on the new wave front is determined by the amplitude (Subscript[ψ, A]) on the preceding wave front, modified by a phase factor (e^(i Θ(A→B))) that encodes the accumulated phase Θ(A→B) during wave propagation.

Interference effect: contributions from different paths must be collected and summed up (integrated over)

The accumulated phase Θ is determined by the propagation time T × the frequency ω of light

proportional to the optical path length L (given that the light propagates with a fixed frequency).

spacetime: ψ(x,t),
space (at fixed-time):ψ(x).

The magnitude (or the absolute amplitude) |ψ| of the wave is related to the intensity of the light, or the probability density to observe a photon at a given position x,

p(x)=((|ψ(x)|)^2).

Normalization: the wavefunction is said to be normalized, if

∫((|ψ(x)|)^2)(d^D)x=∫p(x)(d^D)x=1,

a requirement for the total probability to be 1.

From Fermat to Huygens

S(A->B)=((ℏ ω)/c)L(A->B),

Particle mechanics defines the action S in the variational principle δ S=0.

Wave mechanics defines the phase Θ in the wavelet propagator (e^(i Θ)).

Path Integral Quantization

Quantization of Matter

Classical Mechanics

S[x]=∫L(x,(OverDot[x]),t)dt.

(Subscript[(Subscript[δ S[x]|, x=(OverBar[x])])=δ∫L(x,(OverDot[x]),t)dt|, x=(OverBar[x])])=0.

(d/(d t))(((∂L)/(∂(OverDot[x]))))-((∂L)/(∂x))=0,

Nature may not have sufficient computational resources to carry out the classical mechanics precisely. => Classical mechanics might actually be realized only approximately as a stochastic global optimization, which is computationally more feasible.

Quantum mechanics takes a stochastic approach to optimize the action, which is more natural than the deterministic approach of classical mechanics, if we assume only limited computational resource is available to nature.

Optimization by Interference

f(x)=-(x^2)+2(x^4).

Plot[-x^2+2x^4,{x,-1,1},PlotStyle->Black,FrameLabel->{"\!\( x\)","\!\(f(x)\)"}]

Every point x is a legitimate guess of the solution.

Each point x is associated with an action f(x) (the objective function).

Raise the action f(x) to the exponent (as a phase): (e^(i f(x)/ℏ)) => call it a “probability amplitude” contributed by the point x.

A “Planck constant” ℏ=h/(2π) is introduced as a hyperparameter of the algorithm, to control “how quantum” the algorithm will be.

Contributions from all points must be collected and summed (integrated) up,

Z=(∫_(-∞)^∞)(e^(i f(x)/ℏ))dx.

dz=(e^(i f(x)/ℏ))dx.

dx controls the infinitesimal step size,

(e^(i f(x)/ℏ)) controls the direction to make the displacement,

displacement dz is accumulated to form the partition function,

Z≡∫dz.

For small h (classical limit)

Block[{f,h=0.001,dx,acc,Z1,Z2},
f[x_]:=-x^2+2x^4;
dx=0.01Sqrt[h];
acc=ReIm@Prepend[Accumulate[Table[Exp[2π I f[x]/h]dx,{x,dx/2,1/2+2Sqrt[2h],dx}]],0];
{Z1,Z2}={NIntegrate[Exp[2π I f[x]/h],{x,0,1/(2Sqrt[3])}],NIntegrate[Exp[2π I f[x]/h],{x,1/(2Sqrt[3]),∞}]};
Graphics[{{AbsoluteThickness[1],Line/@{acc,-acc}},{Arrowheads[0.08],
{Red,Arrow/@ReIm/@{{-Z1-Z2,-Z1},{-Z1,Z1},{Z1,Z1+Z2}},Text["\!\(Z\_\("<>#1<>"\)\)",ReIm@#2,#3]&@@@{{"-1/2",-Z1-Z2/2,{1,-1}},{"0",0,{-0.5,-1}},{"1/2",Z1+Z2/2,{-1,1}}}},
{Blue,Arrow/@ReIm/@{{-Z1-Z2,Z1+Z2}},
Text["\!\( Z\)",ReIm@(Z1+Z2)/2,{0,-0.6}]}}},
PlotLabel->Style["\!\(h="<>ToString[h]<>"\)",14,Black],ImageSize->160]
]

Z=(Subscript[Z, -1/2])+(Subscript[Z, 0])+(Subscript[Z, 1/2]),

Z can break up into three smaller contributions, which correspond to the contributions around the three stationary points: x=0,±1/2.

Around the stationary point, phase changes slowly (∂_x)f(x)∼0 => constructive interference => large contribution to the partition function.

The solutions of stationary points (classical solutions) emerge from interference due to their dominant contribution to the probability amplitude.

*More precisely, the partition function is actually evaluated with respect to the momentum k,

Z(k)≡∫dz (e^(i k x))≃(Subscript[Z, -1/2])(e^(-i k/2))+(Subscript[Z, 0])+(Subscript[Z, 1/2])(e^(i k/2)).

Then its Fourier spectrum (OverTilde[Z])(x)=∫dk Z(k)(e^(-i k x)) will reveal the saddle points.

Block[{ϵ=0.0001,ℏ=0.001,a=1/32,L=3,Z},
Z[k_]:=NIntegrate[Exp[(I-ϵ)(-x^2+2x^4)/ℏ]Exp[I k x],{x,-∞,∞}];
ListLinePlot[RotateLeft[#,Ceiling[Length[#]/2]]&[Abs[Fourier[Z/@Range[-π/a,π/a,2π/L]]]^2],DataRange->{-L/2,L/2},FrameLabel->{"\!\( x\)","\!\(\[LeftBracketingBar]Z\&~(x)\[RightBracketingBar]\^2\)"},PlotStyle->Black,PlotRange->All,AspectRatio->1/2,ImageSize->200]]

For intermediate h

Block[{f,h=0.1,dx,acc,Z1,Z2},
f[x_]:=-x^2+2x^4;
dx=0.01Sqrt[h];
acc=ReIm@Prepend[Accumulate[Table[Exp[2π I f[x]/h]dx,{x,dx/2,1/2+Sqrt[2h],dx}]],0];
{Z1,Z2}={NIntegrate[Exp[2π I f[x]/h],{x,0,1/(2Sqrt[3])}],NIntegrate[Exp[2π I f[x]/h],{x,1/(2Sqrt[3]),∞}]};
Graphics[{{AbsoluteThickness[1],Line/@{acc,-acc}},{Arrowheads[0.08],
{Red,Dashed,Arrow/@ReIm/@{{-Z1-Z2,-Z1},{-Z1,Z1},{Z1,Z1+Z2}},Text["\!\(Z\_\("<>#1<>"\)\)",ReIm@#2,#3]&@@@{{"-1/2",-Z1-Z2/2,{-1,-1}},{"0",0,{-1,-1}},{"1/2",Z1+Z2/2,{1,1}}}},
{Blue,Arrow/@ReIm/@{{-Z1-Z2,Z1+Z2}},Text["\!\( Z\)",0.9ReIm@(Z1+Z2),{0,-1}]}}},
PlotLabel->Style["\!\(h="<>ToString[h]<>"\)",14,Black],ImageSize->140]
]

The decomposition of Z into three subdominant amplitudes is not very well defined. => Quantum fluctuations start to smear out nearby stationary points.

For large h (quantum limit)

Block[{f,h=10,dx,acc,Z1,Z2},
f[x_]:=-x^2+2x^4;
dx=0.01Sqrt[h];
acc=ReIm@Prepend[Accumulate[Table[Exp[2π I f[x]/h]dx,{x,dx/2,1/2+Sqrt[h/2],dx}]],0];
{Z1,Z2}={NIntegrate[Exp[2π I f[x]/h],{x,0,1/(2Sqrt[3])}],NIntegrate[Exp[2π I f[x]/h],{x,1/(2Sqrt[3]),∞}]};
Graphics[{{AbsoluteThickness[1],Line/@{acc,-acc}},{Arrowheads[0.08],
{Blue,Arrow/@ReIm/@{{-Z1-Z2,Z1+Z2}},Text["\!\( Z\)",0.9ReIm@(Z1+Z2),{0,-1}]}}},
PlotLabel->Style["\!\(h="<>ToString[h]<>"\)",14,Black],ImageSize->160]
]

Stationary points are indistinguishable if quantum fluctuations are too large. => As if there is only one (approximate) stationary point around x=0.

If there is no sufficient resolution power, fine structures in the action landscape will be ignored by quantum mechanics. In this way, the computational complexity is controlled.

h=6.62607004×(10^(-34))J s.

h is very small (in our everyday unit) => our nature has a pretty high resolution of action => no need to worry about the resolution limit in the macroscopic world => classical mechanics works well.

However, in the microscopic world, nature’s resolution limit can be approached => “round-off error” may occur => one consequence is the quantization of atomic orbitals (discrete energy levels etc.).

Path Integral and Wave Function

Feynman’s principle: The probability (Subscript[p, A→B]) for a particle to propagate from A to B is given by the square modulus of a complex number (Subscript[K, A→B]) called the transition amplitude

The transition amplitude is given by adding together the contributions of all paths x from A to B.

Block[{randpath},
randpath[p1_,p2_]:=BSplineCurve@Table[(1-t)p1+t p2+{0,t(1-t)RandomReal[{-1,1}]},{t,0,1,0.05}];
Graphics[{{Opacity[0.5],AbsoluteThickness[1],Table[{ColorData[0,i],randpath[{0,0},{1,0}]},{i,10}]},Point/@{{0,0},{1,0}},
Text["\!\( A\)",{0,0},{1.2,0}],Text["\!\( B\)",{1,0},{-1.2,0}]}]]

The contribution of each particular path is proportional to (e^(i S[x]/ℏ)), where S[x] is the action of the path x.

In the limit of ℏ→0, the classical path (OverBar[x]) (that satisfies δ S[(OverBar[x])]=0) will dominate the transition amplitude,

(Subscript[K, A→B])∼(e^(i S[(OverBar[x])]/ℏ)).

Quantum mechanics reduces to classical mechanics in the limit of ℏ->0.

(Subscript[K, A→B])=(∫_C)(Subscript[K, A→C])(Subscript[K, C→B]).

(Subscript[K, ((Subscript[x, 0]),(Subscript[t, 0]))->((Subscript[x, N]),(Subscript[t, N]))])=∫d(Subscript[x, 1])…d(Subscript[x, N-1])(Subscript[K, ((Subscript[x, 0]),(Subscript[t, 0]))→((Subscript[x, 1]),(Subscript[t, 1]))])… (Subscript[K, ((Subscript[x, N-1]),(Subscript[t, N-1]))→((Subscript[x, N]),(Subscript[t, N]))]).

ψ((Subscript[x, k+1]),(Subscript[t, k+1]))=∫d(Subscript[x, k])(Subscript[K, ((Subscript[x, k]),(Subscript[t, k]))->((Subscript[x, k+1]),(Subscript[t, k+1]))])ψ((Subscript[x, k]),(Subscript[t, k])).

ψ((Subscript[x, k+1]),(Subscript[t, k+1]))∝∫ⅅ[x]exp((i/ℏ)S[x])ψ((Subscript[x, k]),(Subscript[t, k])),

for each given initial point (Subscript[x, k]), integrate over paths x(t) subject to the boundary conditions x((Subscript[t, k]))=(Subscript[x, k]) and x((Subscript[t, k+1]))=(Subscript[x, k+1]),

finally integrate over choices of initial point (Subscript[x, k]).

Deriving the Schrödinger Equation

Action in a Time Slice

x(t)=(Subscript[x, 0])+v t,

v=(((Subscript[x, 1])-(Subscript[x, 0]))/((Subscript[t, 1])-(Subscript[t, 0])))=(((Subscript[x, 1])-(Subscript[x, 0]))/(δ t)).

S[x]=(1/2)((m((((Subscript[x, 1])-(Subscript[x, 0]))/(δ t))))^2)δ t=(m/(2δ t))((((Subscript[x, 1])-(Subscript[x, 0])))^2).

Path Integral in a Time Slice

ψ((Subscript[x, 1]),t+δ t)∝∫d(Subscript[x, 0])exp((i/ℏ)S[x])ψ((Subscript[x, 0]),t)
=∫d(Subscript[x, 0])exp(((i m)/(2ℏ δ t))((((Subscript[x, 1])-(Subscript[x, 0])))^2))ψ((Subscript[x, 0]),t).

The proportional sign “∝” implies that the normalization factor is not determined yet. (It will be determined later.)

ψ((Subscript[x, 0]),t)=ψ((Subscript[x, 1])+a,t)
=ψ((Subscript[x, 1]),t)+a ψ'((Subscript[x, 1]),t)+((a^2)/(2!))ψ''((Subscript[x, 1]),t)+((a^3)/(3!))(ψ^((3)))((Subscript[x, 1]),t)+…
=(∑_(n=0)^∞)((a^n)/(n!))(∂_(Subscript[x, 1])^n)ψ((Subscript[x, 1]),t).

ψ((Subscript[x, 1]),t+δ t)∝(∑_(n=0)^∞)∫da exp(((i m)/(2ℏ δ t))(a^2))((a^n)/(n!))(∂_(Subscript[x, 1])^n)ψ((Subscript[x, 1]),t).

The idea is that the time-evolved wavefunction can be expressed as the original wavefunction “dressed” by its (different orders of) derivatives.

Plot[{Exp[-x^2/2],-x Exp[-x^2/2],(x^2-1)Exp[-x^2/2]},{x,-5,5},PlotRange->All,Epilog->{{#1,Text["\!\("<>#2<>" ψ(x)\)",#3]}&@@@{{Blue,"",{1.8,0.9}},{Red,"∂\_x",{3.5,-0.5}},{Green,"∂\_x\%2",{2.2,-0.9}}}},FrameLabel->{"\!\(x \)","\!\(∂\_x\%n ψ(x)\)"},ImageSize->200]

For example, ψ(x) is a wave packet.

ψ(x)+λ ψ'(x): shift the wave packet around.

ψ(x)+λ ψ''(x): expand or shrink the wave packet.

Locality of Physics: the time evolution should only involve local modifications of the wavefunction ψ(x) (mostly within the light-cone) in each step.

Computing the Coefficients (Subscript[λ, n])

(Subscript[λ, n])=((1+(((-1))^n))/2)((Sqrt[π])/((2^n)Γ(1+(n/2))))(((-((i m)/(2ℏ δ t))))^(-((1+n)/2))).

Evaluate the integral in Eq. () for (Subscript[λ, n]).

FunctionExpand[Integrate[Exp[I c a^2]a^n/n!,{a,-∞,∞}, GenerateConditions->False]]

The first term (1+(((-1))^n))/2 just discriminates even and odd n.

((1+(((-1))^n))/2)=({{{, ({{1, if n∈even,}, {0, if n∈odd.}})}})

So as long as n∈odd, (Subscript[λ, n])=0. We only need to consider the case of even n.

For even n, the first several (Subscript[λ, n]) are given by

(Subscript[λ, 0])=(Sqrt[π])(((-((i m)/(2ℏ δ t))))^(-1/2)),
(Subscript[λ, 2])=((Sqrt[π])/4)(((-((i m)/(2ℏ δ t))))^(-3/2))=((i )/4)(((2ℏ δ t)/m))(Subscript[λ, 0]),
(Subscript[λ, 4])=((Sqrt[π])/32)(((-((i m)/(2ℏ δ t))))^(-5/2))=-(1/32)(((((2ℏ δ t)/m)))^2)(Subscript[λ, 0]),
…

Compute the first several (Subscript[λ, n]) for even integer n using Eq. ().

Block[{λ},
λ=Association@Table[
n->Sqrt[π]/(2^n Gamma[1+n/2])(-I c)^(-(1+n)/2),{n,0,4,2}];
{λ,λ/λ[0]}]

Determining the Normalization

ψ(x,t+δ t)∝(Subscript[λ, 0])(1+(i/4)(((2ℏ δ t)/m))(∂_x^2)-(1/32)(((((2ℏ δ t)/m)))^2)(∂_x^4)+…)ψ(x,t).

If we take δ t=0, all higher order terms vanishes,

ψ(x,t)∝(Subscript[λ, 0])ψ(x,t).

So obviously, the normalization factor should be such to cancelled out (Subscript[λ, 0]).

ψ(x,t+δ t)=(1+(i/4)(((2ℏ δ t)/m))(∂_x^2)-(1/32)(((((2ℏ δ t)/m)))^2)(∂_x^4)+…)ψ(x,t).

Taking the Limit of δ t->0

(∂_t)ψ(x,t)=(lim_(δ t->0))((ψ(x,t+δ t)-ψ(x,t))/(δ t))
=(lim_(δ t->0))(1/(δ t))((i/4)(((2ℏ δ t)/m))(∂_x^2)-(1/32)(((((2ℏ δ t)/m)))^2)(∂_x^4)+…)ψ(x,t)

Only the first term survives under the limit δ t->0,

(∂_t)ψ(x,t)=((i ℏ)/(2m))(∂_x^2)ψ(x,t).

All the higher order terms will have higher powers in δ t, so they should all vanish under the limit δ t->0.

Adding Potential Energy

Δ S=-V(x)δ t.

ψ(x,t+δ t)=(e^(i Δ S/ℏ))(Subscript[ψ, 0])(x,t+δ t)
=(e^(-i V(x)δ t/ℏ))(Subscript[ψ, 0])(x,t+δ t)
=(1-(i/ℏ)V(x)δ t+…)(Subscript[ψ, 0])(x,t+δ t),

ψ(x,t+δ t)=(1-(i/ℏ)V(x)δ t+…)(1+(i/4)(((2ℏ δ t)/m))(∂_x^2)+…)ψ(x,t)
=(1+(i/4)(((2ℏ δ t)/m))(∂_x^2)-(i/ℏ)V(x)δ t+…)ψ(x,t).

(∂_t)ψ(x,t)=((i ℏ)/(2m))(∂_x^2)ψ(x,t)-(i/ℏ)V(x)ψ(x,t),

Time-Independent Case

Derive Eq. () from Eq. ().

i ℏ(∂_t)ψ(x,t)=i ℏ(∂_t)(ψ(x)(e^(-i E t/ℏ)))
=i ℏ((∂_t)(e^(-i E t/ℏ)))ψ(x)
=i ℏ(-((i E)/ℏ)(e^(-i E t/ℏ)))ψ(x)
=E ψ(x) (e^(-i E t/ℏ)).

-((ℏ^2)/(2m))(∂_x^2)ψ(x,t)+V(x)ψ(x,t)
=-((ℏ^2)/(2m))(∂_x^2)ψ(x)(e^(-i E t/ℏ))+V(x)ψ(x)(e^(-i E t/ℏ))
=(-((ℏ^2)/(2m))(∂_x^2)+V(x))ψ(x)(e^(-i E t/ℏ)).

(-((ℏ^2)/(2m))(∂_x^2)+V(x))ψ(x)=E ψ(x).

(Subscript[E, n]): eigen energies,

(Subscript[ψ, n])(x): corresponding eigen wavefunctions,

Semiclassical Approach

WKB Approximation (General)

Goal: find approximate solution of Eq. (), keeping only the leading quantum effects (i.e., the leading terms of ℏ).

To the leading (0th) order of ℏ, the action function S(x,t) is governed by

also known as the Hamilton-Jacobi equation.

In general, given the Hamiltonian function H(x,p,t) of the system, the Hamilton-Jacobi equation reads

(∂_t)S(x,t)+H(x,∇S(x,t),t)=0.

Eq. () is a special case of Eq. () for a particle moving in 1D with H(x,p,t)=(1/(2m))(p^2)+V(x).

Physical meaning of action derivatives:

Energy: (negative) rate of action accumulation in time

Momentum: action accumulation rate in space (along every direction), or the spatial gradient of action

To the next leading (1st) order of ℏ,

(∂_t)A(x,t)+(1/(2m))(2(∂_x)A(x,t)(∂_x)S(x,t)+A(x,t) (∂_x^2)S(x,t))=0,

which determines A(x,t) from the solution of S(x,t).

Derive Eq. () and Eq. ().

i ℏ(∂_t)(A (e^(i S/ℏ)))=(-((ℏ^2)/(2m))(∂_x^2)+V)A (e^(i S/ℏ)).

i ℏ(∂_t)(A (e^(i S/ℏ)))=i ℏ ((∂_t)A (e^(i S/ℏ))+A (∂_t)(e^(i S/ℏ)))
=i ℏ ((∂_t)A (e^(i S/ℏ))+(i/ℏ)A (∂_t)S (e^(i S/ℏ)))
=(i ℏ (∂_t)A -A (∂_t)S)(e^(i S/ℏ)),

-((ℏ^2)/(2m))(∂_x^2)(A (e^(i S/ℏ)))=-((ℏ^2)/(2m))(∂_x)((∂_x)A (e^(i S/ℏ))+A (∂_x)(e^(i S/ℏ)))
=-((ℏ^2)/(2m))(∂_x)((∂_x)A (e^(i S/ℏ))+(i/ℏ)A (∂_x)S (e^(i S/ℏ)))
=-((ℏ^2)/(2m))(∂_x)(((∂_x)A+(i/ℏ)A (∂_x)S)(e^(i S/ℏ)))
=-((ℏ^2)/(2m))((∂_x)((∂_x)A+(i/ℏ)A (∂_x)S)(e^(i S/ℏ))+((∂_x)A+(i/ℏ)A (∂_x)S)(∂_x)(e^(i S/ℏ)))
=-((ℏ^2)/(2m))(((∂_x^2)A+(i/ℏ)(∂_x)A (∂_x)S+(i/ℏ)A (∂_x^2)S)(e^(i S/ℏ))+(i/ℏ)((∂_x)A+(i/ℏ)A (∂_x)S)(∂_x)S (e^(i S/ℏ)))
=-((ℏ^2)/(2m))((∂_x^2)A+((2i)/ℏ)(∂_x)A (∂_x)S+(i/ℏ)A (∂_x^2)S-(1/(ℏ^2))A ((((∂_x)S))^2))(e^(i S/ℏ))
=(1/(2m))(-(ℏ^2)(∂_x^2)A-2i ℏ(∂_x)A (∂_x)S-i ℏ A (∂_x^2)S+A ((((∂_x)S))^2))(e^(i S/ℏ)).

𝒪((ℏ^0)): to the leading order, we have

-(∂_t)S A (e^(i S/ℏ))=((1/(2m))((((∂_x)S))^2)+V)A (e^(i S/ℏ)).

Eliminating common factors, we obtain

(∂_t)S+(1/(2m))((((∂_x)S))^2)+V=0.

𝒪((ℏ^1)): to the next leading order, we have

i ℏ (∂_t)A (e^(i S/ℏ))=(1/(2m))(-2i ℏ(∂_x)A (∂_x)S-i ℏ A (∂_x^2)S)(e^(i S/ℏ)),

Eliminating common factors, we obtain

(∂_t)A +(1/(2m))(2(∂_x)A (∂_x)S+A (∂_x^2)S)=0.

Solutions of Hamilton-Jacobi Equation

Particle in the Free Space

V(x)=0.

S(x,0)=p x=m (Subscript[v, 0])x,

(Subscript[v, 0]) is the initial velocity of the particle at time t=0,

p=m (Subscript[v, 0]) is the momentum of the particle, which will remain conserved,

E=(1/(2m))(p^2)=(1/2)m (Subscript[v, 0]^2) is the (kinetic) energy of the particle.

Manipulate[DynamicModule[{m=1,V,S,xmin,xmax,tmin,tmax},
V[x_]:=0;
{{xmin,xmax},{tmin,tmax}}={{-5,5},{-5,5}};
Quiet[S=NDSolveValue[{D[S[x,t],t]+(D[S[x,t],x])^2/(2m)+V[x]==0,S[x,0]==m v0 x},S,{x,xmin,xmax},{t,tmin,tmax}]];
Block[{cfun="Rainbow"},
ContourPlot[S[x,t],{x,xmin,xmax},{t,tmin,tmax},PlotRange->All,Contours->Range[-200,200,2],ContourStyle->Directive[AbsoluteThickness[1],Opacity[0.4],Black],ColorFunction->cfun,FrameLabel->{"\!\( x\)","\!\( t\)"},PlotLabel->"\!\(S(x,t)\)",ImageSize->180]]],{{v0,1,"\!\(v\_0\)"},-3,3,ImageSize->Small},Paneled->False]

The contours of S(x,t) are wave fronts (equal-phase surface) in the spacetime.

The solution S(x,t) in Eq. () corresponds to a plane wave solution of the wavefunction

This turns out to be the exact solution of the Schrödinger equation for a free particle (the WKB approximation becomes exact in this case).

Particle under a Constant Force

V(x)=-F x,

(∂_t)S+(1/(2m))((((∂_x)S))^2)+V(x)=0,

Manipulate[DynamicModule[{V,S,xmin,xmax,tmin,tmax},
V[x_]:=-F x;
{{xmin,xmax},{tmin,tmax}}={{-5,5},{-5,5}};
Quiet[S=NDSolveValue[{D[S[x,t],t]+(D[S[x,t],x])^2/(2m)+V[x]==0,S[x,0]==m v0 x},S,{x,xmin,xmax},{t,tmin,tmax}]];
Block[{cfun="Rainbow"},
Show[ContourPlot[S[x,t],{x,xmin,xmax},{t,tmin,tmax},PlotRange->All,Contours->Range[-200,200,2],ContourStyle->Directive[AbsoluteThickness[1],Opacity[0.4],Black],ColorFunction->cfun,FrameLabel->{"\!\( x\)","\!\( t\)"},PlotLabel->"\!\(S(x,t)\)",ImageSize->180],ContourPlot[x==1/2(F/m)t^2+v0 t+m v0^2/(2F),{x,xmin,xmax},{t,tmin,tmax},ContourStyle->Black]]]],
{{m,1},0.5,10,ImageSize->Small},{{F,1},-2,2,ImageSize->Small},
{{v0,0,"\!\(v\_0\)"},-3,3,ImageSize->Small},Paneled->False]

In the m→∞ limit, (∂_t)S=-V(x), such that

S(x,t)=-V(x) t=F x t,

creating a growing spatial gradient of the action

p=(∂_x)S=F t,

corresponding to a momentum that increases in time.

When m is finite, the kinetic energy ((((∂_x)S))^2)/(2m) grows with the momentum p=(∂_x)S, which also contributes to the total energy and alters the rate of action accumulation in time. => This leads to curvature in the constant-action contours, signaling acceleration in the particle’s motion.

The classical trajectory (in black) of the particle corresponds to the family of stationary points of S(x,t) in the spacetime, which turns out to form a parabola

x=(1/2)a (((t+((Subscript[v, 0])/a)))^2).

(Subscript[v, 0])=(OverDot[x])(t=0) is the initial velocity of the particle at time t=0,

a is the acceleration of the particle. It increases with F and decreases with m, and can be verified to follow

a=(F/m),

which recovers Newton’s 2nd law (F=m a).

WKB Approximation (Time-Independent)

S(x,t)=W(x)-E t,

The spatial part W(x) of the action satisfies the stationary Hamilton-Jacobi equation, reduced from Eq. (),

(1/(2m))((((∂_x)W(x)))^2)+V(x)=E.

Derive Eq. () from Eq. ().

-E+(1/(2m))((((∂_x)W(x)))^2)+V(x)=0,

Given a time-independent Hamiltonian H(x,p), a more general form of Eq. () is

H(x,∇W(x))=E.

Eq. () can be solved by introducing the momentum function p(x) — the rate that the action is accumulated in space,

such that Eq. () becomes an algebraic equation

(((p(x))^2)/(2m))+V(x)=E,

with the solution(s) given by

Then the solution of W(x) can be reconstructed by integration

W(x)=(∫^x)p(x')dx'=(∫^x)(Sqrt[2m (E-V(x'))])dx'.

Eq. () also reduces to its stationary form

(∂_x)log A(x)=-(1/2)(∂_x)log p(x),

whose solution is

A(x)=(C/(Sqrt[p(x)])).

Derive Eq. ().

(1/(2m))(2(∂_x)A(∂_x)W+A(∂_x^2)W)=0,

2(∂_x)A(∂_x)W=-A(∂_x^2)W
=>(((∂_x)A)/A)=-(1/2)(((∂_x^2)W)/((∂_x)W))
=>(∂_x)log A=-(1/2)(∂_x)log((∂_x)W).

(∂_x)log A=-(1/2)(∂_x)log p.

Classically allowed regions (V(x)<E):

p(x)∈ℝ, the WKB wavefunction ψ(x) exhibits wavy behavior.

Both ± solutions of p(x) are valid, corresponding to right-moving and left-moving waves.

Classically forbidden regions (V(x)>E):

p(x)∈𝕀, the WKB wavefunction ψ(x) exhibits exponential decay (or grow) behavior

ψ(x)≈(C/(Sqrt[|p(x)|]))exp(∓(1/ℏ)(∫^x)|p(x')|dx').

Only one of the ± solutions of p(x) will be valid, which corresponding to the decaying wave, as the particle’s probability density must diminish as it enters the classical forbidden regions. The invalid solution will correspond to an growing wave.

Transition region (V(x)→E): p(x)→0, the amplitude diverges as ((p(x))^(-1/2)), and the WKB wavefunction is ill-defined. Joining the WKB wave function across the transition region is a rather complicated task, more can be found in Ref. [].

Scattering states

Quantum climbing: the potential top is lower than the energy level E.

Block[{m=1,Ε=2,ℏ=0.05,dx=0.02,V,p,W,ψ,xmin,xmax},
V[x_]:=1.9Exp[-x^2];
p[x_]:=Sqrt[2m(Ε-V[x])];
W[x_?NumericQ]:=Quiet[NIntegrate[p[x1],{x1,-3,x}]];
ψ[x_?NumericQ]:=Exp[I W[x]/ℏ]/Sqrt[Abs[p[x]]];
{xmin,xmax}={-3,3};
Plot[V[x],{x,xmin,xmax},PlotStyle->Black,Axes->False,FrameLabel->{"\!\( x\)","\!\(V(x)\)"},Prolog->{AbsoluteThickness[1],LightGray,InfiniteLine[{#,Ε}&/@{0,1}]},Epilog->{Translate[{Red,Cases[ListLinePlot[Block[{data},data=Table[If[Abs[p[x]]dx>ℏ/8,{x,Re@ψ[x]},{x,Null}],{x,xmin,xmax,dx}];
data[[All,2]]=0.5data[[All,2]]/Abs[data[[1,2]]];data]],_Line,∞],Text["\!\(Re ψ(x)\)",{2,1},{0,0}]},{0,Ε}]},PlotRange->{0,4},PlotRangePadding->{Scaled[0.02],Scaled[0.05]},ImageSize->180]
]

Quantum tunneling: the potential top is higher than the energy level E.

Block[{m=1,Ε=2,ℏ=0.05,dx=0.02,V,p,W,ψ,xmin,xmax},
V[x_]:=2.1Exp[-x^2];
p[x_]:=Sqrt[2m(Ε-V[x])];
W[x_?NumericQ]:=Quiet[NIntegrate[p[x1],{x1,-3,x}]];
ψ[x_?NumericQ]:=Exp[I W[x]/ℏ]/Sqrt[Abs[p[x]]];
{xmin,xmax}={-3,3};
Plot[V[x],{x,xmin,xmax},PlotStyle->Black,Axes->False,FrameLabel->{"\!\( x\)","\!\(V(x)\)"},Prolog->{AbsoluteThickness[1],LightGray,InfiniteLine[{#,Ε}&/@{0,1}]},Epilog->{Translate[{Red,Cases[ListLinePlot[Block[{data},data=Table[If[Abs[p[x]]dx>ℏ/10,{x,Re@ψ[x]},{x,Null}],{x,xmin,xmax,dx}];
data[[All,2]]=0.5data[[All,2]]/Abs[data[[1,2]]];data]],_Line,∞],Text["\!\(Re ψ(x)\)",{2,1},{0,0}]},{0,Ε}]},PlotRange->{0,4},PlotRangePadding->{Scaled[0.02],Scaled[0.05]},ImageSize->180]
]

Bound state: the potential grows higher than the energy level E towards both sides.

Block[{m=1,Ε=2,ℏ=0.05,dx=0.02,V,xmin,xmax,p,W,ψ},
V[x_]:=x^2/2-x^3/6+x^4/24;
{xmin,xmax}=x/.NSolve[V[x]==Ε,x,Reals];
p[x_]:=If[Ε>V[x],Sqrt[2m(Ε-V[x])],Sign[x]Sqrt[2m(Ε-V[x])]];
W[x_?NumericQ]:=Quiet[NIntegrate[p[x1],{x1,0,x}]];
ψ[x_?NumericQ]:=Exp[I W[x]/ℏ]/Sqrt[Abs[p[x]]];
Plot[V[x],{x,-2,3},PlotStyle->Black,Axes->False,FrameLabel->{"\!\( x\)","\!\(V(x)\)"},Prolog->{{AbsoluteThickness[1],LightGray,InfiniteLine[{#,Ε}&/@{0,1}],Table[InfiniteLine[{x,#}&/@{0,1}],{x,{xmin,xmax}}]},{Text["\!\(x\_min\)",{xmin,0},{0,-0.5}],Text["\!\(x\_max\)",{xmax,0},{0,-0.5}]}},Epilog->{Translate[{Red,Cases[ListLinePlot[Block[{data},data=Table[If[Abs[p[x]]dx>ℏ/6,{x,Re@ψ[x]},{x,Null}],{x,-2,3,dx}];
data[[All,2]]=data[[All,2]]/Max[Abs[DeleteCases[data[[All,2]],Null]]];data]],_Line,∞],Text["\!\(Re ψ(x)\)",{0,1}]},{0,Ε}]},ImageSize->180]
]

(Subscript[x, min]), (Subscript[x, max]): two classical turning points, at which E=V(x) and the particle will be bounced back in the classical limit.

Total phase acquired by the wavefunction between the classical turning points is given by

Θ((Subscript[x, min])→(Subscript[x, max]))=(1/ℏ)W((Subscript[x, min])→(Subscript[x, max])),

where W is the corresponding action,

W((Subscript[x, min])→(Subscript[x, max]))=(∫_(Subscript[x, min])^(Subscript[x, max]))p(x)dx=(∫_(Subscript[x, min])^(Subscript[x, max]))(Sqrt[2m (E-V(x))])dx.

Wikipedia, WKB approximation.

Bohr-Sommerfeld Quantization

Intuition: Consider a sine wave, with one node pinned to (Subscript[x, min]), how to pin another node to (Subscript[x, max]) by varying Θ((Subscript[x, min])→(Subscript[x, max]))?

Manipulate[Block[{k,col},{k,col}=If[Abs[Mod[Θ/π,1,-1/2]]<0.1,{Round[Θ/π]π,Green},{Θ,Red}];Plot[Sin[k x],{x,-0.5,1.5},PlotStyle->col,Prolog->{{LightGray,AbsoluteThickness[1],InfiniteLine[{#,0}&/@{0,1}],Table[InfiniteLine[{x,#}&/@{0,1}],{x,{0,1}}]}},Epilog->{{Opacity[0.5],Rectangle[{-1,-2},{0,2}],Rectangle[{1,-2},{2,2}]}},PlotLabel->StringTemplate["\!\(Θ(x\_min→x\_max)=`1`π\)"][If[#==1,"",#]&@(k/π)],Axes->False,FrameLabel->{None,"\!\(ψ(x)\)"},FrameTicks->{{None,None},{{{0,"\!\(x\_min\)"},{1,"\!\(x\_max\)"}},None}},PlotRange->{-1.2,1.2},AspectRatio->1/3,ImageSize->220]],{{Θ,π},0.5π,5π},Paneled->False]

Bohr-Sommerfeld quantization condition: To confine the wave in the region [(Subscript[x, min]),(Subscript[x, max])], we must pin the wave nodes on both turning points, which requires the phase acquired between the turning points to be an integer of π, i.e.

for n=1,2,3,….

Consider the potential

V(x)=(1/2)m (ω^2)(x^2).

ω - angular frequency of the oscillator.

Plot[1/2 x^2,{x,-3,3},GridLines->{{-2,0,2},{2}},Epilog->{Text["\!\(-x\_0\)",{-2,0},{0,-0.6}],Text["\!\(x\_0\)",{2,0},{0,-0.6}],Text["\!\( E\)",{-3,2},{-1,0}]},Axes->False,PlotStyle->Black,FrameLabel->{"\!\( x\)","\!\(V(x)\)"},ImageSize->180]

Let ±(Subscript[x, 0]) be the turning points, at which E=V(x), such that

E=(1/2)m (ω^2)(Subscript[x, 0]^2).

The Bohr-Sommerfeld quantization condition Eq. () requires

((m ω)/ℏ)(∫_(-(Subscript[x, 0]))^(Subscript[x, 0]))(Sqrt[(Subscript[x, 0]^2)-(x^2)])dx=((π m ω (Subscript[x, 0]^2))/(2ℏ))=n π,

which sets (Subscript[x, 0]^2)=2n ℏ/(m ω). By Eq. (), the energy that correspond to such turning points is

(Subscript[E, n])=n ℏ ω,

for n=1,2,3,….

(Subscript[E, n])=(n+(1/2))ℏ ω,