Determination of Gravitational Potential at Ground Using Optical-Atomic Clocks on Board Satellites and on Ground Stations and Relevant Simulation Experiments

Abstract

The general relativity theory provides a potential way to directly determine the gravitational potential (GP) difference by comparing the running rate or vibration frequencies of two optical-atomic clocks located at two stations. Recently we proposed an approach referred to as satellite frequency signal transmission based on the Doppler canceling technique or tri-frequency combination technique to determine the GP difference between a satellite and a ground site via exchanging microwave signals. Here, as an extension of our previous study, we aim to formulate determination of GP at ground stations and establish simulation experiments in different cases, including determining the GP at a ground station via one or more satellites and determining the GP difference between two ground stations via one or more satellites. Concerning each case we made simulating experiments, and results show that the precision of the GP at a ground station and that of the GP difference between two stations, determined via one satellite, are, respectively, about 0.383 and 0.454 m$^{2}$/s$^{2}$, assuming the clocks with inaccuracy of about $1\times 10^{-18}$ (s/s) level are available. If more satellites equipped with ultra-high-precise clocks are available, the precision of the determined GP (difference) at ground stations can be further improved.

1 Introduction

According to the theory of general relativity (GR), an atomic (or optical atomic) clock’s running rate and its vibration frequency will change at different positions with different gravitational potentials (Einstein 1915; Weinberg 1972). Conversely, one can determine the gravitational potential (GP) at a space point or on ground by measuring the change of clocks’ running rates (Bjerhammar 1985) or by measuring the change of electromagnetic signals’ frequencies (Shen et al. 1993). These alternative approaches of determining GP (difference), referred to as the clock transportation comparison (CTC) and frequency signal transmission comparison (FSTC), respectively, require clocks or oscillators with ultra-high precision, say $1\times 10^{-18}$, which is equivalent to 1 cm in height. The time-frequency related confirmation of the GR by various studies (Pound and Rebka 1959; Pound and Snider 1965; Hafele and Keating 1972; Vessot and Levine 1979; Turneaure et al. 1983; Chou et al. 2010) provides a potential and prospective way to directly determine the GP based on the CTC and FSTC.

In recent years, with quick development of high-precision clock manufacturing technology, the optical-atomic clocks (OACs) with relative instability around $10^{-18}$ in several hours and inaccuracy of $10^{-18}$ level have been generated in the laboratory (Hinkley et al. 2013; Bloom et al. 2014; Ushijima et al. 2015), and OACs with such precision level are promising to be installed on satellites in the near future (Schiller et al. 2007; Tino et al. 2007). Since the current precision level of OACs is sufficient for applying the CTC or FSTC approach to determining GP, it attracts more and more attention from geodesy, geoscience and academia (Brumberg and Groten 2001; Pavlis and Weiss 2003; Bondarescu et al. 2015). By far there are generally three kinds of methods that apply the GR to GP determination: (1) transport clocks between two stations on ground and determine the GP difference between the two stations by measuring the accumulated difference of the clocks ticks (Bjerhammar 1985), (2) connect two stations by optical fiber or coaxial cable and transmit frequency signals or time signals between the two stations (Shen and Peng 2012; Shen 2013a, b; Shen and Shen 2015), (3) transmit frequency signals among different stations on ground via GNSS-type (or communication-type) satellites (Shen et al. 1993, 2011; Shen 1998; Shen and Ning 2005).

Although the three kinds of methods mentioned above are all showing potential of determining GP, the first two kinds have obvious drawbacks. For example, the clock transportation comparison approach (Bjerhammar 1985) is laborious and time-consuming, and the errors induced by transportation are difficult to control. The cable time transfer approach (Shen and Shen 2015) or the fiber frequency transfer method (Shen and Peng 2012; Shen 2013a, b) is constrained by the distance between the two stations and increases the complexity especially in the cases that we need to connect stations separated by ocean and mountainous areas. Although the fiber frequency transfer comparison has reached fairly high precision, about $10^{-19}$ level in relative accuracy (Grosche et al. 2009; Calonico et al. 2014), the requirement of fibers limit its application in geodesy because we cannot conveniently determine the GP at an arbitrary position. As contrast, the third kind of method is most flexible and promising, since we can bridge any two places with one or several satellites. It is less laborious, fast and unlimited to geography and distance.

The third method of the GP determination requires FSTC between ground and a satellite. Currently, most relative experiments and researches aim to validate the gravitational redshift effect predicted by the GR. Conversely, if the GR is proved to be reliable or accurate enough, we can determine the GP based on the gravitational redshift effect. The first ground-satellite frequency transfer experiment is the Gravity Probe A (GP-A) experiment in 1976 (Vessot and Levine 1979), which aims to test the GR at $10^{-15}$ level in frequency accuracy. And this experiment is the most precise direct test of the gravitational redshift to date. Delva et al. (2015) proposed to test the gravitational redshift using Galileo satellites with the frequency precision of $10^{-16}$, but when the experiments will be put into practice remains uncertain. The next ground-satellite experiment similar to GP-A is the future Atomic Clock Ensemble in Space (ACES) experiment planned to fly on the International Space Station in 2017 (Cacciapuoti and Salomon 2011). The ACES will carry out both the time and frequency transfer experiment to test the gravitational redshift, and the precision of frequency transfer is supposed to be at $10^{-17}$ level. Furthermore, the STE-QUEST (Altschul et al. 2014) project, planned to launch in 2024, is proposed to equip an optic-atomic clock with the stability of $1\times 10^{-18}$ in few hours to test the gravitational redshift. Hence, the third method is very prospective in the near future.

Compared to time and frequency transfer on ground, the transfer between ground and a satellite confronted much more problems and challenges. For example, the atmosphere and ionosphere will cause signal delay and frequency shift, and the Earth rotation and tidal effect also impose influence to the time and frequency transfer. When the accuracy requirement of frequency transfer reaches $10^{-18}$ level or even higher, instead of the clock stability, the systematic errors might be the dominant error sources. Wolf and Petit (1995) detailedly studies in detail the clock synchronization in the vicinity of the Earth at the accuracy level of $10^{-18}$. They analyzed the influence of various error sources, including tidal effect, Doppler effect, external masses (Sun, Moon and other planets), atmosphere pressure, polar motion and so on, and the error introduced by each of these factors is below $1\times 10^{-18}$ level after correction. But they did not consider the frequency shifts caused by ionosphere and troposphere, which also need to be corrected. Ashby (1998) and Blanchet et al. (2001) reexamined the GP-A test and improved the frequency transfer equation by introducing the $c^{-3}$ terms to the accuracy level of $5\times 10^{-17}$. They also analyzed the influence of Shapiro time delay (Shapiro 1964) in frequency transfer. Considering the GP-A test, Linet and Teyssandier (2002) formulated a frequency shift in a gravitational field generated by an axisymmetric rotating body and provided a one-way frequency transfer equation accurate to $c^{-4}$ terms, equivalent to the relative accuracy level higher than $1\times 10^{-18}$ in frequency. Shen et al. (2016) analyzed the ionosphere and troposphere influences to the frequency links between a ground station and a GNSS-type satellite, and the introduced errors can be reduced to $10^{-19}$ level after proper correction.

Thus, following our previous idea (Shen et al. 1993, 2011; Shen 1998; Shen and Ning 2005), we formulated an approach using frequency signals links based on the Doppler canceling technique (DCT, see Kleppner et al. 1970; Vessot and Levine 1979) or tri-frequency combination (TrFC) technique to practically realize the determination of the GP difference between a satellite and a ground station (Shen et al. 2016), which is referred to as satellite frequency signal transmission (SFST). Based on our theoretical formulation, the SFST can reach 1 m$^{2}$/s$^{2}$ if the clocks’ inaccuracy can achieve $1\times 10^{-17}$ level (Shen et al. 2016). Here we will extend the study of Shen et al. (2016) to the accuracy level of $10^{-18}$, focusing on determining the GP at ground at centimeter level and conduct relevant simulation experiments to show how to realize the GP determination based on the SFST.

2 Gravitational Potential Difference Determination Between a Satellite and a Ground Site

Referring to Fig. 1, the SFST contains three microwave links (Shen et al. 2016). An emitter at a ground station P emits a frequency signal $f_{\mathrm{e}}$ at time $t_1$. When the signal is received by a satellite S at time $t_2$, it immediately transmits the received signal $f_{\mathrm{e}}^{\prime }$ and emits a frequency signal $f_{\mathrm{s}}$ at the same time. These two signals transmitted and emitted from the satellite are received by a receiver at the ground station P at time $t_3$. During the period of the emitting and receiving, the position of the ground station in space has been changed from P to $P^{\prime }$ (see Fig. 1).

Fig. 1

(modified after Shen et al. 2016)

Ground station P emits a frequency signal $f_{\mathrm{e}}$ at time $t_1$. Satellite S transmits the received signal $f_{\mathrm{e}}^{\prime }$ and emits a new frequency signal $f_{\mathrm{s}}$ at time $t_2$. The ground station receives signals $f_{\mathrm{e}}^{\prime \prime }$ and $f_{\mathrm{s}}^{\prime }$ at time $t_3$ at position $P^{\prime }$. $\phi$ is GP, ${\mathbf {r}}$ is position vector, ${{\mathbf {v}}}$ is velocity vector, ${\mathbf {a}}$ is centrifugal acceleration vector

