Formulation to test gravitational redshift based on the tri-frequency combination of ACES frequency links

Atomic Clock Ensemble in Space (ACES) is an ESA mission mainly designed to test gravitational redshift with high-performance atomic clocks in space and on the ground. A crucial part of this experiment lies in its two-way Microwave Link (MWL), which uses the uplink of carrier frequency 13.475 GHz (Ku band) and downlinks of carrier frequencies 14.70333 GHz (Ku band) and 2248 MHz (S band) to transfer time and frequency. The formulation based on the time comparison has been studied for over a decade. However, there are advantages of using frequency comparison instead of time comparison to test gravitational redshift. Hence, we develop a tri-frequency combination (TFC) method based on the measurements of the frequency shifts of three independent MWLs between ACES and a ground station. The potential scientific object requires stabilities of atomic clocks at least $3\times10^{-16}$/day, so we must consider various effects, including the Doppler effect, second-order Doppler effect, atmospheric frequency shift, tidal effects, refraction caused by the atmosphere, and Shapiro effect, with accuracy levels of tens of centimeters. The ACES payload will be launched as previously planned in the middle of 2021, and the formulation proposed in this study will enable testing gravitational redshift at an accuracy level of at least $2\times10^{-6}$, which is more than one order higher than the present accuracy level of $7\times10^{-5}$.


Introduction
General relativity theory [1] concludes three classic predictions: Mercury precession, light deflection and gravitational redshift. The first two have been confirmed by [1] and a group led by [2], but the gravitational redshift was not tested until 1960.
The first direct experimental verifications of gravitational redshift are the series of Pound-Rebka-Snider experiments during 1960-1965 [3,4], who observed the shift using a Mssbauer emitter and absorber at the Jefferson Physical Laboratory tower at Harvard University. Later, there is an around-the-world experiment. Four cesium beam clocks were used to fly around the world on commercial jet flights during several days in October 1971, and they flew in opposite directions while recording the time differences [5]. Additionally, other types of experiments measure the shift of spectral lines in the Sun's gravitational field since 1960 [6]. Typically, a Galileo solar redshift experiment tested the gravitational redshift to 1% accuracy [7]. The most famous test was obtained by the Gravity Probe A (GPA) mission in June 1976, which launched a hydrogen maser onboard a rocket to a height of 10 000 km [8]. During its flight, frequency comparisons were conducted between the maser on the rocket and a corresponding maser on the ground. The consistency of the relativistic frequency shift with the prediction was 7 × 10 −5 [9]. Until now, the most precise indirect tests were performed by eccentric Galileo satellites. The tests are based on the satellites GSAT-0201 and GSAT-0202 of the European Global Navigation Satellite System (GNSS) Galileo, which were accidentally delivered on elliptic instead of circular orbits. Two research teams simultaneously published their results with (0.19 ± 2.48) × 10 −5 [10] and (−0.9 ± 1.4) × 10 −5 [11], respectively.
The Atomic Clock Ensemble in Space (ACES) experiment [12,13,14,15], which was installed onboard the International Space Station (ISS), is an ESA-CNES mission mainly planned to test the gravitational redshift. Equipped with atomic clocks of fractional frequency instability and inaccuracy of (1 − 3) × 10 −16 , it aims to test the gravitational redshift at a level of 2 × 10 −6 [13,14], which is one and a half orders higher than the GPA experiment.
The main onboard instruments are an active hydrogen maser (SHM) and a cold cesium atoms (PHARAO). The PHARAO clock reaches a fractional frequency stability of 1.1 × 10 −13 √ τ , where τ is the integration time in seconds, and an accuracy of a few parts in 10 16 [14]. Meanwhile, SHM demonstrates a fractional frequency instability of 1.5 × 10 −15 after 10 000 s of integration time. Combining the short-term stability of the H-maser with the long-term stability and accuracy of the cesium clock, two clocks will generate an on-board time scale [13,14].
ACES enables frequency/time comparisons between ISS and ground stations by using two independent time & frequency transfer links (Microwave Links (MWL) and European Laser Timing (ELT) optical link) to test general relativity and develop applications in geodesy (relativistic geodesy) and time & frequency metrology [13,14]. These science objectives are closely related to the MWL performance [15], and its performance plays a key role in this study. MWL uses the uplink of carrier frequency 13.475 GHz (Ku band) and downlinks of carrier frequencies 14.70333 GHz (Ku band) and 2248 MHz (S band) to transfer time and frequency. MWL will perform with time deviation better than 0.3 ps at 300 s, 7 ps at 1 day, and 23 ps at 10 days of integration time [16]. These performances, which surpass those of existing techniques (TWSTFT and GPS) by 1-2 orders of magnitude, will enable comparisons (common view and uncommon view) of ground clocks with 10 −17 frequency resolution after a few days of integration [16].
Concerning the ACES mission, some studies have addressed the test of gravitational redshift based on time comparison [12,17,15], but there are almost no publications related to the frequency comparison. Compared with time comparison, frequency comparison has the following advantages: (1) it can weaken the effect of the phase ambiguity because the frequency measurement is irrelevant with ranging and is a consequence of counting during a short time; (2) it can determine the instant gravitational potential, while for time comparison, we must accumulate data to solve the time changing rate to deduce the gravitational redshift value. However, the accuracy of measuring the instant frequency is largely constrained, which implies that we must also accumulate observations to obtain results with higher accuracy.
In our study, we proposed a new formulation, referred to as tri-frequency combination (TFC) to obtain the gravitational potential difference by combining three frequency observations. For the one-way frequency transfer model with a precision requirement of 10 −16 , we adopt a formulation accurate to c −3 order in free space with medium, which was proposed by [18]. For our theoretical contributions, we extended the model of [18] from free space (vacuum) to real space with media (see section 2 and Appendix A) and formulated the approach to eliminate the Doppler frequency shift (the term Doppler effect or Doppler frequency shift mentioned in this paper refers to the firstorder Doppler effect) considering the time offset among three links (see section 3 and Appendix B). Our final TFC model can successfully eliminate all types of shifts to the order of 10 −16 . To verify our model and analyze the demanded magnitude of parameters, we designed simulation experiments considering the real orbit, reliable clocks noises, real atmosphere and real gravity (see section 5).

