Numerical Orbit Integrations Research Articles

Non-linear dynamics of a test particle near the Lagrange points L4 and L5 (Earth-Moon and Sun-Earth case)

The two-bodies problem can be fully solved, and was solved by Kepler (1609) and Newton (1687). The general three-body problem is often given as an example of a mathematical problem that ‘can’t be solved’. So, there is no general analytical solution. This problem can be significant and a special case of this problem is the Circular Restricted Three-Body Problem (CRTBP), which can be applied to the Earth-Moon system with a spacecraft, the Sun-Earth system with an asteroid, etc. In this paper, let’s focus on the motion of a test particle near the triangular Lagrange points L4 and L5 in the Earth-Moon and the Sun-Earth systems. Studying the movement of an object around these points is especially important for space mission design. To generate a trajectory around these points, the non-linear equations of motion for the circular restricted three-body problem were numerically integrated into MATLAB® 2023 software and the results are presented in the plane (x, y) and the phase plane (x, vx) and (y, vy). By numerical orbit integration, it is possible to investigate what happens when the displacement is relatively large or short from the Lagrange points. Then the small astronomical body may vibrate around these points. The results in this paper are shown in the rotating and inertia axes. Various initial positions near the Lagrange points and velocities are used to produce various paths the test particle can take. The same examples of numerical studies of trajectories associated with Lagrange points are shown in the inertial and the rotating coordinates system and are discussed. From the results of the numerical tests performed in MATLAB® 2023, it is possible to saw that there are different types of periodic, quasi-periodic, and chaotic orbits

Pitch angle scattering of fast particles by low frequency magnetic fluctuations

The adiabatic invariance of the magnetic moment during particle motion is of fundamental importance to the dynamics of magnetized plasma. The related rate of pitch angle scattering is investigated here for fast particles that thermally stream through static magnetic perturbations. For a uniform magnetic field with a localized perturbation, it is found that the curvature parameter κ2=min(Rc/ρL) does not predict the level of pitch angle scattering. Instead, based on numerical integration of particle orbits in prescribed magnetic fields, we derive predictions for the particle scattering rates, which can be characterized by the relative perturbation amplitude A≃δB/B, and the particle Larmor radius normalized by the field-aligned wavelength of the perturbations, ρL/λ||. Particles with ρL/λ||≃1 are subject to strong pitch angle scattering, while the level of scattering vanishes for both the limits of ρL/λ||≪1 and ρL/λ||≫1. The results are summarized in terms of a scattering operator, suitable for including the described scattering in basic kinetic models.

Relativistic effects on triple black holes: Burrau’s problem revisited

ABSTRACT We explore, using numerical simulations, the influence of mass and distance on the evolution of triple black hole systems. Following in the direction of Burrau’s famous 3,4,5 problem, black holes are initially placed at the vertices of Pythagorean triangles. Numerical integration of orbits was conducted using relativistic corrections (post-Newtonian) up to the 2.5th order with ARCcode. As a descriptor of the evolution of the systems, the lifetimes, the number of two-body encounters and the number of mergers were all analysed. We found that as the mass unit of the black holes increased, there was strong positive correlation with the fraction of mergers (0.9868), strong negative correlation with the average number of two-body encounters (−0.9016), and the average lifetimes of the triple systems decayed exponentially (determination coefficient of 0.9986). Around the mass unit range of 105.5–105.6M⊙, there was a transition from escape dominated dynamics to merger dominated dynamics. However, in the mass unit range of 106–109M⊙ with 1 pc distance unit, we find that 25 per cent of cases resulted in the escape of a supermassive black hole (SMBH) which may be a cause for wandering SMBH’s found in some galaxies/galactic merger remnants.

Numerical integration of particle orbits in discontinuous fields using VENUS-LEVIS and SPEC

