Fast numerics for the spin orbit equation with realistic tidal dissipation and constant eccentricity

We present an algorithm for the rapid numerical integration of a time-periodic ODE with a small dissipation term that is C 1 in the velocity. Such an ODE arises as a model of spin–orbit coupling in a star/planet system, and the motivation for devising a fast algorithm for its solution comes from the desire to estimate probability of capture in various solutions, via Monte Carlo simulation: the integration times are very long, since we are interested in phenomena occurring on timescales of the order of 10 6 –10 7 years. The proposed algorithm is based on the high-order Euler method which was described in Bartuccelli et al. (Celest Mech Dyn Astron 121(3):233–260, 2015), and it requires computer algebra to set up the code for its implementation. The payoff is an overall increase in speed by a factor of about 7.5 compared to standard numerical methods. Means for accelerating the purely numerical computation are also discussed.


Introduction
Since the discovery that Mercury, uniquely in the solar system, is trapped in a 3:2 spinorbit resonance with the Sun-that is, it turns three times on its own axis for every two revolutions it makes around the Sun-there have been attempts to explain how this came about. Observations going back to the formation of the Solar system are clearly impossible, but simulation of its history is not. To be more specific, what is needed is a 'spin-orbit model', a model of the interaction between a satellite (Mercury) and the body that it orbits (the Sun). The model is in the form of a differential equation (ODE). By solving this ODE starting from many randomly chosen initial conditions, the dynamics of the satellite can then be simulated; the various solutions in which it can become captured are counted, and probabilities of capture can thereby be estimated. This is a classic example of the technique of Monte Carlo simulation, and the use of random initial conditions reflects the fact that we do not know the early history of the Solar system: hence, we choose initial conditions at random so as to cover a wide range of possibilities, and from these we estimate capture probabilities.
The context of such a calculation is wider than that of a single, albeit unusual, case in our Solar system however. Orbiting around other stars, there is evidence for the existence of exoplanets. The resonances in which these are trapped are likely to have consequences for the possibility of life having developed there. This underlines the importance of computing capture probabilities in a wider context.
Clearly, dissipation in the spin-orbit model must play an important role, since this is what leads to eventual capture in one of the available attractors. In this paper, we are concerned with a dissipation model proposed in Noyelles et al. (2014), which we refer to throughout as the NFME model. Here, the authors point out how the model that they advocate describes more realistically the effect of tidal dissipation, in comparison with the much simpler constant time lag (CTL) model. The latter has been used extensively in the literature since the seminal paper by Goldreich and Peale (1966)-see for instance (Correia and Laskar 2004;Chierchia 2008, 2009;Bartuccelli et al. 2015) for recent results on the basins of attraction of the principal attractors. A clear drawback of the CTL model is that, for the parameters of the Sun-Mercury system, the main attractor (with 70% probability of being observed) is quasi-periodic, with an approximate velocity ratio of 1.256-such a solution is usually referred to as supersynchronous in the literature. As we have seen, Mercury is in fact in a 3:2 spin-orbit resonance, and the CTL model indicates that the probability of this is only about 8%. On the other hand, with the NFME model, much more plausibly, all attractors have mean velocity ratio close to rational values, and about 42% of the initial data are captured by an attractor with mean velocity ratio close to 3/2-see Table 4. The shortcomings of the CTL model have also been stressed in Ferraz-Mello (2013, 2015, where the creep tide theory has been developed. According to such a theory, the mean excess angular velocity of the supersynchronous solution turns out to be very small (proportional to the relaxation factor) in the case of very large viscosity.
The greater realism of the NFME model comes at a cost, however: 1. The NFME model is not smooth; in fact, it is only C 1 in the angular velocity,θ , as we shall see; it also involves functions with rapid variations. 2. There is considerable detail to take into account in implementing the NFME model mathematically-see Sect. 2. 3. The functions in the model take significantly more computation time to evaluate, and so numerical calculations are comparatively slow.
The first point represents a substantial obstacle to the application of perturbation theory to the problem, as well as for implementing a fast algorithm to solve the ODE. Both of these objectives are important. For instance, perturbation theory is a useful tool for obtaining any kind of analytical results for the dynamics. For the present work, though, we should now discuss the motivation for developing fast algorithms.
It is important in Monte Carlo estimation to give bounds on the probabilities obtained. In practice, this boils down to specifying, say, the 95% confidence interval (CI) for the probabilities, and these confidence intervals shrink as I , the number of initial conditions investigated, grows-specifically, the width of a given CI decreases as I −1/2 (Walpole et al. 1998). This important result is one motivation for developing fast algorithms for the spin-orbit ODE: the more initial conditions we consider, the better defined our probability estimates become. Additionally, realistic values of dissipation in this kind of problem are very small indeed, and so capture typically takes place over a period measured in millions of years. Both of these facts motivate us to develop fast solution algorithms.
The main thrust of this paper is to show how previous work (Bartuccelli et al. 2015) on fast algorithms can be applied to solving the spin-orbit ODE with dissipation given by the NFME model. To this end, we first concisely set out the model for spin-orbit coupling in, for instance, the Sun-Mercury system, as described by Makarov (2012) and further discussed by Noyelles, Frouard, Makarov and Efroimsky in Noyelles et al. (2014). We then propose a fast numerical algorithm, based on the high-order Euler method (HEM) introduced in Bartuccelli et al. (2015), for solving the spin-orbit ODE incorporating the NFME model. Essentially, we need to deal with the time integration of a trajectory in a different way when it is close to a resonance: this requires a careful partition of the phase space into subsets that we call strips, with a different numerical integrator being used in each strip.
After rigorous testing, we conclude that the proposed algorithm is about 7.5 times faster than simply using a general-purpose numerical ODE solver-a useful gain in speed.
In most of the previous related work, for instance, Correia and Laskar (2004), Celletti and Chierchia (2008), Makarov (2012) and Makarov et al. (2012), only about 1000 initial conditions were considered, almost certainly because of computer time limitations. By contrast, we are able to consider more than 50,000 initial conditions here, with the total computation time being of order 2 weeks.
The dynamics of the solutions of the spin-orbit ODE with the NFME model are also interesting, and we discuss these in a separate paper (Bartuccelli et al. 2017).
Throughout this paper, the units for mass, length and time will be kg, km and Earth years, respectively. For ease of cross-referencing, an equation numbered (N · x) in this paper is equation (x) in Noyelles et al. (2014), modified if necessary, but only for consistency with our notation. When we refer to a 'standard numerical method' or 'numerical method' for ODE solving, we mean one of the family of well-known general-purpose ODE solvers (e.g. Runge-Kutta, Adams); when we mean HEM, which is also numerical, we refer to it explicitly by name.

The spin-orbit ODE
where θ(t) is the angle between the semimajor axis of the planet, assumed to be an ellipsoid, and the major axis of its orbit, ξ is a measure of the inhomogeneity of the planet (with ξ = 2/5 for a homogeneous sphere), M planet is the mass of the planet and R is its radius. Numerical values for these, and indeed all relevant parameters can be found in Table 1.

The triaxiality torque
The triaxiality-caused torque along the z-axis, T TRI z (θ, t), is a torque exerted on the planet by the gravitational field of the star, arising from the fact that the planet is not a perfect sphere. Goldreich and Peale (1966) and the authors of Noyelles et al. (2014) approximate T TRI z by its quadrupole part, which is given by where the principal moments of inertia of the planet, along the x, y and z axes, respectively, are A < B < C = ξ M planet R 2 ; n is its mean motion; r = r (t) is the distance between the centres of mass of the star and the planet at time t; a is the semimajor axis of the orbit of the planet; and f = f (t) is the true anomaly as seen from the star. We further define M(t), the mean anomaly, which is such that n =Ṁ, with a good approximation to n (Noyelles et al. 2014) being given by 1 n ≈ G(M planet + M star )/a 3 , where G is the gravitational constant and M star is the mass of the star. Hence, aside from a constant of integration which we set to zero, M = nt. The standard procedure from this point is to make the following pair of Fourier expansions: where X n,m k (e) are the Hansen coefficients, which depend on the orbital eccentricity e. The authors of Noyelles et al. (2014) compute them via (Duriez 2007) is the k-th order Bessel function of the first kind, and C n,m r,s = n + 1 + m r n + 1 − m s with the binomial coefficients (extended to all integer arguments) being given by For convenience, we follow Noyelles et al. (2014) by defining the functions G lpq (e) = X −l−1, l−2 p l− p+q (e). Then finally, In Noyelles et al. (2014), the sum is over q ∈ Q TRI = {−4, . . . , 6}.

The tidal torque
The authors of Noyelles et al. (2014) forcefully point out the problems associated with the constant time lag (CTL) model, which leads to a dissipation term of the form γ (θ − ω), with γ and ω constants. They argue that this model, among other things, is incompatible with the physically plausible rheology of terrestrial planets, and can also give rise to quasi-periodic solutions for the spin-orbit problem (Correia and Laskar 2004;Celletti and Chierchia 2009;Bartuccelli et al. 2015). Tidal dissipation is modelled in Noyelles et al. (2014) by expanding both the tide-raising potential of the star, and the tidal potential of the planet, as Fourier series. The Fourier modes are written ω lmpq , where l, m, p, q are integers. Appendix B of Noyelles et al. (2014) justifies the simplification l = m = 2, p = 0, which arises from the smallness of the obliquity of the planet (axial tilt), in combination with an averaging argument. This then gives the following approximation for the polar component of the tidal torque: (N.11a) Here, Q TIDE = {−1, . . . , 7}; H 2 (ω 220q ) = k 2 (ω 220q ) sin 2 (ω 220q ) , where k 2 is a positive definite, even function of the mode ω lmpq and 2 is an odd function of the mode; sgn() is the usual signum function. Hence, overall, H 2 (x) sgn(x) is an odd function of x. The reason why H 2 is written in this form is that it will turn out to depend on a fractional power of |ω lmpq |: the oddness of H 2 (x) sgn(x) is therefore retained without raising a negative number to this fractional power. The additive error terms for this expression are O(e 8 ) and O R 7 /a 7 , in which , which is O(1), is a phase lag. For Mercury, e 8 ≈ 10 −6 and (R/a) 7 ≈ 10 −31 . In Noyelles et al. (2014), the additional approximation is made that, since the terms with q = j oscillate, they naturally affect the detailed dynamics of the planet, but nonetheless, the overall capture probabilities are insensitive to them . Hence, the tidal torque is well approximated by the secular part only, i.e. the q = j terms; this gives Neglecting apsidal and nodal precessions, we have with l ≥ 2, m, p = 0, . . . , l and q ∈ Z. For the reasons mentioned earlier, we consider only the mode ω 220q = (q + 2)n − 2θ . Finally, we give the NFME model for the function H l (ω lmpq ). We use the abbreviation χ = ω lmpq , in terms of which with μ being the unrelaxed rigidity; and In these expressions, is the usual gamma function, τ A and τ M are the Andrade and Maxwell times of the mantle, respectively, and α is the Andrade parameter.

Parameter values for the Sun/Mercury system
In Table 1, we give numerical values appropriate to the Sun and Mercury for the parameters introduced so far. We standardise the units to kg for mass, km for length and (Earth) years for time.

The triaxiality and tidal angular accelerations
In order to study the spin-orbit problem further, it is convenient to introduce two quantities with the dimensions of angular acceleration Taken individually, these quantities are not accelerations per se-more correctly, they might be described as the triaxiality/tidal contributions to the angular acceleration, since, when computed on the solution θ(t), their sum isθ , as we shall see from Eq. (3)-but for brevity, we refer to them as the triaxiality and tidal accelerations. We then define and where and R (χ) and I (χ) are defined in Eq. (N.17 mod.).
In terms of these angular accelerations, the spin-orbit ODE becomes and this is our starting point for the rest of the paper.

The NFME model with practical values
We now look at the behaviour and magnitude of the various functions that go to make up the NFME model, using the values for the Sun and Mercury given in Table 1.

Hansen coefficients
In Fig. 1, we plot the Hansen coefficients relevant to the problem, G 20q = X −3,2 q+2 , for e = 0.2056, 0.3 and 0.4 and for q = −12, . . . , 12. These were computed from Eq. (N.9) by truncating the first sum at p = 120 and using 20 significant figures for computation-these values allow the computation of the coefficients to a much higher precision than necessary.

The triaxiality acceleration
We first make a simple estimate of the order of magnitude ofθ TRI (θ, t). Treating θ and t as independent variables, it is clear from Eq In Fig. 2, we give a plot ofθ TRI (θ 0 +θ 0 t, t) for e = 0.2056, θ 0 = 2.0,θ 0 = 29.0, the latter two values being (almost arbitrarily) chosen to approximate a solution starting from an angular velocity slightly greater than n. Note thatθ TRI ∈ [−0.19, 0.19], which is consistent with the estimate, from Eq. (4), of [−0.21, 0.21] when e = 0.2056. As we shall see, for example from Fig. 3, the triaxiality term dominates the tidal term by several orders of magnitude for most of the time.

The tidal acceleration
We now discuss the tidal accelerationθ TIDE (θ). Throughout this section, we take e = 0.2056. Figures 3, 4 and 5 show the tidal acceleration, plotted vertically, versus the relative rate of rotation,θ/n. This function is illustrated in Fig. 3, which showsθ TIDE (θ), plotted over the entire range of interest,θ/n ∈ [−1, 5]. The dominant features here are the 'kinks', which occur aṫ θ/n = 1 + q/2 for q = −1, . . . , 7, this being the range of the sum, Q TIDE , defining the tidal acceleration.
In order to compute probabilities of capture by the different solutions, with small confidence intervals, via Monte Carlo simulation, many solutions to Eq. (3), starting from random Fig. 6 The derivative of the tidal angular acceleration with respect to the relative rate of rotation, that is, n dθ TIDE /dθ , showing the cusp atθ/n = 3/2. Sincë θ TIDE is C 1 , its first derivative exists for allθ initial conditions in the set Q = [0, π] × [θ min ,θ max ], must be computed. Many are needed because, for a given confidence level, the width of the confidence interval is proportional to I −1/2 , where I is the number of solutions computed: more simulations narrow the confidence interval, but rather slowly. We take Q = [0, π] × [0, 5n] in what follows. The triaxiality acceleration depends on 2θ , and so we need only consider θ 0 ∈ [0, π]; and the right-most kink is atθ = 9n/2, so we choose the maximum value ofθ 0 to be somewhat greater than this.
We now describe a fast algorithm for solving Eq. (3), which is based on the high-order Euler method (HEM), described in detail in Bartuccelli et al. (2015). This method has to be adapted for the NFME model because of the discontinuities in the second derivative of θ TIDE (θ)-that is, at the centres of the kinks-that occur atθ/n = 1/2, 1, 3/2, . . . , 9/2. The adaptation we use requires that the subset Q of the (θ,θ) phase plane be split into strips, all with θ ∈ [0, π), but with a variety of ranges ofθ . This splitting is needed in order to meet the error criterion. A strip can be of type H (HEM), when it is far enough removed from the kinks that the HEM can be used, and type N (Numerical), surrounding a kink, where a suitable numerical method has to be used. In a type N strip, there is a possibility that the ODE is stiff, and the numerical method chosen should take this into account.
The details of this splitting are given in Fig. 7, in whichθ/n is plotted horizontally, with type H regions being shown as continuous lines, type N as dashed lines. The strips are grouped into ten bands, one band covering the region between two kinks. For instance, Band 0 consists of four type H strips, centred onθ/n = 0.1, 0.275, 0.39 and 0.45 and with widths 0.2, 0.15, 0.08 and 0.04, respectively; and one type N strip, whereθ/n ∈ [0.47, 0.5], θ = 0.5n being the position of the first kink. The reasoning behind the choice of these values is given in Sects. 4.1 and 4.2. The expansion point θ e for the series solution is always the centre of the strip, and is shown by 'x' in Fig. 7.

Region H: HEM applies
The HEM, which is a fixed timestep implementation of the Frobenius method, can in principle be used to solve Eq. (3)  andθ TIDE are analytic. In other words, the HEM is a numerically implemented analytical continuation of θ(t) on a set of overlapping discs of fixed radius in the complex t-plane (hence, fixed timestep), where derivatives of θ of any order are computed from the differential equation (3). Now, from its definition,θ TRI (θ, t) is an entire function for all θ, t, so this is no bar to using the HEM. However,θ TIDE (θ) is C 1 inθ, with its second derivative being undefined atθ/n = 1/2, 1, 3/2, . . . , 9/2. Hence, we approximateθ TIDE by its Taylor series of degree D TID about a pointθ =θ e , bearing in mind that this will only be good for values oḟ θ far from the kinks. In the work reported here, our error criterion is satisfied with D TID = 25 for all strips. With these provisos, the HEM works well. We briefly review the method here; full details can be found in Bartuccelli et al. (2015).
Let the state vector x(t) = θ(t),θ(t) and define t i = ih, i ∈ N 0 , where h > 0 is a timestep. The ODE allows us to compute, recursively, derivatives of θ(t) of all orders, far from the kinks. Let the j-th time derivative of θ(t) be written as θ ( j) (t). We then write the degree-D s series solution to Eq. (3) as where D s > 1 (D s = 1 gives the Euler method), and f j ( in t ofθ TRI as T 0 = 2π/n (which is one 'Mercury year'), and put x k = x(kT 0 ). Then, the Poincaré map, P, is defined by x k+1 = P(x k ). Iteration of P many times starting from a given initial condition x 0 = (θ (0),θ(0)) enables us to generate a sequence of 'snapshots' of the state variables as the system evolves, from which we can deduce the attractor in which the system is eventually captured. For details of capture detection, see below. Using a large number I of initial conditions, one can then estimate the probability of capture in any of the possible steady-state solutions.
In order to satisfy the error criterion, timestep h needs to be sufficiently small. As with all the parameters mentioned so far, there is a trade-off between speed and accuracy; a choice of h = T 0 /M with M = 24 is a good compromise in practice. Finally, then, the map P can be built up from the functions p i by P(x) = p M • p M−1 • . . . • p 1 (x), and this is the HEM.
The expressions for p i are derived by computer algebra and are initially very large, but most of the terms are negligible and hence can be pruned away. By 'term' we mean here a polynomial inθ -these are then multiplied by powers of sin 2θ, cos 2θ in order to make up p i . In practice, every term whose magnitude is less than 10 −16 at the largest and smallest values ofθ within a strip is deleted. The value of 10 −16 is chosen because the usual double precision arithmetic is carried out to around 16 significant figures. The resulting expression is then converted to Horner form (Press et al. 1992) for efficient evaluation. Typically, after pruning has been carried out, the expressions for p i contain 20-35 terms and are typically of the form: θ -component = θ + r 1 + sr 2 + cr 3 + scr 4 + s 2 r 5 ,θ -component = r 6 + sr 7 + cr 8 + scr 9 + s 2 r 10 + s 2 cr 11 + s 3 r 12 , where the r k are polynomials inθ(t i−1 ) of degree up to about 12, and s = sin 2θ , c = cos 2θ . An explicit example follows; this is the θ -component of p 1 (x, y), where (x, y) = (θ (0),θ(0)). It is valid for the middle strip of Band 3 (see Fig. 7 + 0.542e − 20y)y)y)y)y)y)y)y)y)y)y evaluated to 3 significant figures (in practice, it would need to be evaluated to 18 s.f.). Note that there are 48 strips and 24 integrators in each strip, one for each timestep, and two sets of integrators, one for θ and one forθ . Hence, the total number of integrators is 48 2 = 2304.
We now describe the error criterion used. The final version of the code for computing the Poincaré map via the HEM is compared to a high-accuracy, standard numerical computation of the same thing. The full expression forθ TIDE is used in the accurate numerical computation, not its series approximation. The numerical algorithm used is the standard Runge-Kutta method as implemented in the computer algebra software Maple, computing to 25 significant figures and with absolute and relative error tolerances of 10 −15 . The numerical and HEM computations of P(x) are then compared, for each of the type H strips, using 250 uniformly distributed random values of x in each strip. The comparison gives the maximum value of the modulus of the difference between each component, computed both ways, over the 250 random points. The maximum difference observed over all strips is assumed to be representative of the overall maximum difference. Its values are about 3×10 −14 , 1.4×10 −13 in the θ -andθ-components, respectively.
It is possible that the ODE may be stiff in type N strips, so we choose two numerical methods and compare the results. The methods used are: (1) an explicit, order 7 Runge-Kutta (RK) method due to Dormand and Price, as described in Hairer et al. (1993), and (2) the Adams method/backward differentiation formulae (BDF) (Radhakrishnan and Hindmarsh 1993), with the ability to switch automatically between them. The Adams method is a variable-order (up to 12 in the implementation used here), explicit predictor-corrector method, which, along with RK, is suitable for non-stiff problems, whereas BDF is suitable for stiff problems.
In practice, with the parameters in Table 1, BDF is not used-the problem turns out to be insufficiently stiff-so the comparison between (1) and (2) above comes down to comparing the RK and Adams methods. The implementation of Adams used (Radhakrishnan and Hindmarsh 1993) is approximately 1.6 times slower than RK for this problem, but the probabilities obtained from integration starting from the same 3200 random points in Q are in good agreement-see Table 2-so, for type N strips, we choose RK, for which we fix both the absolute and relative error tolerances to be 2 × 10 −14 . This value is chosen to be comparable with the error entailed by polynomial interpolation.
The probabilities obtained with RK and Adams are not identical, neither should we expect them to be. The long-term fate of a given trajectory depends very sensitively on the details of its computation, implemented for very long integration times. What is important in the end These were computed using the Runge-Kutta and the Adams method in N-type strips, with HEM being used elsewhere. The same 3200 uniformly distributed random initial conditions were used in both cases. Confidence intervals in brackets are unreliable since there are too few data points in the case of these attractors is the probabilities themselves, and Table 2 shows these to be robust against the algorithm used to solve the ODE. Computation time can be saved by efficient calculation ofθ TRI andθ TIDE . The expression forθ TRI can be converted into a polynomial of degree 1 in cos 2θ and sin 2θ , and degree 8 in cos nt/2 and sin nt/2. In Horner form, this polynomial can be evaluated efficiently, using 2 sin/cos evaluations, 16 addition and 36 multiplication operations.
Evaluation ofθ TIDE (θ) in the obvious way is computationally expensive, since it requires the calculation of one fractional power per term in Eq. (N.11b)-nine in all. A more efficient way to evaluate it is: Case 1 ifθ is far from a kink, then use a pre-computed Chebyshev polynomial fit (Press et al. 1992) toθ TIDE -the function is very smooth here; Case 2 ifθ is close to a kink, compute the contribution toθ TIDE from that kink exactly, according to the appropriate single term in the sum (N.11b), with the effect of the remaining kinks being replaced by a Chebyshev polynomial fit.
Hence, at most one fractional power is computed per evaluation ofθ TIDE . In practice, we use a polynomial only (Case 1 above), of degree 25, unless 2θ/n is within ±0.08 of an integer, when Case 2 applies. In Case 2, we use a polynomial of degree 7 to fit the remaining terms. The resulting absolute error is no more than 4 × 10 −14 . For comparison, the ratio of the CPU-time taken to evaluateθ TIDE directly, and via polynomial fitting, is about 5.1. A comparison of the timings in different circumstances is given in Table 3, for which 1000 random initial conditions were used. The CPU-time taken to iterate P 10 5 times starting from each of these was measured, with any data in which the trajectory moved from a type H strip to one of type N, or vice versa, being rejected. As mentioned earlier in this section, the Adams method, for this problem at least, is about 1.6 times slower than RK.

