Causal evolution of wave packets

Numerical investigation of causality in the relativistic-Schrödinger system
A companion to Causal evolution of wave packets arXiv:1610.00764 [quant-ph]

Michal Eckstein, Tomasz Miller ©

michal@eckstein.pl, T.Miller@mini.pw.edu.pl

Definitions

A time counting procedure

Time[x_]:=Module[{t,y},
{t,y}=x;
Print[y];
Print["Time taken: ",t," seconds."]
]

The relativistic-Schrödinger system

ℏ=1; (* Adopting natural units *)
c=1;
m=1; (* Fixing m = 1 *)

Δ[p_]:=Sqrt[p^2+m^2]; (* The dispersion relation *)
ψ[t_?NumericQ,x_?NumericQ,ψ0p_,par_,B_,p0_]:=(2 π ℏ)^(-1/2) NIntegrate[Exp[-I/ℏ (Δ[p] t-p x)] ψ0p[p,par,p0],{p,-B+p0,B+p0}]; (* Evolution of the wave-packet via Fourier transform. *)
(*  ψ0p - initial wave packet in the momentum representation,
	par - parameter of the initial wave packet (for instance width),
	B - numerical cutoff,
    p0 - average momentum of the wave packet. *)

P[T_?NumericQ,a_?NumericQ,b_?NumericQ,ψ0p_,par_,B_,p0_]:=NIntegrate[Abs[ψ[T,x,ψ0p,par,B,p0]]^2,{x,a-c T,b+c T}];(* Probability density integrated over the interval [a-cT,b+cT] *)

Causality breakdown measure

ψ0[x_,ψ0p_,par_,p0_]:=FourierTransform[ψ0p[p,par,p0],p,x,FourierParameters->{0,1/ℏ}]

M[t_,a_,b_,ψ0p_,par_,B_,p0_]:=Max[0,Integrate[((Abs[ψ0[x,ψ0p,par,p0]])^2),{x,a,b}]-P[t,a,b,ψ0p,par,B,p0]]

Discrete plots

DPlot[ψ0p_,par_,B_,p0_,t1_,t2_,x1_,x2_,N_,st_]:=DiscretePlot[M[t,x1,x2,ψ0p,par,B,p0],{t,t1,t2,(t2-t1)/N},PlotStyle->{st},Filling->None,Joined->True] 
(* Yields a plot of the quantity M(t,[(Subscript[x, 1]),(Subscript[x, 2])]) for t∈[(Subscript[t, 1]),(Subscript[t, 2])] with the resolution of N points. *)
DPlotX[ψ0p_,par_,B_,p0_,t_,x1_,x2_,N_,st_]:=DiscretePlot[M[t,-y,y,ψ0p,par,B,p0],{y,x1,x2,(x2-x1)/N},PlotStyle->{st},Filling->None,Joined->True] 
(* Yields a plot of the quantity M(t,[(Subscript[x, 1]),(Subscript[x, 2])]) for t∈[(Subscript[t, 1]),(Subscript[t, 2])] with the resolution of N points. *)

Gaussian wave packets

Definitions

ψ0G[x_]:=(1/(((π d))^(1/4)))(E^((-(x^2))/(2d)))(E^(I p0 x))

Integrate[((Abs[ψ0G[x]])^2),{x,-∞,∞},Assumptions->{d>0,p0∈Reals,x∈Reals}]

FullSimplify[InverseFourierTransform[ψ0G[x],x,p,FourierParameters->{0,1/ℏ}],Assumptions->{d>0,ℏ>0}]

ψ0p[p_,d_,p0_]:=(((d^(1/4)) (E^(-(1/2) d (((p-p0))^2))))/(π^(1/4)))

Integrate[((Abs[ψ0p[p,d,p0]])^2),{p,-∞,∞},Assumptions->{d>0,p0∈Reals}]

The case d = 1, (Subscript[p, 0 ])=0

d=1;p0=0;

m=1;

Setting the value of B - the numerical integration limit

B=(10/(Sqrt[d]));

N@ψ0p[B,d,p0]

The evolution of probability density

Time@Timing@Plot[{((Abs[ψ0[x,ψ0p,d,p0]])^2),((Abs[ψ[1,x,ψ0p,d,B,p0]])^2),((Abs[ψ[2,x,ψ0p,d,B,p0]])^2),((Abs[ψ[3,x,ψ0p,d,B,p0]])^2)},{x,-5,5}]

Time@Timing@Plot3D[((Abs[ψ[t,x,ψ0p,d,B,p0]])^2),{x,-5,5},{t,0,5}]

General view

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,d,B,p0],{t,0,10},{a,0,10},PlotLegends->Automatic]

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,d,B,p0],{t,0,3},{a,0,5},PlotLegends->Automatic]

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,d,B,p0],{t,0,2},{a,2.5,4},PlotLegends->Automatic]

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,d,B,p0],{t,0.3,1},{a,2.5,4},PlotLegends->Automatic]

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,d,B,p0],{t,0.7,0.9},{a,2.8,3},PlotLegends->Automatic]

Time@Timing@FindMaximum[{M[t,-a,a,ψ0p,d,B,p0],0.75<t<0.9,2.85<a<2.93},{{t,0.8},{a,2.89}},PrecisionGoal->8,AccuracyGoal->8]

Checking that M(t,[-3,3]) indeed equals 0 for short times

Time@Timing@Plot[P[t,-3,3,ψ0p,d,B,p0],{t,0,0.2}]

Time@Timing@Show[{
DPlot[ψ0p,d,B,p0,0,3,-3,3,100,Red],
DPlot[ψ0p,d,B,p0,0,3,-3.1,3.1,100,Blue],
DPlot[ψ0p,d,B,p0,0,3,-3.2,3.2,100,Green],
DPlot[ψ0p,d,B,p0,0,3,-3.3,3.3,100,Yellow]},AxesLabel->{"t","d"}]

p2=Show[{DPlot[ψ0p,d,B,p0,0,0.2,-3,3,200,{Black,Thick}],DPlot[ψ0p,d,B,p0,0,0.2,-3.1,3.1,100,{Black,Thick,Dashing[Tiny]}],DPlot[ψ0p,d,B,p0,0,0.2,-3.2,3.2,100,{Black,Thick,DotDashed}],DPlot[ψ0p,d,B,p0,0,0.2,-3.3,3.3,100,{Black,Thick,Dashing[Medium]}]},AxesLabel->{Style[TraditionalForm[t],FontSize->20],Style[TraditionalForm[HoldForm@M[t,((Subscript[ψ, 0])^G),(Subscript[(𝒦^("    ")), a])]],FontSize->20]},TicksStyle->FontSize->20,Ticks->{{0,0.1,0.2},{0,{2 (10^(-6)),"2×(10^(-6))"},{4 (10^(-6)),"4×(10^(-6))"}}},PlotRange->All,ImageSize->300]

