This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Numerics Volume, Chapter 1

Consider a quantum version of the undamped Duffing oscillator (a forced, friction-free anharmonic oscillator [✶], [✶])

i ℏ ((∂ψ(x,t))/(∂t))=-((ℏ^2)/(2m))(((∂^2)ψ(x,t))/(∂(t^2)))+(α (x^4)-β (x^2)-γ x cos(ω t))ψ(x,t)

with minimum uncertainty wave packet initial condition, localized around (Subscript[x, 0]) with mean momentum (Subscript[p, 0])

ψ(x,0)=(1/(Sqrt[(Sqrt[2 π]) (Subscript[σ, x])]))exp(-((((x-(Subscript[x, 0])))^2)/(4 (Subscript[σ, x]^2)))+((i (Subscript[p, 0]) x)/ℏ)).

For the parameters m=1, α=1/2, β=10, γ=10, ω=6, (Subscript[x, 0])=-3, (Subscript[p, 0])=10, (Subscript[σ, x])=1/4, solve the Schrödinger equation over a few periods 2π/ω and visualize the resulting ((|ψ(x,t)|)^2) for some ℏ<1 [✶], [✶], [✶].

Make an animation of the corresponding Wigner distribution [✶]

W( x,p, t)=(1/(2π ℏ))(∫_(-∞)^∞)(e^(i (p/ℏ) y))ψ(x-(y/2),t)(OverBar[ψ(x+(y/2),t)])dy

as a function of t. Color regions of positive and negative W( x,p, t) differently. W( x,p, t) obeys the differential equation [✶]

((∂ W( x,p, t))/(∂t))=-(p/m)((∂ W( x,p, t))/(∂x))-(2β x-4 α (x^3)+γ cos(ω t))((∂ W( x,p, t))/(∂p))+(ℏ^3)(α/8) x(((∂^3) W( x,p, t))/(∂(p^3))).

For the corresponding classical Duffing oscillator with the equation of motion

x''(t)-2β x(t)+4 α ((x(t))^3)=γ cos(ω t)

model the time-evolution for (10^6) particles with initial position-momentum distribution according to the above ψ(x,0). Make a corresponding animation of the time-evolution of the classical ensemble.

MinimumUncertaintyPacket[{x0_, p0_, σx_}, ℏ_] = 
1/Sqrt[Sqrt[2 Pi] σx] Exp[-(x - x0)^2/(4 σx^2) + I p0 x/ℏ];

Next, we solve the Schrödinger equation for the Duffing potential. To obtain a numerically accurate solution, we use an eighth order difference scheme and a couple thousand spatial discretization points. In addition, we will use a precision and accuracy goal of about six.

DuffingSchrödingerSolution[{ℏ_, m_, α_, β_, γ_, ω_}, {x0_, p0_, σx_}, 
                           {L_, tMax_, pag_:5}, ppL_] :=
NDSolve[{(* Schrödinger equation *)
         I ℏ D[ψ[x, t], t] == -ℏ^2/(2 m) D[ψ[x, t], {x, 2}] + 
                              (α x^4 - β x^2 - γ x Cos[ω t]) ψ[x, t],
         (* minimum uncertainty initial wave packet *)
         ψ[x, 0] == MinimumUncertaintyPacket[{x0, p0, σx}, ℏ], 
         (* Dirichlet boundary conditions *) ψ[-L, t] == 0, ψ[L, t] == 0},
         ψ, {x, -L, L}, {t, 0, tMax},
        (* option for accurate solution *)
         MaxSteps -> 10^5, PrecisionGoal -> pag, AccuracyGoal -> pag, 
         Method -> {"MethodOfLines",  Method -> StiffnessSwitching,
                    "SpatialDiscretization" -> 
                    {"TensorProductGrid", "DifferenceOrder" -> 8,
                     "MaxPoints" -> ppL, "MinPoints" -> ppL}}]

For the parameters at hand, at x=±8 the value of the Duffing potential is always greater than 1300. Taking a WKB approximation for a free particle with momentum (Subscript[p, 0]) and approximating the classical turning point to (Subscript[x, t])≈-5, this means that at x=±8, the wave function is attenuated by more than 20 orders of magnitude compared with its maximal value inside the double well. So, we can safely set ψ(x=±8,t)=0. Next, we declare the values of the parameters we will use in the following as globally visible variables. We follow the time-evolution of the packet for seven/eight periods.