One-way frequency transfer between ISS and ground station
For MWLs, the ACES mission uses two different antennas: one Ku-band antenna for uplink and downlink and one S-band antenna for only downward signals. It uses the uplink of carrier frequency 13.475 GHz (Ku band, and the frequency shift will be broadcast to the ground station afterwards) and downlinks of carrier frequency 14.70333 GHz (Ku band) and 2248 MHz (S band) [12]. These three frequencies are denoted by f 1 , f 2 , and f 3 throughout this study. The goal of testing accuracy is 2 × 10 −6 ; thus, we need a frequency transfer model to the level of 1 × 10 −16 , which requires a relativistic model to the order of c −3 .
First, we consider a downlink from satellite A to ground B. The frequency transfer ratio f A /f B between proper frequencies f A and f B is determined by the clocks on the satellite (A) and the ground (B). In practice, this is achieved using the transmission of photons from A to B and the following formula where ν A and ν B are the proper frequencies of the photon at A and B. In a general relativistic framework, the proper frequency shift of the photon from A to B is expressed by [18] The first factor on the right-hand side is the sum of the gravitational redshift and transverse Doppler frequency shift, and U E is Newtonian potential of the Earth in the frame of Earth-Centered Earth-Fixed (ECEF). We denote radial vectors are the coordinate velocities. To the required order of 1/c 3 , the last factor in equation (2) is obtained from [18] with R AB = r B − r A , R AB = | R AB |, and N AB = R AB /R AB . The last terms in equations (3) and (4) are caused by curved geometry in the general relativistic framework and referred to as Shapiro effect. In this formulation, we approximate the Earth as a spherically symmetric body, since the J2 term does not exceed the magnitude of 4 × 10 −17 [18]. Using the simplified notation Figure 1. Principle of ray refraction through the atmosphere for ISS. r is distance from the Earth center, β is the angle between the tangent of the electromagnetic wave and the normal of the layer, θ is the angle between AB line direction and layer normal, and γ is the complementary angle of θ.
we have where A Shap is Shapiro effects given by equations (3) and (4). Models (5)-(6) only hold in vacuum. In a real space with medium, electromagnetic waves experience a change in direction of propagation or refractive bending when they are transmitted through the atmosphere, which is divided into the troposphere (0-60 km) and ionosphere (60-2000 km). Because of the refraction phenomenon, the direction of the refracted ray at the space station slightly differs from the unrefracted line-of-sight direction [19], as figure 1 shows. This phenomenon has significant applications in the GPS/MET (Global Positioning System/Meteorology) experiment [20] and will cause a slight change of Doppler effect in equation (6). Here, by definingĀ dop as Doppler frequency considering the atmosphere, we have The refractive index in the ionosphere is relevant with the carrier frequency, but the refractive index in troposphere is nearly irrelevant with it. Therefore, for all links f 1 , f 2 and f 3 , supposing that they are simultaneously emitted, the bending effects of the troposphere are approximately identical, which makes it easy to wipe out the tropospheric part. For the ionospheric part, to the order of f −2 , we have [21,20] where n E is the electron density per cubic meter; high orders such as f −3 are neglected because they are at least two magnitudes smaller than the order of f −2 [22], which we will later analyze.
In the phase form, Doppler frequency shift is given by phase path P of the radio wave [21,23,24,25]: where the velocity of the source induces a change in λ, and the velocity of the observer changes dP/dt. Thus, in vacuum, the first-order Doppler frequency shift can be derived in terms of A dop , as expressed by equation (5).
More specifically, we have [26] where T is a unit vector of the wave normal, n is the refractive index, and α is the angle between the wave normal and the ray direction. We suppose that the atmosphere is an isotropic medium, the refractive index is nearly independent of ray directions, and α = 0 [27]. From this equation, the atmospheric influence can be explained by two reasons: the refractive index varies with time [23,24], and the wave path varies because the observer moves [28].
Due to the velocity of the source, wavelength λ is expressed as ConsideringĀ dop defined in equation (7), with (10) and (11), we havē With the height of ISS of approximately 400 km [12], the ACES-ground links lie in the middle layer of the ionosphere. The integral term of the refractive index in equation (12) can be expressed as the sum of the ionospheric and tropospheric parts, and if we suppose α = 0, we have [25] B A ∂n ∂t where n e is the electron density along the trajectory, M 1 = 77.6 × 10 −6 p/T , and M 2 = 0.373ε/T 2 with temperature T , total pressure p and partial pressure of water vapor ε along the trajectory. For this expansion, equation (13) can be rewritten as where (referring to Appendix A) where v Ax , v Ay , v Bx , and v By are components of velocities v A and v B of space station A and ground site B projected in the refraction plane in figure 1, and δ A and δ B are deviated angles from the line of sight. More details are referred to Appendix A.
Focusing on the effect of the variation of the refractive index and wave path, we can obtain the expression ofĀ dop ; then, we separate the refraction part (not time-varying parts), ionospheric and tropospheric frequency shift (time-varying parts), and Shapiro effect. Thus, the one-way frequency transfer model should be written as where where δf refr is the bending effect on Doppler frequency shift, which is caused by refraction, δf ion and δf trop are atmospheric effects caused by the time-varying refractive index. For the ACES links, estimates show that the magnitude of δf dop is approximately 10 −5 ; δf rel is approximately 10 −10 [18]; δf ion , δf ion and δf trop will be estimated in section 5.

