Shadowing unstable orbits of the Sitnikov elliptic 3-body problem
aa r X i v : . [ a s t r o - ph . I M ] M a y Mon. Not. R. Astron. Soc. , 1–9 (2002) Printed 15 November 2018 (MN L A TEX style file v2.2)
Shadowing unstable orbits of the Sitnikov elliptic 3-bodyproblem
D. J. Urminsky ⋆ Department of Physics and Centre for Computational Relativity and Gravitation, Rochester Institute of Technology,85 Lomb Memorial Drive, Rochester, NY, 14623, USA
15 November 2018
ABSTRACT
Errors in numerical simulations of gravitating systems can be magnified exponentiallyover short periods of time. Numerical shadowing provides a way of demonstratingthat the dynamics represented by numerical simulations are representative of truedynamics. Using the Sitnikov Problem as an example, it is demonstrated that unstableorbits of the 3-body problem can be shadowed for long periods of time. In addition, itis shown that the stretching of phase space near escape and capture regions is a causefor the failure of the shadowing refinement procedure.
Key words: celestial mechanics, stellar dynamics, methods: numerical
The sensitivity which N -body integrations exhibit to smallchanges in initial conditions and to numerical errors hasbeen an active area of research since Miller’s landmarkstudy. Miller (1964) demonstrated the exponential diver-gence of near-by orbits for systems with N
32 and foundthat the separation of nearby orbits increases rapidly whenclose binary interactions occur. He suggested that the diver-gence of near-by orbits is too rapid to be solely accounted bybinary interactions and suggests that there must be a col-lective effect to account for the results. However, Standish(1968) showed that the divergence rate was reduced if thepotential was replaced with a softened potential and con-cluded that the divergence is mainly due to close binaryinteractions.The dramatic effects of numerical errors on N -body in-tegrations was also demonstrated in an important paper byLecar (1968). After coordinating a study with 11 differentintegrations of the same 25-body problem for 2.5 crossingtimes, Lecar found that quantities such as half mass radiusand the moment of inertia can change by as much as 100percent. In a study with N = 3, Dejonghe and Hut (2001)demonstrated that the amplification of initial errors can in-crease by as much as 10 . In addition, they showed thatthe growth of errors during close encounters can be ampli-fied by as much as 10 , however some of the growth can berecovered after the encounter is over.The sensitivity to small changes in initial condi-tions and numerical errors is a property associated withchaotic systems. A measure of the sensitivity of numer- ⋆ E-mail: [email protected] ical errors can be determined by the Lyapunov expo-nent λ . Early work suggested that the Lyapunov expo-nent is inversely proportional to the crossing time t cr (Kandrup and Smith 1991; Heggie 1991; Goodman et al.1993). However, Goodman et al. (1993) suggest a depen-dence on N of the form λ − = t cr / log N or perhaps λ − = t cr / log(log( N )), implying that as N increases therate of separation decreases and the Lyapunov exponent in-creases. The log( N ) dependence was later numerically veri-fied by Hemsendorf and Merritt (1991).Despite the difficulty calculating solutions to N -bodyintegrations, computers still remain a useful tool to studyself gravitating systems. If numerical errors in numerical so-lutions to the N -body problem cause such drastic changesin the actual positions and velocities of particles how can wetrust the dynamics that these solutions represent? Shadow-ing is a way of proving that a true solution to a dynamicalsystem follows close to a numerical solution. If true orbitscan be found close to numerical orbits then the dynamicsrepresented by the numerical solutions represents true dy-namics.This study will discuss the existence of shadow orbitsfor the gravitational 3-body problem. First, definitions andconcepts related to shadowing of dynamical systems will beintroduced. Next, a refinement procedure which makes cor-rections to numerical orbits to reduce the errors incurred ateach time step will be presented. The Sitnikov problem willthen be presented and used as a simple model to discussescape and capture of orbits. An approximate Poincar´e mapis then presented to model orbits of the Sitnikov problemand will be used in conjunction with the refinement proce-dure to discuss the validity of numerical solutions by way ofshadowing. The failure of the refinement procedure to find c (cid:13) D. J. Urminsky shadow orbits will then be discussed and regions of phase-space where the procedure fails will be delineated. Finally, itwill be demonstrated that the shadow times for this problemcan be modeled as a Poisson process.
Consider the autonomous ordinary differential equation˙ x = f ( x ) , (1)where x ∈ R n and f : R n → R n is a C vector field withthe associated flow represented by ϑ t . A sequence of points { y k } Mk =0 is said to be a pseudo-orbit if there is an associ-ated bounded sequence { h k } Mk =0 of positive time such that, | y k +1 − ϑ h k ( y k ) | < δ, (2)for k = 0 , , ..., M , where δ >
0. An example of a pseudo-orbit is a numerical solution to (1). To show that a pseudo-orbit represents some true dynamics for (1), it would beenough to show that a true orbit follows close to the pseudo-orbit. The pseudo-orbit described above is said to be shad-owed by a true orbit if there is a sequence of points { x k } Mk =0 and positive times { t k } Mk =0 with ϑ t k ( x k ) = x k +1 such that | x k − y k | < ǫ, (3)and | t k − h k | < ǫ, (4)for k = 0 , , .., M and small ǫ >
0. The sequence { x k } Mk =0 is known as a shadow-orbit . The shadow-orbit is a truesolution to (1).The first general contributions made on shadowing fordynamical systems were the shadowing theorems of Anosov(1967) and Bowen (1972). Anosov and Bowen consideredhyperbolic systems and showed that any pseudo-orbit ona hyperbolic invariant set has a shadow-orbit. These the-orems were generalized for pseudo-orbits in the vicinity ofa hyperbolic set (Kato 1991; Nadzieja 1991; Coomes et al.1995). For non-hyperbolic systems or for orbits which arefar from hyperbolic invariant sets these theorems do notapply. Shadowing theorems do exist for pseudo-orbits ofnon-hyperbolic systems and usually rely on numerical ver-ification of a theorem (Coomes et al. 1994; Chow et al.1989; Chow and Palmer 1991; Chow and Van Vleck 1994;Van Vleck 1995). Shadowing N -body simulations was first demonstratedby Quinlan and Tremaine (1992) and Hayes (2001). Boththese studies considered the refinement procedure found inGrebogi et al. (1990) to find numerical shadows for the N -body problem. The refinement procedure is a noise reduc-tion technique which can be used to show the existence ofshadow-orbits. This procedure will be presented for two di-mensional dynamical maps, however the procedure can eas-ily be adapted for flows and has been extended to higherdimensional systems by Quinlan and Tremaine (1992).Consider the pseudo-orbit { p k } Mk =0 of a map f ∈ R .The goal is to find a new less noisy orbit { ˆ p k } Mk =0 close tothe original orbit. Let e k represent the one step error where e k = p k − f ( p k − ) . (5)The refined orbit is constructed byˆ p k = p k + Φ k , (6)where Φ k is the correction at time step k . Combine equa-tions (5) and (6) to obtain Φ k = f (ˆ p k − ) − e k − f ( p k − ) , (7)where ˆ p k = f (ˆ p k − ). Assuming that the correction, Φ k , issmall, expand f (ˆ p k − ) about p k − in a Taylor series to get, f (ˆ p k − ) ≈ f ( p k − ) + L k − Φ k − , (8)where L k is the linearized map at the k th time step. Sub-stitute (8) into (7) to obtain Φ k ≈ L k − Φ k − − e k . (9)It is also assumed that the linearized map has an expandingdirection, u k , and a contracting direction, s k , at each timestep k . With this assumption, the objective is to find thesequences { Φ k } Mk =0 and { e k } Mk =0 in the coordinates { u k } Mk =0 and { e k } Mk =0 by Φ k = α k u k + β k s k (10)and e k = η k u k + ζ k s k . (11)The expanding and contracting directions follow the lin-earized maps, u k +1 = L k u k , (12)and s k +1 = L k s k . (13)For a random | u | = 1, equation (12) gives u k aligned withunstable direction at p k after just a few iterations. Startingwith a random s M and iterating (13) backwards gives s k aligned with the stable direction at p k after a few iterations.Substitute (10) and (11) into (9) to get, α k +1 u k +1 + β k +1 s k +1 = L k ( α k u k + β k s k ) + ( η k +1 u k +1 + ζ k +1 s k +1 ) . (14)Substituting (12) and (13) into (14) yields recursive rela-tionships for { α k } Nk =0 and { β k } Nk =0 where, α k +1 = | L k u k | α k − η k +1 ,β k +1 = | L k s k | β k − ζ k +1 . (15)Equations (15) are made computationally stable by calcu-lating the coefficients α k starting with α M and iteratingbackwards and the coefficients β k are calculated by choos-ing an initial β and iterating forwards. The choice of α M and β are arbitrary and are taken to be α M = β = 0. Thusthe sequence of correction coefficients are given by α M = 0 , α k = ( α k +1 + η k ) / | L k u k | ,β = 0 , β k = | L k s k | β k − ζ k +1 , (16)where the values of η k and ζ k can be determined directlyfrom (5) and (11).Once { ˆ p k } Mk =0 has been found, the refinement procedurecan be iterated. Generally, the number of significant digitsdoubles on each iteration of the process. However, cases have c (cid:13) , 1–9 hadowing unstable orbits mL 3 m1m2 Figure 1.
The Sitnikov Problem been found where the convergence is much slower or does notconverge.The convergence of the refinement procedure does notin itself show the existence of a shadow-orbit. Grebogi et al.(1990) provide a containment procedure in two dimen-sions which rigorously proves the existence of a shadow-orbit. The containment technique was later extended tothree dimensional systems by Hayes (2001). A more practi-cal approach for higher dimensional systems was developedby Sauer and Yorke (1991). They showed that for a givenpseudo-orbit, if the refinement procedure converges - to ma-chine precision - and certain quantities of a theorem remainbounded, then the pseudo-orbit has a shadow-orbit.It has been found (Quinlan and Tremaine 1992; Hayes2003), that one can tell from the convergence of the refine-ment procedure alone whether a given pseudo-orbit can beshadowed. So, if iterations of the refinement procedure con-verge to a new orbit where the one-step errors are the sizeof machine precision, then it is inferred that a shadow-orbitexists for the given pseudo-orbit. The new orbit found bythe refinement procedure is called a numerical shadow . In this study of shadowing for the 3-body problem, a specialconfiguration of the restricted 3-body problem known as theSitnikov problem will be considered. The Sitnikov problemis the problem of the motion of a mass-less particle, m ,on the axis of symmetry of an equal-mass binary (Figure1). Following Moser (1973), units are chosen such that thegravitational constant G = 1 and the total mass M = 1.Under these conditions, the equation of motion for m isgiven by¨ z = − z √ z + r , (17)where z is the position of m and r the distance from thecentre of mass to one of the binary masses. The distance r can be approximated to first order in the eccentricity, e , by r ≈
12 (1 − e cos t ) , (18)and the specific energy of m can be defined by E = 12 | ˙ z | − √ r + z . (19)Taking the plane of motion of the binary ( z = 0) as aSurface Of Section (SOS), consider a map, φ : ( v , t ) → ( v , t ), which takes m from one crossing of the SOS tothe next. If m is on the SOS at time t , φ is a map whichbrings v = ˙ z ( t ) to time t > t where v = ˙ z ( t ) and z ( t ) = 0. The map φ has an open domain D in which everypoint returns to the SOS. As time enters into the problemwith period 2 π , D can be considered in polar co-ordinateswhere the radial variable is v and the angular variable isgiven by t . Alternatively, the domain D can be consideredon the surface of a cylinder where the initial position onthe cylinder is defined by t and E . Figure 2 (a) shows thedomain for φ in cylindrical coordinates. The colour of eachpoint represents the number of periods of the binary beforeescape happens. The green regions represent islands of quasi-periodic motion. In Figure 3 an example of a quasi-periodicorbit which visits the islands of stability in the vicinity of aperiod 7 orbit is provided. Urminsky and Heggie (2009) demonstrated that thePoincar´e map φ with (18) can be approximated by asimplectic map ϕ : ( t , E ) → ( t , E ) where E / = E + a cos( t ) + b sin( t ) ,t / = t + 2 C ( − E / ) − / ,t = t / + 2 C ( − E / ) − / ,E = E / − a cos( t ) + b sin( t ) , (20)and a , b and C are constants. The quantities t / and E / are approximations of the time and energy values of m ata local maximum distance from the SOS. It is clear (Figure2 (a)) that the change in energy of m from one crossing tothe next is periodic in time and the trigonometric terms in(20) can be though of as the lowest order in a Fourier ap-proximation to this change. The change in time is obtainedby approximating the motion of m as Keplerian. The con-stants a and b are proportional to the eccentricity of thebinary whose values can be shown to be, a ≈ . eb ≈ . e (21)and the constant C = π/ (2 √ Through interactions with the binary as it crosses the SOS, m can gain sufficient energy such that it leaves the SOSand does not return. It can be shown that for some positivetime t ∗ and positive ν = (1 − e ) /
2, if12 ˙ z ( t ∗ ) − z ( t ∗ ) + ν > , (22)then | z ( t ) | → t → ∞ . Setting z = 0 in (22) givesa lower bound on the velocity of orbits which escape onthe SOS. The solid black region at the top of Figure 2 (a)demonstrates how this condition over estimates the escapeboundary. All energy and time values in this region do notreturn to the SOS. c (cid:13) , 1–9 D. J. Urminsky (b) 0 1 2 3 4 5 6t -5-4-3-2-1 0 1 E -5-4-3-2-1 0 1 E Figure 2. (a) Domain in ( t , E )-space for the map φ where e = 0 .
61. Each point represents an initial condition and the associated colourrepresents the number of periods of the binary before escape. The green regions towards the bottom of the graph represent quasi-periodicorbits which remain bound for all time. The solid black region at the top of the image are initial conditions outside of the domain of φ .The escape criterion used is effective in determining escape but crude in approximating the escape boundary on the SOS. (b) Domain in( t , E )-space for the map ϕ where e = 0 .
61. The colour associated with each initial condition represents the number of periods of thebinary before escape. In (b), the number of periods of the binary is determined by t M / π where M is the number of iterations of themap. The green regions represent quasi-periodic orbits which do not escape. The map, ϕ , provides an accurate way of determiningescape and capture. From equation (20) it is found that themapping ϕ is defined in a region, E < − a cos( t ) − b sin( t ) := ∂ D , (23)for t ∈ [0 , π ] (24)as time enters into the mapping with period 2 π . The curve ∂ D is the escape boundary. Time and energy values above ∂ D are said to have escaped. The domain, D , can be de-fined by (23) and (24). Initial conditions in D are mappedinto the region, D , defined by E < − a cos( t ) + b sin( t ) := ∂ D , (25)for t ∈ [0 , π ]. Figure 4 shows how the boundaries ∂ D and ∂ D intersect. Orbits are mapped from the region under thecurve ∂ D to the region under the curve ∂ D . The region B = D \D represents energy and time values for which or-bits are captured. In the context of the differential equation,these are orbits which come from infinity and get captured by the binary. Similarly, the initial conditions in the region B = D \D are energy and time values for which ϕ is un-defined. Again, in the context of the differential equation,the region B represents orbits which escape from the sys-tem. Finally, note that initial conditions for the differentialequation are such that ˙ z > z = 0 on the SOS. So from(19), the initial energy can be bounded from below by, E > − / | r ( t ) | , (26)for t ∈ [0 , π ].Initial conditions are chosen in D with (26) for e = 0 . t M / π where M is the number of iterations of the map ϕ .The green regions represent stable motion whose orbits re-main bounded. As energy increases orbits become unstableand escape from the system. Notice the similarities betweenFigure 2 (a) and 2 (b). Both domains have islands repre-senting stable orbits as well as large regions representingunstable orbits. In Figure 5 an example of a quasi-periodicorbit near a period 7 orbit is provided. In addition to the c (cid:13) , 1–9 hadowing unstable orbits -4.5-4-3.5-3-2.5-2-1.5-1 0 1 2 3 4 5 6 E ne r g y time (mod 2 π ) Figure 3.
An example of a quasi-periodic orbit near a period 7orbit for equation (17) on the SOS z = 0. Initial conditions are˙ z (0) = 1 . z (0) = 0 . e = 0 .
61 and the phase of the binary is . ∂ D ∂ D D B B π Figure 4.
The curve ∂ D represents a lower bound of energy andtime values for which ϕ is undefined. Similarly, the curve ∂ D represents a lower bound of energy and time values for which theinverse map ϕ − is undefined. The shaded region is the domain D for the map ϕ . The two regions labeled B and B boundedby the curves ∂ D and ∂ D are the capture and escape regionsrespectively. similarities between Figure 2 (a) and (b), Urminsky (2009)demonstrates that the map ϕ satisfies Lemmas similar toLemmas 1-5 in Moser (1973) (pages 87-91) and that ϕ , like φ , possesses a hyperbolic invariant set on which ϕ is topo-logically equivalent to the shift map. The map ϕ provides a simple way of studying shadowing fororbits like those of the Sitnikov problem. The approximatemap is used to avoid integrating between successive cross-ings of the SOS thus obtaining a tremendous speed up incalculations. In addition, the one step error can more eas-ily be controlled. At each time step uniformly distributednoise | δ k | δ is added to generate the pseudo-orbit. Therefinement procedure is then used to reduce the noise levelto machine precision. Since ϕ is a 2-dimensional mapping, -2.8-2.6-2.4-2.2-2-1.8-1.6-1.4-1.2 0 1 2 3 4 5 6 E ne r g y time (mod 2 π ) Figure 5.
An example of a quasi-periodic orbit near a period7 orbit for the map ϕ . Initial conditions are t = 6 . E = − . e = 0 . the refinement procedure can be directly applied as shownin section 2.1. Using the containment and refinement procedure,Grebogi et al. (1990) successfully demonstrated theexistence of shadows for pseudo-orbits of length 10 ormore. To test the algorithm the refinement procedure isapplied to long lived orbits of the map ϕ . As seen in Figure2 (b), there are regions of stable motion where orbits remainbounded forever. The refinement procedure is applied tothese orbits and it is found that most can be shadowed formany iterations. Some of these are shown in Figure 9.As shown by Dvorak et al. (1998) for the Sitnikov prob-lem, the map ϕ has ‘sticky’ regions where orbits can betrapped for long periods of time before escape. In Figure 6an example of a sticky orbit trapped in the vicinity of is-lands of stable quasi-periodic orbits is shown. The inset plotin Figure 6 is a magnification of the orbit near one of theislands. By sampling the phase space around the islands ofstable motion, one can find many sticky orbits which survivefor long periods of time. In Figure 7 the shadow distance isplotted against the number of iterations of the map for sev-eral sticky orbits where e = 0 .
61. It is shown that as thenumber of iterations increases, the distance of the numer-ical shadow from the pseudo-orbit increases proportionallyto the number of iterations.
Consider uniformly distributed initial values in B (Figure4) for e = 0 .
25. Initial values are iterated forward for a max-imum of 100000 iterations up to the penultimate iterationbefore escaping. For each orbit, the number of iterations, M , the orbit was ‘shadow-able’ for as well as the shadow-distance are recorded. The orbits are binned into bins oflength one iteration and averaged over the bin. The resultsare plotted in Figure 8 where the dots represent the averageshadow-distance at each iteration of the map. Note that as c (cid:13)000
25. Initial values are iterated forward for a max-imum of 100000 iterations up to the penultimate iterationbefore escaping. For each orbit, the number of iterations, M , the orbit was ‘shadow-able’ for as well as the shadow-distance are recorded. The orbits are binned into bins oflength one iteration and averaged over the bin. The resultsare plotted in Figure 8 where the dots represent the averageshadow-distance at each iteration of the map. Note that as c (cid:13)000 , 1–9 D. J. Urminsky -2.8-2.6-2.4-2.2-2-1.8-1.6-1.4-1.2 0 1 2 3 4 5 6 E n t n -1.25-1.225 2.7 3 3.3 3.6 Figure 6.
An example of a ‘sticky’ orbit which remains close toislands of stable orbits for the map ϕ for 500 , S hado w d i s t an c e M (
Figure 7.
Number of iterations verses shadow distance ǫ . M increases, there is increasing variability on the distribu-tion of average shadow distances. The data can be fit withthe curve 7 × M which is similar to the results in Fig-ure 7 where the shadow distance is proportional to the orbitlength. Numerical shadows have been found using the refinementprocedure for orbits whose length exceeds 10 iterations forthe map ϕ . However, what happens when numerical shad-ows are not found? What causes the refinement procedure tofail? First, it should be noted that the failure of the refine-ment algorithm to converge to a numerical shadow does notimply that there is not a shadow-orbit for a given pseudo-orbit. A shadow may still exist but the refinement procedurewas not able to converge towards it. Quinlan and Tremaine(1992) and Hayes (2003) found that shadowing breaks downduring close encounters between particles. This is due to thestretching of the velocity subspace during a close encounter.In the Sitnikov problem, m interacts with the binary onthe SOS and the distance separating m with the binary s hado w d i s t an c e M (
Figure 8.
Each point is average shadowing distance for the as-sociate shadow length M for 10 initial conditions in D where e = 0 . masses is bounded from below (and above) on the SOS. Incontrast to the problems discussed in the above mentionedstudies arbitrarily close encounters do not occur in the Sit-nikov problem. However, escape and capture occur duringclose encounters with the binary as m crosses the SOS.Near the escape and capture boundaries slight changes inthe energy of m as it crosses the SOS can lead to signifi-cant changes in the duration of successive crossings of theSOS. The map provides a simple way of sampling the phasespace on the SOS to find regions where shadowing is morelikely to fail.Figure 9 shows shadowing results of 10 initial condi-tions. The colour of each point represents either success, yel-low, or failure, black, of the refinement procedure. Note thatonly orbits which survived more than three iterations of themap are considered. This is because the choice of u and s M would influence the results for short lived orbits as (12) and(13) may not have had enough time to align u k and s k inthe proper directions. From Figure 9 it can be seen that therefinement procedure tends to fail near the escape boundary ∂ D . Note also that the refinement procedure fails near theboundaries of regions containing orbits which escape afterthree or less iterations.The reason the refinement procedure fails in these re-gions is that there is a stretching of subspace as orbits nearthe boundary ∂ D . At a given iteration k , the distance fromboundary, ∂ D , is given by, d = | E k + a cos( t k ) + b sin( t k ) | . (27)From the Jacobian of (20) it can be shown that | L u k | ∼ d − / . (28)Thus, as d →
0, the correction coefficients α and β go to 0and ∞ respectively making it more difficult for the refine-ment procedure to converge.Figure 10 shows the density of successfully shadowedorbits based on the closest approach to the boundary ∂ D for increasing eccentricity values. For each shown eccen-tricity value, we select 100000 uniformly distributed ini-tial conditions in the region defined by t ∈ ( π, π ) and E ∈ ( ∂ D + 2 b sin( t ) , ∂ D ). These boundaries describea band of initial conditions bounded above by the escape c (cid:13) , 1–9 hadowing unstable orbits
0 1 2 3 4 5 6t -5-4-3-2-1 0 E Figure 9.
The figure shows one million initial conditions for themap ϕ where e = 0 .
61. The map ϕ was applied to each initial con-dition 50,000 times or until the resulting orbit escaped. The colourassociated with each point represents the successful applicationof the refinement algorithm. Black represents initial conditionswhere the refinement procedure failed to converge. Yellow repre-sents the successful application of the refinement procedure. Onlyorbits which were longer than three iterations of ϕ are considered. p r obab ili t y den s i t y d (distance from boundary)e=0.002e=0.102e=0.202e=0.302e=0.402e=0.502e=0.602 Figure 10.
Probability density of close approaches to the escapeboundary for shadow-able orbits. boundary. This band also encompasses the capture region B . The drop in the density to the right of each curve occursat the distance between the lower boundary curve and theescape boundary. Note that the density drops off as initialconditions approach the escape boundary. Data was fittedusing a variable bandwidth kernel density function. It was found above that as orbits approach the escapeboundary the likelihood of an orbit being shadowed de-creases. This has an impact on the shadow-ability of orbitsin the capture region B . The capture region area is directlyproportional to the eccentricity of the binary. As e → B become pushed up against the f r a c t i on s hado w ab l e e (eccentricity)data0.08 log(16948 e) Figure 11.
Fraction of capture orbits shadow-able using the re-finement procedure for increasing eccentricities of the binary. p r obab ili t y den s i t y T (shadow duration)data.0019e -.0019T .025/T1e-09 1e-07 1e-05 1e-031e3 1e4 1e5
Figure 12.
Probability density of shadow durations for the map ϕ where e=0.61. The amplitude of the one step noise was setat 10 − . For shadow durations T < T thedensity is inversely proportional to T . boundary ∂ D . It is expected then that for small eccentric-ities, orbits would be less likely to be shadow-able.To test this hypothesis, 10 uniformly distributed ini-tial conditions are selected in B and iterated forwards untileach orbit escapes. This is performed for a variety of eccen-tricity values and the fraction of shadow-able orbits in eachcase is determined. The results are shown in Figure 11. Thefraction of shadow-able orbits increases as the eccentricity ofthe binary increases. This is because the area of the captureregion increases proportionally to e . As the area increases,initial conditions can be selected at a much further distancefrom the escape boundary making them more likely shadow-able. The failure of the refinement procedure can happen at anypoint along the orbit and not necessarily at a close approachto the escape boundary. The shadow duration is defined asthe number of iterations for which a given orbit can be shad- c (cid:13) , 1–9 D. J. Urminsky p r obab ili t y den s i t y lifetime data.0005*e -.0005x Figure 13.
Probability density for the lifetime, P Tk =0 t k , for themap ϕ where e=0.61. The amplitude of the noise is 10 − . It wasfound that the distribution best fit an exponential distribution. owed. For an orbit { ( t i , E i ) } Mi =0 the shadow duration, T , cantake on positive integer values T < M .Consider the initial conditions for e = 0 .
61 shown inFigure 9. For each resulting orbit, it is determined how longthe orbit is shadow-able. Figure 12 provides some informa-tion on the distribution of shadow lengths. Initial condi-tions are chosen in D and iterated forwards in time us-ing (20). Each orbit is iterated for 50 ,
000 iterations or un-til the solution escapes. The solid line in Figure 12 repre-sents the density of the numerical experiments. The spikeat 50,000 iterations is mostly due to quasi-periodic orbitswhich remain bounded for all time. As shown in Figure12, the density can be approximated, for small iterations,by an exponential density function given by ξ exp( − ξx ) for ξ = 0 . < M < , . /M .The map approximates the time between crossings ofthe SOS by considering the motion of m to be Keplerian.Instead of considering the distribution of the shadow dura-tion in terms of the number of iterations of ϕ we can insteadconsider the distribution of shadow times, t M where M is thenumber of iterations of the orbits for which it was shadow-able. The solid line in Figure 13 represents the probabil-ity density of shadow time for the numerical experiments.Again, the data can best be approximated by the exponen-tial density function for ξ = . N -body systems have an exponential distributionand can be thought of as a Poisson process. The above results confirm, for short lived orbits, previousinvestigations (Hayes 2003) that showed numerical shadowdurations, M , for gravitating systems follow a Poisson pro-cess with a exponential density function. The result foundin this study suggests for longer lived orbits, the densityfunction is better approximated by a function proportionalto 1 /M . This may be because the population of longer lived orbits tends to be dominated by stable orbits, however thishas not been investigated.In section 5, areas of phase-space where the refinementprocedure is more likely to fail are characterized. These areasare near escape boundaries where there is sufficient stretch-ing of phase-space to cause the refinement procedure to failto converge to a less noisy orbit. Interestingly this seems tobe due to the growth of the variational equations over onetime step. This does not rule out the failure of the the re-finement procedure by the accumulative effect of the growthof the variational equation associated with large Lyapunovexponents as discussed by Zhu and Hayes (2009).In Figure 11 it is demonstrated that as the volumeof phase-space representing capture orbits decreases, it be-comes increasingly difficult to shadow capture orbits. Thisis a result of the distribution of failures of the refinementprocedure seen in Figure 10. As the volume of phase-spaceassociated with capture decreases, capture orbits get pushedup against the boundary ∂ D where the one-step growth ofthe variational equations causes the refinement procedure tofail. Finally, it was found that the shadow distance for anorbit is proportional to the number of iterations of the map(Figures 7 and 8). It was noticed that if in addition to t and E , orbits were required to be shadow-able at the half steps t / and E / , then initially shadow-able orbits continuedto be shadow-able. When shadowing at the half step wasrequired, the shadow distance typically increased by abouta factor of two.The Sitnikov problem discussed in this study providesa straight forward way of characterizing a domain of initialconditions as well as regions of stable and unstable motion.Work in progress considers slight changes to the Sitnikovproblem in order to study shadowing of unstable orbits. Forexample, Soulis et al. (2007) consider slight perturbationsto the mass and position (away from the z -axis) of m anddelineate regions of stable and unstable motion. It would beexpected that, like the results found in this study, shadow-ing with the refinement procedure breaks down near bound-aries of escape for unstable orbits. In fact, the break downof the refinement procedure near escape boundaries wouldbe expected for general 3-body configurations. As solutionsapproach parabolic escape boundaries, an orbit can undergoincreasingly long ejections from the left-over binary system.Small changes in the energy of an orbit in this region cancause significant changes in the time of return for the orbit.If the refinement procedure could make changes to the orbitsso as to conserve the energy of the ejected body it might im-prove the success rate of the refinement procedure. Finally,the Sitnikov 4-body problem (Soulis et al. 2008) provides astarting point for examining the relationship between theshadowing distance and the number of bodies. Extra bodiescan be added in circular orbits about the center of mass.Hayes (2003) demonstrates that as the number of movingbodies in a fixed potential increases the shadow durationsdecrease. It would be of interest to determine if a similarrelationship holds for the Sitnikov N -body problem.It should be stressed again that the failure of the refine-ment procedure does not necessarily mean that a shadowdoes not exist for a given pseudo-orbit. It may very well bethat shadows do exist for orbits in regions where the refine-ment procedure fails. We are encouraged that this may be c (cid:13) , 1–9 hadowing unstable orbits the case. Both the Sitnikov problem and the approximatePoincar´e map possess a hyperbolic invariant set, Λ, near theescape boundaries (see Moser (1973) and Urminsky (2009)respectively). Despite the fact that Λ is near the boundary ∂ D , the shadowing theorems by Anosov (1967) and Bowen(1972) guarantee that any pseudo-orbit on Λ has an asso-ciated shadow-orbit. This demonstrates that being in thevicinity on the escape boundary does not necessarily ruleout the existence of shadow-orbits. ACKNOWLEDGMENTS
DU was supported by the National Aeronautics and SpaceAdministration through grant NNX-07AH15G. The authorwould like to thank D. Heggie, D. Merritt and D. Dicken fortheir helpful suggestions. In addition, the author would liketo thank the anonymous referee for his/her careful readingof the manuscript and useful suggestions.
REFERENCES
Anosov D.V, 1967, Trudy Mat. Inst. Steklov, 90, 209Bowen R., 1972, Amer. J. Math., 94, 1Chow S., Lin X., Palmer K.J., 1989, Differential equations,Lecture Notes in Pure and Appl. Math., 118, Dekker, NewYork, 127Chow S., Palmer K.J., 1991, J. Dynam. Differential Equa-tions, 3, 3, 361Chow S., Van Vleck E., 1994, SIAM J. Sci. Comput., 15,4, 959Coomes B.A., Ko¸cak H., Palmer K.J, 1994, J. Comput.Appl. Math., 52, 35Coomes B.A., Ko¸cak H., Palmer K.J, 1995, Z. Angew.Math. Phys., 46, 1, 85Dejonghe H., Hut P., 1986, in Hut P., McMillan S., eds,The Use of Supercomputers in Stellar Dynamics, LectureNotes in Physics, 246, Springer, Berlin, 212Dvorak R., Contopoulos G., Efthymiopoulos C., 1998,Planet. Space Sci., 46, 1567Goodman J., Heggie D.C., Hut P., 1993, ApJ, 415, 213Grebogi C., Hammel S.M., Yorke J.A., Sauer T., 1990,Phys. Rev. Lett., 65, 1527Hayes W.B., 2001, Ph.D. thesis, University of Toronto,TorontoHayes W.B., 2003, Phys. Rev. Lett., 90, 5, 054104Heggie D.C., 1991, in Roy A.E., ed., Predictability, Sta-bility and Chaos in N -Body Dynamical Systems, PlenumPress, New York, 47Hemsendorf M., Merritt D., 2002, Apj, 580, 606Kandrup H.E., Smith H., 1991, ApJ, 374, 255Kato K., 1991, Mem. Fac. Sci. Kochi Univ. Ser. A Math.,21, 43Nadzieja T., 1991, Arch. Math. (Brno), 27A, 65Lecar M., 1968, Bull. Astron., 91, 3, 213Miller R.H., 1964, ApJ, 140, 250Moser J., 1973, Stable and Random Motions in DynamicalSystems, Princeton U. Press, PrincetonQuinlan G.D., Tremaine S., 1992, MNRAS, 259, 505Sauer T., Yorke J.A., 1991, Nonlinearity, 4, 961 Soulis P., Bountis T., Dvorak R., 2007, Celest. Mech. Dyn.Astron., 99, 129Soulis P.S., Papadakis K.E., Bountis T., 2008, Celest.Mech. Dyn. Astron., 100, 251Standish E.M., 1968, Ph.D. thesis, Yale University, NewHavenUrminsky D.J., Heggie D.C., 2009, MNRAS, 392, 1051Urminsky D.J., 2009, Ph.D. thesis, University of Edin-burgh, EdinburghVan Vleck E., 1995, SIAM J. Sci. Comput., 16, 5, 1177Zhu Y.-K., Hayes W.B., 2009, AAS/Division of DynamicalAstronomy Meeting, 40, 11 c (cid:13)000