Solving a time-dependent Schrödinger equation

Question

I want to solve the time-dependent Schrödinger equation:

$$ i\partial_t \psi(t) = H(t)\psi(t)$$

for matrix, time-dependent $H(t)$ and vector $\psi$.

What is an efficient way of doing this so that it efficiently scales to high-dimensional spaces?

Answer 1

Time-dependent case

in the time-dependent case, $[H(t),H(t')]\neq0$ in general and we need to time-order, ie, the operator taking a state from $t=0$ to $t=\tau$ is $U(0,\tau)=\mathcal{T}\exp(-i\int_0^\tau dt\, H(t))$ with $\mathcal{T}$ the time-ordering operator. In practice we just split the time interval into lots of small pieces (basically using the Baker-Campbell-Hausdorff thing).

So, consider the time-dependent Hamiltonian for a two-level system:

$$ H = \left( \begin{array}{cc} \epsilon_1 && b \cos(\omega t) \\ b\cos(\omega t) && \epsilon_2 \\ \end{array} \right) $$

i.e. two level coupled by a time-periodic driving (see here). Even this simplest possible periodically-driven system can't be solved analytically in general.

Anyway, here's a function to construct the hamiltonian:

ham[e1_, e2_, b_, omega_, 
  t_] := {{e1, b*Cos[omega*t]}, {b*Cos[omega*t], e2}}

and here's one to construct the propagator from some initial time to some final time, given a function to construct the Hamiltonian matrix at each point in time (and splitting the interval into $n$ slices--you should try with increasing $n$ until your results stop changing):

ClearAll@constructU;
constructU::usage = "constructU[h,tinit,tfinal,n]";
constructU[h_, tinit_, tfinal_, n_] := 
 Module[{dt = N[(tfinal - tinit)/n], 
   curVal = IdentityMatrix[Length@h[0]]}, 
  Do[curVal = MatrixExp[-I*h[t]*dt].curVal, {t, tinit, tfinal - dt, 
    dt}];
  curVal]

This constructs the operator $U(0,\tau)=\mathcal{T}\exp(-i\int_0^\tau dt\,H(t))$ as $$ U(0,\tau)\approx\prod_{n=0}^{N}\exp\left( -iH(ndt)dt \right) $$ with $N=\tau/dt-1$ (or its ceiling anyway). This is an approximation to the correct $U$.

And now here is how to look at the time-dependent expectation of $\sigma_z$ for different coupling strengths $b$:

ClearAll[cU, psi0];
psi0 = {1., 0};
Manipulate[
 ListPlot[
  Table[
   Chop[#\[Conjugate].PauliMatrix[3].#] &@(constructU[
       ham[-1., 1., b, 1., #] &, 0, upt, 100].psi0),
   {upt, .01, 20, .1}
   ],
  Joined -> True,
  PlotRange -> {-1, 1}
  ],
 {b, 0, 2}
 ]

Alternatively, you could calculate the wavefunction at some time tfinal given the wavefunction at time tinit with this:

propPsi[h_, psi0_, tinit_, tfinal_, n_] := 
 Module[{dt = N[(tfinal - tinit)/n],
   psi = psi0},
  Do[
   psi = MatrixExp[-I*h[t]*dt, psi], {t, tinit, tfinal - dt, dt}
   ];
  psi]

which uses the form MatrixExp[-I*h*t,v]. For large sparse matrices (eg, for h a many-body Hamiltonian), this can be much faster at the cost of losing access to $U$.

Answer 2

Since there hasn't been any discussion of NDSOlve yet, let me point out that for a finite-dimensional Hilbert space where the Schrödinger equation is merely a first-order equation in time, it's easiest to just do this (using the two-dimensional Hamiltonian ham from acl's answer):

ham[e1_, e2_, b_, omega_, 
   t_] := {{e1, b*Cos[omega*t]}, {b*Cos[omega*t], e2}};
Manipulate[Module[{ψ, sol, tMax = 20},
  sol = First@NDSolve[{I D[ψ[t], t] == 
       ham[-1, 1, b, 1, t] .ψ[t], ψ[0] == {1,0}}, ψ, {t, 0, tMax}];
  Plot[Chop[#\[Conjugate].PauliMatrix[3].#] &@(ψ /. sol)[t], 
     {t, 0, tMax}, PlotRange -> {-1, 1}]
  ],
 {{b, 1}, 0, 2}
 ]

I copied the parameters from acl's answer too, to show the direct comparison in the Manipulate. Here the vector $\psi$ is recognized by NDSolve as two-dimensional, so the formulation of the problem is quite concise, and we can leave the time step choice up to Mathematica instead of choosing a discretization ourselves.

Answer 3

Frame it as a set of Linear ODEs and solve it somehow. I usually use Implicit Runge Kutta in the interaction picture.

solver[H_, a_] := 
soln = Module[{d, init, eq, vars, solargs, t, t0, tf}, 
d = Dimensions[H][[1]];
t0 = a[[2]];
tf = a[[3]];
t = a[[1]];
u[t_] := Table[Subscript[u, i, j][t], {i, 1, d}, {j, 1, d}];
init = (u[t0] == IdentityMatrix[d]);
eq = (I u'[t] == H.u[t]);
vars = Flatten[Table[Subscript[u, i, j], {i, 1, d}, {j, 1, d}]];
solargs = LogicalExpand[eq && init];
Return[
 NDSolve[solargs, vars, a, 
  Method -> {"FixedStep", 
    Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 10}}, 
  StartingStepSize -> tf/100, MaxSteps -> Infinity]]];

U[t_] := u[t] /. soln[[1]]

Alternatively, you could solve for the $\psi(t)$ and obtain U as $|\psi(t)\rangle\langle \psi(0)| $ and use an appropriate normalisation to preserve probability.