Tri-frequency combination
Although all of these links are independent and can be synchronized afterwards by data processing, the synchronization error may cause severe residual errors. In this study, we list coordinate time t 1 ∼ t 6 to identify six events and use a combination of three frequency links of ACES to test the gravitational redshift, as figure 2 shows. Later, we will analyze the required synchronization precision to test gravitational redshift at the 2 × 10 −6 level. In our formulation, we use a nonrotating geocentric space-time coordinate system. At coordinate time t 1 , Ku-band signal f 1 is emitted and received by the space station at coordinate time t 2 . Meanwhile, two signals f 2 and f 3 are emitted from the space station at coordinate time t 3 and t 4 , respectively, and received by the ground station at coordinate time t 5 and t 6 , respectively. If we define a coordinate time interval by T ij = t j − t i , T 23 and T 34 will theoretically be synchronized to zero, but in practice, there is a difference between them.
For ACES links f 1 = 13.475 GHz, f 2 = 14.70333 GHz and f 3 = 2248 MHz, the third link is of low frequency, which greatly suffers from ionospheric effects. If we suppose that T 34 is extremely small (< 1 µs), the only different error between link 2 and link 3 is the ionospheric error, because other shifts in link 2 and link 3 are close. We define f ′ 1 , f ′ 2 and f ′ 3 as the received frequencies corresponding to emitted frequencies f 1 , f 2 and f 3 . (5), (16) and (17), the Doppler part, relativistic parts (hereafter, they refer to transverse

ISS trajectory
Ground station trajectory MWL principle (modified after [29]). At time t 1 , the ground station emits signal f 1 to ISS, which is received at time t 2 . Meanwhile, two signals f 2 and f 3 are emitted to the ground station at time t 3 and t 4 (they are approximately t 2 ), which are received at time t 5 and t 6 .
Doppler effects and gravitational redshift) and tropospheric part are cancelled, and we obtain where δ ion A and δ ion B are the refractive angles relevant with the carrier frequency, and they are inversely proportional to the square of the carrier frequency (see Appendix A). All terms in the last middle bracket in equation (18) are inversely proportional to the square of the carrier frequency.
Referring to figure 2, the ground station has moved a certain distance (approximately 1 meter) from t 1 to t 5 ; thus, the first-order Doppler frequency shift cannot be completely cancelled from equation (16). By the same technique, we divide where Bi is the position of the ground station at time t i , and Ai is the position of ISS at time t i . Here, m 1 /m 2 is the solved ionospheric part where the terms in the second bracket are related to ionospheric effects and determined by equation (18), which can give rise to the following relation If the difference in space parameters of the space station and ground station at different time points is ignored, the Doppler effect in equation (19) can be directly eliminated. However, at different time points such as t 2 and t 3 , spatial parameters of ISS will be moderately different, which results in Doppler residuals in equation (19). According to Appendix B and accurate to the order of c −3 , we have where m 1 /m 2 is determined by equation (21); A Shap,link2 is a term of c −3 ; since the magnitudes of T 23 and T 15 are at most c −1 , the difference between A Shap,link1 and A Shap,link2 in equation (21) is at most c −4 , which is negligible, as is the difference between v 2 a2 /c 2 and v 2 a3 /c 2 .

