High Performance Clocks and Gravity Field Determination

TL;DR

Advances in high-precision atomic clocks, reaching 10^{-18} accuracy, are transforming geodesy by enabling new methods for defining height systems, measuring gravity fields, and improving space-borne gravimetry.

Contribution

This paper discusses how ultra-precise clocks can redefine geodetic reference frames, improve gravity field measurements, and address current limitations in space-based gravimetry.

Findings

Clocks with 10^{-18} accuracy can measure height differences of about 1 cm.

Clock data can resolve discrepancies in traditional height systems.

Potential for clocks to enhance regional gravity field recovery.

Abstract

Time measured by an ideal clock crucially depends on the gravitational potential and velocity of the clock according to general relativity. Technological advances in manufacturing high-precision atomic clocks have rapidly improved their accuracy and stability over the last decade that approached the level of 10$^{-18}$. Based on a fully relativistic description of the background gravitational physics, we discuss the impact of those highly-precise clocks on the realization of reference frames and time scales used in geodesy. We discuss the current definitions of basic geodetic concepts and come to the conclusion that the advances in clocks and other metrological technologies will soon require the re-definition of time scales or, at least, clarification to ensure their continuity and consistent use in practice. The relative frequency shift between two clocks is directly related to the difference in the values of the gravity potential at the points of clock's localization. According to general relativity the relative accuracy of clocks in 10$^{-18}$ is equivalent to measuring the gravitational red shift effect between two clocks with the height difference amounting to 1 cm. We show how clock measurements can provide geopotential numbers for the realization of gravity-field-related height systems and can resolve discrepancies in classically-determined height systems as well as between national height systems. Another application of clocks is the direct use of observed potential differences for the improved recovery of regional gravity field solutions. Finally, clock measurements for space-borne gravimetry are analyzed along with closely-related deficiencies of this method like an extra-ordinary knowledge of the spacecraft velocity, etc. For all these applications besides the near-future prospects, we also discuss the challenges that are related to using those novel clock data in geodesy.

Figures

Equations

$$g_{00}$$

$$g_{0i}$$

$$g_{ij}$$

$$w(t,\mathbf{x})$$

$$w^i(t,\mathbf{x})$$

$$V(T,\mathbf{X})$$

$$V^i(T,\mathbf{X})$$

$$w_E(t,\mathbf{x})$$

$$w_E^i(t,\mathbf{x})$$

$$V_E(T,\mathbf{X})$$

$$V_E^i(T,\mathbf{X})$$

$$\Sigma (T,\mathbf{X})=\sigma (t,\mathbf{x})=\rho (t,\mathbf{x})+\mathcal{O}\left(c^{-2}\right),$$

$${\Sigma }^i(T,\mathbf{X})$$

$${\sigma }^i(t,\mathbf{x})$$

$$V_E(T,\mathbf{X})$$

$$V_E^i(T,\mathbf{X})$$

$$\bar{w}(t,\mathbf{x})$$

$${\bar{w}}^i(t,\mathbf{x})$$

$$\bar{V}(T,\mathbf{X})$$

$${\bar{V}}^i(T,\mathbf{X})$$

$$Q_{ij}$$

$$Q_{ijk}$$

$$C_{ijk}$$

$$Q_i=\frac{1}{M_E}\left({\ddot{I}}_E^i-Q_{ij}{I}_E^j-\frac{1}{2}{Q}_{ijk}{I}_E^{jk}\right).$$

$$C_{ij}=-{\bar{w}}^{[i,j]}({\mathbf{x}}_E)+\frac{3}{4}{v}_E^{[i}{\bar{w}}^{,j]}({\mathbf{x}}_E)+\frac{1}{4}{v}_E^{[i}{Q}^{j]},$$

$$X^i=r^i+\frac{1}{c^2}\left[\frac{1}{2}{v}_E^i({\mathbf{v}}_E\cdot \mathbf{r})+\bar{w}({\mathbf{x}}_E)r^i+r^i({\mathbf{a}}_E\cdot \mathbf{r})-\frac{1}{2}{a}_E^i{r}^2\right]+\mathcal{O}(c^{-4}),$$

$$T=t-\frac{1}{c^2}\left[A(t)+\mathbf{v}\cdot \mathbf{r}\right]+\frac{1}{c^4}[B(t)+B^i(t)r^i+B^{ij}(t)r^i{r}^j+C(t,\mathbf{x})]+\mathcal{O}(c^{-5}),$$

$$\frac{dA(t)}{dt}$$

$$\frac{dB(t)}{dt}$$

$$B^i(t)$$

$$B^{ij}(t)$$

$$C(t,\mathbf{x})$$

$$\dot{\bar{w}}\equiv \bar{w}{}_{,t}+v_E^i\bar{w}{}_{,i},$$

$$X_{\mathsf{ITRS}}^i={\Lambda }^{ij}{X}^j,$$

$$X^i={\Lambda }^{ji}{X}_{\mathsf{ITRS}}^j.$$

$${\Lambda }^{ij}=W^{ki}(T)R_3^{kp}\left(\mathrm{GAST}\right)N^{pq}(T)P^{ql}(T)B^{lj},$$

$$P^{ij}$$

$$N^{ij}$$

$$W^{ij}$$

$$v_{\mathsf{GCRS}}^i={ϵ}^{ijk}{\omega }^j{X}_{\mathsf{ITRS}}^k+{\dot{X}}_{\mathsf{ITRS}}^i,$$

Full text

Based on a fully relativistic description of the background gravitational physics, we discuss the impact of those highly-precise clocks on the realization of reference frames and time scales used in geodesy. We discuss the current definitions of basic geodetic concepts and come to the conclusion that the advances in clocks and other metrological technologies will soon require the re-definition of time scales or, at least, clarification to ensure their continuity and consistent use in practice.

The relative frequency shift between two clocks is directly related to the difference in the values of the gravity potential at the points of clock’s localization. According to general relativity the relative accuracy of clocks in $10^{-18}$ is equivalent to measuring the gravitational red shift effect between two clocks with the height difference amounting to 1 cm. This makes the clocks an indispensable tool in high-precision geodesy in addition to laser ranging and space geodetic techniques.

We show how clock measurements can provide geopotential numbers for the realization of gravity-field-related height systems and can resolve discrepancies in classically-determined height systems as well as between national height systems. Another application of clocks is the direct use of observed potential differences for the improved recovery of regional gravity field solutions. Finally, clock measurements for space-borne gravimetry are analyzed along with closely-related deficiencies of this method like an extra-ordinary knowledge of the spacecraft velocity, etc. For all these applications besides the near-future prospects, we also discuss the challenges that are related to using those novel clock data in geodesy.

††journal: Space Science Reviews

1 Introduction

The development of ultra-precise clocks opens the possibility for enhancing several geodetic applications, including an improved definition of timescales, datum connecting and the observation of Earth’s gravity field. This paper addresses the possible impact of high-precision clocks on these applications.

The definition of timescales has wide-ranging consequences for several applications. For example, there are numerous questions associated with realizing a global reference closest to the mean sea level, see e.g. a discussion in (Sánchez, 2012). Work is still under way in working groups of the International Association of Geodesy (IAG) to provide conventional definitions and procedures for the realization of a Vertical Datum Standardization.

A major objective of geodesy is the determination of physical heights, i.e. those heights related to the gravity field. The use of ultra-precise clocks to this aim has been addressed in Mai and Müller (2013), indicating that clocks with an accuracy of $10^{-18}$ allow for example datum connections at a level of 1 cm if combined with precise positioning by Global Navigation Satellite Systems (GNSS). Combination of high-precision terrestrial clock and GNSS data with e.g. space-borne gravimetric data might enhance spatial resolution and accuracy of existing global gravity field models and associated definition of vertical datums.

Finally, this paper will address the possible use of high-precision space-borne clocks for global gravimetry. This could be realized with clocks that have an ultra-high stability at short integration times (Mayrhofer and Pail, 2012). An assessment will be made of the achievable gravity field retrieval performance by an efficient error propagation tool, where the performance for ultra-precise clocks will be compared with other observation techniques including the already proved techniques of low-low Satellite-to-Satellite Traking (ll-SST) and Satellite Gravity Gradiometry (SGG).

This paper is organized as follows. After introducing the relevant background on relativistic reference systems (Section 2) and their relation to the commonly used International Terrestrial Reference System (ITRS) in Section 3, a recap is provided of the required fundamental equations of time dilation (Section 4). The impact of high-precision clocks on the definition of timescales will be addressed in Section 5. The expected benefit of using high-precision clocks for terrestrial gravimetry will be discussed in Section 6, including their use for datum connection and height comparison (Section 6.3) as well as for enhancing regional gravity field solutions (Section 6.4). This will be followed by an assessment of using high-precision clocks for space-borne gravimetry (Section 7). Finally, conclusions will be drawn in Section 8.

