# Animated Wave propagation using split operator method

Here is the state of the art: By solving the time-dependent Schrödinger equation, one obtains that the propagated wave function after time step $\Delta t$ can be calculated by applying the time-evolution operator on the wave packet at any instant $t$ \begin{equation} \Psi({\vec r},t+\Delta t)=e^{-i H\Delta t}\Psi({\vec r},t), \end{equation} where we assumed that the Hamiltonian $H$ is time-independent.

Following the split-operator method, the time evolution operator for the Hamiltonian can be calculated as a series of matrix multiplications \begin{equation} \Psi({\vec r},t+\Delta t)=M1.M2.M1\Psi({\vec r},t)+O(\Delta t^3), \end{equation}

I would like to do iteratively:I start with psi(r,0) to calculate psi(r,Δt); then I use the psi(r,Δt), multiply it by the matrices and obtain psi(r,Δt+Δt); then I take psi(r,Δt+Δt), multiply it by the matrices and obtain psi(r,Δt+Δt+Δt), and so on.

 ass = {Δt>0, v > 0, b > 0, \hbar > 0, kx >= 0,ky >= 0};

η[ℏ_, b_, x_, y_] =b/(2 Sqrt[π]) (x + I*y)*E^(-b*(x^2 + y^2)/(4*ℏ) )*{{1}, {1}};

M1[Δt_,v_,b_,y_,ℏ_] = {{1, 0}, {0, 1}}*Cos[(Δt*v*b*y)/(2 ℏ)]-
Sin[(Δt*v*b*y)/(2 ℏ)]/((Δt*v*b*y)/(2 ℏ))*{{0, -b y}, {-b y, 0}};

M2[Δt_,v_,b_,y_,ℏ_,kx_,ky_] = {{1, 0}, {0, 1}}*Cos[(Δt*v Sqrt[kx^2 + ky^2])/(2 ℏ)]-Sin[Δt*v Sqrt[kx^2 + ky^2]]/(Δt*v Sqrt[kx^2 + ky^2])*{{0, kx - I*ky}, {kx + I*ky, 0}};


The Time Evolution of a Wave Packet can be calculated as a series of matrix multiplication

Ψ[Δt_,v_,b_,x_,y_,ℏ_, kx_,ky_] := M1[Δt, v, b,y, ℏ].M2[Δt, v, b, y, ℏ, kx,ky].M1[Δt, v, b, y, ℏ].η[ℏ, b,x, y];


This is the same expression as in Eq. (16) here

• I answered you question here: mathematica.stackexchange.com/questions/34334/… – Peltio Apr 30 '15 at 12:42
• I used the name for the general state, instead of the specific instance for your example. I have fixed it, and given you a hint on how to plot the evolution of the state. If only this question was reopened... (and frankly I don't understand why it was put on hold at all - apart for not using mathematica syntax, it was pretty clear what you wanted). I hope, once the question is reopened, you will edit it to provide meaningful values for the parameters and add a picture of the resulting plot. – Peltio Apr 30 '15 at 17:38
• @Peltio This method is supposed to keep the normalization of the wave function ψ[0 + n*Δt] – user0322 Apr 30 '15 at 22:23
• The formulation of the problem is incorrect. The matrices are infinite-dimensional either in real or Fourier space and need truncation. You can't get the time evolution of a spinor with position dependence by just using 2D matrix multiplication. The main issue here is unrelated to Mathematica. – Jens May 4 '15 at 0:56
• @Jens, you solved it nicely. Many thanks. – user0322 May 4 '15 at 14:51

The goal of the question is to solve a system of coupled linear differential equations representing a two-component wave function depending on two spatial coordinates x, y, and one time variable t. The system is governed by a time-independent Hamiltonian that acts on the two components and the spatial coordinates.

The question contains a matrix version of the time-evolution operator that is two-dimensional, i.e., acts only on the wave function components. It is translated from a paper in which the treatment of the spatial degrees of freedom is not written down explicitly. This caused the misunderstanding that the evolution across a time step can be done by a multiplication with the two-dimensional matrix. In practice, what you really have to do in order to implement this formalism is to form a KroneckerProduct between the matrix form you have, and the matrices representing the discretization of the position-dependence, or equivalently its Fourier transform.

But implementing this from scratch isn't really necessary, because the built-in NDSolve is able to handle this kind of problem, too. Here is an example:

Clear[x, y, ψ1, ψ2, b, eqn, eqnWithInitial, aX, aY, v];

{σx, σy} = 1/2 PauliMatrix[{1, 2}];

I D[{ψ1[x, y, t], ψ2[x, y, t]},
t] == σx.(-I D[{ψ1[x, y, t], ψ2[x, y, t]}, x] +
aX {ψ1[x, y, t], ψ2[x, y, t]}) +
σy.(-I D[{ψ1[x, y, t], ψ2[x, y, t]}, y] +
aY {ψ1[x, y, t], ψ2[x, y, t]}) +
v[x, y] {ψ1[x, y, t], ψ2[x, y, t]}];

eqnWithInitial =
Join[eqn,
Thread[{ψ1[x, y, 0], ψ2[x, y, 0]} == {1,
1} (x + I*y) Exp[-(x^2 + y^2)]],
Thread[{ψ1[-5, y, t], ψ2[-5, y, t]} == {ψ1[5, y,
t], ψ2[5, y, t]}],
Thread[{ψ1[x, -5, t], ψ2[x, -5, t]} == {ψ1[x, 5,
t], ψ2[x, 5, t]}]];