P3=Show[{DPlot[ψ0p,d,B,p0,0,3,-3,3,200,{Black,Thick}],DPlot[ψ0p,d,B,p0,0,3,-3.1,3.1,100,{Black,Thick,Dashing[Tiny]}],DPlot[ψ0p,d,B,p0,0,3,-3.2,3.2,100,{Black,Thick,DotDashed}],DPlot[ψ0p,d,B,p0,0,3,-3.3,3.3,100,{Black,Thick,Dashing[Medium]}]},AxesLabel->{Style[TraditionalForm[t],Large],Style[TraditionalForm[HoldForm@M[t,((Subscript[ψ, 0])^G),(Subscript[(𝒦^("    ")), a])]],Large]},TicksStyle->Large,Ticks->{Automatic,{0,{(10^(-5)),"(10^(-5))"},{2 (10^(-5)),"2×(10^(-5))"},{3 (10^(-5)),"3×(10^(-5))"}}}, Epilog->{Inset[p2,{2.4,2.8 (10^(-5))}],Inset[Framed[Column[{Style["  a",Large],
LineLegend[{Black},{Style[3,FontSize->20]},LegendMarkerSize->{30,10}],LineLegend[{Dashing[Tiny]},{Style[3.1,FontSize->20]},LegendMarkerSize->{30,10}],
LineLegend[{DotDashed},{Style[3.2,FontSize->20]},LegendMarkerSize->{30,10}],
LineLegend[{Dashing[Medium]},{Style[3.3,FontSize->20]},LegendMarkerSize->{30,10}]
}],RoundingRadius->5,FrameMargins->10],Scaled[{0.9,0.4}]]}]

The time-scale (Subscript[t, 0])

𝒦 = [-4,4]

Time@Timing@DPlot[ψ0p,d,B,p0,0,0.1,-4,4,100,Blue]

Time@Timing@DPlot[ψ0p,d,B,p0,0,0.01,-4,4,100,Blue]

𝒦 = [-5,5]

Time@Timing@DPlot[ψ0p,d,B,p0,0,0.01,-5,5,100,Blue]

Time@Timing@DPlot[ψ0p,d,B,p0,0,0.001,-5,5,100,Blue]

Time@Timing@DPlot[ψ0p,d,B,p0,0,0.0001,-5,5,100,Blue]

The time-scale (Subscript[t, 1])

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,d,B,p0],{t,0.7,0.9},{a,2.8,3},PlotLegends->Automatic]

Choosing the adequate accuracy for FindMaximum

DPlot[ψ0p,d,B,p0,0.79,0.82,-2.9,2.9,50,Blue]

DPlot[ψ0p,d,B,p0,0.81,0.815,-2.9,2.9,50,Blue]

Time@Timing@FindMaximum[{M[t,-2.9,2.9,ψ0p,d,B,p0],0.5<t<1},{t,0.8},PrecisionGoal->8,AccuracyGoal->8]

Time@Timing[t1=Table[{a,FindMaximum[{M[t,-a,a,ψ0p,d,B,p0],0.5<t<1},{t,0.8},PrecisionGoal->8,AccuracyGoal->8]},{a,2.6,5,0.1}]]

t1//TableForm

The time-scale (Subscript[t, 2])

𝒦 = [-3,3]

Time@Timing@DPlot[ψ0p,d,B,p0,0,10,-3,3,100,Blue]

𝒦 = [-4,4]

Time@Timing@DPlot[ψ0p,d,B,p0,0,10,-4,4,100,Blue]

Time@Timing@DPlot[ψ0p,d,B,p0,6,10,-4,4,50,Blue]

𝒦 = [-5,5]

Time@Timing@DPlot[ψ0p,d,B,p0,0,10,-5,5,100,Blue]

Time@Timing@DPlot[ψ0p,d,B,p0,6,10,-5,5,50,Blue]

Time@Timing@DPlot[ψ0p,d,B,p0,6,12,-5,5,50,Blue]

Time@Timing@DPlot[ψ0p,d,B,p0,6,15,-5,5,50,Blue]

Time@Timing@DPlot[ψ0p,d,B,p0,6,20,-5,5,50,Blue]

Asymmetric 𝒦

For[i=-4,i≤4,i++,
For[j=-4,j≤4,j++,Print[(Subscript[p, 0])," = ",p0,", ",𝒦," = [",-4+i/2,",",4+j/2,"]"];
Time@Timing@DPlot[ψ0p,d,B,p0,0,3,-4+i/2,4+j/2,50,Blue];
Print["================================================"]
]
]

Time@Timing@FindMaximum[{M[t,a,b,ψ0p,d,B,p0],0<t<1,-6<a<b<6},{{t,0.8},{a,-3},{b,3}},PrecisionGoal->8,AccuracyGoal->8]

Thus, for (Subscript[p, 0]) = 0, the maximum of M is attained for symmetric K = {0} x [-a,a].

Dependence on (Subscript[p, 0])

p0 = 1/2

p0=1/2;

The evolution of probability density

Time@Timing@Plot[{((Abs[ψ0[x,ψ0p,d,p0]])^2),((Abs[ψ[1,x,ψ0p,d,B,p0]])^2),((Abs[ψ[2,x,ψ0p,d,B,p0]])^2),((Abs[ψ[3,x,ψ0p,d,B,p0]])^2)},{x,-5,5}]

Time@Timing@Plot3D[((Abs[ψ[t,x,ψ0p,d,B,p0]])^2),{x,-5,5},{t,0,5}]

General view

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,d,B,p0],{t,0,2},{a,2.5,4},PlotLegends->Automatic]

Time@Timing@FindMaximum[{M[t,-a,a,ψ0p,d,B,p0],0.7<t<0.9,3<a<3.1},{{t,0.8},{a,3.05}},PrecisionGoal->8,AccuracyGoal->8]

Asymmetric 𝒦

For[i=-4,i≤4,i++,
For[j=-4,j≤4,j++,Print[(Subscript[p, 0])," = ",p0,", ",𝒦," = [",-4+i/2,",",4+j/2,"]"];
Time@Timing@DPlot[ψ0p,d,B,p0,0,3,-4+i/2,4+j/2,50,Blue];
Print["================================================"]
]
]

Time@Timing@FindMaximum[{M[t,a,b,ψ0p,d,B,p0],0<t<1,-6<a<6,a<b<6},{{t,0.8},{a,-7/2},{b,5/2}},PrecisionGoal->8,AccuracyGoal->8]

For (Subscript[p, 0]) = 1/2, the maximum of M is attained for asymmetric K = {0} x [a,b].

p0 = 1

p0=1;

The evolution of probability density

Time@Timing@Plot[{((Abs[ψ0[x,ψ0p,d,p0]])^2),((Abs[ψ[1,x,ψ0p,d,B,p0]])^2),((Abs[ψ[2,x,ψ0p,d,B,p0]])^2),((Abs[ψ[3,x,ψ0p,d,B,p0]])^2)},{x,-5,5}]

