Evolution equations for two-planet solar system

For the moment, let us consider a simplified solar system that consists of the Sun and two planets. See Figure 10.1. Let the Sun be of mass $M$ and position vector ${\bf R}_s$. Likewise, let the two planets have masses $m$ and $m'$ and position vectors ${\bf R}$ and ${\bf R}'$, respectively. Here, we are assuming that $m,\,m'\ll M$. Finally, let ${\bf r} = {\bf R}-{\bf R}_s$ and ${\bf r}'={\bf R}'-{\bf R}_s$ be the position vector of each planet relative to the Sun. Without loss of generality, we can assume that $r'>r$.

Figure 10.1: A simplified model of the solar system.

In an inertial reference frame, the equations of motion of the various elements of our simplified solar system are

$$M\,\ddot{\bf R}_s = G\,M\,m\,\frac{\bf r}{r^{\,3}} + G\,M\,m'\,\frac{{\bf r}'}{r'^{\,3}},$$ (10.1)

$$m\,\ddot{\bf R} = G\,m\,m'\,\frac{{\bf r}'-{\bf r}}{\vert{\bf r}'-{\bf r}\vert^{\,3}} - G\,m\,M\,\frac{{\bf r}}{r^{\,3}},$$ (10.2)

$$m'\,\ddot{\bf R}' = G\,m'\,m\,\frac{{\bf r}-{\bf r}'}{\vert{\bf r}-{\bf r}'\vert^{\,3}} - G\,m'\,M\,\frac{{\bf r}'}{r'^{\,3}}.$$ (10.3)

It thus follows that

$$\ddot{\bf r} + \mu\,\frac{\bf r}{r^{\,3}} =G\,m'\left(\frac{{\bf r}'-{\bf r}}{\vert{\bf r}'-{\bf r}\vert^{\,3}} - \frac{{\bf r}'}{r'^{\,3}}\right),$$ (10.4)

$$\ddot{\bf r}' + \mu'\,\frac{{\bf r}'}{r'^{\,3}} = G\,m\left(\frac{{\bf r}-{\bf r}'}{\vert{\bf r}-{\bf r}'\vert^{\,3}} - \frac{{\bf r}}{r^{\,3}}\right),$$ (10.5)

where $\mu= G\,(M+m)$ and $\mu'=G\,(M+m')$. The right-hand sides of these equations specify the interplanetary interaction forces that were neglected in our previous analysis. These right-hand sides can be conveniently expressed as the gradients of potentials:

$$\ddot{\bf r} + \mu\,\frac{\bf r}{r^{\,3}} =\nabla {\cal R},$$ (10.6)

$$\ddot{\bf r}' + \mu'\,\frac{{\bf r}'}{r'^{\,3}} = \nabla' {\cal R}',$$ (10.7)

where