v[x_, y_] := 0;
b = 0;
{aX, aY} = {-b y, 0};
tMax = 8;
solution = First@NDSolve[
eqnWithInitial, {ψ1[x, y, t], ψ2[x, y, t]}, {x, -5,
5}, {y, -5, 5}, {t, 0, tMax},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"DifferenceOrder" -> "Pseudospectral"}}];
Ψ1[x_, y_, t_] = ψ1[x, y, t] /. solution;
Ψ2[x_, y_, t_] = ψ2[x, y, t] /. solution;

pl = Table[
Plot3D[{Re[Ψ1[x, y, t]] - 2,
2 + Re[Ψ2[x, y, t]]}, {x, -5, 5}, {y, -5, 5},
PlotRange -> {-3, 3},
PlotStyle -> {Red, Directive[Opacity[.9], Orange]},
BoxRatios -> 1], {t, 0, tMax, tMax/20}];

ListAnimate[pl] This is the time evolution without potential and without magnetic field. Shown are the two components of the wave function in a plot that represents their real parts versus x and y, offset for clarity.

Next, turn on a magnetic field:

b = 1;
{aX, aY} = {-b y, 0};
solution = First@NDSolve[
eqnWithInitial, {ψ1[x, y, t], ψ2[x, y, t]}, {x, -5,
5}, {y, -5, 5}, {t, 0, tMax},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"DifferenceOrder" -> "Pseudospectral"}}];
Ψ1[x_, y_, t_] = ψ1[x, y, t] /. solution;
Ψ2[x_, y_, t_] = ψ2[x, y, t] /. solution;

pl = Table[
Plot3D[{Re[Ψ1[x, y, t]] - 2,
2 + Re[Ψ2[x, y, t]]}, {x, -5, 5}, {y, -5, 5},
PlotRange -> {-3, 3},
PlotStyle -> {Red, Directive[Opacity[.9], Orange]},
BoxRatios -> 1], {t, 0, tMax, tMax/20}];

ListAnimate[pl] In the NDSolve command, I use a set of differential equations created from the Dirac Hamiltonian in eq. (12) of the paper linked in the question. This way, there is no need to re-invent the wheel by implementing our own time-stepping algorithm. I'm leaving that work to Mathematica. In order to make this work, I chose periodic boundary conditions in space, and a Gaussian initial condition with orbital angular momentum 1 (as it also appears in the question) at t = 0. All inessential parameters were given some arbitrary values. The main parameters are the magnetic field b which enters the vector potential ax, ay through the Landau gauge, and the potential v[x, y]. I set the latter to zero because it doesn't appear in the question at all. But you could define it to be a nonzero function of position, as long as it's consistent with the periodic boundary conditions.

Edit: checking probability conservation

In response to the comment, I would suggest something like this if you want to check whether the numerical integration has been able to preserve the norm of the solution, as it should for a hermitian Hamiltonian:

gridNormsB = With[{delta = .1},
Table[
Total@Flatten@
Table[
Abs[Ψ1[x, y, t]]^2 + Abs[Ψ2[x, y, t]]^2,
{x, -5, 5, delta}, {y, -5, 5, delta}],
{t, 0, tMax, 1}]]

(*
==> {156.28745, 156.2309, 156.23151, 156.18101, 156.1363, 156.12584, 156.08589, 156.00673, 155.99424}
*)

ListLinePlot[gridNormsB/gridNormsB[],
PlotRange -> {0, 1.1}, DataRange -> {0, 8}] Here I did not do the norm using NIntegrate because that actually crashed in Mathematica version 8 (maybe I ran out of memory - anyway, the following is faster). Instead, I chose a discrete grid of spatial points and simply summed the quantity Abs[Ψ1[x, y, t]]^2 + Abs[Ψ2[x, y, t]]^2 over these points, varying t from the initial t=0 to tMax = 8, using the solutions of the last calculation (for b = 1). As you can see, the norm is indeed preserved.

• ,it is a classic solution. This is what I wanted. But I have some unclear points like 'Re[Ψ1[x, y, t]] - 2, 2 + Re[Ψ2[x, y, t]]'. I don't understand the term 2. And I would like to plot the integral of the squared absolute value of the wavefunction over space (The probablity density function |ψ1[x, y, t]|^2+|ψ2[x, y, t]|^2 at any time t=5). Later on I will calculate the expectation value of the angular momentum. – user0322 May 4 '15 at 14:57
• The 2 is added and subtracted only to offset the 3D plots vertically, so you can see them better. You can change the plot to Plot3D[ Abs[\[CapitalPsi]1[x,y,t]]^2+Abs[\[CapitalPsi]2[x,y,t]]^2, ...  and that will produce the modulus squared. For more information on the styling, also have a look at the options of the Plot3D command. – Jens May 4 '15 at 18:26
• I tryed this but couldn't give the plot Plot3D[Integrate[ Abs[[CapitalPsi]1[x, y, 1]]^2 + Abs[[CapitalPsi]2[x, y, 1]]^2, {x, -10, 10}, {y, -10, 10}], ColorFunction -> "Rainbow"] – user0322 May 4 '15 at 18:52
• You can't make a plot of that because there are no more variables to vary when you integrate. Also, the integration has to be within the NDSolve domain. – Jens May 4 '15 at 19:17
• thank you so much! In the initial state you forgot the term: b/(2 Sqrt[[Pi]]) – user0322 May 4 '15 at 19:22