Time@Timing@Plot3D[((Abs[ψ[t,x,ψ0p,d,B,p0]])^2),{x,-5,5},{t,0,5}]

General view

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,d,B,p0],{t,0,2},{a,2.5,4},PlotLegends->Automatic]

Time@Timing@FindMaximum[{M[t,-a,a,ψ0p,d,B,p0],0.7<t<0.9,3<a<3.1},{{t,0.8},{a,3.05}},PrecisionGoal->8,AccuracyGoal->8]

Asymmetric 𝒦

For[i=-4,i≤4,i++,For[j=-4,j≤4,j++,Print[(Subscript[p, 0])," = ",p0,", ",𝒦," = [",-4+i/2,",",4+j/2,"]"];
Time@Timing@DPlot[ψ0p,d,B,p0,0,3,-4+i/2,4+j/2,50,Blue];
Print["================================================"]]]

Time@Timing@FindMaximum[{M[t,a,b,ψ0p,d,B,p0],0<t<1,-6<a<6,a<b<6},{{t,0.8},{a,-7/2},{b,5/2}},PrecisionGoal->8,AccuracyGoal->8]

For (Subscript[p, 0]) = 1, the maximum of M is attained for asymmetric K = {0} x [a,b].

The dependendence on d

p0=0;

d=1/10;B=(10/(Sqrt[d]));

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,d,B,p0],{t,0.5,0.9},{a,0.6,0.7},PlotLegends->Automatic]

Time@Timing@FindMaximum[{M[t,-a,a,ψ0p,d,B,p0],0.6<t<0.85,0.6<a<0.66},{t,0.68},{a,0.63},PrecisionGoal->6,AccuracyGoal->6]

d=(10^(-2));B=(10/(Sqrt[d]));

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,d,B,p0],{t,0.5,0.9},{a,0.15,0.18},PlotLegends->Automatic]

Time@Timing@FindMaximum[{M[t,-a,a,ψ0p,d,B,p0],0.55<t<0.75,0.16<a<0.17},{t,0.64},{a,0.165},PrecisionGoal->6,AccuracyGoal->6]

d=(10^(-3));B=(10/(Sqrt[d]));

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,d,B,p0],{t,0.3,0.9},{a,0.04,0.06},PlotLegends->Automatic]

Time@Timing@FindMaximum[{M[t,-a,a,ψ0p,d,B,p0],0.45<t<0.75,0.045<a<0.051},{t,0.58},{a,0.05},PrecisionGoal->6,AccuracyGoal->6]

d=(10^(-4));B=(10/(Sqrt[d]));

Time@Timing@DPlot[ψ0p,d,B,p0,0.5,0.55,-0.015,0.015,20,Blue]

Time@Timing@DPlotX[ψ0p,d,B,p0,0.53,0.014,0.016,20,Blue]

Time@Timing@M[0.53,-0.015,0.015,ψ0p,d,B,p0]

Time@Timing@M[0.528,-0.0148,0.0148,ψ0p,d,B,p0]

d=(10^(-5));B=(10/(Sqrt[d]));

Time@Timing@DPlot[ψ0p,d,B,p0,0.46,0.52,-0.0047,0.0047,20,Blue]

Time@Timing@DPlotX[ψ0p,d,B,p0,0.48,0.0045,0.0055,20,Blue]

Time@Timing@M[0.48,-0.0047,0.0047,ψ0p,d,B,p0]

d=(10^(-6));B=(10/(Sqrt[d]));

Time@Timing@M[0.43,-0.0015,0.0015,ψ0p,d,B,p0]

The massless case

m=0;

Re-defining the procedures

ℏ=1; (* Adopting natural units *)
c=1;

Δ[p_]:=Sqrt[p^2+m^2]; (* The dispersion relation *)
ψ[t_?NumericQ,x_?NumericQ,ψ0p_,par_,B_,p0_]:=(2 π ℏ)^(-1/2) NIntegrate[Exp[-I/ℏ (Δ[p] t-p x)] ψ0p[p,par,p0],{p,-B+p0,B+p0}]; (* Evolution of the wave-packet via Fourier transform. *)
(*  ψ0p - initial wave packet in the momentum representation,
	par - parameter of the initial wave packet (for instance width),
	B - numerical cutoff,
    p0 - average momentum of the wave packet. *)

P[T_?NumericQ,a_?NumericQ,b_?NumericQ,ψ0p_,par_,B_,p0_]:=NIntegrate[Abs[ψ[T,x,ψ0p,par,B,p0]]^2,{x,a-c T,b+c T}];(* Probability density integrated over the interval [a-cT,b+cT] *)

ψ0[x_,ψ0p_,par_,p0_]:=FourierTransform[ψ0p[p,par,p0],p,x,FourierParameters->{0,1/ℏ}]

M[t_?NumericQ,a_?NumericQ,b_?NumericQ,ψ0p_,par_,B_,p0_]:=Max[0,NIntegrate[((Abs[ψ0[x,ψ0p,par,p0]])^2),{x,a,b}]-P[t,a,b,ψ0p,par,B,p0]]

DPlot[ψ0p_,par_,B_,p0_,t1_,t2_,x1_,x2_,N_,st_]:=DiscretePlot[M[t,x1,x2,ψ0p,par,B,p0],{t,t1,t2,(t2-t1)/N},PlotStyle->{st},Filling->None,Joined->True] 
(* Yields a plot of the quantity M(t,[(Subscript[x, 1]),(Subscript[x, 2])]) for t∈[(Subscript[t, 1]),(Subscript[t, 2])] with the resolution of N points. *)
DPlotX[ψ0p_,par_,B_,p0_,t_,x1_,x2_,N_,st_]:=DiscretePlot[M[t,-y,y,ψ0p,par,B,p0],{y,x1,x2,(x2-x1)/N},PlotStyle->{st},Filling->None,Joined->True] 
(* Yields a plot of the quantity M(t,[(Subscript[x, 1]),(Subscript[x, 2])]) for t∈[(Subscript[t, 1]),(Subscript[t, 2])] with the resolution of N points. *)
(* DCPlot[t1_,t2_,x1_,x2_,N_]:=ListContourPlot@Flatten[Table[{t,y,M[t,-y,y]},{t,t1,t2,(t2-t1)/N},{y,x1,x2,(x2-x1)/N}],1] *)

d=1;p0=0;B=(10/(Sqrt[d]));

The evolution of probability density

Time@Timing@Plot[{((Abs[ψ0[x,ψ0p,d,p0]])^2),((Abs[ψ[1,x,ψ0p,d,B,p0]])^2),((Abs[ψ[2,x,ψ0p,d,B,p0]])^2),((Abs[ψ[3,x,ψ0p,d,B,p0]])^2)},{x,-5,5},PlotRange->All]

Time@Timing@Plot3D[((Abs[ψ[t,x,ψ0p,d,B,p0]])^2),{x,-5,5},{t,0,5},PlotRange->All]

General view

