Secular Dynamics for Curved Two-Body Problems

Abstract

Consider the dynamics of two point masses on a surface of constant curvature subject to an attractive force analogue of Newton’s inverse square law, that is under a ’cotangent’ potential. When the distance between the bodies is sufficiently small, the reduced equations of motion may be seen as a perturbation of an integrable system. We take suitable action-angle coordinates to average these perturbing terms and describe dynamical effects of the curvature on the motion of the two-bodies.

Appendices

Action-angle coordinates for the curved Kepler problem

1.1 Curved Conics

The Kepler problem on a sphere of radius $\rho$ with a fixed ’sun’, $q_s$, of mass M and particle, $q\in S^2$, of mass m is given by the Hamiltonian flow of:

$$\begin{aligned} Kep_\kappa := \frac{\Vert p\Vert _\kappa ^2}{2m} - \frac{mM}{\rho }\cot \varphi \end{aligned}$$

on $T^*(S^2\backslash \{|\cot \varphi | = \infty \})$. Here $\Vert \cdot \Vert _\kappa$ is the norm induced by the metric of constant curvature $\kappa = 1/\rho ^2$ on $S^2$ and $\varphi$ is the angular distance from $q_s$ to q (see Fig. 2). For a negatively curved space, one replaces $\cot \varphi$ with $\coth \varphi$.

Because the force is central, letting $\mathbf {q}\in \mathbf {R}^3$ be the position of the particle from the center, c, of the sphere and ${\hat{k}}$ a unit vector from c to $q_s$, the angular momentum:

$$\begin{aligned} G:=m(~\mathbf {q}\times \dot{\mathbf {q}}~)\cdot {\hat{k}} \end{aligned}$$

is a first integral. It corresponds to the rotational symmetry about the $\mathbf {cq}_\mathbf{s }$ axis.

Remark 11

In spherical coordinates, $\mathbf {q} = \rho (\sin \varphi \cos \theta , \sin \varphi \sin \theta , \cos \varphi )$, we have: $\Vert p\Vert _\kappa ^2 = \frac{1}{\rho ^2}(p_\varphi ^2 + \frac{p_\theta ^2}{\sin ^2\varphi })$, and $G = m\rho ^2\sin ^2\varphi \dot{\theta }= p_\theta$. In these coordinates, it is not hard to show analytically that the curved Kepler trajectories centrally project to flat Kepler orbits (Fig. 2). Indeed, setting $R = 1/r := \frac{\cot \varphi }{\rho }$, one finds $\frac{d^2R}{d\theta ^2} + R = m^2M/G^2$, so that the $\theta$-parametrized orbits are: $r = \frac{G^2/m^2M}{1 + e\cos (\theta - g)}$, with e and g being constants of integration. Here $\nu := \theta - g$ is the true anomaly and g is the argument of pericenter.

The curved Kepler trajectories are in fact conic sections on the sphere having a focus at $q_s$ (see Fig. 5). One may express energy and momentum in terms of geometric parameters of such spherical conics.

Fig. 5

A spherical ellipse, with foci at $q_s$ and f is the set of points q for which $|qq_s| + |qf| = 2\alpha$ is constant. We call $\alpha$ its semi-major axis, the midpoint, $q_c$, between the foci its center, and the length $|bq_c| =:\beta$ its semi-minor axis. The number, $\epsilon$, for which $|q_sq_c| = \alpha \epsilon$ is called the eccentricity. When $\epsilon \ne 0$, the closest point to $q_s$ is called the pericenter, pc, and furthest, ac, the apocenter. The angle $\nu := \angle (pc,q_s,q)$ is called the true anomaly of q

Proposition 4

Consider an orbit of the curved Kepler problem along a spherical conic with a focus at $q_s$ (Fig. 5). Let $\alpha$ be the semi-major axis and $\beta$ the semi-minor axis of this conic. The orbit has energy and momentum:

$$\begin{aligned}&Kep_\kappa = -\frac{mM}{\rho }\cot \frac{2\alpha }{\rho },\\&G^2 = m^2M\rho \tan ^2\frac{\beta }{\rho }\cot \frac{\alpha }{\rho }. \end{aligned}$$

Proof

Consider the spherical triangle $\varDelta (f,q_s,q)$ with sidelengths $2\alpha \epsilon , \rho \varphi , 2\alpha - \rho \varphi$ and interior angle $\pi - \nu$ opposite to side ${\mathop {fq}\limits ^{\frown }}$ (see Fig. 5). The cosine rule of spherical trigonometry yields: $r = \rho \tan \varphi = \frac{p^2}{1 + e\cos \nu },$ where $p^2 = \rho \frac{\cos \frac{2\alpha \epsilon }{\rho } - \cos \frac{2\alpha }{\rho }}{\sin \frac{2\alpha }{\rho }}$ and $e = \frac{\sin \frac{2\alpha \epsilon }{\rho }}{\sin \frac{2\alpha }{\rho }}$, so that the orbits are indeed curved conic sections. Comparing with the expression in Remark 11 yields $G^2 = m^2Mp^2$. The expression for G in the proposition follows by using the relation: $\cos \frac{\beta }{\rho } = \frac{\cos \frac{\alpha }{\rho }}{\cos \frac{\alpha \epsilon }{\rho }}$ (consider the right spherical triangle $\varDelta (q_s,q_c,b)$), to simplify.

To obtain the expression for $Kep_\kappa$, following [1], observe that $\varphi _p := \frac{\alpha (1-\epsilon )}{\rho }$ and $\varphi _a := \frac{\alpha (1+\epsilon )}{\rho }$ are maximal and minimal values of $\varphi$ over the trajectory having energy $Kep_\kappa =: h$. Consequently $p_\varphi = 0$ at $\varphi _{a,p}$ so we have two solutions of the equation: $h\cos 2\varphi = \frac{mM}{\rho }\sin 2\varphi + h - \frac{G^2}{\rho ^2}.$ Adding and subtracting the above equation evaluated at $\varphi _{a,p}$ yields: $Kep_\kappa = h = \frac{mM}{\rho }\frac{\sin 2\varphi _a - \sin 2\varphi _p}{\cos 2\varphi _a - \cos 2\varphi _p} = - \frac{mM}{\rho }\cot (\varphi _a + \varphi _p).$ $\square$

Remark 12

The sign of G represents the orbits orientation, with positive G for counterclockwise motion around $q_s$. The same arguments apply to bounded motions when the curvature is negative, replacing trigonometric functions with their hyperbolic counterparts. The energy for bounded motions in a space of curvature $\kappa = - 1/\rho ^2$ is always less than $-mM/\rho$.

Remark 13

The orbits of the curved Kepler problem colliding with $q_s$ may be regularized similarly to the usual Kepler regularization, via an elastic bounce [2].

1.2 Delaunay and Poincaré Coordinates

We have the following analogues of the Delaunay and Poincaré symplectic coordinates for the bounded motions of the curved Kepler problem. Note that to each point of $p\in T^*(S^2\backslash \{\pm q_s\})$, one may associate a pointed spherical conic and conversely, parametrizations of pointed conics may serve as local coordinates for $T^*(S^2\backslash \{\pm q_s\})$.

Proposition 5

Given a pointed spherical conic (Fig. 5), let G be its angular momentum, g its argument of pericenter,

$$\begin{aligned} L^2 := m^2M\rho \tan \frac{\alpha }{\rho },~~L>0, \end{aligned}$$

and $\ell$, the mean anomaly, be proportional to time along the conic to pericenter, scaled so that $\ell \in \mathbf {R}/2\pi \mathbf {Z}$. Then $(L, \ell , G, g)$ are symplectic (Delaunay) coordinates for bounded non-circular motions. The variables:

$$\begin{aligned} \Lambda = L, ~~\lambda = \ell + g, ~~\xi = \sqrt{2(L-|G|)}\cos g, ~~\eta = \sqrt{2(L-|G|)}\sin g, \end{aligned}$$

