10
$\begingroup$

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

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ I answered you question here: mathematica.stackexchange.com/questions/34334/… $\endgroup$
    – Peltio
    Commented Apr 30, 2015 at 12:42
  • $\begingroup$ 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. $\endgroup$
    – Peltio
    Commented Apr 30, 2015 at 17:38
  • $\begingroup$ @Peltio This method is supposed to keep the normalization of the wave function ψ[0 + n*Δt] $\endgroup$
    – user0322
    Commented Apr 30, 2015 at 22:23
  • $\begingroup$ 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. $\endgroup$
    – Jens
    Commented May 4, 2015 at 0:56
  • $\begingroup$ @Jens, you solved it nicely. Many thanks. $\endgroup$
    – user0322
    Commented May 4, 2015 at 14:51

1 Answer 1

Reset to default
16
$\begingroup$

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}];

eqn = Thread[
   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]

v0

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]

b1

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[[1]], 
 PlotRange -> {0, 1.1}, DataRange -> {0, 8}]

norm1

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.

Improve this answer
$\endgroup$
18
  • $\begingroup$ ,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. $\endgroup$
    – user0322
    Commented May 4, 2015 at 14:57
  • $\begingroup$ 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. $\endgroup$
    – Jens
    Commented May 4, 2015 at 18:26
  • $\begingroup$ 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"] $\endgroup$
    – user0322
    Commented May 4, 2015 at 18:52
  • $\begingroup$ 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. $\endgroup$
    – Jens
    Commented May 4, 2015 at 19:17
  • $\begingroup$ thank you so much! In the initial state you forgot the term: b/(2 Sqrt[[Pi]]) $\endgroup$
    – user0322
    Commented May 4, 2015 at 19:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.