%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%%%%%%%%%%%%%%%%%%%%%%% \newcommand{\bbC}{\Bbb C } \newcommand{\bbN}{\Bbb N } \newcommand{\bbR}{\Bbb R } \newcommand{\bbZ}{\Bbb Z } \newenvironment{pf}{{\it Proof}.}{$\Box$} \begin{document} % \draft command makes pacs numbers print %\draft \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.} } \date{\today} \title{Gravitational Ionization: A Chaotic Net in the Kepler System} \maketitle \begin{abstract} The long term nonlinear dynamics of a Keplerian binary system under the combined influences of gravitational radiation damping and external tidal perturbations is analyzed. Gravitational radiation reaction leads the binary system towards eventual collapse, while the external periodic perturbations could lead to the ionization of the system via Arnold diffusion. When these two opposing tendencies nearly balance each other, interesting chaotic behavior occurs that is briefly studied in this paper. It is possible to show that periodic orbits can exist in this system for sufficiently small damping. Moreover, we employ the method of averaging to investigate the phenomenon of capture into resonance. \end{abstract} %\pacs{} \section{Introduction} The Kepler problem traditionally involves two point masses $m_1$ and $m_2$ moving under the influence of their mutual gravitational interaction. The incorporation of relativistic gravitational effects in the Kepler system brings about three main post-Newtonian modifications in the traditional picture: (i) There exist gravitoelectromagnetic post-Newtonian effects that need to be taken into account. The most important of these is the periastron precession that has played an important role in the historical development of general relativity. (ii) Moreover, the system emits gravitational radiation and hence there exist {\it radiative} post-Newtonian effects in the orbit that are characteristic of gravitational radiation damping. This radiative damping is consistent with the observed rate of inward spiraling in the Hulse-Taylor binary pulsar system \cite{hulse,taylor}. (iii) The system is expected to be affected by gravitational waves that have been emitted by other systems as well as by their post-Newtonian tidal perturbations. Gravitational radiation has not yet been directly observed; however, the Hulse-Taylor binary pulsar data appear to be consistent with the notion that gravitational radiation is emitted by binary systems \cite{taylor}. Half of all stars are expected to be members of binary or multiple systems that could emit gravitational waves; therefore, there might be a cosmic background of gravitational radiation that has been generated by various sources throughout the history of the universe. Moreover, there could also be primordial gravitational waves left over from the epoch at which the Hubble expansion began. The inclusion of all of these effects in the Kepler problem is impractical; in fact, the two-body problem in general relativity is intractable. To render the problem amenable to mathematical analysis, it is therefore necessary to replace the actual problem by a model that contains the main physical effects of interest. We consider a Kepler system in which the main effects of emission and absorption of gravitational radiation are taken into consideration. Gravitation is a spin-2 field; therefore, the main radiative effects first occur in general at the quadrupole level in emission as well as in absorption. At this level, the wavelength of the radiation is large compared to the size of the system. In this work, we analyze the radiative perturbations of a Keplerian binary system in the quadrupole approximation. In our previous papers \cite{cmr,cmr1}, we considered the problem of ionization of a Keplerian binary system by a normally incident periodic gravitational wave. To render the absorption problem tractable, we simply ignored the emission of gravitational radiation by the binary system. Therefore, the influence of the gravitational radiation reaction on the binary orbit was not taken into account. In the absence of this dissipative effect, our model of gravitational ionization turned out to be a Hamiltonian system to which the basic results of the Kolmogorov-Arnold-Moser~(KAM) theory could be applied under certain circumstances. The main purpose of the present work is to take gravitational radiation damping into account. Thus, we study in this paper the long term nonlinear evolution of a Keplerian binary system under the combined action of perturbations due to the emission and absorption of gravitational radiation. The nonlinear dynamics of a Keplerian binary that emits and absorbs gravitational radiation is given in our model by \begin{equation}\label{BasicEq} \frac{d^2r^i}{dt^2} + \frac{kr^i}{r^3} + {\cal R}^i = -\epsilon {\cal K}_{ij}(t)\: r^j, \end{equation} where ${\bf r}:= {\bf x}_1-{\bf x}_2$ denotes the relative two-body orbit, $r$ is the length of $\bf r$, and $k = G_0(m_1+m_2)$. Here the radiation reaction term ${\cal R}$ is given by \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 r}\right], \end{equation} where ${\bf v}$ is the relative velocity, $G_0$ is Newton's constant, $c$ is the speed of light in vacuum, and an overdot represents differentiation with respect to time, i.e. $\dot{r} = dr/dt$. We emphasize here that equation (\ref{BasicEq}) represents the simplest equation for the relative motion of two point masses $m_1$ and $m_2$ that contains the main dynamical effects of emission and absorption of gravitational radiation that we wish to study here; in fact, other post-Newtonian contributions have been simply ignored. The justification for this approach emerges from our analysis; that is, the substitution of the actual intractable problem by our simple model permits us to arrive at interesting results of possible astrophysical interest. Let the background spacetime metric be given by $g_{\mu\nu} = \eta_{\mu\nu} + \epsilon \chi_{\mu\nu}$, where $\eta_{\mu\nu}$ denotes the Minkowski metric, $\epsilon$ is a small parameter indicative of the strength of the external radiation field $(0 < \epsilon << 1)$, and $\chi_{\mu\nu}$ represents the gravitational radiation field pervading the space. We impose the transverse-traceless gauge condition such that $\chi_{0\mu}=0$, ${\partial}_j \chi_{ij} = 0$, and $(\chi_{ij})$ is traceless; moreover, $\chi_{ij}$ satisfies the wave equation $\Box \chi_{ij} = 0$. Imagine now that the mutual gravitational interaction between the two masses is turned off and they are therefore following geodesics of the background spacetime manifold. The spatial separation between the two masses is assumed to be small compared to the wavelengths of the background (external) waves; therefore, the relative motion of the two masses can be expressed via the Jacobi equation. It is useful to express the deviation equation with respect to a Fermi coordinate system established along the worldline of the center of mass of the system. The Fermi system is the most natural extension of a local inertial frame along the path of an observer in spacetime. Hence, the relative acceleration of the two bodies in the Fermi system takes on a Newtonian form that is derivable from a quadrupole potential given by $\frac{1}{2} \epsilon {\cal K}_{ij}(t) x^ix^j$ with \begin{equation}\label{QuadPot} {\cal K}_{ij}(t) = -\frac{1}{2}\frac{\partial^2\chi_{ij}}{\partial t^2} (t,{\bf x}_{\mbox{\rm \tiny CM}})\;, \end{equation} where $(m_1+m_2){\bf x}_{\mbox{\rm \tiny CM}} = m_1{\bf x}_1 + m_2{\bf x}_2$. In equation (\ref{BasicEq}), we have neglected terms of order $\epsilon^2$, $\epsilon\delta$, $\delta^2$, or higher, where $\delta$ is the strength of the radiation reaction term as defined in Appendix~\ref{DEQM}; that is, in conformity with the Newtonian appearance of equation (\ref{BasicEq}) the various forces have been linearly superposed. A detailed treatment of equation (\ref{BasicEq}) with ${\cal R} = 0$ is contained in our previous papers \cite{cmr,cmr1}, which considered only normally incident external waves for the sake of simplicity. We shall also choose here an external periodic monochromatic gravitational wave that is perpendicularly incident on the orbital plane. Let this plane be characterized by the vanishing of the third component of ${\bf r}$. Then, the nonzero elements of the tidal matrix $(\epsilon {\cal K}_{ij})$ are given by \begin{eqnarray}\label{TidalMatrix} {\cal K}_{11} & = & - {\cal K}_{22} = \alpha\Omega^2\cos(\Omega t), \nonumber \\ {\cal K}_{12} & = & {\cal K}_{21} = \beta\Omega^2\cos(\Omega t + \rho), \end{eqnarray} where $\alpha$ and $\beta$ are of the order of unity and represent the constant amplitudes of the two independent linearly polarized components of the incident monochromatic wave of frequency $\Omega$ and $\rho$ is the constant phase difference between the two components (cf. \cite{cmr}). It is interesting to note the partial reciprocity between the emission and absorption of gravitational radiation. In the quadrupole approximation, an elliptical Keplerian binary system emits gravitational waves only at frequencies $m\omega$, $m= 1, 2, 3, \cdots$, where $\omega$ is the Keplerian frequency of the orbit \cite{peters}. It follows from a linear perturbation analysis \cite{mashoon1} that when a Newtonian binary system absorbs gravitational energy from an incident wave of frequency $\Omega$, resonances occur at $\Omega = m\omega$, $m= 1, 2, \cdots$. However, it is important to point out that while the binary monotonically loses energy in the form of gravitational radiation, the absorption of gravitational radiation energy is not monotonic. {\it The system can gain or lose energy in absorption}. In our previous papers \cite{cmr,cmr1}, we examined the long term nonlinear evolution of the dynamical system with ${\cal R} = 0$ and proved the possibility of existence of periodic orbits for which the net flow of energy between the binary and the external periodic gravitational radiation field must vanish. In the present work, we extend our previous results by taking due account of gravitational radiation damping; in fact, periodic orbits are shown to exist for sufficiently small $\delta/\epsilon$ The issue of {\it gravitational ionization} provided the original motivation for our work. The possibility of ionization of a Keplerian binary by incident gravitational waves had been first discussed within the framework of {\it linear} perturbation analysis \cite{mashoon1}. However, the linear perturbation treatment breaks down over time as a consequence of the appearance of secular terms in the analysis. On the other hand, the results of the linear analysis can be applied to the possibility of detection --- via the solar system --- of gravitational pulses \cite{mashoon2} or a stochastic cosmological background of gravitational waves \cite{mashoon3}; a review of this general topic can be found in \cite{grishchuk} and further work in this direction is contained in \cite{Linet,Chau,Petukhov,Ivashchenko,Sbytov} and the references cited therein. Furthermore, the tidal influence of very long wavelength gravitational waves on galaxies has been studied in \cite{Matsuda}. Our {\it nonlinear} perturbation treatment \cite{cmr,cmr1} has led to the conjecture that gravitational ionization is tantamount to Arnold diffusion. Moreover, there are interesting circumstances --- other than periodic orbits --- for which ionization does not occur. For instance, it follows from the KAM theorem that when the incident wave is circularly polarized (i.e. $\alpha=\beta$ and $\rho=\pm \pi/2$) the perturbed binary is forever confined and ionization does not occur if $\epsilon$ is below a certain limit $(\epsilon < \epsilon_{\mbox{\rm \tiny KAM}})$. It is the purpose of this paper to determine how all of these results are affected by the presence of gravitational radiation damping. While our previous work \cite{cmr,cmr1} has involved time-dependent (i.e. nonconservative) Hamiltonian systems with chaotic behavior, the introduction of radiation damping leads to a dissipative dynamical system that is not Hamiltonian. At each instant of time $t$, the state of relative motion is determined by the relative position and velocity $({\bf r}, {\bf v})$; however, it is useful to employ instead the six orbital parameters of the osculating ellipse. That is, the unperturbed bounded motion is in general a Keplerian ellipse; therefore, it is interesting to describe the state of relative motion at each instant by the ellipse that the system would follow if the perturbations were turned off at that instant. The orbital elements of the osculating ellipse are closely related to the Delaunay elements employed in this paper. They are particularly useful in a Hamiltonian system \cite{cmr,cmr1}, since the transformation to the Delaunay action-angle variables is canonical and the corresponding generating function is time-independent so that the magnitude of the Hamiltonian is unchanged under the transformation. A limiting case of the relative osculating ellipse is a circle. In this case, the system emits gravitational waves with a single frequency $2\omega$. On the other hand, linear perturbation theory reveals that resonant absorption occurs at wave frequencies $\omega$, $2\omega$, and $3\omega$. However, it should be remarked that after the passage of time the orbit will no longer be circular. The nonlinear perturbation theory employed in our work is based on the ideas and methods originally developed by Poincar\'{e} \cite{poincare}. An important aspect of this theory is its {\it lack of emphasis} on the precise form of the {\it initial conditions}; clearly, this characteristic of the theory is rather useful in celestial mechanics. In this approach, the dynamical behavior of a physical system is investigated after transients have died away and a steady state has been established. In connection with the issue of initial conditions, it is necessary to remark that the external periodic perturbation in our model is meant to represent the dominant component of an initial wave packet composed mainly of long-wavelength Fourier components consistent with the quadrupole approximation under consideration here. On the other hand, high-frequency waves with wavelengths much smaller that the semimajor axis of the system are expected to have negligible influence on the dynamics. That is, along any significant portion of the relative orbit the waves are expected to go through many oscillations and their net influence is on average expected to be vanishingly small. In our model, the strength of the external perturbation is $\epsilon$; at present, it is expected that $\epsilon \sim 10^{-20}$, though gravitational waves have not yet been detected in the laboratory. On the other hand, the strength of radiative damping in our model is given by $\delta$ (cf. Appendix \ref{DEQM}); for instance, for the binary pulsar PSR 1913+16, discovered by Hulse and Taylor~\cite{hulse}, $\delta \sim 10^{-15}$, while for the Earth-Sun system $\delta \sim 10^{-26}$. The more relativistic a binary system is, the larger the magnitude of $\delta$, which for a physical system must always be less than $(20\sqrt{2})^{-1}$ as discussed in Appendix \ref{DEQM}. Throughout this paper, the external perturbing force is assumed to be due to background gravitational waves. However, as explained in detail in our first paper \cite{cmr}, this is not necessary since any tidal influence on the binary orbit by a far-away mass would contribute a forcing term with the same functional form as in the equation of relative motion (\ref{BasicEq}). In general, the binary system cannot be isolated from the other masses in the universe; this is simply a result of the universal character of gravitational attraction. Thus equation (\ref{BasicEq}), which is essentially the damped Newton-Jacobi equation, should contain the tidal influence of distant masses as well as background gravitational waves in the forcing function. That is, the background spacetime curvature --- represented in equation (\ref{BasicEq}) by the tidal matrix --- is due to the combined gravitational influence of background waves as well as masses. The analysis of the general problem would involve the three-dimensional Kepler system. In our {\em planar} model, we neglect the other masses and choose a normally incident monochromatic plane wave for the sake of simplicity. An interesting feature of equation (\ref{BasicEq}) must be noted here: The vector cross product of this equation with the relative position vector results in an equation for the rate of change of the orbital angular momentum with time. In this result, the torque due to the gravitational radiation reaction is proportional to the orbital angular momentum according to equation (\ref{RadiationReaction}); therefore, the action of this torque is only to change the angular momentum normal to the orbital plane. The torque due to the external force is similar in character, since the external wave is transverse and is assumed to be normally incident on the orbital plane. The net result is that the relative motion is {\it planar}, and this fact considerably simplifies the analysis of the relative motion. Moreover, in conformity with the quadrupole approximation the net velocity of the center of mass (i.e. the velocity of the binary system as a whole) is uniform. Therefore, we choose the ``inertial" reference frame $(t, {\bf x})$ for our model such that it is comoving with the whole binary system and that ${\bf x}_{\mbox{\rm \tiny CM}}=0$. In this frame, the orientation of the coordinate system is so chosen that the perturbed Keplerian motion has a positive initial orbital angular momentum. The motion in our model is planar; therefore, it is interesting to express equation (\ref{BasicEq}) in terms of polar coordinates $(r, \theta)$ in the orbital plane. Thus, \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{k}{r^2} + \frac{p_{\theta}^2}{r^3} + \frac{16G_0^2m_1m_2}{5c^5}\frac{p_r}{r^3}\left(p_r^2+6\frac{p_{\theta}^2}{r^2}+ \frac{4}{3}\frac{k}{r}\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{8G_0^2m_1m_2}{5c^5}\frac{p_{\theta}}{r^3} \left( 9p_r^2 - 6 \frac{p_{\theta}^2}{r^2} + 2\frac{k}{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} Here $p_{\theta}$ is the orbital angular momentum and the initial conditions are chosen such that $p_{\theta}>0$. It is clear from these equations that the corresponding vector field in phase space is periodic with frequency $\Omega$. The main dynamical result of our work may now be stated qualitatively in terms of the characteristic forms of the perturbing functions. The gravitational radiation reaction force is proportional to $r^{-q}$ with $q \geq 3$. Thus if $r$ becomes small, radiative damping takes over and leads the system inexorably to collapse. On the other hand, the external potential is proportional to $r^2$ --- i.e. the energy exchange is proportional to the orbital area through which the normally incident wave passes (cf. \cite{mashoon2}) --- so that if $r$ becomes large the external force takes over and the results of our previous work \cite{cmr,cmr1} imply that Arnold diffusion might take place leading to the ionization of the binary system. Thus it appears that the long time dynamical behavior of the system under consideration here is characterized by two possibilities: collapse to the origin and unbounded growth; in fact, our results indicate that the collapse scenario is prevalent in most cases. The plan of this paper is as follows: In Section 2, we discuss the radiation reaction term ${\cal R}$ in equation (\ref{BasicEq}). Section 3 is devoted to the determination of the bifurcation function and the proof of the existence of periodic orbits in the damped system. In Section 4, the average behavior of the dynamical system (\ref{BasicEq}) with $\epsilon = 0$ is determined for all times. Section 5 describes the average properties of the system (\ref{BasicEq}). Section 6 considers the case of circularly polarized incident waves. Section 7 describes the chaotic net using numerical analysis. A number of detailed computations are relegated to the Appendices. The paper relies on our previous work \cite{cmr,cmr1} for background material regarding the topic of {\it gravitational ionization}; however, we have tried to make this paper essentially self-contained. \section{Gravitational Radiation Damping} In the quadrupole approximation under consideration here, gravitational waves carry energy and angular momentum away from the system but not linear momentum. Similarly, in absorption the orbit exchanges energy and angular momentum with the incident gravitational wave but not linear momentum. This implies that the motion of the system as a whole --- i.e. the center-of-mass motion --- is not affected by the emission and absorption of radiation in the quadrupole approximation. Taking proper account of the emission of gravitational radiation, the conservation laws of energy and momentum require that these losses be reflected in the motion of the system. To satisfy this requirement, it is sufficient to include in the equations of motion of a particle in the system a gravitational radiation reaction force per unit inertial mass given by ${\cal A}$, \begin{equation}\label{RadReaForce1} {\cal A}_i = -\frac{2 G_0}{15 c^5}\frac{d^5 {\cal D}_{ij}}{dt^5}x^j, \end{equation} where $({\cal D}_{ij})$ is the quadrupole moment of the system (cf. \cite{ww} and the references cited therein). We may look upon this force as the spin-2 analog of the standard spin-1 radiation reaction force in electrodynamics \cite{lorentz}. The radiation reaction term in the equations of motion cannot be derived from a potential; otherwise, the dissipative dynamical system would be Hamiltonian. Let ${\cal E}$ denote the energy radiated as gravitational waves; then, \begin{equation}\label{EnergyEq} \frac{d{\cal E}}{dt} = \frac{G_0}{45c^5}\stackrel{...}{{\cal D}}_{ij} \stackrel{...}{{\cal D}}^{ij}, \end{equation} where \begin{displaymath} {\cal D}_{ij} = \int \tilde\rho(3x^ix^j-\delta_{ij}x^2) dV. \end{displaymath} Here the quadrupole moment of the system is traceless by definition, $\tilde\rho$ is the mass density, and $dV$ is the volume element. Equation~(\ref{EnergyEq}) is a standard result of general relativity in the quadrupole approximation; that is, it can be obtained from a detailed treatment of the linearized form of the gravitational field equations together with the physical interpretation of the Landau-Lifshitz pseudotensor. Following this approach, a similar expression can be derived in the quadrupole approximation for the angular momentum carried away by the gravitational waves. The energy and angular momentum radiated away via gravitational waves are lost by the orbit; therefore, the equation of orbital motion should reflect this loss. For both energy and angular momentum, this must be accomplished by the introduction of the radiation reaction force. By analogy with electrodynamics, we expect that the radiation reaction force per unit inertial mass is of the form given by equation (\ref{RadReaForce1}), so that for the binary system under consideration we can express the equations of motion as \begin{equation}\label{RadEq1} \frac{d^2x^i_1}{dt^2} + \frac{G_0m_2({\bf x}_1-{\bf x}_2)^i}{r^3} = -\frac{2G_0}{15c^5}\frac{d^5{\cal D}_{ij}}{dt^5} x^j_1 -\epsilon {\cal K}_{ij}(t) x_1^j, \end{equation} \begin{equation}\label{RadEq2} \frac{d^2x^i_2}{dt^2} + \frac{G_0m_1({\bf x}_2-{\bf x}_1)^i}{r^3} = -\frac{2G_0}{15c^5}\frac{d^5{\cal D}_{ij}}{dt^5} x^j_2 -\epsilon {\cal K}_{ij}(t) x_2^j, \end{equation} where ${\bf r} = {\bf x}_1-{\bf x}_2$. Multiplying the first equation by $m_1$ and the second equation by $m_2$ and adding them results in \[ \frac{d^2 x^i_{\mbox{\rm\tiny CM}}}{dt^2}=-\left[ \frac{2G_0}{15c^5}\frac{d^5{\cal D}_{ij}}{dt^5} +\epsilon {\cal K}_{ij}(t) \right] x^j_{\mbox{\rm \tiny CM}}. \] In the absence of perturbations, ${\bf x}_{\mbox{\rm \tiny CM}}$ is at the origin of the coordinate system; therefore, this equation of motion for ${\bf x}_{\mbox{\rm \tiny CM}}$ implies that it is effectively unchanged in our approximate treatment. It is thus possible to set ${\bf x}_{\mbox{\rm \tiny CM}} = 0$, so that the emission or absorption of gravitational waves does not affect the total linear momentum of the system in the quadrupole approximation. The relative motion is obtained by subtracting the two equations: \begin{equation}\label{RadReaEq} \frac{d^2 r^i}{dt^2} + \frac{G_0(m_1+m_2)r^i}{r^3} = -\frac{2 G_0}{15 c^5}\frac{d^5{\cal D}_{ij}}{dt^5} r^j -\epsilon {\cal K}_{ij}(t) r^j. \end{equation} It follows that the full content of the equations of motion is simply contained in the expressions ${\bf x}_1 = m_2{\bf r}/M$ and ${\bf x}_2 = -m_1{\bf r}/M$, where $M = m_1 + m_2$ is the total mass and ${\bf r}$ is given by the equation of relative motion (\ref{RadReaEq}). The quadrupole moment of the system is then given by \begin{equation}\label{QuadMom} {\cal D}_{ij} = m_1(3x^i_1x^j_1-\delta_{ij}x^2_1) + m_2(3x^i_2x^j_2-\delta_{ij} x^2_2) = \mu(3r^ir^j-\delta_{ij}r^2), \end{equation} where $\mu$ is the reduced mass, i.e. $\mu = m_1m_2/M$. It is clear from the form of equations (\ref{RadReaEq}) and (\ref{QuadMom}) that the relative motion must be determined from an equation in which the order of differentiation with respect to time exceeds two. On the other hand, the higher derivative terms only appear in the radiation reaction acceleration that has an overall dimensionless strength given by a small parameter $\delta$ that is particularly small for Keplerian (i.e. nonrelativistic) binaries as discussed in Appendix \ref{DEQM}. It follows that the problem can be avoided in this case if we proceed iteratively and substitute for $\ddot{r}^i$ from equation~(\ref{RadReaEq}) every time it appears in the calculation of gravitational radiation damping. Then, keeping only terms of the first order in the perturbation parameters $\epsilon$ and $\delta$ in equation (\ref{RadReaEq}) we find that \begin{eqnarray}\label{RadDamping} \frac{1}{\mu}\frac{d^5{\cal D}_{ij}}{dt^5} &=& -2\frac{k\dot{r}}{r^4} \left( 9v^2-15\dot{r}^2-\frac{8k}{r}\right) \delta_{ij} -90\frac{k\dot{r}}{r^6}(3v^2-7\dot{r}^2)r^ir^j \nonumber \\ && + 24\frac{k}{r^5}\left(3v^2-15\dot{r}^2-\frac{k}{r}\right)(r^iv^j+r^jv^i) +180\frac{k\dot{r}}{r^4}v^iv^j, \end{eqnarray} which is the approximate expression to be used in the radiative damping term \begin{equation}\label{RadRea2A} {\cal R}^i = \frac{2 G_0}{15 c^5} \frac{d^5{\cal D}_{ij}}{dt^5} r^j \end{equation} in order to obtain equation (\ref{RadiationReaction}). In deriving equation~(\ref{RadDamping}), we have repeatedly used the fact that ${\bf r}\cdot{\bf v} = r\dot{r}$. In the following sections, we shall use equations~(\ref{BasicEq}) and~(\ref{RadiationReaction}) to study the dynamical behavior of the Kepler system over a long period of time. Let us now return to equation~(\ref{RadReaEq}) and point out that this equation of relative motion --- in which the effects of emission and absorption of gravitational radiation have been taken into account in the quadrupole approximation and all other retardation effects have been neglected --- is consistent with the conservation laws of energy and momentum. We will illustrate this point explicitly for the energy of the system; the case of angular momentum is similar but will not be treated here. Let us therefore derive an energy expression for the equation of relative motion (\ref{RadReaEq}). Multiplying both sides of that equation by the reduced mass and the relative velocity results in \begin{equation}\label{RelEng1} \frac{d}{dt}\left[ E + \frac{G_0}{45 c^5} \left(\stackrel{....}{\cal D}_{ij}\dot{\cal D}^{ij}- \stackrel{...}{\cal D}_{ij}\ddot{\cal D}^{ij}\right)\right] = -\frac{G_0}{45 c^5} \stackrel{...}{\cal D}_{ij}\stackrel{...}{\cal D}^{ij} -\frac{1}{6}\epsilon{\cal K}_{ij}\dot{\cal D}^{ij}, \end{equation} where we have used $\dot{\cal D}_{ij} = \mu\left[3(v^ir^j+r^iv^j)-2{\bf r}\cdot{\bf v}\:\delta_{ij}\right]$ and the fact that $({\cal D}_{ij})$ is symmetric and traceless. The first term in the square brackets in equation (\ref{RelEng1}) is the Keplerian energy of the relative orbit $E$, $2E=\mu v^2-2k\mu/r$, and the second term is a relativistic contribution. The rate at which this latter term varies averages out --- over the period of the unperturbed Keplerian motion --- to a quantity that is negligible at the level of approximation under consideration in this paper. It follows that the average rate of loss of Keplerian energy of the dynamical system by radiation damping is equal to the average rate at which energy leaves the orbit via gravitational waves, as expected. This latter rate is given by averaging equation (\ref{EnergyEq}) over one period of elliptical Keplerian motion and the result is \cite{peters} \begin{equation} \left<\frac{d{\cal E}}{dt}\right> = \frac{32G_0^4m_1^2m_2^2(m_1+m_2)}{5c^5a^5(1-e^2)^{7/2}} \left(1+\frac{73}{24}e^2+ \frac{37}{96}e^4\right), \end{equation} where $a$ and $e$ are the semimajor axis and the eccentricity of the Keplerian ellipse, respectively. At this average rate, gravitational radiation energy permanently leaves the system and goes off to infinity. No energy is actually lost; the orbital energy is simply converted into gravitational radiation energy. It is sometimes stated (cf. \cite{ww} and references therein) that the force due to gravitational radiation damping could be derived from a quadrupole potential. If so, the dynamical system under consideration here would be Hamiltonian. That is, a simple comparison of the two perturbing influences in equation (\ref{RadReaEq}) reveals that this system would be Hamiltonian if ${d^5{\cal D}_{ij}}/{dt^5}$ were simply a function of time. However, this cannot be the case since ${d^5{\cal D}_{ij}}/{dt^5}$ is a function of ${\bf r}$ and ${\bf v}$; moreover, the system is manifestly dissipative. \section{Continuation of Periodic Orbits}\label{conpo} The investigation of the dynamics of the perturbed Kepler system (\ref{BasicEq}) can be simplified considerably if the equations of relative motion are expressed in terms of action-angle variables that are suited to the unperturbed system \cite{Kov,Stern}. If the perturbations of the Kepler system (\ref{MathEQM}) are turned off at a given instant, the relative motion will follow an osculating elliptical orbit with semimajor axis $a$, eccentricity $e$, eccentric anomaly $\hat{u}$, and true anomaly $\hat{v}$. The orientation of the underlying coordinate system can be so chosen that the relative orbital angular momentum is positive. The action-angle variables appropriate to the Kepler problem are the Delaunay elements \cite{Kov,Stern} that are given for the planar problem under consideration by \begin{eqnarray}\label{PProblem} L := (ka)^{1/2},\qquad & & G:= p_{\theta} = L(1-e^2)^{1/2}, \nonumber \\ \ell:= \hat{u} - e \sin{\hat{u}}, \qquad & & g:=\theta - \hat{v}. \end{eqnarray} To analyze the dynamical behavior contained in equation (\ref{BasicEq}) in terms of Delaunay variables, we first put these equations in dimensionless form. This is done in Appendix \ref{DEQM}. Then, the dimensionless equations of motion are expressed in terms of Delaunay elements; this transformation is presented in Appendix \ref{EQMDelaunay}. This results in the following form of the dynamical equations: \begin{eqnarray}\label{D2EQM1} \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 \begin{eqnarray}\label{feq1} 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} and $f_L$,$f_G$, $f_{\ell}$, $f_g$ can be expressed in terms of Delaunay elements using classical methods of celestial mechanics involving the Bessel functions (cf.\ Appendix \ref{EQMDelaunay}). The right hand side of the system (\ref{D2EQM1}) is periodic in time with frequency $\Omega$. For $\delta = 0$, the system (\ref{D2EQM1}) has periodic orbits as proved in \cite{cmr}; moreover, this Hamiltonian system exhibits dynamical behavior \cite{cmr,cmr1} that appears to be characteristic of Arnold diffusion \cite{arnold,liberman}. In the following, we show that when radiation reaction is taken into account the periodic orbits persist for sufficiently small radiative damping. For $\epsilon=0$, on the other hand, the system~(\ref{D2EQM1}) continuously loses energy to gravitational radiation and eventual collapse is inevitable; moreover, this dynamical behavior is structurally stable and is therefore expected to persist for sufficiently small $\epsilon\ne 0$. Let us now assume that there are relatively prime positive integers $m$ and $n$ such that $m\omega = n\Omega$, i.e. the unperturbed Keplerian orbit is in resonance with the disturbing function. Then, it turns out that the three-dimensional resonance manifold (period manifold) \[{\cal Z}^L = \{(L, G, \ell, g): m\omega = n\Omega\} \] is a normally nondegenerate manifold~\cite{cmr}. Moreover, the solution of the unperturbed system starting at $(L_0, G_0, \ell_0, g_0)$ is given by \[t \mapsto (L_0, G_0, \omega \bar{t} + \ell_0, g_0),\]where $\bar{t} = t - t_0$ and $t_0$ is the initial time. It turns out that we can set $t_0 = 0$ with no loss of generality. It follows from \cite{cmr,cmr1,ccc} that we must project the partial derivative of the Poincar\'{e} map with respect to $\epsilon$ onto the complement of the range of the infinitesimal displacement. The partial derivative of the Poincar\'{e} map with respect to $\epsilon$ at $\epsilon = 0$ is obtained from the solution of the second variational initial value problem \begin{eqnarray}\label{IVP1} \dot{L}_{\epsilon} & = & - \left(\frac{\partial{\cal C}}{\partial \ell} \phi(t) + \frac{\partial{\cal S}}{\partial \ell} \psi(t)\right) + \Delta f_L, \nonumber \\ \dot{G}_{\epsilon} & = & - \left(\frac{\partial{\cal C}}{\partial g} \phi(t) + \frac{\partial{\cal S}}{\partial g} \psi(t)\right) + \Delta f_G , \nonumber \\ \dot{\ell}_{\epsilon} & = & -\frac{3}{L^4}L_{\epsilon} + \left(\frac{\partial{\cal C}}{\partial L} \phi(t) + \frac{\partial{\cal S}}{\partial L} \psi(t)\right) + \Delta f_{\ell}, \nonumber \\ \dot{g}_{\epsilon} & = & \left(\frac{\partial{\cal C}}{\partial G} \phi(t) + \frac{\partial{\cal S}}{\partial G}\psi(t)\right) + \Delta f_g, \end{eqnarray} with \begin{eqnarray}\label{IVP2} L_{\epsilon}(0, L_0, G_0, \ell_0, g_0) = 0, & \hspace*{.25in} & G_{\epsilon}(0, L_0, G_0, \ell_0, g_0) = 0, \nonumber \\ \ell_{\epsilon}(0, L_0, G_0, \ell_0, g_0) = 0, & \hspace*{.25in} & g_{\epsilon} (0, L_0, G_0, \ell_0, g_0) = 0, \end{eqnarray} where $\Delta = \delta/\epsilon$ so that $0<\Delta<\infty$. In fact, the partial derivative is given by the vector \[\left[L_{\epsilon}\left(m\frac{2\pi}{\Omega}, L_0, G_0, \ell_0, g_0 \right), G_{\epsilon}\left(m\frac{2\pi}{\Omega}, L_0, G_0,\ell_0, g_0 \right),\right. \hspace*{1in}\] \[ \hspace*{1in}\left. \ell_{\epsilon}\left(m\frac{2\pi}{\Omega}, L_0, G_0, \ell_0, g_0 \right), g_{\epsilon}\left(m\frac{2\pi}{\Omega}, L_0, G_0,\ell_0, g_0 \right)\right] .\] Here the independent variables are $(t,L,G,\ell,g)$ and the last four variables of each component function give the initial point for the original solution of the unperturbed system, i.e. equations (\ref{D2EQM1}) with $\epsilon = 0$ and $\delta=0$. The appropriate projection is onto the first, second and fourth coordinates. In fact, the bifurcation function is this projection restricted to ${\cal Z}^L$, and is given by \begin{eqnarray}\label{BifurcationEq} {\cal B} & = & \left[L_{\epsilon}\left(m\frac{2\pi}{\Omega}, L_0, G_0, \ell_0, g_0 \right), G_{\epsilon}\left(m\frac{2\pi}{\Omega}, L_0, G_0,\ell_0, g_0 \right), \right. \nonumber \\ & & \hspace*{.5in} \left. g_{\epsilon}\left(m\frac{2\pi}{\Omega}, L_0, G_0, \ell_0, g_0 \right)\right], \end{eqnarray} where $L=L_0$ is fixed by the choice of resonance, i.e. $m/L_0^3 = n \Omega,$ and $G_0>0$ by our choice of spatial coordinate axes. Dropping subscripts indicating initial values of Delaunay elements for the sake of simplicity, the components of the bifurcation function (\ref{BifurcationEq}) can now be written as \begin{eqnarray}\label{BifrucationExp} {\cal B}^L(G, \ell, g) & = & -\frac{\partial {\cal I}}{\partial \ell} +\Delta\int_0^{2\pi m/\Omega}f_L(L,G,\omega t+\ell,g)\,dt, \nonumber \\ {\cal B}^G(G, \ell, g) & = & -\frac{\partial {\cal I}}{\partial g} +\Delta\int_0^{2\pi m/\Omega}f_G(L,G,\omega t+\ell,g)\,dt, \nonumber \\ {\cal B}^g(G, \ell, g) & = & \frac{\partial {\cal I}}{\partial G} +\Delta\int_0^{2\pi m/\Omega}f_g(L,G,\omega t+\ell,g)\,dt, \end{eqnarray} where ${\cal I}$ is defined by \begin{equation}\label{IEq1} {\cal I} = \int_0^{2\pi m/\Omega} [{\cal C}(L, G, \omega t + \ell, g)\phi(t) + {\cal S}(L, G, \omega t + \ell, g)\psi(t)]\: dt, \end{equation} and $\cal C$ and $\cal S$ are given in Appendix~\ref{EQMDelaunay}. The integral in (\ref{IEq1}) has been calculated in \cite{cmr} with the result that ${\cal I} = 0$ for $n > 1$, while for $n = 1$, \begin{eqnarray}\label{Ieqn3} {\cal I} & = & \frac{1}{2}\pi ma^2\Omega \{ \alpha(A_m \cos{m \ell} \cos{2g} - B_m \sin{m \ell} \sin{2g}) \hspace*{1in} \nonumber \\ & & \hspace*{.5in} + \beta [A_m \cos{(m\ell - \rho)} \sin{2g} + B_m \sin{(m\ell- \rho)} \cos{2g}] \}. \end{eqnarray} The periodic orbits that continue are determined by the simple zeros of the bifurcation function (\ref{BifrucationExp}). The integrals multiplying $\Delta$ in equations (\ref{BifrucationExp}) are given in Appendix \ref{RadForceApp}; therefore, periodic orbits continue to exist when the following equations are satisfied by simple zeros: \begin{eqnarray}\label{ZerosBif1} \frac{\partial {\cal I}}{\partial \ell} + 2\pi \Delta \frac{L^3}{G^7}\left( 8 + \frac{73}{3} e^2 + \frac{37}{12} e^4 \right) & = & 0, \nonumber \\ \frac{\partial {\cal I}}{\partial g} + 2\pi \Delta \frac{1}{G^4} ( 8 + 7 e^2) & = & 0, \nonumber \\ \frac{\partial {\cal I}}{\partial G} & = & 0. \end{eqnarray} For $n > 1$, ${\cal I} = 0$ and hence these equations have no solution as the expressions multiplying $\Delta$ are manifestly positive. For $n = 1$ and $\Delta = 0$, we have proved in \cite{cmr} that there are simple zeros of the bifurcation function for all positive integers $m$. It follows that equation (\ref{ZerosBif1}) must have simple zeros as well for sufficiently small $\Delta$. It is important to remark here that this conclusion would not be altered by the inclusion of terms of order $\epsilon^2$, $\epsilon\delta$, $\delta^2$, or higher in the original system (\ref{BasicEq}). In fact, the incorporation of such higher order effects could only affect the shape of a periodic orbit of the type investigated here but not its existence. The Kepler system under investigation in this paper may in general have other periodic orbits than those revealed by our first order method. The investigation of the stability of the periodic orbits under consideration is beyond the scope of this work. If there is a stable periodic orbit for $\epsilon>0$ and $\delta>0$, then its (open) basin of attraction consists of trajectories that are permanently captured into resonance. This same issue will be discussed in connection with the averaging methods that we introduce in the next section. \section{Averaged Dynamics}\label{avedyn} In this section, we wish to consider our system in Delaunay elements when no external gravitational waves are present. The motivation for this discussion is the experimental observation of the inward spiraling of the members of the Hulse-Taylor binary pulsar \cite{hulse}. Thus we consider the traditional Kepler problem except that the post-Newtonian effect of gravitational radiation damping is taken into account in the quadrupole approximation. We expect that the system remains planar and loses energy and angular momentum such that eventually $L\to 0$, $G \to 0$, and the system collapses. The behavior of the system in the infinite past is less trivial. The precise manner in which the orbit behaves on average for $t\to\pm\infty$ can be determined using the method of averaging. In this case, we have \begin{eqnarray}\label{HulseEQM} \dot{L} = \delta f_L, & &\qquad \dot{G} = \delta f_G, \nonumber \\ \dot{\ell} = \omega + \delta f_{\ell}, & & \qquad \dot{g} = \delta f_g. \end{eqnarray} Equations (\ref{HulseEQM}) are of the general form \begin{equation}\label{HulseEQM1} \dot{I} = \epsilon f(I, \varphi, \epsilon), \qquad \dot{\varphi} = \omega(I) + \epsilon g(I, \varphi, \epsilon), \end{equation} where both $f$ and $g$ are $2\pi$ periodic in the angle variable $\varphi$. In fact, in our case $I = (L, G, g)$ and $\varphi = \ell$. The Averaging Theorem~\cite{arnold2} asserts the following: Suppose that $t \mapsto (I(t), \varphi(t))$ is a solution of equations (\ref{HulseEQM1}) and $\omega (I(t)) > 0$ for all $t$. If $t \mapsto J(t)$ is the solution of the averaged system \[ \dot{J} = \epsilon \frac{1}{2\pi} \int_0^{2\pi} f(J, \varphi, 0) d\varphi \] with $J(0) = I(0)$, then, for sufficiently small $\epsilon$, there is a constant $C$ independent of $\epsilon$ such that \[ |I(t) - J(t)| < C\epsilon \qquad \mbox{\rm for} \qquad 0 \leq t \leq \frac{1}{\epsilon}. \] That is, the ``action variables" of the original system (\ref{HulseEQM1}) remain close to the solution of the averaged system over a long time-scale. It is important to note that the Averaging Theorem, in the generality stated above, is only true for the case of a single frequency, i.e. there is only one angle variable and its frequency does not vanish. If there are more than one angle present, then resonances among the frequencies must be taken into account. Though $g$ is an angle variable, its evolution is slow and its average rate of variation vanishes; hence, it may be regarded as an ``action'' variable for the purposes of this section. Therefore in our case, the averaged system obtained from (\ref{HulseEQM}), where we use $L$, $G$, and $g$ for the averaged variables, $e = \sqrt{L^2-G^2}/L$, and $\Delta=\delta/\epsilon$, is given by (cf. Appendix~\ref{RadForceApp}) \begin{eqnarray}\label{AveEQM} \dot{L} & = & -\epsilon \frac{\Delta}{G^7} \left(8 + \frac{73}{3}e^2 + \frac{37}{12}e^4\right) , \nonumber \\ \dot{G} & = & -\epsilon \frac{\Delta}{L^3G^4}(8 + 7 e^2), \nonumber \\ \dot{g} & = & 0. \end{eqnarray} To determine the dynamics of the averaged system (\ref{AveEQM}) for the variables $L$ and $G$, we note that equations (\ref{AveEQM}) are autonomous. Of course, the system is in this case singular on the (invariant) $L$-axis. Also, recall that it suffices to consider only the case where $G>0$. After multiplication by $12L^5G^7/\epsilon\Delta$, we obtain the dynamically equivalent system of equations \begin{eqnarray}\label{AveEQM2} \dot{L} & = & -L(37G^4-366G^2L^2+425L^4), \nonumber \\ \dot{G} & = & -G^3(180L^2-84G^2). \end{eqnarray} The Delaunay elements are defined only for $L \geq 0$ and $0 \epsilon_{\mbox{\rm \tiny KAM}}$, the situation is less clear and --- as discussed in previous sections --- there could be an interesting interplay between collapse and Arnold diffusion. It must be pointed out that these results are independent of the frequency of incident radiation $\Omega$ and only require that $\alpha = \beta$ and $\rho = \pm\pi/2$. The dynamical system can be viewed from a frame rotating with frequency $\Omega /2$; in this frame, there exists an energy integral once the radiative damping is neglected. That is, the incident gravitational wave stands completely still as a consequence of the gravitational helicity-rotation coupling \cite{mashhoon3} and hence the dynamical system is time independent in the absence of damping. The KAM theorem would then imply that for sufficiently small $\epsilon$ the system is forever bounded. This confinement brings into rather sharp focus the nonreciprocal nature of emission and absorption of gravitational radiation. Clearly, the continuous loss of energy by the binary would eventually result in the collapse of the system. \section{Transient Chaos} The equations of motion~(\ref{MathEQM}) can be integrated numerically given a set of initial conditions at $t = 0$. For this purpose it is convenient to express these equations in dimensionless form using the scale transformations discussed in Appendix~\ref{DEQM}. The actual system that we use is given by \begin{eqnarray}\label{NumExpEq} \frac{dr}{dt}& = &p_r, \qquad \frac{d\theta}{dt} = \frac{p_{\theta}}{r^2},\nonumber\\ \frac{dp_r}{dt} & = & -\frac{1}{r^2} + \frac{p_{\theta}^2}{r^3}+4\delta \frac{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} & = & 2\delta \frac{p_{\theta}}{r^3}\left(9 p_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} We remark that this system is numerically ill-conditioned near collision, i.e. near $r=0$. Unfortunately, there does not seem to be a way to regularize the system for numerical integration as can be done for regular perturbations of Kepler motion as discussed in~\cite{Stiefel}. The problem is that the two perturbation terms in our system enter with negative and positive powers of the variable $r$. A possible set of initial conditions for~(\ref{NumExpEq}) that we use for our numerical experiments in this work is $(p_r, p_{\theta}, r, \theta) = (e, 1, 1, 0)$, which corresponds to an osculating ellipse at $t = 0$ with unit orbital angular momentum, eccentricity $e$, semimajor axis $a=(1-e^2)^{-1}$, and Keplerian frequency $\omega = (1-e^2)^{3/2}$. The true anomaly of this osculating ellipse is given by $\hat{v} = \pi/2$, the Delaunay element $g = -\pi/2$, and the periastron occurs at $\theta = -\pi/2$. The main conclusion of the analysis in this paper is that for sufficiently small $\epsilon$, the majority of initial conditions lead {\em on average} to collapse. However, from basic results in nonlinear dynamics, we expect that near resonances there would be transverse intersections of stable and unstable manifolds of some of the perturbed periodic orbits. This suggests the existence of transient chaos --- chaotic invariant sets that are not attractors. We will see in a moment that numerical experiments are consistent with this observation. However, since we do not have estimates on the size of the perturbation parameters for which the expected behavior is valid, there does not seem to be a mathematical obstruction to the existence of attracting chaotic sets (strange attractors) for some choices of the perturbation parameters. In this regard, we mention the connection between the planar Kepler system and a two-dimensional anharmonic system with an external periodic force and damping~\cite{Stiefel}. The existence of transient chaos for small amplitude forcing and the possibility of the existence of strange attractors for larger amplitude forcing and damping is an interesting and well documented dynamical feature of the systems in the latter class. In fact, the existence of transient chaos and strange attractors for the case of a single oscillator is well known (see, for example,~\cite{wig}). It should also be mentioned that the connection between the Kepler system and oscillators is quite general; in this paper, we have restricted our attention to the planar Kepler problem for the sake of simplicity. It is possible that the three-dimensional Kepler system given by the general form of equations (\ref{BasicEq}) and (\ref{RadiationReaction}) would exhibit novel features not found in the planar case. We have conducted a {\it preliminary} numerical search for a strange attractor in our perturbed Kepler system, but we have not succeeded in finding such an attractor. This is an interesting problem for further research. As a numerical example consistent with transient chaotic behavior, we refer to Figure~\ref{trans} where the evolution of the action variables, strobed at each cycle of the incident wave, is depicted. While both action variables tend to eventual collapse, complicated transient dynamics is indicated in this case. It is interesting to note that here the orbit passes through a dense set of $(m:n)$ resonances that may be responsible for the complicated behavior of these action variables. The evolution of these variables in the bottom panel should be compared and contrasted with that in the bottom panel of Figure~\ref{capfig}. Finally, we remark that while energy continuously leaves the system via the emission of gravitational waves, the definite trend toward collapse that results may --- under special circumstances --- be counter-balanced by the input of energy into the system by the external periodic perturbations. If the net flow of energy into the system is positive on average, gravitational ionization will take place---the semimajor axis of the osculating ellipse will grow, albeit not necessarily monotonically. This balance seems to depend very sensitively on the choice of parameters and initial conditions. For instance, let the initial conditions be $(p_r, p_{\theta}, r, \theta) = (0.5, 1, 1, 0)$ in the system~(\ref{NumExpEq}) and set the parameters such that $\alpha = 2$, $\beta = 2$, $\rho = 0$, and $\Omega = 1.299038106$; for $\epsilon = 0.005$ and $\delta/\epsilon = 0.00198992$ the system appears to ionize over a long time-scale, while for $\delta/\epsilon = 0.00198993$ the system collapses. \appendix \section{Scale Transformations}\label{DEQM} To facilitate the analysis of the dynamical behavior of the Kepler problem with both incident gravitational waves and gravitational radiation damping, we transform the equations of motion for this system to dimensionless form. To this end, let $r^i = R_0 \hat{r}^i$ and $t = T_0 \hat{t}$, where $R_0$ and $T_0$ are scale parameters and $\hat{r}^i$ and $\hat{t}$ are dimensionless. This will be true for all quantities with a hat. Substituting for $r^i$ and $t$ in equations~(\ref{BasicEq}) and~(\ref{RadiationReaction}) leads to \begin{eqnarray}\label{HatEqOfMotion} \frac{d^2\hat{r}^i}{d\hat{t}^2} + \hat{k}\frac{\hat{r}^i}{\hat{r}^3} + \frac{\delta}{\hat{r}^3}\left[\left(12 \hat{v}^2 -30\dot{\hat{r}}^2- 4\frac{\hat{k}}{\hat{r}}\right)\frac{d\hat{r}^i}{d\hat{t}} - \frac{\dot{\hat{r}}} {\hat{r}}\left(36\hat{v}^2 \right. \right. & & \nonumber \\ \left. \left. -50\dot{\hat{r}}^2 + \frac{4}{3}\frac{\hat{k}}{\hat{r}} \right) \hat{r}^i\right] + \epsilon \hat{{\cal K}}_{ij}\hat{r}^j & = & 0, \end{eqnarray} where \begin{displaymath} \hat{k} := k\frac{T_0^2}{R_0^3}, \qquad \delta := \frac{4}{5} \frac{G_0^2m_1m_2}{c^5T_0R_0}, \qquad \mbox{\rm and} \qquad \hat{\cal K}_{ij} := T_0^2{\cal K}_{ij}. \end{displaymath} Let us now consider the physical interpretation of such a transformation. If we take a binary system with initial period $2\pi T_0$ and semimajor axis $R_0$, then by Kepler's law \[ \omega^2 = \frac{G_0(m_1+m_2)}{R_0^3} = \frac{1}{T_0^2}, \] i.e. $\hat{k} = 1$. In general, one can define $\hat{\omega}=\omega T_0$ and $\hat{a}=a/R_0$; then, Kepler's law becomes $\hat{\omega}^2=\hat{k}/\hat{a}^3$, where $\hat{a}$ is the dimensionless semimajor axis. We also note that $\delta$ is dimensionless and turns out to be a small parameter, $0<\delta<<1$, for all realistic binary systems. By setting $\hat{k} = 1$ and dropping all hats in equation (\ref{HatEqOfMotion}), we obtain \begin{eqnarray}\label{FinalEqOfMotion} \frac{d^2 r^i}{dt^2} + \frac{r^i}{r^3} + \frac{\delta}{r^3} \left[\left(12v^2 -30\dot{r}^2-\frac{4}{r}\right)\frac{dr^i}{dt} -\frac{\dot{r}}{r}\left(36v^2 \right. \right. & & \nonumber \\ \left. \left. -50\dot{r}^2 + \frac{4}{3r} \right)r^i\right] + \epsilon {\cal K}_{ij}r^j & = & 0, \end{eqnarray} which is the form of the equation of motion that we investigate in this paper. Now everything is dimensionless in equation (\ref{FinalEqOfMotion}); there are two small parameters $\epsilon$ and $\delta$ and because $\hat{k}=1$, the underlying scales $R_0$ and $T_0$ are connected such that $R_0^3 = k T_0^2$. Furthermore, if we take the direction of incidence of the external gravitational waves to be perpendicular to the initial orbital plane, then it follows directly from~(\ref{FinalEqOfMotion}) that the damped motion is planar. Let us now estimate the magnitude of $\delta$ for binaries of physical interest. Expressing $\delta$ in terms of the total mass $M = m_1+m_2$ and the reduced mass $\mu = m_1m_2/M$, we find \[\delta = \frac{4}{5}\left(\frac{\mu}{M}\right)\left(\frac{R_0}{cT_0}\right)\left( \frac{G_0M}{c^2R_0}\right)^2.\] Furthermore, our assumption that $\hat{k}=1$ implies \[\frac{R_0}{cT_0} = \left(\frac{G_0M}{c^2R_0}\right)^{1/2};\] hence, we have \[\delta = \frac{1}{5}\left(\frac{4\mu}{M}\right)\left(\frac{G_0M}{c^2R_0} \right)^{5/2}.\] Here $4\mu/M \leq 1$ by definition, and $G_0M/c^2R_0$ is always less than $1/2$. The quantity $2G_0M/c^2R_0$ is the ratio of the Schwarzschild radius of the system to the semimajor axis of the binary; thus, the system is nearly a black hole for $2G_0M/c^2R_0 \sim 1$, while for a Keplerian binary $2G_0M/c^2R_0<<1$, which indicates that the semimajor axis of the orbit is much larger than the gravitational radius of the binary system. It follows from these remarks that for a physical system $\delta < 0.04$. For the Earth-Sun system, for instance, $\mu/M \simeq 3 \times 10^{-6}$ and $G_0M/c^2R_0 \simeq 10^{-8}$, so that $\delta \simeq 10^{-26}$; by comparison, $\delta$ is negligibly small for an artificial satellite in orbit about the Earth. On the other hand, for the Hulse-Taylor binary pulsar, PSR 1913+16, we find that $\delta \simeq 10^{-15}$. \section{Equations of Motion in Delaunay \newline Elements}\label{EQMDelaunay} Let us now transform the equations of motion (\ref{FinalEqOfMotion}) from Cartesian coordinates and velocities to Delaunay elements. The standard derivation of the equations of motion in Delaunay variables is based on the assumption that the system under consideration is Hamiltonian. However, the dynamical system in our treatment is dissipative; therefore, in this appendix we present a general derivation for the planar Kepler problem. This transformation is performed in two steps. First equation (\ref{FinalEqOfMotion}) is expressed in polar coordinates and then in Delaunay elements. Let \[ f_i = -\frac{1}{r^3}\left[\left(12v^2-30\dot{r}^2-\frac{4}{r}\right) \frac{dr^i}{dt} -\frac{\dot{r}}{r}\left(36v^2-50\dot{r}^2+\frac{4}{3r}\right)r^i \right] \] and ${\bf r} = (r\cos{\theta}, r\sin{\theta}, 0)$; then, equation (\ref{FinalEqOfMotion}) can be written in polar coordinates as \begin{eqnarray}\label{PolarEQM} \ddot{r} - r\dot{\theta}^2 + \frac{1}{r^2} &=&\delta f_r - \epsilon {\cal K}_{r}, \nonumber \\ r\ddot{\theta} +2\dot{r}\dot{\theta} & = & \delta f_{\theta} - \epsilon {\cal K}_{\theta}, \end{eqnarray} where \begin{displaymath} f_r = f_1 \cos{\theta} + f_2 \sin{\theta},\qquad f_{\theta} = -f_1 \sin{\theta} + f_2 \cos{\theta}, \end{displaymath} \begin{equation}\label{frftheta} {\cal K}_r = r({\cal K}_{11}\cos{2\theta}+{\cal K}_{12}\sin{2\theta}), \qquad {\cal K}_{\theta} = r(-{\cal K}_{11}\sin{2\theta}+{\cal K}_{12} \cos{2\theta}). \end{equation} In determining ${\cal K}_{r}$ and ${\cal K}_{\theta}$, we have utilized the fact that the matrix $({\cal K}_{ij})$ is symmetric and traceless. In what follows, we let \begin{equation}\label{Hatdefn} F_r: = \delta f_r - \epsilon {\cal K}_{r},\qquad F_{\theta}:= \delta f_{\theta} - \epsilon {\cal K}_{\theta} \end{equation} for the sake of simplicity. The Delaunay elements are defined by the ellipse tangent to the path at time $t$. The state of relative motion is specified by $(r,\theta,p_r,p_\theta)$ in polar coordinates; therefore, using the standard formulae for elliptic motion we have $p_rp_\theta=e\sin\hat{v}$ and $p_\theta^2/r=1+e\cos\hat {v}$, which uniquely specify the eccentricity $e$ and the true anomaly $\hat{v}$ of the osculating ellipse. Next, the semimajor axis of this ellipse $a$ is obtained from $p_\theta^2=a(1-e^2)$ and the eccentric anomaly $\hat {u}$ can be determined from \begin{displaymath} \cos{\hat{v}} = \frac{\cos{\hat{u}}-e}{1-e\cos{\hat{u}}}\:, \qquad \sin{\hat{v}}=\frac{\sqrt{1-e^2}\sin{\hat{u}}}{1-e\cos{\hat{u}}}. \end{displaymath} The Delaunay elements are then given by \begin{eqnarray}\label{DelaunayElements} L = \sqrt{a},\: & & \: G = p_{\theta}, \nonumber \\ \ell =\hat{u}-e\sin{\hat{u}},\: & & \: g=\theta-\hat{v}. \end{eqnarray} To obtain the equations of motion in Delaunay elements, we observe that the energy and orbital angular momentum of the relative motion per unit mass are \begin{eqnarray}\label{EllipseEnergy} \frac{E}{\mu} & = & \frac{1}{2}(\dot{r}^2+r^2\dot{\theta}^2)-\frac{1}{r} = -\frac{1}{2a}, \nonumber \\ p_{\theta} & = & r^2\dot \theta= \sqrt{ a(1-e^2)}, \end{eqnarray} where the orientation of the spatial coordinate system is chosen such that $p_\theta>0.$ Now \begin{displaymath} \frac{1}{\mu}\frac{dE}{dt} = {\bf F} \cdot {\bf v}, \end{displaymath} where $ {\bf F}\cdot {\bf v} = F_r \: \dot{r} + F_{\theta} \: r\dot{\theta} $. It follows that \begin{equation}\label{Leqn} \frac{dL}{dt} = \frac{a}{\sqrt{1-e^2}} \left[F_r e \sin{\hat{v}} + F_{\theta} \frac{a(1-e^2)}{r}\right], \end{equation} where we have replaced $\dot{r}$ by $e \sin{\hat{v}}/p_{\theta}$, which follows from equation (\ref{EllipseEnergy}). It is a direct consequence of equations (\ref{PolarEQM}) and (\ref{Hatdefn}) that \begin{equation}\label{Geqn} \frac{dG}{dt} = r F_{\theta}. \end{equation} To determine $d\ell/dt$ and $dg/dt$, we first note the following useful relationships that are obtained from equations (\ref{Leqn}) and (\ref{Geqn}), the definitions of the Delaunay elements, and the Kepler law $\omega = 1/L^3$, \begin{eqnarray}\label{adotedot} \frac{da}{dt} & = & \frac{2}{\omega\sqrt{1-e^2}}\left[F_r \: e \sin{\hat{v}} + F_{\theta} \: \frac{a(1-e^2)}{r} \right], \nonumber \\ \frac{de}{dt} & = & \frac{\sqrt{1-e^2}}{a\omega}\left[F_r \: \sin{\hat{v}} + F_{\theta} \:\left(e+ \frac{r+a}{a}\cos{\hat{v}}\right) \right]. \end{eqnarray} Consider first $dg/dt = d\theta/dt - d\hat{v}/dt$, where \begin{equation}\label{thetadot} \frac{d\theta}{dt} = \frac{\sqrt{a(1-e^2)}}{r^2} \: . \end{equation} Differentiating the relationship \[\ln{r} = \ln{a} + \ln{(1-e^2)} - \ln{(1+e\cos{\hat{v}})}\] with respect to time and using the definition of $g$ together with equations~(\ref{adotedot}) and~(\ref{thetadot}), we have \begin{equation}\label{geqn} \frac{dg}{dt} = \frac{\sqrt{a(1-e^2)}}{e}\left[-F_r \: \cos{\hat{v}} + F_{\theta}\, \Big( 1 + \frac{r}{a(1-e^2)} \Big)\sin{\hat{v}} \right]. \end{equation} Let us now compute $d\ell/dt$ using the definition of $\ell$: The Kepler equation $\ell=\hat{u}-e\sin{\hat {u}}$ expresses the mean anomaly $\ell$ in terms of the eccentric anomaly and the eccentricity of the orbit. Thus, \begin{equation}\label{ell1eqn} \frac{d\ell}{dt} = \frac{r}{a}\: \frac{d\hat{u}}{dt}-\frac{r}{a\sqrt{1-e^2}}\sin{\hat{v}} \: \frac{de}{dt}. \end{equation} Using the identity \begin{displaymath} \dot{r} - \dot{a}\frac{r}{a}+\dot{e}\: \frac{a-r}{e} = \frac{er\sin{\hat{v}}}{\sqrt{1-e^2}}\: \frac{d\hat{u}}{dt}, \end{displaymath} which is obtained by differentiating $r=a(1-e\cos{\hat{u}})$ with respect to time, and equation (\ref{adotedot}), we find \begin{equation}\label{elleqn} \frac{d\ell}{dt} = \omega + \frac{r}{e\sqrt{a}}\left[F_r\: (-2e+\cos{\hat{v}}+e\cos^2{\hat{v}})-F_{\theta}\, (2+e\cos{\hat{v}})\sin{\hat{v}} \right]. \end{equation} Therefore, the equations of motion in Cartesian coordinates are equivalent to Delaunay's system given by equations~(\ref{Leqn}),~(\ref{Geqn}),~(\ref{elleqn}), and~(\ref{geqn}). These equations for the planar problem are closely related to corresponding equations for the variation of the orbital elements of the osculating ellipse. In the general, i.e. three-dimensional, Kepler problem, the osculating ellipse is characterized by six orbital elements and the variation of these orbital elements due to perturbing forces is given by a set of equations that are usually associated in celestial mechanics with the name of Lagrange (e.g. ``Lagrange's planetary equations''). In general, the path of the osculating ellipse in the six-dimensional manifold of orbital elements encounters singularities at $e=0$ and $e=1$. These singularities, corresponding to a circle ($G=L$) and a straight line ($G=0)$, respectively, are also evident in the Delaunay equations. The right hand side of the set of Delaunay equations should also be represented in Delaunay elements. To accomplish this, recall equations~(\ref{Hatdefn}) and let $\delta f_D-\epsilon {\cal K}_D$, where $D$ denotes any one of the Delaunay variables, be the net contribution of radiation reaction and external forces to the rate of variation of Delaunay element $D$ with respect to time. We have shown in a previous paper \cite{cmr} that \begin{eqnarray} {\cal K}_L & = & \frac{\partial{\cal C}}{\partial \ell} \phi(t) + \frac{\partial{\cal S}}{\partial \ell} \psi(t), \nonumber \\ {\cal K}_{G} & = & \frac{\partial{\cal C}}{\partial g} \phi(t) + \frac{\partial{\cal S}}{\partial g} \psi(t), \nonumber \\ {\cal K}_{\ell} & = &- \left(\frac{\partial{\cal C}}{\partial L} \phi(t) + \frac{\partial{\cal S}}{\partial L} \psi(t)\right), \nonumber \\ {\cal K}_g & = & -\left(\frac{\partial{\cal C}}{\partial G} \phi(t) + \frac{\partial{\cal S}}{\partial G} \psi(t)\right), \end{eqnarray} where \begin{equation}\label{phipsieqn} \phi(t) = \frac{1}{2}\alpha\Omega^2\cos{(\Omega t)}, \: \: \psi(t) = \frac{1}{2}\beta\Omega^2\cos{(\Omega t+\rho)}, \end{equation} and \begin{eqnarray}\label{CSfs} {\cal C}(L,G,\ell,g) & = & \frac{5}{2}a^2e^2\cos{2g} \nonumber \\ && +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} \nonumber \\ && +a^2\sum_{\nu =1}^\infty (A_\nu\sin{2g}\cos {\nu\ell} +B_\nu\cos{2g}\sin{\nu \ell}). \end{eqnarray} Here we have introduced $A_{\nu}$ and $B_{\nu}$, which are given in terms of the Bessel functions by \begin{eqnarray}\label{ABnu} A_\nu & = &\frac{4}{\nu^2e^2} (2\nu e(1-e^2)J_\nu'(\nu e)-(2-e^2)J_\nu(\nu e)), \nonumber \\ B_\nu & = & -\frac{8}{\nu^2e^2}\sqrt{1-e^2}\, (e J_\nu'(\nu e)-\nu (1-e^2)J_\nu(\nu e)). \nonumber \end{eqnarray} Let us now compute the contribution of radiation damping terms to the equations of motion. Using equation (\ref{frftheta}) and the relationships \[v^2 = \frac{2}{r} -\frac{1}{a}, \: \: \dot{r}^2 = -\frac{1}{a}+\frac{2}{r} -\frac{G^2}{r^2}\:, \] we find that \begin{eqnarray}\label{frfthetadel} f_r & = & -\frac{4 e\sin{\hat{v}}}{GL^2r^3} \left(1- \frac{10}{3}\frac{L^2}{r}-5\frac{L^2G^2}{r^2}\right), \nonumber \\ f_{\theta} & = & -\frac{18 G}{L^2 r^4}\left(1-\frac{20}{9}\frac{L^2}{r} +\frac{5}{3}\frac{L^2G^2}{r^2}\right). \end{eqnarray} The substitution of equation~(\ref{frfthetadel}) into the Delaunay equations via equation~(\ref{Hatdefn}) results in the set of equations presented in equation~(\ref{feq1}), where $e = \sqrt{1-G^2/L^2}$; moreover, these relations involve powers of $1/r$ as well as $\sin{\hat{v}}$ times powers of $1/r$. Therefore, to complete the discussion we need Fourier expansions of these expressions in terms of $\cos{\nu\ell}$ and $\sin{\nu\ell}$ similar to those given in equation (\ref{CSfs}). This can be accomplished by extending the classical methods of celestial mechanics described, for instance, in the definitive monograph of Watson~\cite{watson}. The final expressions for $f_D$, $D\in\{L,G,\ell,g\}$, are not used explicitly in this paper; therefore, we do not present them here for the sake of brevity. Finally, let us note that the Delaunay equations can now be expressed in the form given in equation~(\ref{D2EQM1}). \section{Average Rate of Damping}\label{RadForceApp} The average contribution of radiation damping to the equations of motion in Delaunay variables is employed in this paper in the calculation of the bifurcation function (Section~\ref{conpo}) as well as in the description of averaged dynamics (Section~\ref{avedyn}). In the light of the results of the previous appendix, the quantity $f_D$, where $D$ is any one of the Delaunay variables, is given by equation~(\ref{D2EQM1}) and can be expressed as a Fourier series \begin{equation}\label{AveBasicEq} f_{D}(L, G, \ell, g) = a_0 +\sum_{\nu = 1}^{\infty}(a_{\nu}\cos{\nu\ell}+b_{\nu}\sin{\nu\ell}). \end{equation} We are interested in $\bar f_D$, which is the average of $f_D$ with respect to the ``fast'' angular variable $\ell$, i.e. \[ \bar f_D:=\frac{1}{2\pi}\int_0^{2\pi}f_D(L,G,\ell,g)\,d\ell. \] It is clear that $\bar f_D=a_0$, and this is the quantity that is needed in Section~\ref{avedyn}. The contribution of radiation damping to the bifurcation function $\cal B$ is given by $\Delta{\cal F}_{D}$, where \begin{equation}\label{BifAveFun} {\cal F}_{D} = \int_0^{2\pi m/ \Omega} f_{D}(L, G, \omega t + \ell, g) \: dt. \end{equation} Using the change of variable $\tilde t = \Omega t/m + \ell/n$, the resonance condition $\omega = n\Omega /m$, and the periodicity of $f_{D}$ with respect to $\ell$, equation (\ref{BifAveFun}) can be written as ${\cal F}_{D}=2\pi n L^3 \bar{f}_D$. It remains to compute $\bar{f}_D$. Using the relation $d\ell = r^2 d\hat{v}/(a^2\sqrt{1-e^2})$, we can write \begin{equation}\label{87} \bar{f}_{\cal D} = \frac{1}{2\pi a^2\sqrt{1-e^2}} \int_0^{2\pi} r^2{\cal F}_{D}\,d\hat{v}= \frac{1}{L^3 G}, \end{equation} where the angular brackets denote the average over the true anomaly $\hat{v}$. The various averages may be evaluated using \[ G^2 \left< \frac{1}{r} \right> = 1,\quad G^4 \left< \frac{1}{r^2} \right> = 1+\frac{1}{2} e^2,\quad G^6 \left< \frac{1}{r^3} \right> = 1+\frac{3}{2} e^2, \] \[ G^8 \left< \frac{1}{r^4} \right> = 1+3 e^2 + \frac{3}{8} e^4, \quad G^{10} \left< \frac{1}{r^5} \right> = 1+ 5 e^2 + \frac{15}{8} e^4, \] \[ G^2 \left< \frac{\sin^2{\hat{v}}}{r} \right> = \frac{1}{2}, \quad G^4 \left< \frac{\sin^2{\hat{v}}}{r^2} \right> = \frac{1}{2}\left( 1+\frac{e^2}{4} \right), \] \begin{equation}\label{AveValues} G^6 \left< \frac{\sin^2{\hat{v}}}{r^3} \right> = \frac{1}{2}\left( 1+\frac{3}{4}e^2 \right). \end{equation} We find that \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} where the last two equations simply follow from the fact that \[=0\] for every continuous function $h$. For the bifurcation function $\cal B$, only ${\cal F}_L$, ${\cal F}_G$, and ${\cal F}_g$ are needed and these are then given by \begin{eqnarray}\label{fqAve2} {\cal F}_L& =& - 2\pi n\frac{L^3}{G^7}\left( 8+\frac{73}{3} e^2+ \frac{37}{12} e^4 \right), \nonumber \\ {\cal F}_G &=& -2\pi n\frac{1}{G^4} (8+7e^2), \nonumber \\ {\cal F}_g & = & 0. \end{eqnarray} %\begin{thebibliography}{99} \begin{references} \bibitem{hulse} Hulse R A and Taylor J H 1975 Discovery of a pulsar in a binary system {\em Astrophys.\ J.} {\bf 195} L51--L53 \bibitem{taylor} Taylor J H, Wolszczan A, Damour T and Weisberg J M 1992 Experimental constraints on strong-field relativistic gravity {\em Nature} {\bf 355} 132--136 \bibitem{cmr} Chicone C, Mashhoon B and Retzloff D G 1996 Gravitational ionization: periodic orbits of binary systems perturbed by gravitational radiation {\em Annales de l'Institut Henri Poincar\'{e}, Physique Th\'{e}orique} {\bf 64} 87--125 \bibitem{cmr1} Chicone C, Mashhoon B and Retzloff D G 1996 On the ionization of a Keplerian binary system by periodic gravitational radiation {\em J. Math.\ Phys.} {\bf 37} 3997--4016 \bibitem{peters} Peters P C and Mathews J 1963 Gravitational radiation from point masses in a Keplerian orbit {\em Phys. Rev.} {\bf 131} 435--440 \bibitem{mashoon1} Mashhoon B 1978 On tidal resonance {\em Astrophys.\ J.} {\bf 223} 285--298; 1977 Tidal radiation {\em ibid.} {\bf 216} 591--609 \bibitem{mashoon2} Mashhoon B 1979 On the detection of gravitational radiation by the Doppler tracking of spacecraft {\em Astrophys.\ J.} {\bf 227} 1019--1036; 1980 Absorption of gravitational radiation {\em GRG Abstracts} (Jena Meeting) 1505--1507 \bibitem{mashoon3} Mashhoon B, Carr B J and Hu B L 1981 The influence of cosmological gravitational waves on a Newtonian binary system {\em Astrophys.\ J.} {\bf 246} 569--591 \bibitem{grishchuk} Grishchuk L P and Polnarev A G 1980 Gravitational waves and their interaction with matter and fields {\em General Relativity and Gravitation (Vol 2)} ed A Held (New York: Plenum Press) 393--434 \bibitem{Linet} Linet B 1982 Absorption of gravitational waves by nearly Newtonian systems {\em Gen.\ Rel.\ Grav.} {\bf 14} 479--492 \bibitem{Chau} Nelson L A and Chau W J 1982 Orbital perturbations of a gravitationally bound two-body system with the passage of gravitational waves {\em Astrophys.\ J.} {\bf 254} 735--747 \bibitem{Petukhov} Gal'tsov D V and Petukhov V I 1982 Amplification of gravitational waves interacting with a binary system {\em Phys.\ Lett.\ A} {\bf 90} 351--353 \bibitem{Ivashchenko} Ivashchenko A V 1987 Variation of the Keplerian elements of a planetary orbit under the action of a gravitational wave {\em Sov.\ Astron.} {\bf 31} 76--79 \bibitem{Sbytov} Kochkin R A and Sbytov Yu G 1987 Influence of a periodic gravitational wave on the parameters of a binary system {\em Sov.\ Astron.} {\bf 31} 537--540 \bibitem{Matsuda} Futamase T and Matsuda T 1979 Resonance between primordial gravitational waves and gravitationally bound systems {\em Prog.\ Theor.\ Phys.} {\bf 61} 86--93; 1979 Resonance between primordial gravitational waves and a galaxy {\em ibid}. {\bf 61} 1320--1326 \bibitem{poincare} Poincar\'{e} H 1892-99 {\em Les M\'{e}thodes Nouvelles de la M\'{e}canique C\'{e}leste (Vol 1-3)} (Paris: Gauthier-Villars) \bibitem{ww} Walker M and Will C M 1979 Relativistic Kepler problem. I. Behavior in the distant past of orbits with gravitational radiation damping {\em Phys. Rev. D} {\bf 19} 3483--3494; 1979 Relativistic Kepler problem. II. Asymptotic behavior of the field in the infinite past {\it ibid}. {\bf 19} 3495--3508 \bibitem{lorentz} Lorentz H A 1952 {\it The Theory of Electrons} (New York: Dover) \bibitem{Kov} Kovalevsky J 1967 {\em Introduction to Celestial Mechanics,} {\em Astrophysics and Space Science Library (Vol~ 7)} (New York: Springer-Verlag) \bibitem{Stern} Sternberg S 1969 {\em Celestial Mechanics (Vol 1--2)} (New York: Benjamin) \bibitem{ccc} Chicone C 1995 A geometric approach to regular perturbation theory with an application to hydrodynamics {\em Trans.\ Amer.\ Math.\ Soc.} {\bf 347} 4559--4598 \bibitem{arnold} Arnold V I 1963 Small denominators and problems of stability of motion in classical and celestial mechanics {\em Russian Math.\ Surveys} {\bf 18} 85--191 \bibitem{liberman} Lichtenberg A J and Lieberman M A 1992 {\em Regular and Chaotic Dynamics} (New York: Springer-Verlag) \bibitem{arnold2} Arnold V I (ed) 1988 {\em Dynamical Systems III}, {\em Encyclopedia of Mathematical Sciences (Vol 3)} (New York: Springer-Verlag) \bibitem{lochak} Lochak P and Meunier C 1988 {\em Multiphase Averaging for Classical Systems with Applications to Adiabatic Theorems} (New York: Springer Verlag) \bibitem{mashhoon3} Mashhoon B 1993 Accelerated observers and conformal invariance {\em Quantum Gravity and Beyond: Essays in Honor of Louis Witten}, ed F Mansouri and J J Scanio (Singapore: World Scientific) 257--271 \bibitem{Stiefel} Kustaanheimo P and Stiefel E 1965 Perturbation theory of Kepler motion based on spinor regularization {\em J. Reine Angew.\ Math.} {\bf 218} 204--219; Stiefel E L and Scheisele G 1971 {\em Linear and Regular Celestial Mechanics} (New York: Springer-Verlag) \bibitem{wig} Wiggins S 1990 {\em Introduction to Applied Nonlinear Dynamical Systems and Chaos} (New York: Springer-Verlag) \bibitem{watson} Watson G N 1966 {\em A Treatise on the Theory of Bessel Functions} (New York: Cambridge University Press) %\end{thebibliography} \end{references} %\newpage %\centerline{Figure Captions} %\vspace{2in} \centerline{Hard copy or PostScript files for figures are available from} \begin{verbatim} carmen@chicone.math.missouri.edu \end{verbatim} \begin{figure}[h] \vspace*{1in} \caption[]{ The plots are for system~(\ref{NumExpEq}) with parameter values $\epsilon=10^{-4}$, $\delta/\epsilon=10^{-3}$, $\alpha=2$, $\beta=2$, $\rho=0$, and $\Omega= .6495190528$. The top panel shows $L=\sqrt{a}$ versus time for the initial conditions $(p_r,p_\theta,r,\theta)$ equal to $(.38927, .61085 , 2.70900, -2.45677)$. The middle panel shows $L$ versus time and the bottom panel shows $G$ versus $L$ for the initial conditions $(.3893, .6109 , 2.709, -2.4568)$. Here, $L\approx 1.1547$ corresponds to $(1:1)$ resonance, $1/L^3=\Omega$. \label{capfig}} \end{figure} \begin{figure}[h] \caption[]{ The plot of $G$ versus $L$ for system~(\ref{NumExpEq}) with initial conditions $(p_r,p_\theta,r,\theta)$ equal to $(.5, 1 , 1, 0)$ and parameter values $\epsilon=10^{-3}$, $\delta/\epsilon=10^{-5}$, $\alpha=2$, $\beta=2$, $\rho=0$, and $\Omega=1.299038106$. The top panel depicts 157000 points, one after each time interval $2 \pi/\Omega$. The middle panel is a blow up for iterates 110000--130000, the bottom panel for 140000--155000. \label{trans}} \end{figure} \end{document}