(* system phase space sizes *) 
L = 8; P = 20; 
(* ℏ, m, and potential parameters *) 
ℏ = 1/3; m = 1; α = 1/2; β = 10; γ = 10;
(* wave packet parameters *) 
x0 = -3; p0 = 10; σx = 1/4; σp = ℏ/(2 σx);
(* external force and time parameters *) 
ω = 6; T = 2 Pi/ω; 𝓅 = 8; tMax = 𝓅 T;
(* animation and density plot parameters *)
framesPerPeriod = 10; pp = 256;
Δt = tMax/(𝓅 framesPerPeriod);

Using now 3600 spatial discretization points gives a solution less in than 30 minutes on a 2GHz computer.

(dsNDSol = DuffingSchrödingerSolution[{ℏ, m, α, β, γ, ω}, 
                                      {x0, p0, σx}, {L, tMax, 5.5}, 3600]) // 
                                                                      Timing

The time-dependent Duffing potential is a double well potential in which the depth/height of the two minima oscillates periodically. As a result, a wave function initially localized in one of the wells and having an average momentum comparable to the height of the middle wall of the potential, will quickly spread out over both wells. The following graphic shows this quite clearly. After each reflection as the outer potential walls, the wave function acquires a complicated interference pattern with the adiabatic bound state structure of the double well imprinted. (For the time-development of localized states in a static double-well potential, see [✶].)

DensityPlot[Evaluate[Abs[ψ[x, t]]^(1/4) /. dsNDSol[[1]]], 
            {x, -L, L}, {t, 0, tMax}, 
            PlotRange -> All, Mesh -> False, PlotPoints -> 600]

As a check for quality of solution, we check wether the wavefunction ψ is still normalized to (∫_(-L)^L)ψ(x,t)(OverBar[ψ(x,t)]) dx=1 at (Subscript[t, max]). At t=7 T, the deviation of the norm is less than 0.03%. Using precision and accuracy goal settings of 6, this error could be reduced to 0.01%. The visualizations to be created later would be visually identical, so we use the faster-to-process solution dsNDSol.

NIntegrate[Evaluate[Abs[ψ[x, tMax]]^2 /. dsNDSol[[1]]], {x, -L, L},
           MaxRecursion -> 24, PrecisionGoal -> 10] // InputForm

For comparison, here is a time evolution that starts in one well and only slowly leaks out into the neighbor well.

DensityPlot[Evaluate[Abs[ψ[x, t]]^(1/4) /. 
(* solve Schrödinger equation *)
(DuffingSchrödingerSolution[
  {0.326, 1, 0.929, 18.896, 1.767, 1.420},  
  {-1.503, 5.397, 0.476}, {7, 6, 3}, 1200])[[1]]], 
            {x, -7, 7}, {t, 0, 6}, 
            PlotRange -> All, Mesh -> False, PlotPoints -> 600];

To calculate the Wigner functions, we do not use directly its 2+1D differential equation, but rather use ψ(x,t) from above and perform a numerical Fourier transform. Given a wave functions ψ(x,t) for a fixed t, the function ΨToW generates an array of dimension d×d that approximates the Wigner functions W(t,x,p) on an equidistant grid in [-L,L]×[-P,P]. To perform the numerical Fourier transform sufficiently accurate, we use d λ points for the spatial discretization. With the notations (Subscript[y, r])=-2L+r Δy, (Subscript[p, s])=-2(Subscript[P, h])+s Δp, r=1,…,n and Δy=4L/n, Δp=4(Subscript[P, h])/n, (Subscript[P, h])=P/2, we have the following approximation of the Fourier transform through a discrete fractional Fourier transform (Subscript[ℱℱ𝒯, {α,β}])[(Subscript[{(Subscript[f, r])}, r=1,…,n])].

F(p)=(∫_(-∞)^∞)exp(i((p )/ℏ)y)f(y)dy≈(∫_(-2L)^(2L))exp(i((p )/ℏ)y)f(y)dy≈Δy(∑_(r=1)^n)exp(i((p )/ℏ)(Subscript[y, r]))f((Subscript[y, r]))
F((Subscript[p, s]))=Δy (Subscript[A, s])(Subscript[ℱℱ𝒯, {1,((4 L P)/(n π ℏ))}])[(Subscript[{(Subscript[B, r])f((Subscript[y, r]))}, r=1,…,n])]

Here (Subscript[A, s])=exp(i/ℏ( (4 L (Subscript[P, h])-Δp Δy)- Δp (2 L-Δy)s)) and (Subscript[B, r])=exp(-i/ℏ (2 (Subscript[P, h])-Δp) Δy r).

Off[Fourier::fpopt2];