2 Relativistic Celestial Reference Systems

Current theory of relativistic reference systems in the solar system has been formulated by Kopejkin (1988) and Brumberg and Kopejkin (1989a, b) and further developed by Damour et al. (1991, 1992, 1993). It was adopted by the IAU 2000 as a standard for processing high-precise astronomical and geodetic observations (Soffel et al., 2003). This theory provides precise theoretical definitions of the barycentric celestial reference system (BCRS) and the geocentric celestial reference system (GCRS) as well as the relations between them (Kopeikin et al., 2011). The important distinction of the BCRS and GCRS from their Newtonian counterparts comes from general relativity which predicts the effects of relativistic contraction of (coordinate) spatial distances and time dilation caused by the relative orbital motion of GCRS with respect to BCRS and the presence of the gravitational field of Sun, Moon, and other planets. Because of these effects the spatial coordinates of the two reference systems experience periodic relative variations due to the Lorentz and Einstein contractions of the relative order of $10^{-8}$. The basic coordinate time scales for BCRS and GCRS (called TCB and TCG respectively) are also different due to the time dilation and gravitational redshift (Soffel et al., 2003).

In what follows the small letters refer to BCRS, whereas the capital letters do to GCRS. The Greek letters $\alpha,\beta,...$ denote spacetime indices taking values $0,1,2,3$ with the index [math] belonging to time coordinate. Roman letters will denote the indices taking values $1,2,3$ corresponding to the spatial coordinates only. Bold letters denote spatial vectors, ${\bm{x}}=\{x^{i}\}=\{x^{1},x^{2},x^{3}\}$, ${\bm{X}}=\{X^{i}\}=\{X^{1},X^{2},X^{3}\}$, etc. The symbol ${\cal O}(c^{-n})$ means that all residual terms of order $c^{-n}$ are neglected.

The theory of astronomical reference systems that is outlined below, is formulated in the first post-Newtonian approximation (PNA) of Einstein’s theory of gravity though it can be extended to high-order approximations if necessary. The PNA is a weak-field, slow-motion approximation with three small parameters: $\epsilon=(GM/c^{2}R)^{1/2}$ with $M$ and $R$ being a characteristic mass and size of a body; $\epsilon_{\omega}=\omega R/c$ with $\omega$ being a rotational angular velocity of the body, and $\epsilon_{v}=v/c$ where $v$ is the orbital velocity of the body. According to the virial theorem of gravitational physics (Goldstein et al., 2002) these small parameters in the solar system are ordered as follows: $\epsilon_{\omega}<\epsilon_{v}<\epsilon$, and for simplicity one is using a common denominator of the small parameters, $c^{-1}$, as a book-keeping parameter for a post-Newtonian expansion of the metric tensor, although it is not dimensionless. The first post-Newtonian approximation neglects terms of the order of $c^{-6}$ in the time-time component of the metric tensor $g_{00}$, terms of the order of $c^{-5}$ in the time-space components $g_{0i}$, and terms of the order of $c^{-4}$ in the space-space components $g_{ij}$.

In Einstein’s theory of gravity the gravitational field is described by a metric tensor denoted as $g_{\alpha\beta}=g_{\alpha\beta}(x^{\mu})$ in BCRS and $G_{\alpha\beta}=G_{\alpha\beta}(X^{\mu})$ in CGRS, that also characterizes the geometry of spacetime. The metric tensor is found by solving the Einstein field equations which preserve their form irrespectively of the choice of coordinates. Thus, they can determine the components of the metric tensor only up to four degrees of freedom. This gauge freedom of the metric tensor components corresponds to the free choice of the coordinate system. We will use the harmonic gauge condition (Fock, 1964) imposed on the metric tensor in every coordinate system. The harmonic coordinates are convenient from a mathematical point of view and are often used in other branches of gravitational physics.

Post-Newtonian form of the metric tensor in harmonic coordinates reads (Brumberg and Kopejkin, 1989a, b; Damour et al., 1991)

[TABLE]

where the gravitational potentials are functions of time and spatial coordinates: $w\equiv w(t,{\bm{x}})$, $w^{i}\equiv w^{i}(t,{\bm{x}})$ in BCRS, and $V\equiv V(T,{\bm{X}})$, $V^{i}\equiv V^{i}(T,{\bm{X}})$ in GCRS, respectively. Notice that we have used a letter $V$ to denote the gravitational potentials of the geocentric metric tensor $G_{\alpha\beta}$. The paper by Damour et al. (1991) used $W$ instead of $V$. However, the letter $W$ is ubiquitously used in geodesy to denote the potential of the gravity force in rotating geocentric coordinates system and is represented as an algebraic sum of the gravitational potential $V$ and the centrifugal potential, $Z=1/2{\bm{v}}^{2}$ where ${\bm{v}}=\{v^{i}\}$ is velocity of observer with respect to GCRS (see, for example, Eq. (51) below.)

The gravitational potentials are expressed in terms of the integrals taken over three-dimensional volumes occupied by matter of the bodies comprising the solar system. We single out the gravitational potentials associated with Earth and those produced by the external bodies (Moon, Sun, planets, etc.). In the post-Newtonian approximation this decomposition reads

[TABLE]

and

[TABLE]

where the potentials with sub-index $E$ belong to Earth, and those with the bar to the external bodies.

Explicit expressions of the geopotentials are given in the form of a particular solution of the Poisson equations that are volume integrals,

[TABLE]

where $\sigma(t,{\bm{x}})$, $\Sigma(T,{\bm{X}})$ and $\sigma^{i}(t,{\bm{x}})$, $\Sigma^{i}(T,{\bm{X}})$ are the post-Newtonian mass and mass-current densities of Earth’s matter in the BCRS and GCRS respectively. The densities are directly related to the model of the energy-momentum tensor of matter distribution inside Earth. In most of the practical applications considered below it is sufficient to assume that the mass densities are approximately equal

[TABLE]

where $\rho(t,{\bm{x}})$ is the baryon rest mass density. The current densities

[TABLE]

where $v^{i}_{E}=v^{i}_{E}(t)$ and ${\bm{\omega}}=\{\omega^{i}(t)\}$ are the orbital velocity of Earth and the instantaneous angular velocity of Earth’s rotation, respectively, and the cross between two vectors denote the standard Euclidean cross product.

Earth’s gravitational potentials $V_{E}$ and $V_{E}^{i}$ admit multipolar expansion in the exterior space (outside the Earth),

[TABLE]

where $r=|{\bm{X}}|$ is the radial distance in GCRS, $M_{E}$ is the total (relativistic) mass of Earth, ${\bm{S}}_{E}=\{S^{i}_{E}\}$ is the angular momentum (spin) of Earth, $I_{E}^{i}$ and $I^{ij}_{E}$ are dipole and quadrupole moments of Earth’s, and the overdot denotes a time derivative with respect to $T$=TCG. Usuall assumption is that the Earth’s center of mass is located at the origin of GCRS which makes $I^{i}_{E}=0$, and all terms depending on the dipole moment $I^{i}_{E}$ vanish from equations. We prefer to leave the dipole term explicitly in the equations because it is used as a vector parameter in tracking down the motion of the geocenter with respect to the origin of the International Terrestrial Reference System (ITRS) - for more detail see (Wu et al., 2012; Kuzin et al., 2010; Rebischung et al., 2013).

Explicit expressions of the external gravitational potentials in BCRS are given in terms of the integrals performed over the volumes of the external bodies,

[TABLE]

where $\sigma(t,{\bm{x}}^{\prime}{})$ and $\sigma^{i}(t,{\bm{x}}^{\prime}{})$ describe the distribution of mass and mass-current densities inside the volume of the external body $A$.

External gravitational potentials in GCRS are found as general solutions of the Laplace homogeneous equation that are given in terms of polynomials

[TABLE]

where $Q_{i}=Q_{i}(T)$, $Q_{ij}=Q_{ij}(T)$ and $Q_{ijk}=Q_{ijk}(T)$ are the dipole, quadrupole and octupole moments of the tidal gravitoelectric field, $C_{ij}=C_{[ij]}(T)$, $C_{ijk}=C_{[ij]k}(T)$ are the quadrupole and octupole moments of the tidal gravitomagnetic field where the square parentheses around indices indicate the antisymmetry. The external multipole moments can be expressed in terms of the partial derivatives from the external gravitational potentials which are found by making use of the asymptotic matching of the metric tensor in BCRS and GCRS (Kopeikin et al., 2011)

[TABLE]