For[i=1,i≤12,i++,
Print[𝒦," = [",-i/2,",",i/2,"]"];
Time@Timing@DPlot[ψ0p,d,B,p0,0,3,-i/2,i/2,50,Blue];
Print["================================================"]
]

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,d,B,p0],{t,0,10},{a,0,10},PlotLegends->Automatic]

Time@Timing[r00=ContourPlot[M[t,-a,a,ψ0p,d,B,p0],{t,0,10},{a,0,6},PlotPoints->40,PlotLegends->BarLegend[Automatic,LegendMarkerSize->200,LegendMargins->5,LegendLabel->Style[TraditionalForm[HoldForm@M[t,((Subscript[ψ, 0])^G),(Subscript[(𝒦^("    ")), a])]],FontSize->15],LabelStyle->FontSize->14],ColorFunction->ColorData[{"GrayTones","Reverse"}],ContourStyle->Automatic,FrameLabel->{t,a},FrameStyle->FontSize->15]]

The time-scale (Subscript[t, 0])

𝒦 = [-3,3]

Time@Timing@DPlot[ψ0p,d,B,p0,0,0.1,-3,3,100,Blue]

Time@Timing@DPlot[ψ0p,d,B,p0,0,0.01,-3,3,100,Blue]

𝒦 = [-4,4]

Time@Timing@DPlot[ψ0p,d,B,p0,0,0.01,-4,4,100,Blue]

Time@Timing@DPlot[ψ0p,d,B,p0,0,0.001,-4,4,100,Blue]

Time@Timing@DPlot[ψ0p,d,B,p0,0,0.0001,-4,4,100,Blue]

𝒦 = [-5,5]

Time@Timing@DPlot[ψ0p,d,B,p0,0,0.0001,-5,5,100,Blue]

The time-scale (Subscript[t, 2])

𝒦 = [-(3/2), (3/2)]

Time@Timing@DPlot[ψ0p,d,B,p0,0,10,-3/2,3/2,50,Blue]

Time@Timing@DPlot[ψ0p,d,B,p0,0,100,-3/2,3/2,50,Blue]

𝒦 = [-3, 3]

Time@Timing@DPlot[ψ0p,d,B,p0,0,10,-3,3,50,Blue]

Time@Timing@DPlot[ψ0p,d,B,p0,0,100,-3,3,50,Blue]

𝒦 = [-3/2, 3/2]

M[100,-3/2,3/2,ψ0p,d,B,p0]

M[200,-3/2,3/2,ψ0p,d,B,p0]

M[100,-1.4,1.4,ψ0p,d,B,p0]

M[100,-1.6,1.6,ψ0p,d,B,p0]

M[100,-1.55,1.55,ψ0p,d,B,p0]

M[100,-1.51,1.51,ψ0p,d,B,p0]

M[100,-1.49,1.49,ψ0p,d,B,p0]

M[100,-1.48,1.48,ψ0p,d,B,p0]

M[100,-1.47,1.47,ψ0p,d,B,p0]

M[100,-1.46,1.46,ψ0p,d,B,p0]

M[100,-1.45,1.45,ψ0p,d,B,p0]

M[200,-1.46,1.46,ψ0p,d,B,p0]

Exponential wave packets

Definitions

Re-defining the procedures

ℏ=1; (* Adopting natural units *)
c=1;
m=1; (* Fixing m = 1 *)

Δ[p_]:=Sqrt[p^2+m^2]; (* The dispersion relation *)
ψ[t_?NumericQ,x_?NumericQ,ψ0p_,par_,B_,p0_]:=(2 π ℏ)^(-1/2) NIntegrate[Exp[-I/ℏ (Δ[p] t-p x)] ψ0p[p,par,p0],{p,-B+p0,B+p0}]; (* Evolution of the wave-packet via Fourier transform. *)
(*  ψ0p - initial wave packet in the momentum representation,
	par - parameter of the initial wave packet (for instance width),
	B - numerical cutoff,
    p0 - average momentum of the wave packet. *)

P[T_?NumericQ,a_?NumericQ,b_?NumericQ,ψ0p_,par_,B_,p0_]:=NIntegrate[Abs[ψ[T,x,ψ0p,par,B,p0]]^2,{x,a-c T,b+c T}];(* Probability density integrated over the interval [a-cT,b+cT] *)

ψ0[x_,ψ0p_,par_,p0_]:=FourierTransform[ψ0p[p,par,p0],p,x,FourierParameters->{0,1/ℏ}]

M[t_,a_,b_,ψ0p_,par_,B_,p0_]:=Max[0,Integrate[((Abs[ψ0[x,ψ0p,par,p0]])^2),{x,a,b}]-P[t,a,b,ψ0p,par,B,p0]]

DPlot[ψ0p_,par_,B_,p0_,t1_,t2_,x1_,x2_,N_,st_]:=DiscretePlot[M[t,x1,x2,ψ0p,par,B,p0],{t,t1,t2,(t2-t1)/N},PlotStyle->{st},Filling->None,Joined->True] 
(* Yields a plot of the quantity M(t,[(Subscript[x, 1]),(Subscript[x, 2])]) for t∈[(Subscript[t, 1]),(Subscript[t, 2])] with the resolution of N points. *)
DPlotX[ψ0p_,par_,B_,p0_,t_,x1_,x2_,N_,st_]:=DiscretePlot[M[t,-y,y,ψ0p,par,B,p0],{y,x1,x2,(x2-x1)/N},PlotStyle->{st},Filling->None,Joined->True] 
(* Yields a plot of the quantity M(t,[(Subscript[x, 1]),(Subscript[x, 2])]) for t∈[(Subscript[t, 1]),(Subscript[t, 2])] with the resolution of N points. *)
(* DCPlot[t1_,t2_,x1_,x2_,N_]:=ListContourPlot@Flatten[Table[{t,y,M[t,-y,y]},{t,t1,t2,(t2-t1)/N},{y,x1,x2,(x2-x1)/N}],1] *)

ψ0E[x_]:=(Sqrt[(α/2)])Sech[α x](E^(I p0 x))

Integrate[((Abs[ψ0E[x]])^2),{x,-∞,∞},Assumptions->{α>0,p0∈Reals,x∈Reals}]

FullSimplify[InverseFourierTransform[ψ0E[x],x,p,FourierParameters->{0,1/ℏ}],Assumptions->{α>0,p0∈Reals}]

ψ0p[p_,α_,p0_]:=(((Sqrt[π]) Sech[(((p-p0) π)/(2 α))])/(2 (Sqrt[α])))

Integrate[((Abs[ψ0p[p,α,p0]])^2),{p,-∞,∞},Assumptions->{α>0,p0∈Reals}]

Clear[B]

p0=0;
B[α_]:=30 α

ψ0p[B[α],α,0]//N

α = 3

α=3;

The evolution of probability density

Time@Timing@Plot[{((Abs[ψ0[x,ψ0p,α,p0]])^2),((Abs[ψ[1,x,ψ0p,α,B[α],p0]])^2),((Abs[ψ[2,x,ψ0p,α,B[α],p0]])^2),((Abs[ψ[3,x,ψ0p,α,B[α],p0]])^2)},{x,-5,5},PlotRange->All]

