Instabilities in the Sun-Jupiter-Asteroid Three Body Problem

Transcription

1 Noname manuscript No. (will be inserted by the editor) Instabilities in the Sun-Jupiter-Asteroid Three Body Problem John C. Urschel Joseph R. Galante Received: date / Accepted: date Abstract We consider dynamics of a Sun-Jupiter-Asteroid system, and, under some simplifying assumptions, show the existence of instabilities in the motions of an asteroid. In particular, we show that an asteroid whose initial orbit is far from that of the orbit of Mars can be gradually perturbed into one that crosses Mars orbit. Properly formulated, the motions of the asteroid can be described as a Hamiltonian system with two degrees of freedom, with the dynamics restricted to a large open region of the phase space reduced to an exact area preserving map. Instabilities arise in regions where the map has no invariant curves. The method of MacKay and Percival is used to explicitly rule out the existence of these curves, and results of Mather abstractly guarantee the existence of diffusing orbits. We emphasize that finding such diffusing orbits numerically is quite difficult, and is outside the scope of this paper. Keywords Hamiltonian Systems Restricted Problems Instabilities Asteroid Belt 1 Introduction We consider the restricted circular planar three body problem (RCP3BP) with two massive primaries, which we call the Sun and Jupiter, that perform uniform John C. Urschel Department of Mathematics, Penn State University University Park, State College, PA Tel.: Fax: Joseph Galante Unaffiliated, formerly Department of Mathematics, University of Maryland

2 2 John C. Urschel, Joseph R. Galante circular motion about their center of mass. The system is normalized to mass one so the Sun has mass 1 µ and Jupiter mass µ. We further normalize so that Jupiter rotates with unit speed, resulting in a period of 2π for the primaries. The distance from the Sun to Jupiter is constant and also normalized to one. Our goal is to understand the behavior of the massless asteroid, whose position in polar coordinates is denoted by (r, ψ). It is convenient to consider the system in a rotating frame of reference which rotates with unit speed in the same direction as Jupiter. In this system, the Sun and Jupiter are fixed points on the x-axis corresponding to ψ = 0. We let (r, ϕ) = (r, ψ t) denote the motion of the asteroid in this rotating frame of reference. Our system has a Hamiltonian of the form: H P olar = H 2BP (SA) + H(r, ϕ) := P r P ϕ 2 2r 2 P ϕ 1 + H(r, ϕ; µ) r where P r and P ϕ are the momenta variables conjugate to r and ϕ respectively (Arnol d et al. 2006) and H is the µ-small perturbation of the associated Sun- Asteroid two body problem (2BP(SA)) by the presence of Jupiter. This system arises by initially considering the planar 3BP where the asteroid has mass m, and letting m 0. From this, we have the following equations of motion to describe the motion of our asteroid: ϕ = H P olar P ϕ = P ϕ r 2 1 P ϕ = H P olar ϕ ṙ = H P olar = P r P r = H P olar P r r = P 2 ϕ r 3 = H ϕ 1 r 2 H r In addition, the RCP3BP has a conserved quantity known as the Jacobi constant. J(r, ϕ, ṙ, ϕ) = r2 2 + µ + 1 µ ṙ2 + r 2 ϕ 2 =: U(r, ϕ) ṙ2 + r 2 ϕ 2 d J d S 2 2 where d J and d S are distances from the asteroid to Jupiter and the Sun, respectively. d J (r, ϕ) = ( r 2 2(1 µ)r cos(ϕ) + (1 µ) 2) 1 2 d S (r, ϕ) = ( r 2 + 2µr cos(ϕ) + µ 2) 1 2 (1) The Jacobi constant can be thought of as the total energy of our massless asteroid, with respect to our rotating frame. For a derivation of this conserved quantity, refer to sections and of the text by Arnol d, Kozlov, and Neishtadt (2006). Denote by H(J 0 ) = {(r, ϕ) : U J 0 } a set of points on the plane of motion (configuration space). The connected components of this set are called the Hill regions associated with the Jacobi constant J 0. These regions are the locations in the (r, ϕ) plane (shaded regions in figure 6)

3 Instabilities in the Sun-Jupiter-Asteroid Three Body Problem 3 Fig. 1 Disjoint Hill Regions for µ = 10 3 and J (Galante et at. 2011) Fig. 2 Eccentricity vs. Aphelion and Perihelion Radii for J 0 = 1.55 where the asteroid is allowed to move. Fixing the Jacobi constant restricts dynamics to an invariant energy surface, denoted S(J 0 ) := {(r, ϕ, ṙ, ϕ) : J(r, ϕ, ṙ, ϕ) = J 0 } Most of these surfaces are smooth 3-dimensional manifolds. Let us denote by RCP 3BP (µ, J 0 ) the RCP3BP with Sun-Jupiter mass ratio µ and dynamics restricted to the surface S(J 0 ). For µ 10 3 and J , the set H(J 0 ) consists of three disjoint connected components: a region around the Sun called the inner Hill region, a region around Jupiter called the lunar Hill region, and a noncompact region called the outer Hill region. The boundary of these regions can be found by considering the zero velocity curves, given by ṙ 2 + r 2 ϕ 2 = 0, which are on the boundary of the Hill regions (Arnol d et al. 2006). In this paper, we consider only orbits contained in the inner Hill region, denoted by H in (J 0 ). For convenience, denote S in (J 0 ) = H in (J 0 ) S(J 0 ). When dynamics S in (J 0 ) is considered, we refer exclusively to the case when the inner Hill region is disjoint from the other two. For J 0 greater than 1.52, asteroids stay uniformly bounded away from Jupiter for all time by the energy surface constraint. However, for high eccentricities, the asteroid can make near collisions with the Sun (referred to as sun-grazers; see figure 2). A general result states that there are KAM tori near the Sun which prevent