where the dot over function denotes a total derivative with respect to time, and $v^{i}_{E}=dx^{i}_{E}/dt$ is the velocity of the geocenter with respect to BCRS. Notice that the derivatives from the external potentials are calculated on the world line of the geocenter, ${\bm{x}}_{E}={\bm{x}}_{E}(t)$, so that, for example, $\bar{w}_{,ij}({\bm{x}}_{E})=\bar{w}_{,ij}({\bm{x}})|_{{\bm{x}}={\bm{x}}_{E}}$, etc.

Vector quantity $Q_{i}$ is a small acceleration of the geocenter’s world line with respect to a geodesic world line that is caused by the coupling of the dipole, $I^{i}_{E}$, and quadrupole, $I^{ij}_{E}$, moments of Earth’s gravitational field to the tidal gravitational field of the external bodies (Kopeikin et al., 2011),

[TABLE]

The quadrupole moment, $C_{ik}$, is a post-Newtonian matrix of the dynamic rotation of the GCRS spatial axes with respect to the BCRS ones which appears in the metric tensor because of the IAU 2000 resolution demanding to keep the spatial axes of the BCRS and GCRS aligned to make both coordinate system kinematically-nonrotating (Brumberg and Kopejkin, 1989a; Damour et al., 1991). This demand, however, makes GCRS dynamically rotating with the post-Newtonian angular velocity corresponding to the infinitesimal matrix of rotation $C_{ij}$. With a sufficient accuracy the mathematical expression for the matrix of the post-Newtonian dynamic rotation reads (Kopeikin et al., 2011)

[TABLE]

where the first term in the right side is the gravitomagnetic (Schiff or Lense-Thirring) precession, the second term is the geodetic (de Sitter) precession, and the third term is the Thomas precession caused by the non-geodesic motion of Earth’s geocenter with acceleration $Q_{i}$.

The asymptotic matching technique allows us to derive the transformation law between time and spatial coordinates of BCRS and GCRS (Kopeikin et al., 2011). The post-Newtonian transformation between the spatial coordinates of GCRS and BCRS reads (Soffel et al., 2003; Petit and Luzum, 2010)

[TABLE]

where ${\bm{r}}\equiv{\bm{x}}-{\bm{x}}_{E}(t)$, and ${\bm{a}}_{E}$ is the coordinate acceleration of the Earth’s geocenter, ${\bm{a}}_{E}=d{\bm{v}}_{E}/dt=d^{2}{\bm{x}}_{E}/dt^{2}$.

Post-Newtonian transformation of time coordinates, TCG=$T$ and TCB=$t$, is more complicated. According to the IAU 2000 Resolutions it reads (Soffel et al., 2003; Petit and Luzum, 2010)

[TABLE]

with

[TABLE]

Here, the dot stands for the total time derivative with respect to time $t$, i.e.,

[TABLE]

and the round parentheses denote the symmetrization with respect to two spatial indices.

3 International Terrestrial Reference System

The science of geodesy primarily deals with observations and measurements conducted by terrestrial observers located on Earth. Therefore, it is practically convenient to work in the rotating geocentric coordinate system - the international terrestrial reference system (ITRS) which time coordinate is TCG coordinate time $T$ (the same as the time coordinate in GCRS) and the spatial coordinates are denoted $X_{\scriptscriptstyle\sf{ITRS}}^{i}=\left(X_{\scriptscriptstyle\sf{ITRS}}^{1},X_{\scriptscriptstyle\sf{ITRS}}^{2},X_{\scriptscriptstyle\sf{ITRS}}^{3}\right)$. Transformation from GCRS to ITRS is given by the IERS Conventions (Petit and Luzum, 2010; Kopeikin et al., 2011)

[TABLE]

where $\Lambda^{ij}\equiv\Lambda^{ij}(T)$ is the orthogonal matrix of rotation depending on time $T$. Due to the property of the orthogonal matrices the inverse matrix $\Lambda^{-1}$ of the transformation coincides with the transpose matrix $\left(\Lambda^{-1}\right)^{ij}=\Lambda^{ji}$ so that the inverse transformation between GCRS and ITRS is

[TABLE]

According to the IERS theory of the Earth rotation (Petit and Luzum, 2010) the matrix $\Lambda^{ij}$ can be represented in two equivalent forms corresponding to the, so called, CIO-based transformation and equinox-based transformation (see (Petit and Luzum, 2010, Chapter 5.9) and (Kopeikin et al., 2011, Equation 9.75)). For analytic consideration the equinox-based transformation is more convenient and therefore used in the following. The matrix of the equinox-based transformation is represented as a product of four matrices

[TABLE]

where the time $T$=TCG, $W^{ki}$ is the matrix of the polar wobble, $R_{3}^{kp}\left(\rm{GAST}\right)$ is the matrix of the diurnal rotation depending on the Greenwich Astronomical Sidereal Time (GAST) that is a function of time $T$=TCG, $N^{pq}$ is the matrix of nutation, $P^{ql}$ is the matrix of precession, and the constant matrix $B^{lj}$ describes the, so-called, frame bias.

These matrices have the following structure (Petit and Luzum, 2010)

[TABLE]

where $\varepsilon_{0}$ and $\varepsilon=\varepsilon(T)$ are correspondingly a constant and instantaneous obliquity of the celestial equator to ecliptic, $\chi_{A}=\chi_{A}(T)$, $\omega_{A}=\omega_{A}(T)$, $\psi_{A}=\psi_{A}(T)$ are secular variations in precession, $\Delta\epsilon=\Delta\epsilon(T)$ and $\Delta\psi=\Delta\psi(T)$ are periodic variations in nutation, $x_{\rm p}=x_{\rm p}(T)$ and $y_{\rm p}=y_{\rm p}(T)$ describe the polar wobble, and $s^{\prime}=s^{\prime}{}(T)$ is the small secular variation describing the shift between the ITRS origin of longitude and the terrestrial intermediate origin (TIO) (Kopeikin et al., 2011; Petit and Luzum, 2010).

The components of the rotational matrix are used to calculate the GCRS velocity of motion of a terrestrial observer (clocks) located on Earth’s surface,

[TABLE]

where the overdot denotes a time derivative, $\epsilon^{ijk}$ is a fully anti-symmetric symbol of Levi-Civita with $\epsilon^{123}=+1$, $X_{\scriptscriptstyle\sf{ITRS}}^{i}$ is the ITRS coordinate of the observer (clock), $\dot{X}_{\scriptscriptstyle\sf{ITRS}}^{i}$ is the residual velocity of the observer (clock) with respect to ITRS due to various geophysical reasons or simply because the observer is moving in a car or aircraft, and $\omega^{i}=\omega^{i}(T)$ is the instantaneous angular velocity of the rotation of the ITRS with respect to GCRS,

[TABLE]

In what follows it will be convenient to split the ITRS position of the observer (clock) in two components

[TABLE]

where $n^{i}(h)$ is a unit vector along the direction of the plumb line passing through the point of observation, $X_{\scriptscriptstyle\sf{ITRS}}^{i}$ located on Earth surface, $X_{\scriptscriptstyle\sf{geoid}}^{i}$ is the point on geoid connected by the plumb line to the position of the observer on surface, $dh$ is the element of length along the plumb line, and $H$ is the orthometric height at the observer’s position (Torge and Müller, 2012).

Precise analytic expression for the angular velocity, $\omega^{i}$, is too complicated and we develop it only up to a linear approximation with respect to the variations of the parameters entering (35)–(37). For this purpose we use the following approximations

[TABLE]

where the very last terms in (37)-(41) are functions of time, $T$=TCG, which are changing due to the systematic (both periodic and secular) variations in the spatial orientation of the Earth’s rotational axis, wobble, and tidal friction (Kopeikin et al., 2011; Petit and Luzum, 2010).

Then, in the linear approximation the components of the angular velocity of Earth’s rotation are

[TABLE]

where $\Omega$ is a fixed value of the angular velocity of the Earth rotation adopted by IAU (Petit and Luzum, 2010).

4 Time dilation fundamental equation

We consider a network of atomic clocks located on Earth’s surface at different geographic positions. Each clock moves with respect to the GCRS along world line $X^{i}\equiv X^{i}(T)$. According to general relativity each clock measures its own proper time $\tau$ which is defined by the equation $-c^{2}d\tau^{2}=ds^{2}$ where the interval $ds$ must be calculated along the world line of the clock. In terms of the GCRS metric tensor, the interval of the proper time reads

[TABLE]

where $v^{i}$ is velocity of clock with respect to GCRS, and $T$=TCG. In case of clocks in space $v^{i}$ is the orbital velocity of the spacecraft carrying the clok. If the clock is located on Earth’s surface, the GCRS velocity of the clock is given by (38).

After replacing the GCRS metric (2) in (49) and extracting the root square, we get the fundamental time delay equation

[TABLE]

where time-dependent function $W=W(T)$ is given by (Kopejkin, 1991; Kopeikin et al., 2016)

[TABLE]

