Using the functions in the “2D_TDSE_Chebyshev” notebook

(* A guide to getting started using the functions in 2D_TDSE_Chebyshev.nb.  First that notebook should be opened and all the cells in it evaluated. *)

(* Overriding values set in 2D_TDSE_Chebyshev.  Runs in seconds with 16x16, which is good for playing around. Can get high accuracy and good plots with 64x64, which might take a minute. *)
$Nx=64 ;$Ny=64;

(* Generating a Laguerre-Gaussian wavefunction with n=0, ℓ=1, and width ρ0=1. *)
Ψ0=LG[0,0,0,0,0,1,1];

(* Generating a cylindrically symmetric vector potential describing a uniform magnetic field along the z-axis.  Here the field has strength given by B0[1], which means that a Laguerre-Gaussian wavefunction of width ρ0=1 will be a Landau state (an energy eigenstate of the system).  *)
B=B0[1];
A=CylindricalA[0,0,B];

(* Settings for colour scales to be used in density plots. *)
probDensScaleMax=((Nnℓ[0,0])^2);
currentScaleMax=Max[Absjxy[Ψ0,A]];

(* Plotting probability density and probability current density.  This is done by overlaying one plot on top of the other, and note that at low resolutions these will appear misaligned.  It looks fine with 64x64. *)
Show[PlotProbDens[Ψ0,probDensScaleMax,4,4,0],PlotCurrentDens[Ψ0,A,currentScaleMax,4,4,0],ImageSize->Small]

(* Calculating time-evolution of Ψ0 over a quarter of a cyclotron period. *)
Ψ=ChebyTDSEuniformB[Ψ0,B,T[B]/4];

(* Plotting new probability density and probability current density.  Since Ψ0 is an energy eigenfunction, the probability density and current density are unchanged. *)
Show[PlotProbDens[Ψ,probDensScaleMax,4,4,0],PlotCurrentDens[Ψ,A,currentScaleMax,4,4,0],ImageSize->Small]

(* Can check the accuracy with which the norm of the wavefunction is conserved: *)
Max[Abs[Abs[Ψ]-Abs[Ψ0]]]

(* The expectation value of the kinetic orbital angular momentum can be calculated as follows.  (Note that at the beginning of the 2D_TDSE_Chebyshev notebook I set $ℏ=1.) *)
ExpectationLkinz[Ψ0,0,0,A]
ExpectationLkinz[Ψ,0,0,A]

(* For multiple time steps, it is convenient to apply ChebyTDSEuniformB iteratively using the built-in Mathematica function NestList.  Here this is done for an initial wavefunction which is the same as Ψ0 above expect it also has net momentum in the x-direction given by kx0=2.  This results in cyclotron motion of the centre of mass.  Here there are 4 propagation steps, each equal to a quarter of a cyclotron period.  *)
Ψlist=NestList[Function[Ψi,ChebyTDSEuniformB[Ψi,B,T[B]/4]] ,LG[0,0,2,0,0,1,1],4];

(* Plotting the probability density and probability current density for Ψlist.  Note that the plot area is shifted along the y-axis by an amount σ[B,2$ℏ], which is the radius of the cyclotron orbit.  (Here equal to 0.5.)  This means that the plot area is centred on the axis of the cyclotron motion.  *)
ThreeDplot=Show[PlotEvolution[Ψlist,A,probDensScaleMax,currentScaleMax,4,4,σ[B,2$ℏ]],ImageSize->Large, BaseStyle->{FontSize->18}]