7.5 Variable Time Steps: Better Collision Criteria

If we sprinkle stars into space, roughly in the shape of a star cluster, with random positions and velocities, the time step criterion which we just introduced works fine. The chance that the relative velocities are all zero or extremely small is itself extremely small. However, during a long run, with millions of time steps, highly unlikely cases are bound to occur occasionally, and we should worry about them. In addition, it may happen that we start from initial conditions in which all velocities are zero. This is termed a `cold collapse' type of initial condition, since the particle `temperature', measured by their kinetic energies, is zero, and the particles will start of by falling roughly towards the center of the particle distribution. In such a case, hermite5.C would crash right at the beginning, during the attempt to calculate an infinitely long time step by dividing by zero.

We can make our code more safe, and able to handle cold starts, by adding an extra criterion. In addition to the time scale $t_{1ij} = \vert r_{ij}\vert/\vert v_{ij}\vert$, we introduce a second time scale $t_{2ij} = \sqrt{\vert r_{ij}\vert/\vert da_{ij}\vert}$, where $da_{ij}$ is the relative acceleration between particles $i$ and $j$, due to their mutual attraction. Since acceleration is the second derivative of position, we have to take a square root in order to wind up with a quantity that has the dimension of time.

Note the notation used here: in analogy with the definitions of $r_{ij} = \vert{\bf r}_i - {\bf r}_j\vert$ and $v_{ij} = \vert{\bf v}_i - {\bf v}_j\vert$, we could introduce $a_{ij} = \vert{\bf a}_i - {\bf a}_j\vert$. However, this last quantity would denote the difference between the total accelerations felt by particles $i$ and $j$ due to all other particles in the system. In contrast, we are only interested in that part of the difference in their accelerations due to their own accelerations, $da_{ij}$, since it is that part that is most likely to be important for neighboring particles. In other words, we neglect the gravitational tidal field contributions of all other particles while examining a particle pair. The alternative, of using $a_{ij}$ instead of $da_{ij}$, would be far more expensive. We would have to compute all total accelerations first, and then do a pairwise comparison, which would entail a double pass over all particles, each with a computational cost of order $N$.

Here is the implementation of hermite6.C. With respect to hermite5.C, the differences are all limited to the first function, acc_and_jerk.

Exercise 7.2 (Pythagorean problem: variable time steps)
Repeat the previous exercise, starting again with the Pythagorean initial conditions, but now use hermite6.C. In order to get good energy conservation, roughly how much time do you gain when using a variable integrator?