Function $W$ is the post-Newtonian gravity potential taken at the point of localization of the clock. Notice that it includes both the effects of the gravitational field of the Earth and the external bodies (Sun, Moon, planets) in the form of the tidal terms as shown in equation (4) of the present paper. The clock’s proper time can be calculated by integrating (50) along the world line of the clock,

[TABLE]

where $T_{0}$ is the initial epoch of the integration.

5 Implications for the definition of timescales

The 13th General Conference on Weights and Measures (CGPM) decided in its Resolution 1 in 1967 111The CGPM reports are available at http://www.bipm.org/fr/worldwide-metrology/cgpm/resolutions.html that the second is the duration of 9,192,631,770 periods of the radiation corresponding to the transition between two hyperfine levels of the ground state of the caesium 133 atom. This definition of the time unit called for the adoption of a time scale built by cumulating atomic seconds and International Atomic Time TAI was defined in 1970 by the International Committee for Weights and Measures as the time reference established by the BIH on the basis of the readings of atomic clocks operating in various establishments in accordance with the definition of the second, a definition recognized by the $14^{th}$ CGPM in its Resolution 1 in 1971. In 1980 the definition of TAI was completed by the Consultative Committee for the Definition of the Second, adding TAI is a coordinate time scale defined in a geocentric reference frame with the SI second as realized on the rotating geoid as the scale unit. This definition explicitly refers to TAI as a coordinate time, hence needing a relativistic approach. In 1988 the responsibility of establishing TAI was transferred to the International Bureau for Weights and Measures (BIPM) in Sèvres (France).

TAI is not disseminated directly and Coordinated Universal Time UTC, which was designed to approximate UT1 (a timescale derived from the rotation of the Earth), was chosen as the practical world time reference. Since 1972, UTC differs from TAI by an integral number of seconds, changed when necessary by insertion of a leap second, as predicted and announced by the International Earth Rotation and Reference System Service (IERS). It is not the purpose of this paper to discuss this issue; the point is just to remind that changes in the definition of timescales may have wide-ranging consequences irrespective of the practical implications.

As mentioned in previous sections, TCB and TCG are the time coordinates of the BCRS and GCRS, respectively. The 1991 Recommendation 3 of Resolution A4 of the International Astronomical Union (IAU) 222The IAU Resolutions are available at http://www.iau.org/administration/resolutions/general_assemblies/ defined the scale unit of TCB and TCG to be consistent with the SI second. This means that if readings of the proper time of an observer, expressed in SI seconds, are recomputed into TCB or TCG using the formulas from the IAU Resolutions, without any additional scaling, one gets corresponding values of TCB or TCG in the intended units. The Recommendation also defines the origin of TCB and TCG by the following relation to TAI: TCB (resp. TCG) = TAI + 32.184 s on 1977 January 1st, 0 h TAI, at the geocenter.

In the following, we assume that the temporal variations in the physical quantities (e.g. potential, equipotential surface) due to tides are taken into account through appropriate reductions so that the quantities are considered as quasi-static, i.e. constant or with a slowly varying secular change. For a clock at rest on the Earth’s surface, the relation between proper time and coordinate time is given by Eq. (50) that is

[TABLE]

where for practical purposes, it is sufficient to retain in $W$ only the first two terms of the gravity potential (rotational plus gravitational) given by Eq. (51). Any time differing from TCG by a constant rate may also be chosen as a coordinate time in the geocentric system, and this is the case of Terrestrial Time TT which differs from TCG by a constant rate:

[TABLE]

In the original definition (IAU’1991 Recommendation 4 of Resolution A4), $L_{G}$ is such as the unit of measurement of TT is chosen so that it agrees with the SI second on the geoid, i.e. $L_{G}=W/c^{2}$, where $W$ is the latest estimate of the gravity potential on the geoid. By IAU Resolution B9 (2000) TT is explicitly re-defined with respect to TCG by turning $L_{G}$ in Eq. (54) into a defining constant where $L_{G}~{}=~{}6.969290134\times 10^{-10}$. This constant has been computed as $L_{G}~{}=~{}W_{0}/c^{2}$ taking for $W_{0}$ the value $62636856$ m2/s2 recommended by Special Commission 3 of the International Association of Geodesy (IAG) in the year 1999, before the new definition. This value will be noted below $W_{0}$(2000).

TAI is a realization of coordinate time TT, to within the constant offset of 32.184 s, and is subject to uncertainties in the realization. Timescales derived from TAI, i.e. UTC and other realizations steered to UTC such as GNSS reference times, are also realizations of TT with the same remarks. Because TAI is defined with respect to the rotating geoid, the transformation from proper time to TAI requires, in principle, the actual value, noted $W_{0}$, of the gravity potential on the geoid. Relation (53) transforms to

[TABLE]

where $W_{0}=W({\bm{X}}_{\scriptscriptstyle{\sf geoid}})$ and ${\bm{X}}_{\scriptscriptstyle{\sf geoid}}=\{X^{i}_{\scriptscriptstyle{\sf geoid}}\}$ is the ITRS coordinate on the geoid. In practice, the potential $W$ is expanded in the Taylor series around its value on the geoid (Torge and Müller, 2012)

[TABLE]

so that Eq. (55) is generally applied as

[TABLE]

where $H$ is the orthometric height of the clock and $\bar{g}$ the average value of the acceleration of gravity between the geoid and the clock.

Because a change in height of one meter causes a change in rate of about $1\times 10^{-16}$, Eq. (57) can be somewhat loosely applied for present-day Caesium primary frequency standards which have accuracy of order $1\times 10^{-16}$, assuming that the geoid, the reference level of the height system and the height are all correctly realized or measured with an uncertainty well below one metre. This is no more the case when considering a clock accuracy of the order $1\times 10^{-17}$ and below for which an uncertainty of order 1 cm is needed. For better accuracy, with validity about $5\times 10^{-19}$, Eq. (53) must include tidal terms $V_{\sf tide}\equiv\bar{V}$ given in (17) to obtain

[TABLE]

so that

[TABLE]

The tidal potential has been considered in full detail, for example, in (Agnew, 2007). In the quadrupole gravitational field approximations it reads $V_{\rm tide}\equiv V_{2}$, i.e., we take into account in (17) merely the term depending on the external quadrupole moment $Q_{ij}$. This expression is well known in scientific literature (see, for example, (Agnew, 2007; Torge and Müller, 2012; Simon et al., 2013)), and can be presented in the form

[TABLE]

where

[TABLE]

where

[TABLE]

is the Doodson’s constant for the body A, $\bar{r}_{\scriptscriptstyle{\sf EA}}$ is the reference (mean) value of the distance $r_{\scriptscriptstyle{\sf EA}}$ during one revolution of the body A around Earth, $\delta_{A}$, $h_{A}$ are respectively the declination and the hour angle of the external body A in the ITRS, and the summation goes over all the external bodies of the solar system but Earth.

Note that all developments in this section are valid to an accuracy level of order $1\times 10^{-18}$ and should be re-examined to reach an accuracy level of order $10^{-19}$ ($\equiv$ mm in height). Then, e.g., the full post-Newtonian gravity potential given in Eq. (51) as well as the geoid in its relativistic definition should be considered, cf. (Kopeikin et al., 2011).

Consider the transformation to TCG (the case of TAI will be considered further below) with Eq. (58), one needs to determine the gravity potential at the location of the clocks. This used to be a difficult task at the $10^{-17}$ level worldwide (Petit and Wolf, 1997), however recent global satellite geopotential models allow this determination at the level of a few parts in $10^{18}$ (Denker, 2013). On the other hand the centrifugal potential as well as other time-varying tidal effects due to external masses (direct effect and indirect effect from Earth and ocean tides) can reach a few $10^{-17}$ but can be computed at the $1\times 10^{-18}$ accuracy level (Voigt et al., 2016). In addition, other non-tidal effects due to mass redistributions in the Earth, the hydrosphere or the atmosphere are to be considered at the $10^{-18}$ accuracy level. It is not the purpose of this paper to study these effects but, contrary to tidal effects, it is not straightforward to model them and to define a conventional reference for these effects. Such $10^{-18}$ accuracy in the computation of the relativistic frequency shift will be necessary for the realization of timescales when new clocks with such frequency accuracy will provide the definition of the second.