Error sources
Generally, errors are divided into systematic errors and random errors. All of the aforementioned errors are systematic errors, including Doppler frequency shift, atmospheric frequency shift (including ionospheric, tropospheric and refractive frequency shift), relativistic frequency shift and Shapiro frequency shift. These frequency shifts can be eliminated using our TFC model, but residuals remain. These residuals and tidal effects will be discussed in section 4.1.
Random errors are caused by devices (e.g., atomic clocks, cables and emitters) and measurements (e.g., velocities and accelerations in equation (22) and (23)). Literatures [12,14] have shown the Allan Deviation performance of ACESs clocks (SHM and PHARAO), which shows that PHARAO has better long-term stability. However, these studies did not show the noise components of the clocks of ACES, and we can only simulate the clock data to approach their performance. For measurement noises, the accuracy for parameters r A , v A , a A and T 23 must be carefully controlled. Section 5 will discuss the parameter demands using our simulation data.

Residual errors
Although section 3.1 provided a practical calculation model to test gravitational redshift using the TFC method, there are residual errors, which are mainly reflected in two aspects: First, we have performed many approximations in the model derivation; Second, there are other types of errors in nature that we have not considered, such as tidal effects.
The model of the TFC method can be summarized by equations (21)- (23). In step one, frequency shifts of downlinks 2 and 3 are divided to obtain the ionospheric part; then, frequency shifts of uplink 1 and downlink 2 are divided. In step two, we substitute this result with the calculated residual Doppler effect and ionospheric part, so the gravitational potential difference can be calculated.

Doppler residual errors
For Doppler residuals in section 3.1, we only consider Doppler shift difference between link 1 and link 2 while ignoring that of link 2 and link 3. Since the time difference of link 2 and link 3 is less than 1 µs (T 34 < 1 µs), and they are both downlinks, Doppler shift difference between link 2 and link 3 is much smaller. Supposing that T 34 is 100 ns, we take similar notes as section 3.1 and Appendix B: A denotes time t 3 , B denotes time t 5 , B ′′ denotes time t 4 , and A ′′ denotes time t 6 . Based on the spatial relation, we have With a numerical calculation with equation (24), we obtain that the numerical difference between T 34 and T 56 is tens of nanosecond, which implies that both T 34 and T 56 are in the order of magnitude of c −2 and must be corrected. Through calculation, we neglected the intermediate process and obtained the relation of the Doppler shift difference between link 2 and link 3 With an approximate numerical calculation, and the largest part in equation (25) is cR AB T 34 , which will achieve 3 × 10 −14 . Hence, assuming that T 34 = 100 ns, the maximal Doppler residual errors caused by link 2 and link 3 are 3 × 10 −14 . Since in equation (21) is approximately −0.0053, this error will affect gravitational redshift equation (22) by a magnitude of 5 × 10 −16 .
In Appendix B, we show the relevant approximations. We neglect the terms of c −4 , whose value are less than 10 −17 .

Ionospheric residual errors
In the expansion of the ionosphere refractive index, the terms of f −3 cannot be eliminated by our TFC method and have not been considered in the former sections. According to the study of [30], the f −3 terms of GNSS signals (L band) is approximately one in a hundred of f −2 terms. By reckoning of our simulations, final errors caused by the f −3 terms will be approximately 1 × 10 −15 . Because of its variability and randomness, these errors will be effectively weakened to a considerably small value by averaging with a mass of data. The simulations in section 5 will demonstrate that this effect can be lower than 10 −16 .

Relativistic residual errors
The calculation of relativistic effects will also cause residuals. There are differences among relativistic effects of links 1, 2 and 3; however, in our analysis, we consider it a constant value. The relativistic differences between link 2 and link 3 are much smaller than those between link 1 and link 2 (the time difference is much smaller), so we will only consider the differences between link 1 and link 2. Relativistic effects include gravitational redshift effect and transverse Doppler effect. Regardless of ground displacements such as the solid Earth tide, the gravitational redshift effect only varies with ISS, but the transverse Doppler effect varies with both ISS and ground station. Based on simulated data, the gravitational potential of ISS varies by 100 m 2 /s 2 per second. However, T 23 is at most 1 µs, so the gravitational redshift difference between link 1 and link 2 is approximately 10 −20 . The transverse Doppler frequency shift is proportional to the square of velocity, and its differential expression is Assuming that T 23 = 1 µs and T 15 = 2 ms, the maximal transverse Doppler frequency shift errors caused by ISS and the ground station are 7.6×10 −19 and 3×10 −19 , respectively. Thus, the maximal transverse Doppler frequency shift error is 1 × 10 −18 .

