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.

This decomposition into four distinct entities is used by the C++ Wavepacket package. Now let us demonstrate how to build the skeleton of a propagation scheme out of that. Let us assume we have a one-dimensional grid set up in a Representation instance called "rep" and let us further assume that the Hamiltonian is just the kinetic energy (i.e., we have a free particle with mass "m"). Then we would set up the equations as follows:

:::c++
auto H = std::make_shared<CartesianKineticEnergy>(rep, 0, mass); // 0 - first degree of freedom

auto equation1 = std::make_shared<SchroedingerEquation>(H);
auto solver1 = std::make_shared<ChebychevPrimitive>(equation1);

auto equation2 = std::make_shared<CommutatorLiouvillian>(H);
auto solver2 = std::make_shared<ChebychevPrimitive>(equation2);

auto equation3 = std::make_shared<SchroedingerEquation>(H);
auto solver3 = std::make_shared<ChebychevRelax>(equation3);

auto equation4 = std::make_shared<OneSidedLiouvillian>(H);
auto solver4 = std::make_shared<ChebychevRelax>(equation4);

auto combinedEquations = std::make_shared<EquationSystem>(2);
combinedEquations->setEntry(0, 0, equation1);
combinedEquations->setEntry(1, 1, equation2);
auto combinedSolver = std::make_shared<ChebychevPrimitive>(combinedEquations);
    // note: this will actually not work, because the Chebychev solver needs the spectrum
    // and refuses to guess it, but that is a technical detail here.