Considering that the best Caesium primary frequency standards have an accuracy of order $1-2\times 10^{-16}$, the present requirements are less stringent. In past years, it was claimed an uncertainty for the relativistic frequency shift of $2-3\times 10^{-17}$ for the location of the NIST in Boulder (Colorado, USA) (Pavlis and Weiss, 2003) and of $1\times 10^{-17}$ for the INRIM in Torino (Italy) (Calonico et al., 2007). A review of the frequency standards contributing to TAI in 2015 333BIPM report on Time activities available at http://www.bipm.org/en/bipm/tai/annual-report.html shows that a good fraction adopts a conservative approach in evaluating the relativistic frequency shift e.g. $1\times 10^{-16}$ is used for the LNE-SYRTE in Paris (France) and for the NIM in Beijing (China), $0.5\times 10^{-16}$ for the NPL in Teddington (UK) and for the VNIIFTRI in Mendeleevo (Russia), while the NIST adopts $0.3\times 10^{-16}$ as mentioned above. Some metrological centers, however, like the INRIM or the PTB in Braunschweig (Germany) report an uncertainty at or slightly below $0.1\times 10^{-16}$.

The geoid appearing in the definition of TAI is classically defined as the level surface of the gravity potential closest to the topographic mean sea level. Therefore the value of the potential on the geoid $W_{0}$ depends on the global ocean level which changes with time due to geophysical reasons. Several authors have considered the time variation of $W_{0}$, see e.g. (Burša et al., 2007; Dayoub et al., 2012), but there is some uncertainty in what is accounted for in such a linear model. A recent estimate (Dayoub et al., 2012) over 1993-2009 is $dW_{0}/dt=-2.7\times 10^{-2}~{}\rm m^{2}s^{-2}yr^{-1}$, mostly driven by the sea level change of +2.9 mm/yr. The rate of change of the global ocean level could vary during the next decades, nevertheless, to state an order of magnitude, considering a systematic variation in the sea level of order 3 mm/yr, different definitions of a reference surface for the gravity potential could yield differences in frequency of order $3\times 10^{-18}$ in a decade. In addition, there are numerous questions associated with realizing a global reference closest to the mean sea level, see e.g. a discussion in (Sánchez, 2012), so that work is still under way in IAG working groups to provide conventional definitions and procedures of realization for a Vertical Datum Standardization.

Therefore, there is evidence of inconsistency in the present set of concepts and/or definitions as summarized hereafter. On the geodetic side, there is inconsistency between the three propositions:

• •

The geoid is the equipotential surface closest to the mean sea level;

• •

The value of the gravity potential on the geoid is $W_{0}$;

• •

$W_{0}$ is a defining (fixed) constant.

On the time metrology side, there is inconsistency between the three propositions:

• •

TT is a coordinate time defined by IAU Recommendation B9 (2000);

• •

TAI is a coordinate time with the SI second as realized on the rotating geoid as the scale unit;

• •

TAI is a realization of TT.

The inconsistency is at a level of $10^{-17}$ (equivalent to 9 cm in height) when considering that such fundamental concepts should be valid for decades.

The choice of the value $W_{0}$(2000) provides a formal definition of a surface where clocks run at the same rate as TT, which has been named the chronometric geoid (Wolf and Petit, 1995). As the classical geoid remains linked to the mean sea level, these two surfaces differ and the difference evolves with time. One possible solution on the geodetic side is to update the $W_{0}$ value as necessary to reflect the actual mean sea level as well as improvements in the determination of $W_{0}$, so that $W_{0}$ is no more used as a constant but is varying with time. On the time metrology side, one solution is to explicitly define TAI as a realization of TT so that it is no more linked to the time-varying geoid. Indeed, if the adopted value of the gravity potential on the geoid $W_{0}$ is updated, this has no implication on the $L_{G}$ value which remains a fixed conventional value relating TCG and TT. The $W_{0}$ value and the definition of the geoid are the responsibility of the IAG which addressed the issue in its Resolution 1 (2015) for the definition and realization of an International Height Reference System (IHRS). As mentioned in (Ihde et al., 2015), For a global height reference system, any value $W_{0}$ within a range of a few decimeters can be defined as conventional without affecting the task of defining and realizing a global height reference system. Nevertheless, IAG Resolution 1 (2015) recommends both a new value of the gravity potential at the geoid $W_{0}=62636853.4~{}\rm m^{2}s^{-2}$, here noted $W_{0}$(2015), and to use this value a reference for the IHRS. Should such IHRS-compatible geopotential values become available for vertical reference points, the difference between $W_{0}$(2015) and $W_{0}$(2000), equivalent to $2.9\times 10^{-17}$ in rate shift, has to be noted when realizing TT from an atomic clock. The reformulation of the TAI definition as well as its procedure of realization is under the responsibility of the Consultative Committee for Time and Frequency (CCTF) which has formed a task group to propose, at its 2017 meeting, a new definition of TAI to be adopted by the CGPM in 2018. Changing the wording in the definition of TAI will be one task but it may be the easy one. Applying Eq. (59) will still be a difficult task on a global scale at the $10^{-18}$ level. E.g. while height networks can be directly linked when on the same continent, larger uncertainties persist in the global determination of the geoid and of vertical references from effects affecting the sea surface topography (Fu, 2010).

Placing ultra-accurate clocks on Earth is inconvenient, because of the large number of small geophysical effects that have to be taken into account to compute the relativistic frequency shift. However it is also an opportunity because their accuracy gives access to study these effects. As has long been realized and dreamed of, it might be an optimal solution to have some ultra-accurate clocks in space, to be used as a reference to generate TAI and to which Earth-based clocks could be compared. Comparisons of ultra-accurate clocks could therefore help in the future to establish a worldwide vertical datum.

6 Height systems and terrestrial gravimetry with clocks

Here, we discuss a novel method of the direct use of clock measurements to derive physical heights and to resolve discrepancies in classical height systems. We address time-variable effects that have to be considered when using clock data in practice. These effects include both temporal variations of the gravity field and instabilities in the rotation of the Earth. A particular application is the direct use of observed potential differences for the improved recovery of the regional gravity field in the Massif Central region in the middle of southern France.

6.1 Height systems and clock measurements

To the major objectives of geodesy belong the determination of physical heights (i.e. those heights related to the gravity field) and the determination of the corresponding height reference surface (i.e. the geoid or quasi-geoid) (Torge and Müller, 2012), see Figure 1. For example, the orthometric height is defined by

[TABLE]

where $C_{P}$ is the so-called geopotential number which is the difference of the gravity potential between the geoid and the surface point $P$, and $\bar{g}$ is mean gravity along the (curved) plumbline. A similar relation holds for normal height $H_{N}$ (where the geopotential number $C_{P}$ has to be divided by a mean normal gravity value $\bar{\gamma}$) and quasi-geoid. Geopotential numbers $C_{P}$ are classically obtained via geometrical levelling and terrestrial gravimetry. The latter is required to correct the raw levelling results to consider the non-parallelism of the equipotential surfaces of which the geoid is just that one surface corresponding to the mean ocean surface at rest. This classical approach has some drawbacks such as increasing errors with the length of the levelling loops or systematic errors by combining data from different measurement periods. This method is very time consuming, especially if large areas have to be covered. For detecting vertical deformations, repeated levelling measurements are required which further increases the mentioned complexity.

When establishing a physical height system, normally one refers to a datum point like a selected tide gauge, e.g. NAP (Normaal Amsterdams Peil), and then only differences of potential numbers $\Delta C_{P}$ are used further on. The basic network of physical heights for a certain area, like a country, is then computed at the level of geopotential numbers before they are converted into physical heights by dividing the $C_{P}$ values by $\bar{g}$ or $\bar{\gamma}$.

If the geoid (or quasi-geoid) shall be determined point-wisely, one can use the physical heights together with ellipsoidal GNSS (Global Navigation Satellite System) heights

[TABLE]

where $N$ is the geoid height above the ellipsoid, $h$ is the ellipsoidal height of the surface point $P$. Similar relations hold for the quasi-geoid height: $\zeta=h-H_{N}$. Alternatively, one could apply the Bruns formula

[TABLE]

where the geoid height $N$ is derived from the disturbing potential $T_{0}$ at the geoid and the normal gravity value $\gamma_{0}$ taken at the ellipsoid. A similar expression holds again for the quasi-geoid height $\zeta=T_{P}/\gamma_{Q}$, now the disturbing potential is needed at the surface point $P$ and the normal gravity value at the telluroid point $Q$444The telluroid is an approximation to the Earth surface that is consistent to the theory of Molodensky (Heiskanen and Moritz, 1967).. The disturbing potential $T_{P}$ is defined as difference between the gravity potential of the Earth $W_{P}$ and the normal potential $U_{P}$ which is computed from the values of a mean rotational ellipsoid approximating the Earth: $T_{P}=W_{P}-U_{P}$, cf. (Torge and Müller, 2012). The disturbing potential can be computed by solving a boundary-value problem where various variants exist. Often the long-wavelength contributions of the gravity field are taken from satellite solutions that are combined with terrestrial gravity and topographic data for the shorter wavelengths, cf. (Denker, 2013).