4 4 John C. Urschel, Joseph R. Galante collision (Chenciner et al. 1988). For small µ and away from collisions, the RCP3BP is nearly integrable and can be approximated with the Sun-Asteroid two body problem (2BP(SA)), corresponding to µ = 0. Elliptic motions of a 2BP have two special points where the radial velocity ṙ of the asteroid is zero. The perihelion is the closest point to the Sun 1, denoted r perih, and the aphelion is the farthest point from the Sun, denoted r aph. We define the osculating (or instantaneous) eccentricity e(t) for the RCP3BP to be the eccentricity of the asteroid in the unperturbed 2BP(SA) system, with initial conditions taken to be those of asteroid in the RCP3BP at time t. Theorem 1 Consider the restricted circular planar three body problem with Sun- Jupiter mass ratio µ = Fix a Jacobi Constant J 0 = 1.55, so that there are three disjoint Hill regions and consider dynamics in the inner Hill region. Moreover, assume that ϕ(t ) P ϕ0 > 0 for all trajectories in a set Ω twist {e [0.09, 0.8]} which start on the surface Σ = {ṙ = 0, r 0, ϕ > 0} and return to Σ after time T. Then there exist constants e and e + where e 0.2 and e + 0.6, and trajectories of an asteroid with initial eccentricity e that increase to eccentricity e +. Remark: In section 2.1 we state that under these conditions, there exist orbits that become Mars crossing. We build a mathematical framework and obtain a sufficient condition for Theorem 1 to hold. This condition guarantees applicability of so-called Aubry-Mather theory. Then we run numerical tests and evaluate the region of applicability of this theory in terms of eccentricity. Convincing numerical data shows for the Jacobi constant J 0 = 1.55 there existences of orbits whose eccentricity is varying from 0.2 to 0.6. For applicability of Aubry-Mather theory for the outer Hill region of the RCP3BP, refer to the works of Galante and Kaloshin (2011). For a treatment of diffusion in the Elliptic Restricted 3BP, refer to Liao and Saari (1998). In addition, for a treatment of Mars crossing orbits for the elliptic case, refer to the works of Wisdom (1982, 1983, 1985). 2 Numerical Aspects, Aubry-Mather Theory, and Regions of Instability Let us begin to describe the mathematical framework we employ by first noting that a Jacobi energy surface is a 3-dimensional manifold. Fix µ 10 3 and Jacobi constant J The section Σ = {ṙ = 0, r 0, ϕ > 0} is a well-defined Poincaré section in the inner Hill region. This leads to a well-defined (Poincaré) 1 To be precise, the perihelion is the point where the asteroid is at the closest point to the center of mass of the system, and the Sun is within µ of the center of mass. However, in our Solar System, the radius of the Sun is approximately the Sun-Jupiter distance, so we allow this slight abuse in terminology for small µ.

5 Instabilities in the Sun-Jupiter-Asteroid Three Body Problem 5 return map on an open set Ω. One can show that Ω is homeomorphic to a 2- dimensional cylinder and can be parametrized by angle ϕ and conjugate momenta P ϕ, or, alternatively, by ϕ and eccentricity e = e(p ϕ, J 0 )). Suppose T 2 S in (J 0 ) is an invariant set of the RCP3BP that is diffeomorphic to a 2-dimensional torus. Call T 2 rotational if it cannot be continuously deformed inside S in (J 0 ) into a closed curve or a single point. When µ = 0 (i.e., when there is no perturbation), the problem reduces to the 2BP(SA) system and every such rotational 2-torus is defined by {e = e 0 0}. Bounded motions correspond to e 0 [0, 1). In general, for e bounded away from 1 and µ sufficiently small, many of these rotational 2-tori survive due to KAM (Siegel et al. 1971). Celletti and Chierchia gave a computer assisted proof using µ 10 3 and J in the inner Hill region to show that near e = 0.3 there is a rotational 2-torus T 2 separating S in (J 0 ) into a compact Below T 2 component and a noncompact Above T 2 component (Celletti et al. 2007). We present a complementary method for a specific value of J 0 = 1.55; however, the method works for any µ 10 3 and J Define a Region of Instability (RI) as an open invariant set in S in (J 0 ) which is homeomorphic to an annulus, and has no rotational 2-dimensional tori inside. If there is a rotational 2-torus, then it separates S in (J 0 ) into above and below regions. This provides a topological obstruction to instability. To construct regions of instability, one must know about the existence of invariant curves in a given region. Theorem 2 In the setting of Theorem 1, the RCP3BP, restricted to the inner Hill region of a Jacobi energy surface, has a well-defined Poincaré map F : Ω( Σ) Σ. Its restriction to Ω twist is an exact area-preserving twist map, and there is a subregion Ω AM Ω twist with [e, e + ] T Ω AM, such that for any rotation number ω [ω, ω + ] (ω ± = ω ± (e ± )) there exists a corresponding Aubry-Mather set Σ ω Ω twist with Σ ω Ω AM. Moreover, none of the Aubry-Mather sets Σ ω with ω [ω, ω + ] are ever invariant curves. Abstractly, it is not clear that e < e + (or that ω < ω + ). We fix µ = 10 3 and J 0 = 1.55 for the sake of working with concrete numbers. It is shown numerically for these parameters that on the subset T {e [0.09, 0.8]} Ω twist Ω, the Poincaré map, denoted F = F µ,j0, is an exact, area-preserving twist (EAPT) map. EAPT maps of this region can be studied using Aubry-Mather theory. However it is important to emphasize that Aubry-Mather (AM) theory does not apply directly as Ω may not be an invariant set. Only through a refined study of properties of F can it be shown that AM theory can be applied. A brief review of AM theory is found in section 6 for unfamiliar readers. The cone crossing condition introduced by MacKay and Percival (1985) is used to obtain a sufficient condition to rule out invariant curves, as well as establish the existence of a range of Aubry-Mather invariant sets. In particular, it is used to 2 For J 0 near or less than 1.52 collisions with Jupiter are hard to exclude.