Time@Timing@Plot3D[((Abs[ψ[t,x,ψ0p,α,B[α],p0]])^2),{x,-5,5},{t,0,5},PlotRange->All]

General view

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,α,B[α],p0],{t,0,3},{a,0,5},PlotLegends->Automatic]

For[i=0,i≤10,i++,
Print[(Subscript[p, 0])," = ",p0,", ",𝒦," = [",-1-i/5,",",1+i/5,"]"];
Time@Timing@DPlot[ψ0p,α,B[α],p0,0,2,-1-i/5,1+i/5,20,Blue];
Print["================================================"]
]

M[0.8,-7/5,7/5,ψ0p,α,B[α],p0]

M[0.8,-1.3,1.3,ψ0p,α,B[α],p0]

M[0.8,-1.5,1.5,ψ0p,α,B[α],p0]

M[0.81,-1.4,1.4,ψ0p,α,B[α],p0]

M[0.79,-1.4,1.4,ψ0p,α,B[α],p0]

M[0.78,-1.4,1.4,ψ0p,α,B[α],p0]

α = 2

α=2;

The evolution of probability density

Time@Timing@Plot[{((Abs[ψ0[x,ψ0p,α,p0]])^2),((Abs[ψ[1,x,ψ0p,α,B[α],p0]])^2),((Abs[ψ[2,x,ψ0p,α,B[α],p0]])^2),((Abs[ψ[3,x,ψ0p,α,B[α],p0]])^2)},{x,-5,5},PlotRange->All]

Time@Timing@Plot3D[((Abs[ψ[t,x,ψ0p,α,B[α],p0]])^2),{x,-5,5},{t,0,5},PlotRange->All]

General view

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,α,B[α],p0],{t,0,3},{a,0,5},PlotLegends->Automatic]

For[i=0,i≤10,i++,
Print[(Subscript[p, 0])," = ",p0,", ",𝒦," = [",-3-i/5,",",3+i/5,"]"];
Time@Timing@DPlot[ψ0p,α,B[α],p0,0,2,-3-i/5,3+i/5,20,Blue];
Print["================================================"]
]

M[0.8,-3.2,3.2,ψ0p,α,B[α],p0]

M[0.8,-3.3,3.3,ψ0p,α,B[α],p0]

M[0.8,-3.1,3.1,ψ0p,α,B[α],p0]

M[0.81,-3.2,3.2,ψ0p,α,B[α],p0]

M[0.82,-3.2,3.2,ψ0p,α,B[α],p0]

M[0.83,-3.2,3.2,ψ0p,α,B[α],p0]

M[0.84,-3.2,3.2,ψ0p,α,B[α],p0]

α = 5/3

α=5/3;

The evolution of probability density

Time@Timing@Plot[{((Abs[ψ0[x,ψ0p,α,p0]])^2),((Abs[ψ[1,x,ψ0p,α,B[α],p0]])^2),((Abs[ψ[2,x,ψ0p,α,B[α],p0]])^2),((Abs[ψ[3,x,ψ0p,α,B[α],p0]])^2)},{x,-5,5},PlotRange->All]

Time@Timing@Plot3D[((Abs[ψ[t,x,ψ0p,α,B[α],p0]])^2),{x,-5,5},{t,0,5},PlotRange->All]

General view

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,α,B[α],p0],{t,0,3},{a,0,8},PlotLegends->Automatic]

For[i=0,i≤10,i++,
Print[(Subscript[p, 0])," = ",p0,", ",𝒦," = [",-5-i/5,",",5+i/5,"]"];
Time@Timing@DPlot[ψ0p,α,B[α],p0,0,2,-5-i/5,5+i/5,20,Blue];
Print["================================================"]
]

Time@Timing@DPlot[ψ0p,α,B[α],p0,0,0.1,-5.2,5.2,20,Blue]

The apparent bump for t ∈ [0.0.1] is an artefact of numerics.

M[0.8,-5.2,5.2,ψ0p,α,B[α],p0]

M[0.81,-5.2,5.2,ψ0p,α,B[α],p0]

M[0.82,-5.2,5.2,ψ0p,α,B[α],p0]

M[0.83,-5.2,5.2,ψ0p,α,B[α],p0]

M[0.84,-5.2,5.2,ψ0p,α,B[α],p0]

M[0.85,-5.2,5.2,ψ0p,α,B[α],p0]

M[0.84,-5.1,5.1,ψ0p,α,B[α],p0]

M[0.84,-5.3,5.3,ψ0p,α,B[α],p0]

α = 3/2

α=3/2;

The evolution of probability density

Time@Timing@Plot[{((Abs[ψ0[x,ψ0p,α,p0]])^2),((Abs[ψ[1,x,ψ0p,α,B[α],p0]])^2),((Abs[ψ[2,x,ψ0p,α,B[α],p0]])^2),((Abs[ψ[3,x,ψ0p,α,B[α],p0]])^2)},{x,-5,5},PlotRange->All]

Time@Timing@Plot3D[((Abs[ψ[t,x,ψ0p,α,B[α],p0]])^2),{x,-5,5},{t,0,5},PlotRange->All]

General view

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,α,B[α],p0],{t,0,3},{a,0,10},PlotLegends->Automatic]

For[i=0,i≤10,i++,
Print[(Subscript[p, 0])," = ",p0,", ",𝒦," = [",-7-i/5,",",7+i/5,"]"];
Time@Timing@DPlot[ψ0p,α,B[α],p0,0,2,-7-i/5,7+i/5,20,Blue];
Print["================================================"]
]

Time@Timing@DPlot[ψ0p,α,B[α],p0,0,0.1,-7.4,7.4,20,Blue]

The apparent bump for t ∈ [0.0.1] is an artefact of numerics.

M[0.8,-7.4,7.4,ψ0p,α,B[α],p0]

M[0.81,-7.4,7.4,ψ0p,α,B[α],p0]

M[0.82,-7.4,7.4,ψ0p,α,B[α],p0]

M[0.83,-7.4,7.4,ψ0p,α,B[α],p0]

M[0.84,-7.4,7.4,ψ0p,α,B[α],p0]

M[0.85,-7.4,7.4,ψ0p,α,B[α],p0]

M[0.845,-7.4,7.4,ψ0p,α,B[α],p0]

M[0.846,-7.4,7.4,ψ0p,α,B[α],p0]

M[0.847,-7.4,7.4,ψ0p,α,B[α],p0]

M[0.848,-7.4,7.4,ψ0p,α,B[α],p0]

M[0.849,-7.4,7.4,ψ0p,α,B[α],p0]

M[0.848,-7.5,7.5,ψ0p,α,B[α],p0]

M[0.848,-7.3,7.3,ψ0p,α,B[α],p0]

α = 1

α=1;

The evolution of probability density

