General relativistic observables of the GRAIL mission

This paper describes a testable amendment to the general theory of relativity. The actual abstract contains two equations, so it cannot be properly displayed here per Academia.edu interface limitations. Download the paper to read the abstract.

One of the most famous classical tests of General Relativity is the gravitoelectric secular advance of the pericenter of a test body in the gravitational field of a central mass. To day, this effect has been detected, in an astrophysical context, with the well known measurements of the binary pulsars, and, in the Solar System, by measuring the perihelion shift of Mercury and other planets in the gravitational field of the Sun with the radar ranging technique. The claimed relative accuracy of 2 × 10 of the pulsars measurements is still object of debate and, in fact, according to some people, it should amount, instead, to almost 10. The relative accuracy in measuring the perihelion shift of Mercury is of the order of 10 or better. However, some authors believe that the measurement of the perihelion shift of Mercury may be affected by the uncertainty in the value of the quadrupole mass moment of the Sun. In this paper we explore the possibility of performing a complementary measurement...

Given the extreme accuracy of modern space science, a precise relativistic modeling of observations is required. In particular, it is important to describe properly light propagation through the Solar System. For two decades, several modeling efforts based on the solution of the null geodesic equations have been proposed but they are mainly valid only for the first order Post-Newtonian approximation. However, with the increasing precision of ongoing space missions as Gaia, GAME, BepiColombo, JUNO or JUICE, we know that some corrections up to the second order have to be taken into account for future experiments. We present a procedure to compute the relativistic coordinate time delay, Doppler and astrometric observables avoiding the integration of the null geodesic equation. This is possible using the Time Transfer Function formalism, a powerful tool providing key quantities such as the time of flight of a light signal between two point-events and the tangent vector to its null-geode...