ΨToW[Ψ_, {L_, P_}, ℏ_, d_, λ_] :=
Module[{𝒹 = λ d, Ph, Δy, Δp, b, As, Br, WMat, ΨL, ΨΨC, 
        tpa = Developer`ToPackedArray},
 (* fractional Fourier transform parameters *)
 Ph = P/2; Δy = 4L/𝒹; Δp = 4Ph/𝒹; b = 8. L Ph/(𝒹 π ℏ);
 Br = tpa[Table[Exp[-I/ℏ (2. Ph - Δp) Δy r], {r, 𝒹}]];
 As = tpa[Δy Table[Exp[I/ℏ ((4. L Ph - Δp Δy) - Δp (2 L - Δy) s)], {s, 𝒹}]];
 (* d×d matrix of Wigner function values *)
 Table[(* Fourier transform integrand *)
       ΨL = tpa[Table[Ψ[N[x - y/2]], {y, -2L + 4L/𝒹, 2L, 4L/𝒹}]];
       ΨΨC = ΨL Conjugate[Reverse[ΨL]];
       (* Wigner transform for fixed x *) 
       Wx = As Fourier[Br ΨΨC, FourierParameters -> {1., b}];    
       (* average over λ p-values *)
       1/(2.Pi ℏ) Plus @@@ Partition[Wx, λ],             
       {x, -L, L, 2L/(d - 1)}] // Transpose]

To display the Wigner function, we use a density plot. The Wigner function is purely real, and we color positive parts in blue and negative parts in red [✶]. The function rbColor yields the color value for 0<xM<x<xP>0. For both the quantum and the classical distribution, we will use a rescaled version to obtain a better contrast.

rbColor[x_, {xM_, xP_}] := 
With[{ξ = x/Max[Abs[{xM, xP}]]},
(* negative values in red, positive values in blue; 
   scale with largest positive or negative value *)
RGBColor[Min[If[x <= 0, 1, 1 - ξ], 1], 
         Min[If[x <= 0, 1 - Abs[ξ], 1 - Abs[ξ]], 1], 
         Min[If[x >= 0, 1, -ξ + 1], 1]]]
         
(* use a scaling to show details better *)
scaleW[W_] := Sign[W] Abs[W]^(1/4);

DuffingWignerGraphics[Ψ_, ℏ_, {L_, P_}, {d_, α_}, opts___] := 
Module[{WMat, xM, xP, WScaled},
  (* form Wigner distribution from wave function *)
  WMat = ΨToW[Ψ, {L, P}, ℏ, d, α];
  {xM, xP} = {Min[#], Max[#]}&[WScaled = scaleW[Re[WMat]]];
  (*display colored scaled Wigner distribution *)
  ListDensityPlot[WScaled, Mesh -> False, PlotRange -> All, 
                  MeshRange -> {{-L, L}, {-P, P}}, opts,
                  ColorFunction -> (rbColor[#, {xM, xP}]&),
                  ColorFunctionScaling -> False]]

Later, we will show a side-by-side animation of the time-evolution of the Wigner function and the classical distribution function. For now, we show a detailed plot of the Wigner function after three periods of the external force. The intermingled domains of red and blue show that the Wigner functions acquired a complicated inner structure already after thee periods of the driving force [✶], [✶], [✶], [✶], [✶], [✶], [✶], [✶], [✶], [✶]. The structures exist at a scale smaller than ℏ [✶], [✶], [✶], [✶], [✶], [✶].

Ψ[x_] = If[Abs[x] > L, 0. + 0. I, Evaluate[dsNDSol[[1, 1, 2]][x, 3 T]]];

Show[GraphicsArray[
Block[{$DisplayFunction = Identity},
(* density and contour W == 0 graphics *)
{#, ListContourPlot[(* center part only *)
                    Map[If[Abs[#] < 0.02, 0, #]&, #[[1]], {2}],
                    Contours -> {0}, ContourShading -> False,
                    FrameTicks -> None]}&[
(* the Wigner function density plot *)
DuffingWignerGraphics[Ψ, ℏ, {8, 18}, {600, 8}, FrameTicks -> None]]]]]

Now, we come to the classical counterpart. Here are the equations of motions for the Duffing oscillator.

NewtonDuffingEquations[{m_, α_, β_, γ_, ω_}, {x_, p_}, 
                       {x0s_, p0s_, t0_}] :=
  {x'[t] == p[t]/m, p'[t] == m (2 β x[t] - 4 α x[t]^3 + γ Cos[ω t]), 
   x[t0] == x0s, p[t0] == p0s}

The overall complexity of the motion of the Duffing oscillator becomes visible by following a set of particles with initial conditions on ellipses around the origin in phase space over only 1(1/2) periods. The next graphic shows the resulting shapes of the initial ellipses [✶].

Module[{o = 24000, 𝓅 = (* 3/2 periods only *) 3/2, x0s, p0s, ndsol},
Show[Graphics[
Table[(* solve equations of motion starting with 
         initial conditions on an ellipse around {x, p} = {0, 0} *)
 {x0s, p0s} = Transpose[Table[μ {3. Cos[φ], 6. m Sin[φ]}, 
                              {φ, 0, 2Pi, 2Pi/o}]];               
 ndsol = NDSolve[NewtonDuffingEquations[{m, α, β, γ, ω}, 
                                        {x, p}, {x0s, p0s, 0}] , 
                 {x, p}, {t, 𝓅 T, 𝓅 T}, 
                 PrecisionGoal -> 8, AccuracyGoal -> 8]; 
 {Hue[0.78 (μ - 1/4)/(7/4 - 1/8)], 
  (* polygon bounding the evolved ellipse *)
  Polygon[Transpose[{x[𝓅 T], p[𝓅 T]} /. ndsol[[1]]]]},
{μ, 7/4, 1/8, -1/8}]], AspectRatio -> 1]]

We will use matrices of size 256×256 for the classical and the quantum Wigner distribution. Larger ones are of course possible, but it will take longer to calculate their values. We will generate ten frames per period T, which allows for a reasonable smooth animation. Again, more frames are possible for the price of larger calculation timings.

As the classical analog of the Wigner function time-evolution, we use (10^6) particles with initial positions and momenta according to the initial Wigner distribution. The function RandomArray from the package Statistics`ContinuousDistributions`generates arrays of normal distributed numbers quite quickly.

Needs["Statistics`ContinuousDistributions`"];

{xDist, pDist} =  NormalDistribution @@@ {{x0, σx}, {p0, σp}};

To avoid the storage of the complete trajectory interpolants, we will solve the equations of motions over a time of T/8 for 1000 particles at once, extract their final positions and momenta and keep a binned version of them in the matrices M[τ] before continuing the solution of the differential equations over the next T/8 time interval. The function counterMatrix generates a matrix of dimension pp×pp of the binned position-momentum values of a list of values xpList within the square [-L,L]×[-P,P].

counterMatrix[xpList_, L_, P_, pp_] := 
Module[{xpDiscrete = Ceiling[pp ((# + {L, P})/{2L, 2P}& /@ xpList)],
        M0 = Table[0, {pp}, {pp}], i, j}, 
        Do[If[Min[xpDiscrete[[k]]] >= 1 && Max[xpDiscrete[[k]]] <= pp,
              (* bin numbers *) {j, i} = xpDiscrete[[k]]];
           M0[[i, j]] = M0[[i, j]] + 1, {k, Length[xpDiscrete]}]; M0]

Following now one million trajectories over 8 periods of the driving force takes about 4 hours on a 2GHz computer. We will follow 1000 trajectories in parallel. This avoid the overhead of repeatedly processing the equations but does not use too much memory at one time.

Module[{o = 1000, μ = 1000, q = 100, x0s, p0s, xpList,
        M0 = Developer`ToPackedArray[Table[0, {pp}, {pp}]]},
Do[(* initialize array for time τ *) M[τ] = M0, {τ, 0, tMax, Δt}];
(* bag for {x, p}-points *)
xpBag = Table[{}, {j, μ}, {τ, 0, tMax, Δt}];  
Do[(* list of initial {x, p}-points *)
   xpList = Transpose[RandomArray[#, {o}]& /@ {xDist, pDist}];                    
   (* add values for j-th batch *) 
   M[0] = M[0] + counterMatrix[xpList, L, P, pp];   
   xpBag[[j, 1]] = Take[xpList, q];         
   (* solve ode for t-interval {0, tMax} in steps of T/framesPerPeriod *)
   Do[{x0s, p0s} = Transpose[xpList];
      ndsol = NDSolve[NewtonDuffingEquations[{m, α, β, γ, ω}, {x, p}, 
                      {x0s, p0s, τ - Δt}], {x, p}, {t, τ, τ},
                      PrecisionGoal -> 8, AccuracyGoal -> 8];
      (* update matrices *)
      xpList = Transpose[{x[τ], p[τ]} /. ndsol[[1]]];
      M[τ] = M[τ] + counterMatrix[xpList, L, P, pp];
      (* collect unbinned *)
      xpBag[[j, τ/Δt + 1]] = Take[xpList, q], {τ, Δt, tMax, Δt}], 
   {j, μ}]]; // Timing

Here are the first 50 particles of each run shown in one graphic. To color indicates the time, ranging from red (t=0) to violet (t=(Subscript[t, max])=8T). The originally localized cloud of particles follows an extremely complicated dynamics with positive Lyapunov coefficient. After a few periods, the resulting distribution has a complicated, repeatedly folded-over structure.

Show[Graphics[Reverse @
  Table[{PointSize[0.002], Hue[0.78 τ/tMax], 
               Map[Point, Transpose[xpBag][[τ/Δt + 1]], {-2}]},
        {τ, Δt, tMax, Δt}]]]

Now, we have all data together to display the animation of the time-development of the classical and quantum distribution in parallel. The function classicQuantumDistributionsGraphicsArray generates and displays a graphics array of the two distribution functions for a given τ.

classicQuantumDistributionsGraphicsArray[τ_] := 
Module[{L = 8, P = 20, λ = 12, Ψ},
Show[GraphicsArray[
Block[{$DisplayFunction = Identity},
 {(* classical  phase space distribution *)
  ListDensityPlot[(* scale *) M[τ]^(1/4), Mesh -> False, PlotRange -> All,
                  MeshRange -> {{-L, L}, {-P, P}},  FrameTicks -> None,
                  (* density proportional "blueness" *)
                  ColorFunction -> (RGBColor[1 - #, 1 - #, 1]&)],
 (* Wigner distribution *)    
 (* re-interpolate for fixed τ for faster value access *)   
 ΨT = FunctionInterpolation[Evaluate[If[Abs[x] > L, 0. + 0. I, 
                                       Evaluate[dsNDSol[[1, 1, 2]][x, τ]]]], 
                            {x, -2L, 2L}, InterpolationPrecision -> 8,
                            InterpolationPoints -> 8192, 
                            InterpolationOrder -> 6];
 (ℳ[τ] = (Sign[#] #^4)&[#[[1]]]; #)& @
 DuffingWignerGraphics[ΨT, ℏ, {L, P}, {pp, λ}, FrameTicks -> None]}]]]]

Using now the τ-values for which previously a snapshot of the classical distribution had been calculated, we arrive at the following animation.

Do[If[(* show only every tenth step *) IntegerQ[τ/Δt/10], 
      classicQuantumDistributionsGraphicsArray[τ],
      Block[{$DisplayFunction = Identity}, 
            classicQuantumDistributionsGraphicsArray[τ]]], 
   {τ, 0, tMax, Δt}];

Do[classicQuantumDistributionsGraphicsArray[τ], {τ, 0, tMax, Δt}];

The animation shows how increasingly complicated structures form in the Wigner distribution. As mentioned above, the resulting structures appear at a scale much smaller than ℏ. While the phase space distributions look quite different after already a few periods, the corresponding marginal distributions for x and p are much more similar to each other. The following two graphic pairs show the marginal distributions for the position and the momentum for the classical and the quantum case.

Function[Mℳ, Show[GraphicsArray[
ListPlot3D[Table[Plus @@ #1[Mℳ[τ]], {τ, 0, tMax, Δt}], 
           Mesh -> False, PlotRange -> All, Axes -> False,
           PlotLabel -> #2, DisplayFunction -> Identity]& @@@
           {{Identity, "x"}, {Transpose, "p"}}]]] /@ {M, ℳ}

And the averages of the positions and momentums do agree even better with each other. Only for times larger than the calculated ones will they start to deviate substantially. This means that the larger times displayed are in the order of the Ehrenfest time of the system [✶], [✶], [✶], [✶], [✶], [✶].

With[{T = Range[256]},
Show[GraphicsArray[ Function[{f, l},
Graphics[{{Hue[#1], Line[MapIndexed[{#2[[1]], #1}&, #/#[[1]]&[
                          Table[T.(Plus @@ f[#2[τ]]), {τ, 0, tMax, Δt}]]]]}& @@@
                                {{0.78, M}, {0.00, ℳ}}}, 
         Frame -> True, FrameTicks -> None, PlotLabel -> l]] @@@ 
         {{Identity, "x"}, {Transpose, "p"}}]]]

We could now continue and analyze the value distribution of the Wigner function [✶], but we end here.

If we would add a decoherence term to Wigner evolution, the classical and the quantum distributions would be much more similar [✶], [✶], [✶], [✶], [✶], [✶], [✶], [✶], [✶], [✶].

For other visualizations of time-dependent Wigner functions, see [✶], [✶], [✶], [✶], [✶], [✶], [✶], [✶], [✶], [✶], [✶], [✶], and [✶]. For the use of the Wigner function in the time–frequency space for the classical trajectories, see [✶], [✶].