Based on the procedures as described above (also see Fig. 1), we can extract the gravity frequency shift signals (or equivalently gravitational frequency shift signals). Suppose we set a basic frequency $f_0$ and $f_{\mathrm{e}} = f_{\mathrm{s}} = f_0$, then the frequency shift signals can be determined as depicted in Fig. 2. The frequencies of the signals emitted from ground oscillator and satellite oscillator are $f_0$. The microwave link 1 and link 2 consist a go-return link by a phase-coherent microwave transponder equipped at the satellite and provide two-way frequency shift data as a beat frequency $f_0^{\prime \prime }-f_0$ (Shen et al. 2016). Similarly, the microwave link 3 provides one-way frequency shift data as a beat frequency $f_0^{\prime }-f_0$ (Shen et al. 2016). The output frequency $\Delta f$ is expressed as (Kleppner et al. 1970; Vessot and Levine 1979; Vessot et al. 1980; Shen et al. 2016):

$$\begin{aligned} \Delta f = f_0^{\prime }-f_0-\frac{f_0^{\prime \prime }-f_0}{2}. \end{aligned}$$

(1)

Fig. 2

(modified after Vessot and Levine 1979; Shen et al. 2016)

Ground oscillator emits a frequency signal $f_0$ to the spacecraft (or satellite), then the spacecraft transmits the received signal to ground and emits a frequency signal $f_0$ from spacecraft oscillator to the ground at the same time

In free space, the output frequency expressed by Eq. (1) implies that it can completely cancels the first-order Doppler effect. This is the reason that this procedure is referred to as Doppler canceling technique (DCT, see, e.g., Vessot and Levine 1979). Hence, the GP difference between the satellite and the ground site can be obtained from the following equation (Vessot and Levine 1979):

$$\begin{aligned} \frac{\Delta f}{f_0} = \frac{\phi _{\mathrm{s}} - \phi _{\mathrm{e}}}{c^2} - \frac{\left| {{\mathbf {v}}}_{\mathrm{e}} - {{\mathbf {v}}}_{\mathrm{s}} \right| ^2}{2c^2} - \frac{{\mathbf {r}}_{\mathrm{se}}\cdot {\mathbf {a}}_{\mathrm{e}}}{c^2} \end{aligned}$$

(2)

Here an Earth-centered inertial coordinate frame has been applied in Eq. (3), where $\phi _{\mathrm{s}}$ and $\phi _{\mathrm{e}}$ are Newtonian GPs at spacecraft (or satellite) and ground station, respectively, ${{\mathbf {v}}}_{\mathrm{e}}$ and ${{\mathbf {v}}}_{\mathrm{s}}$ are velocities of ground station and spacecraft, respectively, ${\mathbf {r}}_{\mathrm{se}}$ is vector from spacecraft to ground station, ${\mathbf {a}}_{\mathrm{e}}$ is centrifugal acceleration vector of ground station, c refers to the speed of light in vacuum.

Equation (2) has omitted the terms higher than $c^{-2}$, and it holds only at the accuracy level or a little better than $10^{-15}$ (Cacciapuoti and Salomon 2011). For a higher precision requirement, Ashby (1998) and Blanchet et al. (2001) appended the $c^{-3}$ terms to Eq. (2) and established an equation suitable for an accuracy level of $5\times 10^{17}$. However, to achieve an accuracy level of $1\times 10^{-18}$, terms up to $c^{-4}$ should be considered. A theoretical formula of one-way frequency transfer in free space accurate to $1\times 10^{-18}$ was given by Linet and Teyssandier (2002). Based on the study of Linet and Teyssandier (2002), with three-link frequency transmission as described in Fig. 1, Eq. (1) can be expressed as (derived in the Appendix in detail)

$$\begin{aligned} \frac{\Delta \phi _{\mathrm{es}}}{c^2}\equiv \frac{\phi _{\mathrm{s}} - \phi _{\mathrm{e}}}{c^2}= \frac{\Delta f}{f_0} - \frac{v_{\mathrm{s}}^2 - v_{\mathrm{e}}^2}{2c^2} - \sum ^4_{i=1} q^{(i)} \end{aligned}$$

(3)

where the frequency shift “output” $\Delta f$ is given by expression (1), $\phi$ is the Earth’s Newtonian GP, $\Delta \phi _{\mathrm{es}} = \phi _{\mathrm{s}} - \phi _{\mathrm{e}}$ is the GP difference between satellite and ground station, $q^{(i)} (i=1,2,3,4)$ are referred to Eqs. (50) and (29)–(32), and their explanations are provided thereafter. Equation (3) takes a different form from Eq. (2) because the latter aims only to the accuracy level of $c^{-2}$ (see Vessot and Levine 1979), while the former keeps all terms accurate to $c^{-4}$.

We note that, based on the Doppler canceling technique (see Fig. 2), an oscillator (clock) with stability of $10^{-18}$ is necessary to control its emitting frequency. Then, by tri-frequency combination we may cancel out the Doppler effect and precisely draw out the GP difference between a ground station and the spacecraft. This is the reason why a precise clock on board a spacecraft is needed.

Equation (3) has included the effect of Shapiro delay and the effect of the axisymmetric rotating body of the Earth. It can reach the accuracy of $10^{-19}$ level, but holds only in free space. In real space outside the Earth, a signal’s frequency will be contaminated by ionospheric and tropospheric effects and other influences (Shen et al. 2016), which means that certain corrections should be further appended to the equation. In addition, the potential difference $\phi _{\mathrm{se}}$ contains the influences of other factors (such as tidal effect, gravitational potential fields of celestial bodies). If we aim to obtain the GP difference caused by the Earth, such influences should also be removed. The residual errors after all of the corrections, together with the systematic errors (such as equipment errors, orbit uncertainty), are the main factors that determine the final accuracy of the frequency transfer based on the tri-frequency combination (TrFC) technique.

When various corrections and influences are taken into consideration, Eq. (3) is modified as the following equation

$$\begin{aligned} \frac{\Delta \phi _{\mathrm{es}}}{c^2} \equiv \frac{\phi _{\mathrm{s}}-\phi _{\mathrm{e}}}{c^2} = \frac{\Delta f}{f_0} - \frac{v_{\mathrm{s}}^2 - v_{\mathrm{e}}^2}{2c^2} - \sum ^4_{i=1} q^{(i)} + \Lambda f + \delta f \end{aligned}$$

(4)

where $\Lambda f$ is the sum of all correction terms, $\delta f$ is the sum of all error terms.

The correction term $\Lambda f$ in Eq. (4) is expressed as

$$\begin{aligned} \Lambda f = \Lambda f_{\mathrm{ion}} + \Lambda f_{\mathrm{tro}} + \Lambda f_{\mathrm{tide}} + \Lambda f_{\mathrm{celes}} \end{aligned}$$

(5)

where $\Lambda f_{\mathrm{ion}}$ and $\Lambda f_{\mathrm{tro}}$ are, respectively, the corrections of ionospheric and tropospheric effects (after tri-frequency combination), $\Lambda f_{\mathrm{tide}}$ is the contribution of the additional potential associated with the Earth’s deformation caused by tidal effects, $\Lambda f_{\mathrm{celes}}$ is caused by the GP generated by the main celestial members in our solar system, including the Sun, the Moon and other planets. Accordingly, after the corrections as expressed as (5), the total residual errors are expressed as

$$\begin{aligned} \delta f_{\mathrm{cor}} = \delta f_{\mathrm{ion}} + \delta f_{\mathrm{tro}} + \delta f_{\mathrm{tide}} + \delta f_{\mathrm{celes}} \end{aligned}$$

(6)

Thus the error terms $\delta f$ in Eq. (4) can be expressed as

$$\begin{aligned} \delta f = \delta f_{\mathrm{cor}} + \delta f_{\mathrm{sys}} \end{aligned}$$

(7)

where $\delta f_{\mathrm{sys}}$ is the sum of all relevant systematic errors, which will be discussed later.

The correction terms $\Lambda f_{\mathrm{ion}}$ and $\Lambda f_{\mathrm{tro}}$ have been studied in detail in Shen et al. (2016), expressed as

$$\begin{aligned} \Lambda f_{\mathrm{ion}} = \frac{78{,}570{\bar{\rho }} \left| {\mathbf {r}}_{\mathrm{se}}\right| \left( {{\mathbf {v}}}_{\mathrm{s}} - {{\mathbf {v}}}_{\mathrm{e}} \right) \cdot {\mathbf {a}}_{\mathrm{e}}}{c^2f_0^2H\left| {{\mathbf {v}}}_{\mathrm{s}} - {{\mathbf {v}}}_{\mathrm{e}} \right| } \end{aligned}$$

(8)

and

$$\begin{aligned} \Lambda f_{\mathrm{tro}} = -\frac{60\left( {\bar{M}}_1 + {\bar{M}}_2 \right) \left| {\mathbf {r}}_{\mathrm{se}}\right| \left( {{\mathbf {v}}}_{\mathrm{s}} - {{\mathbf {v}}}_{\mathrm{e}} \right) \cdot {\mathbf {a}}_{\mathrm{e}}}{c^2H\left| {{\mathbf {v}}}_{\mathrm{s}} - {{\mathbf {v}}}_{\mathrm{e}} \right| }, \end{aligned}$$

(9)