The orbit code VENUS-LEVIS is adapted to follow particles in plasma equilibria with discontinuous fields generated by the Stepped Pressure Equilibrium Code (SPEC). The latter is an implementation of the Multi-Region relaxed MHD model, which efficiently computes Taylor states in a series of nested toroidal volumes and supports the formation of magnetic islands and chaotic field regions. To adapt VENUS-LEVIS, an event location procedure is implemented in the existing numerical integrator, which ensures the particle sees the correct field along its trajectory, regardless of the discontinuities present in the Stepped Pressure Equilibrium model. The algorithm is tested in the case where the magnetic field is uniform in the upper and lower half-spaces but has a discontinuity in its direction (shear) on the plane z=0. Particle drifts due to the discontinuity are studied. The convergence properties of the numerical scheme are highlighted by the numerical accuracy, and conservation of the system's invariants, such as energy and momentum. Simulations and convergence studies using the SPEC-LEVIS interface in axisymmetric geometry are then presented. Finally, illustrative particle drifts due to the discontinuity are studied and explained: we examine drifts associated with the change in Larmor radius of passing particles with small excursion from flux surfaces.

Finite gyroradius multidimensional electron hole equilibria

Finite electron gyroradius influences on the trapping and charge density distribution of electron holes of limited transverse extent are calculated analytically and explored by numerical orbit integration in low to moderate magnetic fields. Parallel trapping is shown to depend upon the gyroaveraged potential energy and to give rise to gyroaveraged charge deficit. Both types of average are expressible as convolutions with perpendicular Gaussians of width equal to the thermal gyroradius. Orbit-following confirms these phenomena but also confirms for the first time in self-consistent potential profiles the importance of gyrobounce resonance detrapping and consequent velocity diffusion on stochastic orbits. The averaging strongly reduces the trapped electron deficit that can be sustained by any potential profile whose transverse width is comparable to the gyroradius rg. It effectively prevents equilibrium widths smaller than ∼rg for times longer than a quarter parallel-bounce-period. Avoiding gyrobounce resonance detrapping is even more restrictive, except for very small potential amplitudes, but it takes multiple bounce periods to act. Quantitative criteria are given for both types of orbit loss.

HST Proper Motions of NGC 147 and NGC 185: Orbital Histories and Tests of a Dynamically Coherent Andromeda Satellite Plane

Abstract We present the first proper-motion (PM) measurements for the dwarf elliptical galaxies NGC 147 and NGC 185, two satellite galaxies of M31, using multiepoch Hubble Space Telescope (HST) imaging data with time baselines of ∼8 yr. For each galaxy, we take an error-weighted average of measurements from HST Advanced Camera for Surveys/WFC and Wide Field Camera 3/UVIS to determine the PMs. Our final results for the PMs are (μ W, μ N)N147 = (−0.0232, 0.0378) ± (0.0143, 0.0146) mas yr−1 for NGC 147 and (μ W, μ N)N185 = (−0.0242, 0.0058) ± (0.0141, 0.0147) mas yr−1 for NGC 185. The 2D direction of motion for NGC 147 about M31 is found to be aligned with its tidal tails. The 3D positions and velocities of both galaxies are transformed into a common M31-centric coordinate system to study the detailed orbital histories of the combined M31+NGC 147+NGC 185 system via numerical orbit integration. We find that NGC 147 (NGC 185) had its closest passage to M31 0.3–0.5 Gyr (≳1.6 Gyr) within the past 6 Gyr at distances of ∼70 kpc (70–260 kpc). The pericentric times of NGC 147/NGC 185 correlate qualitatively well with the presence/absence of tidal tails seen around the galaxies. Our PMs show that the orbital poles of NGC 147, as well as NGC 185, albeit to a lesser degree, agree within the uncertainties with the normal of the Great Plane of Andromeda (GPoA). These are the first measurements of the 3D angular momentum vector of any satellite identified as original GPoA members. Our results strengthen the hypothesis that the GPoA may be a dynamically coherent entity. We revisit previous claims that NGC 147 and NGC 185 are binary galaxies and conclude that it is very unlikely that the two galaxies were ever gravitationally bound to each other.

Particle Trapping in Axisymmetric Electron Holes

