Reconstructing the cruise-phase trajectory of deep-space probes in a general relativistic framework: An application to the Cassini gravitational wave experiment

Einstein’s theory of general relativity is playing an increasingly important role in fields such as interplanetary navigation, astrometry, and metrology. Modern spacecraft and interplanetary probe prediction and estimation platforms employ a perturbed Newtonian framework, supplemented with the Einstein-Infeld-Hoffmann n-body equations of motion. While time in Newtonian mechanics is formally universal, the accuracy of modern radiometric tracking systems necessitate linear corrections via increasingly complex and error-prone post-Newtonian techniques—to account for light deflection due to the solar system bodies. With flagship projects such as the ESA/JAXA BepiColombo mission now operating at unprecedented levels of accuracy, we believe the standard corrected Newtonian paradigm is approaching its limits in terms of complexity. In this paper, we employ a novel prototype software, General Relativistic Accelerometer-based Propagation Environment, to reconstruct the Cassini cruise-phase trajectory during its first gravitational wave experiment in a fully relativistic framework. The results presented herein agree with post-processed trajectory information obtained from NASA’s SPICE kernels at the order of centimetres.


Introduction
The accuracy of interplanetary navigation, planetary sciences, and solar system tests of general relativity has greatly benefited from modern improvements in multi-link telecommunication configurations and the use of highfrequency microwave signals between ground stations and interplanetary probes. Precise radiometric tracking data obtained during planetary flybys and dedicated body-centric missions have provided ample opportunity for scientific exploration of the mass, gravitational field, atmosphere, and orbits of celestial bodies; see Ref. [1] and references therein for an overview of notable space missions dedicated to planetary sciences.
The use of joint Ka/X band microwave links coupled with next-generation multi-frequency technologies [2] has yielded unprecedented accuracy for solar system tests of general relativity and facilitated the accurate determination of post-Newtonian parameters [3]. Of particular interest for weak-field tests of gravitation is the determination of the parameter γ (equal to one in general relativity), which measures the amount of space curvature produced by a massive body [3]. The deflection and time delay of photons owing to spacetime curvature can be expressed in terms of γ, which is currently constrained to within (2.1 ± 2.3) × 10 −5 of the predicted unity value [4]. By exploiting the multi-frequency transmission available at both ground station and spacecraft, the solar conjunction experiment performed by the Cassini probe in 2002 [5] provided the first complete solar plasma calibration for both uplink and donwlink signals [6]. The Cassini radioscience experiment improved the estimation of γ by approximately two orders of magnitude, which was previously constrained to within ∼10 −3 of the general relativistic prediction using data from the Viking Martian lander relativity experiment [7].
Through a series of resolutions, the International Astronomical Union (IAU) [8] suggests that all problems in the field of astronomy or astrodynamics should be formulated within the framework of Einstein's theory of general relativity. These include: (i) the derivation of dynamical equations of motion; (ii) the definition of observables; and (iii) the transformations between reference frames and time scales. The high precision range and Doppler observables combined with data from onboard accelerometers [9][10][11][12] necessitate that modern orbitography software employs high-fidelity dynamical and measurement models for spacecraft propagation and nonlinear state and parameter estimation, respectively [13,14].
Popular orbitography software such as the French Space Agency Géodésie par Intégrations Numériques Simultanées (GINS) [15] or NASA's Mission Analysis, Operations, and Navigation Toolkit Environment (MONTE) [16] currently describe the motion of interplanetary spacecraft via a classical Newtonian framework linearly corrected with effective forces, accounting for the effects of general relativity via the socalled n-body Einstein-Infeld-Hoffmann (EIH) equations of motion [17]. Given the stringent accuracy requirements associated with fields such as astrometry, geodesy, celestial mechanics, and metrological sciences [18], spacecraft and interplanetary probe orbit propagation tools based on the so-called "Newton + correction" framework need to include subtle relativistic effects in both state and measurement models to reflect the rapid improvements in modern measurement technologies. The aforementioned ad-hoc process of reducing relativistic phenomena to simple perturbations or relativistic corrections introduces a neo-Newtonian [19] interpretation of the theory of general relativity, which can lead to spurious results when equivalence between Newtonian and relativistic concepts is assumed. Moreover, the "Newton + correction" approach creates the likelihood of neglecting potentially important relativistic contributions if care is not taken. Indeed, it has been discussed at length in relevant literature [20][21][22] whether the dynamical equations and the corresponding measurement models used in the data reduction process of the celebrated Cassini solar-conjunction experiments [23] correctly account for the relative motion of the Sun with respect to the solar system's barycenter. The expected theoretical difference affects the obtained value of the post-Newtonian parameter γ = (2.1 ± 2.3) × 10 −5 by the amount ±1.2 × 10 −4 owing to the barycentric velocity of the Sun [2,21]. Nevertheless, the legacy software used in the Cassini radioscience experiment, i.e., Orbit Determination Program [13,14], employed an approximate model for the motion of the Sun, with a negligible effect on the final result [20].
It is expected that the current best estimate of γ will be improved in the near future by the Mercury Orbiter Radioscience Experiment (MORE) onboard the ESA/JAXA BepiColombo mission [24]. The radioscience instrumentation onboard BepiColombo will provide plasma-free range information in addition to stable Ka/X band microwave links [25], while simultaneously benefiting from the direct measurement of non-gravitational forces in the local frame of the spacecraft by an onboard accelerometer [26]. The novel 24 Mcps pseudo-noise ranging instrumentation onboard BepiColombo [27] will provide two-way range information at the centimeter level, so that secondorder post-Newtonian contributions in the light-time solution need to be considered in the data reduction process [28]. In addition, future space-based gravitational wave detectors, e.g., the Laser Interferometer Space Antenna (LISA), are expected to detect the propagation of low-frequency gravitational waves (10 −5 -1 Hz), complementing the scientific operations of groundbased detectors such as the Laser Interferometer Gravitational wave Observatory (LIGO), VIRGO, and KAGRA collaborations, which focus on higher frequencies in the gravitational wave frequency spectrum. LISA employs time-delay interferometry (TDI) to reduce the effects of laser frequency and optical bench noise [29]. TDI analysis employs solar system barycentric coordinate time while each LISA spacecraft records data according to the individual spacecraft proper time [30]. Hence, sophisticated four-dimensional general relativistic spacetime transformations are required for TDI data reduction to detect gravitational waves using space-based interferometry. Bearing this in mind, we believe that the classical "Newton + correction" paradigm is now reaching its limits in terms of complexity.
In a recent article [31], we presented the first results of a prototype software, namely General Relativistic Accelerometer-based Propagation Environment (GRAPE), which describes the perturbed motion of Parker Solar Probe [32], Mercury Planetary Orbiter [33], and Molniya-like [34] interplanetary probes and spacecraft subject to a radiation-like non-gravitational four-force using the complete framework of general relativity. GRAPE employs extended relativistic equations of motion that account for non-gravitational forces using either end user supplied accelerometer data or approximate dynamical models in the local frame of the spacecraft. We exploit the tetrad formalism [35] to compare dynamical quantities and non-gravitational forces in both co-moving and global frames. The novel approach adopted by GRAPE assumes that the local structure of spacetime is unperturbed by external non-gravitational forces, and that all associated gravitational and relativistic contributions are contained and accounted for at the order of the spacetime metric tensor. In this paper, we discuss the preliminary steps taken to extend GRAPE to the operational planning of future missions within the solar system. Further, we describe the initial approach adopted to interface it with NASA's SPICE kernels [36], and verify our results by comparing the pre-fit Doppler residuals using GRAPE and SPICE with high-precision Ka/Ka band Doppler data obtained from the first Cassini probe gravitational wave experiment (GWE1) [5]. The remainder of this paper is organized as follows. In Sections 2 and 3, we summarize the equations of motion and associated numerical integration procedures implemented in GRAPE, respectively. In Section 4, we present new results that extend GRAPE for practical mission planning using radiometric observables and determine that GRAPE recovers the Cassini trajectory at the order of centimetres. We conclude the paper in Section 5, providing a brief discussion on future plans for the platform.
Unless otherwise stated, we assume the Einstein summation convention for repeated indices; the use of Latin indices (i, j = 1, 2, 3) is reserved for spatial coordinates, and Greek indices (α, β = 0, 1, 2, 3) denote spacetime coordinates, with the 0th component reserved for time with x α = (ct, x i ). Finally, we assume a spacetime metric signature according to (+, −, −, −), a denotes the derivative of a with respect to time, and c denotes the speed of light in vacuum.