Tidal effects
With regard to ground displacements, the relativistic effects have other residuals. According to IERS convention 2010 [22], displacements of reference points are divided into three categories: (1) Tidal motions (mostly near diurnal and semidiurnal frequencies) and other accurately modeled displacements of reference markers (mostly at longer periods); (2) other displacements of reference markers including nontidal motions associated with the changing environmental loads; (3) displacements that affect the internal reference points in the observing instruments. We are interested in the first two types. Tidal motions include solid, ocean and polar tides. The solid tide has the largest magnitude, which will be tens of centimeters [31,22]. Other displacements include nontidal mass redistributions in the atmosphere, oceans, sea-level variations, etc., but their amplitude is much smaller (at most several centimeters) [32]. Crustal deformation and plate motions also contribute to the parameters ( r B , v B , a B ), which can be predicted by the tectonic model NUVEL-1 NNR of [33]. However, the numerical values of the kinematic vectors shows that the rates are only approximately several centimeters per year or tens of micrometers per day, which is 4 orders of magnitude smaller than the solid Earth tide (∼ tens of centimeters per day). Therefore, we will take the solid Earth tide as the dominant contribution.
Based on the theory as proposed in IERS convention [32,34], we estimated the equatorial displacement caused by solid tides, and the maximum is 0.32 m. If a sine function is used to estimate the change in displacement with time, we assume that all tides are diurnal tides, and the effect of the tides on the velocity is 4.65 × 10 −5 m/s. According to those estimates, the magnitude of the position and velocity vectors affects the residual Doppler terms in (22) and (23) by at most 6.76 × 10 −17 , which is much less than the Doppler residuals analyzed in section 3.1.
The solid Earth tide causes displacements and gravitational potential disturbance. The maximum tidal potentials caused by the moon and the sun are 4.41 m 2 /s 2 and 1.60 m 2 /s 2 , respectively [34]. With regard to Love numbers h, l and k, the gravitational potential tide will be 3.74 m 2 /s 2 , which is equivalent to a gravitational redshift of 4.2 × 10 −17 , which is mainly caused by the semidiurnal tide.
Based on the above discussion in section 4.1, the residual amount of each frequency shift is shown in table 1. Among them, the Doppler frequency shift and ionosphererelated frequency shift are relatively large and can be manually eliminated if possible. The ionospheric part is difficult to manually eliminate, as demonstrated in further research (section 5).

Test of gravitational redshift
In a static gravitational field, suppose that two atomic clocks are located at different positions. Then, we compare their frequencies by a certain frequency transfer method. Based on general relativity, gravitational redshift ∆ν between clocks is proportional to their gravitational potential difference ∆U as ∆ν To test the gravitational redshift by a standard convention, parameter α is introduced via the following expression [35] where α vanishes when Einstein equivalence principle (EEP) is valid. In our study, we develop a similar equation using another parameter β where ∆U m is the measured gravitational potential difference by equations (22) and (23), and ∆U is the standard gravitational potential difference developed by the Earth gravity field model. There are testing errors in both ∆U m and ∆U, and the corresponding uncertainties should be calculated using the following equation where u ∆Um and u ∆U are the uncertainties of ∆U m and ∆U, respectively.

Simulation experiments and results
At present, there are no real data. To test our theory and formulations, we present simulation experiments. In these experiments, we use the data of the ISS real orbit, ionosphere, troposphere, calculated gravitational potential by the widely used gravity field model EGM2008 [36], solid Earth tide [34], and simulated clock data by a conventionally accepted stochastic noises model [37,38].

Simulation setup and experiments
In our simulations, we select the station Observatoire de Paris (OP) as the ground station with geographical parameters as shown in table 2. In section 3.2, we discussed that the magnitude of the tidal effect on gravitational redshift is approximately 10 −17 ; nonetheless, in addition to the accuracy of the ACES program, we still added the tidal effect to the simulation experiment. When tidal effects are neglected, the relevant parameters in ECEF related to the ground station OP are considered constant, as shown in table 2. Other parameter settings in our experiment are listed in table 3. Figure 3 shows the procedures of the simulation experiments. The observations in our simulation experiments are carrier frequency values of ACES, which are denoted as