AbstractElectron orbits are calculated in solitary two‐dimensional axisymmetric electrostatic potential structures, typical of plasma electron holes, in order to establish the conditions for the particles to remain trapped. Analytic calculations of the evolution of the parallel energy caused by the perturbing radial electric field (breaking magnetic moment invariance) are shown to agree well with full numerical orbit integration Poincaré plots. The predominant mechanism of detrapping is resonance between the gyrofrequency in the parallel magnetic field and harmonics of the parallel bounce frequency. A region of phase space adjacent to the trapped‐passing boundary in parallel energy is generally stochastic because of island overlap of different harmonics, but except for very strong radial electric field perturbation, more deeply trapped orbits have well‐defined islands and are permanently confined. A simple universal quantitative algorithm is given, and its results plotted as a function of magnetic field strength and hole radial scale length, determining the phase space volume available to sustain the electron hole by depression of the permanently trapped distribution function.

The use of Kepler solver in numerical integrations of quasi-Keplerian orbits

ABSTRACT A Kepler solver is an analytical method used to solve a two-body problem. In this paper, we propose a new correction method by slightly modifying the Kepler solver. The only change to the analytical solutions is that the obtainment of the eccentric anomaly relies on the true anomaly that is associated with a unit radial vector calculated by an integrator. This scheme rigorously conserves all integrals and orbital elements except the mean longitude. However, the Kepler energy, angular momentum vector, and Laplace–Runge–Lenz vector for perturbed Kepler problems are slowly varying quantities. However, their integral invariant relations give the quantities high-precision values that directly govern five slowly varying orbital elements. These elements combined with the eccentric anomaly determine the desired numerical solutions. The newly proposed method can considerably reduce various errors for a post-Newtonian two-body problem compared with an uncorrected integrator, making it suitable for a dissipative two-body problem. Spurious secular changes of some elements or quasi-integrals in the outer Solar system may be caused by short integration times of the fourth-order Runge–Kutta algorithm. However, they can be eliminated in a long integration time of 108 yr by the proposed method, similar to Wisdom–Holman second-order symplectic integrator. The proposed method has an advantage over the symplectic algorithm in the accuracy but gives a larger slope to the phase error growth.

A Numerical Method for Determining the Elements of Circumbinary Orbits and Its Application to Circumbinary Planets and the Satellites of Pluto-Charon

Abstract Planets and satellites orbiting a binary system exist in the solar system and extrasolar planetary systems. Their orbits can be significantly different from Keplerian orbits, if they are close to the binary and the secondary-to-primary mass ratio is high. A proper description of a circumbinary orbit is in terms of the free eccentricity e free at the epicyclic frequency κ 0, forced eccentricity e forced at the mean motion n 0, and oscillations at higher frequencies forced by the non-axisymmetric components of the binary’s potential. We show that accurate numerical values for the amplitudes and frequencies of these terms can be extracted from numerical orbit integrations by applying fast Fourier transformation (FFT) to the cylindrical distance between the circumbinary object and the center of mass of the binary as a function of time. We apply this method to three Kepler circumbinary planets and the satellites of Pluto-Charon. For the satellite Styx of Pluto-Charon, the FFT results for κ 0 and e free differ significantly from the first-order analytic value and the value reported by Showalter & Hamilton, respectively. We show that the deviation in κ 0 is likely due to the effect of the 3:1 mean-motion resonance and discuss the implications of the lower value for e free.

Hybrid-precision arithmetic for numerical orbit integration towards future satellite gravimetry missions

