I need to numerically evolve a density matrix using this formula(Actually I have more terms but right nows I am starting with this and facing problems):

$$\dot\rho(t) = -i[H(t), \rho(t)]$$

$H(t)$ is the Hamiltonian which is of a form $H_o + \epsilon(t)H_1$ where both $H_o$ and $H_1$ are constant $20\times20$ matrices and $\epsilon(t)$ is a piece-wise constant function.

So, my numerical scheme is basically

$$\rho(t + dt) = \rho(t) -i[H(t), \rho(t)]dt$$

and substituting $\rho(t + dt)$ to the right hand side to get $\rho(t + 2dt)$ and so on.

The problem I am facing is that the entries of the density matrix quickly blow up.

I have tried checking whether the density matrix remains a valid density matrix at each iteration: namely if it is Hermitian, if it's trace is one and if it's eigenvalues sum to one.

Now the first two conditions are always met but it's eignevalues are becoming larger and larger. In some sense it seems reasonable too as I could have changed the time-step $dt$ to say 2 times it's value and the eigenvalues of $\rho(t+dt)$ would change. So is it necessary to employ some sort of normalization procedure to the density matrix after each iteration?

Note: I have tried to use very small values of $dt$ and there is no improvement so the problem is not in the value of the time-step.

  • $\begingroup$ The trace can't be zero and the eigenvalues sum to one, it's the same thing... $\endgroup$
    – kηives
    Feb 26, 2018 at 15:36
  • $\begingroup$ Sorry I meant trace being one. And then the trace does not equal the sum of eigenvalues if the matrix is not triangular. $\endgroup$ Feb 28, 2018 at 5:52

1 Answer 1


As always with numerical problems, it's a good idea to simplify things as much as possible to isolate what's really going wrong. You can simplify the problem by considering the evolution of a pure state $\psi$ instead of an entire density matrix, so solving

$$i\dot\psi = H\psi$$

for the vector $\psi$ instead of the matrix $\rho$. Additionally, you can forget the perturbation $\epsilon$ and just assume that $H = H_0$.

The numerical method you're using is called the explicit Euler scheme. The explicit Euler scheme is guaranteed to give a sequence of solutions that, for any instant in time, will converge to the true solution in the limit as $\delta t$ goes to 0. However, it is not guaranteed to preserve every property that might be useful to you, for example the fact that the probability amplitude should be conserved.

To see this in more detail, let's compare the exact solution $\Psi$ with the explicit Euler solution $\psi$. At time $t + \delta t$, you can write down an expression for $\Psi$ using the matrix exponential:

$$\Psi(t + \delta t) = e^{-i\delta t\cdot H}\Psi(t).$$

Since the Hamiltonian $H$ is Hermitian, the operator $e^{-i\delta t\cdot H}$ is unitary. One consequence is that the norm of $\Psi(t + \delta t)$ is the same as the norm of $\Psi(t)$, so the probability amplitude is conserved for the exact system. With the explicit Euler approximation, however, we can write the value $\psi_{n + 1}$ of the wave function at time $(n + 1)\delta t$ in terms of the previous value $\psi_n$ as:

$$\psi_{n + 1} = \psi_n - i\delta t H \psi_n = (I - i\delta t\cdot H)\psi_n.$$

The operator $I - i\delta t\cdot H$ is not unitary, so you can only expect conservation of probability in the limit as $\delta t \to 0$ and $n \to \infty$ such that $n\delta t$ stays fixed. This is really unfortunate because our approximation scheme doesn't preserve one of the most important properties of the physical system.

Luckily, there are many more timestepping schemes besides explicit Euler. Another approach might be to approximate $e^{-i\delta tH}$ by the 1, 1-Pade approximant:

$$e^{-i\delta tH} \approx \left(1 + \frac{i\delta tH}{2}\right)^{-1}\left(1 - \frac{i\delta tH}{2}\right).$$

You can show with a little linear algebra that this approximation is a unitary operator, in which case all of your discrete approximations to $\Psi$ will conserve probability up to floating-point error. This approximation is equivalent to using the Crank-Nicholson discretization of the problem instead of explicit Euler.

Leveque's book is a good introductory reference on finite difference methods. This will give you a lot of the necessary background on constructing stable discretizations of ODE and PDE. It's possible to design numerical methods explicitly for preserving things like the total probability amplitude; these are called symplectic methods. The best reference I know is Hairer, Lubich and Wanner, but unfortunately this book is quite advanced. You could start by looking at the Verlet method, which is the simplest symplectic method for Hamiltonian systems.

  • $\begingroup$ Thanks but I am curious if the Euler solution does converge to the true solution, why doesn't it have the necessary properties that the true solution is supposed to have? In theory the equation is supposed to have a unique solution which has the properties preserved. $\endgroup$ Feb 15, 2018 at 16:28
  • 2
    $\begingroup$ I've added some more to my answer. If you're just getting started with solving ODEs numerically, I highly recommend Leveque's book. You'll want to have a pretty good understanding of consistency and stability of ODE discretizations to be able to decide what is or is not a good timestepping method for a given ODE. $\endgroup$ Feb 15, 2018 at 22:41

Your Answer

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

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