  Dark matter (DM) haloes can be subject to gravothermal collapse if DM is not
collisionless but has strong self-interactions. When the scattering can
efficiently transfer heat from the centre to the outskirts, the central region
of the halo collapses and reaches densities much higher than those for
collisionless DM. This phenomenon is potentially observable in studies of
strong lensing. Current theoretical efforts are motivated by observations of
surprisingly dense substructures. A comparison with them requires accurate
predictions. One method to obtain such predictions is to use N-body
simulations. The collapsed haloes are extreme systems that pose severe
challenges to state-of-the-art codes used to model self-interacting dark matter
(SIDM). We investigate the root of such problems with a focus on energy
non-conservation. We run N-body simulations with and without DM
self-interactions of an isolated DM-only halo and change numerical parameters
relevant to the accuracy of the simulation. We find that not only the numerical
scheme for the DM self-interactions can lead to energy non-conservation but
also the modelling of gravitational interaction and the time integration are
problematic. The issues we found are: (a) particles changing their time step in
a non-time-reversible manner; (b) the asymmetry in the tree-based gravitational
force evaluation; (c) SIDM velocity kicks break the symplectic nature and time
symmetry. Tuning the parameters of the simulation to achieve a high accuracy
allows for conserving energy not only at early stages of the evolution but also
at late ones. However, the cost of the simulations becomes prohibitively large.
Some problems making the simulations of the gravothermal collapse phase
inaccurate, can be overcome by choosing appropriate numerical schemes. However,
others remain challenging. Our findings motivate further work on addressing
these challenges.