6 6 John C. Urschel, Joseph R. Galante establish that there are e ±, e < e + with T [e, e + ] Ω AM. It is the set Ω AM that AM theory is ultimately applied to in order to obtain the existence of orbits making large deviations in eccentricity. The cone crossing condition exploits the fact that Aubry-Mather sets are Lipschitz graphs. The MacKay-Percival sufficient condition is obtained in two steps: Establish bounds on the Lipschitz slope of potential invariant curves as well as Aubry-Mather sets. This gives rise to a conefield in the tangent space to Ω twist, and a range of rotation numbers such that the corresponding Aubry-Mather sets are contained in Ω AM. Show that there is a vertical strip of initial conditions such that there are tangent vectors crossing this conefield. This rules out invariant curves in Ω twist. Once a region free of invariant curves is established inside the region of twist, an application of Mather s variational method provides the existence of an orbit crossing Ω AM, completing the proof of Theorem 1 (Mather et al. 1994). It is important to emphasize that the orbits in Theorem 1 are not constructed by means of numerical integration, but through abstract variational principles found in Aubry- Mather theory. Convincing numerical data is presented to obtain concrete bounds. To make the proof fully rigorous, our initial numerical tests must be verified, for example by using the machinery of interval arithmetic. 2.1 Mars crossing orbits On S in (1.55), an asteroid with e = 0.2 has radius r [0.533, 0.800]. The semimajor axis of the orbit of Mars is 1.5AU and the semi-major axis of the orbit of Jupiter is 5.2AU which places Mars at position r = in normalized coordinates. Also on S in (1.55), an asteroid with e = 0.6 has radius r [0.204, 0.818]. Thus we have the following: Corollary 3 Under the hypothesis of Theorem 1, there exist orbits of the asteroid which become Mars crossing. In the Asteroid belt there are approximately 1.7 million asteroids of radius of at least 1km. The orbits with J 0 = 1.55 and e = 0.2 are at the boundary of the main Asteroid belt. Some of them might have had significant oscillations of eccentricity as this theorem suggests. Even though an asteroid whose orbit is Mars crossing has a small chance of being captured by it (see section of Arnol d, Kozlov, and Neishtadt (2006)), enough attempts could have led to a capture. However, we stress that Mars is a not a part of our model, and, therefore, this claim cannot be completely justified by the numerics in this paper. For larger eccentricity asteroids, the perturbative effects of Mars become less negligible. Letting Jupiter have positive eccentricity leads to the so-called Restricted Planar Elliptic Three Body Problem. This is a Hamiltonian system of two and a half degrees of freedom. It has more room for instabilities. The well-known instabilities occur in Kirkwood gaps, when the period of Jupiter and of an asteroid are

7 Instabilities in the Sun-Jupiter-Asteroid Three Body Problem 7 given by a rational with both numerators and denominators small, e.g. 1/3, 2/5. Mathematically, these instabilities have been recently studied by Fejoz, Guardia, Kaloshin, and Roldan (2011). For a review of diffusion in asteroid belt resonances, refer to the work of Ferraz-Mello (1999). 3 Plan of the Proof The overall plan of our proof is as follows: we convert the dynamics to a return map, eliminate the existence of invariant curves over a section of the return map, and, finally, use existing theory to conclude that this non-existence of invariant curves results in an orbit which diffuses over the interval [e, e + ]. These three steps combined form the numerical proof of the main theorem. More explicitly, the general flow of the remainder of the paper is as follows: Step 1: The dynamics is formulated as that of an exact area preserving twist map (EAPT) F µ. In the case µ = 0, it is easy to show that F µ is a twist map. For µ > 0, numerics are used to confirm this for a large region of the phase space. Step 2: Approximate the EAPT F µ by G. The approximated map G shares the same qualitative behaviors as F µ, has some desired properties that F µ lacks, and is much easier to deal with numerically. However, G is not guaranteed to be area preserving. Step 3: Rule out the the presence of invariant curves for G. We use the aforementioned cone crossing condition. By design, if this condition holds for an area-preserving twist map F µ in a region Ω, then F µ has no rotational invariant curves in Ω. This condition relies on a certain conefield. We construct a conefield for the map G and verify that the cone crossing condition holds for G in a certain region Ω. Then we show that this condition is robust and also holds for F µ. This implies that F µ has no rotational invariant curves in Ω. Step 4: Apply the Mather Connecting Theorem We utilize the Mather Connecting Theorem and Aubry-Mather theory to conclude there exists an orbit which diffuses in the manner we prescribed. 4 Formulation of a Twist Map Recall that motions of the asteroid in rotating polar coordinates (r, ϕ) can be viewed as the solutions to Hamilton s equations of motion with a Hamiltonian of

