A&A 366, 699-707 (2001)
DOI: 10.1051/0004-6361:20000267

Resonance trapping of planetesimals subject to gas drag:
The location of corotation points

D. E. Santos Jr. - R. S. Gomes

Observatório Nacional, Rua General José Cristino, 77 20910-211 Rio de Janeiro RJ, Brazil

Received 14 August 2000 / Accepted 7 November 2000

Abstract
We present two methods to determine corotation equilibrium points for the elliptical, restricted three body problem with a dissipative force. We show that our first order method is valid for high planetary eccentricity and determine corotation points for interior resonances, an improvement in respect with previous analytical first order methods, however it fails to determine the equilibrium points for the case of the extended corotation resonances. We develop another generic method which determines the solutions for this case. This method gives the corotation points with extreme accuracy for all cases (exterior, interior and extended). We compare our results with numerical integrations of the complete N-body equations in order to determine the accuracy of the points determined by both methods.

Key words: solar system: formation - planets and satellites: general

1 Introduction

Resonance trapping occurs when the gain of energy and angular momentum of a particle due to resonant gravitational perturbations from a planet exactly compensates for the loss of energy and angular momentum due to an external dissipation (gas drag for instance). Since the discovery of this phenomenon for the case of planetesimals subject to gas drag by Weidenschilling & Davis (1985), several works have been done in an attempt to understand the dynamics as well as the implications of this process in the formation of the Solar System (see for instance, Beaugé & Ferraz-Mello 1993; Beaugé et al. 1994; Gomes 1995; Kary & Lissauer 1995). Weidenschilling & Davis suggested that fragments of collisions between planetesimals would pass through resonances and under the influence of gas drag be accreted by a growing planetary embryo. Patterson (1987) differently suggested that planetary formation could have processed at two-body exterior resonances due to capture and global accretion of planetesimals at resonance sites starting with Venus for the inner planets and Jupiter for the outer planets. Beaugé et al. (1994) did a numeric simulation envisaging the formation of a proto-Saturn by considering the joint effects of planetary perturbations, gas drag, mutual gravitation and collisions among a swarm of 1000 equal-radius planetesimals. They suggested that the mechanism proposed by Patterson could explain the quasi-commensurability of the planets for the outer solar system. During its evolution in inward spiraling orbits, a planetesimal can pass through all resonances and collide with the Sun or a planet, be sent away from the solar system by a close approach with the perturber or be trapped into a resonance with the planet. When a mean motion and a secular resonance occur simultaneously, we say that the body is trapped into a corotation resonance. Otherwise, a libration resonance takes place. There are important differences between the geometry of these two types of resonances. In the libration case, trapped planetesimals share the same average semi-major axis, but other orbital elements present a noticeable dispersal. This happens due to the circulation of the longitude of the perihelion. In other words, planetesimals trapped in libration resonance are distributed over all longitudes for a given time. On the other hand, equal-sizes planetesimals trapped in corotation resonances accumulate into a single point. In this case, the longitude of the perihelion librates around a fixed value (in respect with the perturber's perihelion direction). These points of accumulation have been understood as sites for accretion of small planetesimals during the process of planetary formation. In a recent work, Mothé-Diniz & Gomes (2000) studied the collisional dynamics of planetesimals subject to a gas drag dissipation. They observed that the relative collisional speed between the planetesimals at these points was very low, reinforcing the idea that corotation points are favorable places for accretion of small planetesimals.

There is usually a maximum radius of the planetesimal and a minimum one for which corotation trapping takes place. Above the maximum planetesimal radius, trappings always occur in libration resonance. Although this fact limits the range of applicability of corotation resonance trapping in planetesimals accretion, the very efficient mass accumulation process thus induced motivates the present study on a more comprehensive determination of corotation equilibrium points for any resonance, planetary orbit and mass. It must also be emphasized that the planetesimals sizes for which corotation resonances are most effective range from a few meters to a few kilometers (depending on the planet's mass, nebula's physical properties, etc.) Accretion at this size range is hardly explained through physical coalescence (valid for much smaller grains) or mutual gravitation (valid for large bodies). Moreover, this size range is that for which radial shift by drag is most effective, yielding a drag time scale of just a few thousand years (Lissauer 1993). A mechanism that can both halt the fast radial migration and provide a process of effective accretion may have had an important influence in the initial process of planetary formation.

The first work to systematically determine corotation points with gas drag was done by Beaugé & Ferraz-Mello (1993) who found these equilibrium points by developing a first order method based on mean equations, considering the case of a Stokes drag dissipation. There are however some limitations related to the first order averaging method used by Beaugé & Ferraz-Mello (1993). Due to an approximation in the classical Laplace expansion for the perturbing function, this method does not give good results for high (planetary or planetesimal) eccentricities. Also, they did not find these points for interior resonances. Before that, Weidenschilling & Davis (1985) claimed that trappings could only occur for exterior resonances. They explained that only for exterior resonance the effect of planetary resonance (move away the planetesimal semi-major axis from the planet) could counterbalance the effect of gas drag. Patterson (1987) also did not find them, but more recently, some numerical simulations have shown trappings of planetesimals for interior resonances (see for instance, Kary & Lissauer 1995). In particular, for the case of corotation resonances, these trappings are stable.

For libration points, in particular, Beaugé (1999) showed that some effects of short-period terms in the location of the equilibrium solutions that are not considered in the first order averaging method brings significant difference to the location of the points. More recently, Beaugé et al. (1999) developed a second order method to explain the discrepancy between the corotation equilibrium points given by the analytic methods and numerical simulations. The results obtained by this new model showed a considerable improvement as compared with the first order method. However, even using the second order method, Beaugé et al. (1999) did not found some types of trappings that are found in numerical simulations, as the extended corotation resonance (Mothé-Diniz & Gomes 2000) and trappings of particles into an interior resonance with the planet (Kary & Lissauer 1995; Mothé-Diniz & Gomes 2000).

In this paper, we use an iterative method to determine equilibrium points in corotation resonance. We show, by using a generic method (described in Sect. 3), that our results are in a very good agreement with those from numerical integration of the complete N-body equations. We also found, with the same method, the corotation points for all cases pointed out above. The paper is organized as follows. In Sect. 2, we present a first order method valid for any planet or planetesimal eccentricity. We show some results obtained by comparing the points given by the method with the ones given by numerical integration of the complete N-body equations. In Sect. 3, we present a generic method valid not only for any eccentricity, but also for any planetary mass. With this method we obtain accurate corotation points for both cases: extended resonances and interior resonances. We also vary some parameters as eccentricity and mass of the planet in order to know the effects of high eccentricities and large masses in the location of the points. Stability and instability of these points are also checked in order to know the limits of stable points. In the last section, we discuss the results obtained and draw some conclusions.

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
Open with DEXTER

	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$
Open with DEXTER

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.

3 Generic method

Our first order method, although giving a new result (corotation points for interior resonances) and being valid for high eccentricities of the planet and planetesimal, failed to determine the equilibrium points for the case of exterior extended resonances. Thus, in order to determine these correct solutions and also to generally improve the accuracy, we developed a generic method. We might opt for developing a second order method as done by Beaugé et al. (1999) but, as this would probably keep the discrepancies and would anyway spend a large computer time to solve the equations, we chose to develop a complete generic method which spends less time than the second order one.

For this method we do a complete (all variables varying point to point) numerical integration of the Eqs. (7), considering also the dissipative terms $\dot \varpi_{ {\rm f.d.}}$ and $\dot \Phi_{ {\rm f.d.}}$ in the equations. This demands a correction in the expression for the temporal derivative of the eccentric anomaly (given by Eq. (16)) to

$$\dot u={n+\dot \sigma+\dot e \sin u\over (1-e\cos u)}$$ (20)

where now $\dot l= n + \dot \sigma$.

Another important change for this method as compared with the first order method refers to the definition of F_i (Eq. (14)) and its derivatives (Eqs. (18)). In the first order method, F_i is defined as function of the average variables and the derivatives of F_i with respect to the variables $\left(\partial F_{i}\over \partial (a,e,\varpi,\Phi)\right)$ are obtained analytically but, for the generic method the complete integration of the Eqs. (7) requires a new definition for F_i and its derivatives. Now F_i is defined as a function of the initial variables, i.e., $F_{i}=(\xi_{i}(T)-\xi_{i}(0))/T$, where $\vec \xi(0)=(a_{0},e_{0},\varpi_{0},\Phi_{0})$ and we obtain the derivatives numerically, by taking this initial point $(a_{0},e_{0},\varpi_{0},\Phi_{0})$ and a small enough increment $\delta a$, $\delta e$, $\delta \varpi$ and $\delta \Phi$ for each variable in the neighborhood of this point. Thus, we calculate all the derivatives at this point numerically by

$${\partial F_{i}\over \partial a}\mid_{\vec \xi_{0}}= {{F_{i}(a_{0... ...,\varpi_{0},\Phi_{0})- F_{i}(a_{0},e_{0},\varpi_{0},\Phi_{0})} \over \delta a},$$

$${\partial F_{i}\over \partial e}\mid_{\vec \xi_{0}}= {{F_{i}(a_{0... ...,\varpi_{0},\Phi_{0})- F_{i}(a_{0},e_{0},\varpi_{0},\Phi_{0})} \over \delta e},$$ (21)

$${\partial F_{i}\over \partial \varpi}\mid_{\vec \xi_{0}}= {{F_{i}... ...varpi,\Phi_{0})- F_{i}(a_{0},e_{0},\varpi_{0}, \Phi_{0})} \over \delta \varpi},$$

$${\partial F_{i}\over \partial \Phi}\mid_{\vec \xi_{0}}= {{F_{i}(a... ...}+\delta \Phi)- F_{i}(a_{0},e_{0},\varpi_{0},\Phi_{0})} \over \delta \Phi}\cdot$$

In the first order method previously presented, we defined F_i as function of average variables thus the iterated convergence process is considered with respect to the average value for each variable but, for the generic method, F_i is defined as a function of initial variables and consequently the convergence process now is taken with respect to the initial value of each variable and when convergence is achieved we perform the average for each variable.

For this method, we first start with our basic example, calculating the equilibrium points for the exterior 2:3 resonance considering Stokes drag as the dissipative force and a planet with Jupiter mass $\mu=\mu_{\rm J}$ and a fixed orbit with elements $a=5.2~{\rm AU}$, e=0.05 and $\varpi=0.0$. Figure 6 shows the corotation points for this case. It can be observed that these points are in a better agreement than those obtained by the first order method (Fig. 1). Following the examples considered in the previous method, we determine the points for the 5:7 resonance. Comparing the result shown in Fig. 7 with the one obtained with the first order method (Fig. 3) we observe that this method gives the correct solutions and this is indeed a case of extended resonance.

Next, we determined the points for the 7:9 resonance considering now three planetary masses: $\mu={1\over 3150.0}$, $\mu={1\over 3500.0}$ and $\mu={1\over 4500.0}$ and the same orbit as considered above. It can be observed in Fig. 8 that this is also a case of the extended resonance for planet's mass equal to ${1\over 3500.0}$. Observing this plot, we also point out that the points found in Fig. 4 (obtained by the first order method) for planetary mass equal to $\mu={1\over 3150.0}$ are wrongly determined and that the extended resonance occurs only for limited values of planet's mass. We determined the points for another case of extended resonance. Figure 9 shows these points for the 4:7 resonance with also three planetary masses considered: $\mu={1\over 400.0}$, $\mu=\mu_{\rm J}$ and $\mu={1\over 2000.0}$ and the same orbit. We have not found any extended corotation for first order resonances even considering high planet's mass or eccentricity. For all examples above, Stokes drag was considered as the dissipative force.

For interior resonances considering this method we took both Stokes and v² gas drag laws as dissipative forces. Figure 10 shows the corotation points for the 3:2 resonance considering both Stokes and v² gas drag. Although we had stopped the numerical points in C=10^-6 (year^-1, AU^-1) they extended for all values below this point. We also determined the points for other interior resonances, for instance, 2:1 (Fig. 11), 5:2 (Fig. 12) and 5:3 (Fig. 13) each of them considering both dissipative forces. For all cases above we considered the perturber's mass: $\mu=\mu_{\rm J}$ and planetary orbit with elements: $a=5.2~{\rm AU}$, e=0.4 and $\varpi=0.0$.

We also study the effect of planet's mass and eccentricity variation in the location of the points. For the first case we consider the interior 3:2 resonance with perturber's masses: $\mu={1\over 10}\mu_{\rm J}$, $\mu={1\over 2}\mu_{\rm J}$, $2~M_{\rm J}$ and $4~\mu_{\rm J}$ and fixed orbit with elements: $a=5.2~{\rm AU}$, e=0.4 and $\varpi=0.0$. For the second case, we consider the same resonance and a planet with Jupiter's mass and the same orbit except for the eccentricity which is varying from 0.1 to 0.6. It can be observed in Fig. 14 that mass variation causes a small influence in the location of the points for small values of the drag coefficient $(C<10^{-5}~{\rm years}^{-1})$. On the other hand, Fig. 15 shows that the planet's eccentricity variation brings a large variation in the planetesimal's eccentricity for small values of drag coefficient $(C<10^{-4}~{\rm years}^{-1})$. Also as eccentricity increases, the range of drag coefficients that enable resonance trapping shrinks. Although it seems that for any non null planetary eccentricity there is some $\Delta C$ left that enables resonance trappings into these important first order resonances, for a low enough planetary eccentricity, this C range becomes too narrow for any practical meaning (note the logarithm scale of Fig. 15). Thus the choice of a high planetary eccentricity for the example in Fig. 10 was not incidental. It is not without reason that trappings into important interior resonances like 2:1 and 3:2 are not presented in previous works. If we assume a Jupiter-like planet, including its eccentricity equal to 0.05, these trappings are virtually impossible. However, considering a higher planetary eccentricity as here presented, trappings into these low order/high "j'' resonances are possible with a non negligible probability.

At last, we determined the points for two high j interior resonances 11:10 and 9:8 (Figs. 16 and 17), considering a planet with $\mu={1\over 1000}\mu_{\rm J}$ and fixed orbit with elements: $a=5.2~{\rm AU}$, e=0.05 and $\varpi=0.0$. We observe for these two examples that for high j the interior resonances are not extended anymore. These interior resonance trappings, for small planetary masses and modest eccentricity, have been found in previous works (Kary & Lissauer 1995). They take place for resonances closer to the planet (high j). For large planetary mass, this region is chaotic and no trapping can take place, however for small j (further from the planet), trapping is again possible but just for high planetary eccentricity as Fig. 15 suggests.

4 Conclusions

Two methods have been developed to determine corotation equilibrium points for the elliptical restricted 3-body problem. The simpler one is a first order (in the small parameter) method which has however no limitation for the planet or planetesimal eccentricities since it does not assume any truncation of the disturbing function. In this way, corotation points can be determined for any planetary eccentricity and any resonance and the results certainly show an improvement with respect to previous analytical methods (Beaugé & Ferraz-Mello). Moreover, corotation points for interior resonances are also determined by this first order method. Notwithstanding the improvement we obtain through this method, some corotation points found in numerical simulations cannot be reproduced by it. This is the case of the (exterior) extended corotation points (Mothé-Diniz & Gomes 2000) which are only obtained with higher order (in the small parameter) theory. A generic method, which is basically a numerical iterative algorithm, is thus developed to determine corotation points not only for any planetary orbit but also for any planetary mass. In this way, extended corotation points could be reproduced. It is noteworthy that both methods work for any nonconservative force acting on the massless body and we did tests with drag models both proportional to the velocity absolute value and to its square. By presenting both methods, we aimed at showing all the limitations imposed by a first order theory, although not truncated in the perturbing function. This point was not very clear from previous analytical works on this subject (Beaugé 1998; Beaugé et al. 1999). We must note however that the generic method must be in the future preferred to the first order one since, besides being much more accurate, it also presented a computing time comparable to that of the first order method.

Our results basically confirm that corotation resonances occur usually for a limited range of drag coefficients or, other things being constant, for a limited range of planetesimal sizes. For any resonance, given a non null planetary eccentricity, there is a minimum size for which trapping is possible and this will occur in the corotation regime. There is also a maximum size below which trappings always occur in the corotation regime. Above this limit, trapping will take place always in the libration regime. The range of sizes for which corotation trapping is possible will thus be a function of planetary eccentricity (increasing with it) and the specific resonance (usually increases for resonances closer to the planet). For a standard nebular model, planetesimals sizes for which corotation trapping occurs, considering a Jupiter-like planet (distance from the star and mass) extends from a few meters to near 1 km. It is interesting to note that for this size range gas drag is very efficient in inducing a fast orbital decay for planetesimals, thus a resonance trapping by a proto-planet can halt this fast decay. Moreover, because corotation resonance induces a fast accumulation of equal sized particles into common equilibrium points, with approximations taking place with very low relative speed, these points may work to produce larger planetesimals from smaller ones, which will then suffer a much slower drag decay.

For some resonances and planetary mass, the so called extended corotation resonance can take place. Its main characteristic is that there is no lower limit for drag rate (or upper limit for planetesimal size) below (above) which the regime is changed to libration. In this sense, the accretion process could be in principle carried on until a planetary embryo was formed. The problem is that this kind of resonance was only observed for some higher order resonances (4:7, 7:9), which are associated to a generally lower capture probability and also overlaps with more important lower order resonances inducing possible high velocity impacts among objects in different resonances.

Important (low j) first order interior resonance trappings have been found only for high planetary eccentricities with Jupiter sized planets. For modest (Jupiter-like) planetary eccentricity, only interior resonance closer to the planet (high j) are possible. In this case, planetary mass must be increasingly lower so as to avoid chaos. In all cases, however, trapping probability was found in the numerical examples to be very low, what limits the applicability of interior resonance trapping for mass accretion.

The importance of corotation resonance trapping seems thus limited in stopping the decay and promoting the accretion of meter-1 km sized planetesimals in exterior resonances with protoplanets in order to form larger objects that are subject to strong enough mutual gravitation forces to help in their further accretion process. In this sense, this mechanism is particularly useful if we consider the recently ressurected gravitational instability model (Boss 1997) for the formation of Jupiter and possibly Saturn. If the first planet is to be formed by the more conventional accretion theory no resonance trapping can be claimed to start meter sized planetesimals accretion. Now one question naturally rises: what is the fate of planetesimals that would be formed at corotation points and after gaining mass pass to the libration regime? For this case, there is no accumulation point. Since trapped planetesimals have higher eccentricities than nontrapped ones, a natural consequence, first noted by Weidenschilling & Davis (1985), would be the high velocity disrupting impacts of trapped planetesimals producing smaller fragments that would naturally escape resonance (larger drag rate) and eventually accrete to the planet. Yet it is not obvious that high eccentricities for planetesimals trapped in (the same) resonance will generally produce high velocity collisions. Because orbits of trapped planetesimals cannot be considered random, this point deserves a deeper investigation, which we are presently addressing. In advance it seems now clear that collisions between bodies of similar sizes (suffering similar drag forces) will occur in relatively low velocities, increasing however for higher planetary eccentricity.

The methods developed in this paper helped us to determine precisely the real limits of influence of a corotation resonance trapping. Our main goal, for which these methods work as a tool, is to determine the real influence of resonance trapping of planetesimals with proto-planets in a forming planetary system. Does this phenomenon favor accretion or hamper it? We plan to have a fairly precise answer to this question with our ongoing research.

Acknowledgements

D. E. Santos Jr. is grateful to CAPES for the scholarship.

References

 
• Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Prog. Theor. Phys., 56, 1756 In the text NASA ADS  
• Beaugé, C., & Ferraz-Mello, S. 1993, Icarus, 103, 301 In the text NASA ADS  
• Beaugé, C., Aarseth, S. J., & Ferraz-Mello, S. 1994, MNRAS, 270, 21 In the text NASA ADS  
• Beaugé, C. 1998, Celest. Mech. Dynam., Astron., 68, 283 In the text  
• Beaugé, C., Lemaître, A., & Jancart, S. 1999, Planet. Space Sci., 47(5), 643 In the text NASA ADS  
• Boss, A. 1997, Science, 276, 1836 In the text NASA ADS  
• Brouwer, D., & Clemence, G. M. 1961, Methods of Celestial Mechanics (Academic Press, New York) In the text  
• Everhart, E. 1985, in Dynamics of Comets: Their Origin and Evolution, ed. A. Carusi & G. B. Valsech, 195-202 In the text  
• Ferraz-Mello, S., Tsuchida, M., & Klafke, J. C. 1993, Celest. Mech. Dynam. Astron., 55, 25 In the text NASA ADS  
• Gomes, R. S. 1995, Icarus, 115, 47 In the text NASA ADS  
• Kary, D. M., & Lissauer, J. J. 1995, Icarus, 117, 1 In the text NASA ADS  
• Lissauer, J. J. 1993, ARA&A, 31, 129 In the text NASA ADS  
• Mothé-Diniz, T., & Gomes, R. S. 2000, Planet. Space Sci., 48, 305 In the text  
• Patterson, C. W. 1987, Icarus, 70, 319 In the text NASA ADS  
• Shampine, L. F., & Gordon, M. K. 1975, Computer Solution of Ordinary Differential Equations: the Initial Value Problem (W. H. Freeman and Co., San Francisco) In the text  
• Weidenschilling, S. J. 1977, MNRAS, 180, 57 In the text NASA ADS  
• Weidenschilling, S. J., & Davies, D. R. 1985, Icarus, 62, 16 In the text NASA ADS

 

Online Material

	Figure 2: Same as Fig. 1 for planet's eccentricity equal to 0.1
Open with DEXTER

	Figure 3: Same as Fig. 1 for the 5:7 resonance
Open with DEXTER

	Figure 4: Same as Fig. 1 for the 7:9 resonance and different planetary masses. Dashes stand for planet's mass equal to $\mu ={1\over 3150}$ and dots stand for planet's mass equal to $\mu ={1\over 4500}$
Open with DEXTER

	Figure 6: Same as Fig. 1, with continuous line standing for the generic method
Open with DEXTER

	Figure 7: Same as Fig. 6, for the 5:7 extended resonance
Open with DEXTER

	Figure 8: Corotation points for the 7:9 resonance, for the generic method, planet as in Fig. 1, except for planetary masses. Continuous line stands for planet's mass equal to $\mu ={1\over 400}$, dashes stand for planet's mass equal to $\mu ={1\over 3500}$ (extended corotation) and and dots stand for planet's mass equal to $\mu ={1\over 4500}$
Open with DEXTER

	Figure 9: Corotation points for the 4:7 resonance, for the generic method, planet as in Fig. 1, except for planetary masses. Dots stand for planet's mass equal to $\mu ={1\over 400}$ (extended corotation), continuous line stands for planet's mass equal to $\mu ={1\over 1047.355}$ and dashes stand for planet's mass equal to $\mu ={1\over 2000}$
Open with DEXTER

	Figure 10: Corotation points for the 3:2 resonance, for the generic method, planet as in Fig. 5. Continuous line and circles stand for points determined considering Stokes drag (generic and numerical methods) while dashes stand for v2 drag
Open with DEXTER

	Figure 11: Corotation points for the 2:1 resonance, planet as in Fig. 5 and symbols as in Fig. 10
Open with DEXTER

	Figure 12: Corotation points for the 5:2 resonance, planet as in Fig. 5 and symbols as in Fig. 10
Open with DEXTER

	Figure 13: Corotation points for the 5:3 resonance, planet as in Fig. 5 and symbols as in Fig. 10
Open with DEXTER

	Figure 14: Variation of corotation points for the case of Stokes drag, with planet's mass for the 3:2 resonance, generic method with planet, except for masses, as in Fig. 5
Open with DEXTER

	Figure 15: Variation of corotation points for the case of Stokes drag, with planet's eccentricity, for the 3:2 resonance, generic method with planet, except for eccentricities, as in Fig. 5
Open with DEXTER

	Figure 16: Corotation points for the 11:10 interior resonance, Stokes drag, generic method, with $\mu_{\rm p} = 1/1000~\mu_{\rm J}$, and planetary orbital elements as in Fig. 1. Continuous line stands for the generic method and points stand for numerical integrations
Open with DEXTER

	Figure 17: Same as Fig. 16 for the 9:8 resonance, only generic method
Open with DEXTER