The precision of numerical orbit integration is crucial in precise orbit determination and gravity field recovery, especially in the context of future satellite gravimetry missions where nanometer-level inter-satellite measurements are provided by laser ranging interferometer. In this paper, we propose a hybrid-precision arithmetic in both Cowell’s and Encke’s formulations for numerical orbit integration after comprehensive analyses on the computational time costs and numerical errors by dividing the integration process into two parts, namely the increment calculation and the orbit propagation. Based on such a separation, we find that the round-off errors basically occur in the least time-consuming orbit propagation part while the more time-demanding increment calculation part is less sensitive to the round-off due to its small numerical magnitude. Therefore, it is natural in the hybrid-precision arithmetic to use the double- and quadruple-precision arithmetic in the increment calculation part and orbit propagation part, respectively. Numerical results demonstrate that the hybrid-precision arithmetic is as efficient as the double-precision method, while far better than its double-precision counterparts in terms of both accuracy and relative precision. In consideration of the ultra-high precision inter-satellite range and range rate measurements of future satellite gravimetry missions, the currently widely-used double-precision arithmetic cannot meet the numerical requirements of one-day arc orbit integration for the inter-satellite range measurement, in which numerical errors should be at least the same level of the measurement precision and had better be one-order magnitude smaller to fully avoid their contaminations, while for the inter-satellite range rate the double-precision Gauss-Jackson integrator is qualified. Instead, by implementation of hybrid-precision arithmetic the range rate precision is easily achieved at 10−12 m/s in terms of the maximum absolute error of the one-day arc by either Cowell’s or Encke’s formulation. Furthermore, the sub-nanometer-level range precision is obtainable in the Encke’s formulation with reference orbit selected as the best-fit one. Therefore, it is recommended to adopt the hybrid-precision arithmetic for orbit integration to fully exploit the ultra-high precise measurement from the pure numerical data processing point of view.

On the triple-star origin of the planetary nebula Sh 2-71

ABSTRACT Recent studies have indicated that triple-star systems may play a role in the formation of an appreciable number of planetary nebulae, however, only one triple central star is known to date (and that system is likely too wide to have had much influence on the evolution of its component stars). Here, we consider the possibility that Sh 2-71 was formed by a triple system that has since broken apart. We present the discovery of two regions of emission, seemingly aligned with the proposed tertiary orbit (i.e. in line with the axis formed by the two candidate central star systems previously considered in the literature). We also perform a few simple tests of the plausibility of the triple hypothesis based on the observed properties (coordinates, radial velocities, distances, and proper motions) of the stars observed close to the projected centre of the nebula, adding further support through numerical integrations of binary orbits responding to mass loss. Although a number of open questions remain, we conclude that Sh 2-71 is currently one of the best candidates for planetary nebula formation influenced by triple-star interactions.

Fast Estimation of Orbital Parameters in Milky Way-like Potentials

Orbital parameters, such as eccentricity and maximum vertical excursion, of stars in the Milky Way are an important tool for understanding its dynamics and evolution, but calculation of such parameters usually relies on computationally expensive numerical orbit integration. We present and test a fast method for estimating these parameters using an application of the Stäckel fudge, used previously for the estimation of action-angle variables. We show that the method is highly accurate, to a level of <1% in eccentricity, over a large range of relevant orbits and in different Milky Way-like potentials, and demonstrate its validity by estimating the eccentricity distribution of the RAVE-TGAS data set and comparing it with that from orbit integration. Using this method, the orbital characteristics of the ∼7 million Gaia DR2 stars with radial velocity measurements are computed with Monte Carlo sampled errors in ∼116 hours of parallelized cpu time, at a speed that we estimate to be ∼3 to 4 orders of magnitude faster than using numerical orbit integration. We demonstrate using this catalog that Gaia DR2 samples a large range of orbits in the solar vicinity, down to those with rap ≲ 2.5 kpc, and out to rperi ≳ 13 kpc. We also show that many of the features present in orbital parameter space have a low mean zmax, suggesting that they likely result from disk dynamical effects.

Modeling approaches for precise relativistic orbits: Analytical, Lie-series, and pN approximation

Accurate orbit modeling plays a key role in contemporary and future space missions such as GRACE and its successor GRACE-FO, GNSS, and altimetry missions. To fully exploit the technological capabilities and correctly interpret measurements, relativistic orbital effects need to be taken into account.Within the theory of General Relativity, equations of motion for freely falling test objects, such as satellites orbiting the Earth, are given by the geodesic equation. We analyze and compare different solution methods in a spherically symmetric background, i.e.for the Schwarzschild spacetime, as a test bed. We investigate satellite orbits and use direct numerical orbit integration as well as the semi-analytical Lie-series approach. The results are compared to the exact analytical reference solution in terms of elliptic functions. For a set of exemplary orbits, we determine the respective accuracy of the different methods.Within the post-Newtonian approximation of General Relativity, modified orbital equations are obtained by adding relativistic corrections to the Newtonian equations of motion. We analyze the accuracy of this approximation with respect to the general relativistic setting. Therefore, we solve the post-Newtonian equation of motion using the eXtended High Performance Satellite dynamics Simulator. For corresponding initial conditions, we compare orbits in the Schwarzschild spacetime to those in its post-Newtonian approximation. Moreover, we compare the magnitude of relativistic contributions to several typical perturbations of satellite orbits due to, e.g., solar radiation pressure, Earth’s albedo, and atmospheric drag. This comparison is done for our test scenarios and for a real GRACE orbit to highlight the importance of relativistic effects in geodetic space missions. For the considered orbits, first-order relativistic contributions give accelerations of about 20 nm/s2 and are dominant in the radial direction.