10 10 John C. Urschel, Joseph R. Galante eccentric anomaly, is given implicitly by the Kepler equation u e sin(u) = l. The variable ( ) ϕ = g + f, where the variable f, known as the true anomaly, is given by tan f 1+e 2 = 1 e tan ( ) u 2 or, alternatively, r = L 2 (1 e 2 ) 1+ecos(f). Because the return map F µ is defined using trajectories corresponding to a full revolution of the asteroid, the variable u is periodic with period 2π. When u = 2πk, k Z, then f = 2πk. When µ = 0, there is no precession of the ellipse of orbit, i.e. g + t is constant in rotating coordinates, and hence ϕ + t = ψ in nonrotating coordinate is periodic. (This a restatement of the fact that the 2BP(SA) has a periodic solution, an ellipse.) Due to the 2π periodicity of the equations of motion in ϕ, ϕ 1 simplifies to ϕ 1 = ϕ 0 + T (ϕ 0, P ϕ0 ). In the 2BP(SA) case (µ = 0), it holds that T (ϕ 0, P ϕ0 ) = 2π(2J 0 2P ϕ0 ) 3/2. Hence, it is possible to write in general F µ : ( ϕ0 P ϕ0 ) ( ϕ1 P ϕ1 ) = ( ) ϕ0 + 2π(2J 0 2P ϕ0 ) 3/2 + µf ϕ (ϕ 0, P ϕ0 ; µ) P ϕ0 + µf Pϕ (ϕ 0, P ϕ0 ; µ) (3) where F ϕ (.; µ), F Pϕ (.; µ) are smooth functions in a domain bounded away from e = 1. We note that this map is similar in many ways to the Keplerian map describing highly eccentric cometary motion as described by Broucke and Petrovsky (1987). 4.2 F µ is an EAPT Because F µ arises as a return map of a Hamiltonian system, the map is area preserving. Additionally, due to the fact that the Hamiltonian has two degrees of freedom, the map is also exact. We claim the map is also a twist map. (See section 6 for precise definitions of exact, area preserving, and twist for an abstract map.) We say that F µ is a twist map in a region Ω if ϕ 1 P ϕ0 > 0 (ϕ 0, P ϕ0 ) Ω It is easy to see that this holds everywhere for µ = 0. Lemma 4 The unperturbed map F 0 is a twist map in the inner Hill Region for J 0 > 1.5. ϕ Proof: 1 P ϕ0 = 6π(2J 0 2P ϕ0 ) 5/2 > 0 since in the inner Hill region P ϕ0 < J 0 when J 0 > 1.5. For µ 0, careful estimates of the perturbation terms are needed to prove twisting over a domain. 3 Let Ω twist be the largest domain in Ω where the inequality ϕ 1 P ϕ0 > 0 holds. In practice, we are only able to find a subdomain of Ω twist 3 Notice that the angle ϕ enters into the perturbation H (see (2)). As the Poincaré map F µ is defined for approximately one revolution of the asteroid, then the change in the ϕ component should average out for higher order terms in the µ expansion of H.

11 Instabilities in the Sun-Jupiter-Asteroid Three Body Problem 11 which is selected experimentally to be the largest domain in which the twist inequality can be proven to hold using the current best estimates on the perturbation terms. The twist property is expected to fail near e = 1 because of close encounters with the Sun. The map F µ is not even defined near this region because of this. Current best estimates show the following: Claim 5 For µ = 10 3 and J 0 = 1.55, the map F µ is a twist map in the inner Hill Region for e [0.09, 0.8]. Numerical Proof: A computer is used to compute the map F µ. The annulus (ϕ, P ϕ ) [0, 2π] [0.39, 0.83] is divided up into boxes of size π 16 by The region e [0.09, 0.8] corresponds to P ϕ [0.39, 0.83], using conservation of the Jacobi constant 4. To compute the twist term, the partial derivative ϕ 1 P ϕ0 is estimated using a difference quotient with initial conditions varying by 10 7 in the P ϕ direction. The computer finds the required inequality holds in the specified region at each grid point using this approximation. The dominant factor in the computation of the partial derivatives for twist arises from the 2BP(SA) system (see proof of Lemma 4). Because we are away from singularities at e = 1 (which can cause near collisions with the sun), the approximation made is acceptable and the behavior is close to that of the unperturbed system. Remark: Currently this proof is verified using only numerics. To make this numerical evidence into a proof, one can use rigorous numerical integration to compute the partial derivative in question, and use interval arithmetic to make a uniform bound for each box in the annulus. This has been carried out for a few boxes in the grid using the CAPD package. Implementing rigorous numerical integration for cases when the asteroid makes a close approach to the sun requires extraordinarily good estimates on the perturbation term, meaning a high order integrator with small step size, and a very fine grid is needed. This can be quite costly in terms of computing power, and the authors were unable to show the result for every box in the grid covering the given annulus. 5 The Cone Crossing Condition First we consider an abstract setting and then apply it to the RCP3BP. This section follows the work found in MacKay and Percival s paper (1985). Suppose F (θ 0, I 0 ) = (θ 1, I 1 ) is an EAPT on an invariant domain Ω diffeomorphic to a cylinder T R (θ, I). A (rotational) invariant curve C Ω (in this paper we study only rotational invariant curves, so we omit rotational for brevity) is an invariant set of F which is diffeomorphic to a circle and cannot be continuously deformed to a single point in Ω. A theorem of Birkhoff states that invariant curves are graphs over T. 4 More precisely, it follows from the geometry of ellipses that e = 1 G2 L 2 in Delauney coordinates; a quick conversion to polar coordinates yields the formula e = 1 2J 0 Pϕ 2 + 2Pϕ. 3

12 12 John C. Urschel, Joseph R. Galante Theorem 6 (Birkhoff) Suppose C is an invariant curve of an EAPT. Then there exists a Lipschitz continuous function P such that C {(θ, P (θ)) : θ T}. Let D(x 1, x 2 ) = P (x 2) P (x 1 ) x 2 x 1, x 1 x 2. Recall that P is Lipschitz if there exist D ± such that for all x 1, x 2, D D(x 1, x 2 ) D +. Obtaining the Lipschitz property in Birkhoff s theorem is not hard. Simply consider the first image of the vertical line segment of height di. It can be shown ( ) I1 I ( 0 θ1 that the image of the line segment under F approximately has slope ). Taking sups over all such partial derivatives in Ω yields D +. Doing the same for F 1 I 0 yields D. Recall, that each Aubry-Mather set is a Lipschitz graph (Mather et al. 1994). We obtain the following corollary. Corollary 7 Each invariant curve in Ω, at every point, has a slope inside [D, D + ]. The same is true about the Lipschitz constant for each Aubry-Mather set contained in Ω. Consider tangent space orbits (F, DF ) : (x, v) (F (x), D x F (v)), where x = (θ, I) T R and v = (δθ, δi) T R. Theorem 8 (Cone Crossing Condition) Suppose {x i = F i (x 0 )} is an orbit and {v i } is the tangent component to the orbit. Additionally suppose v 0 is in the upper half cone (i.e. above lines of slope D ± originating at x 0 ), and suppose there exists an n > 0 so that v n is in the lower half cone (i.e. below the lines of slope D ± originating at x n ). Then the orbit does not belong to an invariant curve. Fig. 4 Cone Crossing Condition Proof: By continuity in initial conditions, a nearby orbit would cross the invariant curve. This is a contradiction (see fig. 4). MacKay and Percival include exposition regarding the equivalence of the cone crossing condition and action-minimizing orbits (MacKay et al. 1985). In short,