General relativistic observables of the GRAIL mission Slava G. Turyshev1,3 , Viktor T. Toth2 , and Mikhail V. Sazhin3 1 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109-0899, USA 2 arXiv:1212.0232v4 [gr-qc] 18 Dec 2012 3 Ottawa, ON K1N 9H5, Canada and Sternberg Astronomical Institute, Lomonosov Moscow State University, Moscow, Russia (Dated: December 19, 2012) We present a realization of astronomical relativistic reference frames in the Solar System and its application to the GRAIL mission. We model the necessary spacetime coordinate transformations for light-trip time computations and address some practical aspects of the implementation of the resulting model. We develop all the relevant relativistic coordinate transformations that are needed to describe the motion of the GRAIL spacecraft and to compute all observable quantities. We take into account major relativistic effects contributing to the dual one-way range observable, which is derived from one-way signal travel times between the two GRAIL spacecraft. We develop a general relativistic model for this fundamental observable of GRAIL, accurate to 1 µm. We develop and present a relativistic model for another key observable of this experiment, the dual one-way range-rate, accurate to 1 µm/s. The presented formulation justifies the basic assumptions behind the design of the GRAIL mission. It may also be used to further improve the already impressive results of this lunar gravity recovery experiment after the mission is complete. Finally, we present transformation rules for frequencies and gravitational potentials and their application to GRAIL. I. INTRODUCTION Several past, present and planned space missions utilize a pair of spacecraft orbiting a celestial body in a tight formation. Continuous high-precision range and range-rate measurements between the spacecraft yield detailed information about the gravity field of the target body. Missions of this type include the Gravity Recovery and Climate Experiment (GRACE) mission [1] in orbit around the Earth; the Gravity Recovery and Interior Laboratory (GRAIL) mission, which comprises two spacecraft in orbit around the Moon [2–5]; and planned missions such as a GRACE Follow-on mission or a proposal for a GRAIL-like mission in orbit around Mars. Of these, the mission of particular current interest is GRAIL, as the two GRAIL spacecraft are presently (2012) orbiting the Moon. In this paper, we therefore focus on the GRAIL mission and its science observables. However, the lessons learned are also applicable to other, similar experiments. To reach its science objectives, the GRAIL mission relies on precision navigation of both spacecraft and accurate range measurements between the two lunar orbiters performed with their on-board Ka-band ranging (KBR) system. The instantaneous one-way range measurements performed at each spacecraft are time-tagged and processed on the ground to form dual one-way range (DOWR) measurements [6]. The mission relies on precision timing of all critical events (using the on-board ultra-stable oscillator, or USO) related to the transmission and reception of various microwave signals used on GRAIL for formation tracking and navigation. The resulting time series of highly accurate radio-metric data will allow for a major increase in accuracy when studying the gravity field of the Moon. The differential nature of the science measurements allows for the removal of a number of measurement errors introduced in the process. In particular, the approach compensates for errors due to long-term instabilities of the on-board USOs. This allows for an improvement in accuracy by about two orders of magnitude when compared to other techniques. In fact, the anticipated accuracies are of the order of 1 µm in range and 1 µm/s in range rate. It was recognized early on during the mission development that due to the expected high accuracy of ranging data on GRAIL, models of its observables must be formulated within the framework of Einstein’s general theory of relativity. In fact, a naive application of the observable models developed for the GRACE mission [1] may have led to significant model discrepancy (as emphasized in Ref. [6]), as these models do not take into account relativistic contributions that are critical for GRAIL. The ultimate observable model for GRAIL must correctly describe all the timing events occurring during the science operations of the mission, including both the navigation observables (Sand X-band, ∼ 2 GHz and ∼ 8 GHz correspondingly) and inter-spacecraft tracking (Ka-band, ∼ 32 GHz) data. The model must represent the different times at which the events are computed, involving the time of transmission of the Ka-band signal at one of the spacecraft, say GRAIL-A, at tA0 , and the reception of this signal by its twin, GRAIL-B, at time tB . In addition, the model must include a description of the process of transmitting S-band and X-band navigation signals from either spacecraft and reception of this signal at a Deep Space Network (DSN) tracking station at time tC . 2 FIG. 1: Representative geometry (not to scale) of the vectors involved in the computation of the GRAIL observables. “SSB” is the Solar System barycenter, “E” is the center of the Earth, “M” is the center of Moon, “EMB” is the Earth-Moon Barycenter. “A” and “B” are the positions of the GRAIL-A and GRAIL-B spacecraft, respectively, and “C” is the position of the DSN tracking antenna on the surface of the Earth. We model the range RAB = |RAB | between the two spacecraft A and B (see Fig. 1 for geometry and notations) as: RAB = |RAB | = |xB − xA | = |(xEM + xM + yB ) − (xEM + xM + yA )|, (1) where xEM is the vector connecting the Solar System barycenter (SSB) with the Earth-Moon barycenter (EMB), xM is the vector from the EMB to the Moon’s (M) center of mass, xA and xB are vectors connecting the SSB with the positions of the two GRAIL orbiters and vectors yA and yB connect the Moon’s center of mass with the orbiters. For navigation purposes, both orbiters maintain communication links with a ground-based DSN antenna. The range RAC = |RAC | between a GRAIL spacecraft (GRAIL-A, for instance) and a ground-based antenna can be modeled as: RAC = |RAC | = |xC − xA | = |(xEM + xE + yC ) − (xEM + xM + yA )|, (2) where xE is the vector from the EMB to the geocenter (E), xC is the vector connecting the SSB with the ground antenna whereas the vector yC determines the geocentric position of the ground antenna’s reference point. For actual computations, we use several different reference systems1 . The Solar System Barycentric Coordinate Reference System (BCRS) has its origin at the SSB. The origin of the Geocentric Coordinate Reference System (GCRS) is the Earth’s center of mass. Positions of DSN ground stations are given with respect to another terrestrial coordinate system, the Topocentric Coordinate Reference System (TCRS; see also Ref. [8]). We also consider the Lunicentric Coordinate Reference System (LCRS; for additional discussion, see Ref. [9]), the origin of which is fixed at the Moon’s center of mass. Finally, we attach to each spacecraft its Satellite Coordinate Reference System (SCRS; for a similar approach aimed to construct a reference frame for the GAIA project, see Ref. [10]). (We discuss these reference frames and their relationships in depth in Sec. II.) Equations (1) and (2) offer a good starting point to develop an appropriate relativistic formulation for the experiment. The six vectors involved in Eqs. (1)–(2) can be expressed in terms of their respective points of origin: e.g., xEM would be expressed in the BCRS, yA and yB in the LCRS, RAB and RAC in the SCRS of GRAIL-A, etc. Each of these coordinate systems has a corresponding time coordinate. To compute the vector sums and differences, all vectors involved must be converted to a common relativistic space-time reference system. Although in general relativity one can introduce any reference frame to describe the experiment, the best practical choice is offered by some realization 1 Following Refs. [7, 27], we use the term “reference system” to describe a purely mathematical construction, while a “reference frame” is a physical realization of such. 3 of the BCRS. We will use a realization of the BCRS that is called the SSB reference frame. The coordinate time associated with the BCRS is TCB (Barycentric Coordinate Time). For practical applications, it is often preferable to use another time scale, the TDB (Barycentric Dynamical Time). Currently published planetary ephemerides are provided using TDB. TDB and TCB differ only by a linear scaling. The advantage of using TDB is that the difference between it and terrestrial timescales (e.g., TT, defined in Sec. II D) is as small as possible and periodic. The choice of the TDB as the SSB time coordinate is realized by the appropriate linear scaling of space coordinates and planetary masses (see [11–13] for review). The vectors xE , xM , and xEM are readily available in the SSB reference frame, obtained by numerical integration and from Solar System ephemerides [14]. The vectors yA , yB and yC have to be transformed to the SSB frame from geocentric and lunicentric reference systems, respectively. Clearly, the required conversion between reference systems also involves conversion of the relativistic time coordinate. The equations of motion of the Moon and Earth, including all the relativistic effects at an accuracy even exceeding that of the GRAIL experiment, have already been discussed elsewhere [15]; here we concentrate on the computation of observables. In this paper, we focus on the formulation of a relativistic model for computing the observables of the GRAIL mission, with results that are applicable to other past and planned missions with similar observables. We address some practical aspects of the implementation of these computations. In Sec. II we discuss all relevant relativistic four-dimensional reference systems and the transformations that are required to make the vector sums in Eqs. (1) and (2) computable. In Sec. III we discuss the process of forming the inter-satellite Ka-band range (KBR) observables of GRAIL and derive a model for the dual one-way range (DOWR) observable. We also develop a relativistic model for another fundamental observable on GRAIL: the dual one-way range-rate (DOWRR). We conclude with a set of recommendations and an outlook in Sec. IV. In order to keep the main body of the paper focused, we chose to present some calculational details in the form of appendices. In Appendix A we present some important derivations: In Appendix A 1 we derive the solution for the post-Minkowskian space-time in general relativity, in Appendix A 2 we derive analytic expressions to describe the phase of an electromagnetic signal in gravitational field, and in Appendix A 3 we discuss the coordinate gravitational time delay. In Appendix B contains a discussion on the evaluation of the integral that is needed to assess the full accuracy of the DOWR observable. Finally, in Appendix C we briefly address the transfer of a precision frequency reference between the spacecraft and a ground station. The notational conventions used in this paper are as follows. Latin indices from the beginning of the alphabet, a, b, c, ..., are used to denote Solar System bodies. Latin indices from the second half of the alphabet (m, n, ...) are space-time indices that run from 0 to 3. Greek indices α, β, ... are spatial indices P that run from 1 to 3. In case of 3 repeated indices in products, the Einstein summation rule applies: e.g., am bm = m=0 am bm . Bold letters denote spatial (three-dimensional) vectors: e.g., a = (a1 , a2 , a3 ), b = (b1 , b2 , b3 ). The dot is used to indicate the Euclidean inner product of spatial vectors: e.g., (a ·b) = a1 b1 + a2 b2 + a3 b3 . Latin indices are raised and lowered using the metric gmn . The Minkowski (flat) space-time metric is given by γmn = diag(1, −1, −1, −1), so that γµν aµ bν = −(a · b). We use powers of the inverse of the speed of light, c−1 , and the gravitational constant, G as bookkeeping devices for order terms: in the low-velocity (v ≪ c), weak-field (GM/r ≪ c2 ) approximation, a quantity of O(c−2 ) ≃ O(G), for instance, has a magnitude comparable to v 2 /c2 or GM/c2 r. The notation O(ak , bℓ ) is used to indicate that the preceding expression is free of terms containing powers of a greater than or equal to k, and powers of b greater than or equal to ℓ. II. SPACE-TIME REFERENCE FRAMES AND TRANSFORMATIONS The theory of general relativity is generally covariant. In the Riemannian geometry that underlies the theory, coordinate charts are merely labels. One may choose an arbitrary coordinate system to describe the results of a particular experiment. Space-time coordinates have no direct physical meaning and it is essential to construct physical observables as coordinate-independent quantities. On the other hand, some of the available coordinate systems have important practical advantages. These systems are usually associated with a particular celestial body, ground-based facility or spacecraft, thereby yielding a material realization of a reference system to be used to describe the results of precision experiments. In order to interpret the results of observations or experiments, one picks a specific coordinate system that is chosen for the sake of convenience and calculational expediency, formulates a coordinate picture of the measurement procedure, and then derives the observable. It is also known that an ill-defined reference frame may lead to the appearance of non-physical terms that may significantly complicate the interpretation of the data. Therefore, in practical problems involving relativistic reference frames, choosing the right coordinate system with clearly understood properties is of paramount importance, even as we recognize that in principle, all (non-degenerate) coordinate systems are created equal [7]. In a recent study [15], we presented a new approach to investigate the dynamics of an isolated, gravitationally 4 bound astronomical N -body system in the weak field, slow-motion approximation of the general theory of relativity. Celestial bodies are described using an arbitrary energy-momentum tensor and assumed to possess any number of internal multipole moments. Using the harmonic gauge conditions together with a requirement for preserving conservation laws, we were able to construct the relativistic proper reference frame associated with a particular body. We also were able to determine explicitly all the terms of the resulting coordinate transformations and their inverses. In this paper we rely on the results obtained in Refs. [15, 16] and develop a set of coordinate reference frames for GRAIL. To reach its scientific objectives, in addition to the BCRS, GRAIL will have to utilize a set of several fundamental coordinate reference frames. These include terrestrial reference systems, namely the GCRS and the TCRS, and lunar reference systems, the LCRS and SCRS. In Ref. [15], we presented the detailed structure of the representations of the metric tensor corresponding to the various reference frames involved, the rules for transforming relativistic gravitational potentials, the coordinate transformations between the frames and the resulting relativistic equations of motion. The accuracy that is achievable by these calculations is sufficient to accommodate modern-day experiments in the Solar System and exceeds that needed for GRAIL. Here, we present the essential part of these transformations between various coordinate systems involved and dealing with transformations of relativistic time scales and position vectors, at the level of accuracy required by GRAIL. A. Barycentric Coordinate Reference Frame (BCRS) The Barycentric Celestial Reference System (BCRS) is defined with coordinates {xm } ≡ (ct, x = xα ), where t is TCB. The BCRS is a particular implementation of a barycentric reference system in the Solar System. The metric tensor gmn (x) of the BCRS satisfies the harmonic gauge condition. It can be written [15] as g00 = 1 − 2 2 w + 4 w2 + O(c−6 ), 2 c c g0α = − 4 γαλ wλ + O(c−5 ), c3 gαβ = γαβ + γαβ 2 w + O(c−4 ), c2 (3) where w and wλ are harmonic gauge potentials that can be presented, at the level of accuracy suitable for the purposes of the GRAIL mission (i.e., neglecting higher order mass- and current-multipole moments), in the form [7, 15, 16]: w = X GMb b w = rb X GMb b rb 1+ o 1 n 2 X GMc 2 1 1 2v − − (n · v ) − (r · a ) + O(c−3 ), b b b b 2 2 b c2 rcb (4) c6=b vb + O(c−2 ), (5) where rb = x − zb , rb = |rb |, and nb = rb /rb , with zb being the barycentric position of body b, and we use rab = rb − ra to denote the vector separating two bodies a and b. Also, the overdot denotes ordinary differentiation with respect to t, vb = żb (vb = |vb |) and ab = z̈b (ab = |ab |) are the barycentric velocity and acceleration of body b, and Mb is its rest mass. Lastly, the summation in (4)–(5) is being performed over all the bodies b = 1, 2..., N in the Solar System. The metric tensor (3) and the gravitational potentials (4)–(5) have sufficient accuracy for modern precision experiments in the Solar System [7, 15]. From Fig. 1, we can read off the barycentric positions of the Earth, zE , and the Moon, zM : zE = xEM + xE and zM = xEM + xM , respectively. Both of these vectors, and corresponding velocities vE = żE and zM = żM can be computed in the first post-Newtonian approximation using the Einstein-Infeld-Hoffmann (EIH) equations in the coordinates of the BCRS [15, 17–20]: z̈a = o X GMb rab n X GMc X GMc 1 1+ 2 −4 − + ṙa2 + 2ṙb2 − 4(ṙa · ṙb ) − 23 (nab · ṙb )2 + 12 (rab · r̈b ) + 3 rab c rac rbc b6=a c6=a c6=b X GMb r̈b 1 X GMb rab · (4ṙa − 3ṙb ) ṙab + 27 + O(c−4 ), (6) + 2 3 c rab rab b6=a b6=a where rab = |rab | and nab = rab /rab . When describing the motion of spacecraft in the Solar System, the models also include forces of attraction between the zonal harmonics of the bodies of interest and forces from asteroids and planetary satellites (see details in Ref. [21]). To determine the orbits of planets and the spacecraft, one must also describe the propagation of electromagnetic signals between any two points in space. The light-time equation corresponding to the metric tensor (3) and written 5 to the accuracy sufficient for GRAIL has the form (see also Ref. [22, 23]): X GMb rb + rb + rb |r2 − r1 | 2 12 + O(c−5 ), t2 − t1 = ln 1b + (1 + γ) b − rb c c3 r + r 1 2 12 b (7) where t1 refers to the signal transmission time and t2 refers to the reception time, while r1,2 are the barycentric b positions of the transmitter and receiver. Also, r1,2 are the distances of the transmitter and receiver from the body b b and r12 is their spatial separation [18, 20]. The logarithmic contribution in (7) is the Shapiro gravitational time delay that, in the case of GRAIL, is mostly due to the Moon, the Earth, and the Sun. (Note that the O(c−5 ) terms are beyond GRAIL’s sensitivity; see the analysis in Sec. III B.) The general relativistic equations of motion (6) and light-time equation (7) are used to produce numerical codes for the purposes of constructing Solar System ephemerides, spacecraft navigation [18, 21] and analysis of gravitational experiments in the Solar System [20, 24]. GRAIL also relies on these equations to compute its range and rangerate observables between the two spacecraft in lunar orbit. The numerical algorithm developed for this purpose [3] iteratively solves the light-time equation (7) in the SSB frame in terms of the instantaneous distance between the two spacecraft. Our objective is to develop an explicit analytical model for all the quantities involved in these highprecision computations. For this purpose, we need a clearly defined set of astronomical reference frames, which we discuss next. B. Relativistic coordinate transformations between various reference frames To describe the dynamics of an N -body system (such as the Solar System) in general relativity, one may choose to introduce N + 1 reference frames, each with its own coordinate chart. We need one global coordinate chart defined for the inertial reference frame that covers the entire system under consideration (e.g., BCRS). In the immediate vicinity of each of the N bodies in the system we can also introduce a set of local coordinates defined in the frame associated with this body (body-centric system). In the remainder of this paper, we use {xm } to represent the coordinates of the global inertial frame and {yam } to be the local coordinates of the accelerated proper reference frame of body a. In Ref. [15], we showed that the transformations between the harmonic coordinates of the BCRS {xm } and nonrotating body-centric reference systems {yam } (such as the GCRS or LCRS) may be written in the following form: Z ya0h n o i ′0 a 1 2 + O(c−4 ), (8) dy x0 = ya0 + c−2 c(va · ya ) + v + U a ext 2 a 0 ya 0 x = ya + za + c−2 n 1 2 va (va o a · ya ) − ya Uext + [ωa × ya ] + 21 aa ya2 − ya (ya · aa ) + O(c−4 ), (9) where za is the vector that connects the origin of the {xm } reference system with the origin of the {yam } reference system. Note that the accuracy of timing for GRAIL is limited by the performance of the on-board USO, which have an error of O(10−13 ) for 103 s of integration time [6]. Therefore, the c−4 terms in Eq. (8), which are at most of order ∼ v 4 /c4 ≃ 10−16 , are negligible for GRAIL, even in the absolute sense. The differential nature of the observables on GRAIL further reduces the sensitivity of the mission to such small terms in the transformations. For a complete post-Newtonian form of these transformations, including the terms c−4 and their explicit derivation, consult Ref. [15]. The inverses of the transformations (8)–(9) can be written as Z x0h i n o a ′0 1 2 ya0 = x0 − c−2 c (va · ra ) + v + U dx + O(c−4 ), (10) ext 2 a x00 ya = ra + c−2 n 1 2 va (va o a · ra ) + ra Uext + [ω a × ra ] + ra (ra · aa ) − 12 aa ra2 + O(c−4 ), (11) a where ra = x−za . The quantity Uext in Eqs. (8)–(9) and (10)–(11) is the Newtonian gravitational potential (including, if necessary, multipole corrections) due to all bodies in the Solar System other than body a, at the location of body a. Furthermore, aa is the Newtonian acceleration of body a due to the combined gravity of all other bodies. Later in a this section, we will present the expressions for Uext and aa for each of the chosen reference frames. µν Finally, ω a in Eqs. (9) and (11) is the vector associated with the relativistic precession given as ωaα = 21 ǫα µν ωa , 1 αβ with ǫα µν being the fully antisymmetric Levi-Civita symbol, normalized as ǫ23 = 1, and the matrix ωa having the form [15]: o X GMb n β β 3 α α −2 3 β nα ). (12) ω̇aαβ = − ba ( 2 va − 2vb ) − nba ( 2 va − 2vb ) + O(c 2 rba b6=a 6 The expression for the relativistic precession matrix is given here only for the sake of completeness. Because of their small magnitude (∼ 10−15 m), these terms will not affect the GRAIL measurements (see discussion in Sec. II F). In the rest of this section, we discuss four fundamental body-centric reference frames that are useful for collection and interpretation of GRAIL data. C. Coordinate systems used in the vicinity of the Earth In the vicinity of the Earth, two standard coordinate systems are utilized: the Geocentric Coordinate Reference System (GCRS), centered at the Earth’s center of mass is used to track orbits in the vicinity of the Earth. The positions of objects on the surface of the Earth, such as DSN ground stations, are usually given in the Topocentric Coordinate Reference System (TCRS). 1. Geocentric Coordinate Reference System (GCRS) When constructing a body-centric coordinate reference frame for a body a at the level of accuracy anticipated for a GRAIL, it is sufficient to consider only monopole contributions to the external potential Uext of all the bodies in the Solar System (for the Earth it is mostly the Sun and the Moon) excluding body a itself [15]. Thus, for the GCRS, a the Newtonian potential of the external bodies (i.e., excluding the Earth) Uext and the corresponding acceleration aa that are present in the coordinate transformations Eqs. (8)–(9), have the form [15]: E Uext = X GMb + O(c−2 ), rbE b6=E E aE = −∇Uext =− X b6=E GMb rbE −2 ), 3 + O(c rbE (13) where summation is performed over all the bodies excluding the Earth (symbolically, b 6= E), the vector that connects body b with the Earth’s center of mass is represented by rbE = xE − xb and the contributions of the higher multipole moments of mass distribution within the bodies are neglected due to their smallness. a The transformations given by Eqs. (8)–(9), together with the potential Uext and the acceleration aE given by E Eq. (13), determine the metric tensor gmn of the non-rotating GCRS [15]. We denote the coordinates of this reference m 0 E frame as {yE } ≡ (yE , yE ) and present the metric tensor gmn in the following form: E g00 = 1− 2 2 2 + O(c−6 ), w[E] + 4 w[E] 2 c c E g0α = −γαλ 4 λ w + O(c−5 ), c3 E E gαβ = γαβ + γαβ 2 wE + O(c−4 ), c2 (14) λ where wE and wE are the scalar and vector harmonic potentials that are given by wE = UE + utidal + O(c−4 ), E (15) wE = − (16) G −2 ), 3 [yE × SE ] + O(c 2yE where SE in Eq. (16) is the Earth’s angular momentum. The scalar potential wE is formed as a linear superposition of the gravitational potential UE of the isolated Earth and the tidal potential utidal produced by all the Solar System E bodies (excluding the Earth itself, b 6= E) evaluated at the origin of the GCRF. The Earth’s gravitational potential UE at a location defined by spherical coordinates (yE , φ, θ) is given by UE = +ℓ ∞ X o X GME n R0E ℓ E E 1+ Pℓk (cos θ)(Cℓk cos kφ + Sℓk sin kφ) + O(c−4 ), yE yE (17) ℓ=2 k=0 E E where R0E being the Earth’s radius, Pℓk are the Legendre polynomials, while Cℓk and Sℓk are spherical harmonic coefficients that characterize the Earth. At the level of sensitivity of GRAIL, only the lowest order spherical harmonic coefficients need to be accounted for, and time-dependent contributions due to the elasticity of the Earth can be ignored. Insofar as the tidal potential utidal is concerned, for GRAIL it is sufficient to keep only its Newtonian E contribution (primarily due to the Moon and the Sun) which can be given as usual: utidal = E X b6=E X GM b 3 −2 3(nbE · yE )2 − y2E + O(yE , c ), Ub (rbE + yE ) − Ub (rbE ) − yE · ∇Ub (rbE ) ≃ 3 2rbE b6=E (18) 7 where Ub is the Newtonian gravitational potential of body b, rbE is the vector connecting the center of mass of body b with that of the Earth, and ∇ denotes the divergence with respect to yE . Note that in Eq. (18) we omitted relativistic tidal contributions of O(c−2 ) that are produced by the external gravitational potentials. These are of the order of 10−16 compared to UE and, thus, completely negligible for GRAIL. In addition, we present only the largest term in 2 the tidal potential of the order of ∼ yE ; however, using the explicit form of the tidal potential Eq. (18), one can easily evaluate this expression to any order needed for a particular problem. The proper time at the origin of the GCRS is called the Geocentric Coordinate Time (TCG), denoted here as tTCG . It relates to the barycentric time TCB t as X GMb i 1 h v2 dtTCG = 1− 2 E + + O(c−4 ) ≈ 1 − 1.48 × 10−8 . dt c 2 rbE (19) b6=E The Earth’s barycentric velocity vE and position zE can be computed from Eq. (6). 2. Topocentric Coordinate Reference System (TCRS): proper and coordinate times To obtain the metric of the topocentric coordinate reference system, the TCRS, one can transform the metric E C gmn of the GCRS using coordinate transformations given by Eqs. (8)–(9), where the “external” potential Uext is the gravitational potential wE given by Eq. (15) and evaluated at the surface of the Earth: X C = wE (yC ) = UE (yC ) + Uext Ub (rbE + yC ) − Ub (rbE ) − yC · ∇Ub (rbE ) + O(c−2 ), (20) b6=E aC = C −∇Uext = −∇UE (yC ) − X b6=E ∇Ub (rbE + yC ) − ∇Ub (rbE ) + O(c−2 ), (21) where yC is the position vector of the DSN station in the GCRS. Note that UE (yC ) must be treated as the potential of an extended body and include a multipolar expansion with sufficient accuracy, taking into account time-dependent terms due to tidal effects on the elastic Earth. The proper time τC , kept by a clock located at the GCRS coordinate position yC (t), and moving with the coordinate velocity vC0 = dyC /dtTCG = [ΩE × yC ], where ΩE is the angular rotational velocity of the Earth at C, is determined by i X GMb dτC 1h 3 3(nbE · yC )2 − y2C + (aE · yC ) + O(yC , c−4 ), = 1 − 2 12 v2C0 + UE (yC ) + 3 dtTCG c 2rbE (22) b6=E where aE is the Earth’s acceleration in the BCRS, Eq. (13), and nbE is a unit spatial vector in the body-Earth direction, i.e., nbE = rbE /|rbE |, where rbE is the vector connecting body b with the Earth. The term within the square brackets in Eq. (22) is the sum of Newtonian tides due to the Sun, the Moon, and other bodies at the clock location yC . These terms are small for Earth stations (of order 2 × 10−17 ) and are negligible for GRAIL. The last term is due to non-inertiality of the GCRS and accounts for the Earth’s finite size. This term is evaluated to be of the order of 4.2 × 10−13 , which is about 10−3 smaller compared to the gravity potential on the surface of the Earth and, thus, it can be omitted. Therefore, at the accuracy required for GRAIL, it is sufficient to keep only the first two terms in Eq. (22) when defining the relationship between the proper time τC and the coordinate time tTCG : i 1h dτC = 1 − 2 21 v2C0 + UE (yC ) + O(c−4 ). dtTCG c (23) At the level of accuracy required for GRAIL, it is important to account in Eq. (23) for the oblateness (non-sphericity) of the Earth’s Newtonian potential, which is given in the form of Eq. (17). In fact, when we model the Earth’s gravity potential, we need to take into account quadrupole and higher moments, time-dependent terms due to tides as well as the tidal displacement of the DSN station. For example, for a clock situated on the surface of the Earth, the relativistic correction term appearing in Eq. (23) is given at the needed precision by v2C0 + UE (yC ) = W0 − 2 Z 0 hC gdh, (24) 8 where W0 = 6.2636856 × 107 m2 /s2 is the Earth’s potential at the reference geoid while g denotes the Earth’s acceleration (gravitational plus centrifugal), and where hC is the clock’s altitude above the reference geoid. Finally, we present the relation of the proper time read by the clock on the surface of the Earth at point C with respect to the TCB. Expressing dτC /dt = (dτC /dtTCG )(dtTCG /dt), with the help of Eq. (8) together with Eqs. (19) and (23), at the level of accuracy sufficient for GRAIL, we have: i X dτC 1h Ub (rbE + yC ) + O(10−17 ), = 1 − 2 12 (vE + [ΩE × yC ])2 + UE (yC ) + dt c (25) b6=E where the first term in the brackets, vE + [ΩE × yC ] ≡ vE + vC0 = vC , is the barycentric velocity of the DSN station. For details on the recommended relativistic formulation of GCRS consult Refs. [11, 13, 18]. Coordinate transformations (in particular, transformations involving topocentric coordinates) are discussed extensively in the IERS Conventions2 . D. Relativistic timekeeping in the Solar System Spacecraft radio science observations are clock and frequency measurements made at Earth stations [18]. For this purpose, the time coordinate called Terrestrial Time (TT) is defined. TT is related to TCG linearly by definition: dtTT = 1 − LG , dtTCG (26) where LG = 6.969290134 × 10−10 by definition. This definition accounts for the secular term due to the Earth’s potential when converting between TCG and the time measured by an idealized clock on the Earth geoid [11–13, 18]. Using Eq. (23), we also have i dτC dtTCG 1h dτC = = 1 + LG − 2 21 v2C0 + UE (yC ) + O(c−4 ). dtTT dtTCG dtTT c (27) On the other hand, equations of motion in the Solar System are often evaluated using another defined time scale, TDB. TDB time (tTDB ) is also related to TCB time t linearly: dtTDB = 1 − LB , dt (28) where LB = 1.550519768 × 10−8 by definition, accounting for all secular terms due to the solar gravitational field, the Earth’s orbital velocity, and the Earth potential on the geoid. The relationship between TCB and TCG is nonlinear; these are the coordinate times of two coordinate systems related to one another by the space-time transformations (8)–(9) and their inverses (10)–(11). The relationship between TT and TDB, therefore, is also nonlinear. The difference is dominated by an annual periodic term with an amplitude of ∼ 1.6 × 10−3 s. The definition of TT and TDB ensures the absence of a significant linear term. For accurate computations in the SSB reference frame, observed times of transmission and reception need to be converted from TT to TDB. E. Coordinate reference frames in the vicinity of the Moon In the vicinity of the Moon, once again we consider two coordinate systems. The lunicentric LCRS is a coordinate system used, for instance, to represent lunar orbits. To describe experiments carried out on board the GRAIL spacecraft, we use the SCRS. 2 All software, technical specification and other relevant materials associated with the IERS Conventions (2010) can be found at http://www.iers.org/ 9 1. Lunar Coordinate Reference System (LCRS) In complete analogy to the formulation of the GCRS (discussed in Sec. II C 1) and similarly to the approach advocated in Ref. [9], the metric tensor of the LCRS may be obtained by transforming the metric (3) and the potentials (4)–(5) of the BCRS using the coordinate transformations given by Eqs. (8)–(9) where it is sufficient to a consider only the monopole contribution to Uext from all the bodies of the Solar System excluding the Moon: M Uext = X GMb + O(c−2 ), rbM M aM = −∇Uext =− b6=M X b6=a GMb rbM + O(c−2 ), 3 rbM (29) where summation is performed over all the bodies excluding the Moon (b 6= M) and the contributions due to the higher multipole moments of the mass distributions within the bodies are neglected. Applying the coordinate transformations given by Eqs. (8)–(9) together with the external potential and acceleration M (29), one can derive the metric tensor gmn of the non-rotating lunar coordinate reference system (LCRS). Denoting m 0 the coordinates of the LCRS as {yM } ≡ (yM , yM ), this tensor may be presented in the following form (which is identical to the expressions in Sec. II C 1 after making the substitution E → M): M g00 = 1− 2 2 2 wM + 4 wM + O(c−6 ), c2 c M g0α = −γαλ 4 λ w + O(c−5 ), c3 M M gαβ = γαβ + γαβ 2 w + O(c−4 ), c2 M (30) λ where the scalar and vector potentials wM and wM are given as: wM = UM + utidal + O(c−4 ), M (31) wM = − (32) G −2 ), 3 [yM × SM ] + O(c 2yM where SM in Eq. (32) is the Moon’s angular momentum. Similarly to the GCRS, the scalar potential wM (31) is a linear superposition of the proper gravitational potential of the Moon (with R0M being the Moon’s radius, MM its M M mass, while Cℓk and Sℓk are the Moon’s spherical harmonic coefficients): UM = ∞ X +ℓ o X R0M ℓ GMM n M M 1+ Pℓk (cos θ)(Cℓk cos kφ + Sℓk sin kφ) + O(c−4 ) yM yM (33) ℓ=2 k=0 plus tidal contributions utidal produced by all the Solar System bodies (excluding the Moon itself, b 6= M) evaluated M at the origin of the LCRF. For the GRAIL’s accuracy, it is sufficient to keep only the Newtonian contribution to the tidal potential produced by the external bodies which can be presented as X GM X b 2 2 3 utidal = Ub (rbM + yM ) − Ub (rbM ) − yM · ∇Ub (rbM ) ≃ 3(n · y ) − y + O(yM , c−2 ), (34) M bM M M 3 2rbM b6=M b6=M where nbM is a unit spatial vector in the body-Moon direction, i.e., nbM = rbM /|rbM |, where rbM is the vector connecting body b with the Moon. Note that the relativistic tidal contributions of 1/c2 order that are due to external potentials have a magnitude of 10−16 when compared to UM and, thus, they were omitted in Eq. (34). In addition, 2 we present only the largest term in the tidal potential of the order of ∼ yM ; however, using the explicit form of the tidal potential Eq. (34), one easily evaluate this expression to any order needed for a particular problem. M At the same time, we must account for deformations of the elastic Moon, expressed in the form of corrections ∆Cℓk M and ∆Sℓk to the lunar spherical coefficients, due to the tidal potential of body b, located at lunicentric spherical coordinates (rbM , φbM , θbM ) [19, 25]: ℓ+1 s M R (ℓ + 2)[(ℓ − k)!]3 ∆Cℓk cos kφbM b 0M M = 4kℓ Pℓk (cos θbM ) . (35) ∆Sℓk sin kφbM MM rbM [(ℓ + k)!]3 The lunar Love number k2M ≃ 0.025 [26] thus introduces a significant time-dependent contribution to the spherical M M harmonic coefficients Cℓk and Sℓk , which must be written as the sums M M0 M Cℓk = Cℓk + ∆Cℓk and M M0 M Sℓk = Sℓk + ∆Sℓk , M0 M0 where we used Cℓk and Sℓk to denote the constant part of the lunar spherical harmonic coefficients. (36) 10 2. Lunar Coordinate Time (TCL) There are several different time coordinates to be considered for GRAIL. In addition to the terrestrial time scales defined in Sec. II D, GRAIL also relies on the timing events reported at proper times measured by clocks on board the lunar orbiters. Thus, one would need to introduce a realization of lunar coordinate time (TCL) and a spacecraft proper time (ST). The lunicentric orbits of the GRAIL spacecraft are coupled to the orbit of the Moon mostly through the difference between the acceleration of the probe and that of the Moon due to the gravitational pull of the Earth and the Sun (the Earth’s and the Sun’s tidal terms). This coupling is weak because the Earth and Sun tides are, respectively, just 2.5 × 10−7 and 4.7 × 10−8 times the monopole acceleration due to the Moon. Relativistic perturbations containing the mass of the Moon are small (∼ 4.6 × 10−11 m/s2 ) to the point that they are not measurable, being easily absorbed into the much larger non-gravitational perturbations (for instance, solar radiation pressure). Should we conclude that general relativity does not matter in the computation of the the lunicentric orbit of the spacecraft? The answer is negative, but the main relativistic effect does not appear in the equation of motion. According to Eq. (8), the differential equation that gives the local proper time tTCL at the origin of the LCRS as it relates to the barycentric time TCB t is i X dtTCL 1 h v2 = 1− 2 M + Ub (rbM ) + O(c−4 ), dt c 2 (37) b6=M where the Moon’s barycentric velocity vM and position zM can be computed from the EIH equations (6). Equation (37) establishes the relationship between the TCL (tTCL ) and TCB (t) time scales. Truncated to the first post-Newtonian (1PN) order (we put the clock at the origin of its proper reference system, yM = 0 and drop on the right hand side O(c−4 ) terms that are in principle known, but certainly not needed for our purposes), it is given by a differential equation X GMb i dtTCL 1 h v2 + O(c−4 ) ≈ 1 − 1.48 × 10−8 , =1− 2 M + dt c 2 rbM (38) b6=M which can be solved by a quadrature formula once the orbits of the Moon, the Sun and the other planets are known. 3. Satellite Coordinate Reference System (SCRS) To determine the metric tensor for the satellite coordinate reference system, the SCRS, we perform the coordinate transformation given by Eqs. (8)–(9), where the “external” potential and acceleration determined by the potential wM given by Eq. (31), taken at the lunicentric position yA of the spacecraft GRAIL-A are (the equations are identical for GRAIL-B, except for the substitution A → B): X A Uext = wM (yA ) = UM (yA ) + Ub (rbM + yA ) − Ub (rbM ) − yA · ∇Ub (rbM ) + O(c−2 ), (39) b6=M where yA is the solution of the equations of motion of the GRAIL-A spacecraft in the lunicentric frame. This equation can be obtained from equations of geodesics and the metric tensor of the LCRS (30) with relativistic gravitational potentials given by (31)–(32), (33), and (39). Including all the terms of the order of ∼ 10−12 m/s2 and larger, the equation of spacecraft motion of the GRAIL spacecraft in the LCRS takes the form: X aA0 = −∇UM (yA ) − ∇Ub (rbM + yA ) − ∇Ub (rbM ) + b6=M + GMM 4GMM 2 2 + aNG + O(10−13 m/s ), n − v n + 4(n · v )v A A A A0 A0 A0 2 c2 y A yA (40) where aA0 ≡ d2 yA /dt2TCL and vA0 = dyA /dtTCL are the lunicentric acceleration and velocity of the spacecraft in a non-rotating LCRS, also nA = yA /yA is the unit vector in the direction of the spacecraft, yA = |yA | and aNG is the contribution of non-gravitational forces affecting the motion of a spacecraft (e.g., solar radiation pressure, thermal imbalance, outgassing, etc). The first term in Eq. (40) is the contribution of the lunar gravity potential UM (yA ), given by Eq. (33). Note that the potential UM (yA ) must be treated as the potential of an extended body and include a multipolar expansion 11 with sufficient accuracy. The second term in Eq. (40) is due to the tidal potential at the location of the spacecraft produced by external bodies utidal (mostly the Earth and the Sun) which is given by Eq. (34). To reach GRAIL’s M accuracy requirement of ∼ 10−12 m/s2 , one would have to account for several terms in the expansion beyond the 5 2 second order one ∼ yM given in Eq. (34). In fact, terms up to ∼ yM in the tidal potential utidal are needed. The M group of terms on the second line of Eq. (40) is the relativistic Schwarzschild perturbation due to the spherically symmetrical component of the Moon’s gravitational field. The first two terms in this group are of the order of ∼ 1.86 × 10−10 m/s2 and ∼ 4.63 × 10−11 m/s2 , respectively. These are large enough to be in the equations of motion. Given the nearly-circular orbit of the GRAIL spacecraft, the magnitude of the last term in this group is reduced by the orbital eccentricity, which is eA ∼ 0.018. This fact reduces the contribution of this term by nearly two orders of magnitude when compared to the first two terms, making it barely observable with GRAIL. Note that the lunar angular momentum SM present in the relativistic vector gravity potential Eq. (32) produces contribution to Eq. (40) of the order of ∼ 10−14 m/s2 , which makes it negligible for GRAIL. A m 0 A a result, the metric tensor gmn representing space-time in the coordinates {ŷA } ≡ (ŷA , ŷA ) of the proper nonrotating spacecraft coordinate reference system (SCRS) may be given in the following form: A g00 = 1− 2 wA + O(c−6 ), c2 A g0α = O(c−5 ), A gαβ = γαβ + γαβ 2 w + O(c−4 ), c2 A where wA is the tidal contribution produced by the Moon on the world-line of the spacecraft: GMM 2 3 −4 2 ). 3(n · ŷ ) − ŷ wA = A A + O(ŷA , c A 3 2ŷA (41) (42) We can also determine the differential equation that relates the rate of the spacecraft proper τA time, as measured by an on-board clock in lunar orbit, to the time in LCRS, tTCL : i X dτA 1 h v2 = 1 − 2 A0 + UM (yA ) + Ub (rbM + yA ) − Ub (rbM ) − yA · ∇Ub (rbM ) + (aM · yA ) + O(c−4 ), (43) dtTCL c 2 b6=M where aM is the barycentric acceleration of the Moon, Eq. (29), see Ref. [15]. As a result, we can establish the rate of the spacecraft proper time with respect to the time of the BCRS, t = tTDB : i X 1h dτA = 1 − 2 12 (vM + vA0 )2 + UM (yA ) + Ub (rbM + yA ) + O(c−4 ). (44) dt c b6=M This result summarizes the relationship of the proper time of an on-board clock in lunar orbit τA and the TDB. 4. Transformation of gravitational potentials To complete the description of the LCRS, we present the transformation rules for the relativistic gravitational potentials. In Ref. [15], we obtained the structure of the metric tensors corresponding to the local space-times in the reference frames relevant for GRAIL, expressed in terms of harmonic gauge potentials w and w. We also derived the rules for transforming relativistic gravitational potentials, the coordinate transformations between the frames and resulting relativistic equations of motion. Applying these results to GRAIL, we see that the scalar and vector gravitational potentials of the Moon wM (yM ) and wM (yM ) (as measured at the LCRS) relate to those measured in the coordinates of the BCRS wM (rBCRS ) and wM (rBCRS ) as wM (yM ) = 2v 2 4 1 + 2M wM (rBCRS ) + 2 vM · wM (rBCRS ) + O(c−4 ), c c wM (yM ) = wM (rBCRS ) − vM wM (rBCRS ) + O(c−2 ). (45) (46) We estimate the magnitude of wM given by Eq. (32) as ∼ 2.4 × 106 m3 /s3 . This results in a value of ∼ 3.2 × 10−6 m2 /s2 for the third term in (45), which is four orders of magnitude too small compared to the second term in 2 that expression, the scalar potential of the Moon multiplied by 2vM /c2 that was evaluated to be ∼ 5.5 × 10−2 m2 /s2 . Thus, to determine the relationship between the LCRS-defined mass of the Moon and its barycentrically defined mass, 2 we must multiply the latter by the factor (1 + 2vM /c2 ) ≈ (1 + 2 × 10−8 ). Such a transformation results in a small, but observable effect. As far as the GRAIL’s accuracy in concerned, contributions to other multipoles of the lunar gravity field are not sensitive to such a small correction factor. 12 As we know, GRAIL determines the lunar gravity field relying on the EIH equations of motion (6) for the bodies of the Solar System, including the Moon, the Earth, and the GRAIL spacecraft. For this, Eq. (6) must also includes the gravitational potential of the extended Moon (33) and tidal potentials due to the Earth and the Sun (34). Therefore, one would have to transform the resulting barycentrically defined gravitational potential of the Moon from coordinates of the BCRS to those of the proper lunicentric frame LCRS. A concern was that such a procedure may lead to some unwanted biases in the determination of the lunar gravity field. Our approach allows one to evaluate the general relativistic effects on the largest coefficients to the lunar gravity potential corresponding to this transformation. Substituting Eq. (45) into Eq. (40), we essentially modify the equation 2 of motion of the GRAIL spacecraft by accounting for the (1 + 2vM /c2 ) factor. Now, we can represent the barycentric ′ velocity of the Moon as vM = vEM + vM , where vEM is the barycentric velocity of the Earth-Moon barycenter and 2 2 ′2 ′ v′M is the velocity of the Moon in the EMB frame. Therefore, vM = vEM + vM + 2vEM vM cos D, where D = θ − θ′ is the the difference between the longitudes of the mean Moon and the mean Sun with a period of 29.531 days. The constant part in the barycentric velocity of the Moon vM may be easily absorbed in the determination of the lunar mass as a bias with magnitude of 2 × 10−8 MM . The variability in vM , if not properly removed in accordance with Eq. (45), may introduce an additional time-dependent bias with a magnitude of 2.2 × 10−9 MM cos D, which can be removed in the data analysis. Finally, given the nearly circular obit of the GRAIL spacecraft around the Moon, the 1/c2 terms in the second line of Eq. (40) can be seen and a modification of the Newtonian point-mass acceleration 2 −10 N of the Moon aN aA and A = GMM nA /yA . These terms are nearly constant, have combined magnitude of ∼ 1 × 10 −10 would be easily absorbed in the determination of the lunar mass as a small bias of 1 × 10 MM . For spacecraft with lesser sensitivity, these corrections are irrelevant. However, at the micron-level sensitivity of the GRAIL mission, they become noticeable. It is, of course, possible to absorb small constant or periodic terms into constants such as MM during data analysis, with no impact on mission objectives or the quality of the mission’s results. Nonetheless, pursuing these small corrections is worthwhile, demonstrating that a spacecraft with GRAIL’s sensitivity is already a practical instrument for relativistic geodesy in the lunar environment, and paving the way for future missions that will operate at even greater accuracy. F. Transformations of position vectors Equation (9) establishes the relationship between the coordinates of the local body-centric coordinate reference frame and the coordinates of the global BCRS [15]: o n a + [ω a × ya ] + 12 aa ya2 − ya (ya · aa ) + O(c−4 ). (47) ra = ya + c−2 − 12 va (va · ya ) − ya Ūext Considering the anticipated accuracy of the GRAIL experiment, we can simplify Eq. (47) by noting that the last three terms in this expression are much smaller than needed for GRAIL. Indeed, we can evaluate the magnitude of the third term [15] as [ω a × ya ] ≃ ya GM⊙ vE ∆t/AU2 , where ∆t is the signal propagation time, M⊙ is the mass of the Sun and AU ≃ 1.5 × 1011 m is the astronomical unit. The Moon-Earth radio-signal propagation time is ∆t ≃ 1.3 s. a ; therefore, the third term within the curly braces However, even for ∆t = 103 s, we have [ω a × ya ] . 2 × 10−4 ya Uext in Eq. (47) is negligible for GRAIL. The last two terms within the curly braces in Eq. (47) are dependent on the acceleration of the planet center, and are ignored for the reasons that we discuss at the end of this subsection. The required space-time transformations relate the position of the ground-based antenna and that of the spacecraft. Using superscript indices to indicate explicitly the dependence on the various time scales TT, TDB, etc., the terrestrial (geocentric) coordinates yTT of the antenna must be transformed into TDB-compatible (barycentric) coordinates C rTDB = zC − zE + O(c−4 ), where zC being the barycentric position of the DSN station zC = zE + yC . Similarly to C the approach developed in Ref. [27] for the BepiColombo mission, this transformation is expressed by TDB 1 X 1 TT rTDB = y 1 − vE , · yTT U (r ) − 2 vTDB b bE C C E C 2 c 2c (48) b6=E P where b6=E Ub (xE ) is the Newtonian gravitational potential due to bodies other than the Earth at the geocenter, xE is the SSB position and vTDB = dzE /dt is the SSB velocity of the Earth (as defined in the paragraph before Eq. (6)). E The time coordinate must also be changed consistently together with the spatial coordinates. The effect of this change on velocities is given by: h 1 X 1 TDB TT TDB i dτC TDB TT 1 − v · vC0 vE , U (r ) − vTDB − v = v b bE C E C0 c2 2c2 E dt b6=E (49) 13 TT where vTT C0 = dyC /dτC . Note that Eq. (49) contains the factor dτC /dt, identical to Eq. (19), that deals with time transformation: τC is the local time for a ground-based antenna, that is, TT, and t is the corresponding TDB time. Similar to Eq. (48), the lunicentric coordinates of the orbiter yTCL are transformed into BCRS coordinates rTDB = A A −4 zA − zM + O(c ), where zA being the barycentric position of the spacecraft zA = zM + yA : 1 TDB TCL TDB 1 X TCL vM , v · yA U (r ) − 1 − rTDB = y b bM A A c2 2c2 M (50) b6=M P where b6=M Ub (xM ) is the gravitational potential due to bodies other than the Moon at the Moon’s barycenter, xM is the Moon’s position in SSB coordinates and vTDB = dzM /dt is the Moon’s SSB velocity. M The corresponding velocity transformation is given by i dt h 1 X 1 TCL vTDB 1− 2 · vTCL vTDB − vTDB = vTCL , Ub (rbM ) − 2 vTDB M M A0 A M A0 c 2c dt (51) b6=M TCL with vTCL /dtTCL and dtTCL /dtTDB given by Eq. (38). The relations for the other spacecraft is obtained by A0 = dyA replacing A → B. Note that in all the coordinate transformations presented in this section, we neglected terms that contain the SSB acceleration of the planet center, as these have an additional small parameter (Rb /zb ), where Rb is distance from the planetary barycenter and zb is the distance of the planet from the SSB. Even for the Earth-Moon distance, these acceleration-dependent terms are at most of the order of 10−3 compared to the other 1/c2 terms, both velocity- and Newtonian potential-dependent ones, making the acceleration-dependent terms negligible for the results above. III. FORMING Ka-BAND RANGE (KBR) OBSERVABLES FOR GRAIL One can demonstrate that in the post-Minkowskian approximation, appropriate for most Solar System experiments including GRAIL, as seen from BCRS, the phase of an electromagnetic wave that is passing by a gravitating body with the mass M can be presented as (see a detailed derivation in Appendix A 2, with the result given by Eq. (A22)): 2GM h rA + r + RA i + O(G2 , c−3 ), ln ϕ(t, x) = k0 ct − RA − c2 rA + r − RA (52) where k m = k 0 (1, k) is a constant null vector directed along the trajectory of propagation of the unperturbed electromagnetic wave such that γmn k m k n = 0, also k 0 = ω/c where ω is the constant frequency of the unperturbed wave. We also use the following notations: k= RA , RA RA = x − xA , RA = |RA |, also r = x − z, r = |r| rA = xA − z, rA = |rA |, (53) where z is the time-dependent spatial coordinate of the massive body. The phase ϕ of an electromagnetic wave that was emitted at the point xm A0 = (ctA0 , xA0 ) and received at the point m xB = (ctB , xB ) remains constant along the path of this wave [28, 29]. In particular, if λ is an affine parameter along the wave’s path, the derivative of the phase satisfies the equation dϕ ∂ϕ dxm = = Km K m = 0, dλ ∂xm dλ (54) which suggests that ϕ[xm (λ)] = const. In other words, along the signal’s world-line the phase stays constant and equal to its initial value ϕ(t, x) = k0 ctA0 = ωA0 tA0 . Equating the values of the phase given by Eq. (52) at two points A0 and B as ϕ(tA0 , xA0 ) = ϕ(tB , xB ), we can determine the gravitational delay of the signal moving through a particular space-time. Indeed, up to O(c−3 ) the coordinate time transfer, which is defined as tB − tA0 = TAB , is given by [15]: TAB = tB − tA0 = 2GM h rA0 + rB + RAB i RAB + O(c−4 ), ln + c c3 rA0 + rB − RAB (55) where the logarithmic term represents the Shapiro time delay. Also, rB = xB − z and RAB = xB − xA0 . With these results, we can now formulate a relativistic model for the fundamental timing observables on GRAIL. 14 FIG. 2: Timing events on GRAIL: Depicted (not to scale) are the world-lines of the GRAIL-A and GRAIL-B spacecraft with corresponding proper times τA and τB . Ka-band signals that are emitted from points A0 and B0 are received at B and A, respectively. The times tA0 , tB0 , tA and tB are the coordinate times at these points, as measured in the BCRS. A. The inter-spacecraft ranging observables Consider a clock with proper frequency fA0 , located at moving point A0 , that emits a signal with frequency fA0 at an instant of proper time τA0 measured on the world-line of the clock. This signal is received by the moving point B at an instant of the proper time τB taken at the world-line body B and the instantaneous phase of this signal is compared with the phase of the local oscillator with proper frequency fB0 of the clock located at point B. The measurable quantity is the difference between the instantaneous phases of the two signals compared at point B. Instrumentally, at point B one measures the fractional difference dnAB in the number of cycles dnB A0 received from the clock at point A0 and the number of the locally generated cycles dnB . Mathematically, this quantity may be expressed at the point B at an instance of the proper time dτB as: B dnAB = dnB − dnB A0 = fB0 dτB − fA0 dτB , (56) B where fA0 is the frequency of the oscillator A as detected at B. Assuming that the number of pulses sent from spacecraft A, dnA0 , and received on spacecraft B, dnB A0 , are the B = dn , we can express f via its value at the proper time of emission on spacecraft A (see also same, or dnB A0 A0 A0 Eq. (C3) below): B fA0 dnB dτA0 A0 dτA0 = = . fA0 dτB dnA0 dτB (57) Furthermore, the instantaneous difference of the number of cycles measured on spacecraft B, as given by Eq. (56), takes the form: dnAB = fB0 dτB − fA0 dτA0 . (58) Eq. (58) is the difference in the number of cycles generated by the two oscillators during the given proper time intervals along the world-lines of the two clocks. We can now express Eq. (58) in terms of the coordinate time: dnAB = fB0 dτ B dtB dtB − fA0 dτ A0 dtA0 dtA0 . (59) Note that Eq. (59) cannot be integrated in general case if the two time variables tA0 and tB are treated as independent. However, in our case the points A0 and B are connected by a time-like geodesic, and therefore, the coordinate 15 times tA0 and tB are connected by the light-time equation (55) that reads: tB − tA0 = TAB (tA0 , tB ) = 1 2GM h rA + rB + RAB i . ln |rB (tB ) − rA (tA0 )| + c c3 rA + rB − RAB (60) We note that Eq. (60) can be used to express either tA0 as a function of tB or vice versa. Observables on GRAIL are time-stamped using the time of reception (that is, tB ). We therefore have more direct access to xA (tB ) rather than xA (tA0 ), and the first term on the right hand side of Eq. (60) gets modified by Sagnac correction terms (as observed in Ref. [30]) consistently to the order 1/c3 : RAB = dAB + (dAB · vA ) dAB 2 + 2 vA + (nAB · vA )2 − (dAB · aA ) + O(c−3 ), c 2c (61) where dAB = xB (tB ) − xA (tB ) is the coordinate distance between A and B at the moment of reception at B (we have dAB = |dAB | and nAB = dAB /dAB ), where vA = vA (tB ) denotes the coordinate velocity of spacecraft A at that instant, and where aB is the acceleration of A (in all the order 1/c3 terms we can use quantities at tA0 or tB ). In this case, Eq. (60) becomes TAB (tB ) = 2GM h rA + rB + dAB i dAB (dAB · vA ) dAB 2 2 + O(c−4 ), (62) + ln + v + (n · v ) − (d · a ) + AB A AB A A c c2 2c3 c3 rA + rB − dAB where all quantities here are taken at the instant of reception tB . In the case of the GRAIL mission, when signal transmission between the two spacecraft is concerned the first term in Eq. (61), which is of order 1/c, is ∼6.67 µs. The second term in Eq. (61) represents the Sagnac term of order 1/c2 and can amount to ∼3.67 ns at GRAIL’s orbit around the Moon; the third Sagnac term, of order 1/c3 , is ∼0.02 ps (comparable to the lunar Shapiro term, which is ∼0.04 ps). Note that expression for TBA (tA ) may be obtained from Eq. (62) by interchanging A ↔ B. B. Dual One-Way Range (DOWR) observables on GRAIL To develop an analytical form for the DOWR observable, we note that Eq. (60) could be used to express either tA0 as a function of tB or vice versa. As observables on GRAIL are time-stamped using the time of reception (that is, tB ), in the following we treat tA0 as a function of tB , i.e., tA0 = tA0 (tB ). Furthermore, we can write TAB (tA0 , tB ) = TAB (tB ). This allow us to present Eq. (59) as dτ dτ dnAB = fB0 dtB − fA0 d tB − TAB (tB ) . (63) dt B dt A To integrate Eq. (63), we rely on Eq. (44) and introduce a function uB (tB ) that allows us to write dτB /dt as: 1 dτB = 1 − 2 uB (tB ) + O(c−4 ), dt c where uB (tB ) = v2A X GMb + O(c−2 ). + 2 rbA (64) b Using this definition of uB (tB ) given in Eq. (64) allows us to integrate the first term in Eq. (63) as Z tB Z tB dτ dτ 1 uB (t′B )dt′B + O(c−4 ), dtB = fB0 tB − t0B + 2 fB0 uB (tB ) tB − t0B − fB0 dt B dt B c t0B t0B where t0B is used to denote the (for now, arbitrary) start of the integration interval. Similarly, we have the following expression for the second term of Eq. (63): Z tB dτ dτ fA0 tB − TAB (tB ) − t0B − TAB (t0B ) + d tB − TAB (tB ) = fA0 dt A dt A t0B Z tB −TAB (tB ) 1 + 2 fA0 uA (tB ) tB − TAB (tB ) − t0B − TAB (t0B ) − uA (t′B )dt′B + O(c−4 ). c t0B −TAB (t0B ) (65) (66) Transforming from proper to coordinate frequencies, as fB = fB0 (dτ /dt)B and fA = fA0 (dτ /dt)A , and using all the results developed in this section, we can integrate Eq. (63) and present the result in terms of the phase difference as (67) ∆nAB (tB ) = fB tB − fA tB − TAB (tB ) + ǫAB + δnAB + O(c−4 ), 16 where ǫAB ≡ ǫAB (t0B , tB ) is given by ǫAB (t0B , tB ) = Z tB 1n ′ ′ 0 u (t )dt u (t )(t − t ) − f B B B B B B0 B − B c2 t0B Z − fA0 uA (tB ) tB − TAB (tB ) − t0B − TAB (t0B ) − tB −TAB (tB ) t0B −TAB (t0B ) uA (t′B )dt′B o + O(c−4 ), (68) and δnAB ≡ δnAB (t0B ) is an integration constant determined by the initial conditions: δnAB (t0B ) = −fB t0B + fA t0B − TAB (t0B ) + O(c−4 ), (69) The second observable ∆nBA (tA ) that deals with the signal propagation from B0 to A is derived in an analogous way. To formulate the relativistic model for the dual one-way range (DOWR) observables on GRAIL, we need the expressions derived above for the instantaneous phase differences measured at both spacecraft, nAB (tB ) and nBA (tA ), which are given by Eqs. (67), together with the instantaneous delays measured at the points of signal reception at both spacecraft, TAB (tB ) and TBA (tA ), as given by Eqs. (61) and (62). As a result, Eq. (67) becomes: ∆nAB (tB ) = (fB − fA ) tB + fA TAB (tB ) + ǫAB + δnAB + O(c−4 ) d 2GMM h rA + rB + dAB i AB = (fB − fA ) tB + fA + ln + c c3 rA + rB − dAB (d · v ) d AB A AB 2 2 + fA + ǫAB + δnAB + O(c−4 ), + v + (n · v ) − (d · a ) AB A AB A A c2 2c3 (70) where MM is the mass of the Moon. An expression for ∆nBA (tA ) may be obtained from (70) by interchanging A ↔ B. The authors of Ref. [6] discuss the interpolation algorithm realized on GRAIL to synchronize the LGRS clocks on both spacecraft in coordinate time. Here we just assume that synchronization is achieved, so that tA = tB = t and t0A = t0B = t0 . We now can form a quantity, that is called dual one-way range (DOWR): Rdowr (t) = c ∆nAB (t) + ∆nBA (t) . fA + fB Substituting Eq. (70), we have the following result for Rdowr : (fA vA − fB vB ) · dAB dAB (fB aB − fA aA ) · dAB + + 2 Rdowr (t) = dAB + c(fA + fB ) 2c fA + fB dAB fA v2A + (nAB · vA )2 + fB v2B + (nAB · vB )2 2GMM h rA + rB + dAB i + 2 + + ln 2c fA + fB c2 rA + rB − dAB c δnAB + δnBA c ǫAB + ǫBA + + O(c−3 ), + fA + fB fA + fB (71) (72) where ǫAB and δnAB depend on the choice of the start of the integration intervals, i.e., t0A and t0B . We now discuss each of the seven terms present in Eq. (72) and evaluate their magnitudes and relevance for GRAIL. To develop numerical estimates for the magnitude of the various terms that we consider, we use mission parameters that are provided in Table I. The first term in Eq. (72) is the instantaneous Euclidean distance dAB ≃ 200 km (see Table I) between the two lunar orbiters. The next three terms are the first (∼ 1/c) and the second (1/c2 ) order Sagnac effects. These terms are due to the fact that representing the observables only in terms of the received times tA and tB on the two spacecraft is equivalent to a rotation of the reference system. To evaluate the second term in Eq. (72), we use the identity (fA vA − fB vB ) = − 21 (fA + fB )(vB − vA ) − 21 (fB − fA )(vA + vB ). Given ∆fAB = |fB − fA | ∼ 103 Hz, we get: (vA + vB ) · dAB (fA vA − fB vB ) · dAB (vAB · dAB ) fB − fA = − − = c(fA + fB ) 2c fA + fB 2c (73) 17 TABLE I: Select parameters of the GRAIL mission (some taken from Ref. [3]) and the Earth-Moon system, along with corresponding symbols and approximate formulae used in the text. Parameter Symbol(s) Values used GRAIL Mission Inter-spacecraft range dAB 200 km Inter-spacecraft range-rate d˙AB = (nAB · vAB ) 2 m/s Lunar altitude hG 55 km Lunicentric velocity vA0 = |vA0 | ≃ |vB0 | 1.65 km/s Relative spacecraft velocity vAB ≃ vA dAB /(RM + hG ) 185 m/s Lunicentric acceleration aA0 = |aA0 | ≃ |aB0 | 1.53 m/s2 Relative spacecraft acceleration aAB ≃ aA dAB /(RM + hG ) 0.17 m/s2 Ka-band frequency fA ≃ fB 32 GHz Frequency difference ∆fAB = fB − fA ∼ 103 Hz Earth-Moon system Moon’s geocentric velocity — 1 km/s EMB orbital velocity — 30 km/s DSN geocentric velocity — 465 m/s Earth mass parameter GME 3.98 × 1014 m3 /s2 Moon mass parameter GMM 4.90 × 1012 m3 /s2 Earth radius RE 6.371 × 106 m Moon radius RM 1.737 × 106 m = (−0.061423 + 3 × 10−7 ) m, (74) where vAB = vB − vA . Therefore, the second term in Eq. (74) is less than 1 µm and it can be omitted. The third term in Eq. (72) is the second order (∼ 1/c2 ) acceleration-dependent Sagnac effect. We evaluate this term in a manner similar to Eq. (74) and obtain the magnitude: (aA + aB ) · dAB fB − fA (aAB · dAB ) dAB (fB aB − fA aA ) · dAB = dAB − dAB 2c2 fA + fB 4c2 fA + fB 4c2 = (2 × 10−8 + 1.2 × 10−14 ) m. (75) Thus, the entire third term in Eq. (72) may be safely omitted. The fourth term on the right-hand side of Eq. (72) is the second order (∼ 1/c2 ) Sagnac effect. As a result this term may be evaluated as dAB fA v2A + (nAB · vA )2 + fB v2B + (nAB · vB )2 dAB 2 2 2 2 = v + (n · v ) + v + (n · v ) + AB A AB B A B 2c2 f +f 4c2 A B dAB fB − fA vAB · (vB + vA ) + (nAB · vAB ) nAB · (vB + vA ) = (0.002 + 2 × 10−13 ) m. (76) + 2 4c fA + fB Thus, the first term in (76) must be kept in the model. One can further evaluate this term by representing the barycentirc velocities of the GRAIL twins as vA = vM + vA0 and vB = vM + vB0 , where vM is the barycentric velocity of the Moon and vA0 and vB0 are the lunicentric velocities of the two orbiters. By doing this, one can see that there will be three terms, each of which is important for the GRAIL model. The term ∼ dAB (vM /c)2 contributes up to 2 mm to the DOWR. The term ∼ dAB (vM vA0 /c2 ) contributes up to 110 µm to the DOWR, and the last term ∼ dAB (vA0 /c)2 , also frequency-dependent, contributes up to 6 µm to this observable. Thus, each of these terms must be accounted for in the relativistic model of GRAIL observables. The fifth term in Eq. (72) is the Shapiro gravitational time delay. Assuming a spacecraft altitude hG = 55 km, this term contributes (4GMM /c2 )(dAB /(rA + rB )) = 12 µm to the DOWR and, thus, it may be accounted for in the range model in the following approximated form, keeping just the largest (12 µm) term: 2GMM h rA + rB + dAB i 4GMM dAB d3AB 4GMM ≈ ln + = 12 µm + 1.3 × 10−8 m. 2 2 2 c rA + rB − dAB c rA + rB 3c (rA + rB )3 (77) Concerning the sixth term in Eq. (72), in Appendix B we show that, for the times-scales of signal propagation 2 realized on GRAIL (dAB /c ≃ 1 ms), this term is of the order of 1/c4 and contributes less than 1 × 10−15(t − t0 )2 m/s to the DOWR. An acceleration error of this magnitude yields a range error of less than 1 µm over the course of 6 hours, and it is thus completely negligible. 18 The last term in Eq. (72) is of O(c−2 ). This term represents the phase ambiguity in the DOWR observable at t0 . A method dealing with this term was outlined in Ref. [6]. We denote this term as δn0 ≡ δn0 (t0 ) and keep it in the model. As a result, Eq. (72) can be presented in the following simplified form: n (vAB · nAB ) 4GMM o 1 + O(0.5 µm). (78) Rdowr (t) = dAB 1 − + 2 v2A + (nAB · vA )2 + v2B + (nAB · vB )2 + 2 2c 4c c (rA + rB ) Up to this point, we treated the start t0 of the integration interval in Eq. (65) as arbitrary. We now see that after negligible contributions are omitted, the start of the integration interval enters Eq. (78) only in the form of the definition of the phase ambiguity δn(t0 ). As we indicated above, dealing with this term is discussed in Ref. [6]. Once the effects of this phase ambiguity are accounted for, our formulation of the instantaneous DOWR observable, in the form of Eq. (78), becomes independent of the choice of the start of the integration interval, and thus t0 is truly arbitrary, even as we maintain an instantaneous range accuracy better than 1 µm, as needed for the GRAIL mission. From Fig. 1 we can see that the vectors RA and RB are given as RA = xEM + xM + yA and RB = xEM + xM + yB . (79) These vectors are measured simultaneously with the signal reception in TBD and are needed to compute Eq. (78). C. Dual One-Way Range-Rate (DOWRR) observables on GRAIL To develop an analytical form for the DOWRR observable, we use Eq. (59) to express it as ṅAB (tB ) = dτ dτ dt dnAB B A A0 = fB0 . − fA0 dtB dtB dtA dtB (80) Using the notation fB = fB0 (dτB /dtB ) and fA = fA0 (dτA /dtA ) for the coordinate frequencies of the two clocks, we can present Eq. (80) as ṅAB (tB ) = fB − fA dtA0 . dtB (81) As with the DOWR, the second observable DOWRR deals with the signal propagation from B0 to A is derived in an analogous way. Similarly to Eq. (81), at the time tA on the spacecraft A we have ṅBA (tA ) = fA − fB dtB0 , dtA (82) where the time of signal’s emission tB0 may be presented as a function of signal reception tA as tB0 = tB0 (tA ). Following the procedure outlined in Ref. [6], we assume that the LGRS clocks on both spacecraft are synchronized, such that tA = tB = t. We now can form a quantity that is called dual one-way range-rate (DOWRR): vdowrr (t) = c fA (dtA0 /dtB ) + fB (dtB0 /dtA ) ṅAB (t) + ṅBA (t) . =c 1− fA + fB fA + fB (83) The ratio of coordinate times dtA0 /dtB (and similarly dtB0 /dtA ) can be computed by differentiating the coordinate time transfer equation (60) for tB − tA0 = TAB (tA0 , tB ) (and similarly for tA − tB0 = TBA (tA , tB0 )) with respect to the reception time tB . This procedure was already performed in Appendix A 3, resulting in Eq. (A33) for dtA0 /dtB . From this equation, the ratio dtB0 /dtA is obtained by interchanging A ↔ B. Substituting these results for dtA0 /dtB and dtB0 /dtA into Eq. (83) we obtain the following expression for vdowrr : (fB aB − fA aA ) · dAB 1 (fB vB − fA vA ) · vAB + + vdowrr (t) = (nAB · vAB ) − c fA + fB fA + fB 1 fA (nAB · vA )(v AB · vA ) + fB (nAB · vB )(v AB · vB ) + 2 − c fA + fB 4GMM dAB (n · v ) + (n · v ) + O(c−2 ). (84) − A A B B c2 (rA + rB )2 19 The first term in Eq. (84) is the first order (∼ 1/c) Doppler term, which may be as high as 2 m/s (see Table I); it clearly must be kept in the model. The second term on the right hand side of Eq. (84) can be evaluated using Eq. (73) as (fB vB − fA vA ) · vAB (vA + vB ) · vAB v2AB fB − fA = + = (5.2 × 10−5 + 6 × 10−10 ) m/s. (85) c(fA + fB ) 2c fA + fB 2c Therefore, the second term in Eq. (85) can be dropped, but the v2AB /2c term must be kept in the model. The third term in Eq. (84) is the second order (∼ 1/c2 ) acceleration-dependent Sagnac effect. To evaluate this term, we note that the acceleration vectors of the spacecraft point in different directions due to the ∼200 km separation between the two craft. The vector difference can be calculated as aAB ≃ 0.17 m/s2 . We evaluate this term in a manner similar to Eq. (75) to obtain a magnitude of (fB aB − fA aA ) · dAB (aA + aB ) · dAB (aAB · dAB ) fB − fA = + = (6 × 10−5 + 4 × 10−11 ) m/s. (86) c(fA + fB ) 2c fA + fB 2c We see that the second term in Eq. (86) is less than the needed accuracy of 1 µm/s and it can be omitted; the (aAB · dAB )/2c term, however, must be kept in the model. The (1/c2 ) term on the second line of (84) can be presented as o 1 n 1 n fA (nAB · vA )(vAB · vA ) + fB (nAB · vB )(vAB · vB ) o = 2 (nAB · vA )(vAB · vA ) + (nAB · vB )(v AB · vB ) + 2 c f + fB 2c A o 1 fB − fA n + 2 (nAB · vB )(vAB · vB ) − (nAB · vA )(vAB · vA ) = (2 × 10−6 + 6 × 10−14 ) m/s. (87) 2c fB + fA Therefore, only the first of the two terms on the right-hand side of Eq. (87) must be retained. As a result, given the strict formation configuration implemented on the GRAIL mission, the model for DOWRR on GRAIL given by Eq. (84) has the following form: o o 1 n 1n 2 vAB + (aAB · dAB ) + 2 (nAB · vA )(v AB · vA ) + (nAB · vB )(vAB · vB ) − vdowrr (t) = (nAB · vAB ) − 2c 2c dAB 4GMM (n · v ) + (n · v ) + O(0.1 µm/s). (88) − A A B B c2 (rA + rB )2 Equation (88) represents the instantaneous DOWRR observable for the GRAIL mission, developed to a level of accuracy better than 1 µm/s. One can verify that the result given in Eq. (88) may be obtained directly from Eq. (78) by simply differentiating Eq. (78) with respect to time and retaining terms to the appropriate order. IV. CONCLUSIONS AND RECOMMENDATIONS We considered the formulation of a relativistic model for the observables of the GRAIL mission. We addressed some practical aspects of implementing the relevant computations. We derived an analytic expression that characterizes the process of forming the Ka-band ranging observables of GRAIL and developed a model for the dual one-way range (DOWR) observable. We also briefly addressed the transformation of relativistic gravitational potentials. This material can be used to improve the accuracy of modeling of the GRAIL fundamental observables. We presented a hierarchy of relativistic coordinate reference frames that are needed to GRAIL. In this respect, we introduced the barycentric (BCRS), geocentric (GCRS), topocentric (TCRS), lunicentric (LCRS) and spacecraft (SCRS) coordinate reference systems, together with the structure of the corresponding metric tensors in each of these systems and the form of the proper relativistic gravitational potentials—all presented at the accuracy required for GRAIL. We advocate a definition for the LCRS with its proper time, which we call the TCL. We presented the rules for transforming time and position measurements between the reference frames involved. The formula given by Eq. (78) is the main result of this paper. It is derived for the first time at this high level of accuracy including the terms of the 1/c2 order. The final expression (78) is relatively simple and easy to utilize in practice. The equations we provide for time and frequency transfers are accurate to the level of 1 µm when used to analyze GRAIL ranging data. Modeling the DOWR observable at this level of accuracy is the most important priority for the mission and must be taken into account for the science data analysis. Most of the relativistic computations for GRAIL are done implicitly and are based on the models and tools available within the framework of JPL’s Multiple Interferometric Ranging Analysis and GPS Ensemble (MIRAGE) software 20 [3]. General relativistic equations of motion form the “back-bone” of the entire suite of models in MIRAGE and rely on the formulation given in Ref. [18]. To navigate the GRAIL spacecraft, the code transforms the proper time of each of the GRAIL spacecraft to the time based on the SSB frame and integrates the spacecraft’s barycentric equations of motion. To determine the inter-spacecraft range, the code then iteratively solves the barycentric light-time equations in terms of instantaneous distance (by recomputing the transmitter’s position bearing in mind the elapsed light-time) in the presence of the Shapiro term. The analytical closed-form solution for DOWR (78) is not only more elegant, it allows for direct investigation of the observables and possible error terms under various circumstances in data analysis. We also developed a similarly accurate formulation for the DOWRR observable. Equation (88) allows us to calculate the value of this observable with an accuracy that is significantly better than 1 µm/s. In the form presented, Eqs. (78) and (88) can be readily incorporated into computer code used to model orbits and radio-metric observables. All the quantities in these equations are directly computable once the numerical positions and velocities of the spacecraft and Solar System bodies are known. It is also relatively straightforward to compute partial derivatives of these equations with respect to the the GCRS, which facilitates their use in fast numerical integration codes and optimizing solvers. Lastly, although we presented ideal, noise-free solutions, one can add relevant noise sources, including those in Ref. [6]. In a practical sense, the small relativistic terms that we calculated are easily absorbed into constant and periodic ad-hoc biases that are introduced during data analysis, with no impact whatsoever on mission objectives or the quality of the mission’s results. Yet the existence of these terms, and the fact that they are observable at the level of sensitivity of the GRAIL mission demonstrate that GRAIL is already a practical instrument for relativistic geodesy, especially after the mission is complete and all data sets for the primary and extended mission phases are assembled [4, 5]. For future spacecraft that operate at even greater accuracy, accounting for these relativistic terms will be essential. Although this paper was aimed specifically at discussing the range and range-rate observables of the GRAIL mission, we note that the solutions presented here are also applicable to other, similar missions. Foremost comes to mind the GRACE with a mission design very similar to that of GRAIL. Indeed, the calculations presented here may help shed light on the origin of small residual terms that were seen in the GRACE range and range rate observables [31]. We will further investigate this possibility with results to be reported elsewhere. Clearly, the model to be developed for the GRACE Follow-on mission [32, 33] must include similar higher-order terms to reach the anticipated DOWR and DOWRR at the level of few nm and nm/s correspondingly. We will address these issues in a subsequent publications. Acknowledgments We thank Sami W. Asmar, William M. Folkner, Nathaniel E. Harvey, Alexander S. Konopliv, Gerhard L. Kruizinga, Ryan S. Park, Michael M. Watkins, James G. Williams, Dah-Ning Yuan of JPL and Maria T. Zuber of MIT for their interest and support during the work and preparation of this manuscript. We also thank Sergey M. Kopeikin and Sergey A. Klioner for their insightful comments and suggestions. We also thank the anonymous referee for valuable comments on this manuscript. This work was performed at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. [1] J. R. Kim, Ph.D. thesis, University of Texas at Austin (2000). [2] R. Roncoli and K. Fujii, AIAA Guidance, Navigation, and Control Conference, Toronto, AIAA Paper 2010-8383 (2010). [3] R. S. Park, S. W. Asmar, E. G. Fahnestock, A. S. Konopliv, W. Lu, and M. M. Watkins, J. of Spacecraft and Rockets 49, 390 (2012). [4] M. T. Zuber, D. E. Smith, D. H. Lehman, T. L. Hoffman, S. W. Asmar, and M. M. Watkins, Space Sci. Rev., submitted (2012). [5] S. W. Asmar, A. S. Konopliv, S. R. Park, M. M. Watkins, G. Kruizinga, M. T. Zuber, D. E. Smith, J. G. Williams, M. Paik, D.-N. Yuan, et al., Space Sci. Rev., submitted (2012). [6] G. L. H. Kruizinga and W. I. Bertiger, GRAIL Project Memorandum (2009). [7] M. Soffel, S. A. Klioner, G. Petit, S. M. Kopeikin, P. Bretagnon, V. A. Brumberg, N. Capitaine, T. Damour, T. Fukushima, B. Guinot, et al., Astron. J. 126, 2687 (2003). [8] D. D. McCarthy, Tech. Rep., U.S. Naval Observatory (2010), IERS Conventions, URL http://www.iers.org/IERS/EN/Publications/TechnicalNotes/tn36.html. [9] S. Kopeikin and Y. Xie, Celestial Mechanics and Dynamical Astronomy 108, 245 (2010). [10] S. A. Klioner, Phys. Rev. D 69, 124001 (2004), arXiv:0311540 [astro-ph]. [11] S. A. Klioner, Astron. Astrophys. 478, 951 (2008). [12] S. A. Klioner, N. Capitaine, W. Folkner, B. Guinot, T. Y. Huang, S. Kopeikin, G. Petit, E. Pitjeva, P. K. Seidelmann, and M. Soffel, in Relativity in Fundamental Astronomy: Dynamics, Reference Frames, and Data Analysis, edited by S. Klioner, P. K. Seidelmann, and M. Soffel (Cambridge University Press, 2010). 21 [13] S. M. Kopeikin, M. Efroimsky, and G. Kaplan, Relativistic Celestial Mechanics of the Solar System (Wiley-VCH, 2011). [14] W. M. Folkner, J. G. Williams, and D. H. Boggs, IPN Progress Report 42-178, 1 (2009), http://ipnpr.jpl.nasa.gov/ progress report/42-178/178C.pdf. [15] S. G. Turyshev and V. T. Toth, in preparation (2012). [16] S. G. Turyshev, O. L. Minazzoli, and V. T. Toth, J. Math. Phys. 53, 032501 (2012), arXiv:1109.1796 [gr-qc]. [17] A. Einstein, L. Infeld, and B. Hoffmann, The Annals of Mathematics 39, 65 (1938). [18] T. D. Moyer, Formulation for Observed and Computed Values of Deep Space Network Data Types for Navigation, JPL Deep-Space Communications and Navigation Series (Wiley-Interscience, 2003). [19] O. Montenbruck and B. Gill, Satellite Orbits (Springer, 2005), 3rd ed. [20] S. G. Turyshev, Ann. Rev. Nucl. Part. Sci. 58, 207 (2008), arXiv:0806.1731 [gr-qc]. [21] E. M. Standish and J. G. Williams (Mill Valley: University Science Books, in press, 2012), Explanatory Supplement to the American Ephemeris and Nautical Almanac, P. K. Seidelmann, ed., chap. 8. [22] P. Teyssandier and C. Le Poncin-Lafitte, Classical and Quantum Gravity 25, 145020 (2008), arXiv:0803.0277 [gr-qc]. [23] N. Ashby and B. Bertotti, Classical and Quantum Gravity 27, 145013 (2010), arXiv:0912.2705 [gr-qc]. [24] C. M. Will, Theory and experiment in gravitational physics (Cambridge University Press, 2000), 2nd ed. [25] J. G. Williams, D. H. Boggs, C. F. Yoder, J. T. Ratcliff, and J. O. Dickey, J. Geophys. Res. 106, 27933 (2001). [26] S. Goossens and K. Matsumoto, Geophys. Res. Lett. 35, L02204 (2008). [27] G. Tommei, A. Milani, and D. Vokrouhlicky, Celest. Mech. Dynam. Astron. 107, 285 (2010). [28] L. D. Landau and E. M. Lifshitz, The Classical Theory of Fields (in Russian) (Nauka, Moscow, 1988), 7th ed. [29] C. W. Misner, K. S. Thorne, and J. A. Wheeler, Gravitation (W. H. Freeman & Co. (San Francisco), 1973). [30] L. Blanchet, C. Salomon, P. Teyssandier, and P. Wolf, Astron. Astrophys. 370, 320 (2001). [31] W. Bertiger, Y. Bar-Sever, S. Bettadpur, S. Desai, C. Dunn, B. Haines, G. Kruizinga, D. Kuang, S. Nandi, L. Romans, et al., Proceedings of ION GPS 2002, Portland OR (2002). [32] B. Loomis, Ph.D. thesis, University of Texas at Austin (2009). [33] M. M. Watkins, F. Flechtner, P. Morton, M. A. Gross, and S. V. Bettadpur, AGU Fall Meeting Abstracts p. A7 (2011). [34] S. M. Kopeikin, Mon. Not. R. Astron. Soc. 399, 1539 (2009). [35] M. V. Sazhin, I. Y. Vlasov, O. S. Sazhina, and V. G. Turyshev, Astronomy Reports 54, 959 (2010). Appendix A: The phase of an electromagnetic signal in gravitational field In this Appendix, we present derivations of the formulas use for time and frequency transfer in the GRAIL experiment. The derivation presented in this appendix is based on material that can be found in standard textbooks such as Refs. [29] and [24]. A general solution is presented to the problem of light propagation in a gravitational field in the linearized approximation. 1. General-relativistic post-Minkowskian space-time To develop the solution to the equations of the general theory of relativity in the post-Minkowskian approximation, we introduce the post-Minkowskian decomposition of the metric tensor gmn as gmn = γmn + hmn + O(G2 ), (A1) where hmn denotes the post-Minkowskian perturbation of the Minkowski metric tensor γmn . Following [29], we impose the harmonic gauge condition on the metric tensor gmn , given in the form √ ∂m (A2) or ∂m hmn − 21 ∂ n h = O(G2 ), −gg mn = 0, where h = γkl hkl + O(G2 ). In the first post-Minkowskian approximation of general relativity [20, 29], Einstein’s equations Rmn = 16πG/c4 T mn − 12 g mn T take the following form in arbitrary harmonic coordinates {xm } = (ct, x): 16πG mn 1 mn 1 ∂2 2 − ∇ hmn = T − 2 γ T + O(G2 ), c2 ∂t2 c4 (A3) where T mn is the stress-energy tensor describing a body that deflects a light ray and T = γkl T kl + O(G). In the linearized approximation and neglecting the higher multipole moments, this tensor is given as [28]: p T mn (t, x) = M um un 1 − β 2 δ 3 (x − z(t)) + O(G), (A4) 22 where M is the rest mass of the body, z(t) its time-dependent spatial coordinate, β = c−1 dz/dt and δ 3 (x) is the three-dimensional Dirac delta function. The body’s normalized four-velocity um , such that um um = 1, is given by vα v α /c dt 1 + O(G), uα = p + O(G). (A5) , u0 = um = (u0 , uα ) = u0 1, = p c ds 1 − β2 1 − β2 √ Note that in Eq. (A4) we have neglected the factor of −g. This is done because in the linearized approximation √ 2 −g = 1 + O(G) and the quadratic terms ∝ G are irrelevant in T mn since the corresponding time-dependent terms of the second post-Minkowskian order are currently unobservable in measurements made in the Solar System. We can now write down the Green’s function solution to Eqs. (A3): Z p 4GM hmn (t, x) = 1 − β ′2 u′m u′n − 21 γ mn δ 3 (x′ − z(t′ ))G(t, x; t′ , x′ )d3 x′ dt′ + O(G2 ), (A6) 2 c where G(t, x; t′ , x′ ) is the Green’s function G(t, x; t′ , x′ ) = G(t − t′ ; x − x′ ) = 1 1 1 ′ ′ δ t − t − |x − x | . 4π |x − x′ | c (A7) Integrating Eq. (A6), one obtains the post-Minkowskian metric tensor perturbation in terms of retarded LiénardWiechert tensor potentials [29, 34]: hmn (t, x) = 4GM um un − 21 γ mn + O(G), c2 um r m (A8) where rm = xm − z m (tret ) = [c(t − tret ), (x − z(tret )]. In Eq. (A8), all time-dependent quantities are taken at a retarded time tret (defined by Eq. (A10) below): um ≡ um (tret ) = c−1 dz m (tret )/dtret is the body’s four-velocity, β(tret ) = c−1 dz(tret )/dtret is the body’s coordinate velocity normalized to the speed of light c. In the solution to Eq. (A3), given in terms of the retarded Liénard-Wiechert potentials, all quantities involved including the distance rm = xm − z m (tret ), the body’s world-line z m (tret ) = [ctret , z(tret )], and the four-velocity um (tret ) are functions of the retarded time tret . It is known (see, for instance, [34]) that the retarded time in the first post-Newtonian approximation may be found from the null-cone equation γmn rm rn ≡ γmn [xm − z m (tret )][xn − z n (tret )] = 0, (A9) suggesting that the retarded time tret = tret (t, x) is established as a solution to the equation 1 tret = t − |x − z(tret )|. c (A10) Note that Eq. (A10) has an analytic solution only in the case of uniform motion of the gravitating body along a straight line. 2. Phase of the electromagnetic wave The phase of an electromagnetic wave is a scalar function that is invariant under a set of general coordinate transformations. In the geometric optics approximation, the phase is found as a solution to the eikonal equation [28, 29, 34]: g mn ∂m ϕ∂n ϕ = 0, (A11) with g mn = γ mn − hmn + O(G2 ). Equation (A11) is a direct consequence of Maxwell’s equations. Its solution describes the front of an electromagnetic wave propagating in curved space-time. The solution’s geometric properties are defined by the metric tensor (A1), where hmn (A8) is the solution of the linearized Einstein equations (A3) with stress-energy tensor (A4). To solve Eq. (A11), we introduce a covector of the electromagnetic wavefront in curved space-time, Km = ∂m ϕ. We use λ to denote an affine parameter along the trajectory of a light ray being orthogonal to the wavefront ϕ. Vector K m = dxm /dλ = g mn ∂n ϕ is tangent to the light ray. Equation (A11) states that the vector K m simply is null or gmn K m K n = 0. Therefore, the light rays are null geodesics [28] described by 1 dKm = ∂m gkl K k K l . dλ 2 (A12) 23 Since eikonal and light-ray equations, given by Eqs. (A11) and (A12) respectively, have equivalent physical content in the general theory of relativity, one can use either of them to study the properties of an electromagnetic wave. However, the eikonal equation offers a more straightforward way to study the propagation of a wave. To find a solution of Eq. (A11), we expand the eikonal ϕ with respect to the gravitational constant G assuming that the unperturbed solution of Eq. (A11) is a plane wave. The expansion may be given as Z ϕ(t, x) = ϕ0 + km dxm + ϕG (t, x) + O(G2 ), (A13) where ϕ0 is an integration constant and k m = k 0 (1, k) is a constant null vector (i.e., γmn k m k n = 0) along the direction of propagation of the unperturbed electromagnetic wavefront. Furthermore, k 0 = ω/c where ω is the constant frequency of the unperturbed wave, and ϕG is the perturbation of the eikonal to O(G), which is yet to be determined. Substituting Eqs. (A8) and (A12) into (A11) and keeping only first order terms in G, we obtain an ordinary differential equation to determine ϕG : dϕG 1 2GM (km um )2 + O(G2 ), = hmn km kn = dλ 2 c2 um r m (A14) which alternatively can be obtained as a first integral of the null geodesic equation (A12). We can now integrate (A14) while keeping in mind that dk m = 0 and employing the exact relationship [34]: ds 1 dλ = = d ln [km rm ] . m m um r km r km um (A15) Neglecting the body’s acceleration (or dum = 0), a plane-wave solution of Eq. (A14) has the form ϕG (t, x) = 2GM (km um ) ln [km rm ] , c2 (A16) where all quantities on the right-hand side are taken at the retarded instant of time tret in agreement with (A10). Therefore, we can now write the post-Minkowskian expansion for the phase of the electromagnetic wave as: Z 2GM ϕ(t, x) = ϕ0 + km dxm + (km um ) ln [km rm ] + O(G2 ), (A17) c2 which can be presented in the following form 2GM 1 − (v · k)/c 0 p ϕ(t, x) = ϕ0 + k0 ct − k · x + ln k (r − k · r) + O(G2 ), 2 c 1 − v 2 /c2 (A18) where all the quantities in the last term are taken at the retarded time tret defined by Eq. (A10). Let us now consider signal propagation from a point (ctA , xA ) to a point (ct, x). Then along the signal’s path the phase (A18) will change according to 2GM 1 − (v · k)/c h r − k · r i p + O(G2 ), (A19) ϕ(t, x) = k0 ct − RA + ln c2 rA − k · rA 1 − v 2 /c2 where we used the following notations: k = RA /RA , RA = x − xA , RA = |RA | and r = x − z(tret ), r = |r|. One can further simplify the argument of the logarithmic term in Eq. (A19) as r−k·r rA + r − RA = . rA − k · rA rA + r + RA Thus, Eq. (A19) takes the form: 2GM 1 − (v · k)/c h rA + r + RA i p ϕ(t, x) = k0 ct − RA − + O(G2 ). ln 2 2 c2 r + r − R 1 − v /c A A (A20) (A21) Effects of order of vG/c3 are very small and one can neglect them in Eq. (A21). As a result, the post-Minkowskian phase of the electromagnetic wave, with an accuracy appropriate for modern-day Solar System experiments [30, 34], can be presented as: 2GM h rA + r + RA i + O(G2 , c−3 ). (A22) ln ϕ(t, x) = k0 ct − RA − c2 rA + r − RA Along the signal’s path the phase stays constant and equal to ϕ(tA , xA ) = k0 ctA . 24 3. Coordinate gravitational time delay Consider the case of one-way signal transmission. Let A be the emitting station, with BCRS position xA (t), and B the receiving station, with position xB (t). Also, z(t) is the vector connecting the SSB to a gravitating body. We denote by tA the coordinate time at the instant of emission of a radio signal, and by tB the coordinate time at the instant of reception. We put rA = xA (tA ) − z(tA ), rB = xB (tB ) − z(tB ), also RAB = xB (tB ) − xA (tA ), NAB = RAB /RAB and RAB = |RAB |, rA = |rA |, rB = |rB | are the Euclidean norms of these vectors. We know that along the signal’s path the phase stays constant. Thus, equating the eikonal of the wave given by Eq. (A22) at the two points A and B as ϕ(tA , xA ) = ϕ(tB , xB ), we determine the gravitational delay of the signal moving through a stationary space-time. Indeed, up to O(c−3 ), the coordinate time transfer TAB = tB − tA is given by [15, 30]: TAB = tB − tA = 2GM h rA + rB + RAB i RAB , ln + c c3 rA + rB − RAB (A23) where the logarithmic term represents the Shapiro time delay. The ratio of coordinate times dtA0 /dtB (and similarly dtB0 /dtA ) can be determined directly by differentiating the coordinate time transfer equation (60) for tB − tA0 = TAB (tA0 , tB ) (and similarly for tA − tB0 = TBA (tA , tB0 )) with respect to the reception time tB . In other words, we must evaluate d d n1 2GM rA + rB + RAB o . (A24) (tB − tA0 ) = ln |rB (tB ) − rA (tA0 )| + dtB dtB c c3 rA + rB − RAB While performing the differentiation, we must account for the fact that the coordinate distance between A0 and B depends on both the times of emission and reception, i.e. RAB = |rB (tB ) − rA (tA0 )|. For instance, we have dRAB dtA0 = NAB · vB − vA . dtB dtB (A25) qB dtA0 = , dtB qA0 (A26) For notational convenience, we henceforth denote the vectors by rA = xA (tA ) − z(tA ) and rB = xB (tB ) − z(tB ). As before, we have rA = |rA | and rB = |rB |, as well as the coordinate velocities vA = ẋA (tA ) − ż(tA ) and vB = ẋB (tB ) − ż(tB ). Performing the differentiation in Eq. (A24) to order 1/c3 , the last factor in Eq. (C1) dtA0 /dtB can be given by the ratio (in agreement with the recent work on the ACES [30] and RadioAstron [35] missions; for convenience, we adopt notations similar to those introduced in Ref. [30]): where qA0 and qB are derived from 1 4GM (rA + rB )(NAB · vA ) + RAB (rA · vA )/rA qA0 = 1 − (NAB · vA ) − , 2 c c3 (rA + rB )2 − RAB 4GM (rA + rB )(NAB · vB ) − RAB (rB · vB )/rB 1 . qB = 1 − (NAB · vB ) − 2 c c3 (rA + rB )2 − RAB (A27) (A28) In an experiment, the position of the transmitter A0 may be recorded at the time of reception tB rather than at the time of emission tA0 , i.e. we may have more direct access to xA (tB ) rather than xA (tA0 ), and the formulae (A27)–(A28) get modified by Sagnac correction terms consistently to the order 1/c3 : 1 1 qA0 = 1 − (nAB · vA ) − 2 v2A − (nAB · vA )2 − (aA · dAB ) + c c o 1 n 2 + 3 (vA − (nAB · vA )2 )(nAB · vA ) + dAB 3(aA · vA ) − (nAB · aA )(nAB · vA ) − (ȧA · dAB ) − 2c 4GM (rA + rB )(nAB · vA ) + dAB (rA · vA )/rA − + O(c−4 ), (A29) c3 (rA + rB )2 − d2AB 1 1 qB = 1 − (nAB · vB ) − 2 (vA · vB ) − (nAB · vA )(nAB · vB ) + c c o 1 n 2 + 3 (vA − (nAB · vA )2 )(nAB · vB ) + dAB (aA · vB ) − (nAB · aA )(nAB · vB ) − 2c 25 − 4GM (rA + rB )(nAB · vB ) − dAB (rB · vB )/rB + O(c−4 ). c3 (rA + rB )2 − d2AB (A30) For the GRAIL spacecraft in lunar orbit, in the lunicentric reference frame the first, 1/c term in Eqs. (A29)–(A30) is ∼ 10−4 in the barycentric frame and ∼ 5.5 × 10−6 in the lunicentric frame. The second term is the 1/c2 Sagnac term, which is ∼ 10−8 for the BCRS and ∼ 10−11 for the LCRS. The 1/c3 terms are ∼ 8 × 10−17 for the Sagnac term and ∼ 1.4 × 10−16 for the Shapiro term. Note that the result given in Eqs. (A29)–(A30) has been obtained assuming that the field of the Earth is spherically symmetric. Indeed, the J2 -terms in the factor (qA /qB ) do not exceed 10−15 . As a result, the expression dtA0 /dtB from Eq. (A26) has the form dtA0 qB 1 1 = = 1 − (nAB · vAB ) − 2 (vA · vAB ) + (aA · dAB ) − dtB qA0 c c 1 n − 3 2(vA · vAB )(nAB · vA ) + (v2A − (nAB · vA )2 )(nAB · vAB ) + 2c o + dAB 2(aA · vA ) + 2(nAB · aA )(nAB · vA ) − (aA · vAB ) − (nAB · aA )(nAB · vAB ) − (ȧA · dAB ) − 4GM (rA + rB )(nAB · vAB ) − dAB (rB · vB )/rB ) + (rA · vA )/rA ) − 3 + O(c−3 ). (A31) c (rA + rB )2 − d2AB We can evaluate the magnitude of each of the seven terms in Eq. (A31) using the values from Table I. We will use these basic formation parameters to evaluate the terms in (84), keeping only those that contribute to the range-rate model more than ∆vdowrr = 1 µm/s. The error terms in these expressions must be less than ∆vdowrr /c = 3 × 10−15 . Keeping these numbers in mind we see that the (1/c) term in (A31) is of the order 7 × 10−9 and must be kept in the model. The two 1/c2 terms were evaluated to be of the order of (vA · vAB ) + (aA · dAB ) /c2 ≈ 6 × 10−11 + 3 × 10−12. Thus, both of these terms must be in the model. Among the 1/c3 terms, the first one is ∼ 6 × 10−15 and must be kept. Next is the term that contains (nAB · vAB ), which was evaluated to be at most 3 × 10−17 and is too small to be in the model. The remaining 1/c3 term, which is ∝ dAB , is at most ∼ 5 × 10−16 and thus, this entire term may be neglected. Lastly, the Shapiro term in Eq. (A31) was evaluated to be 4GM (nAB · vAB ) 4GM dAB (nA · vA ) + (nB · vB ) − ≈ 8 × 10−19 + 3 × 10−15 . (A32) + c3 (rA + rB ) c3 (rA + rB )2 Therefore, we will keep only the second part of the Shapiro term in the model. As a result, the expression (A31) may be presented in the following simplified form: 1 1 dtA0 = 1 − (nAB · vAB ) − 2 (vA · vAB ) + (aA · dAB ) − dtB c c 4GM 1 dAB − 3 (nAB · vA )(vAB · vA ) + (nA · vA ) + (nB · vB ) + O(5 × 10−16 ), 3 2 c c (rA + rB ) (A33) where nA = rA /rA and nB = rB /rB . This expression accounts for all the terms that may have magnitudes larger than 5 × 10−16 and are thus relevant to the GRAIL mission configuration. Appendix B: Evaluating the DOWR integral expression In Sec. III B we defined the quantities ǫAB (tB ) and ǫBA (tA ), which were given in Eq. (68). We now construct the quantity ǫAB (tB ) + ǫBA (tA ): 1n ǫAB (tB ) + ǫBA (tA ) = 2 fB0 uB (tB )(tB − t0B ) − uB (tA ) tA − TBA (tA ) − t0A − TBA (t0A ) − c Z tB −tA +TBA (tA ) − uB (t′ )dt′ + t0 −t0 +TBA (t0 ) B A A 0 + fA0 uA (tA )(tA − tA ) − uA (tB ) tB − TAB (tB ) − t0B − TAB (t0B ) − Z tA −tB +TAB (tB ) o + O(c−4 ). uA (t′ )dt′ − t0A −t0B +TAB (t0B ) (B1) 26 To estimate the magnitude of this quantity, we can assume that the clocks A and B are perfectly synchronized. In reality, nothing is perfect and the GRAIL mission’s design relies on a synchronization procedure (see Ref. [6] for details) that, of course, leaves synchronization errors that percolate in the data analysis. One can study the impact of synchronization errors on the measurement accuracy, using the approach outlined here. However, presently we concerned with the evaluation of the magnitude of the error that may arise in the ideal case, if we were to drop the ǫAB and ǫBA terms. Thus, we assume that tA = tB = t and t0A = t0B = t0 , so that the quantity given by Eq. (B1) is reduced to Z TAB (t) 1n uA (t′ )dt′ + ǫAB (t) + ǫBA (t) = 2 fA0 uA (t) TAB (t) − TAB (t0 ) − c TAB (t0 ) Z TBA (t) o + O(c−4 ). (B2) uB (t′ )dt′ + fB0 uB (t) TBA (t) − TBA (t0 ) − TBA (t0 ) We series expand the integrals in Eq. (B2) around the start t0 of the integration interal. For the integral term multiplied by fA0 inside the second set of square brackets in Eq. (B2), we have uA (t′ ) = uA (t0 )+u̇A (t0 )(t′ −t0 )+O(∆t2 ) (where ∆t = t − t0 ), and obtain: uA (t) TAB (t) − TAB (t0 ) − Z TAB (t) TAB (t0 ) 2 uA (t′ )dt′ = 12 u̇A (t0 ) TAB (t) − TAB (t0 ) + O(∆t3 ). (B3) The integral term in the first set of square brackets in Eq. (B2), which is multiplied by fB0 , can be evaluated in a similar way. Keeping just the leading terms (∼ dAB /c) in TAB and TBA , we present, for instance, TAB (t) − TAB (t0 ) = c−1 d˙AB (t − t0 ) + O(c−2 ), so that Eq. (B2) becomes: ǫAB (t) + ǫBA (t) = 1 fA0 u̇A (t0 ) + fB0 u̇B (t0 ) d˙2AB (t − t0 )2 + O(c−4 ) = O(c−4 ). 4 2c (B4) Remembering the form of uB (t) in Eq. (64), and denoting by f0 the typical GRAIL radio frequency (fA0 ≃ fB0 ≃ f0 , see Table I), we see that the term above is of the order of ∼ (vA aA d˙2AB /2c4 )(t − t0 )2 f0 ≪ 7.7 × 10−24 (t − t0 )2 f0 at most. This term is multiplied by ∼ c/2f0 as its contribution to the DOWR measurement is calculated in Eq. (72): 2 the magnitude of this contribution is therefore less than 1 × 10−15 (t − t0 )2 m/s , which is negligible for GRAIL. Appendix C: Relativistic frequency transfer between spacecraft and a DSN antenna The DOWR observable is defined in terms of the coordinate frequencies fA and fB of both signals generated on-board the two GRAIL spacecraft. However, these frequencies are measured by a ground-based DSN station. Specifically, in the case of the GRAIL spacecraft, while the spacecraft-to-ground transmission takes place using frequency bands that are different from the frequencies used for inter-spacecraft communication, the two frequencies are synthesized using the same timing source (USO) on board. Thus, measuring the frequency of the spacecraft-to-ground transmission is, in effect, a measurement of the inter-spacecraft frequency as well. Below, we simply assume that the spacecraft-to-ground transmission does serve as a means to measure fA and fB without going into further detail. We discuss how one can use these measurements, taken using the TT time coordinate, and transform them to the TDB time coordinate. We focus only on relativistic frequency transformations between the frames involved and will neglect frequency-dependent media effects, which are easy to reinstate when needed. Although communication between the spacecraft is conducted using Ka-band microwave signals and spacecraft-toDSN is done relying on X-band signals, these frequencies are related by simple numeric factors. Therefore, in a slight abuse of notation, we use the same symbol for the inter-spacecraft and spacecraft-to-ground frequencies. 1. One-way frequency transfer We consider the situation when a signal with frequency f is measured by an electronic counter whose register is incremented by 1 each time the magnitude of the signal changes from minus to plus. The number of cycles dn measured by this counter in the interval of proper time dτ is then dn = f dτ . Thus, the frequency transfer between transmitter (at point A) and receiver (at point C) requires the determination of the ratio fAC /fA0 between the proper frequencies fA0 transmitted by satellite (A) and fAC on the ground (C). The infinitesimal proper time intervals dτA 27 and dτC correspond to the infinitesimal number of cycles, dn, at the transmission and reception points A and C, so that dnC A = dnA . Therefore, the one-way frequency shift during the transfer from A to C is (dτ /dt)A dtA dnC dτA fAC A dτA = = = . fA0 dτC dnA dτC (dτ /dt)C dtC (C1) For (dτ /dt)C /(dτ /dt)A , we get 1 − c−2 U (rA ) + 12 v2A (dτ /dt)A + O(c−4 ), (C2) = (dτ /dt)C 1 − c−2 U (rC ) + 21 v2C P P where U (rA ) = b Ub (rbM + yA ) and U (rC ) = b Ub (rbE + yC ) are the Newtonian potentials at the points A and C and vectors vA and vC are the barycentric velocities of the spacecraft and the DSN station correspondingly. The factor given by Eq. (C2) consists of the Einstein gravitational red-shift and second-order Doppler effects, both of order 1/c2 . Therefore, Eq. (C1) becomes 1 − c−2 U (rA ) + 12 v2A dtA fAC dτA = = + O(c−4 ). (C3) fA0 dτC 1 − c−2 U (rC ) + 21 v2C dtC The dtA /dtC term on the right-hand-side of Eq. (C3) is the ratio of coordinate periods of the same signal at A and C. It contains both the (∼ 1/c) Doppler effect and the (∼ 1/c3 ) terms that we seek. To compute this factor at the accuracy of 3 × 10−15 , it is sufficient to treat the Earth and Moon gravitational potentials as monopole potentials, i.e., the approximation U = GM/r and consequently, the use of Eq. (44) is sufficient. Similarly to the discussion in Sec.A 3, the ratio dtA /dtC can be computed by a direct differentiation of the coordinate time transfer TAC = tC − tA with respect to the emission time tA . In Appendix A 3 we computed all the relevant terms, accurate to the ∼ 1/c3 order. Using the time of reception tC (and expressing the time of emission as a function of the time of reception, tA = tA (tC )), the factor dtA /dtC in Eq. (C1) can be given (similar to Eq. (A26)) by qC dtA = dtC qA (C4) where qA and qC are then given to O(c−3 ) from Eqs. (A29)–(A30) as 1 1 qA = 1 − (nAC · vA ) − 2 v2A − (vA · nAC )2 − (aA · dAC ) + O(c−3 ), (C5) c c 1 1 qC = 1 − (nAC · vC ) − 2 (vA · vC ) − (vA · nAC )(vC · nAC ) + O(c−3 ), (C6) c c where dAC = xC (tC ) − xA (tC ) is the coordinate distance between A and C at the moment of reception at C (we have dAC = |dAC | and nAC = dAC /dAC ), where vA (tC ) denotes the coordinate velocity of the station A at that instant, and where aC is the acceleration of A. Therefore, the ratio qC /qA can be given as qC 1 1 = 1 − (nAC · vAC ) + 2 21 v2A − 21 v2C + 12 v2AC − (aA · dAC ) + O(c−3 ), (C7) qA c c where all the quantities involved are given at the time of reception tC . Therefore, for the one-way frequency transfer, we have 1 − c−2 U (rA ) + 21 v2A qC fAC · = , fA0 1 − c−2 U (rC ) + 21 v2C qA with the ratio qA /qC given to sufficient accuracy by Eq. (C7). Up to O(c−3 ), this expression becomes: X fAC 1 1 = 1 − (nAC · vAC ) + 2 21 v2AC + Ub (rbE + yC ) − Ub (rbM + yA ) − (aA · dAC ) + O(5 × 10−14 ). fA0 c c (C8) (C9) b Note that all the quantities in Eq. (C9) are taken at the reception time. This relation determines the frequency displacement due to the combined motion of the emitter and the receiver (intrinsic Doppler effect) and the difference in the gravitational potentials at the points of emission and reception of the signal (gravitational displacement of the frequency). In the case of GRAIL, the one-way frequency transfer given by Eqs. (C5)–(C6) can be evaluated numerically as follows. The first-order Doppler effect is |(nAC · vA )/c| = 5.5 × 10−6 ; for the ground |(nAC · vC )/c| = 1.6 × 10−6 . The second-order Doppler effect is v2A /2c2 = 3.4 × 10−10 for the satellite; v2C /2c2 = 1.3 × 10−12 for the ground. The gravitational red-shift (Einstein) effect is given by UA /c2 = UE (rA )/c2 = 6.5 × 10−10 ; UC /c2 = 6.9 × 10−10 . The O(c−3 ) terms are less than 3.6 × 10−14 for the spacecraft and 2.2 × 10−15 for the Earth station, thus, they are omitted. 28 2. Integrated Doppler Effect To measure the frequency of the received signal, DSN receivers count the number of cycles received from the spacecraft over intervals of time measured by high precision clocks located at the receiver station. We consider the integrated (one-way) Doppler effect that is at the basis of forming the Doppler observable for spacecraft such as GRAIL that are equipped with a precision on-board frequency source. Consider a clock (Fig. 2) with proper frequency fA0 , located at point A (GRAIL-A), which emits light signals during the proper time interval (τA1 , τA2 ). These signals are received with frequency fAC by a ground-based DSN station, located at C (see Fig. 1), over a proper time interval (τC1 , τC2 ). The measured quantity over the interval (τC1 , τC2 ) is the number of cycles received at point C from the transmitter at point A. The number of cycles received at point C must equal the number of cycles emitted at point A. Hence, Z τA2 Z τC2 C fA0 dτA , (C10) fA dτC = τA1 τC1 where the frequency fA0 is the constant (with respect to τA ) precision oscillator frequency of the transmitter, and may be moved outside the integral sign. The remaining integral can be transformed to the proper time τC with the use of Eq. (C3): Z τC2 C Z t2 C Z τC2 Z τA2 fA fA dτC dτA dτC = dτC = dt, (C11) dτA = dτ f f C A0 A0 dt τC1 t1 τC1 τA1 with ratio dτC /dt given by Eq. (25) as i X 1h dτC = 1 − 2 12 v2C + Ub (rbE + yC ) + O(10−17 ). dt c (C12) b Also, in Eq. (C9) we expressed the ratio fAC /fA0 at the adopted accuracy using terms that are functions of t. Therefore, Z τA2 1 Z t2 h i X 1 1 2 vAC − 12 v2C − Ub (rbM +yA )−(aA ·dAC ) dt+O(c−3 ), (C13) dτA = t2 −t1 − dAC (t2 )−dAC (t1 ) + 2 2 c c t1 τA1 b where we relied on the identity (nAC · vAC ) = d˙AC . Using these intermediate results, we find that the total number of cycles NAC , received during the count interval ∆t = t2 − t1 , is given by: Z τC2 1 fAC dτC = fA0 (t2 − t1 ) − fA0 dAC (t2 ) − dAC (t1 ) + NAC = c τC1 Z t2 h i X 1 1 2 + 2 fA0 vAC − 12 v2C − Ub (rbM + yA ) − (aA · dAC ) dt + O(c−3 ). (C14) 2 c t1 b Dividing the number of cycles NAC by the count interval ∆t yields the Doppler observable fˆAC : NC fˆAC = A . ∆t (C15) Therefore, the spacecraft frequency fˆA that is measured at a DSN receiver during the count interval of proper time ∆t can be modeled as Z t2 h i X dAC (t2 ) − dAC (t1 ) 1 fˆAC 1 2 = 1− + 2 vAC − 12 v2C − Ub (rbM + yA ) − (aA · dAC ) dt + O(c−3 ), (C16) 2 fA0 c∆t c ∆t t1 b where the range dAC is defined as dAC = |xC − xA |. 3. Transmitted spacecraft frequency, as seen from the BCRS In a 1-way Doppler transmission between a GRAIL spacecraft and a DSN station, the measured frequency received at the DSN station, fˆAC , is related to the transmitted frequency fA0 as given in Eq. (C16). Our objective is to express 29 the frequency received at the DSN using the TDB time coordinate of the BCRS. To do this, we imagine a hypothetical receiver that is at rest with respect to the BCRS, and co-located with the DSN receiver at the moment of receiving the same signal. Using Eq. (C9) we can relate the instantaneous frequency fA that is received by this hypothetical receiver to the transmitted frequency fA0 : X fA 1 1 Ub (zb ) − Ub (rbM + yA ) + (aA · dA ) + O(c−3 ), = 1 − (nA · vA ) + 2 12 v2A + fA0 c c (C17) b where P b Ub (zb ) is the gravitational potential at the origin of the SSB (see discussion after Eqs. (4)–(5)) defined as X b µ∗b zb = 0 and n X GMc o 1 µ∗b = GMb 1 + 2 vb2 − + O(c−4 ), 2c rbc (C18) c6=b also, zb is the barycentic position vector of the body b, vb = |żb | is its barycentic velocity, and all times are measured in the TDB time. Analogously to the discussion in the Section C 2, the Doppler observable fˆA that corresponds to the instantaneous frequency fA is established during the count interval of TDB time of ∆t = t2 − t1 and is described as Z t2 h i X fˆA dA (t2 ) − dA (t1 ) 1 1 2 = 1− + 2 v + U (z ) − U (r + y ) + (a · d ) dt + O(c−3 ), (C19) b b b bM A A A fA0 c∆t c ∆t t1 2 A b where we used the fact that dA = |xA | and d˙A = (nA · vA ). Finally, with the help of Eqs. (C16) and (C19), we obtain the following relation between two Doppler observables— the spacecraft’s frequency fˆA reported at the BCRS and the same frequency fˆAC as measured at the DSN: h fˆ ih fˆC i−1 A A + O(c−3 ) fˆA = fˆAC fA0 fA0 n 1 = fˆAC 1 + [dAC (t2 ) − dAC (t1 )] − [dA (t2 ) − dA (t1 )] + c∆t o X Z t2 1 Ub (zb )dt + O(∆t−2 , c−3 ). (vA · dC )(t2 ) − (vA · dC )(t1 ) + + 2 c ∆t t1 (C20) b Note that the vectors RAC (given in Eq. (2)), RA and RB , (given in Eq. (79)) are measured simultaneously at the time of signal reception. The result of Eq. (C20) is what we need in order to derive the DOWR and DWORR observables of the GRAIL mission given by Eqs. (78) and (86). With the X-band navigation on GRAIL with frequency ∼ 8.4 MHz, the result is accurate to ∼0.5 mHz, which is sufficient for the purposes of navigating the GRAIL twins around the Moon. The expression for frequency transformation, Eq. (C20), relates the frequency of the signal transmitted by the GRAIL spacecraft as reported at the BCRS to that measured by the DSN station. This result is complete up to the 1/c2 order, sufficient for GRAIL.

RELATED PAPERS

RELATED TOPICS

Log In

General relativistic observables of the GRAIL mission

General relativistic observables of the GRAIL mission

Related Papers

RELATED PAPERS

RELATED TOPICS