Triaxiality Inhibitors in N-body Simulations

Abstract Numerous previous studies have investigated the phenomenon wherein initially spherical N-body systems are distorted to triaxial shapes. We report on an investigation of a previously described orbital instability that should oppose triaxiality. After verifying the instability with numerical orbit integrations that extend the original analysis, we search for evidence of the instability in N-body systems that become triaxial. Our results highlight the difficulty of separating dynamical process from finite-N effects. While we argue that our analysis points to the presence of the instability in simulated triaxial systems, discreteness appears to play a role in mimicking the instability. This suggests that predicting the shapes of real-world systems, such as dark matter halos around galaxies, based on such simulations, involves more uncertainty than previously thought.

Space Motions of the Dwarf Spheroidal Galaxies Draco and Sculptor Based on HST Proper Motions with a ∼10 yr Time Baseline

Abstract We present new proper motion (PM) measurements of the dwarf spheroidal galaxies (dSphs) Draco and Sculptor using multiepoch images obtained with the Hubble Space Telescope ACS/WFC. Our PM results have uncertainties far lower than previous measurements, even those made with the same instrument. The PM results for Draco and Sculptor are , , and , , . The implied Galactocentric velocity vectors for Draco and Sculptor have radial and tangential components: , , and , . We study the detailed orbital histories of both Draco and Sculptor via numerical orbit integrations. Orbital periods of Draco and Sculptor are found to be 1–2 Gyr and 2–5 Gyr, respectively, accounting for uncertainties in the Milky Way (MW) mass. We also study the influence of the Large Magellanic Cloud (LMC) on the orbits of Draco and Sculptor. Overall, the inclusion of the LMC increases the scatter in the orbital results. Based on our calculations, Draco shows a rather wide range of orbital parameters depending on the MW mass and inclusion/exclusion of the LMC, but Sculptor’s orbit is very well constrained, with its most recent pericentric approach to the MW being 0.3–0.4 Gyr ago. Our new PMs imply that the orbital trajectories of both Draco and Sculptor are confined within the “Disk of Satellites,” better so than implied by earlier PM measurements, and likely rule out the possibility that these two galaxies were accreted together as part of a tightly bound group.

New members of Datura family

Asteroid families are groups of minor planets which have a common origin in catastrophic disruptions. Young asteroid families are very interesting because they represent the product of their parent body's fragmentation before orbital and physical evolutionary processes could have changed them. A group of minor asteroids associated with the largest body (1270) Datura is of particular interest because it has enough known members and resides in the inner part of the main asteroid belt and is easier to observe than families (with similar physical characteristics on their surfaces) at further distances. Up to now, 7 members of this family are known. Here we discuss three new members of the Datura Family: (338309) 2002 VR17, 2002 RH291 and 2014 OE206. To prove that these recently-discovered members belong to the Datura family, we conducted numerical orbit integrations with all gravitational perturbation over the last 800kyrs. In the results, we have found that (338309) 2002 VR17 and 2002 RH291 are very close to the mean orbit of this family throughout the calculation. In the case of 2014 OE206, it has a strongly chaotic orbit.

Formation dynamics in geostationary ring

A relative motion model for a satellite formation composed of two Earth-orbiting spacecraft located in the geostationary ring is developed taking into account major gravitational and non-gravitational forces. A previously existing model featuring perturbation due to \(J_2\) is enhanced by the perturbations due to solar radiation pressure arising from unequal area-to-mass ratios, as well as the secular and long-periodic gravitational perturbations due to the Sun and the Moon. The extended relative motion model is validated using several typical formation geometries against a reference generated by numerical integration of the absolute orbits of the two spacecraft. The results of this work can find application in future on-orbit servicing and formation flying missions in near-geostationary orbit.