13 Instabilities in the Sun-Jupiter-Asteroid Three Body Problem 13 orbits which satisfy the cone crossing condition are not action-minimizing. This is because invariant curves are action-minimizing. See the associated paper for additional comments on other types of action minimizing (Aubry-Mather) sets; also see section Some preliminary numerics In order to use the cone crossing condition in practice, tight bounds on the slopes of the cones must be obtained. It is inefficient to use the uniform upper bound on the entire phase space Ω, especially when the objects considered lay only in a smaller subset. The following algorithm attempts to rule out invariant curves by obtaining better bounds using more localized information about potentially invariant objects. 1. Input an initial condition x 0 = (θ 0, I 0 ) and a number of iterates n. In practice, n is chosen so that θ n θ 0 > 2π. 2. Compute a bound I(x 0, n) so that x i = F i µ(x 0 ) T [I 0 I(x 0, n), I 0 + I(x 0, n)] for i = 0,..., n. The interval is often referred to as a localization interval because it localizes the curve Compute D ± on the annulus T [I 0 I(x 0, n), I 0 + I(x 0, n)]. 4. Use the cone crossing condition to rule out invariant curves on the annulus T [I 0 I(x 0, n), I 0 + I(x 0, n)] using n iterates. It suffices to use the vector (1, D + ) for initial conditions of the equations of variation and observe whether images in the tangent space drop below the vector (1, D ). One does not need to apply the cone crossing condition at every point in the annulus to rule out invariant curves. It suffices to use a vertical strip in the I direction going from the bottom to the top. If all points in the vertical strip satisfy the cone crossing condition, then there are no invariant curves. (If there were, they would have have to pass through the strip by construction of I.) If not all points in strip satisfy the cone crossing condition, the test to rule out invariant curves is inconclusive. In this case, a higher number of iterates may be required to rule out invariant curves, there may be point of a periodic orbit on the vertical slice, or there may be an invariant curve in the domain. The cone crossing condition only applies to invariant domains. While all the points in the region Ω twist by definition satisfy the twist property, this domain might not be invariant. By first constructing bounded regions where potential invariant curves must lay, then using the cone-crossing condition to rule out the possibility of existence invariant curves in those regions, this algorithm no longer requires invariance of the domain of twist. Instead, only the weaker condition that the localization interval remain inside in the domain of twist is required. Precise numerical estimates are required for concrete problems to ensure this condition is met. In practice this creates problems ruling out invariant curves in domains close to the boundary of Ω twist, however for a sufficiently large domain of twist, there 5 This construction actually localizes all Aubry-Mather sets with rotation symbol ω [ n+1 +, 1 ]. See section 6 for definitions. 1 n

14 14 John C. Urschel, Joseph R. Galante will be a sizable region where the presence of invariant curves can be ruled out. Remarks on Application to the RCP3BP: For J 0 = 1.55, µ = 0.001, select T [0.39, 0.83] (ϕ, P ϕ ) 6 to pick points from. It follows from Claim 5 that F µ is an EAPT in this domain. Step 2 of the algorithm, the construction of the localization intervals, can be estimated simply estimating the largest jump in the action component, P ϕ. From the equations of motion for the RCP3BP map, this can be estimated by integrating upper bounds on the perturbation term over one revolution. A priori estimates for this are not hard to obtain and this step can be quick if one is willing to accept a larger localization interval. Step 3 can be implemented using interval arithmetic to bound domains on a vertical strip. MacKay and Percival adopted this approach, though for our purposes this is quite expensive. Step 4 requires bounds on the equations of variation. A computer is programmed to compute the map F µ using µ = and J 0 = The space (ϕ, P ϕ ) [0, 2π] [0.57, 0.8] (corresponding to e [0.2, 0.6]) is divided up into boxes of size π 16 by To compute localization intervals, 10 iterates of the map are computed and the difference from initial conditions is measured. In this region of the phase space the size of the localization interval is no larger than 39µ. Using the localization intervals, estimates of D ± are computed using all points in the grid. Unfortunately, the cone crossing condition is not usually satisfied using this technique. Better estimates are needed. Using either a smaller grid or better estimates on both the perturbation terms (section 7) and on the flow is a possible approach. We pursue the method that MacKay and Percival suggested of refining estimates on the bounds D ± and combine it with approximations specific to our map F µ to obtain the cone crossing condition. Specifically, MacKay and Percival note that the map DF induces a map on cones. Under repeated forward iterates of an initial vertical vector, the images of the vectors decrease uniformly, converging to an eigenvector of DF. Similarly for DF 1 (MacKay et al. 1985). This has the effect of refining the cones. Armed with sufficiently good approximations, a computer is able to estimate D ± in each strip of size 0.01 in P ϕ and find that D Similarly, D This might suggest some type of symmetry with respect to D + and D. However, we do not assume anything of the sort. Refining at each point produces a measurable increase the bounds obtained for D ± (4.5µ vs 39µ). However this is still not sufficient to obtain a large region with no invariant curves. Below we shall make explicit the approximations used to obtain the 4.5µ bound, and provide an additional method of refinement which shall ultimately yield our desired result. 6 From hereon out, results may be stated using the (ϕ, P ϕ) parameterization of annulus, as opposed to (ϕ, e). The former parameterization is easier to work with numerically; the later is better for intuition.