Capture test
We now describe the test used to detect capture of a solution. Figure 8 showsθ k =θ(kT 0 ) plotted against k for 0 ≤ k ≤ 7.8×10 6 . In this case, the initial condition was x 0 = (0, 49) and capture took place after about 7.7 × 10 6 iterations of the Poincaré map, which corresponds to 7.7 × 10 6 T 0 = 1.85 × 10 6 year.
There are several ways that capture could be detected. We choose to divide theθ k dataset into blocks of length L and compute the least squares gradient of each block. As can be seen from Fig. 8, this gradient will be negative pre-capture, and close to zero post-capture.  (0, 1.878n), is plotted against k, for 0 ≤ k ≤ 7.8 × 10 6 . In this case, capture takes place in a solution with θ /n ≈ 3/2. Two small regions are shown on a magnified scale. On the left, we see the pre-capture dynamics and on the right, post-capture. The HEM was used for all computations outside the strip 1.47 ≤θ/n ≤ 1.53, marked 'Type N' In detail the test is as follows. Letȳ j be the mean value ofθ k /n over the j-th block of length L, and let m j be the least squares gradient ofθ k plotted against k in block j. Capture is deemed to have taken place when 2ȳ j − 2ȳ j < ε i and m j < ε m for K successive blocks.
Here, [x] is the nearest integer to x, and the factor of two occurs in the first expression because capture can take place at integer or half-integer values ofθ/n. In practice, we choose L = 10,000, K = 8, ε i = 10 −3 and ε m = 3 × 10 −7 . The speed of a probability computation depends on the choice of these parameters. Figure 9 shows a typical plot of m j versus block number j, from j = 0 until capture, which takes place at j = 770.

