Numerics.TimeStepConvergence
Choosing time steps and truncating the Hamiltonian
Here, we want to discuss how you converge a wavepacket calculation with respect to the time discretization. In particular, you should understand
- how to estimate a reasonable initial time step.
- how to monitor that this time step works (reducing it of course if needed).
- what unexpected issue to look out for.
To restrict the scope, we mostly skip convergence issues with respect to the grid. So either you already converged the grid, or we just take Hamiltonian plus numerical errors from the finite grid as The System and disregard deviations from reality until they come back for revenge.
Furthermore, we will only consider time-independent systems here. Basically, if you can converge a time-independent system, you can also converge a time-dependent system: Guess a slightly smaller time step and adjust until you are done. However, introducing the required formalities for proper reasoning (cavity-dressed states etc.) would make this description a lot more complex.
Theory
Estimating a time step
As so often with knowledge-based work, there is no algorithm that gives you the required time step for a given setup. You need some experience, and some trial and error. Nevertheless, there are best practices.
First of all, note that you do not formally need to estimate a time step if you use an adaptive method, like the Chebychev method or adaptive ODE methods. These take an approximate error bound, and either select the time step or other parameters to stay below this error bound. However, you should always still do a simple estimation to ensure reasonable values. If you blindly use an adaptive ODE method, it may take a lot longer than you expect. If you blindly use the Chebychev method, it may choose broken parameters, in particular polynomials of such large order that they overflow the double-precision arithmetics.
The estimation of the time step generally differs between stepping methods, like ODE methods or the split operator method, and polynomial methods.
- Stepping methods: Let us assume for simplicity that your relevant energy range is [0, Emax]. Then the fastest dynamics are given by the high-energy oscillations with Tmin=2π/Emax. Typically, you need multiple time steps to resolve a single oscillation, with the exact number depending on the exact method and Hamiltonian. As a rule of thumb, a time step of Tmin/10 gives a good initial guess. However, you may still be in for surprises as we will discuss below for a model system, so you may want to truncate your Hamiltonian to values around Emax.
- Polynomial methods: For the Chebychev method, a good parameter is the Kosloff alpha number,
,
in terms of the spectral range ΔE and the time step Δt. For typical parameters, the alpha value is roughly the number of polynomials required, give or take some factor <2. It should generally be chosen in the range 40-60. For values much less than 40, the Chebychev method is inefficient, for values much larger than 60, you risk including numerically unstable high-order polynomials. While other polynomial methods have different approaches and do not give the same guarantees, this range of alpha values should still give a reasonable rule of thumb estimate in the absence of more specific measures. As you can see directly from the formula, you may want to truncate your Hamiltonian and reduce the spectral range ΔE to allow larger time steps.
Quality metrics
So now you have estimated a reasonable time step and want to check if the results are acceptable. How to do that?
The generic way consists of running the simulation with different time steps (or errors for adaptive methods) and comparing the results. What you typically see is that at some point they converge. In principle, you can monitor any observable that seems relevant. An unbiased metric is the distance between the time-dependent states, for example at the end of the simulation. You always have a useful distance induced by the Hilbert space's scalar product. If v1 and v2 are two solutions for different time steps, the distance is
where we have assumed normalized states, and in the last transformation dropped small constants; we are not interested in exact numbers, but in seeing the distance shrink.
At the end of the day, converging results always requires such a procedure. After all, you need at least a rough order-of-magnitude idea of your time discretization error. However, for checking on the fly that the results are ok, there exists a simpler metric.
We first note that in general the norm of a wavefunction is constant over time. This is a direct consequence of unitary time evolution. All numerical schemes formally produce approximations to the time evolution operator, and unless they are carefully crafted, these approximations are not unitary. Hence, most numerical schemes fail to preserve the norm, and we can easily identify a reasonable time step by checking that the norm of our wavepacket stays constant.
However, like all good ideas, there are certain cases where this scheme fails. In particular, as soon as we use absorbing boundary conditions, the norm ceases to be a preserved quantity, and the interpretation of the norm becomes difficult.
Example system
Introduction
As an example for how the theory works out, we take one of the Wavepacket demos as blueprint. The demo, HarmOscillator/Gaussian_1D/1, describes a one-dimensional harmonic oscillator (ω = m = 1 a.u.), where the initial state is a displaced Gaussian wavepacket.
The native time scale is the classical oscillation period, which is 2π a.u. for the given parameters. Apart from an offset of 1/2, the eigenenergies are spaced by the oscillator angular frequency ω, hence the time scale for the n-th quantum eigenstates is approximately the n-th fraction of the classical oscillation period.
We will keep the total simulation time to about one classical period. Longer times would make the simulations slightly slower, but do not contribute qualitatively new features.
Figure 1: Projection |<φn|ψ0>|2 of the initial state ψ0 onto the energy eigenstates φn as a function of the eigenstate's energy. Red crosses: 128 grid points, interval [-10,10] a.u.. Blue dots: 512 points, interval [-20,20] a.u.
The first figure shows the projection of the initial state onto the eigenstates of the harmonic oscillators, that is, their populations as red crosses. The Python script that generated this plot can be found here. The populations decay roughly exponentially down to about 10^-9 at an energy of 50 a.u. and remain on this level. This latter feature is not what we would expect from an analytic solution, instead the populations should decrease monotonically with energy. In fact, this feature is an artefact from the finite grid, which is particularly pronounced for high energies. To demonstrate this, we can use a brute force approach and double the size and density of the grid, thereby increasing the number of points fourfold. The resulting blue dotted line confirms this assumption, the populations decrease further as expected from theory.
Note that the grid is not seriously flawed; we merely have practically unpopulated high-energy states (populations of 10-9) that have such high energies that they are no longer described very well by our finite grid, but experience all the grid artefacts (reflections, periodic boundary conditions etc.). These states are basically "numerical noise" that is acceptable because it should have no significant impact on the simulation due to the low population. Such tiny artefacts always exist.
As a final remark, this system lacks some of the real-world complications for simplicity. The Hamiltonian is Hermitian and time-independent, there are neither open-system terms and no absorbing boundary conditions. None of this poses a serious challenge; usually all of these terms should be small compared to the intrinsic system dynamics, so you just reduce the time step a little more.
Choosing a time step
From figure 1, we can estimate that an energy of 50 a.u. should be enough to faithfully reproduce the relevant dynamics of our wavepacket. This energy corresponds roughly to an oscillation time of 0.1 a.u. If we further aim to represent one oscillation period by 10 time steps, we end up with a time step of 0.01 a.u.
Now let us see what happens with the norm of the wavepacket when we apply a fifth-order Runge-Kutta method with different time steps.
Figure 2: Norm of the wavepacket as a function of time for different time steps.
Figure 2 shows the norm of the wavefunction as a function of time for different time steps ( The Python script that produced the plot can be found here). As a technical side note, the actual propagation was done with an adaptive Runge-Kutta method because the C++ version offers nothing else. The normal Runge-Kutta method was emulated by setting extremely large error bounds.
For a too large time step of 0.1 a.u., the norm diverges essentially immediately, showing the immediate use of the norm as a monitoring tool. However, even for our estimate of 0.01 a.u., the norm diverges rapidly at some point, while a smaller time step of 0.005 a.u. shows no divergence.
So, what has happened here?
| dt = 0.01 a.u. | dt = 0.005 a.u. |
|---|---|
 |
 |