15 Instabilities in the Sun-Jupiter-Asteroid Three Body Problem 15 Remark: This method has the advantage that it does not require us to compute a fixed number of iterates n when running the algorithm to test for invariant curves. Instead, we can run a test trajectory for as long as needed to determine if the cone crossing condition is satisfied. 5.2 Some approximations To compute the improved cones in the case of the RCP3BP, the map DF needs to be precisely computed. This is rather complicated in the case of flows (MacKay and Percival used the standard map which is cheap to compute). Let us make several approximations to F µ. First, in the definition of F µ, assume that µf ϕ (.; µ) 0. This is reasonable, as changing ϕ 1 by a µ small quantity over one revolution has little effect, since the ϕ components enter into motions only in the H term, which is already µ small. The effect of µf Pϕ (.; µ) is significantly more important, as it dictates how diffusion in eccentricity occurs. To compute this quantity exactly requires integration of the equations of motion for the RCP3BP, which can be computationally expensive. Instead, we use the integrable 2BP(SA) to compute this term. With the exception of high eccentricity sun grazing comets, starting from the same initial conditions, flows for the RCP3BP and for the 2BP(SA) are quite similar (in fact O(µ) close) over one revolution, so this is a reasonable approximation. First note the following identities for the 2BP (which can be found on pg of the text by Arnol d, Kozlov, and Neishtadt (2006)): u e sin(u) = L 3 t = l r = L 2( 1 e cos(u) ) sin(ϕ ϕ 0 ) = sin(u) 1 e 2 1 e cos(u) cos(ϕ ϕ 0 ) = cos(u) e 1 e cos(u) where ϕ 0 is the angle the asteroid makes with respect to the x-axis (where Jupiter is in the RCP3BP) when it is at the perihelion, i.e., when u = 0 mod 2π and t = 0. Addition formulas for sine and cosine give that sin(ϕ) = sin(u) 1 e 2 cos(ϕ 0 ) + cos(u) e 1 e cos(u) 1 e cos(u) sin(ϕ 0) cos(ϕ) = cos(u) e 1 e cos(u) cos(ϕ 0) sin(u) 1 e 2 sin(ϕ 0 ) 1 e cos(u) We make the estimate µf Pϕ (ϕ 0, P ϕ0 ; µ) p(ϕ 0, P ϕ0 ; µ) := = T 0 2π 0 H ( ) r(t), ϕ(t) 2BP ϕ (SA) dt H ϕ ( )( du ) 1 2BP r(u), ϕ(u) dt (SA) du.

18 18 John C. Urschel, Joseph R. Galante Definition: A conefield C(X) is a collection ( x, v (x), v + (x) ) for each x X where v ± (x) T X are vectors in the tangent space with base point x. In the present application, the vectors in the conefield are taken to the be eigenvectors of Jacobian DF of the EAPT F at each point in the domain Ω twist. MacKay and Percival (1985) proposed to iteratively refine conefields to prove the non-existence of invariant curves. The flow acts upon a conefield at a point through the action of the Jacobian matrix DF on the eigenvectors. For the RCP3BP, define the following refinements: v new + (x) := ( DG ) ( G (x)[ v+ G inv (x) )] (5) inv v new (x) := ( DG ) inv [ ( )] v G(x) G(x) (6) Remark: Further refinements can be obtained by using higher iterates and composition. The refinements work by flowing the eigenvectors for the map DG forward and backwards under the flow in the tangent space. This can be parlayed into an algorithm to refine conefields over an entire space. MacKay and Percival note that if at any point v lays above 7 v +, then the interior of the cone is empty and there is no invariant curve through the phase space at that point. This can be parlayed into the following algorithm. 1. Input n, m and divide the phase space (ϕ, P ϕ ) T R into blocks of size 1/n 1/m. 2. In each block, compute bounds on the eigenvectors of DG and DG inv. Use these eigenvectors to define the conefield for each point in the block. 3. Use refinements 5 and 6 to refine the conefields in each block, taking upper and lower bounds where appropriate. If v lays above v +, conclude there is no invariant curve inside of the block. If v lays below v +, take another refinement, or alternatively, stop and leave the block as inconclusive. A version of this algorithm is implemented. Specifically, a computer is programmed to compute the maps G, G inv and their derivatives, using µ = and J 0 = The space (ϕ, P ϕ ) [0, 2π] [0.57, 0.8] (corresponding to e [0.2, 0.6]) is divided up into boxes of size π 16 by The initial conefields are computed using the eigenvectors. The refinement algorithm is applied to all points on the grid. This provides strong numerical evidence of the following: Claim: For J 0 = 1.55, µ = 0.001, (ϕ, P ϕ ) [0, 2π] [0.57, 0.8] (corresponding to e [0.2, 0.6]) the map F µ has no invariant curves. Numerical Proof: A computer is programmed to compute the maps G, G inv, and their derivatives, using µ = and J 0 = The space (ϕ, P ϕ ) [0, 2π] [0.57, 0.8] (corresponding to e [0.2, 0.6]) is divided up into boxes of size π 16 7 For v = (a, b), let v = (1, b a ) be a normalization of v. A vector v 1 lays above v 2 iff x y, where v 1 = (1, x) and v 2 = (1, y).