where H is the height (in km) of the spacecraft from the ground, ${\bar{\rho }}$ is the average electron density (in m⁻³), M ₁ and M ₂ are substitutions for simplification, defined as $M_1 = 77.6 \times 10^{-6}p/T, M_2 = 0.373\varepsilon /T^2$, and $\bar{M}_1$ and ${\bar{M}}_2$ are the average value of $M_1$ and $M_2$ along the signals’ propagation paths (see Shen et al. 2016), where $p, T, \varepsilon$ are, respectively, total pressure (in mbar), temperature (in degrees K), and partial pressure of water vapor (in mbar); ${\mathbf {a}}_{\mathrm{e}}$ is the acceleration vector of the ground station. The magnitude of correction terms $\Lambda f_{\mathrm{ion}}$ and $\Lambda f_{\mathrm{tro}}$ and their residual errors after corrections are listed in Table 2 (which is explained later).

The deformation of Earth will cause the potential changes outside the Earth and the position changes on the surface of the Earth. These two effects consist of the correction term $\Lambda f_{\mathrm{tide}}$, which are expressed in spherical harmonics expansion series (Farrell 1972). The tide-induced potential changes in the free space are most conveniently modeled as variations in the standard geopotential coefficients $C_{\mathrm{nm}}$ and $S_{\mathrm{nm}}$ (Eanes et al. 1983), and their contributions can be estimated from some global tide models (e.g., Parke 1982). They can also be calculated by some mature softwares (Tsoft for example, see Camp and Vauterin 2005), and the residual error is at the millimeter level.

Concerning the last correction term $\Lambda f_{\mathrm{celes}}$, the GP influence of each planet [except for the Sun, which is expressed as Eq. (11)] can be expressed as:

$$\begin{aligned} V_i = \frac{{\mathrm{GM}}_i}{\left| {\mathbf {r}}_i + {\mathbf {r}}_{\mathrm{se}} \right| } - \frac{{\mathrm{GM}}_i}{r_i} + O_i, (i = {\mathrm{Moon, Mercury, Venus}}, \ldots ) \end{aligned}$$

(10)

where G is gravitational constant, $M_i$ is the mass of the celestial body, ${\mathbf {r}}_i$ is the vector from a celestial body to the ground station, ${\mathbf {r}}_{\mathrm{se}}$ is the vector from the ground station to the satellite. $O_i$ is the higher-order potentials caused by non-spherical distribution. However, Eq. (10) is not suitable for the Sun because of the equivalence principle (Kleppner et al. 1970; Hoffmann 1961). The GP influence of the Sun should be expressed as (Hoffmann 1961):

$$\begin{aligned} V_{\mathrm{Sun}} = {\mathrm{GM}}_{\mathrm{Sun}}\left( \frac{1}{\left| {\mathbf {r}}_{\mathrm{c}} + {\mathbf {r}}_{\mathrm{sat}}\right| } - \frac{1}{r_{\mathrm{c}}} + \frac{{\mathbf {r}}_{\mathrm{sat}}\cdot \mathbf {i}_{\mathrm{c}}}{r_{\mathrm{c}}^2}\right) - {\mathrm{GM}}_{\mathrm{Sun}}\left( \frac{1}{\left| {\mathbf {r}}_{\mathrm{c}} + {\mathbf {r}}_{\mathrm{grd}}\right| } - \frac{1}{r_{\mathrm{c}}} + \frac{{\mathbf {r}}_{\mathrm{grd}}\cdot \mathbf {i}_{\mathrm{c}}}{r_{\mathrm{c}}^2}\right) + O_{\mathrm{Sun}} \end{aligned}$$

(11)

where ${\mathbf {r}}_{\mathrm{c}}$ is the vector from Sun to the Earth’s mass center, $\mathbf {i}_{\mathrm{c}}$ is the unit vector of ${\mathbf {r}}_{\mathrm{c}}$, ${\mathbf {r}}_{\mathrm{sat}}$ and ${\mathbf {r}}_{\mathrm{grd}}$ are, respectively, the vectors of the satellite and ground station with respect to the Earth’s mass center. Because of the long distances and relatively small gravitational influences, we can omit the $O_{i}$ terms in Eqs. (10) and (11). Estimations show that $O_{i}$ does not exceed $10^{-3}$ m$^{2}$/s$^{2}$, equivalent to $10^{-20}$ in frequency influence (see Table 1). The ephemeris of solar system planets can be obtained from Ephemerides of Planets and the Moon (Pitjeva 2013), and the errors of planets’ orbit determination are negligible in our estimation (even an orbit offset of 1 km for the Moon causes a frequency error at the level of $10^{-20}$, and for the other planets are even smaller). For the potential difference measurement between a GNSS-type satellite and a ground station, the largest correction magnitude of each celestial body (when the body, the satellite and the ground station are located in one straight line) is listed in Table 1. We can see that to achieve the accuracy level of $1\times 10^{-18}$ for measuring the Earth’s GP difference, all of other celestial bodies except for Neptune need to be considered. Then, after the celestial bodies’ corrections, the residual errors $\delta f_{\mathrm{celes}}$ are below $1\times 10^{-20}$ (Table 2).

Table 1 Largest correction magnitude of each celestial body for the GP measurement between a GNSS satellite and a ground station

Table 2 Error magnitudes of different error sources in determining GP difference between a satellite and a ground station

The error term $\delta f_{\mathrm{sys}}$ in Eq. (7) is caused by all the effects that cannot be properly or effectively modeled and corrected, such as the equipment delays, clock errors, satellite’s orbit errors. Here we denote $\delta f_\mathrm{{sys}}$ as

$$\begin{aligned} \delta f_{\mathrm{sys}} = \delta f_{\mathrm{vepo}} + \delta f_{\mathrm{delay}} + \delta f_{\mathrm{osc}} + \delta f_{\mathrm{o}} \end{aligned}$$

(12)

where the relevant terms are explained in what follows.

$\delta f_{\mathrm{vepo}}$ is the position and velocity errors of ground station and satellite. The position error in the precise ephemeris of a GPS satellite is about $10^{-2}$ m (Kang et al. 2006; Guo et al. 2015), and the velocity error can be reduced to below $10^{-5}$ m/s (Sharifi et al. 2013). The position error of a ground station is negligible because it is relatively small compared to a satellite. Then, the errors introduced from position vector can be estimated by applying error propagation to Eq. (3), and the amount of $\delta f_{\mathrm{vepo}}$ is below $3.4\times 10^{-19}$ (see Table 2). In addition, Duchayne et al. (2009) have studied the influence of position difference between reference point and mass center of a satellite and concluded that the introduced error is totally negligible.

$\delta f_{\mathrm{delay}}$ is introduced by all the equipment delays. Notice that the DCT method performs frequency transfer, without involving time transfer, thus the hardware time delays in both ground station and satellite can be neglected. But satellite’s transponder delay must be taken into consideration, because the satellite is in motion, its position when receiving signals is different from that when emitting signals. For a transponder’s delay at about 800 ns (Pierno and Varasi 2013), its introduced error $\delta f_{\mathrm{delay}}$ is at the level of $10^{-19}$ (Shen et al. 2016).

$\delta f_{\mathrm{osc}}$ is oscillator (clock) error. Since our aim is to determine the GP difference to accuracy of centimeter level, which means that clocks at least with stability and accuracy of $10^{-18}$ level are required. In this paper we assume that the clocks meet our requirement (some day in the near future), and their instability can achieve $1\times 10^{-18}$ in an hour ($\delta f_{\mathrm{osc}}<10^{-18}$). It should be noted that such high-precise optical-atomic clock has been realized in the laboratory (Bloom et al. 2014). Although currently the stablest clocks onboard satellites are at $10^{-17}$ stability level (Cacciapuoti and Salomon 2011), the stability of $10^{-18}$ level (onboard satellites) will be achieved in the near future. And since an atomic clock is sensitive to temperature and magnetic field (see, e.g., Rochat et al. 2012), it is also important to stabilize the inner environment of a satellite to minimize the introduced errors. In this paper, we do not discuss these influences because it is a topic of clock manufacturing, and different clocks might vary in sensitivity to temperature and magnetic field.

Finally, the term $\delta f_{\mathrm{o}}$ denotes all of the higher-order contributors (multi-path effects, polar motion, etc.) that can be safely neglected. The magnitudes of all correction terms and error terms are listed in Table 2. We can see that the magnitudes of some of the error sources are different from those given by Wolf and Petit (1995), due to the fact that we focus on frequency comparison between satellite and ground links, while they studied the clock synchronization at ground or satellite. And some of error estimates in Wolf and Petit (1995) are undetailed (for example, the influence of celestial bodies).

It should be noted that the DCT method is designed for microwave links, because the measurements concern frequency comparison. To our knowledge, in free or medium space it is not suitable for optical links, which are mainly used for time transfer (e.g., laser time transfer). In addition, there are some constraints for the frequency of $f_0$. A higher value of $f_0$ helps to reduce the influence of the ionosphere [see Eq. (8)] and reduce the refraction of propagation path. However, if the $f_0$ value is too high (>30 GHz for example), the energy required for sustaining the system greatly increases, and the signal will strongly attenuated by the Earth’s atmosphere and particles contained in it, especially during wet weather. In practice, a frequency band range from 2 GHz (adopted in the GP-A experiment) to 15 GHz (adopted in the ACES mission) is suitable.

3 Determination of GP Difference Between Two Ground Sites

The SFST approach (Shen et al. 2016) was designed for determining the GP difference between a satellite and a ground site, as described in Sect. 2. However, if two ground sites are connected to the same satellite via satellite links simultaneously, the satellite can serve as a “bridge” to connect the two ground sites (Shen et al. 1993, 2016). Thus the GP difference between the two ground sites can be determined. Figure 3 depicts the concept. The ground stations $P_1$ and $P_2$ simultaneously observe a satellite S, which is at the same time visible by these two ground stations.

Fig. 3

