Effect of Earth-Moon’s gravity on TianQin’s range acceleration noise. III.
An analytical model

Lei Jiao Xuefeng Zhang [email protected] MOE Key Laboratory of TianQin Mission, TianQin Research Center for Gravitational Physics & School of Physics and Astronomy, Frontiers Science Center for TianQin, Gravitational Wave Research Center of CNSA, Sun Yat-sen University (Zhuhai Campus), Zhuhai 519082, China.

Abstract

TianQin is a proposed space-based gravitational wave detector designed to operate in circular high Earth orbits. As a sequel to [Zhang et al. Phys. Rev. D 103, 062001 (2021)], this work provides an analytical model to account for the perturbing effect of the Earth’s gravity field on the range acceleration noise between two TianQin satellites. For such an “orbital noise,” the Earth’s contribution dominates above $5\times 10^{-5}$ Hz in the frequency spectrum, and the noise calibration and mitigation, if needed, can benefit from in-depth noise modeling. Our model derivation is based on Kaula’s theory of satellite gravimetry with Fourier-style decomposition, and uses circular reference orbits as an approximation. To validate the model, we compare the analytical and numerical results in two main scenarios. First, in the case of the Earth’s static gravity field, both noise spectra are shown to agree well with each other at various orbital inclinations and radii, confirming our previous numerical work while providing more insight. Second, the model is extended to incorporate the Earth’s time-variable gravity. Particularly relevant to TianQin, we augment the formulas to capture the disturbance from the Earth’s free oscillations triggered by earthquakes, of which the mode frequencies enter TianQin’s measurement band above 0.1 mHz. The analytical model may find applications in gravity environment monitoring and noise-reduction pipelines for TianQin.

^†^†preprint: APS/123-QED

I Introduction

TianQin plans to launch three identical satellites equipped with drag-free control of high precision, which orbit the Earth at an altitude of $\sim 10^{5}\mathrm{km}$ . They form a nearly equilateral triangle constellation facing the white dwarf binary RX J0806.3+1527 and detect gravitational waves (GWs) through measuring distance change between satellites using laser interferometry of picometer-level precision [1]. In addition to GWs, the gravity field of the Earth-Moon’s system also causes distance change between the satellites. The effects are two-fold. First, it leads to the constellation deviating from the nominal equilateral triangle, which must be reduced by orbit optimization and control [2, 3, 4]. Second, it may induce in-band range variations mixing with GW signals, and hence pose a potential noise source to TianQin [5].

In order to assess the impact of the Earth-Moon’s gravity, we have developed a program TQPOP (TianQin Quadruple Precision Orbit Propagator) in the previous work [5]. The simulator takes in optimized initial orbit parameters, and uses detailed force models to propagate the satellites’ pure-gravity orbits with low numerical noise that is not available from commonly used double precision. The resulting ephemerides were converted to the inter-satellite range acceleration $\ddot{\rho}$ , and we computed the amplitude spectral density (ASD) and made the comparison with the acceleration noise requirement of TianQin. According to the estimates, we expect the effect of the Earth-Moon’s gravity to be below the noise requirement in the sensitive frequency band of TianQin ( $10^{-4}-1$ Hz) and hence not constituting a showstopper to the mission [5]. Furthermore, we also explored how the Earth-Moon’s gravity disturbance varies with different orbital radii and orientations for TianQin [6]. A general trend is that the disturbance shifts to lower frequencies as the orbital radius increases, which may set the lower bound of the measurement band. Some useful guidelines were drawn for the orbit and constellation design for future geocentric missions like TianQin. All these previous works rely heavily on numerical calculations, and it would be desirable to also examine the relevant issues by analytical modeling so as to, e.g., gain more understanding.

In the line of analytical work, our team have studied the effect of the Sun’s and the Moon’s point masses on TianQin’s constellation stability [7]. By solving Lagrange’s perturbation equations, the explicit orbit solutions have been derived and linked to the arm-length and breathing-angle variations, which can account for the long-term lunisolar influence with sufficient accuracy. In the range-acceleration ASD, the Moon’s point mass dominates below $5\times 10^{-5}$ Hz [5]. Hence in this work, we will instead focus on modeling the disturbance of the Earth’s gravity, as it dominates above $5\times 10^{-5}$ Hz in the frequency spectrum, and is more relevant to TianQin.

Important lessons can be learned from low-low satellite-to-satellite tracking gravimetry missions represented by GRACE (Gravity Recovery And Climate Experiment) and GRACE Follow-On [8, 9]. The missions mainly consist of two identical satellites with separation of about $220$ km and flying in the same circular polar orbit at an altitude of $300-500$ km. The formation measures inter-satellite distance variation using microwave/laser ranging systems. From the viewpoints of experiment design and data processing, TianQin shares many common aspects with GRACE and GRACE Follow-On, and therefore (semi-)analytical work in space gravimetry can offer valuable reference to TianQin’s modeling.

For instance, Kaula’s linear perturbation approach [10, 11, 12, 13, 14] can be used for analytically calculating the Earth’s gravity field coefficients from satellite observables. An important procedure therein is by expressing the Earth’s gravity field experienced by a satellite as a function of the orbital elements [10, 11]. Also based on Kaula’s representation of geo-potential, an analytical method for the Earth’s gravity-field recovery was provided by [15] (see also [16, 17, 18, 19, 20, 21, 22, 23, 24, 25]). The method shares the same observable with the inter-satellite line-of-sight differential gravitational acceleration approach [26, 27, 28, 29, 30], and contains an analytical expression of the differential gravitational acceleration of two GRACE satellites which fly along a circular reference orbit encircling the Earth in uniform rotation.

Based on the above works, we will construct the analytical model for the influence of the Earth’s static gravity field and consider TianQin’s special case of a high altitude and long measurement baseline. Then we will, taking the Earth’s free oscillations as an example, extend the model to include the Earth’s time-variable gravity field. In this latter regard, the reference [31] has discussed the influence of the Earth’s free oscillations triggered by large earthquakes on the inter-satellite measurement of GRACE. We will follow a similar treatment in this work.

This paper is organized as follows. Section II shows the detailed derivation of the analytical model. Section III tests the model with numerical results obtained from TQPOP. Subsection III.1 presents the result for TianQin, and Subsection III.2 presents the results for other orbital inclinations and radii. Section IV extends the model to include time-variable gravity field from the Earth’s free oscillations. Section V carries out the corresponding numerical verification. At last, Section VI draws the conclusions with outlooks.

II Model Derivation: Static Gravity

In this section, we detail the derivation of an explicit expression for the range acceleration $\ddot{\rho}$ between two circularly orbiting TianQin satellites under the influence of the Earth’s static gravity field in uniform rotation. For readers’ convenience, the table containing the symbols and their meanings in this paper is shown in the Appendix A.

II.1 Basic mathematical setup

As shown in Fig. 1, we have two satellites SC1,2 moving along the same nominal circular orbit around the rotating Earth. The phase difference between the two is $\gamma=120^{\circ}$ , i.e., the angle subtended by the satellites with respect to the Earth.

Refer to caption — Figure 1: TianQin satellites in circular prograde orbits around the Earth in the International Celestial Reference System (ICRS) with a phase difference $\gamma=120^{\circ}$ . The orbit radius is $1\times 10^{5}$ km, and the inclination and longitude of the ascending node are $74.5^{\circ}$ and $211.6^{\circ}$ [2], respectively. The arrows on SC1,2 mark the flight directions.

In the International Terrestrial Reference System (ITRS), the Earth’s gravity field can be modeled by [32]

\begin{split}&V(r,\theta,\lambda)=\frac{GM}{R}\sum_{n=0}^{N}\sum_{m=0}^{n}\left(\frac{R}{r}\right)^{n+1}\\ &\times\left(\bar{C}_{nm}\cos(m\lambda)+\bar{S}_{nm}\sin(m\lambda)\right)\bar{P}_{nm}(\cos\theta),\end{split}

(1)

where $G$ is the universal gravitational constant, $M$ the Earth’s total mass, $R$ the Earth’s average radius, and $n$ , $m$ the degree and order of the spherical harmonic expansion, respectively. Moreover, $N$ denotes the truncation degree, and $r$ , $\theta$ and $\lambda$ represent the radius, co-latitude and longitude, respectively, in the spherical coordinate system. Following the convention in gravimetry, the geo-potential takes on positive values. The spherical harmonic coefficients $\bar{C}_{nm}$ and $\bar{S}_{nm}$ , and the associated Legendre functions $\bar{P}_{nm}(\cos\theta)$ are all normalized, with $\bar{C}_{1m}=0=\bar{S}_{1m}$ .

For modeling satellite motion, we switch to the International Celestial Reference System (ICRS). In order to depict the local gravity field experienced by the satellites, one may substitute the spherical coordinates with the orbital elements and express the geo-potential in terms of of the latter [10, 11]:

\begin{split}&V(a,i,e,\omega,M,\Omega,\Theta)=\frac{GM}{R}\sum_{n=0}^{N}\left(\frac{R}{a}\right)^{n+1}\\ &\times\sum_{m=0}^{n}\sum_{p=0}^{n}\bar{F}_{nmp}(i)\sum_{q=-\infty}^{\infty}G_{npq}(e)S_{nmpq}(\omega,M,\Omega,\Theta),\end{split}

(2)

where the arguments $\left\{a,i,e,\omega,M,\Omega,\Theta\right\}$ denote, respectively, the semi-major axis, the inclination, the eccentricity, the argument of perigee, the mean anomaly, the right ascension of the ascending node, and the Greenwich Apparent Sidereal Time (GAST). The explicit form of the function $S_{nmpq}$ is given by

S_{nmpq}(\omega,M,\Omega,\Theta)=\alpha_{nm}\cos\psi_{nmpq}+\beta_{nm}\sin\psi_{nmpq}.

(3)

Here the coefficients $\left\{\alpha_{nm},\beta_{nm}\right\}$ are rearrangements of $\left\{\bar{C}_{nm},\bar{S}_{nm}\right\}$ :