are symplectic (Poincaré) coordinates in a neighborhood of the circular motions. The energy is given by:

$$\begin{aligned} Kep_\kappa = -\frac{m^3M^2}{2L^2} + \kappa \frac{L^2}{2m}. \end{aligned}$$

(5)

Proof

The construction of these coordinates is almost identical as for the planar Kepler problem (see e.g. the presentation in [10]). We only found some difference in the computation determining L owing to, in the curved case, the period as a function of energy (Kepler’s third law) not having such a simple expression.

For non-circular motions, we have symplectic coordinates $(H = Kep_\kappa , t, G, g)$, where t is the time along the orbit (to pericenter). Since orbits of fixed energy, $Kep_\kappa = H$, all have a common period, T(H), we set $\ell = \frac{2\pi }{T(H)}t$ and seek a conjugate coordinate L(H) to $\ell$, i.e. we want to integrate:

$$\begin{aligned} dL = \frac{T(H)}{2\pi }dH. \end{aligned}$$

To integrate this expression, we make some changes of variable. For H a negative energy value admitting bounded motions, let $\varphi _c(H)$ be the angular distance of the circular orbit having energy H. Then:

$$\begin{aligned} H = - \frac{mM}{\rho }\cot 2\varphi _c,~~G_cT(H) = 2\pi m\rho ^2\sin ^2\varphi _c \end{aligned}$$

where $G_c^2 := \rho m^2M\tan \varphi _c$ is the angular momentum of the circular solution. Note that

$$\begin{aligned} (*)~~~~H = -\frac{m^3M^2}{2G_c^2} + \kappa \frac{G_c^2}{2m}. \end{aligned}$$

We compute $\frac{T}{2\pi }dH = dG_c$. So we take $L = G_c$, and have $L^2 = m^2M\rho \tan \frac{\alpha }{\rho }$ from $(*)$ and Propostion 4. $\square$

Remark 14

The mean anomaly, $\ell$, is related to time by $d\ell = n~ dt$ where $n(L) := \frac{m^3M^2}{L^3} + \kappa \frac{L}{m}$. When the curvature is negative, the same arguments lead to $L^2 = m^2M\rho \tanh \frac{\alpha }{\rho }$, and the same expression, Eq. (5), for the energy.

1.3 Other Anomalies

In averaging functions over curved Keplerian orbits, i.e. determining $\frac{1}{2\pi }\int _0^{2\pi } f(q)~d\ell$, it is often useful to perform a change of variables, as the position on the orbit, q, does not have closed form expressions in terms of $\ell$. We collect here some parametrizations of Keplerian orbits. Although not all are necessary for our main results, they may serve useful in other perturbative studies of the curved Kepler problem.

The position on the conic is given explicitely in terms of the true anomaly (Remark 11). By conservation of angular momentum (recall $d\ell = n~dt$ with $n = \frac{m^3M^2}{L^3} + \kappa \frac{L}{m}$):

$$\begin{aligned} G~d\ell =n~m\rho ^2\sin ^2\varphi ~d\nu . \end{aligned}$$

One may centrally project a curved Kepler conic to the tangent plane at $q_s$ and then parametrize this planar conic by its eccentric anomaly, which we denote here by $u_o$. Letting a, e, b be the semi-major axis, eccentricity and minor axis of this planar conic, the position is given by:

$$\begin{aligned} r&= \rho \tan \varphi = a(1 - e\cos u_o)&x&= r\cos \nu = a(\cos u_o - e)&y&=r\sin \nu = b\sin u_o \end{aligned}$$

(6)

and one has:

$$\begin{aligned} \sqrt{\frac{M}{a}}~d\ell =n~\rho \sin \varphi \cos \varphi ~du_o. \end{aligned}$$

(7)