Time@Timing@Plot[{((Abs[ψ0[x,ψ0p,α,p0]])^2),((Abs[ψ[1,x,ψ0p,α,B[α],p0]])^2),((Abs[ψ[2,x,ψ0p,α,B[α],p0]])^2),((Abs[ψ[3,x,ψ0p,α,B[α],p0]])^2)},{x,-5,5}]

Time@Timing@Plot3D[((Abs[ψ[t,x,ψ0p,α,B[α],p0]])^2),{x,-5,5},{t,0,5}]

General view

Time-scale T = 1

For[i=1,i≤100,i++,
Print[(Subscript[p, 0])," = ",p0,", ",𝒦," = [",-i/10,",",i/10,"]"];
Time@Timing@DPlot[ψ0p,α,B[α],p0,0,1,-i/10,i/10,20,Blue];
Print["================================================"]
]

Time-scale T = 10

For[i=1,i≤100,i++,
Print[(Subscript[p, 0])," = ",p0,", ",𝒦," = [",-i/10,",",i/10,"]"];
Time@Timing@DPlot[ψ0p,α,B[α],p0,0,10,-i/10,i/10,20,Blue];
Print["================================================"]
]

Time-scale T = 0.1

For[i=1,i≤100,i++,
Print[(Subscript[p, 0])," = ",p0,", ",𝒦," = [",-i/10,",",i/10,"]"];
Time@Timing@DPlot[ψ0p,α,B[α],p0,0,0.1,-i/10,i/10,20,Blue];
Print["================================================"]
]

Time-scale T = 0.01

For[i=1,i≤100,i++,
Print[(Subscript[p, 0])," = ",p0,", ",𝒦," = [",-i/10,",",i/10,"]"];
Time@Timing@DPlot[ψ0p,α,B[α],p0,0,0.01,-i/10,i/10,20,Blue];
Print["================================================"]
]

α ⟶ ∞

α=10;

For[i=1,i≤10,i++,
Print["α = ",10,", ",𝒦," = [",-i/10,",",i/10,"]"];
Time@Timing@DPlot[ψ0p,α,B[α],p0,0,2,-i/10,i/10,20,Blue];
Print["================================================"]
]

Time@Timing@FindMaximum[{M[t,-a,a,ψ0p,α,B[α],p0],0<t<1,0<a<2/5},{{t,0.8},{a,1/5}},PrecisionGoal->4,AccuracyGoal->4]

α=100;

For[i=1,i≤3,i++,
Print["α = ",10,", ",𝒦," = [",-i/100,",",i/100,"]"];
Time@Timing@DPlot[ψ0p,α,B[α],p0,0,2,-i/100,i/100,20,Blue];
Print["================================================"]
]

Time@Timing@FindMaximum[{M[t,-a,a,ψ0p,α,B[α],p0],0<t<1,0<a<0.1},{{t,0.5},{a,0.02}},PrecisionGoal->4,AccuracyGoal->4]

α=1000;

For[i=1,i≤3,i++,
Print["α = ",10,", ",𝒦," = [",-i/1000,",",i/1000,"]"];
Time@Timing@DPlot[ψ0p,α,B[α],p0,0,1,-i/1000,i/1000,20,Blue];
(* Time@Timing@FindMaximum[{M[t,-4+i/2,4+j/2,ψ0p,d,B,p0],0<t<1.5},{t,0.8},PrecisionGoal->8,AccuracyGoal->8]; *)
Print["================================================"]
]

M[0.5,-0.002,0.002,ψ0p,α,B[α],p0]

M[0.4,-0.002,0.002,ψ0p,α,B[α],p0]

M[0.5,-0.003,0.003,ψ0p,α,B[α],p0]

M[0.6,-0.002,0.002,ψ0p,α,B[α],p0]

α=10000;

Time@Timing@M[0.4,-0.0002,0.0002,ψ0p,α,B[α],p0]

α=100000;

Time@Timing@M[0.35,-0.00002,0.00002,ψ0p,α,B[α],p0]

Time@Timing@M[0.3,-0.00002,0.00002,ψ0p,α,B[α],p0]

Time@Timing@M[0.4,-0.00002,0.00002,ψ0p,α,B[α],p0]

α=1000000;

Time@Timing@M[0.3,-0.000002,0.000002,ψ0p,α,B[α],p0]

Sin - exp wave packets

Definitions

p0=0;

Re-defining the procedures

Clear[M]

ℏ=1; (* Adopting natural units *)
c=1;
m=1; (* Fixing m = 1 *)

Δ[p_]:=Sqrt[p^2+m^2]; (* The dispersion relation *)
ψ[t_?NumericQ,x_?NumericQ,ψ0p_,par_,B_,p0_]:=(2 π ℏ)^(-1/2) NIntegrate[Exp[-I/ℏ (Δ[p] t-p x)] ψ0p[p,par,p0],{p,-B+p0,B+p0}]; (* Evolution of the wave-packet via Fourier transform. *)
(*  ψ0p - initial wave packet in the momentum representation,
	par - parameter of the initial wave packet (for instance width),
	B - numerical cutoff,
    p0 - average momentum of the wave packet. *)

P[T_?NumericQ,a_?NumericQ,b_?NumericQ,ψ0p_,par_,B_,p0_]:=NIntegrate[Abs[ψ[T,x,ψ0p,par,B,p0]]^2,{x,a-c T,b+c T}];(* Probability density integrated over the interval [a-cT,b+cT] *)

ψ0[x_,ψ0p_,par_,p0_]:=FourierTransform[ψ0p[p,par,p0],p,x,FourierParameters->{0,1/ℏ}]

M[t_?NumericQ,a_?NumericQ,b_?NumericQ,ψ0p_,par_,B_,p0_]:=Max[0,NIntegrate[((Abs[ψ0[x,ψ0p,par,p0]])^2),{x,a,b}]-P[t,a,b,ψ0p,par,B,p0]]

DPlot[ψ0p_,par_,B_,p0_,t1_,t2_,x1_,x2_,N_,st_]:=DiscretePlot[M[t,x1,x2,ψ0p,par,B,p0],{t,t1,t2,(t2-t1)/N},PlotStyle->{st},Filling->None,Joined->True] 
(* Yields a plot of the quantity M(t,[(Subscript[x, 1]),(Subscript[x, 2])]) for t∈[(Subscript[t, 1]),(Subscript[t, 2])] with the resolution of N points. *)
DPlotX[ψ0p_,par_,B_,p0_,t_,x1_,x2_,N_,st_]:=DiscretePlot[M[t,-y,y,ψ0p,par,B,p0],{y,x1,x2,(x2-x1)/N},PlotStyle->{st},Filling->None,Joined->True] 
(* Yields a plot of the quantity M(t,[(Subscript[x, 1]),(Subscript[x, 2])]) for t∈[(Subscript[t, 1]),(Subscript[t, 2])] with the resolution of N points. *)

