2 First order method

We consider the equations of motion of a planetesimal assuming the planar, elliptical restricted three-body problem where the gas drag term is added to the equations. Since drag brings planetesimals to the central plane where the planet is supposed to orbit, the planar model is thus justified. Lagrange's equations are considered for the variation of the orbital elements a, e, $\varpi$ and l where a is the semi-major axis, e is the eccentricity, $\varpi$ is the longitude of the perihelion and l is the mean anomaly. In these variables the equations for the variation of the osculating elements are given by (Brouwer & Clemence 1961)

$$\begin{eqnarray*}\dot a = {{2}\over{na}} {{\partial\Re}\over{\partial l}} +\dot a_{ {\rm f.d.}} \end{eqnarray*}$$

$$\dot e = -{{\sqrt{1-e^2}}\over{na^2e}} {{\partial\Re}\over{\p... ...na^2e}} {{\partial\Re}\over{\partial l}} +\dot e_{ {\rm f.d.}}$$

$$\dot \varpi = {{\sqrt{1-e^2}}\over{na^2e}} {{\partial\Re}\over{\partial e}} +\dot \varpi_{ {\rm f.d.}}$$

$$\dot l=n\!-\! {{2}\over{na}}\! \left({{\partial\Re}\over{\p... ...\rm f.d.}}\!=\! n \!-\! \dot \sigma \!+\! \dot l_{ {\rm f.d.}}$$ (1)

$\Re$ is the disturbing function initially assumed in the set of variables $(a,e,\varpi,f,a_{\rm P},e_{\rm P},\varpi_{\rm P},f_{\rm P})$ for the particle and the planet where f is the true anomaly and the subscript ${\rm P}$ refers to the planet.

The derivatives with respect to the variables in Eqs. (1) are obtained through the chain rule from the expressions

$$\tan \left({f\over2}\right)=\sqrt{{1+e\over1-e}}\tan \left({u\over2}\right)$$ (2)

and

$$u-e\sin u=l.$$ (3)

In the end, the partial derivatives are left as functions of u and $u_{\rm P}$ as fast variables where u is the eccentric anomaly.

We also consider, instead of l, the resonant variable ($\Phi$) given by

$$\Phi = (j+k)(l+\varpi)-j(l+\varpi)_{P}-k\varpi$$ (4)

where, j>0 and k are integers with k>0 defining an exterior resonance and k<0 defining an interior one.

Now we have to change the variable l to $\Phi$ in Eqs. (1). A new perturbing function $\Re^\ast$ is defined as an explicit function of a, e, $\Phi$ and $\varpi$, after replacing l by $\Phi$ and $\varpi$ through Eq. (4). This also changes the explicit dependence of $\Re$ with $\varpi$, thus we define $\Re^\ast=\Re^\ast(a,e,\varpi^\ast,\Phi)=\Re(a,e,\varpi,l(\Phi,\varpi))$. The derivatives of the new perturbing function with respect to the new variables are related to the derivatives with respect to the old variables by

$${{\partial\Re^\ast}\over{\partial \varpi\ast}}= {{\partial\Re... ...partial \varpi}}-{j\over(j+k)}{{\partial\Re}\over{\partial l}}$$ (5)

and

$${{\partial\Re^\ast}\over{\partial \Phi}}= {\partial l \over \... ...artial l}}= {1\over(j+k)}{{\partial\Re}\over{\partial l}}\cdot$$ (6)

Equation (1) written in these new set of variables becomes

$$\dot a = {{2}\over{na}}(j+k){{\partial\Re^\ast}\over{\partial \Phi}} +\dot a_{ {\rm f.d.}}$$ (7)

$$\begin{eqnarray*}\dot e =&-&{{\sqrt{1-e^2}}\over{na^2e}} [j-(j+k)\sqrt{1-e^2}] {... ...rtial\Re^\ast}\over{\partial \varpi^\ast}} +\dot e_{ {\rm f.d.}} \end{eqnarray*}$$

$$\dot \varpi^\ast = {{\sqrt{1-e^2}}\over{na^2e}} {{\partial\Re^\ast}\over{\partial e}} +\dot \varpi_{ {\rm f.d.}}$$

$$\dot \Phi =(j+k)(\dot l+\dot \varpi^\ast)-jn_{2}-k\dot \varpi^\ast+\dot \Phi_{ {\rm f.d.}}$$

where $\dot l$ is given by Eq. (1). Of course $\varpi$ and $\varpi^\ast$have the same geometrical meaning and after the new equations are defined as above, we can hereafter just drop the stars for simplification.