Links of frequency signals among a satellite and two ground stations. The satellite S receives frequency signals from two ground stations $P_1$ and $P_2$ simultaneously and then transmits the signals back to ground stations. The ground stations receives the transmitted frequency signals at $P_1^{\prime }$ and $P_2^{\prime }$ because of Earth rotation

According to Eq. (3), for each of the ground stations $P_1$ and $P_2$, we have the following equations:

$$\begin{aligned} \frac{\Delta \phi _{\mathrm{e}1\mathrm{s}}}{c^2} \equiv \frac{\phi _{\mathrm{s}}-\phi _{\mathrm{e}1}}{c^2}& = \frac{\Delta f_1}{f_0} - \frac{v_{\mathrm{s}}^2 - v_{\mathrm{e}1}^2}{2c^2} - \sum ^4_{i=1}q^{(i)}_1 + \Lambda f_1 + \delta f_1 \end{aligned}$$

(13)

$$\begin{aligned} \frac{\Delta \phi _{\mathrm{e}2\mathrm{s}}}{c^2} \equiv \frac{\phi _{\mathrm{s}}-\phi _{\mathrm{e}2}}{c^2}& = \frac{\Delta f_2}{f_0} - \frac{v_{\mathrm{s}}^2 - v_{\mathrm{e}2}^2}{2c^2} -\sum ^4_{i=1}q^{(i)}_2 + \Lambda f_2 + \delta f_2 \end{aligned}$$

(14)

where the subscripts 1 and 2 denote, respectively, the values related to stations $P_1$ and $P_2$. Combining Eqs. (13) and (14), we obtain the equation which contains the GP difference between the two ground stations:

$$\begin{aligned} \frac{\phi _{\mathrm{e}21}}{c^2} &\equiv \frac{\phi _{\mathrm{e}1}-\phi _{\mathrm{e}2}}{c^2}= \frac{\Delta f_2 - \Delta f_1}{f_0} - \frac{v_{\mathrm{e}1}^2 - v_{\mathrm{e}2}^2}{2c^2} \\&\quad -\left( \sum ^4_{i=1}q^{(i)}_2-\sum ^4_{i=1}q^{(i)}_1\right) + (\Lambda f_2 - \Lambda f_1) + \delta f_{12} \end{aligned}$$

(15)

where $\phi _{\mathrm{e}21}=\phi _{\mathrm{e}1}-\phi _{\mathrm{e}2}$ is the Newtonian GP difference between the two ground stations $P_1$ and $P_2$. The error term $\delta f_{12}$ is the sum of the errors $\delta f_{1}$ and $\delta f_{2}$. From Eq. (15) we can see that since there is a pair of SFST links in determining the GP difference between two ground sites, the error magnitude from most error sources would be larger compared to one SFST link. However, the error sources from satellite, such as the error caused by velocity and position, can be significantly reduced because of their partial cancelations as shown by Eq. (15). They cannot be completely canceled out because although we intend to establish a pair of SFST links to one satellite simultaneously, in reality it is not quite possible to link two different ground stations at exactly the same time, because (1) even if the clocks located at two stations have been a prior synchronized, there may still exist time difference and (2) the signals propagation paths between the satellite and the stations are different (see Fig. 3). If the satellite receives the two signals from ground stations $P_1$ and $P_2$ at slightly different instants $t_1$ and $t_2$, the velocities of the satellite (${{\mathbf {v}}}_{\mathrm{s}}$ and ${{\mathbf {v}}}_{\mathrm{s}}^{\prime }$ at times $t_1$ and $t_2$) as shown in Eqs. (13) and (14) will be different. Thus the error term $\delta f_{12}$ in Eq. (15) contains a new error source which comes from asynchronism:

$$\begin{aligned} \delta f_{12} = \delta f_{\mathrm{cor}12} + \delta f_{\mathrm{sys}12} + \delta f_{\mathrm{asy}} \end{aligned}$$

(16)

where $\delta f_{\mathrm{asy}}$ is the asynchronism error, and the meanings of the other terms are the same as described in Eq. (7). In order to estimate the magnitude of $\delta f_{\mathrm{asy}}$, suppose the time interval between the received two signals from two ground stations is $\Delta t$, and in the time duration $\Delta t$ the satellite’s velocity changed from ${{\mathbf {v}}}_{\mathrm{s}}$ to ${{\mathbf {v}}}_{\mathrm{s}}^{\prime }$. Then, we have:

$$\begin{aligned} {{{\mathbf {v}}}_{\mathrm{s}}^{\prime }} = {{{\mathbf {v}}}_{\mathrm{s}}} + 0.5{{\mathbf {a}}_{\mathrm{s}}}\cdot {\Delta t}^2 \end{aligned}$$

(17)

where ${\mathbf {a}}_{\mathrm{s}}$ is centrifugal acceleration vector of the satellite. Substituting the ${{\mathbf {v}}}_{\mathrm{s}}$ to ${{\mathbf {v}}}_{\mathrm{s}}^{\prime }$ in Eq. (14), and then combining it to Eqs. (13) and (15), we can obtain the expression of $\delta f_{\mathrm{asy}}$:

$$\begin{aligned} \delta f_{\mathrm{asy}} < \frac{0.25a_{\mathrm{s}}^2\cdot \Delta t^4 + {{\mathbf {v}}}_{\mathrm{s}} \cdot {\mathbf {a}}_{\mathrm{s}} \cdot \Delta t^2}{c^2} + O(c^{-3}) \end{aligned}$$

(18)

where $O(c^{-3})$ are small amounts of (and higher than) $c^{-3}$ terms. With Eq. (18) we can estimate the influence of $\delta f_{\mathrm{asy}}$. For example, suppose the satellite-receiving time difference $\Delta t = 1$ ms, and the satellite is a GPS satellite whose centrifugal acceleration $|{\mathbf {a}}_{\mathrm{s}}|$ is about 0.558 m/s$^{2}$, velocity $|{{\mathbf {v}}}_{\mathrm{s}}|$ is about 3000 m/s (Zhang et al. 2006). Notice that ${{\mathbf {v}}}_{\mathrm{s}}$ and ${\mathbf {a}}_{\mathrm{s}}$ are almost orthogonal, then the calculated value of $\delta f_{\mathrm{asy}}$ is below $10^{-19}$, which is negligible in our case.

In order to guarantee that the satellite-receiving time difference $\Delta t\,<\,1$ ms or even smaller, we can employ the following two techniques. (1) According to the orbit of satellite, we can preestimate the distances between the satellite and the two ground stations (e.g., the distances of $P_1S$ and $P_2S$ in Fig. 3). Then, we can determine the suitable time for emitting frequency signals from the two ground stations. (2) When the satellite receives the frequency signals from the two ground stations, it can send a feedback signal that contains the information of the current satellite-receiving time difference. Once the two ground stations receive the feedback signals, they can adjust the signals’ emitting times correspondingly. We note that, as mentioned in Sect. 2, exact time synchronization is not necessary, because we make frequency transfer.

4 Simulation Experiments

Sections 2 and 3 provide the approaches of determining the GP difference in two cases: GP difference determination between one satellite and one ground station (Shen et al. 2016) and GP difference determination between two ground stations. In practical applications, these approaches can be used flexibly. For example, to improve the accuracy of the results, we can determine the GP at a certain ground station or the GP difference between two ground stations via several satellites (spacecrafts) which are equipped with high-precision clock systems. For the purpose of potential applications of the SFST approach in GP measurements in the future, in this section we conducted several simulation experiments as examples.

4.1 The Error Models of Various Error Sources

The reliability of a simulation experiment depends on whether the simulating case is close to the real case. In our experiments, we use GPS satellites whose orbit data are obtained from IGS product Web site (www.igs.org/products), and two ground stations located in China whose coordinates are also given. The key problem is the simulation of various error effects as described in Sect. 2 and 3. Because although the magnitude of each error source has been estimated, it is difficult to predict the value of each error in continuous experiments. In order to solve the problem, we adopt three kinds of error models in accordance with the different natures of the error sources.

First, the state of high-precision atomic clocks should be properly simulated. Galleani et al. (2003) have developed a mathematical model for clock error which can be expressed as:

$$\begin{aligned} X_1(t)& = x(0) + y(0)t + a\frac{t^2}{2} + \sigma _1\phi _1(t) + \sigma _2\int _0^t\phi _2(s)\mathrm{d}s \nonumber \\ X_2(t)& = c_2 + at + \sigma _2\phi _2(t) \end{aligned}$$

(19)

where $t\ge 0$ represents time, $X_1$ represents the phase deviation, $X_2$ represents the frequency deviation, x(0) and y(0) are initial conditions of $X_1$ and ${\dot{X}}_1$, respectively, $\phi _1(t)$ and $\phi _2(t)$ are Wiener processes (Brownian motion) defined by $\mathrm{d}W(t)=\xi (t)\mathrm{d}t$, where $\xi (t)$ is a white Gaussian noise with zero mean. $\sigma _1$ and $\sigma _2$ are constants that represent the diffusion coefficients of the two noises, a is a drift term, $c_2$ is the initial condition of $X_2(t)$. According to Eq. (19), a series of simulated clock data with errors embedded can be generated.

Second, we consider the error model of the satellite orbit. In the local satellite frame, the position errors of a satellite have three scalar components $\Delta x$, $\Delta y$ and $\Delta y$. In this local coordinate system (x, y, z), x points to the normal axis of the orbit plane, y points to the tangential axis, z points to the radial axis. According to Hill model, the velocity errors of a satellite satisfy the following equations (Colombo 1986):

