This project develops a computational approach to the inverse Sitnikov problem. Finite-time Lyapunov exponent (FTLE) analysis is used to characterize the transition between regular and chaotic dynamics for different eccentricities. Symbolic dynamics is then employed to search for initial conditions that realize prescribed crossing-time sequences by iteratively refining regions of phase space. The methodology is demonstrated for representative eccentricity values, illustrating the effectiveness of the computational framework and providing insight into the structure of the inverse problem.
Introduction
Introduction
Sitnikov Problem
Sitnikov Problem
In 1961, K. A. Sitnikov published his research on the existence of oscillatory motion in the elliptic restricted three body problem [1]. The setup involves two equally massed gravitational primaries with similar Keplerian orbits and a third body with negligible mass moving entirely along the normal of the orbits and intersecting the orbital plane at the barycenter. Therefore, the third body has 2 degrees of freedom and so its motion can be described with the z position and the velocity in the z-direction ). The configuration of primaries also depends on the eccentricity (e) and the phase of the primaries in their respective orbits (ϕ) that influences the oscillatory motion of the third body.
(
z
The Sitnikov theorem focuses on the oscillation patterns for every plane crossing/zero crossing (z = 0) and states that for sufficiently small e, there exists unique integer m, such that for any sequence [..... , , , , .....] where ∀ k, >= m, there exists an initial condition for the third primary in the Sitnikov problem. This happens due to the dynamics affected by the phase introduced by a non-zero eccentricity. The definition of a Sitnikov sequence and related annotations are discussed in Sitnikov sequence setup.
There are two important questions answered in this work - how to compute the unique integer m as a function of e (P1), and, for a given e and any sequence of integers, how to compute set of initial conditions satisfying the Sitnikov ODE as an inverse Sitnikov problem (P2). To answer the later question, stepwise mappings are implemented for a given eccentricity to segment the phase space into sets of Sitnikov integer. To answer the first question, the Finite-Time Lyapunov Exponent (FTLE) analysis over all the eccentricities determine the boundary of chaotic behavior providing rough estimates of the integer m.
s
k-1
s
k
s
k+1
s
k
There are two important questions answered in this work - how to compute the unique integer m as a function of e (P1), and, for a given e and any sequence of integers, how to compute set of initial conditions satisfying the Sitnikov ODE as an inverse Sitnikov problem (P2). To answer the later question, stepwise mappings are implemented for a given eccentricity to segment the phase space into sets of Sitnikov integer. To answer the first question, the Finite-Time Lyapunov Exponent (FTLE) analysis over all the eccentricities determine the boundary of chaotic behavior providing rough estimates of the integer m.
The reduced form of the third body’s z-dynamics is dependent on two parameters - eccentricity (e) and phase (ϕ). Interestingly, the phase can also be represented with time because of the primaries’ periodicity (T) simply denoted as t mod T. The normalized dimensionless ordinary differential equation (ODE) for the Sitnikov problem with a time period of 1 is given as:
In[]:=
ode[e_]:=Collectz[t]+z[t]==0,z''[t]
∂
t,t
3/2
(+)
2
z[t]
2
(1/2(1+eSin[2πt]))
Sitnikov Sequence Setup
Sitnikov Sequence Setup
Fundamentally, the Sitnikov sequence is a set of integers which is calculated by tracking the total number of complete revolutions of the primaries between two consecutive zero crossing of the third body. Mathematically, these integers are defined as for the crossing. The sequence is infamously known to be periodic for e=0 and the Sitnikov sequence completely depends on the speed at any zero crossing. This is however not true for e 0 cases where the sequence integers are based on the configuration in its last zero crossing.
⌊-⌋
t
k+1
t
k
th
k
th
k
≠
The zero-crossings can be determined for a phase space defined by ( t, ) by propagating the motion and storing the positions at every z=0 crossing.
z
In[]:=
SitnikovSolveCrossing[tend_,x0_,e_,nmax_,opt:OptionsPattern[NDSolve]]:=Module[{count=0,sol,data},{sol,data}=Reap[NDSolve[{ode[e],z[0]==x0[[1]],z'[0]==x0[[2]],WhenEvent[z[t]==0,Sow[{t,z'[t]}];count++;If[count>=nmax,"StopIntegration"]]},z,{t,0,tend},opt]];<|"Solution"->(z/.First[sol]),"Crossings"->If[data==={},{},First[data]]|>]SitnikovSequence[crossings_]:=Floor/@Differences[Prepend[crossings[[All,1]],0]];
Visually the sensitivity of the sequence with eccentricity intuitively explains the chaotic behavior in the system due to the introduction of nonlinearity in time in the ODE. Note that the initial z-position is taken as 0 which is true for the rest for the work because any such initial condition can be propagated back in time to get range of initial conditions for non-zero z and .
z
Out[]=
P1: Finite-Time Lyapunov Exponent Analysis
P1: Finite-Time Lyapunov Exponent Analysis
According to the Sitnikov theorem, there exists a unique integer m = m(e) for a small enough eccentricity e, such that any inverse Sitnikov function would be able to produce output given all the integers in the sequence are greater or equal to m. Intuitively this is a direct consequence of accumulating eccentricity effects on sufficiently large velocities to enter a chaotic motion (high sensitivity to initial conditions) allowing the third body to explore the space of higher velocities.
The separation of near deterministic to chaotic motion for the same system is called separatrix and it’s existence in the phase space is a rough approximation of the m=m(e) function for a given e.
The initial condition sensitivity is measured using finite-time Lyapunov exponent (FTLE) metric that measures the divergence of the image curves for a huge sequence of crossings. By tracking the divergence in the perturbations within the initial conditions, a positive FTLE value for a given initial condition means higher divergence rate compared to a negative FTLE value that means a higher convergence rate. A value near zero means a well deterministic motion for the given initial conditions. We again take a maximum integration time limit of to avoid errors.
10
10
In[]:=
r[t_,e_]:=(1/2)(1+eSin[2Pit]);varCoeff[t_,dz_,z_,e_]:=dz''[t]-dz[t](2z[t]^2-r[t,e]^2)/(z[t]^2+r[t,e]^2)^(5/2)==0ftlePoint[{t0_,v0_},e_,nCross_Integer]:=Module{tCur=t0,vCur=v0,logSum=0.,tFinal=t0,tNext,vNext,growth,k=0},Whilek<nCross,tNext=Null;NDSolve{ode[e],varCoeff[t,dz,z,e],z[tCur]==0,z'[tCur]==vCur,dz[tCur]==1.,dz'[tCur]==0.,WhenEvent[z[t]==0&&t>tCur+0.01,{tNext=t;vNext=z'[t];growth=dz[t];"StopIntegration"}]},{z,dz},{t,tCur,},Method->"StiffnessSwitching",MaxSteps->;If[!NumericQ[tNext],Break[]];logSum+=Log[Abs[growth]];tFinal=tNext;tCur=tNext;vCur=vNext;k++;If[tFinal>t0,logSum/(tFinal-t0),Missing[]]
10
10
6
10
The FTLE heatmap is created for visualizing the structure in the phase space. The FTLE values are simply calculated over a grid of initial conditions in the phase space and is visualized accordingly:
grid=ParallelTable[{t0,v0,ftlePoint[{t0,v0},ecc,nCross]},{t0,0.,1.,1./(nGrid-1)},{v0,-3,3,6/(nGrid-1)}];ftleData=Select[Flatten[grid,1],NumericQ[#[[3]]]&&#[[3]]>-5.&];ListDensityPlotftleData,
P2: Symbolic Dynamics Representation
P2: Symbolic Dynamics Representation
Stepwise Mapping on Discretized Phase Space
Stepwise Mapping on Discretized Phase Space
The so-called 2D positions in the phase space is defined by the phase and the z-velocity of the third body ( { ( - )mod1, } ). There are essentially two stepwise mapping important to address inverse Sitnikov problem - Poincare mapping () & Integer mapping ().
t
k+1
t
k
z
P
e
γ
e
The Poincare stepwise mapping takes a position from the phase space and maps it to another position by numerically integrating the motion until the next zero crossing.
In[]:=
PoincareOneStep[ti_,xi_,e_]:=FlattenReapNDSolve{ode[e],z[ti]==0,z'[ti]==xi,WhenEvent[z[t]==0,{Sow[{Mod[t-ti,1],z'[t]}],"StopIntegration"}]},z[t],{t,ti,}[[2]];Pestep[x_]:=Module[{p=PoincareOneStep[x[[1]],x[[2]],ecc]},If[p==={},{Infinity,0},p]];
10
10
The Integer stepwise mapping takes a position from the phase space and computes the respective Sitnikov integer for the given eccentricity.
In[]:=
γ[ti_,xi_,e_]:=FlattenReapNDSolve{ode[e],z[ti]==0,z'[ti]==xi,WhenEvent[z[t]==0,Sow[Floor[t-ti]];"StopIntegration"]},z[t],{t,ti,}[[2]];
10
10
These maps are used for segmentation of the phase space by collecting all the positions together that have same Sitnikov integers as a set of points in the variable IntervalSets:
data=Table[Flatten[{{y,dz},If[PoincareOneStep[y,dz,ecc]==={},{Infinity,Infinity},PoincareOneStep[y,dz,ecc]],If[γ[y,dz,0.1]==={},{Infinity},γ[y,dz,ecc][[1]]]}],{y,0,1,1/1000},{dz,-3,3,6/1000}];dd=Transpose@DeleteCases[Transpose[data],{{}..}];IntervalSets=<||>;setintegerlist=Flatten@dd[[All,All,5]];inputvzlist=Flatten@dd[[All,All,2]];inputt0list=Flatten@dd[[All,All,1]];IntervalSets=Merge[Thread[setintegerlist->MapThread[List,{inputt0list,inputvzlist}]],Identity]KeyDropFrom[IntervalSets,Infinity];
Visually these integers/intervals forms unique disjoint strips in the phase space as represented:
Interval Strip Intersections
Interval Strip Intersections
Inverse Sitnikov Function
Inverse Sitnikov Function
Defining a method to determine Sitnikov integer for the images contained in the argument x. This is required since the IntervalSet is based of discrete phase space and therefore images would most definitely always contain positions that are not considered by the variable IntervalSets. Intervalproximity method’s purpose is to output list of Sitnikov integers by comparing the input image positions with the nearest positions within IntervalSets and assigning the corresponding Sitnikov integer.
Defining a method to determine intersection of the images of a strip in the argument curve, with the strip corresponding to the Sitnikov integer m.
Defining the final inverse Sitnikov function by iteratively intersecting the images with the strip corresponding to the next integer in the sequence. This would produce set of initial conditions that would satisfy the input sequence. Note if the phase space contained negative velocities then there would always exist paired solutions to this function since a sequence if satisfied according to Sitnikov theorem would have two solutions, positive and negative velocities with same speed.
eccentricity = 0.1
eccentricity = 0.1
eccentricity = 0.5
eccentricity = 0.5
Conclusion
Conclusion
This work explores a symbolic dynamics computational approach in performing an inverse Sitnikov function based on the Sitnikov theorem for the guaranteed existence of initial conditions for certain satisfying sequences. Since the system quickly becomes non-integrable for non-zero eccentricities, the phase space search seems to be the most optimal direction for a solution however also deeming to be a memory intensive method. The computation of unique integer for a given eccentricity that separates the state divergence region with more deterministic region is quantified using the FTLE values which shows distinct separation regions for low and higher values of eccentricity. Moreover, once this unique integer is determined, the second part of the problem of determining the inverse Sitnikov function is approached with discretizing the phase space and realizing phase space strips corresponding to each Sitnikov integer. The image curves of an initial condition from the input sequence that satisfies strip intersection at every step of the Poincare map is ultimately selected to be the initial conditions for the inverse Sitnikov problem.
For smaller eccentricities, the separatrix is well-structured and the corresponding integer m has a relatively higher value. Moreover, the domain of initial conditions with near deterministic motion is relatively larger and almost independent of the phase of the primaries when compared to larger eccentricity cases. This explains why we see regions corresponding to small Sitnikov integer values to be shrinking with increasing eccentricities. One of the major limitations of this technique is its dependence on how accurate the IntervalSets represent the phase space and because of this the ICgenerator would fail to identify initial conditions for a valid sequence. Considering that, for higher eccentricities the ICgenerator is more prone to fail for a long enough sequence since the strips in the chaotic regime has large FTLE values and therefore their stretching and folding per crossings are more intense than lower eccentricities.
For smaller eccentricities, the separatrix is well-structured and the corresponding integer m has a relatively higher value. Moreover, the domain of initial conditions with near deterministic motion is relatively larger and almost independent of the phase of the primaries when compared to larger eccentricity cases. This explains why we see regions corresponding to small Sitnikov integer values to be shrinking with increasing eccentricities. One of the major limitations of this technique is its dependence on how accurate the IntervalSets represent the phase space and because of this the ICgenerator would fail to identify initial conditions for a valid sequence. Considering that, for higher eccentricities the ICgenerator is more prone to fail for a long enough sequence since the strips in the chaotic regime has large FTLE values and therefore their stretching and folding per crossings are more intense than lower eccentricities.
Future Scope
Future Scope
1
.For FTLE, a more adaptive grid resolution for varying thickness of strips corresponding to higher Sitnikov number will provide a better FTLE heatmap for a larger amount of crossings.
2
.Use an adaptive region determination to fully realize the shape of the strips for more efficient strip intersections
3
.The analysis was done over 1 million positions in the phase space however there are still many valid sequences where multiple applications of Poincare mapping stretches the strips high enough to conclude no intersections after a high finite number of crossings. Therefore due to this dependency, a more granular discretization of phase space would correspond to handling longer sequences as input.
Acknowledgments
Acknowledgments
I would like to express my sincere gratitude to my mentors, Pavel Hájek and Jesús Montesinos, whose guidance, patience, and intellectual generosity were instrumental in shaping this project from its earliest conception to its final form. Their deep expertise and thoughtful direction not only sharpened the mathematical and computational foundations of this work but also inspired a broader appreciation for the beauty of dynamical systems. A special thanks to Pedro Márquez-Zacarías for advising and providing sufficient computing resources for the project. I am equally grateful to Stephanie Bowyer for her exceptional organizational dedication in making the program run seamlessly, ensuring that every participant could focus entirely on their research without distraction. Her efforts, alongside those of all the mentors involved, fostered a remarkably engaging and collaborative research environment that made this experience as enriching as it was intellectually stimulating. Finally, I owe a special debt of gratitude to Stephen Wolfram, who generously took the time to personally help me identify a research direction ideally suited to my background, and whose vision for the program created a professional and inspiring setting throughout. This work would not have been possible without the collective support and encouragement of everyone involved.
References
References
1
.Dvorak, Rudolf and Lhotka, Christoph (2014), “Sitnikov problem” Scholarpedia. 10.4249/scholarpedia.11096.
2
.Moser, Jurgen (2016), Stable and Random Motions in Dynamical Systems: With Special Emphasis on Celestial Mechanics (AM-77). Princeton University Press. 978-1-4008-8269-4
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Inverting the Sitnikov Problem
by Suryansh Aryan
Wolfram Community, STAFF PICKS, July 16, 2026
https://community.wolfram.com/groups/-/m/t/3762762
by Suryansh Aryan
Wolfram Community, STAFF PICKS, July 16, 2026
https://community.wolfram.com/groups/-/m/t/3762762