GRAPE: A brief overview
In this section, we discuss the mathematical framework implemented in GRAPE. First, we provide a high-level overview of the platform and discuss some aspects of the software that are unique to general relativity and not common with classical propagation environments. In Section 2.1, we present the equations of motion and discuss the treatment of non-gravitational perturbations in a fully relativistic framework.
In the framework of general relativity, the unperturbed motion of test particles is described using a geodesic equation parameterized by an arbitrary affine parameter s. For massive test particles such as interplanetary probes and near-Earth spacecraft, it is common to parameterize equations of motion using their associated proper time τ such that s ≡ τ . While proper time parameterization is convenient when working in a theoretical framework, it is inconvenient for practical purposes. Spacecraft operators, tracking stations, and orbit determination as well as data reduction software reference the state, error covariance, and observables with respect to a derived external coordinate time scale t such as terrestrial, barycentric, or ephemeris times. Hence, GRAPE internally adopts a global coordinate time t = x 0 /c during numerical propagation, and the corresponding equations of motion are transformed from an affine to a coordinate time parameterization accordingly.
The incorporation of non-gravitational forces in classical Newtonian propagation environments is a well-recognized problem and remains an active area of research in operational astrodynamics [37]. As discussed at length in Ref. [31], the treatment of non-gravitational perturbations in general relativity has not received the same attention in the literature and requires further elucidation. The nonlinear field equations of general relativity [38,39] can be decomposed according to contributions from the Einstein tensor G µν and the energy momentum tensor T µν . The Einstein tensor is composed of the metric tensor and quantities derived from it, e.g., the Ricci and Riemann tensors, which encode the geometry of spacetime so that all gravitational phenomena are captured entirely at the order of the spacetime metric tensor. Similarly, quantities derived from the energy-momentum tensor encode information related to non-gravitational phenomena. In the parlance of operational astrodynamics and interplanetary navigation, gravitational and nongravitational perturbations are obtained from the metric and energy-momentum tensors, respectively. Hence, the perturbed equations of motion employed by GRAPE are derived from the divergence of the energy-momentum tensor, which is linearly decomposed according to contributions from the energy-momentum of an isolated test particle and the corresponding stresses acting on the system that produce geodesic and perturbed motion, respectively [40,41].