As briefly re-called from geodesy textbooks, the gravity potential and mostly differences of the gravity potential are the central quantities to derive physical heights and the geoid. In Section 4 , the dependency of proper time to the gravity potential $W$ is given. Now we consider two clocks with proper times $\tau_{1}$ and $\tau_{2}$ (see also Figure 2) which are directly related to their proper frequencies $f_{1}$ and $f_{2}$ via

[TABLE]

The corresponding proper frequency difference $\Delta f$ between the two clocks placed at two different locations is now derived by introducing Eq. (50) in Eq. (68)

[TABLE]

It is directly related to the difference of the gravity potential $\Delta W_{21}$ between these locations which again corresponds to the difference of the geopotential numbers $\Delta C_{21}=-\Delta W_{21}$. Thus, in future, those clock-based gravity potential values can straightforwardly be used in physical geodesy. For example, the difference between two orthometric heights can be obtained as (see also (Müller, 2016))

[TABLE]

with the difference of the mean gravity values along the plumbline at both locations $\Delta\bar{g}_{12}=\bar{g}_{1}-\bar{g}_{2}$ and the difference of the geopotential numbers $\Delta C_{21}=C_{2}-C_{1}=W_{1}-W_{2}$, which can now be replaced by clock-based frequency differences $\Delta f_{21}$ according to Eq. (69):

[TABLE]

6.2 Time-variable effects on clock measurements

In the following, we want to discuss some properties that have to be considered when using clock measurements and we will give a few examples for potential applications where clock data might especially be beneficial.

At the accuracy level of mm to cm, it is not sufficient anymore to consider gravity potential values as constant in time. They are affected by deformations (up to 40 cm due to tidal effects) and mass variations. In the following, we always indicate the combined effect for both contributions (deformation plus mass effect) as sensed by the clocks. And to make the effect more illustrative, the gravity potential value is converted into a height value by dividing it by some mean gravity value of about 9.81 m/s2. Even if only potential differences are considered, there are remarkable effects of up to 8 cm between PTB in Braunschweig and NPL in London caused by relative perturbations of the solid Earth tides ((Voigt et al., 2015, 2016)). Those variations with changing amplitudes have periods of 12 hours and longer. In addition, many smaller effects have to be considered: periodic effects due to ocean tides (a few cm), non-tidal oceanic and atmospheric effects (1 cm) or episodic effects, e.g., caused by storm surges or droughts (1 cm). Further effects comprise pole tides, variations due to land hydrology, tectonic and GIA-induced processes, etc., see (Voigt et al., 2016). Another class of time-variable effects on the gravity potential are variations in Earth orientation.

The Earth orientation effects affect the rate of clocks located on the Earth surface as they are produced by the time-dependent changes in the velocity-dependent potential of the centrifugal force (Fateev et al., 2015). The centrifugal force potential consists of three terms as a direct consequence of the velocity decomposition equation (38)

[TABLE]

The first term on the right hand side of (72) is the centrifugal potential caused by the rotation of ITRS in space. We shall denote it as

[TABLE]

and decompose it in several components corresponding to secular changes in the rotational velocity, precesssion, nutation, and the polar wobble. It is convenient to write down the decomposition in terms of geocentric spherical coordinates associated with ITRS,

[TABLE]

The centrifugal potential can be represented as a superposition

[TABLE]

where

[TABLE]

The second term in (72) describes the coupling of the rotational velocity of ITRS and the residual velocity $\dot{\bm{X}}{\scriptscriptstyle{\sf ITRS}}$ of clocks due to their motion with respect to ITRS. Taking the time derivative from the ITRS coordinates (74) of clocks and calculating, we get

[TABLE]

We can observe that the vertical speed of clocks and the motion of the clocks along meridian does not contribute to $Z_{\rm c}$, only the longitudinal component of the velocity matters. The term $Z_{\rm c}$ exceeds $10^{-18}$ if the tangential residual velocity of clocks on the Earth’s equator exceeds 0.2 mm/s. We shall assume that clocks are at rest which allows us to neglect the contribution of $Z_{\rm c}$ to the velocity-dependent potential along with the third term in (72). In important cases of clock transportations on cars or aircrafts the term $Z_{c}$ must be retained. All rotational effects can well be modelled. In future, some of those time-variable effects might even be measured in clock networks.

6.3 Clocks for datum connection and uncertainties in height systems

What is the major benefit of highly-precise clock measurements? Clocks can connect distant areas, i.e. they deliver physical height differences for the observed points without being affected by levelling errors or by some smoothing effect when combined global gravity field models, e.g., from GOCE are used. No (good) ground gravity data from terrestrial gravimetry is needed to fill-in regional to local gravity variations. This might especially be interesting in (relatively) unsurveyed countries where no gravity data are available at all, where terrestrial data are difficult to obtain (e.g. in rain forests) or in areas with rough environment (e.g. in mountain areas).

Clock measurements can be used as independent observations to resolve discrepancies in (classical) realisations of height systems and geoid solutions, see example in the Alps (Denker et al., 2015). Even in Germany between north and south over only 700 km, there is a 4 cm discrepancy when determining the quasi-geoid difference by different methods, such as GNSS/levelling versus gravimetric solutions. A similar comparison between NPL, London and PTB, Braunschweig showed differences between the classical methods of up to 15 cm (Denker, 2016). These differences are most probably caused by geometric levelling that accumulates measurement errors with increasing distance.

Nevertheless in Europe, the connection of national height systems (i.e. of the specific vertical datums of each country) might be possible by applying standard geodetic methods - with the above mentioned uncertainties. But in regions like South America, large differences with large error bars are present when comparing various tide gauges, especially in the southern part (see Figure 3 taken from (Sánchez, 2015)) where differences in the decimeter level are present. Here, immediate improvements could be achieved, even with recent clock technology at the few cm level of accuracy, assuming that these clocks could be placed at those tide gauges and be connected by dedicated fiber links. Besides just taking raw difference measurements, one could carry out a new adjustment for the height differences where clock, GNSS and altimetry data are used as joined input.

Another application could be running clocks on ships to independently control altimetric measurements. The ocean surface roughly approximates an equipotential surface in the Earth gravity field (then clock rates would not vary), but deviate from it because of ocean currents like the gulf stream. These deviations reach up to 1 to 2 m and are observed by satellite altimetry together with the marine geoid, e.g., derived from GOCE data (Pail et al., 2010). One idea now is to use combined clock and GNSS measurements along selected profiles that can be compared to the results of altimetric measurements. Once again, possible systematics might by revealed by the clock measurements. Running a clock on a moving platform without interruption, however, is another challenge that is not discussed here. And, it might even be more challenging when the moving clock has to be (frequently) compared to a remote clocks with the same accuracy.

6.4 Clocks for high resolution geopotential recovery

To illustrate what could be the benefit of clocks for high-resolution geopotential recovery, beyond the resolution of the satellites and as complement to near-surface gravity data, we perform simple simulations where we compare the quality of the geopotential reconstruction with and without adding clock data. We consider the Massif Central area, marked by smooth, moderate altitude mountains and volcanic plateaus, see Figure 4 (left), which leads to variations of the gravitational field over a range of spatial scales, as illustrated in Figure 4 (middle and right).

For the tests, we proceed as follows. We sample synthetic gravity and potential data from a spherical harmonics geopotential model, and we build a geopotential control grid. We estimate the gravity potential by least-squares collocation from the synthetic data. We compare the recovered potential on the control grid, to the reference one to assess the improvements that can be obtained using clock data. In this process, different parameters of the problem can be varied, as the spatial density of clock data for instance, studied hereafter.

In more details, the main steps of this methodology are thus:

High spatial resolution 5-km step grids of gravity anomalies $\overline{\delta{g}}$ and the disturbing potential $\overline{T}$ are generated with the program GEOPOT (Smith, 1998), which allows to compute the gravity field related quantities at given locations by using mainly a geopotential model. These reference solutions, see Figure 4, are obtained by using a state of the art geopotential model (Förste et al., 2014, EIGEN-6C4) up to degree and order 2000 (10 km resolution). The long wavelengths of the gravity field covered by the satellites and longer than the extent of the local area are removed up to 200 km resolution, providing centered data or close to zero for the determination of local covariance function. The terrain effects are removed with the help of the topographic potential model dV_ELL_RET2012 (Claessens and Hirt, 2013); 2. 2.

The distribution of gravimetric data is chosen to be the same than the control grid and contains 6989 location points. Several spatial distributions of clock, more and less dense, are generated by randomly sampling points from this regular grid, see Figure 5 (top); and synthetic measurements ${\delta{g}}^{\ast}$ and ${T}^{\ast}$ are simulated on these points as previously; 3. 3.

A white noise is added to the simulated data ${\delta{g}}^{\ast}$ and ${T}^{\ast}$, with a standard deviation of 0.1 m2 $s^{-2}$ (i.e. 1 cm on the geoid) for clocks and 1 mGal for gravimetric measurements; 4. 4.

