Naming the components of equations

Why this blog post

I have just finished Ticket [#124], which introduced some renaming of very basic classes. And while aksing Burkhard for his opinion, he was not able to dig through my explanations. Plus, the Matlab version started a new blog, so I thought I would try to wrap up all the ideas why the new names are great in a blog post and explain along the way how a good part of the C++ version works.

The problem

Let us look at four different equations for propagating some quantum-mechanical state:

The equations are:

1. An ordinary real-time Schrödinger equation with some Hamiltonian that describes the time-evolution of a wave function.
2. An ordinary Liouville-von-Neumann (LvN) equation with the same Hamiltonian that describes the time evolution of a density operator.
3. The Schrödinger equation in imaginary time, again with the same Hamiltonian. This equation has some special uses. In particular, it can be used to relax a wave function to the Hamiltonian's ground state (whose norm decays the slowest).
4. A LvN equation in imaginary time. This equation is pretty useful in thermodynamics: If you use the unit operator at t=0, the solution of this equation is exp(-Ht), which is exactly the form of the density operator in thermal equilibrium.

The goal with the C++ library was to be able to solve all of these equations with as much reuse as possible. To get there, we first need to decompose the equations into components that we can then combine to suit all our needs. If we look hard at these equations, we notice some common patterns:

• The left-hand side of the equations can look differently. In practice, we only need two specific cases, namely the real-time propagation and the imaginary-time propagation.
  We encode this information implicitly in the choice of the propagation scheme. The reason is that real- and imaginary-time evolution behaves differently. For example, you can propagate a wave function back and forth in real time, but in imaginary time, the backwards propagation is numerically unstable with exponentially growing errors. So you choose the left-hand side of the equation by picking the solver for the differential equation.
• Obviously, we have an operator that appears in various ways in the equations. This is a component that should be available directly in the code.
  An operator must be applicable to a wave function or to the ket or bra side of a density operator.
• The operator itself is not enough. Note that the different equations use the operator in different ways. This information must be encoded somehow as well.
  Hence, we need a separate component that takes one or more operators and turns them into one of the standard equations of motion.
• Finally, we also want to build systems of equations. The background is that various treatments for open quantum systems use a form where you propagate a density operator in time plus additional auxiliary or "ghost" densities that are coupled to the "real" equation and contain the memory of the environment. Examples are the procedure by David Tannor et al. https://doi.org/10.1063/1.479669 or the hierarchical equations of motion.
  As a simple goal, we can attempt to propagate equations (1) and (2) together. While this is makes no sense, it shows that we could also propagate coupled equations of motion.