GRAPE: Equations of motion
To account for non-gravitational forces in Einstein's theory of general relativity, we decompose the energymomentum conservation equation according to where we introduced the stress tensor S µν , and we interpret Θ µν as the energy-momentum of an isolated test particle. In its simplest form, Θ µν is given in terms of mass density ρ and four velocity u µ ≡ dx µ /dτ , with Θ µν ≡ ρu µ u ν . Equation (1) gives rise to two distinct cases: (i) S µν = 0 and (ii) S µν ̸ = 0. In the absence of stresses, it is well established that Eq. (1) yields a geodesic equation of motion. In the case where S µν ̸ = 0, we find where ∇ µ S µν = −ρK ν is the force density associated with an arbitrary non-gravitational four-force f ν . Evaluating the term on the left-hand side of Eq. (2) yields [31]: where it is clear that geodesic motion is recovered when K ν ≡ 0. As previously discussed, Eq. (3) can be equivalently expressed with respect to a non-affine coordinate time t following several applications of the chain rule: so that the final equations of motion employed by GRAPE are given by 1 where we introduced the dimensionless quantity β i = v i /c in terms of the test particle four-velocity v i = dx i /dt, and the Christoffel symbols of the second kind are denoted by Γ α βγ . It is worth noting that the components of the four-force in Eq. (5) are given by Eq. (3), and in the appropriate Newtonian limit, Eq. (5) reduces to Newton's second law perturbed by an external force f i , which can be most easily observed via the explicit calculation of the Christoffel symbols given by Eq. (2) in Ref. [42].
On close inspection, Eq. (3) raises two important practical questions: (i) How are the components of K µ obtained? (ii) Given that non-gravitational forces are measured in the local frame of the spacecraft, how does GRAPE relate the quantities between local and global frames, such as the IAU-recommended [8] geocentric and barycentric celestial reference frames? We address both questions below.
The inner product of two four-vectors is a scalar invariant I. In particular, it can be readily shown that in all reference frames, Differentiating Eq. (6) with respect to proper time and substituting the expression given by Eq. (3), we find that at every point in spacetime, the four-force is orthogonal to the four-velocity with For Eq. (7) to be satisfied in the local (L) co-moving frame of the spacecraft where u µ L = (c, 0, 0, 0), we find that f α L = 0, K i . Hence, in the local co-moving frame, the non-gravitational perturbations are purely spatial and are provided to GRAPE by end users using either known dynamical models or data obtained from an onboard accelerometer. The three-dimensional (Newtonian) nongravitational forces acting on interplanetary probes used in classical orbitography are well known and documented; see Ref. [43] for explicit examples. However, in the framework of general relativity, the corresponding models for the impulsive and continuous four-forces are not readily available, posing a theoretical obstacle for general relativistic orbitography. By modeling the non-gravitational four-forces in the local frame of the spacecraft first, we circumvent this issue, whereupon the global four-forces can be obtained using the so-called tetrad formalism [35].
To compare dynamical quantities between local and global frames, GRAPE employs the tetrad formalism as follows. At every point in spacetime, we construct a pseudo-Cartesian reference system by introducing e µ (ν) which satisfy We note that e µ (ν) consists of four orthonormal vectors as opposed to a single tensor quantity known collectively as a tetrad or vierbein [44]. The vector of interest is indicated by the subscript in parentheses and the component of that vector is denoted by the superscript.
An important feature of the tetrad formalism is that quantities projected onto the tetrad have physical meaning. Therefore, denote locally measured spatial and temporal intervals [35]. From Eq. (9), it is clear that the force components are obtained according to Equation (10) provides a direct transformation between local and global force components when e µ (ν) denotes the components of the co-moving tetrad. If, for example, e µ (ν) denoted the components of the so-called natural tetrad, then a further Lorentz transformation is required to account for the relative motion (see Ref. [

GRAPE: Symplectic integration scheme
In general, classical numerical integration algorithms do not preserve the first integrals (constants of motion) associated with dynamical systems; thus, over long integration timescales, energy dissipation occurs, leading to incorrect qualitative dynamical behavior. Symplectic integration schemes [45] are advanced numerical procedures that preserve the symplectic structure of Hamiltonian systems such that the so-called flow of the system defines a canonical transformation and belongs to a class of numerical methods known as structure-preserving geometric integrators. In the celebrated work of Butcher [46], the author introduces and discusses the properties of implicit Runge-Kutta schemes, presenting the Butcher tableau coefficients for a series of s-stage, nth-order algorithms, known collectively as Gauss collocation methods [45]. Implicit Runge-Kutta schemes are symplectic by definition (see Section 3.2 for further details). Hence, given a dynamical system that admits first integrals or quadratic invariant quantities, symplectic integration procedures ensure that these remain constant (or bounded in some cases) throughout the numerical integration.
The dynamical systems associated with classical orbitography software are both time-and velocitydependent (e.g., atmospheric drag perturbations), so that the resulting system of equations are not readily integrable. Hence, in a classical framework, structurepreserving integration schemes do not provide an advantage over non-symplectic algorithms. This is not the case for general relativistic orbitography. Given that Eq. (6) is valid in every reference frame, regardless of time-or velocity-dependent external non-gravitational four-forces, a unique opportunity to exploit structurepreserving integration schemes is available with general relativistic mission planning platforms. GRAPE employs the entire suite of numerical integration procedures presented in Ref. [46], where the reader is referred to Refs. [31,46,47] for the specific Butcher tableau coefficients required for practical implementation. In this section, we review implicit Runge-Kutta algorithms and present a derivation of their structure-preserving properties for quadratic invariants. Finally, we note that this section is based on Ref. [31,Section 4] and remains largely unchanged from its original presentation.

GRAPE: Implicit Runge-Kutta procedures
Implicit Runge-Kutta algorithms proceed as follows. Given an arbitrary system of N differential equations with associated initial conditions: the solution is propagated x n → x n+1 with an integration step-size h according to where Eq. (13) is defined implicitly. Presently, GRAPE determines the individual k i j using fixed-point iteration, where the initial value is assumed to be zero everywhere so that (k 0 ) i j = 0. Improved strategies will be employed in future work to initialize k i j , which requires further investigation [45].

Quadratic invariants
General relativity gives rise to a unique dynamical quantity (see Eq. (6)) formally known as a quadratic invariant [45]. Thus, the symplectic integration procedures implemented in GRAPE provide an additional layer of verification during the propagation interval. The coefficients b i and a ij identified in Ref. [46] satisfy the constraint: which consequently indicates that Gauss collocation methods [46] are symplectic and completely preserve quadratic invariants. The proof is outlined below, and the reader is referred to Ref. [45,Chapter 4] for the complete details.
A first integral I(x j ) associated with an arbitrary system of differential equations (11) satisfies the condition in Eq. (15): We define a quadratic function of Eq. (11) according to where M ij is an arbitrary symmetric matrix. Hence, Q(x k ) is the first integral of Eq. (11) if Eq. (17) holds [45]: where it is clear that we must have dM ij /dt = 0 which is demonstrated as Eq. (18): and using the symmetry M ij = M ji we find if and only if M ij is a constant. According to Eq. (12), we obtain Eq. (20): where Eq. (20) is obtained by squaring Eq. (12). Furthermore, x i n denotes the ith component of an arbitrary vector at time t = t n , and k i m denotes the mth approximate slope associated with the ith differential equation. For Eq. (16) 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 (16) at time t n and C is an arbitrary constant. Hence, using Eq. (13), we can equivalently express Eq. (20) 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 . According to Eqs. (17) and (14), the second and final terms on the right-hand side of Eq. (21) are zero; therefore, for Runge-Kutta methods satisfying Eq. (14), Thus, with respect to the current problem, we find Q(x α ) = I(η µν (x α ), u α ), where I is defined according to Eq. (6), and the condition given by Eq. (19) corresponds to the orthogonality condition associated with the four-velocity and four acceleration; see Eq. (7). It is important to note that in order to enforce strict symplecticity, we must express I in the local comoving tetrad so that at each numerical integration step, the spacetime is assumed to be flat, which is a formal requirement according to the condition dM ij /dt = 0.
To conclude this section, we reiterate the important point: In classical orbit problems, symplectic integration procedures have been shown to have attractive numerical properties for the long-term numerical integration of Hamiltonian systems [45]. Implicit Runge-Kutta schemes [46] are not restricted to Hamiltonian systems. They are a family of implicitly defined numerical procedures for the numerical solution of any system of differential equations, e.g., stiff and non-stiff problems [45]. The Gauss collocation methods introduced above are a particular example of an implicit Runge-Kutta scheme with the added property that their Butcher tableau coefficients satisfy the constraint given by Eq. (14). Hence, they naturally preserve quadratic invariants, as discussed above.
In many classical problems, the introduction of velocity-and time-dependent perturbations results in a system of differential equations that may not readily admit first integrals, such that the structure-preserving properties of symplectic integration schemes can no longer be exploited. The equivalent problem in the framework of general relativity is different; regardless of the specific non-gravitational four-forces entering Eq. (5), the inner product of two four vectors will always be an invariant. In particular, Eq. (6) is always satisfied, despite the motion of the spacecraft or probe being subject to external forces. Accordingly, we choose to employ an implicit Runge-Kutta scheme for the present analysis, which ensures that Eq. (6) does not diverge from c 2 beyond machine precision, providing an additional accuracy check on modifications to the platform. While the results presented in Section 4 are restricted by the duration of GWE1, i.e., O (days), we note that GRAPE can be used for long-term integration, e.g., generating solar system ephemerides [O (centuries)] [48,49], which is of interest to the pulsar timing community for the detection of gravitational waves using time of arrival residuals [50]. The necessary modifications to GRAPE are minor and will be explored in a future study.

Results
In this section, we demonstrate the modifications made to GRAPE since Ref. [31] in preparation for future practical mission planning, which we verify through three independent strategies: (i) using GRAPE, we generate an ephemeris for the Cassini probe for the duration of GWE1 and ensure that Eq. (6)  Doppler ④ ; and (iii) using the SPICE spkezr ⑤ program, we perform a direct comparison of the GRAPE ephemeris data with the Cassini SPICE kernel and project the residuals onto the radial, along-, and cross-track axes in the local time-dependent spacecraft coordinate system; see Section 3.3.3 in Ref. [37] for further details. In the present analysis, the non-gravitational perturbations entering Eq. (23) are inferred directly from the Cassini SPICE kernel, and we assume that they are held constant throughout the experiment. Given that this is the first application of GRAPE to real radiometric tracking data, we persevere with this assumption and note that it is an approximation. We further note that the GRAPE estimation module is currently being developed, where specific time-dependent non-gravitational force models will be applied to parameter estimation problems for probes during interplanetary cruise. The explicit details and steps required to reproduce the present analysis are discussed below. Finally, we note that by construction, the routines utilized internally by the SPICE toolkit automatically take into account the Earth ephemeris, the location of Earth-based tracking stations, and the Earth's rotation. For approximately 40 days starting November 26, 2001, the Cassini probe was tracked during solar opposition in search of gravitational waves in the millihertz frequency band [52]. Given that gravitational wave detection methods depend on the timescale of the radiation, ground-based detectors based on laser interferometry cannot separate signals from seismic and other noise sources in the frequency range 10 −5 -10 −2 Hz. Consequently, searches for low-frequency gravitational waves must employ space-based detectors [5]. The gravitational wave experiments performed by Cassini were, at the time, the most sensitive Doppler tracking campaigns achieved using interplanetary probes. Moreover, at a distance of approximately ∼8 AU, Cassini became the largest gravitational wave antenna ever used [5]. The success of the Cassini gravitational wave experiments and solar system tests of general relativity was mediated by two major technical upgrades. First, Cassini was equipped with a novel multi-link telecommunication configuration with two uplink carriers in the X and Ka bands as well as three downlink carriers in the X, Ka1 (coherent with the Xuplink), and Ka2 (coherent with the Ka-uplink) bands, considerably reducing the electromagnetic wave phase scintillation owing to solar plasma [6,52]. Second, the radioscience experiments performed by Cassini benefited from a dedicated advanced media calibration system located at the Goldstone Deep Space Communications Complex (DSS-25). The novel system consisted of water vapor radiometers, pressure sensors, and microwave temperature profilers that calibrated the frequency modulation caused by the wet and dry components of the troposphere [5,6,52]. The curved spacetime associated with the celestial bodies is described (up to order c −5 in the time component, up to order c −3 in the spatial components, and up to order c −4 in the off-diagonal components) by the components of the n-body metric tensor (see Ref. [14, Section 2]): where the gravitational parameter of the jth celestial body is given by µ j ; the kth component of its coordinate velocity is given byẋ k j ;ṡ 2 j = ẋ k j 2 ; β and γ denote post-Newtonian parameters; the relative barycentric distance between bodies i and j is given by r ij ; and we have introduced the Kronecker delta function according to δ kl . In the present demonstration, we assume that body i refers to the Cassini probe and bodies j = 1, . . . , n − 1 (j ̸ = i) refer to the Sun, the nine planets, and the Moon, where the initial conditions of all n bodies are obtained from the appropriate SPICE kernels on Nov 26, 2001 at 05:00:00 Barycentric Dynamical Time using the SPICE spkezr program. The extended equations of motion implemented in GRAPE follow directly from the definition of the metric tensor components and Eq. (5), and are given by Eq.  in Ref. [14,Section 4]: where we introduced the notation x i x i = x · x to denote the inner product between two vector quantities, and we assume that the components of x j are given in terms of Cartesian coordinates. Equation (23) is defined implicitly in terms ofẍ j , such that the formal numerical solution is obtained by iteration. However, for the present demonstration, we assume that terms of the formẍ i /c 2 can be replaced by their Newtonian counterparts. It is important to note that the final term appearing in Eq. (23) completely captures the general relativistic contributions of an arbitrary non-gravitational four-force f ν . The calculation of observable quantities in a complete general relativistic framework [53] relies on highly accurate light-time solutions and is profoundly different from classical signal-processing paradigms [54,Chapter 2]. In the present study, light-time is calculated using the approach presented in Ref. [14,Chapter 8]. We note that light-time solutions obtained by numerically solving the null geodesic equations where the effects of solar plasma are contained within a modified metric tensor [55] are currently being explored. Figure 1 displays the temporal evolution of the quadratic invariant I given by Eq. (16). The quadratic invariant I depends explicitly on the components of the numerically integrated four-velocity u α and implicitly on the numerically integrated spacetime coordinates x α through the metric tensor components g µν given by Eq. (22). At each integration step, GRAPE internally checks whether I deviates from machine precision. This is employed as an additional accuracy check on modifications made to GRAPE. It is clear from Fig. 1 that the numerically determined four-velocity components satisfy Eq. (6) at machine precision, which can be further reduced to 10 −32 by selecting 128-bit (quadruple) precision internally at a computational cost. The computational demand for double (64-bit precision) and quadruple precision is O (seconds) and O (minutes), respectively. We note that a thorough analysis benchmarking the GRAPE platform is deferred for future study and is outside the scope of the present paper. Figure 2 compares the pre-fit Doppler residual using GRAPE and Cassini SPICE trajectory kernels, as discussed in Steps (4)- (6). The details of the Cassini kernel and observational data used in the present analysis are presented in Footnotes ② and ④. For clarity, we plotted only the residuals for the first pass of Cassini during GWE1. However, we note that the difference in the residuals remains at the ∼0.2 Hz level, with no noticeable secular variation during the first 1 × 10 5 observations. The pre-fit Doppler residuals are slightly biased owing to uncertainties in the initial state of the spacecraft obtained from SPICE. The bias is corrected in the estimation phase, e.g., by employing a (sequential) batch least-squares or Kalman filter algorithm [37]. The initial-state bias in Fig. 2 is standard and is comparable to that of other modern orbit prediction and determination software; see Fig. 1 in Ref. [1] for example. Typically, in orbit determination problems, the goal is to obtain Gaussian residuals with zero mean. In this regard, the results presented herein are slightly different. In its simplest form, orbit determination software is typically decomposed into two modules: prediction and estimation.
The goal of the present study is not estimation via leastsquares or Kalman filter analysis. Instead, this study aims to present the first results of the GRAPE platform during the prediction phase. As noted above, the estimation module of GRAPE is currently under development and will be the topic of a future paper on general relativistic state and parameter estimation. Figure 3 compares the GRAPE and Cassini ephemerides in the time-dependent (R,Ŝ,Ŵ ) spacecraft frame, as discussed in Step (7). Using GRAPE, the reconstructed Cassini trajectory during cruise agrees with the SPICE ephemeris at the level of 10 −4 km. The projected residuals in the cross-and along-track directions, respectively denoted as ∆W and ∆S, are negligible. The projected residual along the radial direction, denoted as ∆R, grows secularly over the duration of the Cassini GWE1. During cruise, the dominant non-gravitational force is in the radial direction. Given that we assume the non-gravitational forces remain constant during GWE1, the projected residual observed along the ∆R axis is expected. We conjecture that once specific time-dependent non-gravitational four-force models are implemented in GRAPE, the ∆R residual will be significantly reduced ⑥ .
It is important that the results presented in this manuscript are easily reproducible. As a result, and in anticipation of a future public release of the GRAPE software, we outline the specific steps required to reproduce Figs. 1-3: (1) Generate an ephemeris using GRAPE by numerically solving Eq. (23) subject to the condition that f α ≡ 0. This step produces an ephemeris for the Cassini trajectory during the first gravitational wave experiment under the assumption that the motion of the probe is governed entirely by the n-body solar system metric tensor components given by Eq. (22). (2) Using the SPICE program mkspk, create a SPICE kernel from the GRAPE-generated ephemeris obtained in Step (1) The GRAPE (G) and SPICE (S) acceleration vectors are respectively defined by where we have introduced the acceleration owing to the n-body metric tensor components only as a M . Under the assumption that a NG ≈ a 0 during the short duration cruise phase (< 40 days), where a 0 denotes a constant vector, we employ Eq. (10) at each integration time-step and produce a new GRAPE-generated kernel to be used with the SPICE software. The ephemerides associated with the new kernel now take into account the non-gravitational forces acting on the probe under the assumption that they remain constant throughout the integration procedure. Again, we reiterate that the adopted approach is an approximation. The implementation of specific time-dependent models in the local tetrad is postponed to future work.
We project the residual position vector between GRAPE and SPICE, denoted by δr ≡ r G −r S , along each of the axes in the local spacecraft frame so that the residuals plotted in Fig. 3

Discussions and conclusions
We presented the first practical results of a novel general relativistic orbitography platform, GRAPE. While the initial results are in good agreement with NASA's SPICE software, further work is required. It is of utmost importance that GRAPE can accurately provide light-time solutions for the calculation of radiometric observables in a fully relativistic framework. Rather than adopting the standard approach in the literature, where light-time is approximated by a series of linearly superimposed post-Newtonian corrections that gradually increase in complexity and are error-prone for higherorder corrections, we are currently implementing lighttime solutions in GRAPE by numerically solving the null geodesic equation given in Ref. [55]. The latter approach guarantees that all relevant relativistic contributions to the dynamics are appropriately accounted for at the order of the metric tensor and voids the possibility of neglecting potentially important terms [20][21][22]. The motivation for designing a general relativistic orbit prediction and estimation tool is not to replace the current popular orbitography platforms based on the corrected Newtonian approach. Rather, it is to provide a fully consistent and highly accurate benchmark for current and future platforms. With radioscience data preprocessing software now publicly available (see for example Ref. [56]), archival tracking data from a myriad of space missions can be easily reduced, thereby providing ample opportunity to revisit historical missions and verify future modifications to the tool. We expect the current manuscript to be the first in a series of papers describing new approaches to general relativistic orbit determination (see for example Ref. [47,Chapter 7]) applicable to missions requiring high precision, e.g., BepiColombo [12], JUICE [57], and VERITAS [58,59] missions, as well as the space-based gravitational wave observatory LISA [60]. As a concluding remark, we note that GRAPE is currently tailored to process cruise-phase data only. Future modifications will extend the platform to planetocentric missions, which will benefit from the extensive treatment of relativistic reference frames and their associated spacetime transformations, time scales, and relativistic multipole structures readily available in the vast and extensive literature (see Refs. [8,61] for further details).
Note: Some findings of this study were presented at the 28th International Symposium on Space Flight Dynamics (ISSFD-2022-036) [62]. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.
The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.