An application of symplectic integration for general relativistic planetary orbitography subject to non-gravitational forces

Spacecraft propagation tools describe the motion of near-Earth objects and interplanetary probes using Newton’s theory of gravity supplemented with the approximate general relativistic n-body Einstein–Infeld–Hoffmann equations of motion. With respect to the general theory of relativity and the long-standing recommendations of the International Astronomical Union for astrometry, celestial mechanics and metrology, we believe modern orbitography software is now reaching its limits in terms of complexity. In this paper, we present the first results of a prototype software titled General Relativistic Accelerometer-based Propagation Environment (GRAPE). We describe the motion of interplanetary probes and spacecraft using extended general relativistic equations of motion which account for non-gravitational forces using end-user supplied accelerometer data or approximate dynamical models. We exploit the unique general relativistic quadratic invariant associated with the orthogonality between four-velocity and acceleration and simulate the perturbed orbits for Molniya, Parker Solar Probe and Mercury Planetary Orbiter-like test particles subject to a radiation-like four-force. The accuracy of the numerical procedure is maintained using a 5-stage, 10th\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^\mathrm{th}$$\end{document}-order structure-preserving Gauss collocation symplectic integration scheme. GRAPE preserves the norm of the tangent vector to the test particle worldline at the order of 10-32\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-32}$$\end{document}.