A probability of capture computation
The main purpose of this paper has been to show how the computation of capture probabilities can be done relatively fast. We now wish to show how this works in practice, by producing probability of capture data for a particular case, a case chosen to allow the algorithm to exhibit its effectiveness. To this end, we now give the results of a calculation of probability of capture for the parameters in Table 1, which are of intrinsic interest as well as being a useful illustration.
As in Bartuccelli et al. (2015), we use the CPU-sec as a unit of time, which is defined in terms of the following sum: whose evaluation requires 5m multiplications and 6m − 1 additions. We define 1 CPU-sec to be the CPU-time taken to evaluate S(N c ), where N c = 5.96 × 10 7 , which is the CPU-time taken to do this computation on the computer used to do some of the calculations in this paper. The time taken can vary according to circumstances, e.g. the loading, the type of processor and so on. Hence, care has to be taken in the codes to scale the CPU-sec appropriately for the particular hardware used: the computation of S(N c ) is timed after every successful capture, and the CPU-sec scaling factor is updated on each occasion.
One other point to note is that a computation of probability is trivial to parallelise. The initial conditions in Q are generated by a pseudo-random number generator, and a different sequence of pseudo-random numbers can be produced just by changing the seed. Hence, by running N copies of the programme on N separate processors with N different seeds, N times the number of initial conditions can be investigated at the same time.
The CPU-time taken to integrate until capture depends strongly on the initial condition, and in particular, the proportion of the integration that is carried out in type N strips (which is a slow process) as opposed to type H strips (where it is fast). By considering 57,600 random initial conditions, we find the following: Mean time to capture: 1156 CPU-sec, with standard deviation 1190 CPU-sec Maximum time to capture: 5161 CPU-sec; minimum time to capture: 1.188 CPU-sec.
From the above, we see that overall, about 1.4 × 10 12 iterations of the Poincaré map were needed to estimate the probability of capture for 57,600 initial conditions, and that 12% of these were in type N strips, with the remaining 88% being in type H strips. Since capture can only take place in the vicinity of a kink, iteration until capture always requires some iteration in N strips. The values given above have a sensitive dependence on the capture parameters.
We can now estimate the factor by which our approach speeds up a typical probability of capture computation. Let T mix (k) be the CPU-time taken to iterate the Poincaré map k times using HEM in type H strips and the chosen numerical method in type N strips, and let T num (k) be the time taken when using the numerical method everywhere. Then, using the data above and from Our approach therefore speeds up this computation by a factor of about 7.5. A histogram of capture times is given in Fig. 10. Let τ (x 0 ) be the number of CPU-sec required to iterate the Poincaré map, starting from x 0 , until the capture criterion of Sect. 4.3 is met. In order to produce Fig. 10, 57,600 random values of x 0 ∈ Q were generated and τ (x 0 ) was computed for each. This figure is a histogram of the proportions of the values of τ (x 0 ) that lie in the ranges 0-200, 200-400, …, 4800-5000 CPU-sec.
Finally, we give two estimates of the probabilities themselves in Table 4. The 95% CI (Walpole et al. 1998) is defined such that the probability of capture in a given attractor A lies in the interval [p − p,p + p], with 95% confidence. Here, p = 1.96 p(1 −p)/I , in which I = 5.76 × 10 4 is the number of initial conditions for the main results, andp is the proportion of these initial conditions that end up in A. For these values, in the worst case, which isp = 1/2, p ≈ 0.41%. We also give, for comparison, an estimate of the same probabilities but obtained by using RK only, with I = 3600. The smaller number of points is due to the much slower computation speed when HEM is not used, but the probabilities obtained are in agreement. The mean computer time to capture for the slow computation was 9055 CPU-sec-c.f. 1156 CPU-sec when HEM is used, the latter being quicker by a factor of 7.8, in close agreement with the claimed speed-up factor of 7.5.

