The Aarseth (1985) timestep criterion is incorrectly implemented

[WalterDehnen]

Environment version 4.3.2, C code, MacOS

Describe the bug src/integrator_ias15.c line 591 reads

timescale2 = (sqrt(y2*y4)+y3) / (sqrt(y3*y5)+y4); // A85

but correct is

(equation 6 of Pham, Rein, Spiegel 2024), i.e.

timescale2 = sqrt((y2*y4+y3*y3) / (y3*y5+y4*y4));  // Aarseth (1985)

To Reproduce The bug is subtle and will not cause a crash, but may explain some behaviour seen in Fig.7 of Pham et al. (2024). Moreover, it only features in with a non-default time-step criterion for IAS15.

[hannorein]

Thanks for opening an issue on this. It is a bit subtle indeed. First, just to avoid confusing, note that all the quantities involved in this line are squared (timescale2 is the timescale squared, y3 is the third derivative squared, etc, see here). So I believe we are talking about the difference between

$$\tau = \sqrt{ \sqrt{ \frac{y^{(2)} y^{(4)}+y^{(3)}y^{(3)} }{y^{(3)}y^{(5)}+y^{(4)}y^{(4)}} } }$$

and

$$\tau = \sqrt{ \frac{\sqrt{ y^{(2)} y^{(4)}}+y^{(3)} }{\sqrt{ y^{(3)}y^{(5)}}+y^{(4)} } }$$

I agree that they are not the same, but the difference is similar to choosing a different norm, which mathematically and for any physical scaling argument shouldn't make a difference. We comment a bit on this in the paper at the end of Section 2. We've tried out various different norms as well as the two different equations above to combine the derivatives into one scalar quantity. We did not notice any systematic difference. Re-running the test shown in Fig 7 for both of the above equations shows the same: both equations give the very similar results. It is not obvious to me why one would be preferred over the other.