$$\begin{aligned} \Delta {\dot{x}}(t)& = - \Omega \Delta x_0 \sin \Omega t + \Delta {\dot{x}}_0 \cos \Omega t \nonumber \\ \Delta {\dot{y}}(t)& = - 2\Delta {\dot{z}}_0 \sin \Omega t + \left( 4\Delta {\dot{y}}_0 + 6\Omega \Delta z_0 \right) \cos \Omega t - \left( 3\Delta {\dot{y}}_0 + 6\Omega \Delta z_0 \right) \nonumber \\ \Delta {\dot{z}}(t)& = \Delta {\dot{z}}_0 \cos \Omega t + \left( 2\Delta {\dot{y}}_0 + 3\Omega \Delta z_0 \right) \sin \Omega t \end{aligned}$$

(20)

where $t\ge 0$ represents time, the subscript “0” denote the initial condition, $\Omega$ is the orbital angular frequency. It should be noted that Eq. (20) holds in a rotating coordinate system (rotates about the x axis). For a non-rotating frame whose axes coincide with the moving ones at time t, there holds the following transformation:

$$\begin{aligned} \Delta {\dot{x}}(t)& = \Delta {\dot{x}}(t)_{\mathrm{NR}} \nonumber \\ \Delta {\dot{y}}(t)& = \Delta {\dot{y}}(t)_{\mathrm{NR}} - \Omega \Delta z(t) \nonumber \\ \Delta {\dot{z}}(t)& = \Delta {\dot{z}}(t)_{\mathrm{NR}} + \Omega \Delta y(t) \end{aligned}$$

(21)

where the subscript “NR” means “non-rotating”.

Finally, for other error sources (see Table 2), currently there are no mature mathematical models to simulate. Thus we adopted a general function of Wiener process to represent each of the other error sources:

$$\begin{aligned} X(t) = a + b\cdot \phi (t) + c\cdot \int _0^t \xi (s)\mathrm{d}s \end{aligned}$$

(22)

where X(t) is the error value at time t, $\phi (t)$ and $\xi (t)$ are both standard white Gaussian noises, a, b and c are constant coefficients. Clearly Eq. (22) is a simplified model and cannot perfectly simulate the values of various error sources. However, taking into consideration the large number of error sources and the relatively small amount of their magnitudes, this simplification is acceptable in our simulation experiments that aims at testing the precision of the SFST approach.

In summary, for each error source we have assigned an error model: Eq. (19) for the clock errors, Eq. (20) for satellite position and velocity errors and Eq. (22) for each of other error sources. The error models are independent of each other, and the coefficient values in a certain equation are determined in accordance with the error magnitude of the relevant error source. Then, the final error model of our simulation experiment is a combination of all these error models. In the following subsections, we conduct four types of simulation experiments and provide the corresponding results. The parameters and relevant error magnitudes used in our simulation experiments are listed in Tables 2 and 3.

Table 3 Relevant parameters used in simulation experiments

4.2 Determination of the GP at a Ground Station Via a Satellite

One of the presently most accurate Earth gravitational model [e.g., EGM2008 (Pavlis et al. 2012a)] provides an accuracy about 10–20 cm (equivalent to 2 m$^2$/s$^{2}$ in potential) at ground and may achieve at least 0.1 m$^2$/s$^{2}$ level at the (GNSS-type) target satellite altitude, which is around 20,000 km above the geoid. Hence, here in this study we just assume that the gravitational potential at the orbit of a GNSS-type or communication-type satellite is given at the accuracy level of 0.1 m$^2$/s$^{2}$, which is equivalent to 1 cm in height in the domain near the Earth’s surface.

We choose a ground station in Wuhan, China, whose geodetic coordinate is $114.32^{\circ }$E, $30.52^{\circ }$W, 50 m. The observation time period is 2.5 h, from 10:30 a.m. to 1:00 p.m., January 31, 2016. In this relatively short time duration, we choose a nearest GPS satellite (PG27) for our simulation experiments. At the ground station, the angle between the observation sight and zenith is $50.16^{\circ }$ at start point. It first decreases then increases to $30.77^{\circ }$ at end point during our experiments, as schematically shown by Fig. 4. Here we use simulation experiments via one satellite links to determine the GP at the Wuhan ground station based on the SFST and analyze the accuracy of the results. The simulation experiment method and relevant results are described as follows.

Fig. 4

Experiments are conducted at the time duration when satellite moves from position S to position ${S}^{\prime }$ (from 10:30 a.m. to 1:00 p.m., January 31, 2016). Ground station moves from P to $P^{\prime }$

First, the orbit information of the GPS satellite PG27 is obtained from IGS product Web site (www.igs.org/products). The precise ephemeris contains position and clock information, the time interval between two data set being 15 min. However, our simulation experiments are conducted every 10 s, hence the required data set was obtained by interpolation. The orbit data and ground position are regarded as true value, and we use EGM2008 model (Pavlis et al. 2012a) to calculate the GP values at ground station and satellite orbit at different times. These GP values are also regarded as true values. Other parameters such as the velocities of ground station and satellite can be calculated. The electron density (ionosphere influence) can also be obtained from IGS, and the atmosphere condition (troposphere influence) can be obtained from Earth Global Reference Atmospheric Model (Leslie and Justus 2011). Then, the ionosphere and troposphere residual corrections $\Lambda f_{\mathrm{ion}}$ and $\Lambda f_{\mathrm{tro}}$ can be calculated from Eqs. (8) and (9). The obtaining of other correction terms ($\Lambda f_{\mathrm{tide}}$, $\Lambda f_{\mathrm{position}}$ and $\Lambda f_{\mathrm{celes}}$) has been illustrated in Sect. 2. These correction terms are all regarded as true values. Thus according to Eq. (4), we can calculate the true value of the output frequency $\Delta f/f_0$.

The next step is adding noises. In Sect. 4.1, we have discussed the three kinds of error models, and the noises are generated according to these models and then added to the relevant true values. Consequently, we get a new set of “observations” which are used to estimate the value of interest. Then, we use Eq. (4) to calculate the GP difference $\phi _{{Si}} - \phi _{{Ei}}$ at time $t_i$ and denote it as $D_i$. Taking equal weight of each observation (at time $t_i$), we have

$$\begin{aligned} {\bar{\phi }}_E = \frac{\sum \nolimits _{i = 1}^n {\left( \phi _{Si} - D_i \right) }}{n}, \quad \sigma _{\phi _E} = \sqrt{\frac{\sum \nolimits _{i=1}^n {\left( \phi _{Si} - D_i - {\bar{\phi }}_E\right) ^2}}{n}} \end{aligned}$$

(23)

where ${\bar{\phi }}_{E}$ and $\sigma _{\phi _{E}}$ are the mean value of the estimated GPs at the ground station and the corresponding standard deviation (SD), respectively, $\phi _{{Si}}$ is the GP (with noises added) at the satellite orbit at time $t_i$, n is the total number of the “observations”.

By comparing the estimated value and the true value at time $t_i$ ($i=1,2,\ldots ,n$), we may verify the reliability of our proposed SFST approach (Shen et al. 2016). The results are shown in Fig. 5a. We can see that most of the absolute offset values are below the order of 1 m$^{2}$/s$^{2}$. There are 900 observations in total, and the mean value of the differences between the estimated GP and the true value at ground station at time $t_i$ $(i=1, 2,\ldots ,n)$ is −0.383 m$^{2}$/s$^{2}$, and the corresponding standard deviation (SD) is 0.385 m$^{2}$/s$^{2}$. The results in details are listed in Table 4 as Case 1.

Table 4 Setup and results of simulation experiments of four cases

Fig. 5

Gravitational potential (GP) of the ground station in Wuhan determined by: a the GPS satellite PG27, b 5 GPS satellites in combination (PG27, PG26, PG23, PG16 and PG08). And the gravitational potential (GP) differences between the two ground stations in Wuhan and Nanjing, determined by: c the GPS satellite PG27, d 5 GPS satellites in combination (PG27, PG26, PG23, PG16 and PG08). Experiment time period is from 10:30 to 13:00, January 31, 2016. We have an “observation” every 10 s and compare the true value with the estimated value. There are 900 comparisons for each satellite and the offsets between true values and estimated values are drawn as time series

4.3 Determination of GP at a Ground Station Via Observing Several Satellites

If more GNSS-type satellites equipped with high-precise clocks are available, the results of determining the GP at the ground station could be improved.

The setup of our second simulation experiment is similar to the first one as described in Sect. 4.2. The experiment date, time duration and location of ground station remain unchanged. The difference is that we use 5 GPS satellites (PG08, PG16, PG23, PG26 and PG27) to establish SFST links to the ground station at Wuhan. They are the most nearest GPS satellites in the experimental time period, and all of the angles among the observation sights and zeniths do not exceed $65^{\circ }$. For each set of links between one satellite and the ground station, the procedures are as same as described in Sect. 4.2, and we obtain one estimate of the GP at the ground station via every satellite. Taking different weights based on the separated accuracies, we obtain the weighted results, expressed as

$$\begin{aligned} {\bar{\phi }}_E = \frac{\sum \nolimits _{j = 1}^5 \left( \phi _j\cdot {\bar{\phi }}_{Ej}\right) }{\sum \nolimits _{j = 1}^5 \phi _j}, \quad \sigma _{\phi _E} = \sqrt{\frac{\sum \nolimits _{j=1}^5 {{\left( \phi _j\cdot \sigma _{\phi _{Ej}}\right) }^2}}{\sum \nolimits _{j=1}^5 \phi _j^2}} \end{aligned}$$