Some more parametrizations arise naturally when one orthogonally projects the curved Kepler ellipse onto the tangent plane at $q_s$. The time-parametrized motion along this plane curve sweeps out area at a constant rate, however it is now a quartic curve: the locus of a 4th order polynomial in the plane. This quartic may be seen naturally as an elliptic curve and then parametrized by Jacobi elliptic functions. Taking $R = \rho \sin \varphi , X = R\cos \nu , Y = R\sin \nu$, the quartic is the projection to the XY-plane of the intersection of the quadratic surfaces:

$$\begin{aligned} R^2 = X^2 + Y^2, ~~~~ (R + eX)^2 = p^2(1 - \kappa R^2), \end{aligned}$$

where p, e are as in the proof of Propostion 4. Consequently we find the parametrization:

$$\begin{aligned} R&= \rho \sin \frac{\alpha }{\rho }k' \text {nd}_k w - \rho \cos \frac{\alpha }{\rho } k \text {cd}_k w&X&= \rho \sin \frac{\alpha }{\rho }k' \text {cd}_k w - \rho \cos \frac{\alpha }{\rho } k \text {nd}_k w \\ Y&= \rho \tan \frac{\beta }{\rho }\cos \frac{\alpha }{\rho }\text {sd}_k w. \end{aligned}$$

where $k = \sin \frac{\alpha \epsilon }{\rho }, k' = \cos \frac{\alpha \epsilon }{\rho }$. Since the area is swept at a constant rate, one computes:

$$\begin{aligned} \rho \sin \frac{\alpha }{\rho }~d\ell = \rho \sin \varphi ~dw. \end{aligned}$$

which integrates to give a ’curved Kepler equation’, i.e. the relation between position and time through:

$$\begin{aligned} \ell = \arccos \text {cd}_k w - \frac{\cot \frac{\alpha }{\rho }}{2}\log \frac{1 + k\text {sn}_k w}{1 - k\text {sn}_k w} \end{aligned}$$

(8)

It turns out that a geometric definition of eccentric anomaly, u (see Fig. 6), is the Jacobi amplitude of w:

$$\begin{aligned} du = \text {dn}_k w~dw, \end{aligned}$$

which can be established by using spherical trigonometry to give the position in terms of u, and some straightforward, although tedious, simplification.

Fig. 6

A ’curved eccentric anomaly’, u, for a (non-circular) Keplerian conic on the sphere. One circumscribes a (spherical) circle centered at $q_c$ around the conic and to a point q on the conic assigns the angle $u :=\angle (pc, q_c, p)$, where p is the intersection of the circle with the perpendicular dropped from q to the major axis

Remark 15

Projecting the spherical orbits from the south pole leads as well to quartic curves in the tangent plane at $q_s$, however along these quartics the time-parametrization is no longer by sweeping area at a constant rate –as it is for orthogonal projection—which we found only led to complicated expressions. Another parametrization of the orbits is presented in [12] eq. (28), and used to derive a different ’curved Kepler equation’ from our Eq. (8).

Expansion of $\langle Per\rangle$

We will compute an expansion, Eq. (12), in powers of our small parameter $\varepsilon$ for $\langle Per\rangle$ (Propostion 2) upto $O(\varepsilon ^5)$. This expansion may be found by using the ’flat eccentric anomaly’, $u_o$, of Eq. (7) to average the terms. With this approach, it is necessary to make use of the following formulas allowing one to translate between the major axis and eccentricity (a, e) of the centrally projected planar conic and our Delaunay coordinates:

$$\begin{aligned}&a = \frac{L^2}{m^2}\left( \frac{1}{1 - \frac{\kappa }{m^4}L^2(L^2-G^2)}\right) = \rho \varepsilon \frac{{\hat{L}}^2}{m^2}(1 + O(\varepsilon ^2)),\\&e^2 = (1 - \frac{{\hat{G}}^2}{{\hat{L}}^2})(1 + \frac{\varepsilon ^2}{m^4}{\hat{L}}^2{\hat{G}}^2) \end{aligned}$$