The disturbing potential $\widetilde{T}$ is estimated from the 6989 synthetic measurements $\delta{g}$ only, and from the combination ${T}^{\ast}$ and ${\delta{g}}^{\ast}$ on the 5-km step grid using the Least-Square Collocation (LSC) method (Moritz, 1980). In this step, we make an assumption on the gravity field regularity in the target area, using a logarithmic 3D covariance function (Forsberg, 1987). This model has the advantage to provide the auto-covariances (ACF) and cross-covariances (CCF) of the potential ${T}^{\ast}$ and its derivatives in closed-form expressions. Parameters of this model are adjusted on the empirical ACF of ${\delta{g}}^{\ast}$ with the program GPFIT (Forsberg and Tscherning, 2008). Note that this covariance function contains low frequencies that we have removed from the data in step 1; 5. 5.

Finally, the potential recovery quality is evaluated for all the data distribution sets and types of data by comparing the statistics of the residuals between $\widetilde{T}$ and $\overline{T}$.

Figure 5(b) shows the residuals between the original and the reconstructed potential for denser and denser distributions of clocks. Clock number zero means that only gravity data have been used.

With no clock, the geoid is reconstructed with an accuracy of 3 cm in rms (or 0.3 m2 s*-2* in term of potential variations). When clocks with an accuracy of 1-cm for geoid height differences (or 0.1 m2 s*-2*) are added to the network, they improve the rms by several factors and reduce the bias up to 2–3 orders of magnitude. Also, we can see that it is not necessary to have a dense clock coverage, beyond 10 percent of the clock number of gravity points, to improve the determination of the disturbing potential. In this test, a coverage of around 1–3 percent of the gravity data could be sufficient to reach the centimeter-level precision and improve greatly the bias.

This simple example shows that clock-based geodetic observables provide useful information at spatial scales beyond those of the satellites for geopotential determination. The gain would be even larger in areas of higher relief. A more detailed study discussing the role of different parameters, such as the noise level in the data, effects of the resolution of gravity measurements and modeling errors can be found in (Lion et al., 2017), in which they consider realistic spatial samplings. The addition of gravimetric data outside the modeling target area could also be considered in future works.

7 Space-borne gravimetry with clocks

Sustained observation of the time varying gravity field of the Earth with the highest possible precision and spatio-temporal resolution is crucial for understanding mass transport in the Earth system, which contributes to a better understanding of all kinds of geophysical processes and their relation to climate change, geohazards, etc. (ESA, 2015). The US/German Gravity Recovery and Climate Experiment (GRACE) has shown the tremendous potential of observing mass changes due to all kinds of processes, including continental hydrology (Tapley et al., 2004), ice melt in Greenland and West-Antarctica (Velicogna and Wahr, 2006a, b), and depletion of huge aquifer systems (Richey et al., 2015). The importance of sustained observation is reflected by the realization of the upcoming GRACE Follow-On mission (Sheard et al., 2012).

GRACE-type missions rely predominantly on so-called low-low satellite-to-satellite tracking (ll-SST), where the distance between two trailing satellites, typically with a distance of the order of 100-200 km, is measured with a precision of a few microns for microwave systems (Dunn et al., 2003). This technique is especially suitable for observing mass changes with a temporal resolution of the order of weeks to months at spatial resolutions of at best a few hundreds of kms. The European Space Agency (ESA) Gravity field and steady-state Ocean Circulation Explorer (GOCE) used the concept of satellite gravity gradiometry (SGG), where minute acceleration differences between adjacent capacitive accelerometers (baseline $\approx$ 0.5 m) are observed (Drinkwater et al., 2007). This concepts allows the observation of the gravity field with a spatial resolution of around 70-100 km, but is not well-suited for observing temporal gravity field changes.

Sustaining and at the same time improving the observation of mass transport through space-borne gravimetry has led to the selection of GRACE Follow-on. This mission will not only carry a radiowave instrument, but also a laser link for providing range observations with a precision that is claimed to be a factor 20 better than the GRACE microwave system (Jet Propulsion Laboratory, 2017). In addition, studies have been carried out that demonstrate a great potential for enhancing both the temporal and spatial resolution of flying multiple ll-SST missions in parallel, and at the same time also enhancing the global isotropy and homogeneity of gravity field solutions (Thales Alenia Spazio et al., 2010). However, there is a limit of what will be possible in terms of spatial and temporal resolution when using ll-SST or GOCE-type SGG. Therefore, new concepts are under study, including the use of cold atom interferometry (Carraz et al., 2015), and the use of ultra-precise space-borne atomic clocks (Mayrhofer and Pail, 2012).

7.1 Space-borne high-precision clocks

Here, we analyze the use of clock comparison between a LEO satellite and terrestrial ground stations for the purpose of determining the Earth’s gravity field. The frequency of an electromagnetic signal transmitted from the ground and observed at the spacecraft will be shifted because of both the relative velocity and the difference in gravitational potential between the transmitter and receiver. By combining the measured frequencies, the following observable is obtained (Kopeikin et al., 2011):

[TABLE]

where $T$ and $\tau$ denote the TCG coordinate time and observer proper time, respectively. The 2 and 1 subscripts denote the receiver (spacecraft) and transmitter (ground station). The terms $d\tau/dT$ are given to post-Newtonian order by Eq. (50). The ratio of coordinate times (in TCG) can, to first order, be expressed as follows:

[TABLE]

where $\mathbf{v}_{2}$ and $\mathbf{v}_{1}$ denote the velocity vectors of receiver and transmitter (evaluated at $T_{2}$ and $T_{1}$), and $\mathbf{v}_{21}=\mathbf{v}_{2}-\mathbf{v}_{1}$. In the above, $\hat{\mathbf{k}}$ denotes the unit vector from transmitter position $\mathbf{r}_{1}$ to receiver position $\mathbf{r}_{2}$, so:

[TABLE]

The relation in Eq. (83) is the so-called radial (or first-order) Doppler effect: the direct result of a change in transmission coordinate time on the reception coordinate time. In the following, we will assume that $d\tau/dT$ at the ground station can be modelled/measured to an uncertainty that is much better than for the space segment. The influence on the velocity of the spacecraft on its local $d\tau/dT$ is termed the quadratic (or second-order) Doppler effect.

For this application, we wish to retrieve the information on the Earth’s gravitational potential $V_{E}$ from measurements of ${f_{2}}/{f_{1}}$. However, the measured frequency ratio encodes both velocity and potential differences. To map the influence of both spacecraft velocity uncertainties and the strength of the potential signature to uncertainties in the model for $\Delta(f_{2}/f_{1})$, we write for the first-order errors:

[TABLE]

where the $\Delta\left(*\right)$ denotes the uncertainty in the parameter/measurement $*$. By combining Eq. (87), with Eqs. (85), and (86), respectively, we obtain the uncertainty in the potential determination $\Delta\left(V_{E}\right)$ induced by the velocity uncertainty as a result of the radial and quadratic Doppler effect, due to $\Delta\left(\hat{\mathbf{k}}\cdot\mathbf{v}_{21}\right)$ and $\Delta\left(v_{2}\right)$, respectively.

The influence of the radial Doppler effect is by far the greatest. The maximum allowable error in $\Delta\left(\hat{\mathbf{k}}\cdot\mathbf{v}_{21}\right)$ needed to obtain a given uncertainty $\Delta V_{E}$ for the potential is obtained from Eqs. (85) and (87):

[TABLE]

A potential height error $\Delta h_{g}$ corresponds to a potential error of $\Delta h_{g}\cdot g$, with $g$ the gravitational acceleration. Consequently, for $\Delta h_{g}$=1 cm, we get $\Delta\left(\hat{\mathbf{k}}\cdot\mathbf{v}_{21}\right)=\frac{0.01g}{c}\approx 3\times 10^{-10}$ m/s or 0.3 nm/s (assuming a LEO satellite flying at 7.5 km/s with $g\approx$ 9 m/s2). This value is well below the orbit uncertainty of state-of-the-art orbit determination. However, it is possible to largely eliminate the radial Doppler shift by using a coherent two-way optical link (Djerroud et al., 2010; Blanchet et al., 2001). It is not possible, however, to eliminate the influence of the quadratic Doppler effect, since we cannot separate the influence of $v$ and $V_{E}$ in the equation for $(d\tau/dT)_{2}$ by combining a one- and two-way link. Again assuming an influence of maximum 1 cm on the geoid determination, the quadratic Doppler term results in a required spacecraft velocity uncertainty of better than $\frac{0.01\times g}{7.5\times 10^{3}}\approx 12~{}\mu$m/s.

In light of the influence of the radial Doppler effect, it will be crucial to have a two-way link, in addition to the one-way frequency comparison. In the error simulations below, it is therefore assumed that such a two-way link is available.