ψ0SE[x_,α_]:=Nr[α]((Sin[x])/x) Sech[α x](E^(I p0 x))

Nr[α_]:=((NIntegrate[((Abs[((Sin[x])/x) Sech[α x]])^2),{x,((-B[α])/(α^2)),((B[α])/(α^2))}])^(-1/2)) (* The normalisation constant *)

B[α_]:=30α

Normalisation check

Table[NIntegrate[((Abs[ψ0SE[x,α]])^2),{x,((-B[α])/(α^2)),((B[α])/(α^2))}],{α,0.1,4,0.1}]

Clear[α]

Computing the inverse Fourier transform via convolution

InverseFourierTransform[((Sin[x])/x),x,p]//FullSimplify

InverseFourierTransform[Sech[α x],x,p]//FullSimplify

Convolve[UnitBox[r/2],Sech[(r π)/(2 α)],r,p,Assumptions->α>0]

ψ0p[p_,α_,p0_]:=Nr[α] (1/2) (Sqrt[(π/2)])(1/α) ((2α)/π) (Gudermannian[(((1+p-p0) π)/(2α))]+Gudermannian[(((1-p+p0) π)/(2α))])

ψ0[x_,ψ0p_,par_,p0_]:=ψ0SE[x,par]

Normalisation check

Table[NIntegrate[((Abs[ψ0p[p,α,p0]])^2),{p,-B[α],B[α]}],{α,0.1,4,0.1}]

ψ0pN[p_,α_,p0_]:=(1/(Sqrt[2π]))NIntegrate[ψ0SE[x,α] (E^(-I p x)),{x,-5B[α]/(α^2),5B[α]/(α^2)}]

Checking the manual computation of the Fourier transform

Time@Timing@Table[Chop[ψ0p[p,α,0]-ψ0pN[p,α,0]],{p,-5,5,0.1},{α,0.1,4,0.1}]

α = 1

Time@Timing@M[0.9,-2.6,2.6,ψ0p,1,B[1],p0]

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,1,B[1],p0],{t,0,3},{a,0,5},PlotLegends->Automatic]

Time@Timing@NMaximize[{M[t,-a,a,ψ0p,1,B[1],p0],0<t<1,0<a<3},{t,a},PrecisionGoal->4,AccuracyGoal->4]

The dependence on α

Time@Timing[datt=Table[{α,NMaximize[{M[t,-a,a,ψ0p,α,B[α],p0],0<t<1,0<a<3},{t,a},PrecisionGoal->4,AccuracyGoal->4]},{α,1/25,4,1/25}]]

dattM=Table[{datt[[i,1]],datt[[i,2,1]]},{i,1,Length@datt}]

Joining the value at α = 0 -- see next Chapter, Section n = 1, Subsection Pm = 1.

dattM=Join[{{0,0.0010632664367498101}},dattM];

ListPlot[dattM,Joined->True,PlotStyle->{Thick,Black},AxesLabel->{Style[TraditionalForm[α],Large],Style[TraditionalForm[HoldForm@(sup_(t≥0))(OverTilde[M])[t,((Subscript[ψ,0])^SE),α]],Large]},TicksStyle->Large]

NotebookSave[]

Compact momentum wave packets

Re-defining the procedures

ℏ=1; (* Adopting natural units *)
c=1;
m=1; (* Fixing m = 1 *)

Δ[p_]:=Sqrt[p^2+m^2]; (* The dispersion relation *)
ψ[t_?NumericQ,x_?NumericQ,ψ0p_,par_,B_,p0_]:=(2 π ℏ)^(-1/2) NIntegrate[Exp[-I/ℏ (Δ[p] t-p x)] ψ0p[p,par,p0],{p,-B+p0,B+p0}]; (* Evolution of the wave-packet via Fourier transform. *)
(*  ψ0p - initial wave packet in the momentum representation,
	par - parameter of the initial wave packet (for instance width),
	B - numerical cutoff,
    p0 - average momentum of the wave packet. *)

P[T_?NumericQ,a_?NumericQ,b_?NumericQ,ψ0p_,par_,B_,p0_]:=NIntegrate[Abs[ψ[T,x,ψ0p,par,B,p0]]^2,{x,a-c T,b+c T}];(* Probability density integrated over the interval [a-cT,b+cT] *)

ψ0[x_,ψ0p_,par_,p0_]:=FourierTransform[ψ0p[p,par,p0],p,x,FourierParameters->{0,1/ℏ}]

M[t_,a_,b_,ψ0p_,par_,B_,p0_]:=Max[0,Integrate[((Abs[ψ0[x,ψ0p,par,p0]])^2),{x,a,b}]-P[t,a,b,ψ0p,par,B,p0]]

DPlot[ψ0p_,par_,B_,p0_,t1_,t2_,x1_,x2_,N_,st_]:=DiscretePlot[M[t,x1,x2,ψ0p,par,B,p0],{t,t1,t2,(t2-t1)/N},PlotStyle->{st},Filling->None,Joined->True] 
(* Yields a plot of the quantity M(t,[(Subscript[x, 1]),(Subscript[x, 2])]) for t∈[(Subscript[t, 1]),(Subscript[t, 2])] with the resolution of N points. *)
DPlotX[ψ0p_,par_,B_,p0_,t_,x1_,x2_,N_,st_]:=DiscretePlot[M[t,-y,y,ψ0p,par,B,p0],{y,x1,x2,(x2-x1)/N},PlotStyle->{st},Filling->None,Joined->True] 
(* Yields a plot of the quantity M(t,[(Subscript[x, 1]),(Subscript[x, 2])]) for t∈[(Subscript[t, 1]),(Subscript[t, 2])] with the resolution of N points. *)

n = 1

Definitions

Clear[p0]

ψ0S[x_]:=((Sin[Pm x])/((Sqrt[Pm]) (Sqrt[π]) x))(E^(I p0 x))

Integrate[((Abs[ψ0S[x]])^2),{x,-∞,∞},Assumptions->{Pm>0,p0∈Reals,x∈Reals}]

FullSimplify[InverseFourierTransform[ψ0S[x],x,p,FourierParameters->{0,1/ℏ}],Assumptions->{Pm>0,p0∈Reals}]

ψ0p[p_,Pm_,p0_]:=((-Sign[p-p0-Pm]+Sign[p-p0+Pm])/(2 (Sqrt[2]) (Sqrt[Pm])))

Integrate[((Abs[ψ0p[p,Pm,p0]])^2),{p,-∞,∞},Assumptions->{Pm>0,p0∈Reals}]

Clear[B]

B[Pm_]:=Pm;

p0=0;

Pm = 1

Pm=1;

The evolution of probability density

Time@Timing@Plot[{((Abs[ψ0[x,ψ0p,Pm,p0]])^2),((Abs[ψ[1,x,ψ0p,Pm,B[Pm],p0]])^2),((Abs[ψ[2,x,ψ0p,Pm,B[Pm],p0]])^2),((Abs[ψ[3,x,ψ0p,Pm,B[Pm],p0]])^2),((Abs[ψ[4,x,ψ0p,Pm,B[Pm],p0]])^2),((Abs[ψ[5,x,ψ0p,Pm,B[Pm],p0]])^2)},{x,-10,10},PlotRange->All]