(24)

where j denotes the jth satellite, $\phi _j$ denotes the weight of the results based on the satellite j, ${\bar{\phi }}_{E_j}$ and $\sigma _{\phi _{E_j}}$ denote the estimated GP at ground station and the corresponding accuracy based on the jth satellite. Then ${\bar{\phi }}_E$ and $\sigma _{\phi _E}$ are the final results by combining the measurements of 5 satellites. Without obvious difference, here we just take equal weight, then the results are shown in Table 4 as case 2, and Fig. 5b shows the offset between the estimated values and the true values at time $t_i$. We can see that the result is better than (a), because some of the error sources can be reduced by multiple measurements. Here there are 900 observations for each satellite, and 4500 observations in total. The mean value of the differences is −0.272 m$^{2}$/s$^{2}$, and the standard deviation (SD) is 0.216 m$^{2}$/s$^{2}$.

4.4 Determination of GP Difference Between Two Ground Stations

If a satellite is connected with two ground stations via the SFST links simultaneously, as shown in Fig. 6, the GP difference between these two ground stations can be measured according to the results of the two groups of the SFST links. In this case, although the error sources from the satellite, such as the error of velocity and position can be significantly reduced, there exist new error sources stemming from the satellite-receiving simultaneity problem (see Sect. 3). Suppose two ground stations A and B are linked to a same satellite as link $L_A$ and $L_B$. The measurement times of $L_A$ are $t_{Ai}$ ($i=1,2,\ldots ,N$), and the measurement times of $L_B$ are $t_{Bi}$ ($i=1,2,\ldots ,N$). These measurement times are recorded by the clock (time-keeping system) equipped on the satellite; thus, they share the same time standard. We use the error models of Eq. (22) to simulate the asynchronism error, and the magnitude of $\delta f_{\mathrm{asy}}$ is below $10^{-19}$ as shown in Eq. (18).

Fig. 6

Experiments are conducted at the time duration when satellite moves from position S to position ${S^{\prime }}$ (from 10:30 a.m. to 1:00 p.m., January 31, 2016). Two ground stations move from $P_1$ and $P_2$ to $P_1^{\prime }$ and $P_2^{\prime }$

The setup of our experiments here is very similar to the first experiment as described in Sect. 4.2. The difference is that we added another ground station which is located in Nanjing (about 500 km from Wuhan station), with its geodetic coordinate being $118.78^{\circ }$E, $32.05^{\circ }$W, 30 m. The theoretical formulation is referred to Sect. 3. Since the two ground sties (Wuhan station A and Nanjing station B) can be bridged by one satellite or multiple satellites simultaneously, we made experiments corresponding to different cases.

For the mentioned two ground stations connected by one satellite, we have:

$$\begin{aligned} \Delta {\bar{\phi }}_{AB} = \frac{\sum \nolimits _{i = 1}^n {\left( \Delta {\hat{\phi }}_{ASi} - \Delta {\hat{\phi }}_{BSi} \right) }}{n} \end{aligned}$$

(25)

and

$$\begin{aligned} \sigma _{\Delta \phi _{AB}} = \sqrt{\frac{\sum \nolimits _{i=1}^n {\left( \Delta \hat{\phi }_{ASi} - \Delta \hat{\phi }_{BSi} - \Delta \bar{\phi }_{AB} \right) ^2}}{n}} \end{aligned}$$

(26)

where $\Delta {\bar{\phi}} _{AB}$ and $\sigma _{\Delta \phi _{AB}}$ are the mean value and standard deviation (SD) of the estimated GP difference between stations A and B, respectively, ${\Delta }\hat{\phi }_{AS_i}$ and ${\Delta \hat{\phi }_{BS_i}}$ are the estimated (calculated after adding noises) values of the GP differences between satellite and ground stations A and B at time $t_i$, respectively. n denotes the total number of the “observation pairs”, and there are 900 pairs of measurements in this case.

Similarly, for two ground stations connected by multiple satellites (here as an example we use 5 satellites), we have:

$$\begin{aligned} {\Delta \bar{\phi }_{AB}} = \frac{\sum \nolimits _{j = 1}^5 \left( \phi _j\cdot {\Delta \bar{\phi }_{(j)AB}}\right) }{\sum \nolimits _{j = 1}^5 \phi _j}, \quad \sigma _{\Delta \phi _{AB}} = \sqrt{\frac{\sum \nolimits _{j=1}^5 {{\left( \phi _j\cdot \sigma _{\Delta \phi _{(j)AB}}\right) }^2}}{\sum \nolimits _{j = 1}^5 \phi _j^2}} \end{aligned}$$

(27)

where j denotes the jth satellite, $\phi _j$ denotes the weight of the jth satellite. Since there are 5 satellites in total, $\Delta {\bar{\phi}} _{(j)AB}$ and $\sigma _{\Delta \phi _{(j)AB}}$ denote the determined results based on the jth satellite. Then ${\Delta {\bar{\phi }}_{AB}}$ and $\sigma _{\Delta \phi _{AB}}$ are the final results by combining the results based on 5 satellites. Taking equal weight, the results are shown in Table 4 (see cases 3 and 4), and Fig. 5c, d. In case 3 (see Fig. 5c) we use one satellite (PG27) to connect two ground stations; the mean value of the differences is 0.454 m$^{2}$/s$^{2}$, and standard deviation (SD) is 0.567 m$^{2}$/s$^{2}$. In case 4 (see Fig. 5d) we use 5 different satellites (PG08, PG16, PG23, PG26, PG27) to connect the two ground stations simultaneously (these 5 satellites are visible at the same time for the two ground stations in the experiment time duration), the mean value of the differences is 0.290 m$^{2}$/s$^{2}$, and standard deviation (SD) is 0.504 m$^{2}$/s$^{2}$. We can see that the results as shown by Fig. 5d are a little better than those as shown by Fig. 5c, but not very obviously. This is because some kinds of errors (such as the clock errors of ground stations, the tidal correction residual errors) cannot be obviously reduced by simply combining multiple satellites.

5 Conclusions

As a further improvement of the study of Shen et al. (2016), in this paper we formulated an approach for determining the GP of a ground station and the GP difference between two ground stations via one or more satellites and provided various simulating experiments addressing four different cases based on the tri-frequency combination (TrFC) technique. The precisions of determining the absolute GP of a ground station and the GP difference between two ground stations are estimated, reaching the level of 0.1 m$^{2}$/s$^{2}$, as long as the clocks’ stability and inaccuracy achieve the level of $1\times 10^{-18}$. Various influence factors (such as the tidal effects, the potentials of other celestial bodies, the frequency influences along the propagation path) have been considered and estimated. Their introduced errors do not exceed the error magnitude of clocks.

In Sect. 4 we have discussed the SFST approach of determining the GP of a ground station given one or several satellites’ GPs. Inversely, given GPs at ground stations with proper distribution, we can determine the GP at the orbit of a flying satellite. One potential application of the SFST approach is to determine the GP distribution along one or several GOCE-type satellites orbits and provide a potential distribution over a quasi-spherical surface constructed by the GOCE-type satellites. To complete this potential and prospective task, we need to first establish a ground datum station network to cover the whole orbits of the GOCE-type satellites via the SFST approach. The core idea is similar to determining the absolute GP at ground stations as described in this paper. Details are discussed in a separate study. Currently, due to the fact that optical clocks with stability of $10^{-18}$ level have been successfully realized (Poli et al. 2014; Bongs et al. 2015), with very quick development of time and frequency science, in the near-future portable optical clocks with stability of $10^{-18}$ level could be also realized. Consequently, the SFST approach will be prospective and potential for determining GP at any space position, not only providing an alternative approach for directly determining the GP, but also providing a way to realize the unification of the world height datum system.

Appendix: Formulation of Tri-Frequency Combination for Frequency Transfer at Accuracy Level of 10−18 in Free Space

For the purpose of an accuracy level of $10^{-18}$ in frequency transfer, the terms of $c^{-4}$ should be taken into consideration. In this Appendix we will derive a formula for the tri-link frequency transmission measurement [see Fig. 1; Eq. (3) in Sect. 2] in free space, achieving the accuracy requirement of $10^{-18}$ level.

Suppose the signal is emitted at point B (satellite), the one-way frequency shift received at A (ground station) is expressed as (Linet and Teyssandier 2002)

$$\begin{aligned} \frac{\Delta f _{BA}}{f} \equiv \frac{f_A}{f_B} -1 =\frac{1}{c^2} (\phi _A-\phi _B) +\frac{1}{2c^2}\big (v^2_A-v^2_B\big )+\sum ^4_{i=1} q^{(i)}_{(BA)} \end{aligned}$$

(28)

where $q^{(i)}$ is in the order of $1/c^i$ ($i=1,\,2,\,3,\,4$), expressed as

$$\begin{aligned} q^{(1)}_{BA}& = - \frac{1}{c}{\mathbf {N}}_{AB}\cdot ({\mathbf {v}}_A-{\mathbf {v}}_B) \end{aligned}$$

(29)

$$\begin{aligned} q^{(2)}_{BA}& = -\frac{1}{c^2} \left[ {\mathbf {N}}_{AB}\cdot ({\mathbf {v}}_A-{\mathbf {v}}_B) ({\mathbf {N}}_{AB}\cdot {\mathbf {v}}_B)\right] \end{aligned}$$

(30)

