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

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

$\displaystyle {\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},$
$\displaystyle {\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)
$\displaystyle {\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},$
$\displaystyle {\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$

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.

3 Generic method