\begin{array}[]{c}\alpha_{nm}=\left\{\begin{array}[]{cc}\phantom{-}\bar{C}_{nm},&n-m=\mathrm{even},\\ -\bar{S}_{nm},&n-m=\mathrm{odd},\end{array}\right.\\ \\ \beta_{nm}=\left\{\begin{array}[]{cc}\bar{S}_{nm},&n-m=\mathrm{even},\\ \bar{C}_{nm},&n-m=\mathrm{odd},\end{array}\right.\end{array}

(4)

and the angular argument reads

\psi_{nmpq}=(n-2p)\omega+(n-2p+q)M+m(\Omega-\Theta).

(5)

For the function $G_{npq}(e)$ , only the case of $e=0$ is used. The expression for the inclination function $\bar{F}_{nmp}(i)$ will be discussed in the next subsection, with the normalization related to $\bar{P}_{nm}(\cos\theta)$ .

It is worth mentioning that in deriving Kaula’s representation Eq. (2), one has neglected the procession, nutation, and polar motion of the Earth’s rotation axis [11]. These effects are typically slow-varying [32]. For TianQin’s measurement band of $10^{-4}-1$ Hz (roughly 1s to 3 hours in time scales), we deem it reasonable to assume uniform Earth rotation in the derivation and will verify it with the numerical simulation in Sec. III.

II.2 Circular reference orbits

Now we introduce another important simplification by adopting circular reference orbits. The reasons are two-fold. First, the nominal orbit of TianQin is circular and has a large radius of $\sim 1\times 10^{5}$ km by design. After orbit optimization, the deviation from perfect circularity is primarily due to the Moon’s and the Sun’s gravitational perturbations, which increase the eccentricity only to $5\times 10^{-4}$ on average [2]. Second, the approximation also makes the further derivation much more tractable, which is already quite complicated on its own. More specifically, in calculating $\ddot{\rho}$ (see Eq. (17)), we will directly use the circular reference orbits and ignore the coupling between the Earth’s non-spherical and third-body gravitational acceleration and the satellite’s perturbed position from the circular orbit. The validity of this approximation is to be confirmed by numerical simulation in Sec. III.

Assuming $e=0$ , Eq. (2) can be simplified to [15]

\begin{split}&V=\frac{GM}{R}\sum_{n=0}^{N}\left(\frac{R}{a}\right)^{n+1}\sum_{m=0}^{n}\sum_{p=0}^{n}\bar{F}_{nmp}(i)\\ &\times\left(\alpha_{nm}\cos\psi_{nmp}+\beta_{nm}\sin\psi_{nmp}\right)\end{split}

(6)

with the angular argument

\psi_{nmp}=(n-2p)\omega^{o}+m\omega_{e},

(7)

and

\omega^{o}=\omega+M,\quad\omega_{e}=\Omega-\Theta.

(8)

Note that only the $q=0$ components of $G_{npq}(0)$ retain the value $1$ while the other vanish. Now in Eq. (6), the orbit elements $\{a,i,e,\Omega,\omega\}$ are of fixed values, and $M$ increases linearly with time.

In practice, it is more convenient to rearrange the summation indices and their order so that Eq. (6) can be rewritten as [15]

\begin{split}&V=\sum_{m=0}^{N}\sum_{k=-N}^{N}\sum_{n=n_{1}[2]}^{n_{2}}u_{n}(a)\bar{F}_{nm}^{k}(i)\\ &\times\left(\alpha_{nm}\cos\psi_{mk}+\beta_{nm}\sin\psi_{mk}\right),\end{split}

(9)

where one has

u_{n}(a)=\frac{GM}{R}\left(\frac{R}{a}\right)^{n+1},

(10)

and

\psi_{mk}=m\omega_{e}+k\omega^{o},

(11)

with $k=n-2p$ . The summation bounds $\{n_{1},n_{2}\}$ are related to $\{m,k,N\}$ by

\begin{array}[]{l}n_{1}=\left\{\begin{array}[]{ll}|k|,&|k|\geq m,\\ m,&|k|<m,\quad m-k=\mathrm{even},\\ m+1,&|k|<m,\quad m-k=\mathrm{odd},\end{array}\right.\\ \\ n_{2}=\left\{\begin{array}[]{ll}N,&N-k=\mathrm{even},\\ N-1,&N-k=\mathrm{odd},\end{array}\right.\end{array}

(12)

and the notation $[2]$ means that the step size of the summation is 2.

Introducing the auxiliary index

g=\left\{\begin{array}[]{ll}\frac{n-m}{2},&n-m=\mathrm{even},\\ \frac{n-m-1}{2},&n-m=\mathrm{odd},\end{array}\right.

(13)

one can express the inclination function $\bar{F}_{nm}^{k}(i)$ as

\begin{split}&\bar{F}_{nm}^{k}(i)=\sum_{s=0}^{m}\sum_{t=t_{1}}^{t_{2}}\sum_{c=c_{1}}^{c_{2}}\frac{(-1)^{c+g}}{2^{2n-2t}}\frac{(2n-2t)!}{t!(n-t)!(n-m-2t)!}\\ &\times C_{m}^{s}C_{n-m-2t+s}^{c}C_{m-s}^{\frac{1}{2}n-\frac{1}{2}k-t-c}(\sin i)^{n-m-2t}(\cos i)^{s},\end{split}

(14)

with

t_{1}=0,\quad t_{2}=\left\{\begin{array}[]{ll}-\frac{1}{2}k+\frac{1}{2}n,&k\geq n-2g,\\ g,&k\leq n-2g,\end{array}\right.

(15)

and

\begin{array}[]{c}c_{1}=\left\{\begin{array}[]{cc}0,&2s-k-2t\leq 2m-n,\\ s-\frac{1}{2}k-t+\frac{1}{2}n-m,&2s-k-2t\geq 2m-n,\end{array}\right.\\ \\ c_{2}=\left\{\begin{array}[]{cc}s-2t+n-m,&2s+k-2t\leq 2m-n,\\ -\frac{1}{2}k-t+\frac{1}{2}n,&2s+k-2t\geq 2m-n.\end{array}\right.\end{array}

(16)

Examples of explicit forms of $\bar{F}_{nm}^{k}(i)$ are shown in Appendix B.

II.3 Range acceleration

By differentiating the (instantaneous) range $\rho=|\mathbf{r}_{2}-\mathbf{r}_{1}|$ between two TianQin satellites, the range acceleration in the presence of the geo-potential $V$ reads

\ddot{\rho}=\left(\nabla V_{2}-\nabla V_{1}\right)\cdot\ \hat{\mathbf{e}}_{12}+\frac{1}{\rho}\left(\left|\dot{\mathbf{r}}_{2}-\dot{\mathbf{r}}_{1}\right|^{2}-\dot{\rho}^{2}\right).

(17)

The first term on the right-hand side,

\mathcal{DA}:=\left(\nabla V_{2}-\nabla V_{1}\right)\cdot\hat{\mathbf{e}}_{12},

(18)

represents the differential gravitational acceleration along the line of sight between the satellites, where the subscripts $\{1,2\}$ indicate SC1 and SC2, respectively, and $\hat{\mathbf{e}}_{12}$ is the unit vector from SC2 to SC1. The second term on the right-hand side of Eq. (17),

\mathcal{CA}:=\frac{1}{\rho}\left(\left|\dot{\mathbf{r}}_{2}-\dot{\mathbf{r}}_{1}\right|^{2}-\dot{\rho}^{2}\right),

(19)

represents the centrifugal acceleration with $\mathbf{r}_{1}$ , $\mathbf{r}_{2}$ the position vectors of SC1 and SC2, respectively. Based on the circular orbits and geometric relation, this term is given by

\mathcal{CA}=\frac{2GM\sin(\gamma/2)}{a^{2}},

(20)

which is a constant and depends on the phase difference $\gamma$ . This means that in evaluating the ASD of $\ddot{\rho}$ , we neglect the contribution from the centrifugal acceleration, of which the validity can be verified by numerical simulation for TianQin.

In terms of the orbit elements, the expression of $\mathcal{DA}$ can written as (see [15] for more details)

\mathcal{DA}=\frac{1}{a}\left(\frac{\partial V_{2}}{\partial\omega^{o}_{2}}-\frac{\partial V_{1}}{\partial\omega^{o}_{1}}\right)\cos\left(\frac{\gamma}{2}\right)+\left(\frac{\partial V_{2}}{\partial r_{2}}+\frac{\partial V_{1}}{\partial r_{1}}\right)\sin\left(\frac{\gamma}{2}\right).

(21)

Then we evaluate the partial derivatives to have

\begin{split}&\frac{\partial V}{\partial\omega^{o}}=\sum_{m=0}^{N}\sum_{k=-N}^{N}\sum_{n=n_{1}[2]}^{n_{2}}ku_{n}(a)\bar{F}^{k}_{nm}(i)\\ &\times\left(\beta_{nm}\cos\psi_{mk}-\alpha_{nm}\sin\psi_{mk}\right)\end{split}

(22)

and

\begin{split}&\frac{\partial V}{\partial r}=\sum_{m=0}^{N}\sum_{k=-N}^{N}\sum_{n=n_{1}[2]}^{n_{2}}-\frac{n+1}{a}u_{n}(a)\bar{F}^{k}_{nm}(i)\\ &\times\left(\alpha_{nm}\cos\psi_{mk}+\beta_{nm}\sin\psi_{mk}\right).\end{split}

(23)

By substitution into Eq. (21), we arrive at [15]

\begin{split}&\mathcal{DA}=\sum_{m=0}^{N}\sum_{k=-N}^{N}\sum_{n=n_{1}[2]}^{n_{2}}\frac{u_{n}(a)}{a}\bar{F}^{k}_{nm}(i)\times\\ &\left[\left(k\cos\left(\frac{\gamma}{2}\right)\beta_{nm}-(n+1)\sin\left(\frac{\gamma}{2}\right)\alpha_{nm}\right)\cos\psi_{mk}^{2}\right.\\ &\!\!\!-\left(k\cos\left(\frac{\gamma}{2}\right)\alpha_{nm}+(n+1)\sin\left(\frac{\gamma}{2}\right)\beta_{nm}\right)\sin\psi_{mk}^{2}\\ &\!\!\!-\left(k\cos\left(\frac{\gamma}{2}\right)\beta_{nm}+(n+1)\sin\left(\frac{\gamma}{2}\right)\alpha_{nm}\right)\cos\psi_{mk}^{1}\\ &\!\!\!\left.+\left(k\cos\left(\frac{\gamma}{2}\right)\alpha_{nm}-(n+1)\sin\left(\frac{\gamma}{2}\right)\beta_{nm}\right)\sin\psi_{mk}^{1}\right],\end{split}

(24)

where the superscript of $\psi_{mk}^{1,2}$ means SC1,2. Introducing $\omega^{om}$ as the argument of latitude of the midpoint of SC1,2 and

\left\{\begin{array}[]{l}\psi_{mk}^{1}=m\omega_{e}+k\left(\omega^{om}-\frac{\gamma}{2}\right),\\ \psi_{mk}^{2}=m\omega_{e}+k\left(\omega^{om}+\frac{\gamma}{2}\right),\end{array}\right.

(25)

we can rewrite $\mathcal{DA}$ as [15]

\begin{split}\mathcal{DA}=\sum_{m=0}^{N}\sum_{k=-N}^{N}&\left[E_{mk}^{c}\cos\left(m\omega_{e}+k\omega^{om}\right)\right.\\ &\!\!\!\left.+E_{mk}^{s}\sin\left(m\omega_{e}+k\omega^{om}\right)\right],\end{split}

(26)

with

	$\displaystyle E_{mk}^{c}=\sum_{n=n_{1}[2]}^{n_{2}}\Xi_{nmk}(a,i,\gamma)\alpha_{nm},$		(27)
	$\displaystyle E_{mk}^{s}=\sum_{n=n_{1}[2]}^{n_{2}}\Xi_{nmk}(a,i,\gamma)\beta_{nm},$		(28)

and

\begin{split}&\Xi_{nmk}(a,i,\gamma)=-2\frac{u_{n}(a)}{a}\bar{F}_{nm}^{k}(i)\\ &\times\left(k\sin\left(k\frac{\gamma}{2}\right)\cos\left(\frac{\gamma}{2}\right)+(n+1)\sin\left(\frac{\gamma}{2}\right)\cos\left(k\frac{\gamma}{2}\right)\right).\end{split}

(29)

The variables $\omega_{e}$ and $\omega^{om}$ can be expressed in terms of the Earth’s rotation frequency $f_{e}$ and the satellites’ orbit frequency $f_{o}$ , i.e.,

\left\{\begin{array}[]{l}\omega_{e}=-2\pi f_{e}t+\omega_{e0},\\ \omega^{om}=2\pi f_{o}t+\omega_{o0},\end{array}\right.

(30)

where $\omega_{e0}$ and $\omega_{o0}$ are the initial phases, and $f_{o}$ is determined by

f_{o}=\sqrt{\frac{GM}{4\pi^{2}a^{3}}}.

(31)

Inserting Eq. (30) into Eq. (26), we obtain

\mathcal{DA}=\sum_{m=0}^{N}\sum_{k=-N}^{N}\left[E_{mk}^{c}\cos\Upsilon_{mk}(t)-E_{mk}^{s}\sin\Upsilon_{mk}(t)\right],

(32)

with the angular argument

\Upsilon_{mk}(t)=2\pi(mf_{e}-kf_{o})t-m\omega_{e0}-k\omega_{o0}.

(33)

Or equivalently, one can rewrite Eq. (32) as

\mathcal{DA}=\sum_{m=0}^{N}\sum_{k=-N}^{N}E_{mk}\cos\left(\Upsilon_{mk}(t)-\omega_{mk}\right),

(34)

with

		$\displaystyle E_{mk}=\sqrt{(E_{mk}^{c})^{2}+(E_{mk}^{s})^{2}},$		(35)
		$\displaystyle\tan\omega_{mk}=-E_{mk}^{s}/E_{mk}^{c}.$		(36)

At last, one can model the range-acceleration between two TianQin satellites by

\ddot{\rho}=\!\!\sum_{m=0}^{N}\sum_{k=-N}^{N}\!\!\!E_{mk}\cos\left(\Upsilon_{mk}(t)-\omega_{mk}\right)+\frac{2GM\sin(\gamma/2)}{a^{2}},

(37)

with $E_{mk}$ , $\Upsilon_{mk}(t)$ , and $\omega_{mk}$ given by Eqs. (27–36). The formula consists of a series of cosine functions of time $t$ . Note that terms with $m=0$ and opposite $k$ are combined into one. Their amplitudes depend on the harmonic coefficients $\bar{C}_{nm}$ , $\bar{S}_{nm}$ and the satellite orbit elements $a$ and $i$ , as well as the phase difference $\gamma$ . The associated frequencies are integer linear combinations of the Earth’s rotation frequency $f_{e}$ and the satellites’ orbit frequency $f_{o}$ . Hence, through a Fourier-style decomposition in terms of frequency components, the formula can reveal the frequency-domain characteristics of the effect of the Earth’s gravity field { $\bar{C}_{nm}$ , $\bar{S}_{nm}$ } on the inter-satellite range-acceleration $\ddot{\rho}$ . Moreover, the dependence on the orbit elements $a,i,\Omega$ also allows one to study the impact of orbit selection [6]. To see explicit examples of Eq. (37), one can refer to Appendix C, which contains the low-degree case of $N=2$ .

II.4 Comments

In Kaula’s original linear perturbation method [10, 11], the orbital element representation Eq. (2) of the geo-potential for a satellite moving along a reference orbit is plugged into Lagrange’s perturbation equations to solve for orbital element perturbation relative to the reference values. The reference orbit is a slowly precessing circular orbit, i.e. that $\{a,e,i\}$ are of fixed values and that $\{\Omega,\omega,M\}$ are linear functions of $t$ . It is obtained also by solving Lagrange’s equations that contain the dominating $\bar{C}_{20}$ term. Once the orbital element perturbations are acquired, one can establish the relationship between the range or range-rate observables and the geo-potential coefficients [12, 13, 14], which can be used to determine the latter from the former for gravimetry.

Unfortunately, the solutions to Lagrange’s equations containing high-degree geo-potential terms have very complicated forms, cumbersome to use for our purposes. Moreover, what we need for the noise assessment is the range-acceleration variation in the frequency domain rather than the inter-satellite range or range-rate in the time domain. Therefore, instead of solving Lagrange’s equations, we directly use Kaula’s representation Eq. (2) and calculate the range acceleration through differential gravitational acceleration sensed by the satellite pair. The treatment is inspired by [15] and takes into account TianQin’s large orbit radius and long baseline. It should be mentioned that the small orbit precession has been neglected in our analytical model, because of the large orbit radius resulting in much weaker influence of $\bar{C}_{20}$ , as opposed to the case of low-orbit gravimetry missions. This has also been confirmed by the numerical work [2].

III Model Verification

In this section, we verify the analytical model with high-precision numerical orbit simulation. First we discuss the case of TianQin’s orbit and compare the resulting ASD curves. Then we alter the orbital elements to different values [6] so as to test the model in more generic settings.

From the previous work ([5], see Sec. IVB, Fig. 3), we have learned that the total range-acceleration ASD can be viewed as the sum of the lunisolar contribution and the Earth’s contribution. Since the latter dominates above $5\times 10^{-5}$ Hz in the total ASD and hence is more relevant to TianQin, we will focus on the Earth’s contribution alone in the verification. For an analytical model of the Sun and the Moon’s effect on the constellation stability, one may refer to [7].

III.1 TianQin’s orbit

In this case, we set the parameters of Eq. (37) according to Table 1 [32], and use EGM2008 [33] for the Earth’s static gravity field with the maximum degree $N=12$ . Then the range acceleration is calculated every 50s for a duration of 90 days. At the same time, we numerically integrate the orbits using the TQPOP program [5] with the force models listed in Table 2. Note that the Earth’s rotation models are also included to test our assumption of uniform Earth rotation in the model.

Table 1: The parameter setting of the analytical model [32].

Symbols	Parameters	Values
$N$	maximum degree	$12$
$a$	orbit radius	$1\times 10^{5}\mathrm{km}$
$i$	inclination	$74.5^{\circ}$
$\Omega$	longitude of ascending node	$211.6^{\circ}$
$GM$	Earth’s gravity constant	$3.986\times 10^{14}\mathrm{m}^{3}/\mathrm{s}^{2}$
$R$	average Earth radius	$6.378\times 10^{6}\mathrm{m}$
$T_{e}$	Earth rotation period	$86164\mathrm{s}$
$\omega_{o0}$	satellite initial phase	$\pi/2$

Table 2: Force models implemented in the numerical simulation.

Models	Specifications
Earth’s precession & nutation	IAU 2006/2000A [34]
Earth’s polar motion	EOP 14 C04 [35]
Earth’s static gravity field	EGM2008 ( $N=12$ ) [33]

The two ASD curves are compared in Fig. 2. Note that the flattening of the red curve above $\sim 10^{-4}$ Hz is due to numerical error of TQPOP (interpolation of the EOP data, likewise for the figures in Sec. III.2), while the green curve from Eq. (37) can extend below this error. Both curves agree well with each other except at the orbital frequency $3\times 10^{-6}$ Hz, and the much lower peak in the green curve is likely owing to the vanishing eccentricity of the circular reference orbit. However, one can see that the analytical model can indeed capture the main spectral feature of the effect due to the Earth’s static gravity field.

III.2 Other inclinations and radii

Firstly we reset the inclination to $30^{\circ}$ , $90^{\circ}$ , and $150^{\circ}$ in the analytical model while keeping the other parameters the same. The results are shown in Fig. 3-5, all of which show good agreement and are consistent with [6]. Additionally, we have also tested the cases of $60^{\circ}$ and $120^{\circ}$ . They all show similar features like Fig. 3-5, and hence omitted here.

Secondly we reset the radius to $0.8\times 10^{5}$ km and $1.2\times 10^{5}$ km in the analytical model while keeping the other parameters the same. The results are shown in Fig. 6-7. Again one has an evident matching of the two results.

Due to rotational symmetry about the Earth’s axis, we do not need to test different longitudes of ascending node. All these results have shown the validity of our model, as well as the approximations made in the derivation. It prompts us to take a further step and include time-variable gravity in the model.

IV Model Derivation: Time-Variable Gravity

The model Eq. (37) can be extended straightforwardly to incorporate the Earth’s time-variable gravity by adding time-dependent corrections to the harmonic coefficients $\bar{C}_{nm}$ and $\bar{S}_{nm}$ . The temporal gravity variations can be categorized into tidal and non-tidal parts. For the Earth’s tides (solid Earth, ocean, pole, and atmospheric), their analytical models in the form of the harmonic correction terms are well-known (e.g., see [34]) with frequency components typically outside TianQin’s measurement band of $10^{-4}-1$ Hz. As has been shown for TianQin, the tidal contributions to the total rang-acceleration ASD are considerably smaller than the static one [5]. Therefore for TianQin, it is more relevant to instead examine non-tidal gravity variations that may enter the measurement band.

From the classification of non-tidal variations, the Earth’s free oscillations are of particular interest since their frequencies are typically in the mHz range. Free oscillations of the Earth are standing waves at discrete frequencies of the solid Earth like a ringing bell ([36], not to be confused with seismic waves). They are typically excited and observed after a large earthquake, and can last for days. In this section, we model the effect of the Earth’s free oscillations triggered by earthquakes as a demonstrating example. Other types of non-tidal variations can be treated in a similar manner.

An earthquake can be characterized by the point-source double-couple model and associated parameters such as scalar seismic moment, dip, rake, and strike [37, 38]. The free oscillation excited by an earthquake is a superposition of the Earth’s normal modes with amplitudes determined from the parameters of the earthquake [36, 39]. The normal modes of the Earth can be calculated by the Preliminary Reference Earth Model (PREM), which is a non-rotating, elastic, spherically symmetric, one dimensional Earth model widely used in literature [40]. According to [31, 41], the gravity change due to free oscillations following an earthquake can be modeled by

\left\{\begin{array}[]{c}\bar{C}_{nm}(t)=\bar{C}_{nm}+\sum_{l=0}^{L}\phantom{}{}_{l}\Delta\bar{C}_{nm}\,\phantom{}{}_{l}\xi_{n}(t-t_{0}),\\ \bar{S}_{nm}(t)=\bar{S}_{nm}+\sum_{l=0}^{L}\phantom{}{}_{l}\Delta\bar{S}_{nm}\,\phantom{}{}_{l}\xi_{n}(t-t_{0}),\end{array}\right.

(38)

where $l$ marks free-oscillation overtones, $L$ is the maximum overtone, and $t_{0}$ is the earthquake occurrence time. The time dependence is given by

\phantom{}_{l}\xi_{n}(t)=\begin{cases}0,&t<0,\\ 1-\phantom{}_{l}\tau_{n}(t)\cos\left(2\pi\phantom{}_{l}f_{n}t\right),&t\geq 0,\\ \end{cases}

(39)

with the attenuation function

\phantom{}_{l}\tau_{n}(t)=\exp\left(-\frac{2\pi\phantom{}_{l}f_{n}t}{2\phantom{}_{l}Q_{n}}\right),

(40)

for each degree and overtone, where $\phantom{}{}_{l}f_{n}$ is the oscillation frequency, $\phantom{}{}_{l}Q_{n}$ is the attenuation factor, and $\phantom{}{}_{l}\Delta\bar{C}_{nm}$ and $\phantom{}{}_{l}\Delta\bar{S}_{nm}$ denote the permanent changes of the harmonic coefficients. Moreover, to be consistent with the notation of Eq. (4), we introduce

\begin{array}[]{l}\phantom{}{}_{l}\Delta\alpha_{nm}=\left\{\begin{array}[]{cc}\phantom{}{}_{l}\Delta\bar{C}_{nm},&n-m=\mathrm{even},\\ -\phantom{}_{l}\Delta\bar{S}_{nm},&n-m=\mathrm{odd},\end{array}\right.\\ \\ \phantom{}{}_{l}\Delta\beta_{nm}=\left\{\begin{array}[]{cc}\phantom{}{}_{l}\Delta\bar{S}_{nm},&n-m=\mathrm{even},\\ \phantom{}{}_{l}\Delta\bar{C}_{nm},&n-m=\mathrm{odd},\end{array}\right.\end{array}

(41)

The expressions for $\alpha_{nm}(t)$ and $\beta_{nm}(t)$ can be written down accordingly.

Including these time variations, the range-acceleration formula can be written as

\ddot{\rho}=\mathcal{DA}_{1}+\mathcal{DA}_{2}+\mathcal{DA}_{3}+\mathcal{CA}.

(42)

Here the first term $\mathcal{DA}_{1}$ is identical to $\mathcal{DA}$ in the absence of the free oscillations. The second term is given by

\mathcal{DA}_{2}=\sum_{m=0}^{N}\sum_{k=-N}^{N}P_{mk}\cos\left(\Upsilon_{mk}(t)-\varepsilon_{mk}\right),

(43)

with

P_{mk}=\sqrt{\left(P_{mk}^{c}\right)^{2}+\left(P_{mk}^{s}\right)^{2}},

(44)

and

\left\{\begin{array}[]{c}P_{mk}^{c}=\sum_{n=n_{1}[2]}^{n_{2}}\sum_{l=0}^{L}\Xi_{nmk}(a,i,\gamma)\phantom{}_{l}\Delta\alpha_{nm},\\ P_{mk}^{s}=\sum_{n=n_{1}[2]}^{n_{2}}\sum_{l=0}^{L}\Xi_{nmk}(a,i,\gamma)\phantom{}_{l}\Delta\beta_{nm},\end{array}\right.

(45)

and

\tan\varepsilon_{mk}=-P_{mk}^{s}/P_{mk}^{c}.

(46)

The term stands for the permanent change of $\ddot{\rho}$ resulting from an earthquake, where the frequency combinations $\Upsilon_{mk}(t)$ are identical to the ones in $\mathcal{DA}$ with the amplitudes and phases slightly altered. The third term carries the free oscillation information and reads

\begin{split}&\mathcal{DA}_{3}=\sum_{n=0}^{N}\sum_{l=0}^{L}\sum_{m=0}^{n}\sum_{k=-n[2]}^{n}\frac{1}{2}\phantom{}_{l}T_{nmk}\,\phantom{}{}_{l}\tau_{n}(t)\times\\ &\left[\cos\left(\phantom{}{}_{l}\Psi_{nmk}^{+}(t)-\phantom{}_{l}\gamma_{nmk}\right)+\cos\left(\phantom{}{}_{l}\Psi_{nmk}^{-}(t)+\phantom{}_{l}\gamma_{nmk}\right)\right],\end{split}

(47)

with the time dependence

\left\{\begin{array}[]{l}\begin{split}\phantom{}{}_{l}\Psi_{nmk}^{+}(t)=&2\pi\left(\phantom{}{}_{l}f_{n}+(mf_{e}-kf_{o})\right)t\\ &-2\pi\phantom{}_{l}f_{n}t_{0}-m\omega_{e0}-k\omega_{o0},\end{split}\\ \begin{split}\phantom{}{}_{l}\Psi_{nmk}^{-}(t)=&2\pi\left(\phantom{}{}_{l}f_{n}-(mf_{e}-kf_{o})\right)t\\ &-2\pi\phantom{}_{l}f_{n}t_{0}+m\omega_{e0}+k\omega_{o0},\end{split}\end{array}\right.

(48)

and

\phantom{}_{l}T_{nmk}=\sqrt{\left(\phantom{}{}_{l}T_{nmk}^{c}\right)^{2}+\left(\phantom{}{}_{l}T_{nmk}^{s}\right)^{2}},

(49)

where one has

\left\{\begin{array}[]{l}\phantom{}{}_{l}T_{nmk}^{c}=-\Xi_{nmk}(a,i,\gamma)\phantom{}_{l}\Delta\alpha_{nm},\\ \phantom{}{}_{l}T_{nmk}^{s}=-\Xi_{nmk}(a,i,\gamma)\phantom{}_{l}\Delta\beta_{nm},\end{array}\right.

(50)

and

\tan\phantom{}_{l}\gamma_{nmk}=-\phantom{}_{l}T_{nmk}^{s}/\phantom{}_{l}T_{nmk}^{c}.

(51)

The term represents a series of damped oscillations after the earthquake. Note that the resulting oscillation frequencies in $\ddot{\rho}$ are the free oscillation frequencies $\phantom{}{}_{l}f_{n}$ plus or minus the integer linear combinations of the Earth’s rotation frequency $f_{e}$ and the satellites’ orbit frequency $f_{o}$ . This indicates splitting of one mode frequency in the range-acceleration ASD due to the coupling with the Earth’s rotation and the satellites’ orbital motion. Furthermore, the amplitudes of the frequency components depend on the parameters $a,i,\gamma$ and the earthquake-induced geo-potential changes $\phantom{}{}_{l}\Delta\alpha_{nm}$ and $\phantom{}{}_{l}\Delta\beta_{nm}$ . For more details, one can refer to Appendix D, which shows an example of the spherical mode $\phantom{}{}_{0}S_{2}$ ( $[n,l]=[2,0]$ ) oscillation in the presence of the Earth’s static $N=2$ gravity field.

V Model Verification

Again we verify the extended model with high-precision numerical orbit simulation. Here we keep using the parameter setting from Sec. III.1. In addition, the moment magnitude of the test earthquake is 6.0 and the occurrence time $t_{0}$ is 30 days after the initial simulation time. The maximum degree $N$ and overtone $L$ are both set to be $12$ . The mode periods $\phantom{}{}_{l}T_{n}=1/\phantom{}_{l}f_{n}$ , the attenuation factor $\phantom{}{}_{l}Q_{n}$ , and the geo-potential corrections $\left\{\phantom{}{}_{l}\Delta\bar{C}_{nm},\phantom{}_{l}\Delta\bar{S}_{nm}\right\}$ are shown in Tables 3 and 4 (only prominent ones are shown, e.g., $\phantom{}{}_{0}S_{2},\phantom{}_{0}S_{3},\phantom{}_{1}S_{2},\phantom{}_{1}S_{3}$ ), which are obtained by following the procedures of [31].

As in the previous case of static gravity, we calculate the inter-satellite range-acceleration ASD from Eq. (42). At the same time, we obtain the numerical result from TQPOP in which the force models consist of the Earth’s static gravity and the gravity change from the free oscillations mentioned above. Note that the Earth’s polar motion is turned off to lower the numerical error above $\sim 10^{-4}$ due to data interpolation. The two ASD curves are compared in FIG.8, and the peaks above $2\times 10^{-4}$ Hz are due to the free oscillations. Both curves agree well with each other, indicating the validity of our model. The agreement also holds for other earthquakes we have tested but not shown here for succinctness. As an extra comment, it should be pointed out that spectrograms are normally preferred for visualizing non-stationary signals such as damped oscillations. Nevertheless, since our purpose is to compare the frequency content of the two results, we still use ASDs but one should keep in mind that the peak levels above $2\times 10^{-4}$ Hz do not reflect the actual amplitudes of the perturbed $\ddot{\rho}$ due to the free oscillations.

Furthermore, we make the comparison in the time domain as well. We calculate the TianQin’s range-acceleration time series using the numerical procedure and the analytical formula, respectively. Then we process the data with a high-pass filter [42] to extract signals of Earth’s free oscillation induced by earthquake above $10^{-4}$ Hz. The numerical and analytical waveforms are compared in Fig.9. The fitting factor is $>99\%$ , indicating a good consistence.

Table 3: Exemplary mode periods and attenuation factors of the tested free oscillations.

Parameters	Values	Parameters	Values
$\phantom{}{}_{0}T_{2}$	$3233.678\mathrm{s}$	$\phantom{}{}_{0}Q_{2}$	$509.6824$
$\phantom{}{}_{0}T_{3}$	$2134.577\mathrm{s}$	$\phantom{}{}_{0}Q_{3}$	$417.5499$
$\cdots$	$\cdots$	$\cdots$	$\cdots$
$\phantom{}{}_{1}T_{2}$	$1471.996\mathrm{s}$	$\phantom{}{}_{1}Q_{2}$	$310.2749$
$\phantom{}{}_{1}T_{3}$	$1064.875\mathrm{s}$	$\phantom{}{}_{1}Q_{3}$	$282.5018$
$\cdots$	$\cdots$	$\cdots$	$\cdots$

Table 4: Corrections to the geo-potential coefficients due to the earthquake, cf. Eq. (38).

Parameters	Values	Parameters	Values
$\phantom{}{}_{0}\Delta\bar{C}_{20}$	$\phantom{+}2.954\times 10^{-17}$	$\phantom{}{}_{0}\Delta\bar{S}_{20}$	$0$
$\phantom{}{}_{0}\Delta\bar{C}_{21}$	$\phantom{+}1.595\times 10^{-16}$	$\phantom{}{}_{0}\Delta\bar{S}_{21}$	$-2.313\times 10^{-16}$
$\phantom{}{}_{0}\Delta\bar{C}_{22}$	$\phantom{+}2.355\times 10^{-16}$	$\phantom{}{}_{0}\Delta\bar{S}_{22}$	$\phantom{+}1.078\times 10^{-16}$
$\phantom{}{}_{0}\Delta\bar{C}_{30}$	$-2.118\times 10^{-17}$	$\phantom{}{}_{0}\Delta\bar{S}_{30}$	$0$
$\phantom{}{}_{0}\Delta\bar{C}_{31}$	$-9.627\times 10^{-16}$	$\phantom{}{}_{0}\Delta\bar{S}_{31}$	$\phantom{+}6.953\times 10^{-16}$
$\phantom{}{}_{0}\Delta\bar{C}_{32}$	$\phantom{+}6.099\times 10^{-17}$	$\phantom{}{}_{0}\Delta\bar{S}_{32}$	$\phantom{+}6.132\times 10^{-17}$
$\phantom{}{}_{0}\Delta\bar{C}_{33}$	$\phantom{+}3.797\times 10^{-16}$	$\phantom{}{}_{0}\Delta\bar{S}_{33}$	$\phantom{+}8.375\times 10^{-16}$
$\cdots$	$\cdots$	$\cdots$	$\cdots$
$\phantom{}{}_{1}\Delta\bar{C}_{20}$	$-1.896\times 10^{-17}$	$\phantom{}{}_{1}\Delta\bar{S}_{20}$	$0$
$\phantom{}{}_{1}\Delta\bar{C}_{21}$	$\phantom{+}1.104\times 10^{-15}$	$\phantom{}{}_{1}\Delta\bar{S}_{21}$	$-8.167\times 10^{-16}$
$\phantom{}{}_{1}\Delta\bar{C}_{22}$	$\phantom{+}1.269\times 10^{-15}$	$\phantom{}{}_{1}\Delta\bar{S}_{22}$	$\phantom{+}3.188\times 10^{-16}$
$\phantom{}{}_{1}\Delta\bar{C}_{30}$	$\phantom{+}2.772\times 10^{-18}$	$\phantom{}{}_{1}\Delta\bar{S}_{30}$	$0$
$\phantom{}{}_{1}\Delta\bar{C}_{31}$	$\phantom{+}1.392\times 10^{-15}$	$\phantom{}{}_{1}\Delta\bar{S}_{31}$	$-1.043\times 10^{-15}$
$\phantom{}{}_{1}\Delta\bar{C}_{32}$	$-9.928\times 10^{-17}$	$\phantom{}{}_{1}\Delta\bar{S}_{32}$	$-4.501\times 10^{-17}$
$\phantom{}{}_{1}\Delta\bar{C}_{33}$	$-5.365\times 10^{-16}$	$\phantom{}{}_{1}\Delta\bar{S}_{33}$	$-1.237\times 10^{-15}$
$\cdots$	$\cdots$	$\cdots$	$\cdots$

VI Concluding Remarks and Outlook

In this work, we provide an analytical model which can account for the perturbing effect of the Earth’s gravity field on TianQin’s inter-satellite range acceleration. Especially, we have verified that the ASD curves obtained from the analytical model matches those from the numerical simulation. The work helps to confirm the numerical result of [5] and to better understand the underlying mechanism in terms of the frequency composition, as well as the dependence on orbit selection [6]. Particularly, we have also extended the model to capture the time-variable gravity from the Earth’s free oscillations, which enters TianQin’s measurement band.

For future applications in environment monitoring and noise-reduction pipelines, the model is useful in studying TianQin’s detector response to gravity disturbance from the Earth, particularly through time delay interferometry combinations [43]. Additionally, the model may also find usage in fast waveform generation for subtracting unwanted terrestrial gravity disturbance in case of occasional large-scale geo-seismic activities. Some related work will be reported in the future.

Acknowledgements.

The authors thank Hang Li, Xiang’e Lei, and Jun Luo for helpful discussions and comments, and special thanks to Kun Liu for providing earthquake parameters. X. Z. is supported by the National Key R&D Program of China (Grant No. 2020YFC2201202 and and 2022YFC2204600).

Appendix A Table of Symbols

Table LABEL:tab:symbols below lists the symbols used in the paper as well as their meanings for quick look-ups.

Table 5: List of symbols and their meanings

Symbols	Meanings
$V$	Earth’s gravity potential
$G$	universal gravitational constant
$M$	Earth’s total mass
$R$	Earth’s average radius
$n$	degree of spherical harmonic expansion
$m$	order of spherical harmonic expansion
$N$	truncation degree
$l$	Earth’s free-oscillation overtone
$L$	maximum overtone
$p$	index with values $0,1,\cdots,n$
$q$	index with values $-\infty,\cdots,-1,0,1,\cdots,+\infty$
$k$	index with values $-N,\cdots,-1,0,1,\cdots,N$
$n_{1},n_{2}$	see Eq.(12)
$t$	time
$t_{0}$	earthquake occurrence time
$r$	radius
$\theta$	co-latitude
$\lambda$	longitude
$\gamma$	phase difference between two satellites
$a$	semi-major axis
$i$	inclination
$e$	eccentricity
$\omega$	argument of perigee
$M$	mean anomaly
$\Omega$	right ascension of ascending node
$\Theta$	Greenwich Apparent Sidereal Time
$\bar{C}_{nm},\bar{S}_{nm}$	spherical harmonic coefficients
$\alpha_{nm}$ , $\beta_{nm}$	rearrangements of $\bar{C}_{nm}$ , $\bar{S}_{nm}$ defined by Eq.(4)
$\phantom{}{}_{l}\Delta\bar{C}_{nm},\phantom{}_{l}\Delta\bar{S}_{nm}$	coefficients of correction of the harmonic coefficients
$\phantom{}{}_{l}\Delta\alpha_{nm},\phantom{}_{l}\Delta\beta_{nm}$	rearrangement of $\phantom{}{}_{l}\Delta\bar{C}_{nm},\phantom{}_{l}\Delta\bar{S}_{nm}$ defined by Eq.(41)
$\bar{P}_{nm}(\cos\theta)$	normalized associated Legendre functions
$S_{nmpq}(\omega,M,\Omega,\Theta)$	see Eq.(3)
$\psi_{nmpq}$	see Eq.(5)
$\psi_{nmp}$	see Eq.(7)
$\psi_{mk}$	see Eq.(11)
$G_{npq}(e)$	functions of eccentricity
$\bar{F}_{nmp}(i)$	functions of inclination
$\bar{F}_{nm}^{k}(i)$	another form of functions of inclination defined by Eq.(14)
$u_{n}(a)$	functions of radius defined by Eq.(10)
$\omega^{o}$	see Eq.(8)
$\omega_{e}$	see Eq.(8)
$\omega^{om}$	argument of latitude of the midpoint of SC1,2
$f_{e}$	Earth’s rotation frequency
$f_{o}$	satellites’ orbit frequency obtained from Eq.(31)
$\omega_{e0}$	initial phases of Earth’s rotation
$\omega_{o0}$	initial phases of satellites’ movement
$\phantom{}{}_{l}f_{n}$	Earth’s free-oscillation frequency
$\phantom{}{}_{l}Q_{n}$	Earth’s free-oscillation attenuation factor
$E_{mk}$	see Eq.(35)
$E_{mk}^{c},E_{mk}^{s}$	see Eq.(27)
$\omega_{mk}$	see Eq.(35)
$P_{mk}$	see Eq.(44)
$P_{mk}^{c},P_{mk}^{s}$	see Eq.(45)
$\varepsilon_{mk}$	see Eq.(46)
$\phantom{}{}_{l}T_{nmk}$	see Eq.(49)
$\phantom{}{}_{l}T_{nmk}^{c},\phantom{}_{l}T_{nmk}^{s}$	see Eq.(50)
$\phantom{}{}_{l}\gamma_{nmk}$	see Eq.(51)
$\Xi_{nmk}(a,i,\gamma)$	see Eq.(29)
$\Upsilon_{mk}(t)$	see Eq.(33)
$\phantom{}{}_{l}\xi_{n}(t)$	see Eq.(39)
$\phantom{}{}_{l}\tau_{n}(t)$	attenuation function defined by Eq.(40)
$\phantom{}{}_{l}\Psi_{nmk}^{+}(t),\phantom{}_{l}\Psi_{nmk}^{-}(t)$	see Eq.(48)
$\mathcal{CA}$	centrifugal acceleration
$\mathcal{DA}$	projected differential gravitational acceleration
$\mathcal{DA}_{1}$	identical to $\mathcal{DA}$
$\mathcal{DA}_{2}$	permanent change of $\ddot{\rho}$ given by Eq.(43)
$\mathcal{DA}_{3}$	attenuation oscillation change of $\ddot{\rho}$ given by Eq.(47)
$\rho$	inter-satellite range
$\dot{\rho}$	inter-satellite range-rate
$\ddot{\rho}$	inter-satellite range-acceleration
$\mathbf{r}_{1},\mathbf{r}_{2}$	position vectors of SC1 and SC2
$\dot{\mathbf{r}}_{1},\dot{\mathbf{r}}_{2}$	velocity vectors of SC1 and SC2
$\hat{\mathbf{e}}_{12}$	unit vector from SC2 to SC1

Appendix B Explicit forms of inclination functions

A non-normalized version of $\bar{F}_{nmp}$ up to $n=4$ can be found in [11]. The explicit forms of normalized $\bar{F}_{nm}^{k}(i)$ up to $n=4$ used in this work are shown below, following the ascending order of $[m,k,n]$ .

\bar{F}_{4\,0}^{-4}(i)=\frac{105}{128}(\sin i)^{4}

(52)

\bar{F}_{3\,0}^{-3}(i)=\frac{5}{16}\sqrt{7}(\sin i)^{3}

(53)

\bar{F}_{2\,0}^{-2}(i)=-\frac{3}{8}\sqrt{5}(\sin i)^{2}

(54)

\bar{F}_{4\,0}^{-2}(i)=-\frac{105}{32}(\sin i)^{4}+\frac{45}{16}(\sin i)^{2}

(55)

\bar{F}_{1\,0}^{-1}(i)=-\frac{1}{2}\sqrt{3}\sin i

(56)

\bar{F}_{3\,0}^{-1}(i)=-\frac{15}{16}\sqrt{7}(\sin i)^{3}+\frac{3}{4}\sqrt{7}\sin i

(57)

\bar{F}_{0\,0}^{0}(i)=1

(58)

\bar{F}_{2\,0}^{0}(i)=\frac{3}{4}\sqrt{5}(\sin i)^{2}-\frac{1}{2}\sqrt{5}

(59)

\bar{F}_{4\,0}^{0}(i)=\frac{315}{64}(\sin i)^{4}-\frac{45}{8}(\sin i)^{2}+\frac{9}{8}

(60)

\bar{F}_{1\,0}^{1}(i)=\frac{1}{2}\sqrt{3}\sin i

(61)

\bar{F}_{3\,0}^{1}(i)=\frac{15}{16}\sqrt{7}(\sin i)^{3}-\frac{3}{4}\sqrt{7}\sin i

(62)

\bar{F}_{2\,0}^{2}(i)=-\frac{3}{8}\sqrt{5}(\sin i)^{2}

(63)

\bar{F}_{4\,0}^{2}(i)=-\frac{105}{32}(\sin i)^{4}+\frac{45}{16}(\sin i)^{2}

(64)

\bar{F}_{3\,0}^{3}(i)=-\frac{5}{16}\sqrt{7}(\sin i)^{3}

(65)

\bar{F}_{4\,0}^{4}(i)=\frac{105}{128}(\sin i)^{4}

(66)

\bar{F}_{4\,1}^{-4}(i)=-\frac{21}{64}\sqrt{10}(\sin i)^{3}\cos i+\frac{21}{64}\sqrt{10}(\sin i)^{3}

(67)

\bar{F}_{3\,1}^{-3}(i)=\frac{5}{32}\sqrt{42}(\sin i)^{2}\cos i-\frac{5}{32}\sqrt{42}(\sin i)^{2}

(68)

\bar{F}_{2\,1}^{-2}(i)=\frac{1}{4}\sqrt{15}\sin i\cos i-\frac{1}{4}\sqrt{15}\sin i

(69)

\begin{split}\bar{F}_{4\,1}^{-2}(i)&=\frac{21}{16}\sqrt{10}(\sin i)^{3}\cos i-\frac{9}{16}\sqrt{10}\sin i\cos i\\ &-\frac{21}{32}\sqrt{10}(\sin i)^{3}+\frac{9}{16}\sqrt{10}\sin i\end{split}

(70)

\bar{F}_{1\,1}^{-1}(i)=-\frac{1}{2}\sqrt{3}\cos i+\frac{1}{2}\sqrt{3}

(71)

\begin{split}\bar{F}_{3\,1}^{-1}(i)=&-\frac{15}{32}\sqrt{42}(\sin i)^{2}\cos i+\frac{5}{32}\sqrt{42}(\sin i)^{2}\\ &+\frac{1}{8}\sqrt{42}\cos i-\frac{1}{8}\sqrt{42}\end{split}

(72)

\bar{F}_{2\,1}^{0}(i)=-\frac{1}{2}\sqrt{15}\sin i\cos i

(73)

\bar{F}_{4\,1}^{0}(i)=-\frac{63}{32}\sqrt{10}(\sin i)^{3}\cos i+\frac{9}{8}\sqrt{10}\sin i\cos i

(74)

\bar{F}_{1\,1}^{1}(i)=\frac{1}{2}\sqrt{3}\cos i+\frac{1}{2}\sqrt{3}

(75)

\begin{split}\bar{F}_{3\,1}^{1}(i)&=\frac{15}{32}\sqrt{42}(\sin i)^{2}\cos i+\frac{5}{32}\sqrt{42}(\sin i)^{2}\\ &-\frac{1}{8}\sqrt{42}\cos i-\frac{1}{8}\sqrt{42}\end{split}

(76)

\bar{F}_{2\,1}^{2}(i)=\frac{1}{4}\sqrt{15}\sin i\cos i+\frac{1}{4}\sqrt{15}\sin i

(77)

\begin{split}\bar{F}_{4\,1}^{2}(i)&=\frac{21}{16}\sqrt{10}(\sin i)^{3}\cos i-\frac{9}{16}\sqrt{10}\sin i\cos i\\ &+\frac{21}{32}\sqrt{10}(\sin i)^{3}-\frac{9}{16}\sqrt{10}\sin i\end{split}

(78)

\bar{F}_{3\,1}^{3}(i)=-\frac{5}{32}\sqrt{42}(\sin i)^{2}\cos i-\frac{5}{32}\sqrt{42}(\sin i)^{2}

(79)

\bar{F}_{4\,1}^{4}(i)=-\frac{21}{64}\sqrt{10}(\sin i)^{3}\cos i-\frac{21}{64}\sqrt{10}(\sin i)^{3}

(80)

\begin{split}\bar{F}_{4\,2}^{-4}(i)=&-\frac{21}{64}\sqrt{5}(\sin i)^{2}(\cos i)^{2}+\frac{21}{32}\sqrt{5}(\sin i)^{2}\cos i\\ &-\frac{21}{64}\sqrt{5}(\sin i)^{2}\end{split}

(81)

\begin{split}\bar{F}_{3\,2}^{-3}(i)=&-\frac{1}{16}\sqrt{105}\sin i(\cos i)^{2}+\frac{1}{8}\sqrt{105}\sin i\cos i\\ &-\frac{1}{16}\sqrt{105}\sin i\end{split}

(82)

\bar{F}_{2\,2}^{-2}(i)=\frac{1}{8}\sqrt{15}(\cos i)^{2}-\frac{1}{4}\sqrt{15}\cos i+\frac{1}{8}\sqrt{15}

(83)

\begin{split}\bar{F}_{4\,2}^{-2}(i)&=\frac{21}{16}\sqrt{5}(\sin i)^{2}(\cos i)^{2}-\frac{21}{16}\sqrt{5}(\sin i)^{2}\cos i\\ &-\frac{3}{16}\sqrt{5}(\cos i)^{2}+\frac{3}{8}\sqrt{5}\cos i-\frac{3}{16}\sqrt{5}\end{split}

(84)

\begin{split}\bar{F}_{3\,2}^{-1}(i)&=\frac{3}{16}\sqrt{105}\sin i(\cos i)^{2}-\frac{1}{8}\sqrt{105}\sin i\cos i\\ &-\frac{1}{16}\sqrt{105}\sin i\end{split}

(85)

\bar{F}_{2\,2}^{0}(i)=-\frac{1}{4}\sqrt{15}(\cos i)^{2}+\frac{1}{4}\sqrt{15}

(86)

\begin{split}\bar{F}_{4\,2}^{0}(i)=&-\frac{63}{32}\sqrt{5}(\sin i)^{2}(\cos i)^{2}+\frac{21}{32}\sqrt{5}(\sin i)^{2}\\ &+\frac{3}{8}\sqrt{5}(\cos i)^{2}-\frac{3}{8}\sqrt{5}\end{split}

(87)

\begin{split}\bar{F}_{3\,2}^{1}(i)=&-\frac{3}{16}\sqrt{105}\sin i(\cos i)^{2}-\frac{1}{8}\sqrt{105}\sin i\cos i\\ &+\frac{1}{16}\sqrt{105}\sin i\end{split}

(88)

\bar{F}_{2\,2}^{2}(i)=\frac{1}{8}\sqrt{15}(\cos i)^{2}+\frac{1}{4}\sqrt{15}\cos i+\frac{1}{8}\sqrt{15}

(89)

\begin{split}\bar{F}_{4\,2}^{2}(i)&=\frac{21}{16}\sqrt{5}(\sin i)^{2}(\cos i)^{2}+\frac{21}{16}\sqrt{5}(\sin i)^{2}\cos i\\ &-\frac{3}{16}\sqrt{5}(\cos i)^{2}-\frac{3}{8}\sqrt{5}\cos i-\frac{3}{16}\sqrt{5}\end{split}

(90)

\begin{split}\bar{F}_{3\,2}^{3}(i)&=\frac{1}{16}\sqrt{105}\sin i(\cos i)^{2}+\frac{1}{8}\sqrt{105}\sin i\cos i\\ &+\frac{1}{16}\sqrt{105}\sin i\end{split}

(91)

\begin{split}\bar{F}_{4\,2}^{4}(i)=&-\frac{21}{64}\sqrt{5}(\sin i)^{2}(\cos i)^{2}-\frac{21}{32}\sqrt{5}(\sin i)^{2}\cos i\\ &-\frac{21}{64}\sqrt{5}(\sin i)^{2}\end{split}

(92)

Appendix C An explicit example for the Earth’s static gravity field

To help gain more intuition about Eq. (37), we exhibit its explicit form for $N=2$ static gravity field. In this case, one has

\begin{split}\ddot{\rho}&=2GM\sin(\gamma/2)/a^{2}+A_{00}+A_{02}\cos\left(2\pi f_{02}t+\varphi_{02}\right)\\ &+A_{10}\cos\left(2\pi f_{10}t+\varphi_{10}\right)+A_{20}\cos\left(2\pi f_{20}t+\varphi_{20}\right)\\ &+A_{1-2}\cos\left(2\pi f_{1-2}t+\varphi_{1-2}\right)+A_{12}\cos\left(2\pi f_{12}t+\varphi_{12}\right)\\ &+A_{2-2}\cos\left(2\pi f_{2-2}t+\varphi_{2-2}\right)+A_{22}\cos\left(2\pi f_{22}t+\varphi_{22}\right).\end{split}

(93)

As can be seen, besides the constant terms, there are several cosine terms of the time $t$ in Eq. (93). The amplitudes, frequencies and phases of these cosine term are given below.
(1) Term with the frequency 0:

\begin{split}A_{00}=&-2\frac{GM}{a^{2}}\sin\left(\frac{\gamma}{2}\right)\\ &-\frac{3\sqrt{5}}{2}\frac{GM}{R^{2}}\left(\frac{R}{a}\right)^{4}\left(3(\sin i)^{2}-2\right)\sin\left(\frac{\gamma}{2}\right)\bar{C}_{20}.\end{split}

(94)

(2) Term with the frequency $f_{02}$ :

f_{02}=2f_{o},

(95)

\begin{split}&A_{02}=\frac{3\sqrt{5}}{2}\frac{GM}{R^{2}}\left(\frac{R}{a}\right)^{4}(\sin i)^{2}\\ &\times\left(2\sin\gamma\cos\left(\frac{\gamma}{2}\right)+3\sin\left(\frac{\gamma}{2}\right)\cos\gamma\right)\bar{C}_{20},\end{split}

(96)

\varphi_{02}=2\omega_{o0}.

(97)

(3) Term with the frequency $f_{10}$ :

f_{10}=f_{e},

(98)

A_{10}=\sqrt{\left(E_{10}^{c}\right)^{2}+\left(E_{10}^{s}\right)^{2}},

(99)

\left\{\begin{array}[]{c}E_{10}^{c}=-3\sqrt{15}\frac{GM}{R^{2}}\left(\frac{R}{a}\right)^{4}\sin i\cos i\sin\left(\frac{\gamma}{2}\right)\bar{S}_{21},\\ E_{10}^{s}=+3\sqrt{15}\frac{GM}{R^{2}}\left(\frac{R}{a}\right)^{4}\sin i\cos i\sin\left(\frac{\gamma}{2}\right)\bar{C}_{21},\end{array}\right.

(100)

\varphi_{10}=-\omega_{e0}-\omega_{10},

(101)

\tan\omega_{10}=-E_{10}^{s}/E_{10}^{c}.

(102)

(4) Term with the frequency $f_{20}$ :

f_{20}=2f_{e},

(103)

A_{20}=\sqrt{\left(E_{20}^{c}\right)^{2}+\left(E_{20}^{s}\right)^{2}},

(104)

\left\{\begin{array}[]{c}E_{20}^{c}=\frac{3\sqrt{15}}{2}\frac{GM}{R^{2}}\left(\frac{R}{a}\right)^{4}\left((\cos i)^{2}-1\right)\sin\left(\frac{\gamma}{2}\right)\bar{C}_{22},\\ E_{20}^{s}=\frac{3\sqrt{15}}{2}\frac{GM}{R^{2}}\left(\frac{R}{a}\right)^{4}\left((\cos i)^{2}-1\right)\sin\left(\frac{\gamma}{2}\right)\bar{S}_{22},\end{array}\right.

(105)

\varphi_{20}=-2\omega_{e0}-\omega_{20},

(106)

\tan\omega_{20}=-E_{20}^{s}/E_{20}^{c}.

(107)

(5) Term with the frequency $f_{1\mp 2}$ :

f_{1\mp 2}=f_{e}\pm 2f_{o},

(108)

A_{1\mp 2}=\sqrt{\left(E_{1\mp 2}^{c}\right)^{2}+\left(E_{1\mp 2}^{s}\right)^{2}},

(109)

\left\{\begin{array}[]{c}\begin{split}&E_{1\mp 2}^{c}=+\frac{\sqrt{15}}{2}\frac{GM}{R^{2}}\left(\frac{R}{a}\right)^{4}\sin i\left(\cos i\mp 1\right)\\ &\times\left(2\sin\gamma\cos\left(\frac{\gamma}{2}\right)+3\sin\left(\frac{\gamma}{2}\right)\cos\gamma\right)\bar{S}_{21},\end{split}\\ \begin{split}&E_{1\mp 2}^{s}=-\frac{\sqrt{15}}{2}\frac{GM}{R^{2}}\left(\frac{R}{a}\right)^{4}\sin i\left(\cos i\mp 1\right)\\ &\times\left(2\sin\gamma\cos\left(\frac{\gamma}{2}\right)+3\sin\left(\frac{\gamma}{2}\right)\cos\gamma\right)\bar{C}_{21},\end{split}\end{array}\right.

(110)

\varphi_{1\mp 2}=-\omega_{e0}\pm 2\omega_{o0}-\omega_{1\mp 2},

(111)

\tan\omega_{1\mp 2}=-E_{1\mp 2}^{s}/E_{1\mp 2}^{c}.

(112)

(6) Term with the frequency $f_{2\mp 2}$ :

f_{2\mp 2}=2f_{e}\pm 2f_{o},

(113)

A_{2\mp 2}=\sqrt{\left(E_{2\mp 2}^{c}\right)^{2}+\left(E_{2\mp 2}^{s}\right)^{2}},

(114)

\left\{\begin{array}[]{c}\begin{split}&E_{2\mp 2}^{c}=-\frac{\sqrt{15}}{4}\frac{GM}{R^{2}}\left(\frac{R}{a}\right)^{4}\left(\cos i\mp 1\right)^{2}\\ &\times\left(2\sin\gamma\cos\left(\frac{\gamma}{2}\right)+3\sin\left(\frac{\gamma}{2}\right)\cos\gamma\right)\bar{C}_{22},\end{split}\\ \begin{split}&E_{2\mp 2}^{s}=-\frac{\sqrt{15}}{4}\frac{GM}{R^{2}}\left(\frac{R}{a}\right)^{4}\left(\cos i\mp 1\right)^{2}\\ &\times\left(2\sin\gamma\cos\left(\frac{\gamma}{2}\right)+3\sin\left(\frac{\gamma}{2}\right)\cos\gamma\right)\bar{S}_{22},\end{split}\end{array}\right.

(115)

\varphi_{1\mp 2}=-2\omega_{e0}\pm 2\omega_{o0}-\omega_{2\mp 2},

(116)

\tan\omega_{2\mp 2}=-E_{2\mp 2}^{s}/E_{2\mp 2}^{c}.

(117)

The resulting frequencies are integer linear combinations of the Earth’s rotation frequency $f_{e}$ and the satellites’ orbital frequency $f_{o}$ , i.e., 2 $f_{o}$ , $f_{e}$ , 2 $f_{e}$ , $f_{e}\mp 2f_{o}$ , $2f_{e}\mp 2f_{o}$ . The amplitudes are functions of the Earth’s gravity field coefficients $\bar{C}_{20}$ , $\bar{C}_{21}$ , $\bar{S}_{21}$ , $\bar{C}_{22}$ , $\bar{S}_{22}$ and TianQin’s orbital parameters $a$ , $i$ , $\gamma$ . The phases are determined by the initial phases of the Earth’s rotation and the satellite orbits, $\omega_{e0}$ and $\omega_{o0}$ , as well as the parameters $a$ , $i$ , $\gamma$ . The frequencies and associated amplitudes correspond to the peaks of the ASD curves.

Appendix D An explicit example for the Earth’s time-variable gravity field

To help gain more intuition about about Eq. (42), we exhibit its explicit form for $N=2$ static gravity field with the oscillation mode $\phantom{}{}_{0}S_{2}$ ( $[n,l]=[2,0]$ ). After the earthquake’s occurrence, the $N=2$ terms change into

\left\{\begin{array}[]{c}\bar{C}_{20}(t)=\bar{C}_{20}+\phantom{}_{0}\Delta\bar{C}_{20}\,\phantom{}{}_{0}\xi_{2}(t-t_{0}),\\ \bar{C}_{21}(t)=\bar{C}_{21}+\phantom{}_{0}\Delta\bar{C}_{21}\,\phantom{}{}_{0}\xi_{2}(t-t_{0}),\\ \bar{S}_{21}(t)=\bar{S}_{21}+\phantom{}_{0}\Delta\bar{S}_{21}\,\phantom{}{}_{0}\xi_{2}(t-t_{0}),\\ \bar{C}_{22}(t)=\bar{C}_{22}+\phantom{}_{0}\Delta\bar{C}_{22}\,\phantom{}{}_{0}\xi_{2}(t-t_{0}),\\ \bar{S}_{22}(t)=\bar{S}_{22}+\phantom{}_{0}\Delta\bar{S}_{22}\,\phantom{}{}_{0}\xi_{2}(t-t_{0}),\end{array}\right.

(118)

with the oscillation function

\phantom{}{}_{0}\xi_{2}(t-t_{0})=1-\phantom{}_{0}\tau_{2}(t)\cos\left(2\pi\phantom{}_{0}f_{2}(t-t_{0})\right)

(119)

and the attenuation

\phantom{}{}_{0}\tau_{2}(t)=\exp\left(-\frac{2\pi\phantom{}_{0}f_{2}(t-t_{0})}{2\phantom{}_{0}Q_{2}}\right).

(120)

The forms of $\mathcal{CA}$ and $\mathcal{DA}_{1}$ are identical to the ones for $N=2$ static gravity field. The term $\mathcal{DA}_{2}$ is negligible compared to $\mathcal{DA}_{1}$ , and hence omitted here. The term $\mathcal{DA}_{3}$ reads

\begin{split}&\mathcal{DA}_{3}=C_{00}\,\phantom{}{}_{0}\tau_{2}(t)\left[\cos\left(2\pi f_{00}^{+}t+\zeta_{00}^{+}\right)+\cos\left(2\pi f_{00}^{-}t+\zeta_{00}^{-}\right)\right]\\ &+C_{02}\,\phantom{}{}_{0}\tau_{2}(t)\left[\cos\left(2\pi f_{02}^{+}t+\zeta_{02}^{+}\right)+\cos\left(2\pi f_{02}^{-}t+\zeta_{02}^{-}\right)\right]\\ &+C_{10}\,\phantom{}{}_{0}\tau_{2}(t)\left[\cos\left(2\pi f_{10}^{+}t+\zeta_{10}^{+}\right)+\cos\left(2\pi f_{10}^{-}t+\zeta_{10}^{-}\right)\right]\\ &+C_{20}\,\phantom{}{}_{0}\tau_{2}(t)\left[\cos\left(2\pi f_{20}^{+}t+\zeta_{20}^{+}\right)+\cos\left(2\pi f_{20}^{-}t+\zeta_{20}^{-}\right)\right]\\ &+C_{1-2}\,\phantom{}{}_{0}\tau_{2}(t)\left[\cos\left(2\pi f_{1-2}^{+}t+\zeta_{1-2}^{+}\right)+\cos\left(2\pi f_{1-2}^{-}t+\zeta_{1-2}^{-}\right)\right]\\ &+C_{12}\,\phantom{}{}_{0}\tau_{2}(t)\left[\cos\left(2\pi f_{12}^{+}t+\zeta_{12}^{+}\right)+\cos\left(2\pi f_{12}^{-}t+\zeta_{12}^{-}\right)\right]\\ &+C_{2-2}\,\phantom{}{}_{0}\tau_{2}(t)\left[\cos\left(2\pi f_{2-2}^{+}t+\zeta_{2-2}^{+}\right)+\cos\left(2\pi f_{2-2}^{-}t+\zeta_{2-2}^{-}\right)\right]\\ &+C_{22}\,\phantom{}{}_{0}\tau_{2}(t)\left[\cos\left(2\pi f_{22}^{+}t+\zeta_{22}^{+}\right)+\cos\left(2\pi f_{22}^{-}t+\zeta_{22}^{-}\right)\right].\end{split}

(121)

It is a summation of several oscillation functions with the attenuation $\phantom{}{}_{0}\tau_{2}(t)$ . The frequency, amplitude, and phase of each term are listed below in their order of appearance in the equation.
(1) Terms with the frequencies $\left\{f_{00}^{+},f_{00}^{-}\right\}$ :

f_{00}^{+}=f_{00}^{-}=\phantom{}_{0}f_{2},

(122)

C_{00}=\frac{3\sqrt{5}}{4}\frac{GM}{R^{2}}\left(\frac{R}{a}\right)^{4}\left(3(\sin i)^{2}-2\right)\sin\left(\frac{\gamma}{2}\right)\phantom{}_{0}\Delta\bar{C}_{20},

(123)

\zeta_{00}^{+}=\zeta_{00}^{-}=-2\pi\phantom{}_{0}f_{2}t_{0}.

(124)

(2) Terms with the frequencies $\left\{f_{02}^{+},f_{02}^{-}\right\}$ :

\left\{\begin{array}[]{c}f_{02}^{+}=\phantom{}_{0}f_{2}-2f_{o},\\ f_{02}^{-}=\phantom{}_{0}f_{2}+2f_{o},\end{array}\right.

(125)

\begin{split}&C_{02}=-\frac{3\sqrt{5}}{4}\frac{GM}{R^{2}}\left(\frac{R}{a}\right)^{4}(\sin i)^{2}\\ &\times\left(2\sin\gamma\cos\left(\frac{\gamma}{2}\right)+3\sin\left(\frac{\gamma}{2}\right)\cos\gamma\right)\phantom{}_{0}\Delta\bar{C}_{20},\end{split}

(126)

\left\{\begin{array}[]{c}\zeta_{02}^{+}=-2\pi\phantom{}_{0}f_{2}t_{0}-2\omega_{o0},\\ \zeta_{02}^{-}=-2\pi\phantom{}_{0}f_{2}t_{0}+2\omega_{o0}.\end{array}\right.

(127)

(3) Terms with the frequencies $\left\{f_{10}^{+},f_{10}^{-}\right\}$ :

\left\{\begin{array}[]{c}f_{10}^{+}=\phantom{}_{0}f_{2}+f_{e},\\ f_{10}^{-}=\phantom{}_{0}f_{2}-f_{e},\end{array}\right.

(128)

C_{10}=\frac{1}{2}\sqrt{\left(\phantom{}{}_{0}T_{210}^{c}\right)^{2}+\left(\phantom{}{}_{0}T_{210}^{s}\right)^{2}},

(129)

\left\{\begin{array}[]{c}\phantom{}{}_{0}T_{210}^{c}=+3\sqrt{15}\frac{GM}{R^{2}}\left(\frac{R}{a}\right)^{4}\sin i\cos i\sin\left(\frac{\gamma}{2}\right)\phantom{}_{0}\Delta\bar{S}_{21},\\ \phantom{}{}_{0}T_{210}^{s}=-3\sqrt{15}\frac{GM}{R^{2}}\left(\frac{R}{a}\right)^{4}\sin i\cos i\sin\left(\frac{\gamma}{2}\right)\phantom{}_{0}\Delta\bar{C}_{21},\end{array}\right.

(130)

\left\{\begin{array}[]{c}\zeta_{10}^{+}=-2\pi\phantom{}_{0}f_{2}t_{0}-\omega_{e0}-\phantom{}_{0}\gamma_{210},\\ \zeta_{10}^{-}=-2\pi\phantom{}_{0}f_{2}t_{0}+\omega_{e0}+\phantom{}_{0}\gamma_{210},\end{array}\right.

(131)

\tan\phantom{}_{0}\gamma_{210}=-\phantom{}_{0}T_{210}^{s}/\phantom{}_{0}T_{210}^{c}.

(132)

(4) Terms with the frequencies $\left\{f_{20}^{+},f_{20}^{-}\right\}$ :

\left\{\begin{array}[]{c}f_{20}^{+}=\phantom{}_{0}f_{2}+2f_{e},\\ f_{20}^{-}=\phantom{}_{0}f_{2}-2f_{e},\end{array}\right.

(133)

C_{20}=\frac{1}{2}\sqrt{\left(\phantom{}{}_{0}T_{220}^{c}\right)^{2}+\left(\phantom{}{}_{0}T_{220}^{s}\right)^{2}},

(134)

\left\{\begin{array}[]{c}\begin{split}\phantom{}{}_{0}T_{220}^{c}=&-\frac{3\sqrt{15}}{2}\frac{GM}{R^{2}}\left(\frac{R}{a}\right)^{4}\\ &\times\left((\cos i)^{2}-1\right)\sin\left(\frac{\gamma}{2}\right)\phantom{}_{0}\Delta\bar{C}_{22},\end{split}\\ \begin{split}\phantom{}{}_{0}T_{220}^{s}=&-\frac{3\sqrt{15}}{2}\frac{GM}{R^{2}}\left(\frac{R}{a}\right)^{4}\\ &\times\left((\cos i)^{2}-1\right)\sin\left(\frac{\gamma}{2}\right)\phantom{}_{0}\Delta\bar{S}_{22},\end{split}\end{array}\right.

(135)

\left\{\begin{array}[]{c}\zeta_{20}^{+}=-2\pi\phantom{}_{0}f_{2}t_{0}-2\omega_{e0}-\phantom{}_{0}\gamma_{220},\\ \zeta_{20}^{-}=-2\pi\phantom{}_{0}f_{2}t_{0}+2\omega_{e0}+\phantom{}_{0}\gamma_{220},\end{array}\right.

(136)

\tan\phantom{}_{0}\gamma_{220}=-\phantom{}_{0}T_{220}^{s}/\phantom{}_{0}T_{220}^{c}.

(137)

(5) Terms with the frequencies $\left\{f_{1\mp 2}^{+},f_{1\mp 2}^{-}\right\}$ :

\left\{\begin{array}[]{c}f_{1\mp 2}^{+}=\phantom{}_{0}f_{2}+f_{e}\pm 2f_{o},\\ f_{1\mp 2}^{-}=\phantom{}_{0}f_{2}-f_{e}\mp 2f_{o},\end{array}\right.

(138)

C_{1\mp 2}=\frac{1}{2}\sqrt{\left(\phantom{}{}_{0}T_{21\mp 2}^{c}\right)^{2}+\left(\phantom{}{}_{0}T_{21\mp 2}^{s}\right)^{2}},

(139)

\left\{\begin{array}[]{c}\begin{split}&\phantom{}{}_{0}T_{21\mp 2}^{c}=-\frac{\sqrt{15}}{2}\frac{GM}{R^{2}}\left(\frac{R}{a}\right)^{4}\sin i\left(\cos i\mp 1\right)\\ &\times\left(2\sin\gamma\cos\left(\frac{\gamma}{2}\right)+3\sin\left(\frac{\gamma}{2}\right)\cos\gamma\right)\phantom{}_{0}\Delta\bar{S}_{21},\end{split}\\ \begin{split}&\phantom{}{}_{0}T_{21\mp 2}^{s}=+\frac{\sqrt{15}}{2}\frac{GM}{R^{2}}\left(\frac{R}{a}\right)^{4}\sin i\left(\cos i\mp 1\right)\\ &\times\left(2\sin\gamma\cos\left(\frac{\gamma}{2}\right)+3\sin\left(\frac{\gamma}{2}\right)\cos\gamma\right)\phantom{}_{0}\Delta\bar{C}_{21},\end{split}\end{array}\right.

(140)

\left\{\begin{array}[]{c}\zeta_{1\mp 2}^{+}=-2\pi\phantom{}_{0}f_{2}t_{0}-\omega_{e0}\pm 2\omega_{o0}-\phantom{}_{0}\gamma_{21\mp 2},\\ \zeta_{1\mp 2}^{-}=-2\pi\phantom{}_{0}f_{2}t_{0}+\omega_{e0}\mp 2\omega_{o0}+\phantom{}_{0}\gamma_{21\mp 2},\end{array}\right.

(141)

\tan\phantom{}_{0}\gamma_{21\mp 2}=-\phantom{}_{0}T_{21\mp 2}^{s}/\phantom{}_{0}T_{21\mp 2}^{c}.

(142)

(6) Terms with the frequencies $\left\{f_{2\mp 2}^{+},f_{2\mp 2}^{-}\right\}$ :

\left\{\begin{array}[]{c}f_{2\mp 2}^{+}=\phantom{}_{0}f_{2}+2f_{e}\pm 2f_{o},\\ f_{2\mp 2}^{-}=\phantom{}_{0}f_{2}-2f_{e}\mp 2f_{o},\end{array}\right.

(143)

C_{2\mp 2}=\frac{1}{2}\sqrt{\left(\phantom{}{}_{0}T_{22\mp 2}^{c}\right)^{2}+\left(\phantom{}{}_{0}T_{22\mp 2}^{s}\right)^{2}},

(144)

\left\{\begin{array}[]{c}\begin{split}&\phantom{}{}_{0}T_{22\mp 2}^{c}=-\frac{\sqrt{15}}{4}\frac{GM}{R^{2}}\left(\frac{R}{a}\right)^{4}\left(\cos i\mp 1\right)^{2}\\ &\times\left(2\sin(\gamma)\cos(\frac{\gamma}{2})+3\sin\left(\frac{\gamma}{2}\right)\cos(\gamma)\right)\phantom{}_{0}\Delta\bar{C}_{22},\end{split}\\ \begin{split}&\phantom{}{}_{0}T_{22\mp 2}^{s}=-\frac{\sqrt{15}}{4}\frac{GM}{R^{2}}\left(\frac{R}{a}\right)^{4}\left(\cos i\mp 1\right)^{2}\\ &\times\left(2\sin(\gamma)\cos(\frac{\gamma}{2})+3\sin\left(\frac{\gamma}{2}\right)\cos(\gamma)\right)\phantom{}_{0}\Delta\bar{S}_{22},\end{split}\end{array}\right.

(145)

\left\{\begin{array}[]{c}\zeta_{2\mp 2}^{+}=-2\pi\phantom{}_{0}f_{2}t_{0}-2\omega_{e0}\pm 2\omega_{o0}-\phantom{}_{0}\gamma_{22\mp 2},\\ \zeta_{2\mp 2}^{-}=-2\pi\phantom{}_{0}f_{2}t_{0}+2\omega_{e0}\mp 2\omega_{o0}+\phantom{}_{0}\gamma_{22\mp 2},\end{array}\right.

(146)

\tan\phantom{}_{0}\gamma_{22\mp 2}=-\phantom{}_{0}T_{22\mp 2}^{s}/\phantom{}_{0}T_{22\mp 2}^{c}.

(147)

The resulting frequencies are integer linear combinations of the ${}_{0}S_{2}$ mode frequency ${}_{0}f_{2}$ , the Earth’s rotation frequency $f_{e}$ , and the satellites’ orbital frequency $f_{o}$ , i.e., ${}_{0}f_{2}$ , ${}_{0}f_{2}\mp 2f_{o}$ , ${}_{0}f_{2}\pm f_{e}$ , ${}_{0}f_{2}\pm 2f_{e}$ , ${}_{0}f_{2}\pm\left(f_{e}\mp 2f_{o}\right)$ , ${}_{0}f_{2}\pm\left(2f_{e}\mp 2f_{o}\right)$ . The amplitudes are functions of the coefficients $\phantom{}{}_{0}\Delta\bar{C}_{20}$ , $\phantom{}{}_{0}\Delta\bar{C}_{21}$ , $\phantom{}{}_{0}\Delta\bar{S}_{21}$ , $\phantom{}{}_{0}\Delta\bar{C}_{22}$ , $\phantom{}{}_{0}\Delta\bar{S}_{22}$ and TianQin’s orbital parameters $a$ , $i$ , $\gamma$ . The phases are determined by ${}_{0}f_{2}$ , the earthquake occurrence time $t_{0}$ , the initial phases of the Earth’s rotation and the satellite orbits, $\omega_{e0}$ and $\omega_{o0}$ , as well as the parameters $a$ , $i$ , $\gamma$ . The frequencies and associated amplitudes correspond to the peaks of the ASD curves in the presence of the free oscillations.

References

J. Luo et al. [2016] J. Luo et al., TianQin: a space-borne gravitational wave detector, Class. Quantum Grav. 33, 035010 (2016).
B.-B. Ye et al. [2019] B.-B. Ye et al., Optimizing orbits for TianQin, Int. J. Mod. Phys. D 28, 1950121 (2019).
Z. Tan, B. Ye and X. Zhang [2020] Z. Tan, B. Ye and X. Zhang, Impact of orbital orientations and radii on TianQin constellation stability, Int. J. Mod. Phys. D 29, 2050056 (2020).
B.-B. Ye et al. [2021] B.-B. Ye et al., Eclipse avoidance in TianQin orbit selection, Phys. Rev. D 103, 042007 (2021).
X. Zhang et al. [2021] X. Zhang et al., Effect of earth-moon’s gravity on TianQin’s range acceleration noise, Phys. Rev. D 103, 062001 (2021).
Luo and Zhang [2022] C. Luo and X. Zhang, Effect of earth-moon’s gravity on TianQin’s range acceleration noise. II. impact of orbit selection, Phys. Rev. D 105, 102007 (2022).
Ye [2022] B. Ye, Theoretical analysis and optimization of TianQin constellation, Ph.D. thesis, Sun Yet-Sen University (2022).
B. D. Tapley and S. Bettadpur [2004] B. D. Tapley and S. Bettadpur, The gravity recovery and climate experiment: Mission overview and early results, Geophys. Res. Lett. 31, L09607 (2004).
R. P. Kornfeld et al. [2019] R. P. Kornfeld et al., Grace-fo: The gravity recovery and climate experiment follow-on mission, J. Spacecr. Rockets 56, 931 (2019).
Kaula [1961] W. M. Kaula, Analysis of gravitational and geometric aspects of geodetic utilization of satellites, Geophys. J. Int. 5, 104 (1961).
Kaula [1966] W. M. Kaula, Theory of Satellite Geodesy (Blaisdell Publishing Company, 1966).
Rosborough [1986] G. W. Rosborough, Satellite Orbit Perturbations Due to the Geopotential, Ph.D. thesis, The University of Texas at Austin (1986).
Rosborough [1987] G. W. Rosborough, Radial, transverse and normal satellite position perturbations due to the geopotential, Celest. Mech. 40, 409 (1987).
Sharma [1995] J. Sharma, Precise determination of the geopotential with a low-low satellite-to-satellite tracking mission, Ph.D. thesis, The University of Texas at Austin (1995).
Sharifi [2006] M. A. Sharifi, Satellite to Satellite Tracking in the Space-wise Approach, Ph.D. thesis, University of Stuttgart (2006).
Colombo [1984] O. L. Colombo, The global mapping of gravity with two satellites, Publications on Geodesy 7-3 (Netherlands Geodetic Commission, 1984).
M. K. Cheng [2002] M. K. Cheng, Gravitational perturbation theory for inter‐satellite tracking, J. Geodesy 76, 169 (2002).
C. Wagner et al. [2006] C. Wagner et al., Degradation of geopotential recovery from short repeat-cycle orbits: application to grace monthly fields, J. Geodesy 80, 94 (2006).
Schrama [1990] E. Schrama, Gravity field error analysis: Applications of GPS receivers and gradiometers on low orbiting platforms, NASA Technical Memorandum 100769 (National Aeronautics and Space Administration, 1990).
E. J. O. Schrama [1991] E. J. O. Schrama, Gravity field error analysis: Applications of Global Positioning System receivers and gradiometers on low orbiting platforms, J. Geophys. Res. 96, 20041 (1991).
Schrama [1989] E. J. O. Schrama, The role of orbit errors in processing of satellite altimeter data, Ph.D. thesis, Delft University of Technology (1989).
P. N. A. M. Visser, K. F. Wakker, B. A. C. Ambrosius [1994] P. N. A. M. Visser, K. F. Wakker, B. A. C. Ambrosius, Global gravity field recovery from the ARISTOTELES satellite mission, J. Geophys. Res. 99, 2841 (1994).
P. N. A. M. Visser et al. [2001] P. N. A. M. Visser et al., Exploring gravity field determination from orbit perturbations of the european gravity mission GOCE, J. Geodesy 75, 89 (2001).
P. N. A. M. Visser [2005] P. N. A. M. Visser, Low-low satellite-to-satellite tracking: a comparison between analytical linear orbit perturbation theory and numerical integration, J. Geodesy 79, 160 (2005).
Visser [1992] P. N. A. M. Visser, The use of satellites in gravity field determination and model adjustment, Ph.D. thesis, Delft University of Technology (1992).
Keller and Sharifi [2005] W. Keller and M. A. Sharifi, Satellite gradiometry using a satellite pair, J. Geodesy 78, 544 (2005).
Hajela [1974] D. P. Hajela, Improved procedures for the recovery of $5^{\circ}$ mean gravity anomalies from ATS-6/GEOS-3 satellite-to-satellite range-rate observation, Scientific Report 276 (Department of Geodetic Science, The Ohio State University, 1974).
Rummel [1980] R. Rummel, Geoid height, geoid height differences, and mean gravity anomalies from ‘low–low’ satellite-to-satellite tracking—an error analysis, Scientific Report 306 (Department of Geodetic Science, The Ohio State University, 1980).
Garcia [2002] R. V. Garcia, Local geoid determination from GRACE mission, Scientific Report 460 (Department of Civil and Environmental Engineering and Geodetic Science, The Ohio State University, 2002).
Han [2003] S.-C. Han, Efficient global gravity determination from satellite-to-satellite tracking (SST), Scientific Report 467 (Department of Civil and Environmental Engineering and Geodetic Science, The Ohio State University, 2003).
K. Ghobadi-Far et al. [2019] K. Ghobadi-Far et al., Gravitational changes of the earth’s free oscillation from earthquakes: Theory and feasibility study using GRACE inter‐satellite tracking, J. Geophys. Res.-Solid Earth 124, 7483 (2019).
Montenbruck and Gill [2000] O. Montenbruck and E. Gill, Satellite Orbits - Models, Methods, and Applications (Springer, 2000).
Pavlis et al. [2012] N. K. Pavlis, S. A. Holmes, S. C. Kenyon, and J. K. Factor, The development and evaluation of the Earth Gravitational Model 2008 (EGM2008), J. Geophys. Res. 117, B04406 (2012).
Petit and Luzum [2010] G. Petit and B. Luzum, IERS Conventions (2010), Technical Report 36 (BUREAU INTERNATIONAL DES POIDS ET MESURES SEVRES (FRANCE), 2010).
[35] http://hpiers.obspm.fr/iers/eop/eopc04/ .
Z. Alterman, H. Jarosch, C. L. Pekeris [1959] Z. Alterman, H. Jarosch, C. L. Pekeris, Oscillations of the earth, Proc. R. Soc. Lond. A 252, 80 (1959).
H. Kanamori and J. J. Cipar [1974] H. Kanamori and J. J. Cipar, Focal process of the great chilean earthquake may 22, 1960, Phys. Earth Planet. Inter. 9, 128 (1974).
H. Kanamori and G. W. Given [1981] H. Kanamori and G. W. Given, Use of long-period surface waves for rapid determination of earthquake-source parameters, Phys. Earth Planet. Inter. 27, 8 (1981).
Petrov [2021] O. V. Petrov, The Earth’s Free Oscillations-Formulation and Solution of the Fundamental Wave Equation of Nature (Springer, 2021).
A. M. Dziewonski and D. L. Anderson [1981] A. M. Dziewonski and D. L. Anderson, Preliminary reference earth model, Phys. Earth Planet. Inter. 25, 297 (1981).
S.-H. Han et al. [2013] S.-H. Han et al., Source parameter inversion for recent great earthquakes from a decade-long observation of global gravity fields, J. Geophys. Res.-Solid Earth 118, 1240 (2013).
Zheng et al. [2023] L. Zheng, S. Yang, and X. Zhang, Doppler effect in TianQin time-delay interferometry (2023), under review.
Tinto and Dhurandhar [2021] M. Tinto and S. V. Dhurandhar, Time-delay interferometry, Living Rev Relativ 24, 1 (2021).

Effect of Earth-Moon’s gravity on TianQin’s range acceleration noise. III. An analytical model

Abstract

I Introduction

II Model Derivation: Static Gravity

II.1 Basic mathematical setup

II.2 Circular reference orbits

II.3 Range acceleration

II.4 Comments

III Model Verification

III.1 TianQin’s orbit

III.2 Other inclinations and radii

IV Model Derivation: Time-Variable Gravity

V Model Verification

VI Concluding Remarks and Outlook

Acknowledgements.

Appendix A Table of Symbols

Appendix B Explicit forms of inclination functions

Appendix C An explicit example for the Earth’s static gravity field

Appendix D An explicit example for the Earth’s time-variable gravity field

References

Effect of Earth-Moon’s gravity on TianQin’s range acceleration noise. III.
An analytical model