$${\cal R}({\bf r},{\bf r}') = \tilde{\mu}'\left(\frac{1}{\vert{\bf r}-{\bf r}'\vert} - \frac{{\bf r}\cdot{\bf r}'}{r'^{\,3}}\right),$$ (10.8)

$${\cal R}' ({\bf r},{\bf r}') =\tilde{\mu}\left(\frac{1}{\vert{\bf r}-{\bf r}'\vert} - \frac{{\bf r}\cdot{\bf r}'}{r^{\,3}}\right),$$ (10.9)

with $\tilde{\mu} = G\,m$, and $\tilde{\mu}'=G\,m'$. Here, ${\cal R}({\bf r},{\bf r}')$ and ${\cal R}'({\bf r}, {\bf r}')$ are termed disturbing functions. Moreover, $\nabla$ and $\nabla'$ are the gradient operators involving the unprimed and primed coordinates, respectively.

In the absence of the second planet, the orbit of the first planet is fully described by its six standard orbital elements (which are constants of its motion): the major radius, $a$; the mean longitude at epoch, $\skew{5}\bar{\lambda}_0$; the eccentricity, $e$; the inclination (to the ecliptic plane), $I$; the longitude of the perihelion, $\varpi$; and the longitude of the ascending node, ${\mit \Omega }$. (See Section 4.12.) As described in Appendix G, the perturbing influence of the second planet causes these elements to slowly evolve in time. Such time-varying orbital elements are generally known as osculating elements.^10.1 Actually, when describing the aforementioned evolution, it is more convenient to work in terms of an alternative set of osculating elements, namely, $a(t)=a^{(0)}[1+\epsilon'\,a^{(1)}(t)]$, $\skew{5}\bar{\lambda}(t) = \skew{5}\bar{\lambda}_0+n^{(0)}\,t + \skew{5}\bar{\lambda}^{(1)}(t)$, $h=e\,\sin\varpi$, $k=e\,\cos\varpi$, $p=\sin I\,\sin {\mit\Omega}$, and $q=\sin I\,\cos{\mit\Omega}$. Here, $\epsilon'=\tilde{\mu}'/\mu=m'/(M+m)\ll 1$ and $n(t)=n^{(0)}[1-(3/2)\,\epsilon'\,a^{(1)}(t)]$, where $n^{(0)}=(\mu/[a^{(0)}]^3)^{1/2}$ is the unperturbed mean orbital angular velocity. In the following, for ease of notation, $a^{(0)}$ and $n^{(0)}$ are written simply as $a$ and $n$, respectively. Furthermore, $\skew{5}\bar\lambda$ will be used as shorthand for $\skew{5}\bar{\lambda}_0+n^{(0)}\,t$. The evolution equations for the first planet's osculating orbital elements take the form (see Section H.2)

$$\frac{d a^{(1)}}{dt} =\epsilon'\,n\left[2\,\alpha\,\frac{\partial ({\cal S}_0+{\cal S}_1)}{\partial\skew{5}\bar{\lambda}}\right],$$ (10.10)

$$\frac{d\skew{5}\bar{\lambda}^{(1)}}{dt} =\epsilon'\,n\left[- \frac{3}{2}\, a^{(1)}- 2\,\alpha^{\,2}\,\fra... ...al S}_1}{\partial h} + k\,\frac{\partial {\cal S}_1}{\partial k}\right)\right],$$ (10.11)

$$\frac{dh}{dt} = \epsilon'\,n\left[-\alpha\, h\,\frac{\partial {\cal S}_0}{\part... ...\lambda}} + \alpha\,\frac{\partial ({\cal S}_1+{\cal S}_2)}{\partial k}\right],$$ (10.12)

$$\frac{dk}{dt} = \epsilon'\,n\left[-\alpha\,k\, \frac{\partial {\cal S}_0}{\part... ...{\lambda}} -\alpha\,\frac{\partial ({\cal S}_1+{\cal S}_2)}{\partial h}\right],$$ (10.13)

$$\frac{dp}{dt} =\epsilon'\,n\left[ - \frac{\alpha}{2}\, p\,\frac{\partial {\cal ... ...\skew{5}\bar{\lambda}} + \alpha\,\frac{\partial {\cal S}_2}{\partial q}\right],$$ (10.14)

$$\frac{dq}{dt} = \epsilon'\,n\left[- \frac{\alpha}{2}\, q\,\frac{\partial {\cal ... ...l\skew{5}\bar{\lambda}}- \alpha\,\frac{\partial {\cal S}_2}{\partial p}\right],$$ (10.15)

where (see Section H.3)

$${\cal S}_0 =\frac{1}{2} \sum_{j=-\infty,\infty} b_{1/2}^{(j)}\,\cos [j\,(\sk... ...bar{\lambda}')] - \alpha\,\cos(\skew{5}\bar{\lambda} - \skew{5}\bar{\lambda}'),$$ (10.16)

$${\cal S}_1 = \frac{1}{2}\sum_{j=-\infty,\infty}\left\{k\,(-2\,j-\alpha\,D)\,... ...}_{1/2} \,\cos[ (1-j)\,\skew{5}\bar{\lambda}+ j\,\skew{5}\bar{\lambda}']\right.$$

$$\phantom{=}+ h\,(-2\,j-\alpha\,D)\,b^{(j)}_{1/2} \,\sin[(1-j)\,\skew{5}\bar{\lambda}+j\,\skew{5}\bar{\lambda}' ]$$

$$\phantom{=}+k'\,(-1+2\,j+\alpha\,D)\,b^{(j-1)}_{1/2} \,\cos[ (1-j)\,\skew{5}\bar{\lambda}+j\,\skew{5}\bar{\lambda}' ]$$

$$\phantom{=}\left. +h'\,(-1+2\,j+\alpha\,D)\,b^{(j-1)}_{1/2} \,\sin[ (1-j)\,\skew{5}\bar{\lambda}+j\,\skew{5}\bar{\lambda}' ]\right\}$$

$$\phantom{=}+ \frac{\alpha}{2}\left\{-k\,\cos(2\,\skew{5}\bar{\lam... ...')+ 3\,k\,\cos \skew{5}\bar{\lambda}' + 3\,h\,\sin\skew{5}\bar{\lambda}'\right.$$

$$\phantom{=}\left.-4\,k'\,\cos(\skew{5}\bar{\lambda}-2\,\skew{5}\b... ...lambda}')+4\,h'\,\sin(\skew{5}\bar{\lambda}-2\,\skew{5}\bar{\lambda}')\right\},$$ (10.17)

$${\cal S}_2 = \frac{1}{8}\,(h^{\,2}+k^{\,2}+h'^{\,2}+k'^{\,2})\,(2\,\alpha\,D... ...{1/2} - \frac{1}{8}\,(p^{\,2}+q^{\,2}+p'^{\,2}+q'^{\,2})\,\alpha\,b_{3/2}^{(1)}$$

$$\phantom{=}+\frac{1}{4}\,(k\,k'+h\,h')\,(2-2\,\alpha\,D - \alpha^{\,2}\,D^{\,2})\,b_{1/2}^{(1)}+\frac{1}{4}\,(p\,p'+q\,q')\,\alpha\,b_{3/2}^{(1)}.$$ (10.18)

Here, $\alpha=a/a'$, $D\equiv d/d\alpha$, and

$$b_s^{(j)}(\alpha) = \frac{1}{\pi}\int_0^{2\pi}\frac{\cos(j\,\psi)\,d\psi}{[1-2\,\alpha\,\cos\psi+\alpha^{\,2}]^s},$$ (10.19)

where $a'$, $\skew{5}\bar{\lambda}'$, $h'$, $k'$, $p'$, $q'$ are the osculating orbital elements of the second planet. The $b_s^{(j)}$ factors are known as Laplace coefficients (Brouwer and Clemence 1961). In deriving these expressions from Equations (10.6) and (10.8), we have expanded to first order in the ratio of the planetary masses to the solar mass; we have then evaluated the secular terms in the disturbing functions (i.e., the terms that are independent of $\skew{5}\bar{\lambda}$ and $\skew{5}\bar{\lambda}'$) to second order in the orbital eccentricities and inclinations. The nonsecular terms in the disturbing functions are evaluated to first order in the eccentricities and inclinations. (See Appendix H.) This expansion procedure is reasonable because the planets all have very small masses compared with that of the Sun, and also have relatively small orbital eccentricities and inclinations.

There is an analogous set of equations, which can be derived from Equations (10.7) and (10.9), that describe the time evolution of the osculating orbital elements of the second planet due to the perturbing influence of the first. These take the form (see Section H.2)

$$\frac{d a^{(1)'}}{dt} =\epsilon\,n'\left[2\,\alpha^{-1}\,\frac{\partial ({\cal S}_0'+{\cal S}_1')}{\partial\skew{5}\bar{\lambda}'}\right],$$ (10.20)

$$\frac{d\skew{5}\bar{\lambda}^{(1)'}}{dt} =\epsilon\,n'\left[- \frac{3}{2}\, a^{(1)'}+2\,\frac{\partial ({\... ..._1'}{\partial h'} + k'\,\frac{\partial {\cal S}_1'}{\partial k'}\right)\right],$$ (10.21)

$$\frac{dh'}{dt} = \epsilon\,n'\left[-\alpha^{-1}\, h'\,\frac{\partial {\cal S}_0'... ...} + \alpha^{-1}\,\frac{\partial ({\cal S}_1'+{\cal S}_2')}{\partial k'}\right],$$ (10.22)

$$\frac{dk'}{dt} = \epsilon\,n'\left[-\alpha^{-1}\,k'\, \frac{\partial {\cal S}_0'... ...'} -\alpha^{-1}\,\frac{\partial ({\cal S}_1'+{\cal S}_2')}{\partial h'}\right],$$ (10.23)

$$\frac{dp'}{dt} =\epsilon\,n'\left[ - \frac{\alpha^{-1}}{2}\, p'\,\frac{\partial ... ...\bar{\lambda}'} + \alpha^{-1}\,\frac{\partial {\cal S}_2'}{\partial q'}\right],$$ (10.24)

$$\frac{dq'}{dt} = \epsilon\,n'\left[- \frac{\alpha^{-1}}{2}\, q'\,\frac{\partial ... ...}\bar{\lambda}'}- \alpha^{-1}\,\frac{\partial {\cal S}_2'}{\partial p'}\right],$$ (10.25)

where (see Section H.3)

$${\cal S}_0' =\frac{\alpha}{2} \sum_{j=-\infty,\infty} b_{1/2}^{(j)}\,\cos [j\... ...pha^{-1}\,\cos(\skew{5}\bar{\lambda}' - \skew{5}\bar{\lambda}),\displaybreak[0]$$ (10.26)

$${\cal S}_1' = \frac{\alpha}{2}\sum_{j=-\infty,\infty}\left\{k\,(-2\,j-\alpha\... ..._{1/2} \,\cos[ j\,\skew{5}\bar{\lambda}' +(1-j)\,\skew{5}\bar{\lambda} ]\right.$$

$$\phantom{=}+ h\,(-2\,j-\alpha\,D)\,b^{(j)}_{1/2} \,\sin[j\,\skew{5}\bar{\lambda}' +(1-j)\,\skew{5}\bar{\lambda}]$$

$$\phantom{=}+k'\,(-1+2\,j+\alpha\,D)\,b^{(j-1)}_{1/2} \,\cos[ j\,\skew{5}\bar{\lambda}'+(1-j)\,\skew{5}\bar{\lambda}]$$

$$\phantom{=}\left. +h'\,(-1+2\,j+\alpha\,D)\,b^{(j-1)}_{1/2} \,\sin[ j\,\skew{5}\bar{\lambda}' +(1-j)\,\skew{5}\bar{\lambda}]\right\}$$

$$\phantom{=}+ \frac{\alpha^{-1}}{2}\left\{-k'\,\cos(2\,\skew{5}\ba... ...})+ 3\,k'\,\cos \skew{5}\bar{\lambda} + 3\,h'\,\sin\skew{5}\bar{\lambda}\right.$$

$$\phantom{=}\left.-4\,k\,\cos(\skew{5}\bar{\lambda}'-2\,\skew{5}\b... ...,\sin(\skew{5}\bar{\lambda}'-2\,\skew{5}\bar{\lambda})\right\},\displaybreak[0]$$ (10.27)

$${\cal S}_2' =\frac{1}{8}\,(h^{\,2}+k^{\,2}+h'^{\,2}+k'^{\,2})\,\alpha\,(2\,\a... ...- \frac{1}{8}\,(p^{\,2}+q^{\,2}+p'^{\,2}+q'^{\,2})\,\alpha^{\,2}\,b_{3/2}^{(1)}$$

$$\phantom{=}+\frac{1}{4}\,(k\,k'+h\,h')\,\alpha\,(2-2\,\alpha\,D - \alpha^{\,2}\,D^{\,2})\,b_{1/2}^{(1)}$$

$$\phantom{=}+\frac{1}{4}\,(p\,p'+q\,q')\,\alpha^{\,2}\,b_{3/2}^{(1)}.$$ (10.28)

Here, $\epsilon=\tilde{\mu}/\mu' = m/(M+m')\ll 1$, and $n'=(\mu'/a'^{\,3})^{1/2}$.