planetary spacecraft and near-Earth satellites using Newton's theory of gravity supplemented with the so-called Einstein-Infeld-Hoffmann equations of motion (Einstein et al. 1938;Soffel and Han 2019). The dynamical equations of motion characterising spacecraft orbits are given in terms of gravitational and non-gravitational perturbations, denoted, respectively, f G and f N G , which are linearly superimposed as small corrections to Newton's universal law of gravitation. Classically, non-gravitational forces are derived analytically using simplified ideal physical models, for example, the so-called cannonball model, and are the main source of error in the orbit determination process (Vallado 2001). Examples of several gravitational and non-gravitational perturbations, and their respective orders of magnitude can be found in (Pireaux et al. 2006, Table 1), and the reader is further referred to (Vallado 2001, Chapter 9) for a detailed discussion on classical perturbation techniques.
Due to the complexity involved in accurately determining spacecraft orbits, the system of equations considered by modern orbitography software require advanced numerical procedures. There exists a multitude of open-source and commercial software suites (Hughes et al. 2014;Vallado et al. 2006;Evans et al. 2018) designed to determine spacecraft ephemerides subject to a surfeit of physical and dynamical configurations. Although many platforms differ in their initial formulation (special and general perturbation approaches (Vallado et al. 2006), semi-analytical techniques, etc.), we may mathematically formulate the problem as follows. The approximate position of spacecraft subject to gravitational and non-gravitational perturbations is determined through numerical integration of the first coordinate-time derivative of the spatial and velocity components, denoted x i = (x 1 , x 2 , x 3 ) and v i = (v 1 , v 2 , v 3 ), respectively, given by the initial value problem where we denote the initial position and velocity of the spacecraft at a given initial epoch t 0 by x i (t 0 ) = x i 0 and v i (t 0 ) = v i 0 . We note that the aforementioned time parameter is absolute and universal in the framework of Newtonian mechanics. The expressions f i G and f i N G are introduced for brevity and denote the i-th component of the individual gravitational and non-gravitational perturbations which are given explicitly by where the spatial gradient operator is given by ∂ i ≡ ∂/∂ x i , the Newtonian gravitational potential is given by U and the subscripts j, k denote the j-th gravitational and k-th nongravitational perturbation, respectively. Owing to the increased accuracy requirements in fields such as astrometry (Soffel and Han 2019), geodesy (Müller et al. 2008;Soffel and Langhans 2012), the rapid improvements in time and frequency stability (Exertier et al. 2019), and further, the development and utilisation of advanced telecommunication systems for radio tracking of interplanetary probes, we require that Einstein's general theory of relativity (d'Inverno 1992;Weinberg 1972) be taken into account for any mission requiring highly accurate orbit information and practically all astronomical and geodynamical observations (Brumberg et al. 1998). Within the Solar System, the so-called effects of general relativity are accounted for using the first post-Newtonian (PN) approximation (Soffel and Han 2019) and we refer the reader to the seminal papers Damour et al. (1991Damour et al. ( , 1992bDamour et al. ( , 1993Damour et al. ( , 1994 for extensive technical details. The PN framework is a slow-motion, weak-field approximation to the Einstein field equations of gravitation and is a formal expansion of the metric tensor components in powers of order O(c −2 ), where c denotes the isotropic speed of light in vacuum. For an Earthorbiting satellite, the typical relativistic effects considered by modern orbitography software are the contributions due to the Schwarzschild field of the Earth, and the DeSitter and Lense-Thirring effects (Misner et al. 1973). As an example, consider the Earth-orbiting Laser Geodynamics Satellite LAGEOS. 1 With respect to the aforementioned relativistic contributions, the additional perturbing effect of the Schwarzschild field is largest in magnitude, causing a decrease in semi-major axis by 4.4mm and a perigee advance of approximately 9mas/day. The magnitude of the Lense-Thirring and DeSitter effects is considerably smaller and amount to a precession of the longitude of the ascending node by approximately 85 and 53μas/day, respectively, with a further precession of the argument of perigee caused by the Lense-Thirring effect of the order 31.6mas/year (Hugentobler 2010;Combrinck 2012). The magnitude of the DeSitter and Lense-Thirring effects as predicted by the general theory of relativity was experimentally verified by the Gravity Probe B mission using IM Pegasi as a guide star where the reader is referred to Everitt et al. (2011) for further details.
The general relativistic corrections to Newton's law of gravitation, denoted f G R , are derived analytically (to first PN order) and added linearly to Eq. (1) b . Hence, the complete system of equations of motion typically considered by modern orbitography software are given by Again, we remind the reader, the expression f G R is introduced for brevity and denotes the sum of the individual relativistic perturbations under consideration. The recommended analytical form for the general relativistic correction to Newton's law of gravitation f G R is maintained and published by the International Earth Rotation and Reference Systems Services (IERS) in a series of technical notes. 2 For the interested reader, derivations of the PN equations of motion can be readily found in the extensive literature (see, for example Damour et al. (1994); Brumberg (2004); Soffel (1989); Huang et al. (1990)). Although we may observe the correct dynamical contributions due to the general theory of relativity by numerically integrating (3), care must be taken. The naive process of linearly adding f G R to Eq. (1) in order to account for the so-called relativistic effects introduces a Newton-esque, or, what is more commonly referred to in the literature, a neo-Newtonian interpretation of the general theory of relativity (Damour et al. 1991). As discussed by many authors previously (Damour et al. 1991;Damour 1987;Damour et al. 1992a;Soffel 2000), the aforementioned reduction of general relativity to the framework of Newtonian theory can lead to the derivation of erroneous equations of motion and the apparent appearance of illusory defects in material bodies (Damour 1987). A further source of confusion may be linked to the recommendations of the International Astronomical Union (IAU) (Soffel et al. 2003), which, through a series of resolutions, suggest that all problems in the field of astronomy or astrodynamics be formulated within the framework of Einstein's general theory of relativity. Specifically, this relates to the derivation of all dynamical equation of motion including those which model the propagation of electromagnetic waves in curved spacetime, the definition of reference frames and the corresponding four-dimensional spacetime transformations connecting them. In order for the IAU resolutions to be sufficiently reflected in the so-called "Newton + correction" frame-work would not only require additional corrections to standard Newtonian quantities, but would further require that we abolish the classical notion of absolute time resulting in an additional hybrid Newton-Einstein theory of celestial mechanics. Considering the aforementioned shortcomings of the "Newton + correction" approach, and the confusion arising in adopting the recommended IAU resolutions in a Newton-dominated theory, we believe the current approach is now reaching its limits in terms of complexity.
In this paper, we discuss the formulation of a new prototype orbitography software titled General Relativistic Accelerometer-based Propagation Environment (GRAPE), formally extending the proposed work of Pireaux et al. (2006). We suggest a novel technique for the propagation of interplanetary probes and near-Earth spacecraft using the full framework of general relativity by accounting for the spatial components of non-gravitational forces using either end-user supplied data from on-board accelerometers (Lenoir et al. 2011;Touboul et al. 1999) or specific force models with adjustable parameters. Accelerometer-based orbit determination has been successfully exploited by the geodesy community for a number of geopotential recovery missions (see, for example Kang et al. (2006); Van Helleputte and Visser (2008)), where the dynamical models used to approximately account for the nongravitational perturbations f N G are replaced with accelerometer data obtained in the local frame of the spacecraft. Hence, errors introduced by simplified dynamical models, omitted second-order perturbations and unaccounted anomalous forces are automatically included at the order of the measurement sensitivity of the on-board accelerometer.
General relativity provides our currently accepted theory of gravitation. Within the framework of general relativity, quantities constructed from the metric tensor components g μν such as Christoffel symbols, curvature scalars and the Ricci, Riemann and Einstein tensors encode the geometry of spacetime which we associate with classical gravitational phenomenon. GRAPE assumes the local structure of spacetime is unperturbed by non-gravitational forces, and all associated gravitational effects are contained within the choice of spacetime metric. Consequently, the general relativistic analogue of the first two terms on the right-hand side of Eq. (3) is accounted for by a suitable choice of metric g μν . Further, we note the spacecraft equations of motion are always accompanied by a quadratic invariant quantity unique to general relativity; offering an opportunity to exploit sophisticated structure-preserving numerical integration techniques irrespective of whether non-gravitational perturbations are being considered (Hairer et al. 2006). Again, we note such an opportunity to exploit structurepreserving integration schemes is not readily available in classical Newtonian theory and is unique to the general theory of relativity which arises due to the orthogonality of four-velocity and acceleration. GRAPE employs a 5-stage, 10 th order symplectic implicit Runge-Kutta scheme (Butcher 1964) to maintain the constancy of the quadratic invariant providing a consistent accuracy check over the considered numerical propagation arc. Additionally, given that the magnitude of the non-gravitational perturbations is measured in the local frame of the spacecraft, GRAPE utilises the so-called tetrad formalism (Weinberg 1972), introducing at each point in spacetime, a locally Minkowski frame, which we associate with the local frame of the spacecraft (up to a Lorentz boost), allowing for convenient transformations, which relate the dynamical quantities dynamical quantities between local and global frames (Brumberg 2017).
The paper is organised as follows: in Sect. 2, we provide an overview of the mathematics employed by the GRAPE orbitography software. We derive the equations of motion for near-Earth objects and interplanetary probes subject to non-gravitational perturbations and discuss the associated quadratic invariant for spacecraft whose wordline is parameterised by proper time τ . Further, we present a brief overview of the tetrad formalism and discuss a novel method to determine the magnitude of the global force components obtained using both Lorentz and tetrad transformations which does not appear to have been previously discussed in the literature. In Sects. 3 and 4, we derive the general relativistic extension of the classical accelerometer equation using a modified geodesic deviation approach (Pireaux et al. 2006) and give an overview of the structure-preserving symplectic numerical integration scheme implemented in GRAPE, respectively. In Sect. 5, we provide a specific illustrative example of the GRAPE tool and consider the perturbed motion of the NASA Parker Solar Probe (Fox et al. 2016), Mercury Planetary Orbiter (Novara 2002) and Molniya-like spacecraft (Daquin et al. 2021). Additionally, we show the conservation of the quadratic invariant in the local frame of each spacecraft whose motion is perturbed by a constant radiation-like non-gravitational force.
Unless otherwise stated, throughout the remainder of the manuscript we assume the Einstein summation convention for repeated indices, and the use of Latin indices (i, j = 1, 2, 3) is reserved for spatial coordinates while Greek indices (α, β = 0, 1, 2, 3) denote spacetime coordinates with the 0 th coordinate being reserved for time so that x α = ct, x i . Finally, we assume a spacetime metric signature according to (+, −, −, −).

GRAPE: mathematical preliminaries
The worldline of a time-like test particle is described by the spacetime coordinates x α (λ), where λ is an arbitrary affine parameter. For a massive time-like particle, such as an interplanetary probe or near-Earth satellite, a common parameterisation of the spacecraft worldline is achieved using the proper time τ (Poisson and Will 2014). However, spacecraft operators track satellites, interplanetary probes, large debris objects and planets with reference to an external coordinate time such as terrestrial, universal, ephemeris or barycentric times (see Vallado (2001) for further time scales and definitions). Hence, for operational purposes, it is convenient to modify the parameterisation and adopt a global coordinate time t e = x 0 /c common to all spacecraft instead of τ as the integration variable. Each parameterisation offers different insight for perturbed spacecraft motion in the framework of general relativity. Accordingly, GRAPE numerically solves both systems of equations of motion and determines spacecraft orbits with respect to both proper time τ and Solar System ephemeris time t e . We note for the purpose of the present paper, we define ephemeris time as the time read by a clock at spatial infinity, void of gravitational effects, and draw the readers attention to (Soffel and Han 2019, Chapter 8) for formal operational details.

On the nature of non-gravitational forces in general relativity
Studies on the motion of near-Earth spacecraft and interplanetary probes subject to nongravitational forces are ubiquitous in the literature and remain an active area of research for modern operational astrodynamics and celestial mechanics (Vallado 2001). However, the same treatment cannot be said for the general theory of relativity where the practical incorporation of non-gravitational forces for near-Earth objects remains somewhat elusive and requires elucidation. The approach adopted by GRAPE follows from the work of Lichnerowicz and Teichmann (1955); Burcev (1962); Charon (1963) which extend the standard geodesic equation of motion to include non-gravitational forces. While the extended equations of motion derived by Lichnerowicz and Teichmann (1955); Burcev (1962); Charon (1963) appear to have important practical implications for applied general relativity, the utility of their results appears to have gone largely unnoticed in the community which we re-derive herein.
Classical dynamical quantities associated with continuous matter distributions such as stresses and mass, energy and flux densities are coalesced into a single unified object in the general theory of relativity known as the energy momentum tensor T αβ (Kopeikin 2007). To account for non-gravitational forces in the framework of general relativity, we decompose the energy-momentum tensor as follows where we define the stress tensor according to S αβ and we interpret Θ αβ as the energymomentum of an isolated test particle or system of non-interacting test particles defined in terms of mass density ρ and test particle four-velocity u α according to Θ αβ ≡ ρu α u β . It is well established that the conservation equation gives rise to the geodesic equation of motion in the absence of stresses (S αβ ≡ 0). In the case where S αβ = 0, we find where we have introduced the force density K β associated with an arbitrary non-gravitational force f β and we have ∇ α S αβ = −ρK β . Eq. (6) is equivalently expressed as where the second term on the left hand side of Eq. (7) is a continuity-like equation which is defined in terms of the force density according to Thus, the complete system of equations of motion describing the perturbed motion of spacecraft subject to an external perturbation f α employed in GRAPE are given by where it is clear for K β = 0, Eq. (9) reduces to the equation of geodesic motion and we introduce the Christoffel symbols of the second kind according to Γ α βγ ≡ g αδ ∂ β g γ δ + ∂ γ g δβ − ∂ δ g βγ /2.
Following several applications of the chain rule so that and using (9) in expressions (11), it can be readily shown that the perturbed equations of motion with respect to ephemeris time are given by where it is clear that the temporal component of (12) is identically zero. Hence, when integrated with respect to ephemeris time, the motion of spacecraft depends only on the spatial coordinates x i and Eq. (12) may be equivalently expressed as where we have introduced the dimensionless quantity β i ≡ v i /c and the contravariant components of the coordinate velocity are given by v i = dx i /dt e . We remind the reader that f α is defined according to Eq. (9). In addition, we make the observation that in the Newtonian limit of Eq. (13), we recover Newton's universal law of gravitation perturbed by an arbitrary non-gravitational force f i . Finally, it is interesting to note for an arbitrary non-affine parameter such as the proper time of an external Solar System probe or celestial body σ , the resulting equations of motion are implicitly defined according to where and so that Eq. (14) must be solved iteratively.

General relativistic invariant
We note the explicit expression for the final quantity (dτ/dt e ) arising in the right-hand side of Eq. (13) is determined by the differential relationship between the square of the magnitude of the invariant spacetime line element (ds) 2 and proper time namely and is given explicitly by We further note, the components of the spacecraft four-velocity are equivalently expressed in terms of ephemeris time according to Hence, on substituting Eq. (18) in (19), it can be readily shown that the particle tangent vector to the worldline satisfies the normalisation condition given by (Poisson and Will 2014) Equation (20) is an invariant of the motion and is defined implicitly in terms of numerically integrated spacetime coordinates arising in the metric tensor components and explicitly in terms of the spacecraft four-velocity vector (or, coordinate velocity if (19) is used) so that I ≡ I g μν (x α ) , u α . Eq. (20) remains valid when the spacecraft is subject to external perturbations and has no Newtonian analogue, giving rise to a unique opportunity for GRAPE to utilise symplectic integration schemes providing an accuracy check over the numerical integration span (see Sect. 4 for further details).

Tetrad formalism
GRAPE utilises the so-called tetrad formalism in order to relate dynamical quantities between local (known as the natural frame) and global frames. That is, at each point in spacetime, we can construct a pseudo-Cartesian reference system through the introduction of a tetrad 3 e μ (ν) so that (Brumberg 2017) where indices are manipulated using the components of the Minkowski metric η αβ and it is clear for a diagonal metric we have the relations In the natural tetrad, all dynamical quantities have physical meaning, so that for example, the expressions denote locally measured intervals of time and distance, respectively (Brumberg 2007). From Eqs. (24), it is clear that the transformations relating natural and global four-velocity components are given by the relations and the corresponding transformations for the non-gravitational forces are given, respectively, by where we have introduced D/Dτ ≡ u α ∇ α and ∇ α denotes the covariant derivative operator. Hence, the resulting expressions for the natural and global four-forces are given, respectively, by where we have used known as the tetrad postulate (Yepez 2011).

Non-gravitational force transformations
Finally, we describe the method adopted by GRAPE in order to determine the components of the non-gravitational forces in the global frame. The covariant derivative of Eq. (20) is by definition (Poisson and Will 2014) and using (9) simplifies to give the orthogonality constraint The spatial components of the external force are measured by on-board accelerometers or modelled by dynamical equations in the local co-moving frame (L) of the spacecraft and provided to GRAPE by end-users. We have, by definition, in the rest frame of the spacecraft u α L = (c, 0, 0, 0) so that in order to satisfy Eq. (30), we find f α L = (0, f i L ). Hence, there is no temporal component of the non-gravitational force in the co-moving frame of the spacecraft. At this point, we make the important observation that in the local rest frame of the spacecraft we have f α Thus, the components of the non-gravitational force density in the natural tetrad are obtained using a Lorentz transformation where the components of Λ α β v i are given by where the Lorentz factor is given by γ = dt e /dτ , the magnitude of the relative frame velocity is given by v and the inverse transformation is given by Λ β α −v i . Finally, the global components of the non-gravitational force are obtained using Eqs. (27).

General relativistic accelerometer equations
Space-based accelerometer technology or drag-free satellites (Lange 1964) have been successfully utilised in space missions across a number of fields including geodesy (Touboul et al. 1999), planetary science (Lucchesi and Iafolla 2006) and gravitational-wave physics (Rodrigues and Touboul 2003). For near-Earth spacecraft, the analytical models for nongravitational forces such as solar radiation pressure and atmospheric drag depend on the physical characteristics of individual spacecraft (for example the area-to-mass ratio and solar reflectivity parameters) and its attitude (for example the orientation of solar panels) which give rise to a number of technical issues in modern orbit determination and prediction software. Typically, the parameters associated with non-gravitational forces such as solar radiation and atmospheric drag coefficients (Vallado et al. 2006) are prescribed as solve-for parameters whose value minimises the residuals between measurement and predictions in the orbit determination process (Schutz et al. 2004). For the motion of spacecraft and interplanetary probes about central bodies other than the Earth, the accurate modelling of non-gravitational forces becomes increasingly difficult. In this section, we derive the so-called general relativistic accelerometer equations for spacecraft which are subject to non-gravitational forces and, using a post-Newtonian expansion of the resulting expressions, describe the approach adopted herein to determine the components of f α .
Consider the perturbed motion of a spacecraft subject to the gravitational field of an arbitrary central body with spacetime coordinates given by x α and equations of motion prescribed according to Eq. (9). Likewise, consider the motion of a Test Mass (TM) located inside the spacecraft at x α T M ≡ x α + δx α , shielded from all non-gravitational forces so that the motion of the test mass follows a geodesic i.e. f α T M ≡ 0 in Eq. (9). Therefore, the relative acceleration between the spacecraft and test mass is given by where we have introduced ξ α ≡ x α T M − x α as the relative spacetime separation vector and we assume ξ 0 = 0. We note in the absence of non-gravitational forces f α ≡ 0, Eq. (33) reduces to the so-called equation of geodesic deviation (d'Inverno 1992).
The IAU-recommended post-Newtonian metric tensor components g μν are given by (Soffel et al. 2003) where W is a scalar gravitational potential defined in the Geocentric Celestial Reference Frame which generalises the definition of the classical Newtonian potential where the reader is referred to Soffel and Han (2019); Soffel (2000) for further details. In general, there exists an additional general relativistic vector potential W j in the off-diagonal time-space metric tensor components which accounts for the rotation of the central body. However, for the current analyses we assume a static spacetime so that W j ≡ 0. The invariant spacetime line element associated with Eqs. (34) is given explicitly by where it is standard practice to denote spacetime coordinates in the local geocentric frame using upper case symbols. Hence, the differential relations between proper and coordinate time are given (to order O(c −2 )) by which gives rise to the first-and second-order differential operators and respectively, where we have introduced Σ ≈ 1 + 2W /c 2 + V 2 /c 2 . Hence, the spatial components of the post-Newtonian accelerometer equation are given by where The nonzero Christoffel symbols associated with the first post-Newtonian metric tensor components (34) are given by where the reader is reminded that in the present manuscript we assume a static spacetime so that the so-called gravitomagnetic vector potential W j ≡ 0 (cf. O'Leary et al. (2018) for the explicit derivation of the Christoffel symbols when W j = 0). The above expressions are obtained following direct application of the first-and second-order differential operators (37) and (38) to Eq. (33). The classical accelerometer equation is deduced in the weak field Newtonian limit (d'Inverno 1992) of expressions (40) to (43) and is given by where, at the measurement accuracy of space-based accelerometers, the tidal effects contained in the expression ∂ j ∂ i U ξ j are negligible so that Eq. (45) is approximately given by Hence, in the Newtonian limit, the relative acceleration between the test mass and the spacecraft is determined by the non-gravitational forces K i only.

GRAPE: structure-preserving integration scheme
In the framework of general relativity, the equations of motion for test particles are given by a system of four nonlinear, second-order, coupled differential equations. Although exact analytical expressions exist for the geodesic motion of test particles in Schwarzschild and Kerr geometries (Misner et al. 1973), the system of equations (9) are, in general, non-integrable and require numerical integration. Classical numerical procedures for the solution of differential equations (Butcher and Goodwin 2008) such as the family of explicit, single-step Runge-Kutta schemes do not preserve first integrals or phase-space volume (or, more generally, the Poincare integral invariants) and, as a result, are unsuitable for problems in orbitography or celestial mechanics requiring long-term stability. Symplectic integration schemes (Hairer et al. 2006) are advanced numerical procedures which preserve the symplectic structure of Hamiltonian systems so that the so-called flow of the system defines a canonical transformation (Yoshida 1990). Symplectic integration algorithms belong to a class of numerical methods known as structure-preserving geometric integrators and have been successfully applied across a myriad of fields (Donnelly and Rogers 2005;Hairer et al. 2006) where the reader is referred to Kinoshita et al. (1991) for a detailed discussion on the merits of symplectic integration algorithms to problems in astronomy.

GRAPE: Gauss collocation procedure
In the celebrated work of Butcher (1964), the author discusses the properties of implicit Runge-Kutta schemes and presents the Butcher tableau coefficients for a series of s-stage, n th -order algorithms known collectively as Gauss collocation methods (Hairer et al. 2006). The coefficients for the 5-stage, 10 th -order integration scheme implemented in GRAPE are given by where the constants α, β, γ are defined by α = 1/2, β = 32/225, γ = 64/225 and the m = 1, 2, . . . , 7 components ω m , ω m are defined according to (Butcher 1964) where δ = √ 70. Hence, given an arbitrary system of N differential equations the solution x n → x n+1 with integration step size h is given by We make the observation that Eqs. (49) are defined implicitly. At present, GRAPE determines the individual k i j using fixed-point iteration where the initial value is assumed zero everywhere so that (k 0 ) i j = 0. In future studies, improved strategies will be employed to initialise the k i j which requires further investigation (Hairer et al. 2006).

Quadratic invariants
Central to the present paper is the conservation of first integrals associated with the perturbed motion of a test particle described by Eq. (9). General relativity gives rise to a unique dynamical quantity (20) known formally as a quadratic invariant (Hairer et al. 2006). The coefficients b i , a i j appearing in Butcher (1964) which, as a result, indicates that Gauss collocation methods (Butcher 1964) are symplectic and preserve quadratic invariants completely (Cooper 1987;Butcher 2016). The proof is outlined below and the reader is referred to (Hairer et al. 2006, Chapter 4) for complete details.
A first integral I(x j ) associated with an arbitrary system of differential equations (47) satisfies the following condition We define a quadratic function of (47) according to where M i j is an arbitrary square, symmetric matrix. Hence, Q(x k ) is a first integral of (47) if the following holds (Hairer et al. 2006) From the above, it is clear that we must have dM i j /dt = 0 which is demonstrated as follows and using the symmetry M i j = M ji we find if and only if M i j is a constant.
According to Eq. (48) we have the following where it is important to note that Eq. (56) is obtained by squaring Eq. (48), x i n denotes the i th component of an arbitrary vector at time t = t n and k i m denotes the m th approximate slope associated with the i th differential equation. In order for Eq. (52) to be classified as a first integral we must have Q n+1 = Q n = ... = Q 0 = C, ∀n where Q n ≡ Q(x k n ) denotes the value of the quadratic function (52) at time t n and C is an arbitrary constant. Hence, using Eq. (49), we can equivalently express Eq. (56) according to where, for brevity, we have introduced k i j = f i T , X i n with T = t n + c j h and X i n = x i n + h s l=1 a jl k i l . Hence, according to Eqs. (53) and (50), the second and final terms on the right hand side of Eq. (57) are zero so that for Runge-Kutta methods satisfying (50), we have M i j x i n+1 x j n+1 = M i j x i n x j n . Thus, with respect to the current problem which, when expressed according to Eq. (52), we find Q(x k ) = I(η μν (x α ), u α ), where I is defined according to (20) and the condition (55) corresponds to the orthogonality condition (30) associated with the four-velocity and acceleration. It is important to note that in order to enforce strict symplecticity, we must express I in the natural tetrad so that at each numerical integration step the spacetime is assumed flat which is a formal requirement according to the condition dM i j /dt.

GRAPE: first results and discussion
In this section, we demonstrate the utility of the newly developed GRAPE software through a series of numerical simulations. We consider the motion of three test particles with orbital parameters assimilated to Molniya (M), Parker Solar Probe (PSP) and Mercury Planetary Orbiter (MPO)-like objects. We assume each test particle is initially located at apoapsis r a so that the initial spacetime coordinates and four-velocity are given by x α 0 = (0, r a , 0, 0) and u α 0 = (c, 0, v a cos i, v a sin i), respectively, where v a is the magnitude of the spacecraft velocity at apoapsis and the reader is referred to Table 1 for further details. The gravitational field of the respective central body is assumed to be modelled by the metric tensor components associated with the spherically symmetric Schwarzschild solution. While the choice of the underlying coordinate system is arbitrary, integration is performed in (Cartesian) isotropic coordinates which are easy to implement numerically.
For the benefit of the reader, and in order to facilitate future replications of the results presented in the current manuscript, we briefly outline the dynamical quantities used in the numerical simulations. The metric tensor components and corresponding invariant spacetime line element associated with the spherically symmetric Schwarzschild field of the central body are given by and respectively, where ρ s = r s /4, ρ = (x i x i ) 1/2 , the Schwarzschild radius is denoted r s and x i = (x, y, z) denote pseudo-Cartesian coordinates (Müller and Grave 2009).We transform dynamical quantities from global to natural frames using Eqs. (22) with the corresponding inverse transformations following trivially. The equations of motion used in the numerical simulation are given by where in the above equations j = i and we define w 1 = 1 − (ρ s /ρ) 2 , w 2 = (ρ + ρ s ) −7 , w 3 = 1 + ρ s /ρ and the Christoffel symbols of the second kind required to derive Eqs. (61) and (62) are given explicitly in Müller and Grave (2009). Finally, we note, in each simulation, we assume a constant radiation-like four-force (in the direction opposite to central body) acting on the spacecraft measured in the local frame (see Sect. 2 for implementation details) where f α ∼ 1 × 10 −6 m/s 2 . Again, it is important to note that the specific form of the external perturbation is selected arbitrarily and will be provided to GRAPE by end-users using either specific force models or data obtained from an accelerometer. It is interesting to consider the evolution of the individual spacecraft proper time when compared with a derived coordinate time as per Fig. 1. Taking into account the order of magnitude and the nonlinear evolution of δt, we can extract information regarding the orbit type. It is clear that Molniya spacecraft are in a highly eccentric orbit given the periodic oscillation occurring at perigee where special and general relativistic time dilation effects are greatest in magnitude. While the order of magnitude of the combined effects is low for the weak gravitational field of the Earth, we observe that for the PSP, the effects are increased by approximately six orders of magnitude due to the Schwarzschild field of the Sun and the high orbital velocity at perihelion. Complete conservation (at the order 10 −32 ) of the quadratic invariant (20) is achieved for the perturbed motion of the three test particles described in Table 1 and presented in Fig. 2. We simulate three orbits per test particle, selecting different integration time steps h which are typical in magnitude for practical operations. For example, consider the simulated orbit of the perturbed NASA Parker Solar Probe-like object. We assume a time step of approximately five  Fig. 2) where the typical order of magnitude used for generating ephemerides for probes during interplanetary cruise is approximately ten minutes. The standard "Newton + correction" approach for the numerical propagation of spacecraft is in conflict with the resolutions of the IAU which necessitate that all astrodynamical problems be formulated in the framework of Einstein's general theory of relativity. Given the rapid increase in measurement technology, we reiterate the important point: modern orbitography software is approaching its limits in terms of complexity. The standard approach to include relativity as an effective force presents both technical and conceptual problems in modern orbitography. It ignores the fundamental principles of relativity and gives rise to opportunities to duplicate or ignore relativistic effects. Further, as precision increases in measurement technology, propagation tools need to be continuously updated in order to account for additional relativistic effects and changes in conventions; increasing opportunities for human error.
We have presented a novel approach to propagate the motion of test particles subject to non-gravitational forces in the complete framework of general relativity which does not appear to have been exploited in the literature. GRAPE accounts for all relativistic contributions and all non-gravitational forces at the order of the chosen spacetime metric and the measurement sensitivity of the on-board accelerometer, respectively. Further, due to the existence of the unique general relativistic quadratic invariant (20), the integration procedure adopted in GRAPE provides an additional layer of security and an accuracy check during the propagation interval. We note that given the order of magnitude of the relativistic contributions arising in the PN metric tensor components, GRAPE employs 128-bit (quadruple) precision to perform calculations. While 64-bit (double) precision is less computationally demanding, we conjecture that high-order relativistic contributions will appear as noise at this level of precision.
It is of interest to the present authors to investigate and perform comparisons with alternative integration procedures in future implementations of GRAPE. We note the truncation error of integration procedures is of a mathematical origin, where the accumulation of roundoff errors is intrinsically linked with the machine architecture. Hence, the optimal time step required is machine-precision dependent, and is obtained as a trade-off between truncation and round-off errors. Ideally, the optimal time step must be determined using a reference integration scheme built with respect to a representation of floating numbers with a higher number of bits with respect to the method to be assessed in order to reduce the numerical round-off noise to nearly zero (Balmino and Barriot 1990). Given that GRAPE is currently utilising 128-bit precision necessitates a reference integrator with 256-bit (octuple) precision which is rarely used in practice. A typical reference integrator is the classical Bulirsch and Stoer algorithm with Richardson extrapolation (see Stoer and Bulirsch (2013) and Press et al. (2007) for implementation details). More recently, there has been renewed interest in utilising high-order, adaptive integration procedures based on Taylor's method (Biscani and and Parker Solar Probe (bottom)-like objects. The individual orbits are propagated with integration time step h which is calculated with respect to the orbital period T . We define δI according to δI ≡ I − c 2 /c 2 where I is calculated at each numerical integration step. Hence, for the perturbed motion of test particles subject to a locally measured radiation-like four-force, GRAPE preserves I at the order of 10 −32 Izzo 2021) which have been shown to be competitive with and in some cases superior to symplectic and non-symplectic schemes.
The coefficients of the implicit Runge-Kutta integration scheme currently implemented in GRAPE are analytically derived in Butcher (1964) up to order 10. While higher order methods allow for a larger time step h and a reduction in computational demand, their derivation is complex. The construction of higher order implicit Runge-Kutta methods which also satisfy the highly nonlinear constraint (50) is a non-trivial task and the reader is referred to Yoshida (1990) and Hairer et al. (2006) for further technical details which are beyond the scope of the current paper.
As a final remark, we outline an approach to extend GRAPE for practical mission planning and compatibility with NASA's SPICE kernels (Acton 1997). This will be formally addressed in a future manuscript and we expect the important results of the Cassini-Huygens (Bertotti et al. 2003) and BepiColumbo (Serra et al. 2018;Iess et al. 2021) radioscience experiments to serve as important verification tests for future versions of GRAPE. When complete, GRAPE will utilise the PN IAU operational-ready spacetime metrics for both geocentric and barycentric studies and will calculate important light-time observables ) associated with interplanetary probes. The complete n−body barycentric metric tensor components are given by (Damour et al. 1991;Soffel 2000) g 00 = 1 − 2w c 2 + 2w 2 c 4 , g 0 j = 4w j /c 3 , g jk = −δ jk 1 + 2w/c 2 , where, for consistency we have adopted the opposite metric signature to Soffel (2000) and w is defined in terms of the individual contributions due to monopole w 0 and higher order post-Newtonian multipole w L moments of the A = 1, . . . , N Solar System bodies and an additional post-Newtonian contribution so that and the vector potential w i accounts for the motion (rotational and translational) of the respective Solar System bodies. Again, we note that for consistency the above definitions (63) and (64) are defined using the opposite metric signature to Soffel (2000). At this point, we refrain from providing explicit definitions of the individual terms as they are dependent on the accuracy required for the specific mission and refer the reader to Soffel (2000) for further details. However, we make the observation that A is defined in terms of Solar System body position, velocity and acceleration which are provided to GRAPE via SPICE kernels. Hence, the resulting equations of motion will be defined implicitly and will need to be solved by iteration, which, in the monopole approximation, is equivalent to the EIH approach discussed in Standish et al. (1992). Finally, we note the analytical computation of the Christoffel symbols associated with (63) can be calculated using computer algebra software packages to avoid computationally demanding numerical differentiation.