The dissipative terms $\dot a_{ {\rm f.d.}}$, $\dot e_{ {\rm f.d.}}$, $\dot \varpi_{ {\rm f.d.}}$ and $\dot \Phi_{ {\rm f.d.}}$ are obtained from Gauss equations (Brouwer & Clemence 1961)

$$\dot a_{ {\rm f.d.}}={{2}\over{n\sqrt{1-e^2}}} [Re\sin f+S(1+e\cos f)]$$ (8)

$$\dot e_{ {\rm f.d.}}={{\sqrt{1-e^2}}\over{na}} [R\sin f+S{{(e+2\cos f+e\cos ^2f)}\over{(1+e\cos f)}}]$$

$$\dot \varpi_{ {\rm f.d.}}={{\sqrt{1-e^2}}\over{nae}} [-R\cos f+S{{(2+e\cos f)}\over{(1+e\cos f)}}\sin f]$$

$$\dot \epsilon_{ {\rm f.d.}}=-{{2(1-e^2)}\over{na}}R+{{e^2}\over{1+e\sqrt{1-e^2}}} \dot \varpi_{ {\rm f.d.}}.$$

For the resonance angle $\Phi_{ {\rm f.d.}}$, the equation is given by

$$\dot \Phi_{ {\rm f.d.}}= (j+k)\dot \epsilon_{ {\rm f.d.}}-k\dot \varpi_{ {\rm f.d.}}$$ (9)

S and R are the radial and transverse components of the dissipative force. We consider two expressions for the dissipative force per unit mass:

$$\vec f=-C\vec V$$ (10)

or

$$\vec f=-CV\vec V$$ (11)