Figure 3: *Projection of the time-dependent state onto eigenstates of the Hamiltonian, similar to figure 1, for different times. Left plot: dt = 0.01. Right plot: dt = 0.005.
To gain some further insight, we calculate spectra similar to figure 1, but for different times. This is done in figure 3 for the two relevant time steps of 0.01 a.u and 0.005 a.u. (Python script here).
These spectra tell a story. Apparently, the numerical errors that lead to the norm deviation accumulate not in the energy range where most of the wavepacket resides (< 50 a.u.), but at the highest energies. Now if we have an eigenstate for a given energy, we try to approximate the correct time evolution, exp(-iEt), by discrete time steps, which gives an error. An example is shown below for the exceedingly bad but simple Euler method. While the fifth-order Runge-Kutta method is way more accurate, its problems are qualitatively similar. In particular, you need to sample one oscillation by multiple steps, say 5 or so, otherwise, you cannot hope to describe the dynamics accurately.
Figure 4: Approximation of the correct time evolution in the complex plane c(t) = exp(-iEt) (black circle) by the Euler method for different time steps as fractions of the oscillation time T = 2π/E.
Now for the highest energy states, we have an energy E of about 200 a.u., hence an oscillation period of about 0.03 a.u. Such an oscillation is poorly sampled by the time step of 0.01 a.u., leading to a rather large error in the norm. For the half as large time step, we sample an oscillation with about 6 discrete time steps, which is just enough for reasonably accurate results. Note that if you look closely at figure 3, the highest energies also show a norm increase for dt=0.005 a.u., but much smaller, maybe one order of magnitude every 10 a.u. Hence, if we would simulate much longer, say, for 100 a.u., the norm would also diverge for the smaller time step.
Importance of truncation
The main result from the application to the model system is that we must always choose our time step suitable for the whole spectrum of our system, not just for the dynamics that we consider relevant.
This observation has an important consequence: If you are interested in large time steps for performance, you should always truncate the spectrum of your Hamiltonian or its components! There are two reasons why truncating operators should not bother you too much:
- On the purely numerical side, these high-energy states are often garbage anyway. As an example, have a look at figure 1 how essentially all states above 50 a.u. are incorrect. Hence, there is little harm done if you screw them up completely by truncating the operators.
-
On the physics side, you need to recall that most systems are approximations. For our harmonic oscillator example, we observed severe numerical errors at energies larger than E = 50 a.u., which corresponds to the fiftieth harmonic oscillator eigenstate. Recall that harmonic oscillators are low-order Taylor expansions of the true potential. So whatever the real system is, we can safely assume that, say, the fiftieth excited state is not going to look like a harmonic oscillator eigenstate.
No matter how accurately you model your system, you will always have such severe approximations. Maybe you use the Born-Oppenheimer approximation and ignore non-adiabatic couplings: This approximation fails for large energies. Maybe you use a model for your thermal bath: Rest assured that your model is not quantitatively correct at high energies.
So: always know and pick the energy range that you consider important. Then be brave and truncate the spectrum of your operators (typically kinetic energy and potential) outside of this range, knowing full well that you limit the validity of your calculations to sufficiently low energies. What you loose in correctness, you will definitely gain in computation speed and numerical robustness.
Figure 5: *Norm as a function of time for different time steps. The kinetic and potential energy have been truncated at 50 a.u. each, that is, values larger than 50 a.u. have been set to exactly 50 a.u. in the corresponding DVR/FBR grid representations.
For the model system, figure 5 shows the effect of truncating the kinetic energy and potential to 50 a.u. each (Python script here).
We can clearly see that the norm no longer diverges for the time step of 0.01 a.u., giving us a factor of 2 decrease in computation cost.
In practice we might truncate at considerably lower energies (20 a.u. or so); here we chose 50 a.u. because that is where the obvious numerical artifacts start.