$$\begin{aligned} q^{(3)}_{BA}& = - \frac{1}{c^3} {\mathbf {N}}_{AB}\cdot ({\mathbf {v}}_A-{\mathbf {v}}_B) \left[ \left( \frac{1}{2}v^2_A-\frac{1}{2}v^2_B\right) +({\mathbf {N}}_{AB}\cdot {\mathbf {v}}_B)\right] \nonumber \\&\quad +\,\frac{1}{c^3} \left[ \left( \phi _A-\phi _B\right) {\mathbf {N}}_{AB}\cdot \left( {\mathbf {v}}_A-{\mathbf {v}}_B\right) - {\mathbf {l}}^{(2)}_A\cdot {\mathbf {v}}_A +{\mathbf {l}}^{(2)}_B\cdot {\mathbf {v}}_B\right] \end{aligned}$$

(31)

$$\begin{aligned} q^{(4)}_{BA}& = \frac{1}{c^4} \left\{ \frac{3}{8}v^4_A-\frac{1}{4}v^2_Av^2_B-\frac{1}{8}v^4_B - \left[ {\mathbf {N}}_{AB}\cdot ({\mathbf {v}}_A-{\mathbf {v}}_B)\right] \cdot ({\mathbf {N}}_{AB}\cdot {\mathbf {v}}_B)\left[ \left( \frac{1}{2}v^2_A-\frac{1}{2}v^2_B\right) +({\mathbf {N}}_{AB}\cdot {\mathbf {v}}_B)^2\right] \right\} \nonumber \\& \quad+\, \frac{1}{c^4} \left\{ (\gamma +1)(\phi _Av^2_A -\phi _Bv^2_B )+ \frac{1}{2}(\phi _A-\phi _B) \left[ (\phi _A-\phi _B) + 2(1-\beta )(\phi _A+\phi _B)\right. \right. \nonumber \\&\left. \left. \quad+\, v^2_A-v^2_B -2 {\mathbf {N}}_{AB}\cdot ({\mathbf {v}}_A-{\mathbf {v}}_B)\cdot ({\mathbf {N}}_{AB}\cdot {\mathbf {v}}_B )\right] + {\mathbf {N}}_{AB}\cdot \left[ ({\mathbf {l}}^{(2)}_B\cdot {\mathbf {v}}_B) ({\mathbf {v}}_A-2{\mathbf {v}}_B)+ ({\mathbf {l}}^{(2)}_A\cdot {\mathbf {v}}_A) ({\mathbf {v}}_B)\right] \right. \nonumber \\&\left. \quad+\, \left[ {\mathbf {l}}^{(3)}_A -2(\gamma +1) \mathbf {\zeta }_A\right] \cdot {\mathbf {v}}_A - \left[ {\mathbf {l}}^{(3)}_B -2(\gamma +1) \mathbf {\zeta }_B\right] \cdot {\mathbf {v}}_B \right\} +\frac{1}{c^2} (\varPsi _A-\varPsi _B) \end{aligned}$$

(32)

where

$$\begin{aligned} {\mathbf {N}}_{AB} = \frac{\mathbf {x}_B-\mathbf {x}_A}{R_{AB}} \end{aligned}$$

(33)

$\phi$ and $\varPsi$ are, respectively, the first and second Newtonian potentials, defined as

$$\begin{aligned} \phi& = G \int \frac{\rho (\mathbf {x}^\prime )}{|\mathbf {x}-\mathbf {x}^\prime |} d^3 \mathbf {x}^\prime \end{aligned}$$

(34)

$$\begin{aligned} \varPsi& = \frac{G}{c^2} \int \frac{\rho ^*(\mathbf {x}^\prime )}{|\mathbf {x}-\mathbf {x}^\prime |}\left[ \left( \gamma +\frac{1}{2}\right) v^2 + (1-2\beta )\phi + \varPi +3\gamma \frac{p}{\rho ^* }\right] d^3 \mathbf {x}^\prime \end{aligned}$$

(35)

where $\rho$ is the rest mass density, $\varPi$ is the specific energy density (ratio of internal energy density to rest mass density), p is the pressure, $\rho ^*$ is the conserved mass density, given by

$$\begin{aligned} \rho ^*=\rho \left[ 1+\frac{1}{c^2}\left( \frac{1}{2}v^2+3\gamma \phi \right) \right] \end{aligned}$$

(36)

and $\mathbf {\zeta }$ is vector potential, defined as

$$\begin{aligned} \mathbf {\zeta }= G \int \frac{\rho ^*(\mathbf {x^\prime })\mathbf {\omega }\times \mathbf {x}^\prime }{|\mathbf {x}-\mathbf {x}^\prime |}d^3 \mathbf {x}^\prime \end{aligned}$$

(37)

where $\mathbf {\omega }$ is the Earth’s temporal rotation angular velocity. Here in this study, for our purpose we just take $\mathbf {\omega }$ as a constant vector in what follows.

Accurate to our requirement, ${\mathbf {l}}_A (x_A, x_B)$ and ${\mathbf {l}}_B (x_A, x_B)$ are expressed as (Linet and Teyssandier 2002)

$$\begin{aligned} {\mathbf {l}}_A (x_A, x_B)& = - {\mathbf {N}}_{AB} + \left[ {\mathbf {l}}_M(x_A, x_B)+ {\mathbf {l}}_{J_2}(x_A, x_B)\right] + \left\{ {\mathbf {l}}_S(x_A, x_B) + {\mathbf {l}}_{v_r}(x_A, x_B) \right\} \end{aligned}$$

(38)

$$\begin{aligned} {\mathbf {l}}_B (x_A, x_B)& = - {\mathbf {N}}_{AB} - \left[ {\mathbf {l}}_M(x_B, x_A)+ {\mathbf {l}}_{J_2}(x_B, x_A)\right] + \left\{ {\mathbf {l}}_S(x_B, x_A) + {\mathbf {l}}_{v_r}(x_B, x_A) \right\} \end{aligned}$$

(39)

where

$$\begin{aligned} {\mathbf {l}}_M(x_A, x_B)& = -(\gamma +1)\frac{2GM}{c^2}\frac{(r_A+r_B){\mathbf {N}}_{AB}+R_{AB}{\mathbf {n}}_A}{(r_A+r_B)^2-R^2_{AB}}\end{aligned}$$

(40)

$$\begin{aligned} {\mathbf {l}}_{J_2}(x_A, x_B)& = (\gamma +1)\frac{GMJ_2}{c^2}\frac{(r_A+r_B)r^2_{\mathrm{e}}}{[(r_A+r_B)^2-R^2_{AB}]^2} \left\{ {\mathbf {N}}_{AB}\left[ 2(\mathbf {k}\cdot {\mathbf {n}}_A +\mathbf {k}\cdot {\mathbf {n}}_B)^2 \frac{(r_A+r_B)^2+3R^2_{AB}}{(r_A+r_B)^2-R^2_{AB}} \right. \right. \nonumber \\&\left. \left. \quad -\,\left( \frac{1-(\mathbf {k}\cdot {\mathbf {n}}_A)^2}{r_A}+\frac{1-(\mathbf {k}\cdot {\mathbf {n}}_B)^2}{r_B}\right) \frac{(r_A+r_B)^2+R^2_{AB}}{r_A+r_B}\right] +2 {\mathbf {n}}_A \frac{R_{AB}}{r_A+r_B}\left[ (\mathbf {k}\cdot {\mathbf {n}}_A +\mathbf {k}\cdot {\mathbf {n}}_B)^2\frac{3(r_A+r_B)^2+R^2_{AB}}{(r_A+r_B)^2-R^2_{AB}} \right. \right. \nonumber \\&\quad \left. \left. -\,\frac{1}{2}\big (1-3(\mathbf {k}\cdot {\mathbf {n}}_A)^2\big )\frac{(3r_A+r_B)(r_A+r_B)-R^2_{AB}}{r^2_A} +(r_A+r_B)\left(\frac{2(\mathbf {k}\cdot {\mathbf {n}}_A)(\mathbf {k}\cdot {\mathbf {n}}_B)}{r_A}-\frac{1-(\mathbf {k}\cdot {\mathbf {n}}_B)^2}{r_B}\right) \right] \right. \nonumber \\&\quad \left. -\,4 \mathbf {k}\frac{R_{AB}}{r_A} \left[ (\mathbf {k}\cdot {\mathbf {n}}_A)\frac{(3r_A+r_B)(r_A+r_B)-R^2_{AB}}{2r_A(r_A+r_B)} +(\mathbf {k}\cdot {\mathbf {n}}_B)\right] \right\} \end{aligned}$$

(41)

$$\begin{aligned} {\mathbf {l}}_S(x_A, x_B)& = \left( \gamma +1+\frac{1}{4}\alpha _1\right) \frac{2GC\omega }{c^3}\frac{r_A+r_B}{r_A\left[ (r_A+r_B)^2-R^2_{AB}\right] } \bigg \{ \mathbf {k}\times {\mathbf {n}}_B \nonumber \\& \quad+\,\frac{2r_Ar_B \mathbf {k}\cdot ({\mathbf {n}}_A \times {\mathbf {n}}_B) }{(r_A+r_B)^2 -R^2_{AB}} \left[ \frac{(3r_A+r_B)(r_A+r_B)-R^2_{AB}}{2r_A(r_A+r_B)} {\mathbf {n}}_A +{\mathbf {N}}_B\right] \bigg \} \end{aligned}$$

(42)

