%This document is typeset in REVTeX
%\documentstyle[preprint,aps,epsf]{revtex}
\documentstyle[aps,epsf]{revtex}
% To force arabic section heading in REVTeX style we use the
% following command.
\renewcommand{\thesection}{\arabic{section}}
\input amssym.def
\newtheorem{thm}{Theorem}[section]
\newtheorem{prop}[thm]{Proposition}
\newtheorem{hyp}[thm]{Hypothesis}
\newtheorem{cor}[thm]{Corollary}
\newtheorem{ques}[thm]{Question}
\newtheorem{lem}[thm]{Lemma}
\newtheorem{res}[thm]{Result}
\newtheorem{defn}[thm]{Definition}
\newtheorem{remark}[thm]{Remark}
\newtheorem{conj}[thm]{Conjecture}
%%%%%%%%%%%%%%FONTS%%%%%%%%%%%%%%%%%%%%%%%
\newenvironment{pf}{{\it Proof}.}{$\Box$}
\begin{document}
\date{\today}
\title{Evolutionary Dynamics While Trapped in Resonance:
A Keplerian Binary System Perturbed by
Gravitational Radiation}
\author{
C. Chicone\thanks{ Department of Mathematics,
University of Missouri, Columbia, MO 65211.
Supported by
the NSF grant DMS-9531811.} ,\quad
B. Mashhoon\thanks{ Department of Physics and Astronomy,
University of Missouri, Columbia, MO 65211.
} \quad
and D. G. Retzloff\thanks{ Department of Chemical Engineering,
University of Missouri, Columbia, MO 65211.}
}
\maketitle
\begin{abstract}
The method of averaging is used to investigate the
phenomenon of capture into resonance for a model that describes
a Keplerian
binary system influenced by radiation damping and external normally
incident periodic gravitational radiation.
The dynamical evolution of the binary orbit while
trapped in resonance is elucidated using the second order partially averaged
system.
This method provides a theoretical framework that can be used to explain the
main evolutionary dynamics of a physical system that has
been trapped in resonance.
\end{abstract}
\section{Introduction}
In two previous papers~\cite{cmr,cmr1}, we considered the long term nonlinear
evolution of a Keplerian binary system that is perturbed by a normally incident
periodic gravitational wave, and in a recent work~\cite{cmr2} we considered the
additional effect of radiation damping,
which is of interest in connection with
the observed behavior of the Hulse-Taylor
binary pulsar PSR~1913+16~\cite{hulse,taylor}.
These studies have been concerned with the issue
of {\em gravitational ionization}, i.e. the possibility that an external
periodic gravitational wave could ionize a Keplerian binary system over a long
period of time. The impetus for this subject has come from the close conceptual
analogy between gravitational ionization and a fundamental physical process,
namely, the electromagnetic ionization of a hydrogen atom. That is, in these
studies one hopes to learn about the disposition of gravitational radiation to
transfer energy and angular momentum in its interaction with matter. In our
recent investigation~\cite{cmr2}, we encountered an interesting dynamical
phenomenon connected with the passage of the binary orbit through resonance. As
the binary system loses energy and angular momentum by emitting gravitational
waves, its semimajor axis and eccentricity decrease monotonically on the
average; however, this process of inward spiraling could stop if the system is
captured into a resonance. The resonance condition fixes the semimajor axis;
therefore, if the semimajor axis decreases to the resonant value and the orbit
is trapped, it then maintains this value {\it on the average} while the external
perturbation deposits energy into the system to compensate for the radiation
damping {\it on the average}.
It turns out that along with this energy deposition, the external
tidal perturbation can also deposit angular momentum into the binary orbit so that
its eccentricity decreases considerably during the passage through resonance.
This was the situation for the particular instance of resonance trapping
reported in~\cite{cmr2}.
In general, the orbital angular momentum can increase or decrease while the orbit
is trapped in resonance.
Figure~\ref{fig1} depicts passage through a
$(1:1)$ resonance in which the eccentricity decreases.
A similar phenomenon but with an increase in eccentricity
has been reported in the recent numerical study
of the three-body problem by Melita and Woolfson~\cite{melita}.
The same situation can occur over a long time-scale in a $(4:1)$ resonance in our
model as depicted in Figure~\ref{fig3}.
The dynamical phenomena associated with an orbit trapped in a resonance
occur over many Keplerian periods of the osculating ellipse;
therefore, it is natural to average the dynamical equations over the
orbital period at resonance. This partial averaging removes the ``fast'' motion
and allows us to see more clearly the ``slow'' motion during trapping.
It is possible
to provide a description of the slow motion as well as
a theoretical justification for the transfer of angular momentum by
means of the method of second order partial averaging near a resonance. The
purpose of the present paper is to study this phenomenon theoretically using the
method of averaging for the resonances that occur in a Keplerian binary that is
perturbed by the emission and absorption of gravitational radiation.
Let us consider the simplest model involving a perturbed Keplerian (i.e.
nonrelativistic) binary that contains the effects of radiation reaction damping
and external tidal perturbation in the lowest (i.e. quadrupole) order. The
tidal perturbation could be due to external masses and gravitational radiation;
in fact, we choose in what follows a normally incident periodic gravitational
wave for the sake of simplicity \cite{cmr2}. The Keplerian binary orbit under
the combined perturbations due to the emission and absorption of gravitational
radiation in the quadrupole approximation is given in our model
by the equation of relative
motion
\begin{equation}\label{RelativeMotion}
\frac{d^2x^i}{dt^2} + \frac{kx^i}{r^3} + {\cal R}^i =
- \epsilon {\cal K}_{ij}\; x^j,
\end{equation}
where ${\bf x}(t)$ is the relative two-body orbit, $r = |{\bf x}|$, and $k =
G_0(m_1+m_2)$. The binary consists of two point masses $m_1$ and $m_2$,
$m_1+m_2 = M$, that move according to ${\bf x}_1(t) = (m_2/M)\; {\bf x}(t)$ and
${\bf x}_2(t) = -(m_1/M)\; {\bf x}(t)$
, so that the center of mass of the system is at the origin of the
coordinates. In fact, the center of mass of the binary can be taken to be at
rest in the approximation under consideration here. The explicit expressions
for ${\cal R}$ and ${\cal K}$ are
\begin{equation}\label{RadiationReaction}
{\cal R} = \frac{4G_0^2m_1m_2}{5c^5r^3}
\left[\left(12 v^2-30\dot{r}^2-\frac{4k}{r}\right){\bf v}
-\frac{\dot{r}}{r}\left(36v^2-50\dot{r}^2+\frac{4k}{3r}\right){\bf x}\right],
\end{equation}
and
\begin{equation}\label{TidalForce}
{\cal K} = \Omega^2 \left[ \begin{array}{ccc}
\alpha \cos{\Omega t} & \beta \cos{(\Omega t + \rho)} & 0 \\
\beta \cos{(\Omega t + \rho)} & -\alpha \cos{\Omega t} & 0 \\
0 & 0 & 0
\end{array} \right],
\end{equation}
where ${\bf v}=\dot{\bf x}(t)$, an overdot denotes differentiation with respect
to time, $\alpha$ and $\beta$ are of the order of unity and are the amplitudes
of the two independent states of linear polarization of the normally incident wave, $\rho$
is a constant phase, and $\Omega$ is the frequency of the external wave.
It is interesting to transform the dynamical
system~(\ref{RelativeMotion})--(\ref{TidalForce}) to dimensionless form.
To accomplish this, let ${\bf x}\rightarrow R_0 {\bf x}$ and $t \rightarrow T_0 t$
where $R_0$ and $T_0$ are
scale parameters. Under this transformation, $k \rightarrow k T_0^2/R_0^3$ and
${\cal K}$ remains unchanged if we let $\Omega \rightarrow \Omega/T_0$. Let us
further restrict $R_0$ and $T_0$ by the relation $k T_0^2 = R_0^3$, so that we
can set $k=1$ in the dynamical equations; for instance, this condition is
satisfied if $R_0$ is the initial semimajor axis of the unperturbed Keplerian
orbit and $2\pi T_0$ is its period. Furthermore,
the dynamical system~(\ref{RelativeMotion})--(\ref{TidalForce}) is planar;
therefore, it is
convenient to express these dimensionless equations in polar coordinates
$(r, \theta)$ in the orbital plane. The result is
\begin{eqnarray}\label{MathEQM}
\frac{dr}{dt} & = & p_r, \nonumber \\
\frac{d\theta}{dt} & = & \frac{p_{\theta}}{r^2}, \nonumber \\
\frac{dp_r}{dt} & = & -\frac{1}{r^2} + \frac{p_{\theta}^2}{r^3} +
\frac{4\delta p_r}{r^3}\left(p_r^2+6\frac{p_{\theta}^2}{r^2}+
\frac{4}{3r}\right) \nonumber \\
& & \hspace*{.25in} -\epsilon r
\Omega^2 [\alpha\cos{2\theta}\cos{\Omega t}
+\beta\sin{2\theta}\cos{(\Omega t+\rho)}], \nonumber \\
\frac{dp_{\theta}}{dt} & = & \frac{2\delta p_{\theta}}{r^3}
\left( 9p_r^2 - 6 \frac{p_{\theta}^2}{r^2} + \frac{2}{r}\right) \nonumber \\
&& \hspace*{.25in} +\epsilon r^2 \Omega^2 [\alpha\sin{2\theta}\cos{\Omega t}
-\beta\cos{2\theta}\cos{(\Omega t+\rho)}],
\end{eqnarray}
where $\delta$, $0<\delta<<1$,
is the dimensionless strength of radiation reaction and is given by
\begin{equation}\label{DeltaDefn}
\delta = \frac{4G_0^2m_1m_2}{5c^5T_0R_0},
\end{equation}
while $\epsilon$, $0<\epsilon<<1$, is the dimensionless strength of the
external periodic perturbation.
In this paper, we let $\delta = \epsilon \Delta$, where $\Delta$,
$0 < \Delta < \infty$, is a parameter that is fixed in the system;
in this way, we avoid dealing with a two parameter $(\epsilon,\delta)$
perturbation problem. In particular, we consider only perturbations that
correspond to fixed directions from the origin of this parameter space.
The full two parameter problem would require the consideration of perturbations
corresponding to all curves in the parameter space.
In the absence of radiative perturbations $(\epsilon =0$ and $\delta=0)$, the
dynamical system (\ref{MathEQM}) describes a Keplerian ellipse. It is therefore
useful to
express the dynamical equations (\ref{MathEQM}) in terms of Delaunay variables
that are closely related to the elements of the osculating ellipse. This is the
ellipse that the relative orbit would describe at time $t$, if the perturbations
were turned off at $t$.
The osculating ellipse always has the same focus,
which is taken to be the
origin of the (polar) coordinates in the space of relative coordinates.
Let the state of relative motion be described by $({\bf
x}, {\bf v})$, or equivalently $(r, \theta, p_r, p_{\theta})$, at time $t$; then,
the energy of the motion fixes the semimajor axis $a$ of the osculating ellipse
and its eccentricity is subsequently fixed by the orbital angular momentum
$p_{\theta}$. Only two angles are left to determine the osculating ellipse
completely: the orientation of the ellipse in the orbital plane given by $g$ and
the position on the ellipse measured from the periastron given by the true
anomaly $\hat{v}$. The latter is obtained from
$p_rp_\theta=e\sin\hat{v}$ and the equation of the
ellipse
\[ \frac{p_{\theta}^2}{r} = 1+e \cos{\hat{v}}, \]
and the former from
$\theta - \hat{v}$. The relevant Delaunay ``action-angle" variables $(L, G,
\ell, g)$ are thus defined by \cite{Kov,Stern}
\begin{eqnarray}
L := a^{1/2}, & \hspace*{.25in} & G := p_{\theta} = L
(1-e^2)^{1/2}, \nonumber \\
\ell := \hat{u} - e \sin{\hat{u}}, & & g := \theta - \hat{v},
\end{eqnarray}
where $\hat{u}$ is the eccentric anomaly of the osculating ellipse, $r =
a(1-e\cos{\hat{u}})$, and $\ell$ is the mean anomaly.
The dynamical system~(\ref{MathEQM}) in terms of Delaunay variables is given
briefly in Appendix~\ref{appendixa} and used in the following section.
The Delaunay equations of motion are useful for the investigation of periodic
orbits using the Poincar\'{e} surface of section technique~\cite{poincare}. It
has been shown in \cite{cmr2} that nearly resonant periodic orbits exist in system
(\ref{MathEQM}) for sufficiently small $\epsilon$ and $\Delta$. These
correspond to
$(m:1)$ resonances, where $(m:n)$ refers to the resonance condition $m\omega =
n\Omega$. Here $m$ and $n$ are two relatively prime integers and $\omega=1/L^3$
is the Keplerian frequency of the orbit. A linear perturbation treatment
\cite{mashoon1} first revealed resonant absorption at $(m:1)$ resonances. There
could, in principle, be other periodic orbits whose existence is not revealed by
our method \cite{ccc,cmr}.
In our numerical investigation of the simple nonlinear model described above,
we found~\cite{cmr2} instances of resonance trapping during which the
behavior of the osculating orbit could not be inferred in a simple manner
on the basis of equation~(\ref{MathEQM}).
However, the dynamics of the {\em averaged} equations
in resonance is simpler to analyze and it turns out that our numerical
results~\cite{cmr2} can be adequately explained using the second order partially
averaged dynamics. The phenomenon of resonance trapping appears to be
of basic significance for the origin of the solar system; therefore,
it is worthwhile to develop a general theoretical framework for the
study of the evolutionary dynamics while trapped in resonance.
The {\em dynamics} of a system when it is locked in
resonance is interesting in any circumstance involving
more than one degree of freedom; for instance, let us suppose that
the resonance condition
fixes an action variable---say, energy---and for a one dimensional
motion this would then imply that the state of the system at resonance is definite.
However, if other action variables are present, they will not necessarily remain
fixed while the system is trapped in resonance.
Instead, the state of the
system will in general vary and its dynamics at resonance is best
investigated using the method of averaging. This is a generalization
of the simple procedure that is commonly
employed in Hamiltonian dynamics: The Hamiltonian is averaged over certain
``fast'' variables and the resulting averaged Hamiltonian is
used to derive new dynamical equations that presumably describe the ``slow''
motion in a certain averaged sense. The general method is described in
Appendix~\ref{appendixb}, and it is applied to the dissipative
dynamical system under consideration here in the rest of this paper.
Resonance is a general and significant physical phenomenon and the description of the
state of a physical system while trapped in resonance is of intrinsic importance. The inherent
dynamics at resonance is trivial for a one dimensional oscillator, but is rich in
physical consequences for higher dimensional systems.
While the general methods
described here could be applied to a wide variety of physical problems, we confine our
attention to a single model.
Our results may, however,
be of qualitative interest in dealing with the three-body
problems that arise in the discussions of the origin of the structure
in the solar system.
The present paper relies on the results of our recent work \cite{cmr2}. We have
repeated here only what is needed for the discussion of the dynamics at
resonance; for a more precise and complete presentation of the
background material, our papers \cite{cmr,cmr1,cmr2} should be consulted.
Finally, a basic limitation of our model
should be noted. The only damping mechanism that we take into account is
gravitational radiation reaction; this is
consistent with our assumption that the binary consists of Newtonian point
masses moving in vacuum except for the presence of background gravitational
radiation. In this model, a theorem of Neishtadt can be used to show that
resonance trapping is a rather rare phenomenon~\cite{cmr2}. However, taking due
account of the finite size and structure of astronomical bodies and the existence
of an ambient medium, we would have to include in our model---among other things---the
additional damping effects of tidal friction as well as the various drag effects of
the ambient medium and electromagnetic radiation
(cf., for instance,~\cite{melita,henrard,beauge,gomes,lai}).
These additional frictional effects could well combine to violate the
condition $N$ of Neishtadt's Theorem~\cite{cmr2,arnold2}.
Thus, resonance trapping
may not be so rare in astrophysics after all~\cite{bj}. Inclusion of these additional
effects is beyond the scope of our work.
\section{Partial Averaging Near A Resonance}
We will consider the dynamics of the model system~(\ref{MathEQM}) that is
derived in~\cite{cmr2}.
It has the following form when expressed in the Delaunay elements for the
Kepler problem under consideration (cf. Appendix~\ref{appendixa}):
\begin{eqnarray}\label{AAAveraging}
\dot{L} & = & -\epsilon
\frac{\partial {\cal H}_{\mbox{\scriptsize ext}}}{\partial \ell} +
\epsilon\Delta f_L, \nonumber \\
\dot{G} & = & -\epsilon
\frac{\partial {\cal H}_{\mbox{\scriptsize ext}}}{\partial g} +
\epsilon\Delta f_G \nonumber \\
\dot{\ell}& = &\quad \frac{1}{L^3}
+ \epsilon \frac{\partial{\cal H}_{\mbox{\scriptsize ext}}}{\partial L}
+\epsilon\Delta f_{\ell},\nonumber \\
\dot{g} & = & \quad\epsilon
\frac{\partial {\cal H}_{\mbox{\scriptsize ext}}}{\partial G} +
\epsilon\Delta f_g,
\end{eqnarray}
where $\epsilon{\cal H}_{\mbox{\scriptsize ext}}$ is the Hamiltonian
corresponding to the external perturbation and
\begin{equation}\label{Hext}
{\cal H}_{\mbox{\scriptsize ext}}=\frac{1}{2}\Omega^2\left[
\alpha {\cal C}(L,G,\ell,g)\cos{\Omega t}+\beta {\cal S}(L,G,\ell,g)\cos(\Omega t+\rho)
\right].
\end{equation}
Here
\begin{eqnarray}\label{CSfs}
{\cal C}(L,G,\ell,g) & = & \frac{5}{2}a^2e^2\cos{2g}
+a^2\sum_{\nu =1}^\infty (A_\nu\cos{2g}\cos {\nu\ell}
-B_\nu\sin{2g}\sin{\nu \ell}), \nonumber \\
{\cal S}(L,G,\ell,g) & = &\frac{5}{2}a^2e^2\sin{2g}
+a^2\sum_{\nu =1}^\infty (A_\nu\sin{2g}\cos {\nu\ell}
+B_\nu\cos{2g}\sin{\nu \ell}), \nonumber \\
A_\nu & = &\frac{4}{\nu^2e^2}
\big[2\nu e(1-e^2)J_\nu'(\nu e)-(2-e^2)J_\nu(\nu e)\big], \nonumber \\
B_\nu & = & -\frac{8}{\nu^2e^2}(1-e^2)^{1/2}\,
\big[e J_\nu'(\nu e)-\nu (1-e^2)J_\nu(\nu e)\big],
\end{eqnarray}
$J_\nu$ is the Bessel function of the first kind of order $\nu$,
and $e=(L^2-G^2)^{1/2}/L$.
The radiation reaction ``forces" $f_D$, $D\in\{L,G,\ell,g\}$, are certain
complicated functions
of the Delaunay variables given in Appendix~\ref{appendixa}.
In fact, we will only require the
averages of the $f_D$ given by
\[\bar f_{D}:=\frac{1}{2\pi}\int_0^{2\pi} f_{D}(L, G, \ell, g) \,
d\ell.\]
These have been computed in~\cite{cmr2}, and are given by
\begin{eqnarray}\label{fqAve3}
\bar f_{L}&=&-\frac{1}{G^7}\big( 8+\frac{73}{3} e^2+ \frac{37}{12} e^4 \big),
\nonumber \\
\bar f_{G}&=&-\frac{1}{L^3G^4}(8+7e^2),\nonumber \\
\bar f_{\ell}&=&0, \qquad \bar f_{g}=0.
\end{eqnarray}
In order to study the dynamics of the system~(\ref{AAAveraging}) at resonance,
we will apply the method of averaging. We note here that averaging
over the fast angle $\ell$ and the time $t$ gives the correct approximate
dynamics for most initial conditions via Neishtadt's Theorem as explained
in our previous paper~\cite{cmr2}. Here we are interested in the orbits
that are captured into resonance. To study them, we consider partial
averaging at each resonance.
Let us fix the value of $L$ at the $(m:n)$ resonance, say
$L=L_*$ with $m/L_*^3=n\Omega$, and consider the deviation of $L$ from
the resonance manifold.
To measure this deviation, we introduce a new variable ${\cal D}$ given by
\begin{equation}\label{coordtrans}
L = L_* + \epsilon^{1/2}\; {\cal D}
\end{equation}
and a new angular variable $\varphi$ by
\[\ell=\frac{1}{{L_*}^3}\: t+\varphi=\frac{n\Omega}{m}t+\varphi.\]
The scale factor $\epsilon^{1/2}$ ensures that, after changing to
the new variables in~(\ref{AAAveraging}), the resulting equations
for $\dot {\cal D}$ and $\dot \varphi$ have the same order in
the new small parameter $\epsilon^{1/2}$ and therefore the system is in
the correct form for averaging.
These new coordinates are standard choices in the mathematical literature
(for more details see, for example,~\cite{arnold2} or~\cite{wig}).
It is important to emphasize here that the small parameter in the
actual dynamics is $\epsilon$; however, the small parameter turns out to
be $\epsilon^{1/2}$ in this case for the averaged dynamics.
To effect the coordinate transformation, we use the expansion
\begin{equation}\label{tayL}
\frac{1}{{L}^3}=\frac{1}{{L_*}^3}\Big[1-\epsilon^{1/2}\frac{3{\cal D}}{L_*}
+\epsilon\frac{6{\cal D}^2}{L_*^2} +O(\epsilon^{3/2})\Big]
\end{equation}
and find
\begin{eqnarray}\label{AA1Averaging}
\dot{{\cal D}} & = & -\epsilon^{1/2} F_{11}-\epsilon {\cal D} F_{12}+O(\epsilon^{3/2}),
\nonumber \\
\dot{G} & = & -\epsilon F_{22}+O(\epsilon^{3/2}), \nonumber \\
\dot{\varphi} & = & -\epsilon^{1/2}\frac{3}{L_*^4}{\cal D}
+ \epsilon \big( \frac{6}{L_*^5}{\cal D}^2+F_{32}\big)+O(\epsilon^{3/2}),\nonumber \\
\dot{g} & = & \quad\epsilon F_{42}+O(\epsilon^{3/2}),
\end{eqnarray}
where the
$F_{ij}(G, n\Omega t/m+\varphi,g,t)$
are defined such that the first index refers to the equation in which
it appears and the second index refers to the perturbation order in powers
of $\epsilon^{1/2}$ that is employed.
These quantities are given by
\begin{eqnarray}\label{Fij}
F_{11}&:=&\frac{\partial{\cal H}_{\mbox{\scriptsize ext}}}{\partial
\ell}(L_*,G,\frac{n\Omega}{m}t+\varphi,g,t)
-\Delta f_L(L_*,G,\frac{n\Omega}{m}t+\varphi,g), \nonumber \\
F_{12}&:=&\frac{\partial^2{\cal H}_{\mbox{\scriptsize ext}}}{\partial L\partial \ell}
(L_*,G,\frac{n\Omega}{m}t+\varphi,g,t)
-\Delta \frac{\partial f_L}{\partial
L}(L_*,G,\frac{n\Omega}{m}t+\varphi,g),\nonumber \\
F_{22}&:=&\frac{\partial{\cal H}_{\mbox{\scriptsize ext}}}{\partial
g}(L_*,G,\frac{n\Omega}{m}t+\varphi,g,t)
-\Delta f_G(L_*,G,\frac{n\Omega}{m}t+\varphi,g),\nonumber \\
F_{32}&:=&\frac{\partial{\cal H}_{\mbox{\scriptsize ext}}}
{\partial L}(L_*,G,\frac{n\Omega}{m}t+\varphi,g,t)
+\Delta f_\ell(L_*,G,\frac{n\Omega}{m}t+\varphi,g), \nonumber\\
F_{42}&:=&\frac{\partial{\cal H}_{\mbox{\scriptsize ext}}}{\partial
G}(L_*,G,\frac{n\Omega}{m}t+\varphi,g,t)
+\Delta f_g(L_*,G,\frac{n\Omega}{m}t+\varphi,g).
\end{eqnarray}
The system~(\ref{AA1Averaging}) is $2\pi m/\Omega$ periodic in
the temporal variable---since $\cal H_{{\mbox{\scriptsize ext}}}$ is
$2\pi /\Omega$ periodic in time---and is in time-periodic standard form.
Anticipating our intention to average to second order, we will apply an
averaging transformation (for a detailed exposition see~\cite{wig}).
It is the characteristic property of this transformation that it automatically
renders system~(\ref{AA1Averaging}) in a form such that to lowest order
the new system is exactly the first order averaged system and the second order
averaged system can be simply obtained by averaging the new
system (cf. Appendix~\ref{appendixb}).
To obtain the desired transformation, we define
\[
\bar{F}_{ij}:= \frac{\Omega}{2\pi m}
\int_0^{2\pi m/\Omega} F_{ij}(G,\frac{n\Omega}{m}s+\varphi,g,s)\,ds,
\]
and the deviation from the mean for $F_{11}$
\begin{equation}\label{avertran}
\lambda(G,\varphi,g,t):=F_{11}(G,\frac{n\Omega}{m}t+\varphi,g,t)-\bar{F}_{11}.
\end{equation}
Furthermore, we define $\Lambda(G,\varphi,g,t)$ to be the antiderivative of
$\lambda(G,\varphi,g,t)$ with respect to $t$ with the additional property
that the average of $\Lambda$ should vanish, i.e.
\[
\int_0^{2\pi m/\Omega} \Lambda(G,\varphi,g,s)\,ds=0.
\]
Moreover, we note that
both $\lambda$ and $\partial\Lambda/\partial\varphi$ have
zero averages. Our averaging transformation is given by
\[
{\cal D}=\widehat {{\cal D}}
-\epsilon^{1/2}\Lambda(\widehat{G},\widehat{\varphi},\widehat{g},t),\quad
G=\widehat{G},\quad \varphi=\widehat{\varphi},\quad g=\widehat{g}.
\]
The averaging transformation is constructed so that its average becomes
the identity transformation.
Let us observe that if $G$, $\varphi$, and $g$ depend on $t$ as solutions of the
system~(\ref{AA1Averaging}), then
\[
\dot\Lambda=\lambda-\epsilon^{1/2}\Big(\frac{3{\cal D}}{L_*^4}\Big)
\frac{\partial\Lambda}{\partial\varphi} +O(\epsilon).
\]
After applying the averaging transformation, we find that the
system~(\ref{AA1Averaging}) takes the form
\begin{eqnarray}\label{AA2Averaging}
\dot{\widehat{{\cal D}}} & = & -\epsilon^{1/2} \bar{F}_{11}
-\epsilon \widehat{{\cal D}} \Big( F_{12}
+\frac{3}{L_*^4}\frac{\partial\Lambda}{\partial\varphi}\Big)
+O(\epsilon^{3/2}), \nonumber \\
\dot{\widehat{G}} & = & -\epsilon F_{22}+O(\epsilon^{3/2}), \nonumber \\
\dot{\widehat{\varphi}} & = & -\epsilon^{1/2}\frac{3}{L_*^4}\widehat{{\cal D}}
+ \epsilon \Big(\frac{6}{L_*^5}{\widehat{{\cal D}}}^2+F_{32}
+\frac{3}{L_*^4}\Lambda\Big)+O(\epsilon^{3/2}), \nonumber \\
\dot{\widehat{g}} & = &\quad \epsilon F_{42}+O(\epsilon^{3/2}).
\end{eqnarray}
Finally, we drop the $O(\epsilon^{3/2})$ terms in~(\ref{AA2Averaging}) and
average the remaining truncated system
to obtain the {\em second order partially averaged system}
\begin{eqnarray}\label{2ndoave}
\dot{\widetilde{{\cal D}}}
& = & -\epsilon^{1/2} \bar{F}_{11}
-\epsilon\widetilde{{\cal D}}\bar{F}_{12},\nonumber \\
\dot{\widetilde{G}} & = & -\epsilon \bar{F}_{22}, \nonumber \\
\dot{\widetilde{\varphi}} & = & -\epsilon^{1/2}\frac{3}{L_*^4}\widetilde{{\cal D}}
+ \epsilon \big(\frac{6}{L_*^5}{\widetilde{{\cal D}}}^2+\bar{F}_{32}\big),\nonumber \\
\dot{\widetilde{g}} & = & \quad \epsilon \bar{F}_{42}.
\end{eqnarray}
This system is the averaged form of system~(\ref{AA1Averaging})
after dropping its $O(\epsilon^{3/2})$ terms;
however, this coincidence is fortuitous in this case.
In general, one has to employ an averaging
transformation in order to obtain the second order averaged system.
To explain the evolutionary dynamics at resonance, we will
replace the actual dynamical
equations by the second order partially averaged system. As explained in
Appendix~\ref{appendixb}, this is a reasonable approximation over a limited time-scale.
We remark that although the second order partially averaged system will be used
to explain some of the features of our model system near its resonances,
the actual dynamics predicted by our model is certainly much more complex than the
averaged equations reveal. In particular, we expect that near the resonances---perhaps
in other regions of the phase space as well---there are chaotic invariant sets and,
therefore, transient chaotic motions~\cite{cmr2}.
The averaged system is nonlinear and could in general exhibit chaos; however,
we have not encountered chaos in the second order averaged equations obtained
from the model under consideration here.
It remains to compute $\bar{F}_{ij}$, where $F_{ij}$ are defined
in~(\ref{Fij}). For this, we recall equation~(\ref{Hext})
and set $\ell=n\Omega t/m+\widetilde{\varphi}$,
expand the trigonometric terms in the Fourier series using trigonometric
identities, and then average over the variable $t$. The required
averages involving the external perturbation can be computed from the average of
${\cal H}_{\mbox{\scriptsize ext}}$. In fact, after some computation, we find that
\begin{eqnarray*}
\bar{{\cal H}}_{\mbox{\scriptsize ext}}&:=&\frac{\Omega}{2\pi m}\int_0^{2\pi m/\Omega}
{\cal H}_{\mbox{\scriptsize ext}}(L_*,\widetilde{G},
\frac{n\Omega}{m}t+\widetilde{\varphi},\widetilde{g},t)\,dt\\
&=&
T_c(L_*,\widetilde{G},\widetilde{g})\cos{m\widetilde{\varphi}}
+T_s(L_*,\widetilde{G},\widetilde{g})\sin{m\widetilde{\varphi}},
\end{eqnarray*}
where, for $n=1$,
\begin{eqnarray}\label{tcts}
T_c & := & \frac{L^4\Omega^2}{4}\big[
\alpha A_m(e)\cos{2g}+\beta A_m(e)\sin{2g}\cos{\rho}
-\beta B_m(e)\cos{2g}\sin{\rho}\big],\nonumber \\
T_s & := & \frac{L^4\Omega^2}{4}\big[-\alpha B_m(e)\sin{2g}+\beta
A_m(e)\sin{2g}\sin{\rho}
+\beta B_m(e)\cos{2g}\cos{\rho}\big],
\end{eqnarray}
and for $n\ne 1$, $\bar{{\cal H}}_{\mbox{\scriptsize ext}}=0$ so
that in this case we can define $T_c=T_s=0$.
The averages of $f_D$, $D\in\{L,G,\ell,g\}$, are given in~(\ref{fqAve3}).
The terms involving radiation damping
that appear in the partially averaged system are obtained from
these expressions by Taylor expansion about the resonant orbit using
equation~(\ref{coordtrans}). For example, we will use
\[
\Gamma(G):=\left.\bar{f}_L\right|_{L=L_*}=
-\left.\frac{1}{G^7}\left(8 + \frac{73}{3} e^2 + \frac{37}{12} e^4 \right)
\right|_{e^2=1-G^2/L_*^2}
\]
and the average of $\partial f_L/\partial L$ at resonance given by
\[
\left.\frac{\partial}{\partial L} \bar{f}_L \right|_{L=L_*}=
\left.-\frac{1}{3 L_*^3 G^5}(146+37 e^2)\right|_{e^2=1-G^2/L_*^2}.
\]
The second order partially averaged system~(\ref{2ndoave}) is
thus given explicitly by
\begin{eqnarray}\label{ex2ndoave}
\dot {\cal D} &=& -\epsilon^{1/2}\Big[-mT_c\sin{m\varphi}+mT_s\cos{m\varphi}
+\frac{\Delta}{G^7}\Big(8 + \frac{73}{3} e^2 + \frac{37}{12} e^4 \Big)\Big]
\nonumber \\
&&\quad-\epsilon {\cal D} \Big[-m\frac{\partial T_c}{\partial L}\sin{m\varphi}
+m\frac{\partial T_s}{\partial L}\cos{m\varphi}
+\frac{\Delta}{3L_*^3 G^5}\Big(146 + 37 e^2 \Big)\Big], \nonumber \\
\dot G &=& -\epsilon\Big[\frac{\partial T_c}{\partial g}\cos{m\varphi}
+\frac{\partial T_s}{\partial g}\sin{m\varphi}
+\frac{\Delta}{L_*^3 G^4}(8 + 7 e^2 )\Big], \nonumber \\
\dot \varphi &=& -\epsilon^{1/2}\frac{3}{L_*^4} {\cal D}+\epsilon\Big(
\frac{6}{L_*^5}{\cal D}^2+\frac{\partial T_c}{\partial L}\cos{m\varphi}
+\frac{\partial T_s}{\partial L}\sin{m\varphi}\Big), \nonumber \\
\dot g &=& \:\epsilon\Big(\frac{\partial T_c}{\partial G}\cos{m\varphi}
+\frac{\partial T_s}{\partial G}\sin{m\varphi}\Big),
\end{eqnarray}
where we have dropped the tildes.
It is clear that $L$ in the
expressions involving $T_c$, $T_s$, and $e$ must be replaced by $L_*$, the value of
$L$ at resonance.
Having derived the equations for the averaged dynamics~(\ref{ex2ndoave}),
we now turn our attention to the consequences of these equations and the comparison
of predictions based on them with the actual dynamics given by~(\ref{MathEQM}).
\section{First Order Averaged Dynamics}\label{drt}
The first order partially averaged system, obtained from~(\ref{ex2ndoave})
by dropping the
$O(\epsilon)$ terms, is given by
\begin{eqnarray}\label{foav}
\dot{\cal D}&=&-\epsilon^{1/2}\Big[-mT_c\sin{m\varphi}+mT_s\cos{m\varphi}
-\Delta \Gamma(G)\Big],\nonumber\\
\dot G&=&0,\nonumber\\
\dot \varphi&=&-\epsilon^{1/2}\Big(\frac{3}{L_*^4}\Big) {\cal D}, \nonumber\\
\dot g&=&0.
\end{eqnarray}
In this approximation, the variables $G$ and $g$ are constants fixed
by the initial conditions while the remaining system in ${\cal D}$ and $\varphi$
is equivalent to a pendulum-type equation with torque; namely,
\begin{equation}\label{pendulum}
\ddot\varphi+\frac{3\epsilon}{{L_*}^4}(mT_c\sin{m\varphi}-mT_s\cos{m\varphi})
=-\frac{3\epsilon}{{L_ *}^4}\: \Delta\: \Gamma(G).
\end{equation}
We also have a second order differential equation---a harmonic
oscillator with slowly varying frequency---for the deviation
${\cal D}$ given by
\begin{equation}\label{fodeveq}
\ddot {\cal D}+\epsilon w^2 {\cal D}=0,
\end{equation}
where
\begin{equation}\label{w}
w:=\Big[\frac{3m^2}{L_*^4}(T_c\cos{m\varphi}+T_s\sin{m\varphi})\Big]^{1/2}.
\end{equation}
These results can be formally justified if they are recast in terms of a
new temporal variable given be $\epsilon^{1/2}t$; however, we use $t$ in
order to facilitate comparison with the actual dynamics.
To show that~(\ref{fodeveq}) is an oscillator, we must show that $w$ is a
real number.
To this end, we suppose that during capture into the resonance the orbit
is near an elliptic region of the first order partially averaged system~(\ref{foav}).
This system is Hamiltonian with an effective energy of the form
\[
-\epsilon^{1/2}\Big[\frac{1}{2}\Big(\frac{3}{L_*^4}\Big) {\cal D}^2+U(\varphi)\Big],
\]
where $U$ represents the effective potential energy given by
\[
U:=-(T_c\cos{m\varphi}+T_s \sin{m\varphi})+(\Delta\Gamma)\varphi.
\]
If $({\cal D}, \varphi)=(0,\varphi_0)$ is an elliptic rest point of the
first order partially averaged system; then,
$U'(\varphi_0)=0$ and $U''(\varphi_0)>0$, where a prime denotes
differentiation with respect to $\varphi$. It follows that
\[U''(\varphi_0)= m^2(T_c\cos{m\varphi_0}+T_s\sin{m\varphi_0}),\]
so that $w_0$, $w_0:=[3 U''(\varphi_0)]^{1/2}/L_*^2$, must be real.
To show that the frequency $\xi:=\epsilon^{1/2}w$ is slowly varying,
we differentiate this frequency with respect to time to
obtain
\[
\dot\xi=-\epsilon\frac{9m^3}{2L_*^8w}(T_s\cos{m\varphi}-T_c\sin{m\varphi}){\cal D}+
O(\epsilon^{3/2}).
\]
It follows from the first order averaged equations that
$\varphi-\varphi_0$ is expected to be of order unity; thus,
$\xi$ is nearly constant over a time-scale of order
$\epsilon^{-1/2}$ since $\dot\xi$ is proportional to $\epsilon$.
In particular, ${\cal D}$ varies on the time-scale $\epsilon^{1/2}t$.
Inspection of equations~(\ref{pendulum}) and~(\ref{fodeveq})
reveals that while ${\cal D}$
is predominantly a simple harmonic oscillator
with frequency $\xi_0=\epsilon^{1/2}w_0$,
the motion of $\varphi$ in time could be
rather complicated involving essentially all harmonics
of $\xi_0$. Therefore, $\varphi$ cannot in general
be considered a slow variable in time.
Librational motions (periodic motions in the phase plane) of the pendulum-type
averaged equation~(\ref{pendulum}) correspond to orbits of the original system
that are captured into the resonance. On the other hand, if there are no
librational motions, then all orbits pass through the resonance.
Thus, we observe a necessary condition for
capture: the pendulum system must have rest points in its phase plane.
The rest points of the first order averaged system are given by
${\cal D}=0$ and $\varphi=\varphi_0$, with
\begin{equation}\label{resteq}
R\sin(m\varphi_0+\eta)=-\Delta \Gamma,
\end{equation}
where $mT_c=R\cos{\eta}$ and $mT_s=-R\sin{\eta}$.
As an immediate consequence, we have the following proposition:
{\em For the $(m:n)$ resonance, if $n\ne 1$, there is no capture into resonance}.
On the other hand, for $n=1$,
if there are librational motions, then
we must have $\Delta|\Gamma/R| \le 1$.
These observations suggest that in the presence of sufficiently strong
radiation damping, i.e. for $\Delta$ sufficiently large, there are no
librational
motions and, as a result, each orbit will pass through all the resonant
manifolds that it encounters.
In particular, since $L=a^{1/2}$, the semimajor
axis of the corresponding osculating ellipse will collapse to zero.
Moreover, capture is only possible for $(m:n)$ resonances with $n=1$ when
$\Delta|\Gamma/R| \le 1$. In this case, if $\Delta|\Gamma/R| < 1$, then
the rest points in the phase plane of the pendulum-type equation~(\ref{pendulum})
at the resonance are all nondegenerate. This is precisely Neishtadt's
condition $B$.
The quantities $\Gamma$ and $R$
depend (nonlinearly) on the variables $(L, G, g)$ as well as
the parameters of the system; therefore, the precise range of their
ratio $\Gamma/R$ is difficult to
specify in general. However, the value of this ratio
can be determined numerically. In fact, to find orbits that are captured
into resonance as displayed in Figures~1--3, we use the first order averaged
system. A rest point of the first order system corresponds
to a ``fixed'' resonant orbit with ${\cal D}=0$ and constant
$G$, $g$, and $\varphi=\varphi_0$. The main characteristics of these orbits
do not change with respect to the slow variable $\epsilon^{1/2}t$;
that is, the resonant orbits are essentially fixed only over a time-scale of
order $\epsilon^{-1/2}$.
After choosing $G$ and $g$,
we solve for $\varphi_0$ at a rest point,
and then test to see if the resulting resonant orbit corresponds
to a libration point of the first order averaged system (e.g. $w_0^2$
must be positive). If it does, we
convert the point with coordinates
$(L_*, G,\varphi_0,g)$ to polar coordinates and then numerically integrate
our model~(\ref{MathEQM}) backward in time from this starting value.
When the backward integration in time is carried
out over a sufficiently long time interval, the orbit is expected to
leave the vicinity of the resonance.
After this occurs, we integrate forward
in time to obtain an orbit of~(\ref{MathEQM}) that is captured
into the resonance.
The first order partially averaged system~(\ref{foav}) is
Hamiltonian; therefore,
we do not expect that its dynamics would give a complete picture of the average
dynamics near a resonance for our {\em dissipative} system. In particular, the
librational motions of the pendulum-like system are not structurally stable.
Indeed, the phase space for the full four dimensional system~(\ref{foav})
is foliated by invariant two dimensional subspaces parametrized by the
variables $G$ and $g$.
If a librational motion exists in a leaf of the foliation, then
the corresponding rest points, the elliptic rest point and the
hyperbolic rest point (or points) at the boundary of the librational region for the
associated pendulum-like system, are degenerate in the four
dimensional phase space of the first order averaged system;
their linearizations have
two zero eigenvalues corresponding to the directions normal to the leaf.
In the second order partially averaged system, viewed as a small perturbation
of~(\ref{foav}),
we expect that only a finite number of the degenerate rest points survive. These
correspond to the continuable periodic solutions described in~\cite{cmr2}.
We expect the corresponding perturbed rest points
to be hyperbolic. As a result, they persist even
when higher order effects are added. The dimensions and positions of the
stable and unstable manifolds of these rest points determine the average
motion near the resonance.
In particular, if one of these rest points is a hyperbolic sink,
then there is an open set of initial conditions that correspond to orbits
{\em permanently} captured into resonance.
Our numerical experiments have provided no evidence for this behavior.
However, another possibility---that is consistent with our
numerical experiments---is that a rest point of the perturbed system
in phase space is stable
along two directions but unstable in the remaining two directions.
Thus a trajectory with initial point near the corresponding stable manifold will
be captured into resonance and undergo librational motions near the resonant
manifold until it spirals outward along the unstable manifold.
\section{Dynamical Evolution During Resonance}\label{dedr}
The dynamics of system~(\ref{ex2ndoave}) is expected to
give a close approximation of the
near resonance behavior of our original model~(\ref{MathEQM}) over
a time-scale of order $\epsilon^{-1/2}$. In particular,
a basic open problem is the following: Determine the positions of the rest points
and the corresponding stable and unstable manifolds of~(\ref{ex2ndoave}).
A solution of this problem would provide the information needed to analyze
the most important features of the behavior of the orbits that pass through the
resonance as well as those that are captured by the resonance. Unfortunately,
rigorous analysis of the dynamics in
the four dimensional phase space of system~(\ref{ex2ndoave})
seems at present to be very difficult. Thus, in lieu of a complete and rigorous
analysis of the dynamics, we will show how to obtain useful information from
approximations of the second order partially averaged system.
A typical example of the behavior of the orbit as it passes through a
$(1:1)$ resonance is depicted in Fig.~\ref{fig1}. The semimajor axis undergoes
librations of increasing amplitude around the resonant value. Furthermore,
the eccentricity of the orbit generally decreases while the orbit is trapped.
We propose to analyze the oscillations of ${\cal D}$ by using the
$O(\epsilon)$ approximation of the second order differential equation for ${\cal D}$
derived from the second order partially averaged equations~(\ref{ex2ndoave}).
In fact, a simple computation
yields the following differential equation
\begin{equation}\label{devos}
\ddot{\cal D}+\epsilon \gamma \dot{\cal D}+\epsilon w^2 {\cal D}=0,
\end{equation}
where
\begin{equation}\label{gam25}
\gamma:=m\Big(
\frac{\partial T_s}{\partial L}\cos{m\varphi}
-\frac{\partial T_c}{\partial L}\sin{m\varphi}\Big)
+\frac{\Delta}{3L^3G^5}(146+37e^2)
\end{equation}
is evaluated at $L=L_*$ in equation~(\ref{devos})
and $w$ is given by equation~(\ref{w}).
Similarly, the second order partially averaged system~(\ref{ex2ndoave})
can be used to derive an equation for the temporal evolution of $\varphi$
to order $\epsilon$; however, the resulting equation turns out to
be identical with equation~(\ref{pendulum}).
It is clear from inspection of the differential equation~(\ref{devos}) that
$L=L_*+\epsilon^{1/2}{\cal D}$ oscillates about its resonant
value $L_*$ with a libration frequency given approximately by
$\xi_0=\epsilon^{1/2} w_0$.
The magnitude of this frequency is in agreement with our numerical results to
within a few percent for the case of resonance trapping depicted in Fig.~\ref{fig1}.
Moreover,
the amplitude of the oscillations will increase (decrease) if $\gamma<0$
($\gamma>0$).
The sign of $\gamma$ varies with the choice of parameters, variables, and
the order of the resonance. Our numerical experiments---which have
by no means been exhaustive---indicate that the $(1:1)$ resonance for the
linearly polarized incident wave with $\alpha=1$ and $\beta=0$
is special in that $\gamma$ is negative at the
elliptic rest points of the first order averaged system.
However, for a resonance with order $m>1$, $\gamma$ is
not of fixed sign. Nevertheless, our numerical experiments indicate that
the amplitude of librations generally increases over a long time-scale
(cf. Fig.~\ref{fig2}).
Figures~\ref{fig2} and~\ref{fig3} illustrate the apparent richness of the
evolutionary dynamics near resonances with $m>1$.
It should be emphasized that the damped (or anti-damped)
oscillator~(\ref{devos}) approximates the
actual dynamics only over a time-scale of order $\epsilon^{-1/2}$
(cf. Appendix~\ref{appendixb}).
To obtain an approximate equation for the envelope of the oscillations,
we consider a time-dependent change of variables for the oscillator~(\ref{devos})
defined by the relation
\begin{equation}\label{Deq}
{\cal D}=Ve^{-(\epsilon/2)\int_0^t\gamma(s)\,ds}.
\end{equation}
After a routine calculation, we obtain the differential equation
\begin{equation}\label{harmos}
\ddot V+\epsilon w^2 V=O(\epsilon^{3/2}).
\end{equation}
Thus, to order $\epsilon$, $V(t)$ in equation~(\ref{Deq}) is
the solution of a harmonic oscillator equation with a
slowly varying frequency $\xi=\epsilon^{1/2}w$. In fact,
equation~(\ref{harmos}) is identical to the oscillator~(\ref{fodeveq})
with frequency $\xi$ obtained from the first order
averaged system.
The solution ${\cal D}(t)$ of the second order equation~(\ref{devos})
in this approximation
is obtained by modulating the amplitude of $V$.
We note that if $\gamma$ and $w$ are constants, then formula~(\ref{Deq})
reduces to
\[
{\cal D}(t)={\cal A}e^{-\epsilon\gamma t/2} \cos(\epsilon^{1/2}wt+\tau),
\]
where ${\cal A}$ and $\tau$ are constants depending on the initial conditions;
this is the standard result for the damped
(or anti-damped) harmonic oscillator with constant coefficients.
Equation~(\ref{Deq}) certainly agrees qualitatively with the numerical
experiments reported in Figures~\ref{fig1}--\ref{fig3}.
To check the accuracy quantitatively, we numerically integrated
system~(\ref{MathEQM}) and simultaneously evaluated
the integral of $\gamma$ in order
to obtain an approximation for the envelope
of ${\cal D} (t)$.
Over an interval of time of length $691$, which is of order $\epsilon^{-1/2}$,
the value of ${\cal D}$ at its maximum predicted from
equation~(\ref{Deq}) differs from the value obtained by numerical integration
of the differential equations by only a few percent (actually $2.5\%$).
The error in the comparison of the two values for the envelope of ${\cal D}$
grows slowly over time; in fact, it remains less than $20\%$ over a time interval
of length $9000$.
The evolutionary dynamics near resonance
is characterized by librational motions described
above as well as variations in angular momentum.
Our numerical experiments suggest considerable variation in the
angular momentum while the orbit is trapped in resonance.
The remainder of this section is devoted to a theoretical explanation
of this phenomenon.
During the time that an orbit is captured into resonance,
we expect that on average it will reside near an elliptic rest point of the
first order partially averaged system.
Thus, we expect the dynamical variables to satisfy approximately
equation~(\ref{resteq}) while the orbit is trapped.
We will show how to determine the dynamical behavior of the angular
momentum $G$ under this assumption.
In particular, for the $(1:1)$ resonance, we claim that $\dot G>0$; that is,
the orbital angular momentum increases.
It is easy to show that
\begin{eqnarray*}
\frac{\partial T_c}{\partial g}-2 T_s&=&
-\frac{1}{2L_*^2}[\alpha \sin 2g-\beta\cos(2g + \rho)]S_m(e), \\
\frac{\partial T_s}{\partial g}+2 T_c&=&
\frac{1}{2L_*^2}[\alpha \cos 2g+\beta\sin(2g + \rho)]S_m(e),
\end{eqnarray*}
where
\[ S_m(e):=m^2[A_m(e)-B_m(e)].\]
After substituting these identities as well as the condition~(\ref{resteq})
into the expression for
$\dot G$ in~(\ref{ex2ndoave}),
we find that
\begin{equation}\label{gdot}
\dot G= -\epsilon [P_m(e)+L_*^{-2} S_m(e) Q_m],
\end{equation}
where, after some manipulation using the relation $e=(1-G^2/L^2_*)^{1/2}$,
\begin{eqnarray*}
P_m(e)&:=&\frac{\Delta}{G^7}\Big[(1-e^2)^{3/2}(8+7 e^2)
-\frac{2}{m}\big( 8+\frac{73}{3} e^2+ \frac{37}{12} e^4 \big)\Big],\\
Q_m&:=&\frac{1}{2} [\alpha\sin(m\varphi-2g)+\beta\cos(m\varphi-2g-\rho)].
\end{eqnarray*}
The function $S_m(e)\to 0$ as $e\to 0$; therefore,
it is clear that $\dot G>0$ near $e=0$ as long as $P_m(e)<0$.
We have $P_1(e)<-8\Delta/G^7$ for $02$, we have $P_m(e)>0$ near $e=0$.
Our conclusion is that near the $(1:1)$ resonance, we can expect
$\dot G>0$ during capture provided the osculating
ellipse is not too eccentric. The other cases are more delicately balanced.
However, for a specific choice of parameters and for a specific
resonance, one can use equation~(\ref{gdot}) to predict
the behavior of $G$ near the resonance.
Consider, for instance, the system~(\ref{MathEQM}) with the following parameter values:
$\epsilon=10^{-4}$, $\Delta=\delta/\epsilon=10^{-3}$,
$\alpha=1$, $\beta=0$, $\rho=0$, and $\Omega= 1$.
The orbit, as depicted in Figure~\ref{fig1}, with
initial conditions given by
\[
(p_r,p_\theta,r,\theta)=(0.2817, 0.6389, 1.6628, 2.9272)
\]
appears to be trapped in $(1:1)$ resonance with a sojourn time of approximately
$40000$.
Moreover,
during this period of time
the angular momentum appears to increase from approximately
$0.66$ to $0.95$.
To demonstrate that our theoretical scheme using the second order partially
averaged system agrees with our numerical results, let us determine the
time-scale over which the angular momentum changes from $0.78$ to $0.95$.
According to the numerical integration of system~(\ref{MathEQM}) this occurs
over about $32700$ time units, whereas our approximate
analysis of the second order partially
averaged system predicts a time-scale of
\[
-\frac{1}{\epsilon}\int_{0.78}^{0.95}
\big[P_1\big((1-G^2/L^2_*)^{1/2}\big)\big]^{-1}\,dG\approx 34200,
\]
a time that is within less than $5\%$ of the value given by our numerical experiment.
In this integral, we have used the fact that in equation~(\ref{gdot}),
the term $L_*^{-2}S_1(e)Q_1$ is small compared with $P_1(e)$.
It should be clear from the results of this section that---by employing
the second order partially averaged system---we have been able
to provide theoretical explanations for the behaviors of semimajor axis
$a$ and eccentricity $e$ of the orbit that is trapped in resonance.
The basic dynamical equations for our model are only valid to order $\epsilon$,
since physical effects of higher order have been neglected in equation~(\ref{MathEQM});
therefore,
it would be inappropriate to go beyond the second order averaged system in
this case. The averaging method is general, however, and can be carried
through to any order.
\section{Conclusion}\label{c}
Some of the observed structures in the solar system are the direct results of
evolutionary dynamics while trapped in resonance
(see \cite{melita}, \cite{beauge}, and \cite{gomes}).
We have argued in this paper that such dynamics can in general be explained using
the method of averaging. In particular, the main aspects of the
evolution of the dynamical system during the time that it is locked in resonance
can be understood on the basis of the second order partially averaged dynamics.
We have illustrated these ideas using a particular model involving a Keplerian
binary system that emits and absorbs gravitational radiation.
\appendix
\section{Delaunay equations of motion}\label{appendixa}
The Delaunay equations of motion for our system are given in~\cite{cmr2} as
\begin{eqnarray}\label{D2EQM}
\frac{dL}{dt} & = & -\epsilon \left(\frac{\partial{\cal C}}{\partial \ell}
\phi(t) +
\frac{\partial{\cal S}}{\partial \ell} \psi(t)\right) + \delta f_L,
\nonumber \\
\frac{dG}{dt} & = & -\epsilon \left(\frac{\partial{\cal C}}{\partial g}
\phi(t) +
\frac{\partial{\cal S}}{\partial g} \psi(t)\right) + \delta f_G,
\nonumber \\
\frac{d\ell}{dt} & = & \omega + \epsilon \left(\frac{\partial{\cal C}}{\partial
L} \phi(t) +
\frac{\partial{\cal S}}{\partial L} \psi(t)\right)
+ \delta f_{\ell}, \nonumber \\
\frac{dg}{dt} & = & \epsilon \left(\frac{\partial{\cal C}}{\partial G}
\phi(t) +
\frac{\partial{\cal S}}{\partial G}\psi(t)\right) + \delta f_g,
\end{eqnarray}
where
\[
\phi(t):=\frac{1}{2} \alpha\Omega^2\cos\Omega t,\qquad
\psi(t):=\frac{1}{2} \beta\Omega^2\cos(\Omega t+\rho),
\]
${\cal C}$ and ${\cal S}$ are given by equation~(\ref{CSfs}) and
\begin{eqnarray}\label{fequations}
f_L & = &
\frac{4}{Lr^3}\left[1-\frac{16}{3}\frac{L^2}{r} +
\left(\frac{20}{3}L^2- \frac{17}{2}G^2\right)\frac{L^2}{r^2}
+\frac{50}{3}\frac{L^4G^2}{r^3}-\frac{25}{2}\frac{L^4G^4}{r^4}\right],
\nonumber \\
f_G & = & -\frac{18 G}{L^2r^3}\left(1 -\frac{20}{9}\frac{L^2}{r}
+\frac{5}{3}\frac{L^2G^2}{r^2}\right), \nonumber \\
f_{\ell} & = & \frac{2\sin{\hat{v}}}{eL^3Gr^2}\left[4e^2+
\frac{1}{3}
\left(73G^2-40L^2\right)\frac{1}{r}
-2\left(1+\frac{70}{3}L^2-\frac{29}{2}G^2\right)\frac{G^2}{r^2}
\right. \nonumber \\
&& \hspace*{.5in}\left.-\frac{25}{3}\frac{L^2G^4}{r^3}
+25\frac{L^2G^6}{r^4} \right], \nonumber \\
f_g & = &
-\frac{2\sin{\hat{v}}}{eL^2r^3}\left[11+
\left(7G^2-\frac{80}{3}L^2\right)\frac{1}{r}-
\frac{25}{3}\frac{L^2G^2}{r^2}+25\frac{L^2G^4}{r^3}\right].
\end{eqnarray}
The $f_D$ can be expressed, in principle, in terms
of Delaunay elements. In fact, expressions of the form $r^{-q}$ and
$r^{-q}\sin{\widehat{v}}$ can be expanded as Fourier series in terms of the mean
anomaly $\ell$ (see \cite{Kov}, p. 54). For example, the following results are known
\cite{Kov,watson}
\begin{eqnarray}\label{watexpan}
\frac{L^2}{r} & = & 1+2 \sum_{n=1}^{\infty}J_n(ne) \cos{n\ell}, \\
\frac{L^4}{r^2}\sin{\hat{v}} & = &
\frac{(1-e^2)^{1/2}}{e} \sum_{n=1}^{\infty}2n J_n(ne) \sin{n\ell}.
\end{eqnarray}
For the calculation of the other similar
expressions---and eventually the $f_D$---see the methods of
classical celestial mechanics involving the Bessel
functions described in~\cite{Kov} and~\cite{watson}.
Only the {\em averages} of the $f_D$---calculated in~\cite{cmr2} and
given in equation~(\ref{fqAve3})---are
needed for the purposes of this paper.
\section{the method of averaging}\label{appendixb}
In this section, we give a brief introduction to the method of averaging as it is
used in this paper. More detailed accounts can be found in many excellent sources.
We suggest references \cite{arnold2},\cite{murdock}, \cite{chirikov},
\cite{wig}, and~\cite{sv},
as they offer a range of mathematical sophistication
as well as diverse points of view.
The method of averaging is usually applied to differential equations in one
of several standard forms.
In celestial mechanics, perhaps the most common standard form
is ``angular standard form''
\[ \dot I=\epsilon f(I,\theta),\qquad \dot \theta=\omega(I)+\epsilon g(I,\theta),\]
where $I$ is an $m$-dimensional variable, $\theta$ is an $n$-dimensional
angular variable, both $f$ and $g$ are $2\pi$ periodic in $\theta$, and
$\epsilon$ is a small parameter.
This form is often obtained after a model of disturbed two-body motion is expressed
in the action-angle variables appropriate for the
unperturbed Kepler problem as in equation~(\ref{AAAveraging}).
The second form is ``time periodic standard form''
\begin{equation}\label{tpsf}
\dot x=\epsilon F(x,t),
\end{equation}
where $x$ is an $m$-dimensional variable and $F$ is $T$-periodic in time.
This form is obtained for many different physical problems.
In particular, in our case --- i.e. equation~(\ref{AA1Averaging})---
the small parameter is $\epsilon^{1/2}$ and the equation of motion is a Taylor
series in $\epsilon^{1/2}$.
Of course, the time periodic standard form can be viewed as a special
case of the angular standard form by introducing a new one dimensional angular variable:
\[
\dot x=\epsilon F^*(x,\theta),\qquad \dot\theta=\frac{2\pi}{T}.
\]
However, it is useful to separate the two cases.
The Averaging Principle is often applied to differential equations after
they are transformed to one of the above standard forms. It consists of
approximating a solution of the angular standard form by the corresponding solution
(with the same initial conditions) of the averaged system
\[
\dot J=\epsilon\frac{1}{(2\pi)^n}\int_0^{2\pi}f(J,\theta)\, d\theta
\]
or, for the time periodic standard form,
by the equation
\[
\dot y=\epsilon \frac{1}{T}\int_0^{T}F(y,t)\, dt.
\]
The validity of the approximation produced by the Averaging Principle is often left
unverified. However, as many examples show, the Averaging Principle is not
valid in general. Thus, a careful analysis of a physical problem must
always deal with this issue. Fortunately, the problem has been treated in great
detail in the mathematical literature. Many useful and rigorous
results are now available.
The basic result involves only one angular variable; that is, it pertains
either to
the time periodic standard form, or to the case when $\theta$ is one dimensional
in the angular standard form. The Averaging Theorem asserts in this
case (with appropriate smoothness of the functions involved)
that the averaged solution
approximates the original solution with an error of order $\epsilon$ on a
time-scale of order $1/\epsilon$;
for example, if $x(0)=y(0)$, then there are constants
$C_1$ and $C_2$ that do not depend on $\epsilon$ such that
\[ |x(t)-y(t)|0$ and $c_2>0$
such that
\[
|x(t,u_0+\epsilon K(u_0,0),\epsilon)-
[u(t,u_0,\epsilon)+\epsilon K(u(t,u_0,\epsilon),t)]|