(recall we set $m = m_1m_2$). Note that by Eq. (6), $r = \rho \tan \varphi = \rho O(\varepsilon )$, so indeed $\varphi = \frac{r}{\rho } + O(\varepsilon ^3)$ is $O(\varepsilon )$.

By Taylor expansion of Eq. (2) in $\varphi$, and then exchanging $\varphi$ to $\frac{r}{\rho }$, we have:

$$\begin{aligned} Per= & {} \varepsilon ^2\left( \frac{{\hat{C}}^2}{2} - {\hat{G}}^2 + \frac{2}{3}(m_2-m_1){\hat{G}}\sqrt{{\hat{C}}^2 - {\hat{G}}^2} \frac{r\cos \theta }{\rho } \right. \nonumber \\&\left. + (({\hat{C}}^2 - {\hat{G}}^2)\cos ^2\theta - {\hat{G}}^2)\sigma \frac{r^2}{\rho ^2}\right) + O(\varepsilon ^5), \end{aligned}$$

(9)

where we set $\sigma = \frac{1 -(m_1^3 + m_2^3)}{6}$. So, to determine $\langle Per\rangle = \frac{1}{2\pi }\int _0^{2\pi } Per~d\ell$ upto $O(\varepsilon ^5)$, it remains to find:

$$\begin{aligned} \langle r\cos \theta \rangle , ~~~~~~~ \langle r^2\cos ^2\theta \rangle , ~~~~~~~ \langle r^2 \rangle . \end{aligned}$$

By Eq. (7), $d\ell = n\sqrt{a}\frac{r}{1 + \kappa r^2} du_o = n\sqrt{a} r (1 - \kappa r^2 + O(\varepsilon ^4))~du_o$, where

$$\begin{aligned} n\sqrt{a} = \frac{m^2}{\rho \varepsilon {\hat{L}}^2} + O(\varepsilon ). \end{aligned}$$

Hence:

$$\begin{aligned} \langle r\cos \theta \rangle = \frac{n\sqrt{a}}{2\pi } \int _0^{2\pi }r^2\cos \theta - \kappa r^3\cos \theta ~du_o + O(\varepsilon ^3). \end{aligned}$$

Using $\theta = \nu + g$, and the expressions for $r, r\cos \nu , r\sin \nu$ of Eq. (6), we obtain:

$$\begin{aligned} \langle r\cos \theta \rangle = a^2 e n \sqrt{a}\cos g \left( -\frac{3}{2} + \kappa a ( 2 + \frac{e^2}{2}) \right) + O(\varepsilon ^3). \end{aligned}$$

Or, in terms of the (scaled) Delaunay coordinates:

$$\begin{aligned} \langle r\cos \theta \rangle = \rho \varepsilon \cos g~ {\hat{L}}\sqrt{{\hat{L}}^2 - {\hat{G}}^2}\left( - \frac{3}{2} + \frac{\varepsilon }{2\rho m^2}(5{\hat{L}}^2 - \hat{G}^2) \right) + O(\varepsilon ^3). \end{aligned}$$

(10)

Likewise, one computes:

$$\begin{aligned} \langle r^2\cos ^2\theta \rangle&= \frac{\rho ^2\varepsilon ^2}{2m^4}\hat{L}^2\left( 2{\hat{L}}^2 + (3 + 5\cos 2g)({\hat{L}}^2 - {\hat{G}}^2)\right) + O(\varepsilon ^4), \nonumber \\ \langle r^2 \rangle&= \frac{\rho ^2\varepsilon ^2}{2m^4}\hat{L}^2(5{\hat{L}}^2 - 3{\hat{G}}^2) + O(\varepsilon ^4). \end{aligned}$$

(11)

Combining eqs. (10), (11) with the averaged Eq. (9), yields:

$$\begin{aligned} \begin{aligned} \langle Per\rangle&= \varepsilon ^2\left( \frac{{\hat{C}}^2}{2} - {\hat{G}}^2\right) + \varepsilon ^3\left( m_\varDelta {\hat{L}}{\hat{G}}\sqrt{({\hat{C}}^2 - {\hat{G}}^2)({\hat{L}}^2 - {\hat{G}}^2)}\cos g\right) \\&\quad + \varepsilon ^4\left[ {\tilde{m}} {\hat{L}}^2 \left( ({\hat{C}}^2 - {\hat{G}}^2)(2{\hat{L}}^2 + ({\hat{L}}^2 - {\hat{G}}^2)(3 + 5\cos 2g)) - {\hat{G}}^2(5{\hat{L}}^2 - 3{\hat{G}}^2)\right) \right. \\&\quad \left. - \frac{m_\varDelta }{3\rho m^2}{\hat{L}}{\hat{G}}\sqrt{({\hat{L}}^2 - {\hat{G}}^2)({\hat{C}}^2 - {\hat{G}}^2)}(5{\hat{L}}^2 - {\hat{G}}^2)\cos g \right] + O(\varepsilon ^5), \end{aligned}\nonumber \\ \end{aligned}$$

(12)

where we set $m_\varDelta = m_1-m_2, {\tilde{m}} = \frac{1 - m_1^3 - m_2^3}{12m^4}$, and $m = m_1m_2$.

Remark 16

When the masses are equal, one may avoid the need to take expansions as we have done here by using the simplified expression in Remark 3, and making use of the different anomalies presented in Sect. A.3 to obtain explicit, albeit still rather complicated, expressions.

In the secular dynamics, we may view ${\hat{L}}, {\hat{C}}$ as parameters, and then describe the dynamics of non-circular motions in the coordinates $({\hat{G}},g)$. To see the behaviour near circular motions (when ${\hat{G}} = {\hat{L}} \le {\hat{C}}$), one may convert Eq. (12) to Poincaré coordinates, in which one finds the circular motion is an elliptic fixed point. Taking into account higher order terms of $\langle Per\rangle$ leads to finer descriptions of the orbits. For example above we have for the most part worked at order $\varepsilon ^2$, when we see the precession properties described for our continued orbits. In Fig. 7, we plot some level curves of Eq. (12).

Fig. 7

Level sets of $\langle Per\rangle$, with ${\hat{L}} = {\hat{C}} = 1$. The vertical axis is the angular momentum, ${\hat{G}}\in [-1,1]$, and the horizontal axis the argument of pericenter, $g\in [0,2\pi ]$. In (a), we have equal masses, while in (b) non-equal masses. The two stable circular orbits are, in these coordinates, blown up to the lines ${\hat{G}} = \pm 1$

Remark 17

The equilibrium points, $G = 0, 2g\equiv 0~\mod \pi$ for equal masses, in Fig. 7a correspond to (regularized) periodic collision orbits of the reduced dynamics (for non-equal masses, one has high eccentricity ’near collision orbits’).

According to Remark 10, the stable (elliptic) points, $g = 0,\pi$, lift to collision orbits bouncing perpendicularly to the equator, $\mathbf {C}^{\perp }$ (depicted in Fig. 1). The closed curves surrounding these stable points represent invariant (punctured) tori of the reduced dynamics different than those described above as, rather than precessing at a constant rate, the osculating conics experience libration: with g oscillating around $\pi$ or zero.

The unstable (hyperbolic) points, $g = \frac{\pi }{2}, \frac{3\pi }{2}$, lift to collision orbits bouncing along an arc aligned with the equator. In the true dynamics, one expects a splitting of the seperatrices connecting these unstable orbits. Note that, by Propostion 3, such a splitting would be exponentially small, i.e. of $O(\exp (-1/\varepsilon ))$. If one could establish such a transversal intersection, one would obtain random motions near these unstable collision orbits admitting the following description. For any sequence $s_0,s_1,...$ with $s_k\in \{0,1\}$, there would exist a near collision orbit for which, during the time interval $[nT, (n+1)T]$ the osculating conic precesses $s_n$ times around, and where $T>0$ is some sufficiently long time period.