Conclusions
The main purpose of this paper has been to show how the computation of probability of capture in various solutions to the spin-orbit ODE, in the case of constant eccentricity e and with a particular model of dissipation (Noyelles et al. 2014), can be accelerated. Estimation of these probabilities requires us to solve the ODE (3), and bearing in mind that this is periodic in time with period T 0 = 2π/n, our objective is to compute the Poincaré map, which advances the state variables x = (θ,θ) by a time T 0 , as fast as possible. This enables us to generate a sequence of values x k , from which salient information about the dynamics can be deduced. In particular, we have in mind the application of this work to estimation of probability of capture in various orbits, and the confidence intervals of these probabilities shrink as more orbits are investigated: hence the need for a fast algorithm.
The dissipation term,θ TIDE (θ), which is a C 1 function of the angular velocityθ , and additionally varies rapidly in the ranges ofθ around the so-called kinks, complicates the computation, and requires the use of a standard numerical ODE solver for some of the time. We point out ways in which these purely numerical computations can be carried out efficiently, by streamlining the calculation of the triaxiality and tidal accelerations. However, in about 88% of the phase plane the high-order Euler method (Bartuccelli et al. 2015) can be used, and this speeds up the computation significantly. Additionally, the triaxiality term, being a smooth function of θ and t, is taken into account in the HEM at almost no additional computational cost.
The present case should be contrasted with the constant time lag model, in which dissipation is just proportional toθ − ω with ω a constant. In the light of its simplicity, this has been used in many publications, for instance, Bartuccelli et al. (2015) and references therein. In that case, only a single set of integrators p i , defined in Eq. (5), was needed to build up the Poincaré map. In other words, there was only one strip, which was the whole of Q, the subset of the phase plane considered, and the HEM could be used everywhere.
Compare that with the current case using the parameters in Table 1, which relate to Mercury and the Sun. Here, Q has to be divided into 48 strips and an integrator set up for each. Setting up the codes for the HEM in each strip, which is done by computer algebra, itself takes significant time, but the payoff is an increase in speed by a factor of approximately 65 in the strips where HEM can be used, compared to using a standard numerical method.
A probability of capture computation requires a capture detection algorithm, and one has been described. It is based on the fact that after capture, the angular velocity has no underlying decay: a captured solution has a close to constant mean value ofθ k =θ(kT 0 ) for all k greater than the value at which capture takes place. Typically, the mean is taken over 10 4 successive k values.
A test run of our algorithm, in which 57,600 initial conditions were iterated until capture took place, reveals that the overall speed of computation is faster by a factor of about 7.5 compared to using a standard numerical algorithm alone.
An important question left unanswered in this paper is 'What is the nature of the solutions in which capture takes place?' The answer turns out to be more complicated than expected: periodic solutions with meanθ/n ≈ −1, −1/2, 1/2, . . . , 4 have been computed (highaccuracy numerics are needed). There is also numerical evidence for the existence of attracting solutions whose period is not a small integer multiple of T 0 . The dynamics of solutions to Eq. (3), both pre-and post-capture, will be subject of a future publication (Bartuccelli et al. 2017).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.