where $\vec V=\vec v-\alpha \vec v_{ {\rm kep}}$, C being the drag coefficient (it relates physical properties of the gas and planetesimal's radius and mass), $\vec v$ the planetesimal velocity and $\vec v_{ {\rm kep}}$ the Keplerian velocity at the planetesimal location; the parameter $\alpha$ is the ratio of the gas speed to the Keplerian speed at the same point, taken as $\alpha=0.995$ (Adachi et al. 1976). The parameter that defines the use of either Eq. (10) or Eq. (11) is the Reynolds number given by $R_{\rm e}={{2 \rho_{ {\rm gas}} V {\rm s}} \over {\mu}}$where $\rho_{ {\rm gas}}$ is gas density, s is the planetesimal radius and $\mu$ is the viscosity of the gas. For high Reynolds number ( $R_{\rm e}>800$) (Weidenschilling 1977), the gas drag acceleration is given by Eq. (11) and for low Reynolds number ( $R_{\rm e}<1$) the gas drag acceleration is given by Eq. (10). Hereafter we will refer to these expressions as Stokes and v² gas drag for Eqs. (10) and (11), respectively. The choice of the right gas drag law thus depends on several parameters. There is anyway a tendency for the v² law to be valid for the inner region of a planetary system since V and $\rho_{ {\rm gas}}$ are higher. Also, given standard parameters for the primordial solar nebula, in many cases trapped planetesimals will experience drag in an intermediate regime ( $1 < R_{\rm e} < 800$) (Mothé-Diniz & Gomes 2000). Anyway, for illustration purposes, in what follows, we will give examples for Stokes drag in the case of exterior resonances and both Stokes and v² drag for interior resonances.

A corotation point is defined as a solution point of the system of equations: $<\dot a>= <\dot e> = <\dot \varpi>$ $= <\dot \Phi>=0$, where the average is taken over a synodic period T. Defining a generic variable $\vec \xi=(a,e,\varpi,\Phi)=(\xi_{1},\xi_{2},\xi_{3},\xi_{4})$ and using the definition of an average value for any $\xi$ we have

$$<\dot \xi_{i}>={{1}\over{T}}\!\int_0^T \!\dot \xi_{i} {\rm d}... ...=\! {{1}\over{T}}[\xi_{i}(T)-\xi_{i}(0)]\!=\!0 , i\!=\!1,...,4$$ (12)

or

$$\xi_{i}(T)=\xi_{i}(0)$$ (13)

which means that all the orbits around an exact corotation point are periodic on the synodic period T. For a point $\vec \xi_{0}=(a_{0},e_{0},\varpi_{0},\Phi_{0})$ in the vicinity of an exact corotation point, we define:

$$F_{i}(\vec \xi_{0})=<\dot \xi_{i}(\vec \xi_{0})>={{1}\over{T}}\int_0^T\dot \xi_{i} {\rm d}t\neq 0.$$ (14)

Before performing the integration, we change the time variable t to u in the integral yielding

$$F_{i}(\vec \xi_{0})=<\dot \xi_{i}>={{n}\over{2\pi j}}\int_0^{2\pi j} {\dot \xi_{i}\over \dot u}{\rm d}u$$ (15)

where the expression for the temporal derivative of the eccentric anomaly is given by

$$\dot u={n\over (1-e\cos u)}\cdot$$ (16)

The disturbing function and the dissipative terms in Eqs. (7) are calculated without any approximation in their expressions. The integrations of Eqs. (7), over a synodic period, are done numerically using the ODE integrator (Shampine & Gordon 1975).

If $\vec \xi_{j}=(a_{j}, e_{j}, \varpi_{j}, \Phi_{j})$ is an initial guess in the vicinity of an exact corotation point we define f_i by:

$$F_{i}(\vec \xi_{j})=f_{i} , i=1,...,4.$$ (17)

In this way, we use a four-dimensional iterative method to numerically solve the system of equations.

$$\frac{\partial F_{1}}{\partial a} \mid_{(\vec \xi_j)} \Delta a + \displaystyle\frac{\partial F_{1}}{\partial e} \mid_{(\vec \xi_j)} \Delta e$$

$$+ \displaystyle\frac{\partial F_{1}}{\partial \varpi} \mid_{(\vec \xi_j)} \Delta \varpi +\displaystyle\frac{\partial F_{1}}{\partial \Phi} \mid_{(\vec \xi_j)} \Delta \Phi=-f_{1},$$

$$\displaystyle\frac{\partial F_{2}}{\partial a} \mid_{(\vec \xi_j)} \Delta a + \displaystyle\frac{\partial F_{2}}{\partial e} \mid_{(\vec \xi_j)} \Delta e$$

$$+ \displaystyle\frac{\partial F_{2}}{\partial \varpi} \mid_{(\vec \xi_j)} \Delta \varpi +\displaystyle\frac{\partial F_{2}}{\partial \Phi} \mid_{(\vec \xi_j)} \Delta \Phi=-f_{2},$$ (18)

$$\displaystyle\frac{\partial F_{3}}{\partial a} \mid_{(\vec \xi_j)} \Delta a + \displaystyle\frac{\partial F_{3}}{\partial e} \mid_{(\vec \xi_j)} \Delta e$$

$$+ \displaystyle\frac{\partial F_{3}}{\partial \varpi} \mid_{(\vec \xi_j)} \Delta \varpi +\displaystyle\frac{\partial F_{3}}{\partial \Phi} \mid_{(\vec \xi_j)} \Delta \Phi=-f_{3},$$

$$\displaystyle\frac{\partial F_{4}}{\partial a} \mid_{(\vec \xi_j)} \Delta a + \displaystyle\frac{\partial F_{4}}{\partial e} \mid_{(\vec \xi_j)} \Delta e$$

$$+ \displaystyle\frac{\partial F_{4}}{\partial \varpi} \mid_{(\vec \xi_j)} \Delta \varpi +\displaystyle\frac{\partial F_{4}}{\partial \Phi} \mid_{(\vec \xi_j)} \Delta \Phi=-f_{4}.$$

The new solutions, $a_{j+1}=a_{j}+\Delta a$, $e_{j+1}=e_{j}+\Delta e$, $\varpi_{j+1}=\varpi_{j}+\Delta \varpi$ and $\Phi_{j+1}=\Phi_{j}+\Delta \Phi$ replace the previous ones and the procedure is iterated until $\mid { { \xi_{i_{j+1}}- \xi_{i_j}} \over { \xi_{i_{j+1}}}} \mid<10^{-\delta}$, for all i=1...4, where $\delta$ gives the accuracy for the determination of the corotation point. The derivatives of F_i introduce second derivatives in respect with $\Re$ which are computed in the variables $(a,e,\varpi,\Phi)$ as previously described and the result integrated in u over a synodic period (u and $u_{\rm P}$ are related through $\dot u_{\rm P}= {n\over n_{\rm P}}{1-e\cos u\over1-e_{\rm P}\cos u_{\rm P}}\dot u$).

By our method, we can also study the linear stability of a equilibrium point calculating the eigenvalues $\lambda_{1}, \lambda_{2}, \lambda_{3}$ and $\lambda_{4}$ of the characteristic equation ${\rm det}(J-\lambda I)=0$ where, J is the Jacobian matrix given by

$$J=\pmatrix{{\partial F_{1}\over \partial a}&{\partial F_{1}\o... ...\over \partial \varpi}&{\partial F_{4}\over \partial \Phi}\cr}$$ (19)

and I is the unit matrix of rank 4. If all real parts of the eigenvalues are negative the point is stable, but if at least one eigenvalue has a positive real part, the point is unstable.

It must be noted that in the first order method the dissipative terms $\dot \varpi_{ {\rm f.d.}}$ and $\dot \Phi_{ {\rm f.d.}}$ are dropped from the Eqs. (7) because they have null average.

To determine the corotation points we have to take an initial guess for the set of variables for a fixed C (drag coefficient). After finding this point (we can do numerical simulation of the N-body equations or use a first order analytical method to have this initial guess) we determine all other points by varying the parameter C. Also a value already obtained for a particular planet eccentricity is an initial guess for a little higher (or lower) planet eccentricity.

To check the validity and accuracy of the (developed) method, we start by determining the corotation points for the 2:3 resonance considering Stokes drag law for the dissipative force. For the perturber, we consider a Jupiter mass planet with $\mu=\mu_{\rm J}={M_{\rm P}\over M_{\odot}}={1\over 1047.355}$ and a fixed orbit with elements: $a_{\rm J}=5.2~{\rm AU}$, $e_{\rm J}=0.05$ and $\varpi_{\rm J}=0.0$. It can be observed in Fig. 1 that there is a good agreement between the solutions determined by the first order method and those obtained by numerical integrations of the complete (N-body) equations of motion using the RADAU 16th order integrator (Everhart 1985).

Figure 1: Orbital elements a, e, $\varpi$ and $\Phi$ versus C corresponding to the corotation points for the 2:3 resonance, for a Jupiter-like planet with $\mu _{\rm J} = 1/1047.355$, $e_{\rm p}=0.05$, $a_{\rm p}=5.2$ AU and $\varpi _{\rm p}=0$. Continuous line stands for the first order method and circles stand for numerical integrations. The dissipative force is Stokes drag

Figure 5: Corotation points for the 3:2 interior resonance for the case of v2drag with planet as in Fig. 1 except for $e_{\rm p}=0.4$

To obtain these points (indicated with dots) we perform an integration of the complete equations of motion until the planetesimal has remained in a stable resonant orbit for a long enough time (usually $\sim$10⁵ years). After that, we take an average for each variable considering only the variations of the elements for the part of the stable planetesimal's orbit. We see that these numeric points match only the stable branch of the solutions. We checked this by computing the eigenvalues of the Jacobian matrix (Eq. (19)), yielding a stable solution for the branch where numerical points were found and an unstable solution corresponding to the other branch. Next, we tried out the method for a high planetary eccentricity (e=0.1) considering the same mass and semi-major axis as for the planet as above. For this case, the corotation points both by our method and those obtained by numerical simulations are shown in Fig. 2, thus confirming that the first order method works also for high eccentricities. Next, we determined the corotation points for the 5:7 resonance, considering the same conditions as in the first case. We find, by comparing with numerical simulations, that there are high discrepancies between the values obtained by both methods (see Fig. 3). This came about because the 5:7 resonance for a Jupiter-like planet is an example of an extended corotation (mothé- Diniz & Gomes 2000). For this case, the developed method does not give good results due to its limitation as a first order theory in the planet's mass whereas extended corotation is strongly dependent on the planet's mass. This limitation brings about the deviations found in Fig. 3. Extended corotation resonance was first noted by Mothé-Diniz & Gomes (2000) in numerical studies and its main characteristic is that there is a continuous range of drag coefficients up to zero for which corotation resonance takes place, or, in terms of planetesimal size, there is no upper limit of planetesimal radius for which this resonance is effective.

As a last example for an exterior resonance we show in Fig. 4 the points for the 7:9 resonance where two planetary masses were considered: $\mu={1\over 3150.0}$ and $\mu ={1\over 4500}$ and the same planetary orbit as that used in the first example. It can be noticed that the variation of the mass yields only a small shift in the location of the points. In all the cases considered above we assumed Stokes drag as the dissipative force.

We next apply our method for interior resonances. Figure 5 shows the corotation points for the interior 3:2 resonance considering the v² gas drag case as the dissipative force and considering the perturber with mass: $\mu=\mu_{\rm J}$ and fixed orbit with elements: $a=5.2~{\rm AU}$, e=0.4 and $\varpi=0.0$. We notice that the 3:2 interior resonance is also a case of the extended corotation and that for a very low value of drag coefficient (which means large planetesimal radius) the equilibrium points tends to a solution without dissipation (Ferraz-Mello et al. 1993). We have found in our applications the extended case for interior resonances with planetary excentricities in a range from 0.1 to 0.6 so, we choose e=0.4 as mean value for the applications shown.