19 Instabilities in the Sun-Jupiter-Asteroid Three Body Problem 19 by The initial conefields are computed using the eigenvectors for DG. The refinement algorithm given above is used for all points on the grid. Each box in the given domain requires only a finite number of iterates (an average of 15) to satisfy the cone crossing condition. Since the map G is µ close to F µ, then the eigenvectors for DG are of order µ close to those of DF µ. For a modest number of iterations, this difference appears to be negligible. Hence showing that the cone crossing condition holds for points in the domain using G implies it holds for F µ as well. Remark: While this provides strong numerical evidence that there are no invariant curves for F µ in the domain e [0.2, 0.6], the above numerical proof is not rigorous. To make it rigorous, upper and lower bounds for the cone fields and their refinements must be given for each box, not just for a representative point in the box. Interval arithmetic may be used to do this, as suggested by both Mackay and Percival (1985) and Galante and Kaloshin (2011). Alternatively, one may dispense with the map G entirely and compute everything using the map F µ at the cost of increased complexity in the interval arithmetic. We have shown that F µ cannot have invariant curves in the region corresponding to e [0.2, 0.6]. This region is contained inside the twist region Ω twist, however Ω twist is not invariant. It is not obvious that Aubry-Mather theory can be applied to this entire non-invariant region. Abstractly, Aubry-Mather (AM) sets arise as minimizers to a variational principle. To ensure their existence, it must be shown that orbits do not stray far from the domain where the minimization is taking place. One way to see this is the case for the RCP3BP is to recall that in section 5.1, localization intervals were constructed for Aubry-Mather sets. Numerical estimates show that orbits with initial conditions inside of T {e [0.2, 0.6]} remain inside of T {e [0.09, 0.8]} Ω twist after one full revolution around the sun. Independent of these concrete estimates, there are theoretical bounds which arise due to the action-minimizing properties of Aubry-Mather sets that keep the orbits we are interested in from straying too far from AM sets. These issues, along with a general review of Aubry-Mather theory is addressed in section 6. In summary, it turns out that e [0.2, 0.6] is large enough for minimizers and connecting orbits to exist and remain bounded inside the twist region. 6 Aubry-Mather Theory A compact invariant region C is bounded by two rotationally invariant curves C and C + such that there are no rotationally invariant curves in between C and C + is called a Birkhoff Region of Instability (BRI). In such BRIs, Birkhoff showed the existence of orbits coming arbitrarily close to C and C + (see the work of Mather and Forni (1994)). A much stronger result is given by Mather (1990), which allows one to specify neighborhoods of certain invariant sets which the orbit must pass through. Before stating this result, a quick overview of basic facts of Aubry-Mather theory shall be given. This review is primarily drawn from the works of Bangert (1988), Mather and Forni (1994), Gole (2001), Moser (1986),

6. Metric spaces In this section we review the basic facts about metric spaces. Definitions. A metric on a non-empty set X is a map with the following properties: d : X X [0, ) (i) If x, y X are points

Basic Concepts of Point Set Topology Notes for OU course Math 4853 Spring 2011 A. Miller 1. Introduction. The definitions of metric space and topological space were developed in the early 1900 s, largely

Duality of linear conic problems Alexander Shapiro and Arkadi Nemirovski Abstract It is well known that the optimal values of a linear programming problem and its dual are equal to each other if at least

Section 4: The Basics of Satellite Orbits MOTION IN SPACE VS. MOTION IN THE ATMOSPHERE The motion of objects in the atmosphere differs in three important ways from the motion of objects in space. First,

Terms and Definitions Absolute Value the magnitude of a number, or the distance from 0 on a real number line Additive Property of Area the process of finding an the area of a shape by totaling the areas

Mathematics Course 111: Algebra I Part IV: Vector Spaces D. R. Wilkins Academic Year 1996-7 9 Vector Spaces A vector space over some field K is an algebraic structure consisting of a set V on which are

+a 1. R In this and the next section we are going to study the properties of sequences of real numbers. Definition 1.1. (Sequence) A sequence is a function with domain N. Example 1.2. A sequence of real

USING MS EXCEL FOR DATA ANALYSIS AND SIMULATION Ian Cooper School of Physics The University of Sydney i.cooper@physics.usyd.edu.au Introduction The numerical calculations performed by scientists and engineers

Chapter 3 Stochastic Inventory Control 1 In this chapter, we consider in much greater details certain dynamic inventory control problems of the type already encountered in section 1.3. In addition to the

4.B. Tangent and normal lines to conics Apollonius work on conics includes a study of tangent and normal lines to these curves. The purpose of this document is to relate his approaches to the modern viewpoints