as section 3 demonstrates) to calculate the gravitational potential difference between ISS and ground station and test the gravitational redshift.
To obtain the received frequency values (f ′ 1 , f ′ 2 , and f ′ 3 ), the original emitted frequency values and various types of frequency shifts are required: (1) Originally emitted frequency values of f 1 , f 2 , and f 3 , and clock noises (including devices noises); (2) frequency shifts including Doppler frequency shift, relativistic frequency shift (including second-order Doppler shift and gravitational redshift), atmospheric frequency shifts (including ionospheric and tropospheric parts), and tidal effects. To calculate these effects, we must have the position, velocity and acceleration information, which can be derived from the orbits of ISS and moving positions of the ground station. To select observable data, we set the condition that the observation elevation angle should be larger than 15 • .
First, we must solve the emitted frequency values, which are frequency series composed of the given frequencies f 1 , f 2 , f 3 and clock noises. Based on the stochastic noise nature of the clocks, there are five types of clock noises: Random Walk FM (frequency modulation), Flicker FM, White FM, Flicker PM (phase modulation) and White PM [37] with spectral densities of the types f −2 , f −1 , f 0 , f 1 and f 2 , respectively.
Before starting our experiments, we provide samples of the simulated clock frequency noise series in figure 4. We simulated five pure types of noises with similar magnitudes (figures 4(a-e)) and their sum (figure 4f). The data length is 20 000. Figure  4a shows a strong trend, figure 4b shows little trend and some periodic patterns, and figures 4(c-e) show irregular patterns. White and random walk noises are the simple types [38]. Flicker noises were calculated by the AR (autoregressive) model [39].
In simulations, we took White FM as the dominant part because Allan deviation performance of PHARAO is very similar to that of pure white FM, and we simulated 864 000 s of clock data with sample intervals of 1 s. The modified Allan deviation (MDEV, [37]) of our simulated data is shown in figure 5, and its performance is similar   with previous studies [12,14]. Figure 5 shows a potential of long-term stability of 10 −16 .
With these clock data, we succeeded in the first step: the emitted frequency values were generated.
In the experiment, to simulate the orbit, we use the daily orbit elements to calculate. The TEC (Total Electron Content) data are interpolated from the grid data, which are derived from CDDIS (The Crustal Dynamics Data Information System). Using the TEC data and time-varying rate, we can simulate the ionospheric frequency shift. The refractive frequency shift caused by the ionosphere is calculated by the refractive angle, which is integrated from the layered structure adopted model [20] of ionosphere. The troposphere is calculated by the zenith tropospheric delay (ZTD) of the dry and wet components, and the projection function adopted the Vienna Mapping Function (VMF1) mode. The refractive frequency shift caused by the ionosphere is not calculated because this part is independent of the carrier frequency and easy to eliminate. We adopted EGM2008 [36] for gravitational potential models. We considered the effect of the tidal effect on the ground gravitational potential because this is the main source of residual error in gravitational redshift. For this part, we first calculated the Earth tide generated potential from parameters of periodic tides [34] using the method of harmonic analysis. Then, we considered Love numbers and solved the gravitational potential tide.
Combining with the analysis in section 4, we did not consider the indirect tidal effect on the coordinates variations of the ground station because this indirect tide is far from other residual errors. ISS is flying in an orbit with an inclination of 51.6 • and a period of 5400 s, and its track of the subsatellite point is shown in figure 6. In this figure, the red part is where we can observe. Because ISS is flying in a low orbit with a large velocity, each pass over the ground station will last approximately 600 seconds at most [15]; for our elevation threshold of 15 • , we can only observe approximately 300 s per pass.

Results and accuracy level of the test
For our experiment, we had 29-day simulated observations, and ISS flew above OP 130 times in total. As shown in figure 7, we drew a subsatellite point trajectory with the station OP at the center. Because the inclination of ISS is 51.6 • , the maximal latitude of the subsatellite point is that value. After analyzing the orbits of ISS, we selected the time when we could observe with elevations greater than 15 • . Based on the position information of ISS and OP, downloaded atmospheric parameters, gravitational quantities and various models, we calculated all types of frequency shifts as shown in figure 3, whose magnitudes are shown in table 1. The result shows that Doppler frequency shift is the dominant part, while Shapiro frequency shift is the smallest part. Refraction effects cannot be neglected because they have larger magnitude than the ionospheric frequency shift and tropospheric frequency shift. In the experiment, although the residuals of some errors (especially the higher-order term of the ionosphere) are greater than 10 −16 , due to the randomness of the errors, simulation experiments show that after a long-term average, they can be less than 10 −16 .
Furthermore, we analyzed each frequency shift in a one-way transfer (link 3, figure  8) and the residual of various frequency shifts after applying the TFC method ( figure 9). For the time length, we only show one pass. Figure 8 is consistent with table 1, where the relativistic frequency shift appears unchanged because the gravitational potential and velocity norm of ISS slowly change with time. There appears to be singularities in figure 8, but those are not real. When ISS was exactly on top of the station, the velocity became vertical to the line of sight (LoS), which made Doppler and Shapiro effects extremely small and resemble a singularity.
The TFC method can largely eliminate various errors, but the analysis in section 4.1 shows that residual errors remain. According to table 1, we calculated the  Doppler residuals, ionospheric-related residuals (including ionospheric frequency shift and ionospheric refraction), transverse Doppler residuals, and tidal effects. The Doppler residual is caused by the Doppler difference between link 2 and link 3. The ionosphericrelated residual is caused by the higher-order ionospheric term. The transverse Doppler residual is caused by the change in velocity. Figure 9 shows that among the residuals of various frequency shifts, the largest component is the ionospheric residual, which can reach a maximum of 1 × 10 −15 and is often at the level of 10 −16 ; the second largest component is the Doppler residual, whose maximum is 1×10 −16 ; the magnitudes of other frequency shifts are much smaller than 10 −16 , so they are negligible. These two frequency shifts should be corrected in the TFC method model, but simulation experiments show that due to the variability of these two frequency shifts, they can be eliminated below 10 −16 after a long-term average, so they can be ignored. Thus, the TFC method model In addition to the discussed systematic residual errors, random errors will be caused by the parameters of the ground station and space station in the results. We will later analyze the impact of these errors. Figure 8 and figure 9 only show the frequency shift during a single flight of the space station, and figure 10 shows the data for a longer time period. We define the observations of one single pass of ISS as an epoch. Figure 10 . Using our TFC model as section 3 proposed, we can obtain the gravitational potential differences (blue line) and compared them with the real differences calculated by EGM2008 in figure 11(a). The results are obtained by only one single epoch, and the dotted line in the figure is the range of the error in the 1 − σ criterion. Figure 11(b) shows the gravitational potential differences after averaging epoch by epoch. The results show that the precision is not good for data sampled per second ( figure 11(a)), which is equivalent to 250 m in height, which is within our expectation. However, after the epoch-by-epoch averaging, we obtained a gravitational potential difference bias series. The precision of the gravitational potential differences has improved to approximately 218 m 2 /s 2 (approximately equivalent to an elevation of 22.2 m), as shown in figure 11(b). After averaging the entire data set, we obtain a gravitational potential difference bias of 4.0 ± 18.6 m 2 /s 2 . Compared with the real averaged gravitational potential difference of 3.8340 × 10 6 m 2 /s 2 , our testing level achieves (1.04 ± 4.85) × 10 −6 . Supposing that the gravitational potential difference calculated by EGM2008 has an uncertainty of 3 m 2 /s 2 , the testing uncertainty achieves √ 18.6 2 + 3 2 3.8340 × 10 6 ≈ 4.91 × 10 −6 ; therefore,  our testing level is (1.04 ± 4.91) × 10 −6 . Here, we only used 29 days of observations. With longer-period experiments, the results will be better.