TEMPORARY CAPTURE OF ASTEROIDS BY A PLANET: DEPENDENCE OF PROGRADE/RETROGRADE CAPTURE ON ASTEROIDS’ SEMIMAJOR AXES

ABSTRACT We have investigated the dependence of the prograde/retrograde temporary capture of asteroids by a planet on their original heliocentric semimajor axes through analytical arguments and numerical orbital integrations in order to discuss the origins of irregular satellites of giant planets. We found that capture is mostly retrograde for the asteroids near the planetary orbit and is prograde for those from farther orbits. An analytical investigation reveals the intrinsic dynamics of these dependences and gives boundary semimajor axes for the change in prograde/retrograde capture. The numerical calculations support the idea of deriving analytical formulae and confirm their dependence. Our numerical results show that the capture probability is much higher for bodies from the inner region than for outer ones. These results imply that retrograde irregular satellites of Jupiter are most likely captured bodies from the nearby orbits of Jupiter that may have the same origin as Trojan asteroids, while prograde irregular satellites originate from far inner regions such as the main-belt asteroid region.

HIGHLY STABLE EVOLUTION OF EARTH'S FUTURE ORBIT DESPITE CHAOTIC BEHAVIOR OF THE SOLAR SYSTEM

Due to the chaotic nature of the Solar System, the question of its dynamic long-term stability can only be answered in a statistical sense, e.g. based on numerical ensemble integrations of nearby orbits. Destabilization, including catastrophic encounters and/or collisions involving the Earth, has been suggested to be initiated through a large increase in Mercury's eccentricity (eM), with an estimated probability of ~1%. However, it has recently been shown that the statistics of numerical Solar System integrations are sensitive to the accuracy and type of numerical algorithm. Here I report results from computationally demanding ensemble integrations (N=1,600 with slightly different initial conditions) at unprecedented accuracy based on the full equations of motion of the eight planets and Pluto over 5Gyr, including contributions from general relativity. The standard symplectic algorithm produced spurious results for highly eccentric orbits and during close encounters, which were hence integrated with a suitable Bulirsch-Stoer algorithm, specifically designed for these situations. The present study yields odds for a large increase in Mercury's eccentricity that are less than previous estimates. Strikingly, in two solutions Mercury continued on highly eccentric orbits (after reaching eM values >0.93) for 80-100Myr before colliding with Venus or the Sun. Most importantly, none of the 1,600 solutions led to a close encounter involving the Earth or a destabilization of Earth's orbit in the future. I conclude that Earth's orbit is dynamically highly stable for billions of years, despite the chaotic behavior of the Solar System.

MEASUREMENT OF PLANET MASSES WITH TRANSIT TIMING VARIATIONS DUE TO SYNODIC “CHOPPING” EFFECTS

Gravitational interactions between planets in transiting exoplanetary systems lead to variations in the times of transit that are diagnostic of the planetary masses and the dynamical state of the system. Here we show that synodic "chopping" contributions to these transit timing variations (TTVs) can be used to uniquely measure the masses of planets without full dynamical analyses involving direct integration of the equations of motion. We present simple analytic formulae for the chopping signal, which are valid (generally <10% error) for modest eccentricities e <~ 0.1. Importantly, these formulae primarily depend on the mass of the perturbing planet, and therefore the chopping signal can be used to break the mass/free-eccentricity degeneracy which can appear for systems near first order mean motion resonances. Using a harmonic analysis, we apply these TTV formulae to a number of Kepler systems which had been previously analyzed with full dynamical analyses. We show that when chopping is measured, the masses of both planets can be determined uniquely, in agreement with previous results, but without the need for numerical orbit integrations. This demonstrates how mass measurements from TTVs may primarily arise from an observable chopping signal. The formula for chopping can also be used to predict the number of transits and timing precision required for future observations, such as those made by TESS or PLATO, in order to infer planetary masses through analysis of TTVs.