15 Limit sets. Lyapunov functions At this point, considering the solutions to ẋ = f(x), x U R 2, (1) we were most interested in the behavior of solutions when t (sometimes, this is called asymptotic behavior

1 Distance J Muscat 1 Metric Spaces Joseph Muscat 2003 (Last revised May 2009) (A revised and expanded version of these notes are now published by Springer.) 1 Distance A metric space can be thought of

AST1100 Lecture Notes 3 Extrasolar planets 1 Detecting extrasolar planets Most models of star formation tell us that the formation of planets is a common process. We expect most stars to have planets orbiting

Senior Secondary Australian Curriculum Mathematical Methods Glossary Unit 1 Functions and graphs Asymptote A line is an asymptote to a curve if the distance between the line and the curve approaches zero

CHAPTER II THE LIMIT OF A SEQUENCE OF NUMBERS DEFINITION OF THE NUMBER e. This chapter contains the beginnings of the most important, and probably the most subtle, notion in mathematical analysis, i.e.,

MA 323 Geometric Modelling Course Notes: Day 02 Model Construction Problem David L. Finn November 30th, 2004 In the next few days, we will introduce some of the basic problems in geometric modelling, and

DIFFERENTIAL GEOMETRY HW 2 CLAY SHONKWILER 2 Prove that the only orientation-reversing isometries of R 2 are glide reflections, that any two glide-reflections through the same distance are conjugate by

Chapter 1 Metric Spaces Many of the arguments you have seen in several variable calculus are almost identical to the corresponding arguments in one variable calculus, especially arguments concerning convergence

Appendix A Set Theory This chapter describes set theory, a mathematical theory that underlies all of modern mathematics. A.1 Basic Definitions Definition A.1.1. A set is an unordered collection of elements.

1 Local Brouwer degree Let D R n be an open set and f : S R n be continuous, D S and c R n. Suppose that the set f 1 (c) D is compact. (1) Then the local Brouwer degree of f at c in the set D is defined.

Chapter 3 Sequences In this chapter, we discuss sequences. We say what it means for a sequence to converge, and define the limit of a convergent sequence. We begin with some preliminary results about the

POWER SETS AND RELATIONS L. MARIZZA A. BAILEY 1. The Power Set Now that we have defined sets as best we can, we can consider a sets of sets. If we were to assume nothing, except the existence of the empty

MA651 Topology. Lecture 6. Separation Axioms. This text is based on the following books: Fundamental concepts of topology by Peter O Neil Elements of Mathematics: General Topology by Nicolas Bourbaki Counterexamples

Stationary random graphs on Z with prescribed iid degrees and finite mean connections Maria Deijfen Johan Jonasson February 2006 Abstract Let F be a probability distribution with support on the non-negative

arxiv:1205.5492v1 [math.co] 24 May 2012 Partitioning edge-coloured complete graphs into monochromatic cycles and paths Alexey Pokrovskiy Departement of Mathematics, London School of Economics and Political

Chapter 27 Solutions PSS 27.2 The Electric Field of a Continuous Distribution of Charge Description: Knight Problem-Solving Strategy 27.2 The Electric Field of a Continuous Distribution of Charge is illustrated.

Kepler s Laws of lanetary Motion and Newton s Law of Universal Gravitation Abstract These notes were written with those students in mind having taken (or are taking) A Calculus and A hysics Newton s law

Mathematics 31 Pre-calculus and Limits Overview After completing this section, students will be epected to have acquired reliability and fluency in the algebraic skills of factoring, operations with radicals

73 Appendix A.1 Basic Notions We take the terms point, line, plane, and space as undefined. We also use the concept of a set and a subset, belongs to or is an element of a set. In a formal axiomatic approach

Factoring Patterns in the Gaussian Plane Steve Phelps Introduction This paper describes discoveries made at the Park City Mathematics Institute, 00, as well as some proofs. Before the summer I understood

The Calculus of Functions of Several Variables Section. Introduction to R n Calculus is the study of functional relationships and how related quantities change with each other. In your first exposure to

Understanding Basic Calculus S.K. Chung Dedicated to all the people who have helped me in my life. i Preface This book is a revised and expanded version of the lecture notes for Basic Calculus and other

Some Comments on the Derivative of a Vector with applications to angular momentum and curvature E. L. Lady (October 18, 2000) Finding the formula in polar coordinates for the angular momentum of a moving

RANDOM INTERVAL HOMEOMORPHISMS MICHA L MISIUREWICZ Indiana University Purdue University Indianapolis This is a joint work with Lluís Alsedà Motivation: A talk by Yulij Ilyashenko. Two interval maps, applied

Solving Simultaneous Equations and Matrices The following represents a systematic investigation for the steps used to solve two simultaneous linear equations in two unknowns. The motivation for considering

PSEUDOARCS, PSEUDOCIRCLES, LAKES OF WADA AND GENERIC MAPS ON S 2 Abstract. We prove a Bruckner-Garg type theorem for the fiber structure of a generic map from a continuum X into the unit interval I. We

MATH 337 Cardinality Dr. Neal, WKU We now shall prove that the rational numbers are a countable set while R is uncountable. This result shows that there are two different magnitudes of infinity. But we

MATHEMATICS: THE LEVEL DESCRIPTIONS In mathematics, there are four attainment targets: using and applying mathematics; number and algebra; shape, space and measures, and handling data. Attainment target

150 Chapter 7 Dynamical Systems Mathematical models of scientific systems often lead to differential equations in which the independent variable is time. Such systems arise in astronomy, biology, chemistry,

Content Area: Mathematics Grade Level Expectations: High School Standard: Number Sense, Properties, and Operations Understand the structure and properties of our number system. At their most basic level

Linear Programming Linear programming refers to problems stated as maximization or minimization of a linear function subject to constraints that are linear equalities and inequalities. Although the study

Orbital Dynamics with Maple (sll --- v1.0, February 2012) Kepler s Laws of Orbital Motion Orbital theory is one of the great triumphs mathematical astronomy. The first understanding of orbits was published

CARDINALITY, COUNTABLE AND UNCOUNTABLE SETS PART ONE With the notion of bijection at hand, it is easy to formalize the idea that two finite sets have the same number of elements: we just need to verify

Continued Fractions and the Euclidean Algorithm Lecture notes prepared for MATH 326, Spring 997 Department of Mathematics and Statistics University at Albany William F Hammond Table of Contents Introduction

Student Performance Q&A: AP Calculus AB and Calculus BC Free-Response Questions The following comments on the free-response questions for AP Calculus AB and Calculus BC were written by the Chief Reader,

CHAPTER 1 General theory of stochastic processes 1.1. Definition of stochastic process First let us recall the definition of a random variable. A random variable is a random number appearing as a result

Lecture 13 Gravity in the Solar System Guiding Questions 1. How was the heliocentric model established? What are monumental steps in the history of the heliocentric model? 2. How do Kepler s three laws

An example of a computable absolutely normal number Verónica Becher Santiago Figueira Abstract The first example of an absolutely normal number was given by Sierpinski in 96, twenty years before the concept

AP Calculus BC Course Description: Advanced Placement Calculus BC is primarily concerned with developing the students understanding of the concepts of calculus and providing experiences with its methods