12
$\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$ 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$ 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$ 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$ Commented May 4, 2015 at 0:56
  • $\begingroup$ @Jens, you solved it nicely. Many thanks. $\endgroup$ Commented May 4, 2015 at 14:51

1 Answer 1

Reset to default
18
$\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$
20
  • $\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$ 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$ 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$ 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$ Commented May 4, 2015 at 19:17
  • 1
    $\begingroup$ @ayr For this differential equation, you can't specify spatial boundary conditions (that's what you'd do in the time-independent Schrödinger equation which is of elliptical type). Instead, the entire solution is determined by the initial condition - no spatial boundaries at all. $\endgroup$ Commented Mar 12 at 23:47

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.