Accuracy requirements of relevant parameters
For the calculation of our TFC model, we examine the requirements of the parameter precision (table 5). To obtain results with high precision, we need parameters to satisfy these demands. r A and r B are the positions of ISS and the ground station, respectively; v A and v B are their velocities; a A and a B are their accelerations; T 23 is the time offset of link 1 and link 2. All ground station parameters here can be accurately obtained based on the coordinates and Earth rotation model and do not need to be discussed in detail. Thus, we must only study the parameters related to the space station. These parameters are demanded for the calculation of equation (22) and (23); thus, we derived the total differential of these equations and analyzed the differential relationship between parameters ( r A , v A , a A and T 23 ) and gravitational potential difference. The measurement noises are stochastic, so after averaging, the noises will be largely weakened. Due to this principle, based on the MDEV curve of ACES clocks, we set a target precision of 10 −13 at a sampling rate of 1 second.
According to table 5, the requirements of orbit parameters and acceleration parameters are easy to satisfy, and T 23 is easy to control at that level. Only the precision requirement of the velocity is high and must be carefully solved. It is easy for ISS to satisfy these demands. In our parameter analysis, for each analysis of one parameter, we suppose that the other parameters are constant. In this case, we may clearly examine the effect of the error of each parameter on the results. In the real environment, there are stochastic errors in all parameters. Table 5 shows that only the velocity is the principal error source if all errors exist. Therefore, our parameter analysis scheme is feasible.

Conclusion
In this study, first we proposed a new one-way frequency transfer model based on the frequency transfer theory to order c −3 [18] and Doppler shift considering the refractive effect [26]. Our model is established in free space and has an accuracy of c −3 . Then, we proposed the TFC method to derive the gravitational redshift based on three frequency links of ACES. Unlike other combination methods of ACES [17,15], our model addresses the frequency comparison instead of time comparison. This combination method can largely eliminate Doppler frequency shift, atmospheric effects, etc. We also calculated the residual errors and test whether they could be eliminated. Furthermore, we designed simulation experiments to authenticate our model and analyzed the demanded parameters. The simulation results show that our testing level can achieve (1.04 ± 4.91) × 10 −6 . Equipped with ACES atomic clocks (SHM and PHARAO) with a frequency stability of (1 − 3) × 10 −16 , our study shows that the test of gravitational redshift can achieve an accuracy level of 2 × 10 −6 if longer data (say much longer data than 29 days) are obtained, which is one and a half orders higher than the accuracy level given by [9]. This study provides a new approach to test the gravitational redshift and benefits the brand-new field of relativistic geodesy.
ACES limits its testing precision at the 10 −6 level due to the limitation of the accuracy of the atomic clock it carries. To achieve higher scientific goals, more precise space atomic clocks and time-frequency comparison links are required. The Chinese Space Station (CSS) provides this opportunity. Since the CSS will carry optic-atomic clocks on the order of 10 −18 , which is one and a half orders of magnitude higher than the cold atomic cesium clocks of ACES. It is expected to achieve 10-cm to 5-cm level gravitational potential measurements and can test general relativity at a level of 10 −7 .
However, due to the higher accuracy of the clock on the CSS, there are differences in model and actual measurement. In terms of the model, it must consider the timefrequency comparison model under the relativistic framework of c −4 accuracy, and some of the originally negligible residual items in section 4.1 must be carefully considered. In addition, because the magnitude of the high-order ionospheric term is obvious and difficult to eliminate, one may consider a potential solution to generalize the TFC method to a four-frequency combination and remove the high-order ionospheric terms. If the hardware conditions cannot satisfy the requirements of the quad-band link, a more in-depth study of the ionospheric frequency shift is required. If the accuracy of the model must reach the level of 10 −18 , it is also necessary to consider the tidal effect.
In terms of actual measurement, due to the higher accuracy of gravitational potential, according to the model, there are higher-accuracy requirements in measuring the position and speed of the space station, which introduces further requirements on the positioning system of the space station and hardware facilities for time measurement. If the accuracy requirement is increased by one and a half orders of magnitude, the velocity of space station should achieve the centimeter accuracy level. In summary, the formulation of the TFC method proposed in this study is also feasible for the CSS plan.

