Infinite Potential Wall

The energy eigenfunctions and eigenvalues of the infinite potetial wall:

ψ[x_,n_] = (Sqrt[(2/L)])Sin[((n π)/L)x];En[n_] = (((ℏ^2)(π^2))/(2 m (L^2)))(n^2);
nrule = {L-> 1,m-> 1, ℏ-> 1,p-> 10π}; (*Numerical values chosen for plots*)
Ncut = 40; (* We truncate the Hilbert space dimension to this value *)
$Assumptions = n∈Integers && l∈Integers && n>0 && m > 0 && L>0 && ℏ>0 && l>0 && p ∈ Reals && x ∈ Reals && y ∈ Reals && t∈ Reals;

Problem 4:

Consider the particle in the initial state ϕ(x,0) = ((Subscript[c, 1])(Subscript[ψ, 1])(x)+(Subscript[c, 3])(Subscript[ψ, 3])(x)+(Subscript[c, 5])(Subscript[ψ, 5])(x))(e^(ip x/ℏ)) for (Subscript[c, 3])/(Subscript[c, 1])=-1/3, (Subscript[c, 5])/(Subscript[c, 1])=1/5. Compute the state ϕ(x,t) and the corresponding average values of position , momentum etc.
It is similar to problem 3 but we give to the particle initial momentum <p>=p (see nrule above).
Since now the coefficients (Subscript[c, n]) are infinite in number and complex in algebraic form (you can easily try to compute them), we put numerical values from the beginning.

Solution:

Consider the expression of the wave function interms of the energy eigenfunctions:
ϕ(x,0) ≡ϕ(x) =  (∑_(n=1)^∞)(Subscript[c, n])(Subscript[ψ, n])(x), (Subscript[c, n])= (∫_0^L)(Subscript[ψ, n])*(x)ϕ(x)dx
ϕ(x,t) = (∑_(n=1)^∞)(Subscript[c, n])(Subscript[ψ, n])(x) (e^(-i (Subscript[E, n])t/ℏ))
< x > (t) = (∑_(n,m=1)^∞)(Subscript[c, m])*(Subscript[c, n])(e^(-i ((Subscript[E, n])-(Subscript[E, m]))t/ℏ))(∫_0^L)(Subscript[ψ, m])*(x)x (Subscript[ψ, n])(x)dx
< (x^2) > (t) = (∑_(n,m=1)^∞)(Subscript[c, m])*(Subscript[c, n])(e^(-i ((Subscript[E, n])-(Subscript[E, m]))t/ℏ))(∫_0^L)(Subscript[ψ, m])*(x)(x^2) (Subscript[ψ, n])(x)dx

We compute the normalization factor A and define the wave function at t=0

A = (1/(Sqrt[(∫_0^L)(((ψ[x,1]-(1/3)ψ[x,3]+(1/5)ψ[x,5]))^2)\[DifferentialD]x]))

ϕ[x_,0] = A(ψ[x,1]-(1/3)ψ[x,3]+(1/5)ψ[x,5])(E^(I p x /ℏ));

Plot[{Abs[ϕ[x,0]]/.nrule,Re[ϕ[x,0]]/.nrule,Im[ϕ[x,0]]/.nrule},{x,0,L/.nrule}]

We compute the coefficients (Subscript[c, n]):

Table[cn[n] = Chop[NIntegrate[ψ[x,n] ϕ[x,0]/.nrule,{x,0,L/.nrule},MaxRecursion->12,WorkingPrecision->20]],{n,1,Ncut}]

ϕ[x_,t_] = (∑_(n=1)^Ncut)cn[n] ψ[x,n] (E^(-I En[n]t/ℏ)) ;

Absϕ[x_,t_] = Abs[ϕ[x,t]]/.nrule;
Manipulate[Plot[((Absϕ[x,t])^2),{x,0,1},PlotRange->{{0,1},{0,4.5}}],{t,0,4,0.0001}]

Reϕ[x_,t_] = Re[ϕ[x,t]]/.nrule;
Imϕ[x_,t_] = Im[ϕ[x,t]]/.nrule;
Manipulate[Plot[{Reϕ[x,t],Imϕ[x,t]},{x,0,1},PlotRange->{{0,1},{-2,2.2}},PlotStyle->{Red,Blue}],{t,0,4,0.0001}]

Average energy <E>, <(E^2)> and uncetainty ΔE:

avE = (∑_(n=1)^Ncut)((Abs[cn[n]])^2)En[n];avE2 = (∑_(n=1)^Ncut)((Abs[cn[n]])^2)((En[n])^2);ΔE = (Sqrt[avE2-(avE^2)]);
Print[
"<E> = ", avE," = ",avE/En[1]," (Subscript[E, 1])",
"\n\!\(\*Cell[TextData[{
Cell[BoxData[
FormBox[
RowBox[{"<", 
SuperscriptBox["E", "2"]}], TraditionalForm]]],
"> = "
}]]\)", avE2," = ",avE2/((En[1])^2)," ((Subscript[E, 1])^2)",
"\nΔE = ",ΔE," = ",ΔE/En[1]//Simplify," (Subscript[E, 1])"
]