$$\begin{aligned} {\mathbf {l}}_{v_r}(x_A, x_B)& = \alpha _1\frac{GM}{c^3} \left[ \frac{{\mathbf {v}}_r-({\mathbf {v}}_r\cdot {\mathbf {N}}_{AB}) {\mathbf {N}}_{AB}}{2R_{AB}} {\text {ln}} \frac{r_A+r_B+R_{AB}}{r_A+r_B-R_{AB}} \right. \nonumber \\& \quad\left. +\,({\mathbf {v}}_r\cdot {\mathbf {N}}_{AB}) \frac{(r_A+r_B){\mathbf {N}}_{AB}+R_{AB}{\mathbf {n}}_A}{(r_A+r_B)^2 -R^2_{AB}}\right] \end{aligned}$$

(43)

where C is the Earth’s principal moment of inertia around $z-$axis, r denotes the Euclidean norm of the vector $\mathbf {x}$, $\mathbf {R}_{AB}=\mathbf {x}_B-\mathbf {x}_A$, $R_{AB}=|\mathbf {R}_{AB}|$, $\mathbf {k}$ is the unit vector along positive z axis, ${\mathbf {n}}=\mathbf {x}/r$ is the unit vector along $\mathbf {x}$ direction (example, ${\mathbf {n}}_A=\mathbf {x}_A/r_A$). In Eqs. (31) and (32)

$$\begin{aligned} {\mathbf {l}}^{(2)}_A/c^2& = {\mathbf {l}}_M (x_A, x_B)+{\mathbf {l}}_{J_2} (x_A, x_B),\quad {\mathbf {l}}^{(3)}_A/c^3 = {\mathbf {l}}_S (x_A, x_B)+{\mathbf {l}}_{v_r}(x_A, x_B), \nonumber \\ {\mathbf {l}}^{(2)}_B/c^2& = - {\mathbf {l}}_M (x_B, x_A) -{\mathbf {l}}_{J_2} (x_B, x_A),\quad {\mathbf {l}}^{(3)}_B/c^3 = {\mathbf {l}}_S (x_B, x_A)+{\mathbf {l}}_{v_r}(x_B, x_A), \end{aligned}$$

(44)

where the relevant vectors are given by expressions (40)–(43).

Suppose the emitting frequency at ground station A at time $t_1$ is $f _0$, the receiving frequency at satellite at time $t_2$ is $f ^\prime _0$, similar to Eq. (28), we have

$$\begin{aligned} \frac{f ^\prime _0-f _0}{f _0} = \frac{1}{c^2} \left( \phi _{B(2)}-\phi _{A(1)}\right) +\frac{1}{2c^2}\left( v^2_{B(2)}-v^2_{A(1)}\right) +\sum ^4_{i=1} q^{(i)}_{(12)} \end{aligned}$$

(45)

where the subscript “(i)” corresponds to the instants at $t_i (i=1,\,2,\,3)$, the subscript “(ij)” corresponds to the instants at $t_i$ and $t_j$ $(i = 1,\,2; j = i + 1)$ (for instance, $q^{(i)}_{(12)}$ denotes in fact $q^{(i)}_{A(1)B(2)}$). This received signal at satellite (at time $t_2$) is immediately transponded at time $t_2$ to the ground station, and the receiving frequency $f ^{\prime \prime }$ at ground station at time $t_3$ is, based on Eq. (28), expressed as

$$\begin{aligned} \frac{f ^{\prime \prime } _0-f ^\prime _0}{f _0} = \frac{1}{c^2} \big (\phi _{A(3)}-\phi _{B(2)}\big ) +\frac{1}{2c^2}\big (v^2_{A(3)}-v^2_{B(2)}\big )+ \sum ^4_{i=1} q^{(i)}_{(23)} \end{aligned}$$

(46)

Simultaneously at time $t_2$ the satellite emits a new signal with frequency $f _{\mathrm{s}}$, and the receiving frequency at ground is, based on equation (28), expressed as

$$\begin{aligned} \frac{f ^\prime _{\mathrm{s}}-f _{\mathrm{s}}}{f _0} = \frac{1}{c^2} \big (\phi _{A(3)}-\phi _{B(2)}\big ) +\frac{1}{2c^2}\big (v^2_{A(3)}-v^2_{B(2)}\big )+ \sum ^4_{i=1} q^{(i)}_{(23)} \end{aligned}$$

(47)

Concerning Eq. (47), it is worthy to notice that we omitted the difference between transponding the received signal at time $t_2$ and emitting a new signal at time $t^\prime _2$ at satellite, and that of the receiving signals at time $t_3$ and time $t^\prime _3$ at ground, namely we assume that $t^\prime _2=t_2$, $t^\prime _3=t_3$. This is related to time synchronization, which is not so serious in our frequency transfer scheme, because we compare the frequency, not the time elapsed.

Now we formulate a tri-frequency combination based on the emitting frequencies $f _0$ and $f_{\mathrm{s}}$ at ground and satellite, respectively, and the observed receiving frequencies $f^{\prime \prime }_0$ and $f^\prime _{\mathrm{s}}$. We examine the following frequency shift “output”

$$\begin{aligned} \frac{\Delta f}{f_0}\equiv \frac{\Delta f_{\mathrm{se}}}{f_0} -\frac{f^{\prime \prime }_0-f_0}{2f_0} \equiv \frac{f^\prime _{\mathrm{s}}-f_{\mathrm{s}}}{f_0} -\frac{f^{\prime \prime }_0-f_0}{2f_0} = \frac{f^\prime _{\mathrm{s}}-f_{\mathrm{s}}}{f_0} -\frac{(f^{\prime \prime }_0-f^\prime _0)+(f^\prime _0-f_0)}{2f_0} \end{aligned}$$

(48)

Substituting Eqs. (45)–(47) into Eq. (48), we obtain the following expression

$$\begin{aligned} \frac{\Delta f}{f_0}& = \frac{1}{2c^2} \left[ \phi _{A(3)} + \phi _{A(1)}- 2\phi _{B(2)}\right] + \frac{1}{2c^2} \left[ \frac{1}{2}v^2_{A(3)}+\frac{1}{2}v^2_{A(1)}-v^2_{B(2)}\right] + \frac{1}{2} \left[ \sum ^4_{i=1} q^{(i)}_{(23)} - \sum ^4_{i=1} q^{(i)}_{(12)} \right] \nonumber \\& = \frac{1}{c^2} \left[ \big (\phi _A - \phi _B \big ) + \frac{1}{2} \big ( v^2_A-v^2_B\big )\right] +\sum ^4_{i=1} q^{(i)} \end{aligned}$$

(49)

where $q^{(i)}$ is defined as

$$\begin{aligned} q^{(i)}= \frac{1}{2}\left[ q^{(i)}_{(23)}-q^{(i)}_{(12)}\right] \end{aligned}$$

(50)

where $q^{(i)}_{(12)}$ and $q^{(i)}_{(23)}$ (namely $q^{(i)}_{(A1B2)}$ and $q^{(i)}_{(B2A3)}$) are determined by Eqs. (29)–(32). In the last step of derivation in Eq. (49), we assume that $\phi _{A(3)}=\phi _{A(1)}$ and $v^2_{A(3)}=v^2_{A(1)}$, due to the fact that in a very short period (<0.5 s in general cases for our present purpose), we may consider the GP and speed of the ground station hold invariant. However, if necessary at any time, we may just apply the expression following the first equal sign “$=$” of Eq. (49).

We note that all relevant quantities are related to $t_1, t_2$, or $t_3$. For instance, ${\mathbf {l}}^{(i)}_A$ and ${\mathbf {l}}^{(i)}_B$ are related to $t_1, t_2$, or $t_3$. Based on Eq. (49) one can determine the GP difference between A and B.

The value of $q^{(i)}$ can be calculated by the given velocities of A and B and proper model values of the Newtonian potentials $\phi$, $\varPsi$ and the vector potential $\mathbf {\zeta }$. We use the GP model EGM2008 (Pavlis et al. 2012b) to calculate the model value of $\phi$, denoted as $\phi ^{\mathrm{EGM}08}$, which has at least the accuracy levels of tens of centimeters (one meter height is equivalent to 10 m$^2$/s$^2$ potential) at ground station and several centimeters at the satellite altitude. $\phi ^{\mathrm{EGM}08}$ is a harmonic expansion expression of the Earth’s external GP complete to degree/order 2159. Hence, in Eqs. (31) and (32) we use the following model value

$$\begin{aligned} \phi ^M(\mathbf {x})=\phi ^{\mathrm{EGM}08}(\mathbf {x}) \end{aligned}$$

(51)

In practice, for the purpose of model value used here, it is accurate enough to use the terms up to degree/order 20 of the EGM2008 model.

Since the second Newtonian potential $\varPsi$ itself is in the order of $1/c^2$ [see Eq. (35)], we see from equation (32) that $\varPsi$ plays only a role of the order $1/c^4$. Hence, accurate to the level of $c^{-4}$, in Eq. (32) we may take the following model value

$$\begin{aligned} \varPsi ^M = \frac{G}{c^2} \left[ \left( \gamma +\frac{1}{2}\right) \phi ^{\mathrm{EGM08}}v^2 +(1-2\beta )(\phi ^{\mathrm{EGM}08})^2 \right] \end{aligned}$$

(52)

Concerning the vector potential $\mathbf {\zeta }$, from expressions (37) and (32) we see that it also plays only a role of the order $1/c^4$. Hence, to achieve the accuracy level of $1\times 10^{-18}$, in Eq. (32) we may take the following model value

$$\begin{aligned} \mathbf {\zeta }^M= \frac{GC}{2r^3}\mathbf {\omega }\times \mathbf {x} \end{aligned}$$

(53)

In our simulation experiment, we just set $\beta =1/2, \gamma =1, \alpha _1=0$.