Time@Timing@Plot3D[((Abs[ψ[t,x,ψ0p,Pm,B[Pm],p0]])^2),{x,-10,10},{t,0,5},PlotRange->All,PlotPoints->50]

General view

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,Pm,B[Pm],p0],{t,0,3},{a,0,5},PlotLegends->Automatic]

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,Pm,B[Pm],p0],{t,0,1},{a,2.5,3.2},PlotLegends->Automatic]

General

For[i=1,i≤10,i++,
Print["Pm = ",i/10];
Time@Timing@NMaximize[{M[t,-a,a,ψ0p,i/10,i/10,p0],0<t<1,0<a<10},{t,a},PrecisionGoal->2,AccuracyGoal->2,MaxIterations->20,Method->"DifferentialEvolution"]
]

For[i=2.0,i≤10,i++,
Print["Pm = ",i];
Time@Timing@NMaximize[{M[t,-a,a,ψ0p,i,i,p0],0<t<1,0<a<3},{t,a},PrecisionGoal->2,AccuracyGoal->2,MaxIterations->20,Method->"DifferentialEvolution"]
]

Time@Timing@NMaximize[{M[t,-a,a,ψ0p,100,i,p0],0<t<0.1,0<a<0.3},{t,a},PrecisionGoal->2,AccuracyGoal->2,MaxIterations->20,Method->"DifferentialEvolution"]

n = 2

Definitions

Clear[p0,Pm]

ψ0S[x_]:=(Sqrt[(3/(2 π (Pm^3)))])(((Sin[Pm x])^2)/(x^2))(E^(I p0 x))

Integrate[((Abs[ψ0S[x]])^2),{x,-∞,∞},Assumptions->{Pm>0,p0∈Reals,x∈Reals}]

FullSimplify[InverseFourierTransform[ψ0S[x],x,p,FourierParameters->{0,1/ℏ}],Assumptions->{Pm>0,p0∈Reals}]

ψ0p[p_,Pm_,p0_]:=(((Sqrt[3]) (2 (-p+p0) Sign[p-p0]+(p-p0-2 Pm) Sign[p-p0-2 Pm]+(p-p0+2 Pm) Sign[p-p0+2 Pm]))/(8 (Pm^(3/2))))

Integrate[((Abs[ψ0p[p,Pm,p0]])^2),{p,-∞,∞},Assumptions->{Pm>0,p0∈Reals}]

Clear[B]

B[Pm_]:=2Pm;

p0=0;

Pm = 1

Pm=1;

The evolution of probability density

Time@Timing@Plot[{((Abs[ψ0[x,ψ0p,Pm,p0]])^2),((Abs[ψ[1,x,ψ0p,Pm,B[Pm],p0]])^2),((Abs[ψ[2,x,ψ0p,Pm,B[Pm],p0]])^2),((Abs[ψ[3,x,ψ0p,Pm,B[Pm],p0]])^2)},{x,-5,5},PlotRange->All]

Time@Timing@Plot3D[((Abs[ψ[t,x,ψ0p,Pm,B[Pm],p0]])^2),{x,-5,5},{t,0,5},PlotRange->All,PlotPoints->50]

General view

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,Pm,B[Pm],p0],{t,0,3},{a,0,5},PlotLegends->Automatic]

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,Pm,B[Pm],p0],{t,0,0.6},{a,2.4,3},PlotLegends->Automatic]

General

For[i=1.0,i≤10,i++,
Print["Pm = ",i/10];
Time@Timing@NMaximize[{M[t,-a,a,ψ0p,i/10,2 i/10,p0],0<t<1,0<a<10},{t,a},PrecisionGoal->2,AccuracyGoal->2,MaxIterations->10,Method->"DifferentialEvolution"]
]

For[i=2.0,i≤10,i++,
Print["Pm = ",i];
Time@Timing@NMaximize[{M[t,-a,a,ψ0p,i,2i,p0],0<t<1,0<a<3},{t,a},PrecisionGoal->2,AccuracyGoal->2,MaxIterations->10,Method->"DifferentialEvolution"]
]

n = 3

Definitions

Clear[p0,Pm]

ψ0S[x_]:=(Sqrt[(20/(11 π (Pm^5)))])(((Sin[Pm x])^3)/(x^3))(E^(I p0 x))

Integrate[((Abs[ψ0S[x]])^2),{x,-∞,∞},Assumptions->{Pm>0,p0∈Reals,x∈Reals}]

Simplify[InverseFourierTransform[ψ0S[x],x,p,FourierParameters->{0,1/ℏ}],Assumptions->{Pm>0,p0∈Reals}]

ψ0p[p_,Pm_,p0_]:=(1/(8 (Pm^(5/2))))(Sqrt[(5/22)]) (-(((-p+p0+3 Pm))^2) Sign[p-p0-3 Pm]+3 (((-p+p0+Pm))^2) Sign[p-p0-Pm]-3 (((p-p0+Pm))^2) Sign[p-p0+Pm]+(((p-p0+3 Pm))^2) Sign[p-p0+3 Pm])

Integrate[((Abs[ψ0p[p,Pm,p0]])^2),{p,-∞,∞},Assumptions->{Pm>0,p0∈Reals}]

Clear[B]

B[Pm_]:=3Pm;

p0=0;

Pm = 1

Pm=1;

The evolution of probability density

Time@Timing@Plot[{((Abs[ψ0[x,ψ0p,Pm,p0]])^2),((Abs[ψ[1,x,ψ0p,Pm,B[Pm],p0]])^2),((Abs[ψ[2,x,ψ0p,Pm,B[Pm],p0]])^2),((Abs[ψ[3,x,ψ0p,Pm,B[Pm],p0]])^2)},{x,-5,5},PlotRange->All]

Time@Timing@Plot3D[((Abs[ψ[t,x,ψ0p,Pm,B[Pm],p0]])^2),{x,-5,5},{t,0,5},PlotRange->All,PlotPoints->50]

General view

Time@Timing@ContourPlot[M[t,-a,a,ψ0p,Pm,B[Pm],p0],{t,0,3},{a,0,5},PlotLegends->Automatic]

General

For[i=1.0,i≤10,i++,
Print["Pm = ",i/10];
Time@Timing@NMaximize[{M[t,-a,a,ψ0p,i/10,3 i/10,p0],0<t<1,0<a<10},{t,a},PrecisionGoal->2,AccuracyGoal->2,MaxIterations->10]
]

For[i=2.0,i≤10,i++,
Print["Pm = ",i];
Time@Timing@NMaximize[{M[t,-a,a,ψ0p,i,3i,p0],0<t<1,0<a<3},{t,a},PrecisionGoal->2,AccuracyGoal->2,MaxIterations->10]
]