Position and momenta matrix elements:
(Subscript[x, nm])= (∫_0^L)(Subscript[ψ, n])*(x)x (Subscript[ψ, m])(x)dx,(Subscript[(x^2), nm])= (∫_0^L)(Subscript[ψ, n])*(x)(x^2)(Subscript[ψ, m])(x)dx , (Subscript[p, nm])= (∫_0^L)(Subscript[ψ, n])*(x)(OverHat[p]) ψ(x)dx, (Subscript[(p^2), nm])= (∫_0^L)(Subscript[ψ, n])*(x)((OverHat[p])^2) ψ(x)dx

xnm[n_,l_] = (∫_0^L)ψ[x,n]x ψ[x,l]\[DifferentialD]x;xnm[n_,n_] = (∫_0^L)ψ[x,n]x ψ[x,n]\[DifferentialD]x;
x2nm[n_,l_] = (∫_0^L)ψ[x,n](x^2) ψ[x,l]\[DifferentialD]x;x2nm[n_,n_] = (∫_0^L)ψ[x,n](x^2) ψ[x,n]\[DifferentialD]x;
pnm[n_,l_] = (∫_0^L)ψ[x,n](-I ℏ) (∂_x)ψ[x,l]\[DifferentialD]x /. Cos[l π] -> (((-1))^l);pnm[n_,n_] = (∫_0^L)ψ[x,n](-I ℏ) (∂_x)ψ[x,n]\[DifferentialD]x ;
p2nm[n_,l_] = (∫_0^L)ψ[x,n](- (ℏ^2)) (∂_(x,x))ψ[x,l]\[DifferentialD]x ;p2nm[n_,n_] = (∫_0^L)ψ[x,n](- (ℏ^2)) (∂_(x,x))ψ[x,n]\[DifferentialD]x ;

Table[pnm[i,j],{i,1,5},{j,1,5}]//MatrixForm

Here we have to be careful: Due to large arguments in trig/exp functions the numerical errors are large. In particular we find large real parts. Therefore we compute the real part from the beginning. Notice the use of the Refine[] function in order for mathematica to use the fact that t is real and compute the real part.

δ[i_,j_] := KroneckerDelta[i,j];δ2[i_,j_]:=(2-δ[i,j]);
xavϕ[t_] = (∑_(n=1)^Ncut)(∑_(l=n)^Ncut) (δ2[n,l]Refine[Re[cn[n]\[Conjugate] cn[l] xnm[n,l] (E^(-I(En[l]-En[n])t/ℏ))]])/.nrule//N;
x2avϕ[t_] = (∑_(n=1)^Ncut)(∑_(l=n)^Ncut) (δ2[n,l]Refine[Re[cn[n]\[Conjugate] cn[l] x2nm[n,l] (E^(-I(En[l]-En[n])t/ℏ))]])/.nrule//N;
pavϕ[t_] = (∑_(n=1)^Ncut)(∑_(l=n)^Ncut) (δ2[n,l]Refine[Re[cn[n]\[Conjugate] cn[l] pnm[n,l] (E^(-I(En[l]-En[n])t/ℏ))]])/.nrule//N;
p2avϕ[t_] = (∑_(n=1)^Ncut)(∑_(l=n)^Ncut)(δ2[n,l]Refine[Re[ cn[n]\[Conjugate] cn[l] p2nm[n,l] (E^(-I(En[l]-En[n])t/ℏ))]])/.nrule//N;
Δx[t_] = (Sqrt[x2avϕ[t]-((xavϕ[t])^2)]);Δp[t_]=(Sqrt[p2avϕ[t]-((pavϕ[t])^2)]);

{Plot[xavϕ[t]/.nrule,{t,0,1},PlotLabel-> "<x>(t)",PlotStyle->Directive[Red,Thick]],
Plot[(Sqrt[x2avϕ[t]/.nrule]),{t,0,1},PlotLabel-> "<(x^2)((>)^(1/2))(t)",PlotStyle->Directive[Blue,Thick]]}

This should be the same as the value of p in nrule:

pavϕ[0]

{Plot[pavϕ[t]/.nrule,{t,0,1},PlotLabel-> "<p>(t)",PlotStyle->Directive[Red,Thick]],
Plot[(Sqrt[p2avϕ[t]/.nrule]),{t,0,1},PlotLabel-> "<(p^2)((>)^(1/2))(t)",PlotStyle->Directive[Blue,Thick]]}

{Plot[Δx[t]/.nrule,{t,0,1},PlotLabel-> "Δx(t)",PlotStyle->Directive[Red,Thick]],
Plot[Δp[t]/.nrule,{t,0,1},PlotLabel-> "Δp(t)",PlotStyle->Directive[Blue,Thick]],
Plot[{Δx[t]Δp[t]/.nrule,ℏ/2/.nrule},{t,0,1},PlotLabel-> "Δx(t) Δp(t)",PlotStyle->{Directive[Magenta,Thick],Directive[Green,Thick]}]}