Appendix A.
In our modeling, we focus on the down link from A to B, but since the optical path is reversible, we draw figure 1 as if it is the link from B to A. Assuming that our Earth is a sphere, and the atmosphere is composed of many thin spherical shells, in each of which the medium is homogeneous. Based on Bouguers law (Snells law in spherical surface) [40,30], the following equation holds where n is the refractive index, and i is the zenith angle at each layer. For electromagnetic waves that propagate in a medium with refractive index n, the differential bending angle dα that accrue on the ray path over a differential arc length ds is given by [28] dα = n −1 T × ∇n ds where T is its unit tangent vector defined by T = k/k, k(r) is the vector wavenumber of the ray at point r, and ∆n is the gradient vector of the refraction index. Here, differential bending angle dα is written in vector form and related to the value and direction.
The refraction index in the ionosphere is displayed by equation (8) n = 1 − 40.
3 n e f 2 (A. 3) We suppose that the refractive index only changes in the vertical direction. At the lower middle part of the ionosphere layer, electron density n e is at the peak [20]. If we integrate equation (A.2) from the peak height to the ISS height, the gradient of the refractive index can be expressed as a vertical downward vector form, and the integration is where e g is the unit vector whose direction is defined by T × ∇n. According to equation (A.3), ∂n/∂h is a positive value, and we use where R E is the Earths average radius. For the links of ACES, we consider the ionospheric part (discussed in section 2), and we let [20] n e (h) = where h max and N max correspond to the peak height and peak density, respectively; H is the free electron density scale height. Considering (A.4) and (A.5), we have the following for the ionosphere Substituting (A.7) into (A.5) and taking n = 1 for approximation, equation (A.5) gives where integration I is It is weekly correlated to carrier frequency f , and the tropospheric part can be determined in this manner. For the total bending angle, we have Based on the reversible principle of the optical path, for three signals of ACES, the only difference is carrier frequency f . Moreover, we have α ion ∼ 1/f 2 .
To calculate the effect of refraction on the Doppler effect, we must separately calculate the deflection angles at A and B. If we define δ A and δ B to be the deflection angle between the direction of the ray and the direct connection of A and B, which implies that As figure 1 shows, we have According to equations (A.1) and (A.12), we have For equations (A.10), (A.11) and (A.13), δ A and δ B can be solved by the expression of α From equation (A.14), n A and α are relevant with carrier frequency f (ISS is in the ionosphere layer, the term n A − 1 is proportional to f −2 , and the ground station is in the troposphere layer), and the term n A r A cos θ A − n B r B cos θ B is nearly irrelevant to carrier frequency f because the term n A − 1 is small compared to the entire term. Thus, equation (A.14) can be divided into two parts: the first part is irrelevant to carrier frequency f , and the second part is proportional to f −2 .
To calculate T · v in equation (14), we must project these two vectors to the same plane. We set point B to be the origin, direction OB to be the y-axis, and the tangent of the Earth surface to be the x-axis, and we obey the right-hand rule to define the z-axis. This coordinate system is defined as the Local Link Coordinates System (LLCS). We can transform vectors from ECEF to LLCS and select the x and y components in LLCS, which are the projected coordinates. In the xoy plane of LLCS, vectors T A and T B can be expressed as Considering ACES links f 1 and f 2 , we mark the positions of the ground station at t 1 as B ′ , space station at t 2 as A ′ , space station at t 3 as A, and ground station at t 5 as B (as figure 2 shows).
To express the Doppler frequency shift to the order of c −1 , it is easy to obtain the shift of two links as Considering that the time interval of this process is very short (known for the height of ISS, which is approximately T 12 = T 35 = 1 ms; according to [17], we need T 23 < 1 µs, two Doppler frequency shifts are numerically close.
To simplify the overall frequency shift, we must express all terms at time t 1 and t 2 by the terms at time t 3 and t 5 . For this study, we deduced the algorithm adopted by the appendix of [18]. Using this algorithm, we have According to the actual value, we roughly calculated the speed value of ISS, and v/c is approximately 10 −5 . Based on this magnitude, the magnitude of T 23 is approximately nanosecond; thus, they can be considered a c −2 term. However, T 15 is considered a c After squaring the vector and expanding the root of this result into secondary series, we obtain: Substitute T 15 with (B.5), we have With the velocity of the ground station and space station, in terms of velocity, acceleration and derivative of acceleration [18], we have For the approximation of ACES: N AB · v A /c ≤ 2.6 × 10 −5 , N AB · v B /c ≤ 1.6 × 10 −6 , R AB · a A /c 2 ≤ 1.7 × 10 −10 , and R AB · a B /c 2 ≤ 7 × 10 −13 [18]; to the order of 10 −16 , we have