7.2 Accuracy assessment

A hypothetical mission is defined where a satellite is flying in a low Earth orbit carrying an ultra-precise clock that can be compared continuously with perfect reference clocks. It is assumed that this satellite flies in a circular orbit, where its projection on the surface of the Earth (the so-called ground track) repeats after a certain number of nodal days. For such an observation geometry, use can be made of an efficient error propagation method that was introduced in (Colombo, 1984). This error propagation tool has been compared and validated by comparison with robust numerical integration methods and end-to-end simulations for several observation types, including orbit perturbations, SGG and ll-SST (Schrama, 1991; Visser et al., 2001).

The satellite is assumed to fly in a polar orbit (inclination $i=90^{\circ}$) for which the ground track repeats every month (31 days) in which the satellite completes 497 orbital revolutions. The height above the Earth’s surface of the satellite is about 250 km. The space-borne clock is assumed to have a frequency dependent instability equal to $30\times 10^{-18}/\sqrt{\tau}$, where $\tau$ represents the clock integration time. This instability is about equivalent to a stability of $10^{18}$ for an integration time $\tau$ or 1000 s, cf. (Giorgetta et al., 2013). Assuming the clock instability to be the only source of measurement error $\Delta\left(f_{2}/f_{1}\right)$, Eq. (87) leads to the following error for the observation of the gravity potential:

[TABLE]

According to Eq. (86), the velocity then has to be known with a precision of about:

[TABLE]

Determining the LEO satellite velocity with such a precision is very challenging, but might be feasible with future optical links (Chiodo et al., 2013). The requirement for the velocity determination is based on the assumption that the velocity uncertainty causes an error in the potential that is equal or smaller than the measurement error due to clock instability.

The gravitational potential is expanded in terms of spherical harmonics. The uncertainty in the gravity field coefficients is then written as follows, where we neglect the uncertainty in $r$, $\lambda$ and $\phi$:

[TABLE]

where $\mu$ is the gravity parameter (the product of the universal gravitational constant $G$ and Earth’s mass $M$), $a_{e}$ is the mean equatorial radius, $r,\phi,\lambda$ are the spherical coordinates (radius, geocentric latitude, longitude) denoting the location of the space-borne clock, $\bar{P}_{lm}$ is the normalized Legendre polynomial of spherical harmonic degree $l$ and order $m$, and $\bar{\Delta C}_{lm},\bar{\Delta S}_{lm}$ represent the estimated gravity field coefficients.

For a circular satellite orbit with repeating ground track, Eq. (91) can be represented by a Fourier series (Visser et al., 2003) and the method outlined in (Colombo, 1984) can be used for predicting the gravity field retrieval performance. The gravity field can be determined from time series of observed or determined satellite velocities as well. In that case, use is made of a linear orbit perturbation theory, where the associated equations can also be represented by a Fourier series (Schrama, 1991). Then also the method in (Colombo, 1984) can be used for a gravity field retrieval accuracy assessment. Please note that in case of deriving gravity from satellite velocity perturbations, it is assumed that non-gravitational accelerations are known through e.g. observation by precise accelerometers or absent due to a drag-free control system.

The predicted gravity field retrieval performances according to the error spectra in Eq. (89) and (90) are included in Figure 6. It can be observed that the very precise knowledge of the satellite velocity leads to better gravity field retrieval performance through classical orbit perturbation analysis as compared to the post-Newtonian derivation of potential height from clock observations. For comparison, the predicted gravity field performances for the same LEO satellite are included if equipped with a high-quality Global Navigation Satellite System (GNSS) receiver allowing observation of orbit position perturbations with 1 cm precision, and a gradiometer based on cold atom interferometry (0.4 mE flat noise). We note that this assumption on the uncertainty is optimistic and does not reflect the state of the art for a space instrument, which is about one order of magnitude worse. With the used assumptions, the ultra-stable clock does show a much better gravity field performance than the GPS receiver, and also the gradiometer based on cold atom inferometry up to spherical harmonic degree 30. Figure 6 also includes a performance prediction for a clock that would have a flat error spectrum with a stability of $10^{18}$ leading to 1 cm error in potential height. Such a - highly speculative - performance would allow the observation of time variable gravity with a precision that is relatively close to the currently claimed performance level for GRACE Release 5 Level-2 products ((Dahle et al., 2013), see also Figure 6).

Being able to determine the 3-dimensional velocity of the LEO satellite with a precision according to Eq. (90) leads to a performance that is starting to become competitive with the currently claimed GRACE performance and allows the observation of, for example, gravity field changes due to continental hydrology (signal magnitude also indicated in Figure 6, (Gruber et al., 2011)). In this preliminary analysis, we have not analyzed the possibilities to decouple the signature of the spacecraft’s velocity and gravitational potential from the frequency comparison measurements. Such an approach, as is discussed by Dirkx et al. (2016) for planetary laser ranging data, could be used to quantify more robustly the influence of velocity uncertainties on potential measurements using clocks. Conversely, such an analysis would quantify the degree to which the relativistic Doppler shift can be used to dynamically determine the spacecraft’s state when making use of clock comparisons.

8 Discussion and conclusions

Various geodetic applications of highly accurate and stable clocks at the $10^{-18}$ level for the relative frequency shift (corresponding to 1 cm in height) have been addressed. We discussed both the prospects of implementing in practice the new clock-based measurement concepts and their challenges. Our firm conclusion is that the quantum clocks and time metrology will provide us with unique opportunities to support and largely extend the tools for gravity field measurements and physical height determination based on the direct access to the gravity potential differences.

We revisited the relevant time scales used in geodesy and showed the existence of growing inconsistencies in their definition and realization at a $10^{-18}$ accuracy level. This mainly concerns the use of the (non-constant) gravity potential value $W_{0}$ on the geoid, and its relation to the international time scales TT and TAI. In a foreseeable future, clarification and re-definition of the currently-used time concepts will be required to make them consistent with the highly-precise clock measurements.

We have provided a set of formulas to directly use the observed frequency shifts for the determination of physical height differences, where the relative clock measurements replace the differences of geopotential numbers. Applying this conceptually-new approach can facilitate or even fully resolve the problem of the systematic discrepancies existing between various height systems on the regional, national and intercontinental level. Clocks can connect tide gauges through satellite and/or fiber links providing the vertical datum of national height systems and help to determine the present offsets between them. Even discrepancies between independent geoid solutions stemmed from applying different classical techniques (GNSS/levelling versus gravimetric methods) can be determined and eliminated by making use of the precise chronometric measurements. These discrepancies are known to reaching centimeters to decimeters and can well be recovered already with the present-day clock technology. At the envisaged cm level of accuracy, time-variable parts affecting the clock measurements due to observer’s height, position or Earth’s density variations have to be considered. The most prominent effect is caused by the solid Earth tides but smaller variations, such as ocean tides or inaccurate knowledge of the Earth rotation parameters, play a significant role and are to be taken into consideration in the data processing algorithms.

The direct use of observed gravity potential values for a regional gravity field recovery has been studied as a test case in the Massif Central region crossing a number of mountain ranges. We have clearly demonstrated in our simulation scenario that adding the clock-based potential values to the existing data set would notably improve the final gravity field solution. The bias to a reference solution and total rms error could be reduced remarkably well, up to a few orders of magnitude. Another important conclusion stemming from our simulations is that in solving the problem of gravity field recovery it is not required to have a dense clock network, only a very few percent of clock measurements compared to the number of needed gravity data is sufficient.

The possible use of ultra-precise space-borne clocks has been revisited for deriving potential height differences at satellite altitude and using these satellite data for recovering the global gravity field. It has been shown that extremely challenging requirements have to be met in order to use space-borne clocks to successfully fulfill this task. In order to be able to exploit the gravity field information content of clocks with a stability of $10^{-18}$, precise reference clocks are needed along with two-way frequency transfer for eliminating the Doppler effect correction errors. In addition, the position and velocity of the clock need to be known with a very high precision. The current technological level of GNSS receivers allow for determining the position with the required cm-level precision. However, the velocity needs to be determined with such a high precision, that the associated velocity perturbations caused by the gravity field anomalies might be used as well as a fundamental observable for determining Earth’s gravity field. Finally, the $10^{-18}$ clock stability needs to be achievable for very short integration times of the order of seconds. Even if this would be possible and achieved, the existing space geodesy concepts such as GRACE-type ll-SST lead to a better observability of Earth’s time variable gravity. However, our study does show that the frequency transfer may be used to improve the determination of the spacecraft velocity.

Acknowledgements.

The International Space Science Institute (ISSI), is acknowledged for having provided the opportunity of presenting the work described in this paper at the Workshop on High Performance Clocks, 30 November - 4 December 2015 in Bern, Switzerland.