Effects of lunisolar perturbations on TianQin constellation: An analytical model
阴阳扰动对天琴星座的影响：一个解析模型

Bobing Ye yebb5@mail.sysu.edu.cn Xuefeng Zhang zhangxf38@sysu.edu.cn 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, People’s Republic of China

(May 2, 2024) （5月 2， 2024）

Abstract 抽象

TianQin is a proposed space-based gravitational-wave observatory mission that critically relies on the stability of an equilateral-triangle constellation. Comprising three satellites in high Earth orbits of a $10^{5}$ km radius, this constellation’s geometric configuration is significantly affected by gravitational perturbations, primarily originating from the Moon and the Sun. In this paper, we present an analytical model to quantify the effects of lunisolar perturbations on the TianQin constellation, derived using Lagrange’s planetary equations. The model provides expressions for three kinematic indicators of the constellation: arm-lengths, relative line-of-sight velocities, and breathing angles. Analysis of these indicators reveals that lunisolar perturbations can distort the constellation triangle, resulting in three distinct variations: linear drift, bias, and fluctuation. Furthermore, it is shown that these distortions can be optimized to display solely fluctuating behavior, under certain predefined conditions. These results can serve as the theoretical foundation for numerical simulations and offer insights for engineering a stable constellation in the future.
“天琴”是一项拟议的天基引力波天文台任务，它严重依赖于等边三角形星座的稳定性。该星座由三颗位于半径 1 $10^{5}$ 公里的地球高轨道上的卫星组成，其几何结构受到引力扰动的显着影响，引力扰动主要源自月球和太阳。在本文中，我们提出了一个解析模型来量化阴阳扰动对天琴星座的影响，该模型使用拉格朗日行星方程推导。该模型提供了星座的三个运动学指标的表达式：臂长、相对视距速度和呼吸角。对这些指标的分析表明，阴阳扰动会扭曲星座三角形，导致三种不同的变化：线性漂移、偏差和波动。此外，结果表明，在某些预定义条件下，这些失真可以优化为仅显示波动行为。这些结果可以作为数值模拟的理论基础，并为未来设计稳定的星座提供见解。

^†^†preprint: APS/123-QED
^†预印本：APS/123-QED

I Introduction I 引言

The successful detection of gravitational waves (GWs) by the ground-based observatory LIGO [1] has opened up the era of GW astronomy. To detect GWs in the millihertz range (0.1 mHz–1 Hz), known for its rich sources and to circumvent the impact of seismic noise, space-based GW observatories are highly favored [2, 3]. For such observatories, proposed projects include LISA [4, 5], DECIGO [6], TianQin [7], Taiji [8], etc. Among these, TianQin is a geocentric space-based GW observatory mission that consists of three drag-free controlled satellites with an orbital radius of $10^{5}$ km [7]. The three satellites form a nearly equilateral-triangle constellation, standing almost vertical to the ecliptic, and they employ high-precision laser-ranging interferometry to measure distance changes between satellites for the detection of GWs. The mission will bring rich science prospects to GW astronomy [9, 10, 11].
地面天文台 LIGO [1] 成功探测到引力波（GW），开启了 GW 天文学的时代。为了探测以丰富的来源而闻名的毫赫兹范围（0.1 mHz–1 Hz）的 GW 并规避地震噪声的影响，天基 GW 天文台受到高度青睐 [2， 3]。对于此类天文台，拟议的项目包括 LISA [4， 5]、DECIGO [6]、天琴 [7]、太极 [8] 等。其中，“天琴”是一项地心天基GW天文台任务，由三颗轨道半径为 $10^{5}$ 公里的无阻力受控卫星组成[7]。这三颗卫星形成一个几乎等边三角形的星座，几乎垂直于黄道，它们采用高精度激光测距干涉测量法测量卫星之间的距离变化，以检测 GW。该任务将为 GW 天文学带来丰富的科学前景 [9， 10， 11]。

TianQin, as well as other three-satellite GW missions, relies crucially on the stability of an equilateral-triangle constellation [3, 7]. Unequal variations in the three arm-lengths of the constellation prevent the cancellation of laser frequency noise, which has a profound impact on the design of frequency stabilization systems and requires time-delay interferometry (TDI) [12, 13, 14]. The relative line-of-sight velocities between satellites induce Doppler frequency shifts, affecting phase meter bandwidth and ultra-stable oscillator design [15]. Moreover, changes in the three breathing angles of the triangle directly influence the design of telescopes and beam pointing mechanisms [7]. It is crucial to minimize variations in the triangular constellation, as indicated by these three kinematic indicators.
“天琴”以及其他三颗卫星的 GW 任务在很大程度上依赖于等边三角形星座的稳定性 [3， 7]。星座三个臂长的不相等变化阻止了激光频率噪声的抵消，这对稳频系统的设计有深远的影响，需要时延干涉测量法（TDI） [12， 13， 14]。卫星之间的相对视距速度会引起多普勒频移，影响相位计带宽和超稳定振荡器设计 [15]。此外，三角形三个呼吸角的变化直接影响望远镜和光束指向机构的设计 [7]。正如这三个运动学指标所表明的那样，尽量减少三角形星座的变化至关重要。

Analytical analysis of satellite motion and constellation variations holds significant importance [16, 17, 18, 19, 20, 17]. To identify orbits with minimal variations in the constellation, extensive efforts have been dedicated to numerical orbit optimization and analysis (for a review, see Ref. [19]). The use of analytical models, as opposed to numerical simulations, allows for deeper physical insights and often yields clearer solutions for issues related to satellite motion [16, 17]. Moreover, these analytical models provide the basis for further numerical simulations, enhancing orbit optimization efficiency [18, 19]. They also enable theoretical studies on inter-satellite optical links and light propagation [20, 17].
卫星运动和星座变化的分析具有重要意义[16， 17， 18， 19， 20， 17]。为了确定星座中变化最小的轨道，人们投入了大量精力进行数值轨道优化和分析（综述见参考文献 [19]）。与数值模拟相反，使用分析模型可以获得更深入的物理见解，并且通常可以为与卫星运动相关的问题提供更清晰的解决方案[16,17]。此外，这些解析模型为进一步的数值模拟提供了基础，提高了轨道优化效率 [18， 19]。它们还有助于对星间光链路和光传播进行理论研究 [20， 17]。

Concerning analytical efforts, Ref. [21] first presented the analytical coordinates of the TianQin satellites, based on unperturbed Keplerian orbits, which showed that the arm-lengths of the constellation remain constant when orbital eccentricities are ignored. Furthermore, the leading-order effect of the third-body perturbation was considered to derive expressions for both arm-lengths and breathing angles [22]. These expressions were constructed iteratively, assuming circular orbits, and they were used to study the impact of initial orbit errors. Moreover, the effect of the Earth’s non-spherical gravitational perturbation was analyzed in [23], with a particular focus on its influence on inter-satellite range acceleration noise.
在分析工作方面，参考文献[21]首先提出了天琴卫星的解析坐标，该坐标基于未受扰动的开普勒轨道，结果表明，当忽略轨道偏心率时，星座的臂长保持不变。此外，第三体扰动的超前序效应被认为可以推导出臂长和呼吸角的表达式 [22]。这些表达式是迭代构建的，假设圆轨道，它们被用来研究初始轨道误差的影响。此外，在 [23] 中分析了地球非球面引力扰动的影响，特别关注其对星间距离加速度噪声的影响。

The analytical investigation into the influence of gravitational perturbations on the TianQin constellation is incomplete. Existing models have neglected the satellite’s orbital eccentricity, a crucial factor for constellation stability [24, 25]. Moreover, relying solely on the leading-order lunar perturbation is insufficient to address the high-altitude TianQin orbits. These issues highlight the necessity for an analytical study to develop a more explicit and higher-precision model.
关于引力扰动对天琴星座影响的分析研究尚不完整。现有的模型忽略了卫星的轨道偏心率，这是星座稳定性的一个关键因素[24,25]。此外，仅仅依靠超前阶月球扰动不足以解决高空天琴轨道的问题。这些问题凸显了分析研究开发更明确和更高精度模型的必要性。

In the exploration of three-satellite constellations in heliocentric GW missions, such as LISA [26, 20, 27, 28, 29, 30, 18] and Taiji [17, 31], expressions for these three indicators have been derived and analyzed using either Keplerian orbits or perturbation solutions of satellite orbits. Valuable references are also found in geocentric satellite formation missions, including NASA’s four-satellite Magnetospheric Multiscale (MMS) mission [16], and extensive studies on third-body perturbations in general satellites (see [32] and references therein). Perturbation solutions for third-body effects can be derived by solving Lagrange’s planetary equations [33], where the perturbative potential depends on the orbital elements of both the satellite and perturbing bodies. To directly obtain solutions with instantaneous elements, perturbation methods [33, 34, 35, 36], especially the mean element method [34, 35], are utilized. This method employs a slowly precessing elliptical orbit as a reference, effectively reducing errors in analytical solutions.
在日心GW任务中对三星星座的探索中，如LISA [26， 20， 27， 28， 29， 30， 18] 和太极[17， 31]，已经推导出了这三个指标的表达式，并使用开普勒轨道或卫星轨道的扰动解进行了分析。在地心卫星编队任务中也发现了有价值的参考资料，包括NASA的四颗卫星磁层多尺度（MMS）任务[16]，以及对通用卫星中第三体扰动的广泛研究（见[32]和其中的参考文献）。第三体效应的扰动解可以通过求解拉格朗日行星方程 [33] 来推导出，其中扰动势取决于卫星和扰动体的轨道元件。为了直接获得瞬时单元的解，使用了微扰法 [33， 34， 35， 36]，尤其是均元法 [34， 35]。该方法采用慢进椭圆轨道作为参考，有效减少了解析解中的误差。

In this work, we will construct an analytical model for the TianQin constellation. To address its near-circular, high Earth orbits, we utilize singularity-free Lagrange equations while accounting for lunar, solar perturbations, and Earth’s $J_{2}$ perturbation. This model will then be used to analyze and optimize the three kinematic indicators. Additionally, to facilitate the perturbation-inclusive study, the unperturbed Keplerian orbits of TianQin satellites will also be presented.
在这项工作中，我们将为天琴座构建一个解析模型。为了解决其近圆形的高地球轨道问题，我们利用了无奇点的拉格朗日方程，同时考虑了月球、太阳的扰动和地球 $J_{2}$ 的扰动。然后，该模型将用于分析和优化三个运动指标。此外，为了促进包含微扰的研究，还将介绍天琴卫星的未扰动开普勒轨道。

The paper is organized as follows. In Sec. II, we introduce the Keplerian orbits of the satellites and present the design of the nominal equilateral-triangle constellation. The gravitational perturbations on the constellation are studied in Sec. III. In Sec. IV, we make the concluding remarks.
本文的组织结构如下。在第二部分，我们介绍了卫星的开普勒轨道，并介绍了名义等边三角形星座的设计。星座上的引力扰动在第三节中进行了研究。在第四节中，我们做了结束语。

II Fundamentals of stable TianQin constellation
二、稳定天琴座的基本原理

In this section, we describe the motion of TianQin satellites in the geocentric ecliptic coordinate system and present the orbit design of satellites for a stable equilateral-triangle constellation.
本节描述了天琴卫星在地心黄道坐标系中的运动，并介绍了稳定等边三角形星座的卫星轨道设计。

II.1 Keplerian orbits of satellites
II.1 卫星的开普勒轨道

Within the central gravitational field of the Earth, satellite moves in a Keplerian orbit, as illustrated in Fig. 1.
在地球的中心引力场内，卫星在开普勒轨道上移动，如图 1 所示。1.

Refer to caption — Figure 1: Depiction of the TianQin constellation in the geocentric ecliptic coordinate system. The ecliptic plane is spanned by the $x$ and $y$ axes, with the $x$ -axis directed toward the vernal equinox. The figure also illustrates the orbital coordinate system $\{X,Y,Z\}$ for SC1, where the $X$ -axis points toward the perigee of the satellite’s orbit, and the $Z$ -axis (not shown) is perpendicular to the orbital plane. The angles $i$ , $\Omega$ , $\omega$ , and $\nu$ denote the orbital inclination, longitude of ascending node, argument of perigee, and true anomaly, respectively. Specifically, $i=94.7^{\circ}$ and $\Omega=210.4^{\circ}$ are set to orient the TianQin detector plane toward the reference source, the white-dwarf binary RX J0806.3+1527.
图 1：天琴座在地心黄道坐标系中的描述。黄道平面由 $x$ 和 $y$ 轴跨越， $x$ 其中 -轴指向春分。该图还说明了 SC1 的轨道坐标系 $\{X,Y,Z\}$ ，其中 $X$ -轴指向卫星轨道的近地点， $Z$ 而 -轴（未显示）垂直于轨道平面。角度 $i$ 、 $\Omega$ 、 $\omega$ 和 $\nu$ 分别表示轨道倾角、升交点的经度、近地点的参数和真实异常。具体来说， $i=94.7^{\circ}$ 并将 $\Omega=210.4^{\circ}$ 天琴探测器平面定向到参考源，即白矮星双星 RX J0806.3+1527。

$\{X,Y,Z\}$ is the orbital right-handed coordinate system, with the origin at the Earth’s center of mass. The satellite’s orbital plane is same as the $X$ - $Y$ plane, where the $X$ -axis points toward the perigee. In this system, the satellite’s Cartesian coordinates $(X,Y,Z,\dot{X},\dot{Y},\dot{Z})$ can be denoted as [35, 33]:
$\{X,Y,Z\}$ 是轨道右手坐标系，原点位于地球质心。卫星的轨道平面与 $X$ - $Y$ 平面相同，其中 $X$ -轴指向近地点。在这个系统中，卫星的笛卡尔坐标 $(X,Y,Z,\dot{X},\dot{Y},\dot{Z})$ 可以表示为 [35， 33]：

\displaystyle\left\{\begin{array}[]{l}X=r\cos\nu=a(\cos E-e),\vskip 3.0pt plus% 1.0pt minus 1.0pt\\ Y=r\sin\nu=a\sqrt{1-e^{2}}\sin E,\vskip 3.0pt plus 1.0pt minus 1.0pt\\ Z=0,\vskip 3.0pt plus 1.0pt minus 1.0pt\\ \dot{X}=-\frac{\sqrt{\mu a}}{r}\sin E,\vskip 3.0pt plus 1.0pt minus 1.0pt\\ \dot{Y}=\frac{\sqrt{\mu a}}{r}\sqrt{1-e^{2}}\cos E,\vskip 3.0pt plus 1.0pt % minus 1.0pt\\ \dot{Z}=0,\end{array}\right.

(7)

with $r$ representing the geocentric radius, $\nu$ the true anomaly, $a$ the semimajor axis, $e$ the orbital eccentricity, $\mu=GM_{\text{e}}$ the Earth’s gravitational constant, and $E$ the eccentric anomaly. $E$ satisfies Kepler’s equation,
代表 $r$ 地心半径、 $\nu$ 真实异常、 $a$ 长半轴、 $e$ 轨道偏心率、 $\mu=GM_{\text{e}}$ 地球引力常数和 $E$ 偏心异常。 $E$ 满足 Kepler 方程，

\displaystyle E-e\sin E=M,

(8)

where $M$ denotes the mean anomaly. Specifically, $M$ is given by $M=n\,(t-t_{\text{p}})$ in the two-body problem, with the mean motion $n$ and the passing time of the perigee $t_{\text{p}}$ . Equation (8), which is a transcendental equation, can be solved iteratively, resulting in the following expression [21]:
其中 $M$ 表示均值距平。具体来说， $M$ 由 $M=n\,(t-t_{\text{p}})$ 在双体问题中给出，具有近地点的平均运动 $n$ 和通过时间 $t_{\text{p}}$ 。方程（8）是一个超越方程，可以迭代求解，得到以下表达式 [21]：

\displaystyle E=M+e\sin M+e^{2}\cos M\sin M+\mathcal{O}(e^{3}).

(9)

By substituting Eq. (9) into Eq. (7), one can obtain the explicit coordinates $(X,Y,Z,\dot{X},\dot{Y},\dot{Z})$ .
通过将方程（9）代入方程（7），可以得到显式坐标 $(X,Y,Z,\dot{X},\dot{Y},\dot{Z})$ 。

The orbital planes may not be identical for the three TianQin satellites. Thus, the geocentric ecliptic coordinate system $\{x,y,z\}$ is also employed, where the $x$ - $y$ plane is the ecliptic plane. The $x$ -axis points toward the vernal equinox, and the $z$ -axis is normal to the ecliptic plane. The coordinates $(x,y,z)$ and $(\dot{x},\dot{y},\dot{z})$ in this system can be obtained by $(X,Y,Z)$ and $(\dot{X},\dot{Y},\dot{Z})$ through the following transformation [35, 33]:
三颗天琴卫星的轨道平面可能并不相同。因此，还采用了地心黄道坐标系 $\{x,y,z\}$ ，其中 $x$ - $y$ 平面是黄道平面。 $x$ -轴指向春分， $z$ -轴垂直于黄道平面。在这个系统中，坐标 $(x,y,z)$ 和 $(\dot{x},\dot{y},\dot{z})$ 可以通过以下变换 [35， 33] 获得 $(X,Y,Z)$ $(\dot{X},\dot{Y},\dot{Z})$ ：

	$\displaystyle\begin{bmatrix}x\\ y\\ z\end{bmatrix}=R_{z}(-\Omega)R_{x}(-i)R_{z}(-\omega)\begin{bmatrix}X\\ Y\\ Z\end{bmatrix},$		(10)
	$\displaystyle\begin{bmatrix}\dot{x}\\ \dot{y}\\ \dot{z}\end{bmatrix}=R_{z}(-\Omega)R_{x}(-i)R_{z}(-\omega)\begin{bmatrix}\dot{% X}\\ \dot{Y}\\ \dot{Z}\end{bmatrix},$		(11)

where $\Omega$ , $i$ , and $\omega$ denote the satellite’s longitude of the ascending node, inclination, and argument of perigee, respectively. Additionally, $R_{z}(\gamma)$ and $R_{x}(\gamma)$ are the rotation matrices that rotate vectors by an angle $\gamma$ about the $z-$ or $x-$ axis,
其中 $\Omega$ ，、 $i$ 和 $\omega$ 分别表示卫星的升交点经度、倾角和近地点参数。此外， $R_{z}(\gamma)$ 和 $R_{x}(\gamma)$ 是将向量绕 $z-$ or $x-$ 轴旋转一定角度 $\gamma$ 的旋转矩阵，

	$\displaystyle R_{z}(\gamma)=\begin{bmatrix}\cos\gamma&\sin\gamma&0\\ -\sin\gamma&\cos\gamma&0\\ 0&0&1\end{bmatrix},$		(12)
	$\displaystyle R_{x}(\gamma)=\begin{bmatrix}1&0&0\\ 0&\cos\gamma&\sin\gamma\\ 0&-\sin\gamma&\cos\gamma\end{bmatrix}.$		(13)

Combining Eqs. (7) and (9)-(13), the position vector and velocity vector of SC $k$ ( $k=1$ , 2, 3), $\mathbf{r}_{k}=(x_{k},y_{k},z_{k})$ and $\mathbf{\dot{r}}_{k}=(\dot{x}_{k},\dot{y}_{k},\dot{z}_{k})$ , are given by
组合方程。（7）和（9）-（13），则 SC $k$ （ $k=1$ ， 2， 3） $\mathbf{r}_{k}=(x_{k},y_{k},z_{k})$ 和 $\mathbf{\dot{r}}_{k}=(\dot{x}_{k},\dot{y}_{k},\dot{z}_{k})$ 的位置矢量和速度矢量由下式给出

\displaystyle\left\{\begin{array}[]{l}\displaystyle x_{k}=a_{k}[\cos\Omega_{k}% \cos\lambda_{k}-\cos i_{k}\sin\Omega_{k}\sin\lambda_{k}\vskip 3.0pt plus 1.0pt% minus 1.0pt\\ \displaystyle\phantom{x_{k}=}+\frac{1}{2}e_{k}(\cos\Omega_{k}\,f_{\text{c1}k}-% \cos i_{k}\sin\Omega_{k}\,f_{\text{s1}k})]+\mathcal{O}(e^{2}_{k}),\vskip 3.0pt% plus 1.0pt minus 1.0pt\\ \displaystyle y_{k}=a_{k}[\sin\Omega_{k}\cos\lambda_{k}+\cos i_{k}\cos\Omega_{% k}\sin\lambda_{k}\vskip 3.0pt plus 1.0pt minus 1.0pt\\ \displaystyle\phantom{y_{k}=}+\frac{1}{2}e_{k}(\sin\Omega_{k}\,f_{\text{c1}k}+% \cos i_{k}\cos\Omega_{k}\,f_{\text{s1}k})]+\mathcal{O}(e^{2}_{k}),\vskip 3.0pt% plus 1.0pt minus 1.0pt\\ \displaystyle z_{k}=a_{k}(\sin i_{k}\sin\lambda_{k}+\frac{1}{2}e_{k}\sin i_{k}% \,f_{\text{s1}k})+\mathcal{O}(e^{2}_{k}),\vskip 3.0pt plus 1.0pt minus 1.0pt\\ \displaystyle\dot{x}_{k}=\frac{\sqrt{\mu}}{\sqrt{a_{k}}}[-\cos\Omega_{k}\sin% \lambda_{k}-\cos i_{k}\sin\Omega_{k}\cos\lambda_{k}\vskip 3.0pt plus 1.0pt % minus 1.0pt\\ \displaystyle\phantom{\dot{x}_{k}=}+e_{k}(-\cos\Omega_{k}\,f_{\text{s2}k}-\cos i% _{k}\sin\Omega_{k}\,f_{\text{c2}k})]+\mathcal{O}(e^{2}_{k}),\vskip 3.0pt plus % 1.0pt minus 1.0pt\\ \displaystyle\dot{y}_{k}=\frac{\sqrt{\mu}}{\sqrt{a_{k}}}[-\sin\Omega_{k}\sin% \lambda_{k}+\cos i_{k}\cos\Omega_{k}\cos\lambda_{k}\vskip 3.0pt plus 1.0pt % minus 1.0pt\\ \displaystyle\phantom{\dot{y}_{k}=}+e_{k}(-\sin\Omega_{k}\,f_{\text{s2}k}+\cos i% _{k}\cos\Omega_{k}\,f_{\text{c2}k})]+\mathcal{O}(e^{2}_{k}),\vskip 3.0pt plus % 1.0pt minus 1.0pt\\ \displaystyle\dot{z}_{k}=\frac{\sqrt{\mu}}{\sqrt{a_{k}}}(\sin i_{k}\cos\lambda% _{k}+e_{k}\sin i_{k}\,f_{\text{c2}k})+\mathcal{O}(e^{2}_{k}),\end{array}\right.

(24)

where 哪里

\displaystyle\lambda_{k}:=M_{k}+\omega_{k},

(25)

$f_{\text{c1}k}:=\cos(2\lambda_{k}-\omega_{k})-3\cos\omega_{k}$ , $f_{\text{s1}k}:=\sin(2\lambda_{k}-\omega_{k})-3\sin\omega_{k}$ , $f_{\text{c2}k}:=\cos(2\lambda_{k}-\omega_{k})$ , and $f_{\text{s2}k}:=\sin(2\lambda_{k}-\omega_{k})$ . Define $\sigma\in\{a,e,i,\Omega,\omega,\lambda\}$ , $\sigma_{k}(t)$ in Eq. (24) are straightforwardly determined in the two-body problem by
$f_{\text{c1}k}:=\cos(2\lambda_{k}-\omega_{k})-3\cos\omega_{k}$ $f_{\text{c2}k}:=\cos(2\lambda_{k}-\omega_{k})$ 、 $f_{\text{s1}k}:=\sin(2\lambda_{k}-\omega_{k})-3\sin\omega_{k}$ 和 $f_{\text{s2}k}:=\sin(2\lambda_{k}-\omega_{k})$ .定义 $\sigma\in\{a,e,i,\Omega,\omega,\lambda\}$ ， $\sigma_{k}(t)$ 在方程（24）中直接由下式确定

\displaystyle\sigma_{k}(t)=\sigma_{0k}+\delta_{\iota}n_{k}(t-t_{0}),\qquad% \text{for Keplerian orbit},

(26)

where 哪里

\sigma_{0k}:=\sigma_{k}(t_{0}),\qquad\delta_{\iota}=\left\{\begin{array}[]{ll}% 1,&\sigma=\lambda,\vskip 3.0pt plus 1.0pt minus 1.0pt\\ 0,&\sigma\neq\lambda.\end{array}\right.

(27)

Note that Eq. (24) remains valid even when considering gravitational perturbations, with the only change being the replacement of $\sigma_{k}(t)$ in Eq. (26) with the corresponding perturbation solution.
请注意，即使考虑引力扰动，方程（24）仍然有效，唯一的变化是用相应的扰动解替换了方程（26） $\sigma_{k}(t)$ 中的方程。

II.2 Orbit design of the TianQin constellation
二.2 天琴座的轨道设计

The TianQin constellation is composed of three satellites in geocentric orbits, forming a triangular configuration which continuously evolves in geometry over time. The closer the configuration change approaches an equilateral triangle, the more it aids in alleviating design constraints on measurement system instruments. Therefore, it is essential to find a constellation design with minimal variations.
天琴星座由三颗位于地心轨道上的卫星组成，形成一个三角形结构，随着时间的推移，几何结构不断演变。配置更改越接近等边三角形，就越有助于减轻测量系统仪器的设计限制。因此，找到变化最小的星座设计至关重要。

The constellation is considered more stable if it is closer to an equilateral triangle. There are three main kinematic indicators to characterize the stability, namely, the arm-length $L_{ij}$ , relative line-of-sight velocity between satellites $v_{ij}$ , and breathing angle $\alpha_{k}$ ,
如果星座更靠近等边三角形，则认为它更稳定。有三个主要的运动学指标来表征稳定性，即臂长 $L_{ij}$ 、卫星 $v_{ij}$ 之间的相对视线速度和呼吸角 $\alpha_{k}$ ，

$\displaystyle L_{ij}$	$\displaystyle=\left\|\mathbf{r}_{i}-\mathbf{r}_{j}\right\|,$	(28)
$\displaystyle v_{ij}$	$\displaystyle=\dot{L}_{ij},$	(29)
$\displaystyle\alpha_{k}$	$\displaystyle=\arccos\frac{L_{ki}^{2}+L_{kj}^{2}-L_{ij}^{2}}{2L_{ki}L_{kj}},$	(30)

where $i$ , $j$ , and $k$ take values 1, 2, or 3 and $i\neq j\neq k$ . Substituting Eq. (24) into Eqs. (28)-(30), one can obtain the explicit expressions, for these three kinematic indicators, with forms $L_{ij}(\sigma_{i}(t),\sigma_{j}(t))$ , $v_{ij}(\sigma_{i}(t),\sigma_{j}(t))$ , and $\alpha_{k}(\sigma_{k}(t),\sigma_{i}(t),\sigma_{j}(t))$ , respectively.
其中 $i$ ， $j$ ，和 $k$ 取值 1、2 或 3 和 $i\neq j\neq k$ 。将方程（24）代入方程。（28）-（30）中，对于这三个运动指示符，可以分别获得、和 $v_{ij}(\sigma_{i}(t),\sigma_{j}(t))$ $\alpha_{k}(\sigma_{k}(t),\sigma_{i}(t),\sigma_{j}(t))$ 、的 $L_{ij}(\sigma_{i}(t),\sigma_{j}(t))$ 显式表达式。

To maintain the constellation as an equilateral triangle, i.e. $L_{12}(t)=L_{13}(t)=L_{23}(t)$ , the orbits of the three satellites need to be purposefully designed. One intuitive orbit design involves circular orbits for the satellites in the point-mass gravitational field of the Earth:
为了保持星座为等边三角形 $L_{12}(t)=L_{13}(t)=L_{23}(t)$ ，即三颗卫星的轨道需要有目的地设计。一种直观的轨道设计涉及卫星在地球点质量引力场中的圆形轨道：

\displaystyle e_{1}(t)=e_{2}(t)=e_{3}(t)=0,

(31)

while ensuring that they share the same orbit size, lie in the same orbital plane, and are phased 120 degrees apart from each other:
同时确保它们具有相同的轨道大小，位于相同的轨道平面上，并且彼此相距 120 度：

\displaystyle\left\{\begin{array}[]{l}a_{1}(t)=a_{2}(t)=a_{3}(t),\vskip 3.0pt % plus 1.0pt minus 1.0pt\\ i_{1}(t)=i_{2}(t)=i_{3}(t),\vskip 3.0pt plus 1.0pt minus 1.0pt\\ \Omega_{1}(t)=\Omega_{2}(t)=\Omega_{3}(t),\vskip 3.0pt plus 1.0pt minus 1.0pt% \\ \lambda_{k}(t)=\frac{2\pi}{3}(k-1)+\lambda_{1}(t).\end{array}\right.

(36)

The above requirements on the inter-satellite parameters can be achieved in the two-body problem, if the initial orbital elements $\sigma_{0k}$ in Eq. (26) are set to
如果方程（26）中的初始轨道单元 $\sigma_{0k}$ 设置为

\displaystyle\sigma_{0k}=\sigma_{\text{o}}+\delta_{\iota}\frac{2\pi}{3}(k-1),

(37)

where the parameters with subscript “o” are the nominal ones of the TianQin constellation. For instance, these values can be chosen as $a_{\text{o}}=10^{5}$ km, $e_{\text{o}}=0$ , and $i_{\text{o}}=94.7^{\circ}$ , $\Omega_{\text{o}}=210.4^{\circ}$ , respectively establishing the orbit size and orienting the orbital plane perpendicular to J0806 [7, 24]. The initial value $\lambda_{\text{o}}$ associated with the orbit phase is typically selected to be any value within the range of $0^{\circ}$ to $120^{\circ}$ , or it may be specifically designated to avoid Moon eclipses [37].
其中带下标 “o” 的参数是天琴座的名义参数。例如，这些值可以选择为 $a_{\text{o}}=10^{5}$ km、 $e_{\text{o}}=0$ 、和 $i_{\text{o}}=94.7^{\circ}$ 、 $\Omega_{\text{o}}=210.4^{\circ}$ ，分别建立轨道大小和垂直于 J0806 的轨道平面 [7， 24]。与轨道相位相关的初始值 $\lambda_{\text{o}}$ 通常选择为 $0^{\circ}$ 到 $120^{\circ}$ 范围内的任何值，或者可以专门指定以避免月食 [37]。

To analyze additional nominal orbit design allowing for $e\neq 0$ and quantify the impact of eccentricity on the three indicators, the constraint specified by Eq. (31) is relaxed. Subsequently, employing only Eq. (36) or Eq. (37) (for $\sigma\in\{a,i,\Omega,\lambda\}$ ), the variations of these indicators in the two-body problem, up to the first order of $e$ , can be expressed as
为了分析允许 $e\neq 0$ 和量化偏心率对三个指标影响的额外标称轨道设计，式（31）规定的约束是放宽的。随后，仅使用方程（36）或方程（37）（for $\sigma\in\{a,i,\Omega,\lambda\}$ ），这些指标在两体问题中的变化，直到 $e$ 的一阶，可以表示为

	$\displaystyle L_{ij}^{\text{kepl}}(t)=$	$\displaystyle~{}\sqrt{3}a_{\text{o}}+\frac{\sqrt{7}}{2}a_{\text{o}}[-e_{i}\sin% \left(M_{i}+\beta\right)$
		$\displaystyle+e_{j}\sin\left(M_{j}-\beta\right)],$		(38)

	$\displaystyle v_{ij}^{\text{kepl}}(t)=$	$\displaystyle~{}0+\frac{\sqrt{7}}{2}a_{\text{o}}n_{\text{o}}[-e_{i}\cos\left(M% _{i}+\beta\right)$
		$\displaystyle+e_{j}\cos\left(M_{j}-\beta\right)],$		(39)

$\displaystyle\alpha_{k}^{\text{kepl}}(t)=$	$\displaystyle~{}\frac{\pi}{3}+\frac{\sqrt{7}}{3}e_{k}\cos M_{k}\sin\beta$
	$\displaystyle-\frac{\sqrt{7}}{6}e_{i}[\sin\left(M_{i}-\beta\right)+2\sin\left(% M_{i}+\beta\right)]$
	$\displaystyle+\frac{\sqrt{7}}{6}e_{j}[2\sin\left(M_{j}-\beta\right)+\sin\left(% M_{j}+\beta\right)],$	(40)

where $M_{i}(t)=\frac{2\pi}{3}(i-1)+\lambda_{1}(t)-\omega_{i}$ and $\beta:=\arccos\frac{2}{\sqrt{7}}$ , with the indices $i$ , $j$ , and $k$ using cyclic indexing ( $i,\,j,\,k=1\to 2\to 3\to 1$ ). If we further set $e_{1}=e_{2}=e_{3}\equiv e_{\text{o}}\neq 0$ and $\omega_{1}=\omega_{2}=\omega_{3}$ , then it follows that
其中 $M_{i}(t)=\frac{2\pi}{3}(i-1)+\lambda_{1}(t)-\omega_{i}$ 和 $\beta:=\arccos\frac{2}{\sqrt{7}}$ ，带有索引 $i$ ， $j$ ，并使用 $k$ 循环索引（ $i,\,j,\,k=1\to 2\to 3\to 1$ ）。如果我们进一步设置 $e_{1}=e_{2}=e_{3}\equiv e_{\text{o}}\neq 0$ 和 $\omega_{1}=\omega_{2}=\omega_{3}$ ，则

	$\displaystyle L_{ij}^{\text{kepl-A}}(t)=\sqrt{3}a_{\text{o}}-\frac{\sqrt{3}}{2% }a_{\text{o}}e_{\text{o}}\cos M_{k},$		(41)
	$\displaystyle v_{ij}^{\text{kepl-A}}(t)=0+\frac{\sqrt{3}}{2}a_{\text{o}}n_{% \text{o}}e_{\text{o}}\sin M_{k},$		(42)
	$\displaystyle\alpha_{k}^{\text{kepl-A}}(t)=\frac{\pi}{3}-\frac{\sqrt{3}}{2}e_{% \text{o}}\cos M_{k},$		(43)

where $M_{k}(t)=\frac{2\pi}{3}(k-1)+M_{1}(t)$ . Equations (41)-(43) indicate that, at the zeroth order of $e$ , the three TianQin satellites can form a constant equilateral triangle. However, when accounting for eccentricity, as observed in perturbed orbits, the constellation’s evolution deviates from the ideal equilateral triangle, exhibiting periodic variations.
其中 $M_{k}(t)=\frac{2\pi}{3}(k-1)+M_{1}(t)$ . 方程（41）-（43）表明，在的零阶 $e$ 处，三颗天琴卫星可以形成一个恒定的等边三角形。然而，当考虑到在扰动轨道上观察到的偏心率时，星座的演化偏离了理想的等边三角形，表现出周期性变化。

The close-to-circular orbits, as inspired by Eqs. (31) and (36) for Keplerian orbits, are currently employed in TianQin orbit studies (see, e.g., [21, 24, 13, 38, 39, 40, 41, 42]). It is worth noting that, to obtain the nominal equilateral triangle configuration, there is another option: elliptical frozen orbits. From Eqs. (38)-(40), if $e_{1}=e_{2}=e_{3}\equiv e_{\text{o}}\neq 0$ , and $M_{1}=M_{2}=M_{3}$ , namely $\omega_{k}=\omega_{1}+\frac{2\pi}{3}(k-1)$ , then $L_{ij}^{\text{kepl-B}}(t)\equiv\sqrt{3}a_{\text{o}}-\sqrt{3}a_{\text{o}}e_{% \text{o}}\cos E_{1}(t)$ , representing an equilateral-triangle constellation with three arm-lengths that vary synchronously. Preliminary numerical simulation results show that the long-term stability of the constellation, based on this design, is not as favorable as that of the close-to-circular orbits. Furthermore, the impact of this design on other aspects of the mission, such as point-ahead angle variations associated with the finite speed of light, requires further assessment. In this paper, we focus on the study of a three-satellite constellation with close-to-circular orbits.
受 Eqs 启发的接近圆形的轨道。（31）和（36）目前用于天琴轨道研究（参见 [21， 24， 13， 38， 39， 40， 41， 42]）。值得注意的是，要获得名义上的等边三角形构型，还有另一种选择：椭圆冻结轨道。来自 Eqs.（38）-（40），如果 $e_{1}=e_{2}=e_{3}\equiv e_{\text{o}}\neq 0$ ，和 $M_{1}=M_{2}=M_{3}$ ，即 $\omega_{k}=\omega_{1}+\frac{2\pi}{3}(k-1)$ ，则 $L_{ij}^{\text{kepl-B}}(t)\equiv\sqrt{3}a_{\text{o}}-\sqrt{3}a_{\text{o}}e_{% \text{o}}\cos E_{1}(t)$ ，代表一个等边三角形星座，具有三个同步变化的臂长。初步数值模拟结果表明，基于此设计的星座的长期稳定性不如接近圆形的轨道。此外，这种设计对任务其他方面的影响，例如与有限光速相关的前方角度变化，需要进一步评估。在本文中，我们专注于研究具有接近圆形轨道的三卫星星座。

III Effects of lunisolar perturbations on TianQin constellation
三、阴阳扰动对天秦星座的影响

The TianQin constellation is subject not only to the central gravitational attraction but also to gravitational perturbations. These perturbations can distort the carefully designed equilateral-triangle configuration. To gain a more accurate understanding of the TianQin constellation’s variations, it is crucial to account for these gravitational perturbations.
天琴座不仅受中心引力的影响，还受引力扰动的影响。这些扰动会扭曲精心设计的等边三角形构型。为了更准确地了解天琴座的变化，考虑这些引力扰动至关重要。

The primary perturbations originate from the Moon and the Sun, with magnitudes of approximately $4\times 10^{-4}$ and $2\times 10^{-4}$ , respectively [33]. In this section, we collectively address the point mass effects of these two perturbing bodies. Furthermore, we also incorporate the secular perturbation arising from the third most significant perturbation, Earth’s $J_{2}$ perturbation, which has a magnitude of $6\times 10^{-6}$ . Other perturbations, e.g., the higher-degree non-spherical gravity fields of the Earth, have a minor impact on satellite positions and constellation stability. As illustrated in Fig. 2, these perturbations lead to deviations of approximately 3.3 km in satellite positions, $\pm 2.2$ km in arm-lengths, $\pm 0.0020$ m/s in relative velocities, and $\pm 0.0012^{\circ}$ in breathing angles over a 5-year period. In contrast to Refs. [21, 22, 23, 36, 16, 32], the perturbation solutions developed in this study offer explicit expressions with improved precision, enabling a more precise description of the distinctive $10^{5}$ km-radius orbits relevant to space-based GW detection.
主要扰动起源于月球和太阳，震等分别约为 $4\times 10^{-4}$ 和 $2\times 10^{-4}$ [33]。在本节中，我们将共同讨论这两个扰动体的点质量效应。此外，我们还考虑了由第三个最显著的扰动，即地球 $J_{2}$ 的扰动引起的长期扰动，其震级为 $6\times 10^{-6}$ 。其他扰动，例如地球的更高程度的非球形重力场，对卫星位置和星座稳定性影响较小。如图所示。2，这些扰动导致卫星位置偏差约为 3.3 km，臂长偏差约为 km， $\pm 2.2$ $\pm 0.0020$ 相对速度偏差为 m/s， $\pm 0.0012^{\circ}$ 呼吸角偏差为 5 年。与 Refs 相反。 [21， 22， 23， 36， 16， 32]，本研究开发的扰动解决方案提供了明确的表达式，精度更高，能够更精确地描述与天基GW探测相关的不同 $10^{5}$ km半径轨道。

Figure 2: Time evolution of deviations between the simplified (SNM) and high-precision (HPNM) numerical models for satellite positions and three indicators. SNM incorporates the point-mass gravity fields of the Earth, Moon, and Sun, as well as the Earth’s

J_{2}

. HPNM additionally considers the higher-degree non-spherical gravity fields of the Earth and the point-mass gravity fields of other planets (for further details, see Appendix B.3). In these plots, red corresponds to SC1,

v_{23}

L_{23}

, or

\alpha_{1}

; green represents SC2,

v_{31}

L_{31}

, or

\alpha_{2}

; and blue indicates SC3,

v_{12}

L_{12}

, or

\alpha_{3}

. The initial orbital elements used are the same as those presented in Table 6.
图 2：卫星位置和三个指标的简化（SNM）和高精度（HPNM）数值模型之间偏差的时间演变。SNM 包含地球、月球和太阳的点质量重力场，以及地球的

J_{2}

.HPNM 还考虑了地球的更高阶非球形重力场和其他行星的点质量重力场（有关详细信息，请参阅附录 B.3）。在这些图中，红色对应于 SC1、

v_{23}

、

L_{23}

或

\alpha_{1}

;绿色代表 SC2、

v_{31}

、

L_{31}

或

\alpha_{2}

;蓝色表示 SC3、

v_{12}

、

L_{12}

或

\alpha_{3}

。使用的初始轨道元件与表 6 中所示的相同。

III.1 Dynamic model
三.1 动态模型

The gravitational potential $U$ acting on a satellite can be expressed as
作用在卫星上的引力势 $U$ 可以表示为

\displaystyle U=U_{0}+\mathcal{R},\qquad U_{0}=\frac{\mu}{r},

(44)

where $U_{0}$ is the gravitational potential of a pointlike Earth, and $\mathcal{R}$ represents a perturbative potential describing the satellite’s perturbed motion. Under the influence of $\mathcal{R}$ , the evolution of the satellite’s orbital elements is governed by Lagrange’s planetary equations [43],
其中 $U_{0}$ 是点状地球的引力势，表示 $\mathcal{R}$ 描述卫星扰动运动的扰动势。在 $\mathcal{R}$ 的影响下，卫星轨道元件的演化受拉格朗日行星方程 [43] 的控制，

$\displaystyle\frac{da}{dt}=$	$\displaystyle~{}\frac{2}{na}\frac{\partial\mathcal{R}}{\partial\lambda},$	(45)
$\displaystyle\frac{di}{dt}=$	$\displaystyle~{}\frac{1}{na^{2}\check{e}\sin i}\left[\cos i\left(\xi\frac{% \partial\mathcal{R}}{\partial\eta}-\eta\frac{\partial\mathcal{R}}{\partial\xi}% +\frac{\partial\mathcal{R}}{\partial\lambda}\right)\right.$
	$\displaystyle\left.-\frac{\partial\mathcal{R}}{\partial\Omega}\right],$	(46)
$\displaystyle\frac{d\Omega}{dt}=$	$\displaystyle~{}\frac{1}{na^{2}\check{e}\sin i}\frac{\partial\mathcal{R}}{% \partial i},$	(47)
$\displaystyle\frac{d\xi}{dt}=$	$\displaystyle-\frac{\check{e}}{na^{2}}\frac{\partial\mathcal{R}}{\partial\eta}% -\xi\frac{\check{e}}{na^{2}(1+\check{e})}\frac{\partial\mathcal{R}}{\partial% \lambda}+\eta\cos i\frac{d\Omega}{dt},$	(48)
$\displaystyle\frac{d\eta}{dt}=$	$\displaystyle~{}\frac{\check{e}}{na^{2}}\frac{\partial\mathcal{R}}{\partial\xi% }-\eta\frac{\check{e}}{na^{2}(1+\check{e})}\frac{\partial\mathcal{R}}{\partial% \lambda}-\xi\cos i\frac{d\Omega}{dt},$	(49)
$\displaystyle\frac{d\lambda}{dt}=$	$\displaystyle~{}n-\frac{2}{na}\frac{\partial\mathcal{R}}{\partial a}+\frac{% \check{e}}{na^{2}(1+\check{e})}\left(\xi\frac{\partial\mathcal{R}}{\partial\xi% }+\eta\frac{\partial\mathcal{R}}{\partial\eta}\right)$
	$\displaystyle-\cos i\frac{d\Omega}{dt},$	(50)

where $\check{e}:=\sqrt{1-e^{2}}$ , and new variables $\xi$ and $\eta$ are introduced,
其中 $\check{e}:=\sqrt{1-e^{2}}$ ，和新变量 $\xi$ 和 $\eta$ 被引入，

\displaystyle\xi:=e\cos\omega,\qquad\eta:=e\sin\omega,

(51)

to avoid the singularity at $e=0$ . When $\mathcal{R}=0$ , the solutions to Eqs. (45)-(50) revert to the Keplerian case discussed in Sec. II.1.
为了避免处 $e=0$ 的奇点。当 $\mathcal{R}=0$ ，方程的解程.（45）-（50）回到第 II.1 节讨论的 Keplerian 案。

For TianQin orbits with an orbital radius of $10^{5}$ km, the perturbative potential $\mathcal{R}$ predominantly encompasses the perturbation effects arising from the Sun, Moon, and Earth’s $J_{2}$ term, as expressed in the following expressions:
对于轨道半径为 $10^{5}$ km 的天琴轨道，扰动势 $\mathcal{R}$ 主要包括由太阳、月亮和地球 $J_{2}$ 项引起的扰动效应，如以下表达式所示：

\displaystyle\mathcal{R}=\mathcal{R}_{\text{s}}+\mathcal{R}_{\text{m}}+% \mathcal{R}_{J_{2}}

(52)

with [33]
和 [33]

	$\displaystyle\mathcal{R}_{\text{s}}=\frac{\mu_{2}r^{2}}{r_{2}^{3}}\frac{3\cos^% {2}\psi_{2}-1}{2},$		(53)
	$\displaystyle\mathcal{R}_{\text{m}}=\frac{\mu_{3}}{r_{3}}\sum_{N=2}^{\mathcal{% N}}\left(\frac{r}{r_{3}}\right)^{N}P_{N}(\cos\psi_{3}),$		(54)
	$\displaystyle\mathcal{R}_{\!J_{2}}=-\frac{\mu R_{\text{e}}^{2}}{r^{3}}J_{2}% \frac{3\sin^{2}\varphi-1}{2},$		(55)

where $\mu_{2}=GM_{\text{s}}$ and $\mu_{3}=GM_{\text{m}}$ are the gravitational constants of the Sun and the Moon, respectively. $r_{2}$ and $r_{3}$ denote the geocentric distances of the Sun and the Moon. Moreover, $P_{N}(x)$ is the Legendre polynomial of degree $N$ , with $\mathcal{N}=6$ signifying the truncation degree. The derivation of Eq. (54) is presented in Appendix B.1, suggesting that employing Legendre polynomial expansions is more advantageous than the original square root form (Eq. (117)) for solving Lagrange’s equations. Additionally, $R_{\text{e}}$ stands for the equatorial radius of the Earth, and $J_{2}$ represents the second zonal harmonic coefficient. Furthermore, $\psi_{2}$ is the angular separation of the Sun and the satellite as observed from the Earth’s center,
其中 $\mu_{2}=GM_{\text{s}}$ 和 $\mu_{3}=GM_{\text{m}}$ 分别是太阳和月亮的引力常数。 $r_{2}$ 并 $r_{3}$ 表示太阳和月亮的地心距离。此外， $P_{N}(x)$ 是度 $N$ 的勒让德多项式，表示 $\mathcal{N}=6$ 截断度。方程（54）的推导见附录 B.1，表明采用勒让德多项式展开比原始平方根形式（方程（117））更有利于求解拉格朗日方程。此外， $R_{\text{e}}$ 表示地球的赤道半径，并 $J_{2}$ 表示第二纬向谐波系数。此外， $\psi_{2}$ 是从地心观测到的太阳和卫星的角度分离，

\displaystyle\cos\psi_{2}=\mathbf{\hat{r}_{2}}\cdot\mathbf{\hat{r}},

(56)

where $\mathbf{\hat{r}_{2}}$ and $\mathbf{\hat{r}}$ denote the unit position vectors of the Sun and the satellite, respectively,
其中 $\mathbf{\hat{r}_{2}}$ 和 $\mathbf{\hat{r}}$ 分别表示太阳和卫星的单位位置向量，

\displaystyle\mathbf{\hat{r}_{2}}=\begin{bmatrix}\cos\text{u}_{2}\\ \sin\text{u}_{2}\\ 0\end{bmatrix},\quad\mathbf{\hat{r}}=R_{z}(-\Omega)R_{x}(-i)R_{z}(-\omega)% \begin{bmatrix}\cos\nu\\ \sin\nu\\ 0\end{bmatrix},

(57)

with $\text{u}_{2}:=\Omega_{2}+\omega_{2}+\nu_{2}$ representing the Sun’s ecliptic longitude. Similarly, $\psi_{3}$ is given by
，代表 $\text{u}_{2}:=\Omega_{2}+\omega_{2}+\nu_{2}$ 太阳的黄道经度。同样， $\psi_{3}$ 由

\displaystyle\cos\psi_{3}=\mathbf{\hat{r}_{3}}\cdot\mathbf{\hat{r}},

(58)

with $\mathbf{\hat{r}_{3}}$ being the Moon’s unit position vector,
是 $\mathbf{\hat{r}_{3}}$ 月球的单位位置向量，

\displaystyle\mathbf{\hat{r}_{3}}=R_{z}(-\Omega_{3})R_{x}(-i_{3})R_{z}(-\text{% u}_{3})\begin{bmatrix}1\\ 0\\ 0\end{bmatrix},

(59)

where $\Omega_{3}$ , $i_{3}$ , and $\text{u}_{3}:=\omega_{3}+\nu_{3}$ correspond to the Moon’s longitude of ascending node, inclination, and latitude argument, respectively. $\varphi$ of Eq. (55) signifies the satellite’s geocentric latitude in the Earth-fixed coordinate system,
其中 $\Omega_{3}$ ， $i_{3}$ ，和 $\text{u}_{3}:=\omega_{3}+\nu_{3}$ 分别对应于月球的升交点经度、倾角和纬度参数。 $\varphi$ 方程（55）表示卫星在地球固定坐标系中的地心纬度，

\displaystyle\sin\varphi=\sin i\sin(\omega+\nu).

(60)

Substituting Eqs. (56)-(60) into (53)-(55), $\mathcal{R}$ is formulated as a function of the satellite’s orbital elements and those of the Sun and the Moon, $\mathcal{R}(a,e,i,\Omega,\omega,\nu;\sigma_{23})$ , where $\sigma_{23}=\{r_{2},\text{u}_{2},i_{3},\Omega_{3},r_{3},\text{u}_{3}\}$ . Consequently, the singularity-free form of the potential, $\mathcal{R}(a,\xi,\eta,i,\Omega,\lambda;\sigma_{23})$ , can be obtained. Moreover, it is important to note that Eq. (60) is formulated in the equatorial coordinate system. For a unified description of the influence of all three perturbations on satellite orbits, including both solar and lunar perturbations, the transformation to the ecliptic coordinate system must be taken into account (see Appendix C.3 for more details).
代入 eqs.（56）-（60）转换为（53）-（55）， $\mathcal{R}$ 表示为卫星轨道元件与太阳和月亮轨道元件的函数， $\mathcal{R}(a,e,i,\Omega,\omega,\nu;\sigma_{23})$ 其中 $\sigma_{23}=\{r_{2},\text{u}_{2},i_{3},\Omega_{3},r_{3},\text{u}_{3}\}$ 。因此，可以获得势能 $\mathcal{R}(a,\xi,\eta,i,\Omega,\lambda;\sigma_{23})$ 的无奇点形式。此外，需要注意的是，方程（60）是在赤道坐标系中表示的。要统一描述所有三种扰动对卫星轨道的影响，包括太阳和月球的扰动，必须考虑到黄道坐标系的变换（更多细节见附录 C.3）。

III.2 Motion of the Sun and Moon
三.2 太阳和月亮的运动

To solve the Lagrange equations, the coordinates of the Sun and the Moon, relative to $r_{2}$ , $\text{u}_{2}$ , $i_{3}$ , $\Omega_{3}$ , $r_{3}$ , and $\text{u}_{3}$ , are required. While high-precision coordinates are available through numerical integration, such as the JPL ephemerides DE430 [44], they are less suitable for analytical purposes. References like [35, 45] offer analytical formulas with reduced precision, providing geocentric solar coordinates based on a simplified, unperturbed motion of the Earth around the Sun and expressed using appropriate mean orbital elements. In contrast, the Moon’s motion, influenced by strong solar and terrestrial perturbations, is described through linear terms corresponding to its long-term precessing elliptical orbit and numerous trigonometric terms capturing periodic variations.
要求解拉格朗日方程，需要太阳和月亮相对于 $r_{2}$ 、 $\text{u}_{2}$ $i_{3}$ $\Omega_{3}$ $r_{3}$ 和 $\text{u}_{3}$ 的坐标。虽然可以通过数值积分获得高精度坐标，例如 JPL 星历仪 DE430 [44]，但它们不太适合用于分析目的。像 [35， 45] 这样的参考文献提供了精度降低的解析公式，根据地球围绕太阳的简化、不受干扰的运动提供地心太阳坐标，并使用适当的平均轨道元素表示。相比之下，月球的运动受强烈的太阳和地球扰动的影响，通过对应于其长期进动椭圆轨道的线性项和捕捉周期性变化的众多三角项来描述。

Generally, higher precision in these analytical coordinates results in more complex expressions, rendering the analytical solutions of the Lagrange equations challenging. In this study, for a balance between solvability and precision, essential components of these coordinates are retained, and fitting is applied using JPL ephemerides data [46] from around 2035 to reduce discrepancies in the positions of the Sun and the Moon.
通常，这些解析坐标的精度越高，表达式就越复杂，这使得拉格朗日方程的解析解具有挑战性。在这项研究中，为了在可解性和精度之间取得平衡，保留了这些坐标的基本组成部分，并使用2035年左右的JPL星历数据[46]进行拟合，以减少太阳和月亮位置的差异。

The apparent motion of the Sun around the Earth is approximated as a circular orbit on the ecliptic plane with a one-sidereal-year period,
太阳绕地球的视运动近似为黄道平面上的圆形轨道，周期为单侧实年，

	$\displaystyle r_{2}(t)\simeq\overline{r}_{2},$		(61)
	$\displaystyle\text{u}_{2}(t)\simeq u_{2}(t)=n_{2}\,t+u_{2_{0}},$		(62)

where $\overline{r}_{2}$ is the mean Sun-Earth distance, and $u_{2}$ represents the mean longitude of the Sun. Moreover, $n_{2}=2\pi/(\text{one sidereal year})$ is the mean motion, and $u_{2_{0}}$ denotes the initial phase of the Sun’s orbit. The specific parameter values can be found in Table 1. The Moon’s orbit is considered as an inclined and elliptical precessing orbit,
其中 $\overline{r}_{2}$ 是太阳与地球的平均距离， $u_{2}$ 代表太阳的平均经度。此外， $n_{2}=2\pi/(\text{one sidereal year})$ 是平均运动，表示 $u_{2_{0}}$ 太阳轨道的初始阶段。具体参数值见表1。月球的轨道被认为是一个倾斜和椭圆形的前部轨道，

	$\displaystyle i_{3}(t)\simeq\overline{i}_{3},$		(63)
	$\displaystyle\Omega_{3}(t)\simeq\overline{\Omega}_{3}(t),$		(64)
	$\displaystyle r_{3}(t)\simeq\overline{r}_{3}+r_{3}^{A}\cos\overline{M}_{3}(t),$		(65)
	$\displaystyle\text{u}_{3}(t)\simeq u_{3}(t)+u_{3}^{A}\sin\overline{M}_{3}(t),$		(66)

with 跟

	$\displaystyle\overline{\Omega}_{3}(t)=n_{\Omega_{3}}t+\Omega_{3_{0}},$		(67)
	$\displaystyle\overline{M}_{3}(t)=n_{\!M_{3}}t+M_{3_{0}},$		(68)
	$\displaystyle u_{3}(t)=n_{3}\,t+u_{3_{0}},$		(69)

where $\overline{i}_{3}$ is the mean inclination of the Moon’s orbit, $\overline{\Omega}_{3}$ represents the secular variation in the Moon’s longitude of ascending node, $\overline{r}_{3}$ denotes the mean Earth-Moon distance, and $u_{3}$ signifies the secular variation in the argument of latitude. Equations (65) and (66) include trigonometric corrections, related to the Moon’s mean anomaly $\overline{M}_{3}$ , aimed at more accurately describing the Moon’s motion in the radial and transverse directions. The periods of variation for $\overline{\Omega}_{3}$ , $\overline{M}_{3}$ , and $u_{3}$ are approximately 18.6 years, 27.55 days (anomalistic month), and 27.21 days (draconic month), respectively.
其中 $\overline{i}_{3}$ 是月球轨道的平均倾角， $\overline{\Omega}_{3}$ 表示月球升交点经度的长期变化， $\overline{r}_{3}$ 表示平均地月距离，并 $u_{3}$ 表示纬度论证中的长期变化。方程（65）和（66）包括与月球平均距平相关的三角校正 $\overline{M}_{3}$ ，旨在更准确地描述月球在径向和横向上的运动。 $\overline{\Omega}_{3}$ 、 $\overline{M}_{3}$ 和 $u_{3}$ 的变化周期分别约为 18.6 年、27.55 天（异常月）和 27.21 天（龙月）。

Table 1: Parameter settings for the motion of the Sun and Moon. The subscript “0” indicates values taken at the epoch 1 January, 2034, 00:00:00 UTC.
表 1：太阳和月亮运动的参数设置。下标“0”表示在 2034 年 1 月 1 日 00：00：00 UTC 纪元获取的值。

Symbols 符号	Parameters 参数	Values 值
$\overline{r}_{2}$	Mean Sun-Earth distance 平均日地距离	$1.496\,191\times 10^{8}\,\text{km}$
$n_{2}$	Mean motion of the Earth 地球的平均运动	$2\pi/365.2564\,\text{days}$
$u_{2_{0}}$	Sun initial phase Sun 初始阶段	$280.251^{\circ}$
$\overline{i}_{3}$	Mean lunar orbit inclination 平均月球轨道倾角	$5.162^{\circ}$
$n_{\Omega_{3}}$	Rate of change of $\overline{\Omega}_{3}$ 变化 $\overline{\Omega}_{3}$ 率	$2\pi/18.6\,\text{years}$
$\Omega_{3_{0}}$	Initial phase of $\overline{\Omega}_{3}$ 的起步阶段 $\overline{\Omega}_{3}$	$186.988^{\circ}$
$\overline{r}_{3}$	Mean Earth-Moon distance 平均地月距离	$384\,151\,\text{km}$
$r_{3}^{A}$	Amplitude of correction term 校正项的振幅	$-20\,905\,\text{km}$
$n_{\!M_{3}}$	Rate of change of $\overline{M}_{3}$ 变化 $\overline{M}_{3}$ 率	$2\pi/27.55\,\text{days}$
$M_{3_{0}}$	Initial phase of $\overline{M}_{3}$ 的初始阶段 $\overline{M}_{3}$	$22.578^{\circ}$
$n_{3}$	Rate of change of $u_{3}$ 变化 $u_{3}$ 率	$2\pi/27.21\,\text{days}$
$u_{3_{0}}$	Initial phase of $u_{3}$ 的初始阶段 $u_{3}$	$221.970^{\circ}$
$u_{3}^{A}$	Amplitude of correction term 校正项的振幅	$6.289^{\circ}$

III.3 Lunisolar perturbations on the TianQin satellites
三.3 天琴卫星上的阴阳扰动

Let $\sigma=\begin{bmatrix}a&i&\Omega&\xi&\eta&\lambda\end{bmatrix}^{T}$ ; then, the Lagrange perturbation equations (45)-(50) can be reformulated as
让 $\sigma=\begin{bmatrix}a&i&\Omega&\xi&\eta&\lambda\end{bmatrix}^{T}$ ;那么，拉格朗日扰动方程（45）-（50）可以重新表述为

\displaystyle\frac{d\sigma}{dt}=f_{0}(a)+f_{1}(\sigma,t,\varepsilon),

(70)

where the functions $f_{0}$ and $f_{1}$ are both 6-dimensional vector functions,
其中函数 $f_{0}$ 和 $f_{1}$ 都是 6 维向量函数

	$\displaystyle f_{0}(a)=\delta_{\iota}n,\qquad\delta_{\iota}=\begin{bmatrix}0&0% &0&0&0&1\end{bmatrix}^{T},$		(71)
	$\displaystyle(f_{1})_{\zeta}=\mathcal{O}(\varepsilon),\qquad\zeta=1,2,\cdots,6,$		(72)

and $\varepsilon\ll 1$ is a small parameter. Since the perturbing forces are significantly weaker than the Earth’s central gravitational attraction, the solution of Eq. (70) is assumed to be
and $\varepsilon\ll 1$ 是一个小参数。由于扰动力明显弱于地球的中心引力，因此假设方程（70）的解为

\displaystyle\sigma(t)=\sigma^{(0)}(t)+\sigma^{(1)}(t).

(73)

Here, $\sigma^{(0)}(t)$ represents the unperturbed Keplerian orbit (as shown in Eq. (26)),
这里， $\sigma^{(0)}(t)$ 表示未受扰动的开普勒轨道（如方程（26）所示），

\displaystyle\sigma^{(0)}(t)=\sigma_{0}+\delta_{\iota}n(t-t_{0}),

(74)

and $\sigma^{(1)}(t)$ is the perturbation solution,
并且是 $\sigma^{(1)}(t)$ 扰动解，

\displaystyle\sigma^{(1)}(t)=\bigl{(}\sigma_{c}(t)+\sigma_{l[c]}(t)+\sigma_{l[% l]}(t)+\sigma_{s}(t)\bigr{)}\!\large{|}_{t_{0}}^{t}.

(75)

In Eq. (75), $\sigma^{(1)}(t)$ is decomposed into four parts, distinguished by the unique time scales of orbital variations induced by perturbations: the secular term $\sigma_{c}$ , special long-period term $\sigma_{l[c]}$ , general long-period term $\sigma_{l[l]}$ , and short-period term $\sigma_{s}$ . $\sigma_{c}$ signifies the linear change over time, while $\sigma_{l[c]}$ , $\sigma_{l[l]}$ , and $\sigma_{s}$ are associated with periodic variations. These variations are linked to, for instance, $\Omega_{3}$ with an 18.6-year period, $u_{3}$ with a 27.21-day period, and $\lambda$ with a 3.64-day period. The explicit expression for $\sigma^{(1)}(t)$ can be derived using perturbation methods to solve Eq. (70).
在方程（75）中， $\sigma^{(1)}(t)$ 被分解为四个部分，其特点是由扰动引起的轨道变化的独特时间尺度：长期项 $\sigma_{c}$ 、特殊长周期项 $\sigma_{l[c]}$ 、一般长期项 $\sigma_{l[l]}$ 和短期项 $\sigma_{s}$ 。 $\sigma_{c}$ 表示随时间的线性变化，而 $\sigma_{l[c]}$ 、 $\sigma_{l[l]}$ 和与 $\sigma_{s}$ 周期性变化相关联。例如， $\Omega_{3}$ 这些变化与18.6年、27.21天 $\lambda$ 和3.64天的周期 $u_{3}$ 有关。可以使用扰动方法求解方程（70）来推导的 $\sigma^{(1)}(t)$ 显式表达式。

Table 2: The components of the analytical solution

\sigma(t)

for the TianQin orbit describing the perturbing effects from the Sun, Moon, and Earth’s

J_{2}

. The symbol “

-

” denotes that terms do not appear separately due to the joint effects of solar and lunar perturbations. “

\sim

” indicates neglected contributions, considering that the

J_{2}

perturbation induces a negligible eccentricity variation of

10^{-6}

. Terms like

\sigma_{J_{2}}^{s}

are also ignored, and

\sigma_{J_{2}}^{l[l]}=\sigma_{J_{2}}^{l[c]}=0

. Explicit expressions for each component are detailed in Appendix C.
表 2：天琴轨道的解析解

\sigma(t)

的分量，描述了太阳、月亮和地球的

J_{2}

扰动效应。符号 “

-

” 表示由于太阳和月球扰动的共同影响，术语不会单独出现。“

\sim

” 表示被忽略的贡献，考虑到

J_{2}

扰动引起的的偏心率变化可以忽略不计。

10^{-6}

术语 like

\sigma_{J_{2}}^{s}

也被忽略，和

\sigma_{J_{2}}^{l[l]}=\sigma_{J_{2}}^{l[c]}=0

.附录 C 中详细介绍了每个组件的显式表达式。

$\sigma(t)$	$\sigma_{0}+\delta_{\iota}\overline{n}\,t$	$\sigma_{c}(t)$			$\sigma_{l[c]}(t)$		$\sigma_{l[l]}(t)$		$\sigma_{s}(t)$
$\sigma(t)$	$\sigma_{0}+\delta_{\iota}\overline{n}\,t$	s	m	$J_{2}$	s	m	s	m	s	m
$a(t)$	$a_{0}$	0	0	0	0	0	0	0	$a_{\text{s}}^{s}$	$a_{\text{m}}^{s}$
$i(t)$	$i_{0}$	0	0	$i_{J_{2}}^{c}$	0	$i_{\text{m}}^{l[c]}$	$i_{\text{s}}^{l[l]}$	$i_{\text{m}}^{l[l]}$	$i_{\text{s}}^{s}$	$i_{\text{m}}^{s}$
$\xi(t)$	$\xi_{0}$	$-$	$\xi_{c}$	$\sim$	$-$	$\xi_{l[c]}$	0	$\xi_{\text{m}}^{l[l]}$	$\xi_{\text{s}}^{s}$	$\xi_{\text{m}}^{s}$
$\eta(t)$	$\eta_{0}$	$-$	$\eta_{c}$	$\sim$	$-$	$\eta_{l[c]}$	0	$\eta_{\text{m}}^{l[l]}$	$\eta_{\text{s}}^{s}$	$\eta_{\text{m}}^{s}$
$\Omega(t)$	$\Omega_{0}$	$\Omega_{\text{s}}^{c}$	$\Omega_{\text{m}}^{c}$	$\Omega_{J_{2}}^{c}$	0	$\Omega_{\text{m}}^{l[c]}$	$\Omega_{\text{s}}^{l[l]}$	$\Omega_{\text{m}}^{l[l]}$	$\Omega_{\text{s}}^{s}$	$\Omega_{\text{m}}^{s}$
$\lambda(t)$	$\lambda_{0}+\overline{n}\,t$	$\lambda_{\text{s}}^{c}$	$\lambda_{\text{m}}^{c}$	$\lambda_{J_{2}}^{c}$	0	$\lambda_{\text{m}}^{l[c]}$	$\lambda_{\text{s}}^{l[l]}$	$\lambda_{\text{m}}^{l[l]}$	$\lambda_{\text{s}}^{s}$	$\lambda_{\text{m}}^{s}$

To enhance the precision of the analytical solution, we employ a perturbation method known as the mean element method [34, 35], which uses the mean orbital elements $\overline{\sigma}(t)$ as a reference solution (defined in Eq. (122)), rather than the Keplerian orbit $\sigma^{(0)}(t)$ . Additionally, for simplicity, only terms up to the first order of eccentricity in the solution are considered. For more derivation details, one can refer to Appendix B. The components of the perturbation solution are outlined in Table 2, with the average value $\overline{n}=\sqrt{\frac{\mu}{\overline{a}^{3}}}$ replacing $n$ in $\sigma^{(0)}(t)$ . Detailed expressions for the terms in Table 2 can be found in Appendix C. Moreover, the effectiveness of the analytical solution $\sigma(t)$ is evaluated through a comparison with high-fidelity numerical orbit simulations (see Appendix B.3). For the TianQin orbit, the 5-year average deviation in position is approximately 87 km.
为了提高解析解的精度，我们采用了一种称为均元法 [34， 35] 的微扰方法，该方法使用均值轨道元 $\overline{\sigma}(t)$ 作为参考解（在方程（122）中定义），而不是开普勒轨道 $\sigma^{(0)}(t)$ 。此外，为简单起见，仅考虑解中偏心率的一阶项。有关更多推导的详细信息，可以参考附录 B。表 2 概述了扰动解的组成部分，平均值 $\overline{n}=\sqrt{\frac{\mu}{\overline{a}^{3}}}$ 替换为 $n$ $\sigma^{(0)}(t)$ 。表 2 中术语的详细表达式可在附录 C 中找到。此外，通过与高保真数值轨道模拟的比较来评估解析解 $\sigma(t)$ 的有效性（见附录 B.3）。对于天琴轨道，5年的平均位置偏差约为87公里。

Table 2 illustrates the effects of gravitational perturbations on the orbital elements of the TianQin satellite. Variation in $a$ are solely induced by short-period perturbation. However, the other five orbital elements are also influenced by both secular and long-period perturbations, particularly the secular one, leading to cumulative change. In the case of $i$ , its secular variation are not due to lunisolar perturbations (which would occur when considering second-order eccentricity [47]), but instead result from $J_{2}$ perturbation, tied to a coordinate transformation involving the obliquity $\epsilon$ (see Eq. (182)). As for $\Omega$ , $\lambda$ , $\xi$ , and $\eta$ , their secular variations are predominantly driven by lunar and solar perturbations.
表 2 说明了引力扰动对天琴卫星轨道元件的影响。的变化 $a$ 完全是由短周期扰动引起的。然而，其他五种轨道元素也受到长期和长期扰动的影响，尤其是长期扰动，导致累积变化。在 $i$ 的情况下，它的长期变化不是由于阴阳扰动（在考虑二阶偏心率时会发生 [47]），而是由 $J_{2}$ 扰动引起的，与涉及倾斜度的坐标变换有关 $\epsilon$ （见方程（182））。至于 $\Omega$ 、 $\lambda$ 、 $\xi$ 和 $\eta$ ，它们的长期变化主要是由月球和太阳的扰动驱动的。

The two elements, $\Omega$ and $i$ , determine the orientation of the orbital plane. As indicated in Eqs. (137) and (156), the secular variation of $\Omega$ is dependent on the satellite’s mean semimajor axis $\overline{a}$ and mean inclination $\overline{i}$ , which implies that $\Omega$ experiences negligible precession when $i\sim 90^{\circ}$ . Similarly, the secular variation of $i$ for the TianQin satellite is also minimal ( $<0.1^{\circ}$ in five years). As a result, the orientation of the TianQin detector plane remains nearly constant, changing by less than $2.6^{\circ}$ over five years. This is in stark contrast to LISA orbits [4], where the plane undergoes a full 360-degree rotation annually.
这两个单元 $\Omega$ 和 $i$ 确定轨道平面的方向。如 Eqs 所示。（137）和（156）中，的长期 $\Omega$ 变化取决于卫星的平均长半轴 $\overline{a}$ 和平均倾角 $\overline{i}$ ，这意味着 $\Omega$ 当时经历的岁差可以忽略不计。 $i\sim 90^{\circ}$ 同样，天琴卫星的长期变化 $i$ 也很小（ $<0.1^{\circ}$ 五年内）。因此，天琴探测器平面的方向几乎保持不变，变化时间不到 $2.6^{\circ}$ 五年。这与 LISA 轨道 [4] 形成鲜明对比，LISA 轨道上的飞机每年都会进行完整的 360 度旋转。

For the periodic variations, their periods are linked to the motions of the satellite, the Moon, and the Sun. Especially, for the short-period variation $\sigma_{s}(t)$ , the arguments of trigonometric terms take the form of $\kappa\,\lambda+p\,u_{3}+q\,\theta_{3}$ or $\kappa\,\lambda+p\,U_{2}$ (see, e.g., Eqs. (174) and (146)), which indicates that orbital variations occur at multiples of the satellite’s orbital frequency and are modulated by the motions of the Moon and the Sun. This insight aids in understanding the perturbing effects of the Moon’s and the Sun’s gravitational fields on the TianQin inter-satellite range acceleration noise (cf. Fig. 3 in Ref. [38]).
对于周期性变化，它们的周期与卫星、月亮和太阳的运动有关。特别是，对于短周期变化 $\sigma_{s}(t)$ ，三角项的参数采用 $\kappa\,\lambda+p\,u_{3}+q\,\theta_{3}$ or $\kappa\,\lambda+p\,U_{2}$ 的形式（参见，例如，方程。（174）和（146）），这表明轨道变化发生在卫星轨道频率的倍数处，并受到月球和太阳运动的调制。这种洞察力有助于理解月球和太阳引力场对天琴星间距离加速度噪声的扰动影响（参见参考文献 [38] 中的图 3）。

Moreover, $\sigma_{s}(t)$ is correlated with the orbit phase $\lambda_{k}(t)$ of SC $k$ , with a 120-degree phase difference among the three satellites, indicating that short-period perturbation will influence the relative motion between satellites. However, it can be demonstrated that the other three components of Eq. (75) may have little impact on the relative motion. Furthermore, an ideal equilateral-triangle constellation requires zero eccentricity, which is unlikely to hold, as shown by the perturbation solutions $\xi_{k}^{(1)}(t)$ and $\eta_{k}^{(1)}(t)$ . Consequently, the presence of these components underscores the potential to significantly disturb the nominal TianQin triangle constellation.
此外， $\sigma_{s}(t)$ 与 SC $k$ 的轨道相位 $\lambda_{k}(t)$ 相关，三颗卫星之间的相位差为 120 度，表明短周期扰动会影响卫星之间的相对运动。然而，可以证明方程（75）的其他三个分量可能对相对运动影响很小。此外，理想的等边三角形星座需要零偏心率，这不太可能成立，如扰动解 $\xi_{k}^{(1)}(t)$ 和 $\eta_{k}^{(1)}(t)$ 所示。因此，这些组件的存在强调了显著扰乱名义上的天琴三角形星座的可能性。

III.4 Perturbed motion of the TianQin constellation
III.4 天琴星座的扰动运动

Equations (73)-(75) describe the variations in orbital elements for SC $k$ ( $k$ = 1, 2, 3) under the influence of lunisolar perturbations and Earth’s $J_{2}$ perturbation in the geocentric ecliptic coordinate system,
方程（73）-（75）描述了在月球太阳扰动和地球 $J_{2}$ 扰动在地心黄道坐标系的影响下SC（ $k$ = 1， 2， 3）的轨道元素 $k$ 的变化，

	$\displaystyle\sigma_{k}(t)=$	$\displaystyle~{}\sigma_{0k}+\delta_{\iota}\,\overline{n}_{k}(t-t_{0})$
		$\displaystyle+\bigl{(}\sigma_{k}^{c}(t)+\sigma_{k}^{l[c]}(t)+\sigma_{k}^{l[l]}% (t)+\sigma_{k}^{s}(t)\bigr{)}\!\large{\|}_{t_{0}}^{t},$		(76)

with their explicit expressions detailed in Table 2 and Appendix C. Substituting Eq. (76) into Eq. (24) and using Eq. (51), the position vector $\mathbf{r}_{k}$ and velocity vector $\mathbf{\dot{r}}_{k}$ can be obtained. Then, employing the definitions in Eqs. (28)-(30), one can derive analytical expressions for the constellation’s three kinematic indicators, $L_{ij}(t)$ , $v_{ij}(t)$ , and $\alpha_{k}(t)$ . The time evolution of these three quantities is illustrated in Fig. 3 for a set of simulated TianQin orbits, comparing both analytical and numerical models (see Appendix B.4 for more details).
表2和附录C中详细说明了其明确的表达方式。将式（76）代入式（24）中，使用式（51）， $\mathbf{\dot{r}}_{k}$ 可以得到位置矢量 $\mathbf{r}_{k}$ 和速度矢量。然后，采用方程中的定义。（28）-（30），可以推导出星座的三个运动学指标 $L_{ij}(t)$ 、 $v_{ij}(t)$ 和 $\alpha_{k}(t)$ 的解析表达式。这三个量的时间演化如图所示。图3为一组模拟的天琴轨道，比较了分析和数值模型（详见附录B.4）。

Variations in the triangular constellation can be decomposed into two parts,
三角形星座的变化可以分解为两部分，

$\displaystyle L_{ij}(t)=$	$\displaystyle~{}L_{\text{o}}+\delta L_{ij}(t),$	(77)
$\displaystyle v_{ij}(t)=$	$\displaystyle~{}v_{\text{o}}+\delta v_{ij}(t),$	(78)
$\displaystyle\alpha_{k}(t)=$	$\displaystyle~{}\alpha_{\text{o}}+\delta\alpha_{k}(t),$	(79)

where $L_{\text{o}}=\sqrt{3}\times 10^{5}$ km, $v_{\text{o}}=0$ m/s, and $\alpha_{\text{o}}=\frac{\pi}{3}$ represent the desired equilateral-triangle configuration. Conversely, $\delta L_{ij}(t)$ , $\delta v_{ij}(t)$ , and $\delta\alpha_{k}(t)$ signify distortions from the ideal configuration. As previously mentioned, the magnitude of these distortions significantly impacts GW detection missions, including TDI data processing and the design of instruments such as phase meters and telescopes. It is crucial to minimize these distortions as much as possible.
其中 $L_{\text{o}}=\sqrt{3}\times 10^{5}$ km、 $v_{\text{o}}=0$ m/s 和 $\alpha_{\text{o}}=\frac{\pi}{3}$ 表示所需的等边三角形配置。反之， $\delta L_{ij}(t)$ 、、 $\delta v_{ij}(t)$ 和 $\delta\alpha_{k}(t)$ 表示与理想配置的扭曲。如前所述，这些扭曲的程度极大地影响了GW探测任务，包括TDI数据处理和相位计和望远镜等仪器的设计。尽可能减少这些失真至关重要。

These distortions are all zero when the three satellites are solely influenced by the Earth’s point mass and satisfy the conditions (31) and (36). However, these conditions no longer hold when accounting for gravitational perturbations, which result in variations in $e_{k}(t)$ and the inclusion of non-synchronous short-period terms $\sigma_{k}^{s}(t)$ in $\sigma_{k}(t)$ . To achieve an equilateral triangle, $\sigma_{k}(t)$ in Eqs. (31) and (36), instead, can be substituted with $\sigma_{k}^{\text{ref}}(t)$ :
当三颗卫星仅受地球点质量的影响并满足条件（31）和（36）时，这些失真都为零。然而，在考虑引力扰动时，这些条件不再成立，引力扰动导致中的变化 $e_{k}(t)$ 和不同步的短周期项 $\sigma_{k}^{s}(t)$ 的包含 $\sigma_{k}(t)$ 。为了实现等边三角形， $\sigma_{k}(t)$ 在 Eqs.（31）和（36）可以替换为 $\sigma_{k}^{\text{ref}}(t)$ ：

\displaystyle\left\{\begin{array}[]{l}e_{1}^{\text{ref}}(t)=e_{2}^{\text{ref}}% (t)=e_{3}^{\text{ref}}(t)\equiv e_{\text{o}}=0,\vskip 3.0pt plus 1.0pt minus 1% .0pt\\ a_{1}^{\text{ref}}(t)=a_{2}^{\text{ref}}(t)=a_{3}^{\text{ref}}(t)\equiv a_{% \text{o}}=\frac{L_{\text{o}}}{\sqrt{3}},\vskip 3.0pt plus 1.0pt minus 1.0pt\\ i_{1}^{\text{ref}}(t)=i_{2}^{\text{ref}}(t)=i_{3}^{\text{ref}}(t)\equiv i_{% \text{o}}(t),\vskip 3.0pt plus 1.0pt minus 1.0pt\\ \Omega_{1}^{\text{ref}}(t)=\Omega_{2}^{\text{ref}}(t)=\Omega_{3}^{\text{ref}}(% t)\equiv\Omega_{\text{o}}(t),\vskip 3.0pt plus 1.0pt minus 1.0pt\\ \lambda_{k}^{\text{ref}}(t)\equiv\frac{2\pi}{3}(k-1)+\lambda_{\text{o}}(t),% \end{array}\right.

(85)

where 哪里

\sigma_{k}^{\text{ref}}(t):=\left\{\begin{array}[]{ll}0,&\sigma=e,\vskip 3.0pt% plus 1.0pt minus 1.0pt\\ \sigma_{k}(t)-\sigma_{k}^{s}(t),&\sigma\in\{a,i,\Omega,\lambda\},\end{array}\right.

(86)

and $\sigma_{\text{o}}(t)$ serves as the reference orbit for the synchronous motion of the three satellites. Utilizing Eqs. (76), (85), and (86), the form of $\sigma_{\text{o}}(t)$ for $\sigma\in\{a,i,\Omega,\lambda\}$ is
并 $\sigma_{\text{o}}(t)$ 作为三颗卫星同步运动的参考轨道。利用方程。（76）、（85）和（86），for $\sigma\in\{a,i,\Omega,\lambda\}$ 的形式 $\sigma_{\text{o}}(t)$ 是

\displaystyle\sigma_{\text{o}}(t)=\overline{\sigma}_{0\text{o}}+\delta_{\iota}% \overline{n}_{\text{o}}(t-t_{0})+\Delta\sigma_{\text{o}}^{c}(t)+\Delta\sigma_{% \text{o}}^{l[c]}(t)+\sigma_{\text{o}}^{l[l]}(t),

(87)

with $\overline{\sigma}_{0\text{o}}:=\sigma_{0k}-[\sigma_{\text{o}}^{l[l]}(t_{0})+% \sigma_{k}^{s}(t_{0})]-\delta_{\iota}\frac{2\pi}{3}(k-1)=\text{const}.$ , $\Delta\sigma_{\text{o}}^{c}(t):=\sigma_{\text{o}}^{c}(t)-\sigma_{\text{o}}^{c}% (t_{0})$ , and $\Delta\sigma_{\text{o}}^{l[c]}(t):=\sigma_{\text{o}}^{l[c]}(t)-\sigma_{\text{o% }}^{l[c]}(t_{0})$ . In other words, the three satellites move along the same virtual circular orbit $\sigma_{\text{o}}(t)$ with secular and long-period variations, while maintaining a 120-degree phase difference, forming a rotating, precessing equilateral-triangle constellation.
其中 $\overline{\sigma}_{0\text{o}}:=\sigma_{0k}-[\sigma_{\text{o}}^{l[l]}(t_{0})+% \sigma_{k}^{s}(t_{0})]-\delta_{\iota}\frac{2\pi}{3}(k-1)=\text{const}.$ 、 $\Delta\sigma_{\text{o}}^{c}(t):=\sigma_{\text{o}}^{c}(t)-\sigma_{\text{o}}^{c}% (t_{0})$ 和 $\Delta\sigma_{\text{o}}^{l[c]}(t):=\sigma_{\text{o}}^{l[c]}(t)-\sigma_{\text{o% }}^{l[c]}(t_{0})$ 。换句话说，这三颗卫星沿着同一个虚拟圆形轨道移动 $\sigma_{\text{o}}(t)$ ，具有长期和长期变化，同时保持 120 度的相位差，形成一个旋转的进动等边三角形星座。

Correspondingly, $\delta L_{ij}(t)$ , $\delta v_{ij}(t)$ , and $\delta\alpha_{k}(t)$ result from discrepancies $\delta\sigma_{k}(t)$ :
相应地， $\delta L_{ij}(t)$ ， $\delta v_{ij}(t)$ ，和 $\delta\alpha_{k}(t)$ 则由差异产生 $\delta\sigma_{k}(t)$ ：

\displaystyle\delta\sigma_{k}(t):=\sigma_{k}(t)-\sigma_{\text{o}k}(t),

(88)

between the real orbits $\sigma_{k}(t)$ and the reference orbits $\sigma_{\text{o}k}(t)=\sigma_{\text{o}}(t)+\delta_{\iota}\frac{2\pi}{3}(k-1)$ . By expanding $L_{ij}(t)$ , $v_{ij}(t)$ , and $\alpha_{k}(t)$ into a Taylor series along $\sigma_{\text{o}k}(t)$ , the triangle distortions caused by $\delta\sigma_{k}(t)$ can be obtained. For the arm-length distortion, we have
在实轨道 $\sigma_{k}(t)$ 和参考轨道之间 $\sigma_{\text{o}k}(t)=\sigma_{\text{o}}(t)+\delta_{\iota}\frac{2\pi}{3}(k-1)$ 。通过将、 $v_{ij}(t)$ 和 $\alpha_{k}(t)$ 沿 $\sigma_{\text{o}k}(t)$ 展开 $L_{ij}(t)$ 为泰勒级数，可以得到引起的 $\delta\sigma_{k}(t)$ 三角形畸变。对于臂长畸变，我们有

$\displaystyle\delta L_{ij}(t)=$	$\displaystyle~{}L_{ij}(t)-\sqrt{3}a_{\text{o}}$
$\displaystyle=$	$\displaystyle~{}\frac{\sqrt{3}}{2}[{\delta a}_{i}(t)+{\delta a}_{j}(t)]+\frac{% 1}{2}a_{\text{o}}\delta\lambda_{ji}(t)+\frac{1}{2}a_{\text{o}}\cos i_{\text{o}% }(t)\,\delta\Omega_{{ji}}(t)$
	$\displaystyle+\frac{\sqrt{7}}{2}a_{\text{o}}[\sin M_{\text{o}j}^{-}(t)\,{% \delta e}_{j}(t)-\sin M_{\text{o}i}^{+}(t)\,{\delta e}_{i}(t)]+\mathcal{O}(% \delta\sigma(t)^{2}),$	(89)

where 哪里

\displaystyle\delta\sigma_{ji}(t):=\delta\sigma_{j}(t)-\delta\sigma_{i}(t),

(90)

\displaystyle M_{\text{o}k}^{\pm}(t):=M_{\text{o}k}(t)\pm\beta,\qquad M_{\text% {o}k}(t)=\lambda_{\text{o}k}(t)-\omega_{k}(t),

(91)

and the indices $i$ , $j$ , and $k$ follow a cyclic permutation ( $i,\,j,\,k=1\to 2\to 3\to 1$ ). It can be seen from Eq. (89) that, up to $(\delta\sigma)^{1}$ order, $\delta L_{ij}(t)$ is unaffected by inclination deviations $\delta i_{i}$ and $\delta i_{j}$ . Additionally, deviations in $\Omega$ have minimal influence on $\delta L_{ij}(t)$ due to the approximately $90^{\circ}$ inclinations of TianQin orbits, which render the constellation stability insensitive to orbital plane deviations. Combining Eqs. (90), (88), (76), and (87), and defining $i_{\text{o}}(t)=\overline{i}_{\text{o}}+i^{\epsilon}_{\text{o}}(t)$ , where $\overline{i}_{\text{o}}$ is the average, the right-hand side of Eq. (89) can be categorized into distinct types:
和索引 $i$ ， $j$ 并 $k$ 遵循循环排列（ $i,\,j,\,k=1\to 2\to 3\to 1$ ）。从方程（89）可以看出，在阶 $(\delta\sigma)^{1}$ 次下， $\delta L_{ij}(t)$ 不受倾角偏差和 $\delta i_{j}$ 的影响 $\delta i_{i}$ 。此外，由于天琴轨道的近似 $90^{\circ}$ 倾角，星 $\Omega$ 座的偏差对 $\delta L_{ij}(t)$ 轨道平面偏差的影响很小，这使得星座稳定性对轨道平面偏差不敏感。结合方程。（90）、（88）、（76）和（87），定义 $i_{\text{o}}(t)=\overline{i}_{\text{o}}+i^{\epsilon}_{\text{o}}(t)$ ，其中 $\overline{i}_{\text{o}}$ 是平均值，方程（89）的右侧可以分为不同的类型：

\displaystyle\delta L_{ij}(t)=

\displaystyle~{}\delta L_{ij}^{\text{drift}}(t)+\delta L_{ij}^{\text{bias}}+% \delta L_{ij}^{\text{fluc}}(t)+\mathcal{O}(i^{\epsilon}_{\text{o}}(t)^{1})% \delta\Omega_{{ji}}(t)+\mathcal{O}(\delta\sigma(t)^{2}),

(92)

with 跟

$\displaystyle\delta L_{ij}^{\text{drift}}(t)=$	$\displaystyle~{}\frac{1}{2}a_{\text{o}}\delta\overline{n}_{ji}(t-t_{0})+\frac{% 1}{2}a_{\text{o}}[\delta_{\!\Delta}\lambda_{ji}^{c}(t)+\delta_{\!\Delta}% \lambda_{ji}^{l[c]}(t)]+\frac{1}{2}a_{\text{o}}\cos\overline{i}_{\text{o}}\,[% \delta_{\!\Delta}\Omega_{ji}^{c}(t)+\delta_{\!\Delta}\Omega_{ji}^{l[c]}(t)],$	(93)
$\displaystyle\delta L_{ij}^{\text{bias}}=$	$\displaystyle~{}\frac{\sqrt{3}}{2}(\delta\overline{a}_{0i}+\delta\overline{a}_% {0j})+\frac{1}{2}a_{\text{o}}\delta\overline{\lambda}_{0{ji}}+\frac{1}{2}a_{% \text{o}}\cos\overline{i}_{\text{o}}\,\delta\overline{\Omega}_{0{ji}},$	(94)
$\displaystyle\delta L_{ij}^{\text{fluc}}(t)=$	$\displaystyle~{}\frac{\sqrt{7}}{2}a_{\text{o}}[\sin M_{\text{o}j}^{-}(t)\,e_{j% }(t)-\sin M_{\text{o}i}^{+}(t)\,e_{i}(t)]+\frac{\sqrt{3}}{2}[a^{s}_{i}(t)+a^{s% }_{j}(t)]$
	$\displaystyle+\frac{1}{2}a_{\text{o}}[\delta\lambda_{ji}^{l[l]}(t)+\delta% \lambda^{s}_{ji}(t)]+\frac{1}{2}a_{\text{o}}\cos\overline{i}_{\text{o}}\,[% \delta\Omega_{ji}^{l[l]}(t)+\delta\Omega^{s}_{ji}(t)],$	(95)

where $\delta\overline{n}_{ji}=\overline{n}_{j}-\overline{n}_{i}$ , 其中 $\delta\overline{n}_{ji}=\overline{n}_{j}-\overline{n}_{i}$ ，

\displaystyle\delta_{\!\Delta}\sigma_{ij}(t):=\delta\sigma_{ij}(t)-\delta% \sigma_{ij}(t_{0}),

(96)

$\delta\overline{a}_{0{i}}=\overline{a}_{0{i}}-a_{\text{o}}$ , $\delta\overline{\lambda}_{0{ji}}=\overline{\lambda}_{0{j}}-\overline{\lambda}_% {0{i}}-\frac{2\pi}{3}(j-i)$ , and so on. Equation (92) illustrates that arm-length distortion manifests in three possible types: linear drift $\delta L_{ij}^{\text{drift}}(t)$ , constant bias $\delta L_{ij}^{\text{bias}}$ , and periodic fluctuation $\delta L_{ij}^{\text{fluc}}(t)$ . Specifically, (1) $\delta L_{ij}^{\text{drift}}(t)$ consists of five components of inter-satellite deviations, including $\delta\overline{n}_{ji}$ , $\delta_{\!\Delta}\lambda_{ji}^{c}$ , $\delta_{\!\Delta}\Omega_{ji}^{c}$ , $\delta_{\!\Delta}\lambda_{ji}^{l[c]}$ , and $\delta_{\!\Delta}\Omega_{ji}^{l[c]}$ ; (2) $\delta L_{ij}^{\text{bias}}$ comprises initial mean deviations, $\delta\overline{a}_{0i}$ , $\delta\overline{a}_{0j}$ , $\delta\overline{\lambda}_{0{ji}}$ , and $\delta\overline{\Omega}_{0ji}$ ; and (3) $\delta L_{ij}^{\text{fluc}}(t)$ is linked to eccentricity variations $e_{i}(t)$ and $e_{j}(t)$ and short-period variations in semi-major axis $a^{s}_{i}(t)$ and $a^{s}_{j}(t)$ , along with inter-satellite periodic deviations in $\lambda$ and $\Omega$ . Among these types, the drift, which progressively increases over time, emerges as the predominant factor affecting the stability of the constellation.
$\delta\overline{a}_{0{i}}=\overline{a}_{0{i}}-a_{\text{o}}$ 、 $\delta\overline{\lambda}_{0{ji}}=\overline{\lambda}_{0{j}}-\overline{\lambda}_% {0{i}}-\frac{2\pi}{3}(j-i)$ 、等。方程（92）说明了臂长失真表现为三种可能的类型：线性漂移 $\delta L_{ij}^{\text{drift}}(t)$ 、恒定偏置 $\delta L_{ij}^{\text{bias}}$ 和周期性波动 $\delta L_{ij}^{\text{fluc}}(t)$ 。具体来说，（1） $\delta L_{ij}^{\text{drift}}(t)$ 由卫星间偏差的五个分量组成，包括 $\delta\overline{n}_{ji}$ 、、 $\delta_{\!\Delta}\lambda_{ji}^{c}$ $\delta_{\!\Delta}\Omega_{ji}^{c}$ 、 $\delta_{\!\Delta}\lambda_{ji}^{l[c]}$ 和 $\delta_{\!\Delta}\Omega_{ji}^{l[c]}$ ;（2） $\delta L_{ij}^{\text{bias}}$ 包括初始平均偏差 $\delta\overline{a}_{0i}$ 、、 $\delta\overline{a}_{0j}$ $\delta\overline{\lambda}_{0{ji}}$ 和 $\delta\overline{\Omega}_{0ji}$ ;（3） $\delta L_{ij}^{\text{fluc}}(t)$ 与半长轴 $a^{s}_{i}(t)$ 和 $a^{s}_{j}(t)$ 的偏心率变化 $e_{i}(t)$ 和 $e_{j}(t)$ 短周期变化以及和 $\Omega$ 的 $\lambda$ 卫星间周期性偏差有关。在这些类型中，随着时间的推移逐渐增加的漂移成为影响星座稳定性的主要因素。

Regarding the relative velocity, one has
关于相对速度，有

$\displaystyle\delta v_{ij}(t)=$	$\displaystyle~{}v_{ij}(t)-0$
$\displaystyle=$	$\displaystyle-\frac{3}{4}\frac{v_{\text{o}}}{a_{\text{o}}}[{\delta a}_{j}(t)-{% \delta a}_{i}(t)]+\frac{\sqrt{7}}{2}v_{\text{o}}[\cos M_{\text{o}j}^{-}(t)\,{% \delta e}_{j}(t)-\cos M_{\text{o}i}^{+}(t)\,{\delta e}_{i}(t)]+\mathcal{O}(% \delta\sigma(t)^{2}),$	(97)
$\displaystyle=$	$\displaystyle~{}\delta v_{ij}^{\text{bias}}+\delta v_{ij}^{\text{fluc}}(t)+% \mathcal{O}(\delta\sigma(t)^{2}),$	(98)

where $v_{\text{o}}:=\sqrt{\frac{\mu}{a_{\text{o}}}}$ , and 其中 $v_{\text{o}}:=\sqrt{\frac{\mu}{a_{\text{o}}}}$ ，和

	$\displaystyle\delta v_{ij}^{\text{bias}}=$	$\displaystyle-\frac{3}{4}\frac{v_{\text{o}}}{a_{\text{o}}}(\delta\overline{a}_% {0j}-\delta\overline{a}_{0i}),$		(99)
	$\displaystyle\delta v_{ij}^{\text{fluc}}(t)=$	$\displaystyle-\frac{3}{4}\frac{v_{\text{o}}}{a_{\text{o}}}[a^{s}_{j}(t)-a^{s}_% {i}(t)]+\frac{\sqrt{7}}{2}v_{\text{o}}[\cos M_{\text{o}j}^{-}(t)\,{e}_{j}(t)-% \cos M_{\text{o}i}^{+}(t)\,{e}_{i}(t)].$		(100)

Equation (98) illustrates that there is little long-term variation in relative velocity, $\delta v_{ij}^{\text{drift}}(t)\simeq 0$ , consistent with numerical simulation results (cf. Fig. 10 in Ref. [41]). Additionally, the breathing angle within the TianQin triangle experiences three types of distortion akin to those observed in arm-length:
方程（98）表明相对速度的长期变化很小， $\delta v_{ij}^{\text{drift}}(t)\simeq 0$ 这与数值模拟结果一致（参见参考文献 [41] 中的图 10）。此外，天琴三角内的呼吸角经历了三种类型的扭曲，类似于在手臂长度中观察到的扭曲：

$\displaystyle\delta\alpha_{k}(t)=$	$\displaystyle~{}\alpha_{k}(t)-\frac{\pi}{3}$
$\displaystyle=$	$\displaystyle~{}\frac{1}{2\sqrt{3}a_{\text{o}}}[{\delta a}_{i}(t)+{\delta a}_{% j}(t)-2{\delta a}_{k}(t)]+\frac{1}{2}\delta\lambda_{ji}(t)+\frac{1}{2}\cos i_{% \text{o}}(t)\,\delta\Omega_{{ji}}(t)$
	$\displaystyle+\frac{\sqrt{7}}{6}[f_{e}^{k}(t){\delta e}_{k}(t)+f_{e}^{i}(t){% \delta e}_{i}(t)+f_{e}^{j}(t){\delta e}_{j}(t)]+\mathcal{O}(\delta\sigma(t)^{2})$	(101)
$\displaystyle=$	$\displaystyle~{}\delta\alpha_{k}^{\text{drift}}(t)+\delta\alpha_{k}^{\text{% bias}}+\delta\alpha_{k}^{\text{fluc}}(t)+\mathcal{O}(i^{\epsilon}_{\text{o}}(t% )^{1})\delta\Omega_{{ji}}(t)+\mathcal{O}(\delta\sigma(t)^{2}),$	(102)

with 跟

$\displaystyle\delta\alpha_{k}^{\text{drift}}(t)=$	$\displaystyle~{}\frac{1}{2}\delta\overline{n}_{ji}(t-t_{0})+\frac{1}{2}[\delta% _{\!\Delta}\lambda_{ji}^{c}(t)+\delta_{\!\Delta}\lambda_{ji}^{l[c]}(t)]+\frac{% 1}{2}\cos\overline{i}_{\text{o}}\,[\delta_{\!\Delta}\Omega_{ji}^{c}(t)+\delta_% {\!\Delta}\Omega_{ji}^{l[c]}(t)],$	(103)
$\displaystyle\delta\alpha_{k}^{\text{bias}}=$	$\displaystyle~{}\frac{1}{2\sqrt{3}a_{\text{o}}}(\delta\overline{a}_{0ik}+% \delta\overline{a}_{0jk})+\frac{1}{2}\delta\overline{\lambda}_{0{ji}}+\frac{1}% {2}\cos\overline{i}_{\text{o}}\,\delta\overline{\Omega}_{0{ji}},$	(104)
$\displaystyle\delta\alpha_{k}^{\text{fluc}}(t)=$	$\displaystyle~{}\frac{\sqrt{7}}{6}[f_{e}^{k}(t)e_{k}(t)+f_{e}^{i}(t)e_{i}(t)+f% _{e}^{j}(t)e_{j}(t)]+\frac{1}{2\sqrt{3}a_{\text{o}}}[\delta a^{s}_{ik}(t)+% \delta a^{s}_{jk}(t)]$
	$\displaystyle+\frac{1}{2}[\delta\lambda_{ji}^{l[l]}(t)+\delta\lambda^{s}_{ji}(% t)]+\frac{1}{2}\cos\overline{i}_{\text{o}}\,[\delta\Omega_{ji}^{l[l]}(t)+% \delta\Omega^{s}_{ji}(t)],$	(105)

where 哪里

\displaystyle f_{e}^{k}(t)=2\cos M_{\text{o}k}(t)\sin\beta,\qquad f_{e}^{i}(t)% =-\sin M_{\text{o}i}^{-}(t)-2\sin M_{\text{o}i}^{+}(t),\qquad f_{e}^{j}(t)=2% \sin M_{\text{o}j}^{-}(t)+\sin M_{\text{o}j}^{+}(t).

(106)

Equations (104) and (105) show that bias and fluctuation in the breathing angle, as observed from SC $k$ , are associated with deviations in $a$ and $e$ of all three satellites. However, concerning $\lambda$ and $\Omega$ , they are exclusively linked to the relative differences between the other two satellites.
方程（104）和（105）表明，从 SC $k$ 观察到的呼吸角度的偏差和波动与所有三个卫星的 $a$ $e$ 偏差有关。然而，关于 $\lambda$ 和 $\Omega$ ，它们完全与其他两颗卫星之间的相对差异有关。

The variations in the three types, drift, bias, and fluctuation, all impact the constellation’s stability, necessitating optimization. The drift, associated with $\delta\overline{n}_{ji}=-\frac{3n_{\text{o}}}{2a_{\text{o}}}{\delta\bar{a}}_{{% ji}}$ , can be significantly mitigated by aligning the mean semi-major axes $\overline{a}_{j}^{\,\text{optim}}=\overline{a}_{i}^{\,\text{optim}}$ . More generally, one can see from, e.g., Eqs. (103)-(105), that terms within the drift, bias, and long-period fluctuation, are contingent on the mean or initial mean values of parameters $a$ , $i$ , $\Omega$ , and $\lambda$ . Additionally, the fluctuation is also correlated with eccentricities, which exhibit secular variations and serve as the primary factor influencing the amplitude of the fluctuation. Thus, the optimization of constellation variations can be achieved by imposing the following conditions:
漂移、偏置和波动这三种类型的变化都会影响星座的稳定性，因此需要进行优化。与 $\delta\overline{n}_{ji}=-\frac{3n_{\text{o}}}{2a_{\text{o}}}{\delta\bar{a}}_{{% ji}}$ 关联的漂移可以通过对齐平均半长轴 $\overline{a}_{j}^{\,\text{optim}}=\overline{a}_{i}^{\,\text{optim}}$ 来显著缓解。更一般地说，人们可以从 Eqs 等中看到。（103）-（105），即漂移、偏差和长期波动内的项取决于参数 $a$ 、 $i$ 、 $\Omega$ 和 $\lambda$ 的平均值或初始平均值。此外，波动还与偏心率相关，偏心率表现出长期变化，是影响波动幅度的主要因素。因此，可以通过施加以下条件来实现星座变化的优化：

\displaystyle\left\{\begin{array}[]{l}\overline{a}_{1}^{\,\text{optim}}=% \overline{a}_{2}^{\,\text{optim}}=\overline{a}_{3}^{\,\text{optim}},\\ \overline{i}_{1}^{\,\text{optim}}=\overline{i}_{2}^{\,\text{optim}}=\overline{% i}_{3}^{\,\text{optim}},\\ \overline{\Omega}_{1}^{\,\text{optim}}=\overline{\Omega}_{2}^{\,\text{optim}}=% \overline{\Omega}_{3}^{\,\text{optim}},\\ \overline{\lambda}_{k}^{\,\text{optim}}=\frac{2\pi}{3}(k-1)+\overline{\lambda}% _{1}^{\,\text{optim}},\\ e_{k}^{\text{optim}}(t)\simeq 0.\end{array}\right.

(112)

Further, the optimized indicators, up to the leading order, are
此外，优化后的指标，直到前导顺序，都是

\displaystyle\delta L_{ij}^{\text{optim}}(t)=

\displaystyle~{}\frac{\sqrt{7}}{2}a_{\text{o}}[\sin M_{\text{o}j}^{-}(t)\,e_{j% }(t)-\sin M_{\text{o}i}^{+}(t)\,e_{i}(t)]+\frac{\sqrt{3}}{2}[a^{s}_{i}(t)+a^{s% }_{j}(t)]+\frac{1}{2}a_{\text{o}}\delta\lambda^{s}_{ji}(t)+\frac{1}{2}a_{\text% {o}}\cos\overline{i}_{\text{o}}\,\delta\Omega^{s}_{ji}(t),

(113)

\displaystyle\delta v_{ij}^{\text{optim}}(t)=

\displaystyle~{}\frac{\sqrt{7}}{2}v_{\text{o}}[\cos M_{\text{o}j}^{-}(t)\,{e}_% {j}(t)-\cos M_{\text{o}i}^{+}(t)\,{e}_{i}(t)]-\frac{3}{4}\frac{v_{\text{o}}}{a% _{\text{o}}}[a^{s}_{j}(t)-a^{s}_{i}(t)],

(114)

\displaystyle\delta\alpha_{k}^{\text{optim}}(t)=

\displaystyle~{}\frac{\sqrt{7}}{6}[f_{e}^{k}(t)e_{k}(t)+f_{e}^{i}(t)e_{i}(t)+f% _{e}^{j}(t)e_{j}(t)]+\frac{1}{2\sqrt{3}a_{\text{o}}}[\delta a^{s}_{ik}(t)+% \delta a^{s}_{jk}(t)]+\frac{1}{2}\delta\lambda^{s}_{ji}(t)+\frac{1}{2}\cos% \overline{i}_{\text{o}}\,\delta\Omega^{s}_{ji}(t),

(115)

with right-hand functions adopting the conditions given by Eq. (112), i.e., $a_{\text{o}}=\overline{a}_{k}^{\,\text{optim}}$ , $\overline{i}_{\text{o}}=\overline{i}_{k}^{\,\text{optim}}$ , $e_{k}(t)=e^{\text{optim}}_{k}(t)$ , etc. Equations (113)-(115) reveal the optimized TianQin triangle as intrinsic fluctuation variations induced by perturbations, with amplitude dependent on eccentricities and short-period variations in other elements. Notably, conditions in Eq. (112) correspond to the optimal stable configuration, and therefore can provide useful guidelines for numerical optimization and orbit control.
右手函数采用式（112）给出的条件，即、 $a_{\text{o}}=\overline{a}_{k}^{\,\text{optim}}$ 、 $\overline{i}_{\text{o}}=\overline{i}_{k}^{\,\text{optim}}$ $e_{k}(t)=e^{\text{optim}}_{k}(t)$ 等。方程（113）-（115）揭示了优化后的天琴三角是由扰动引起的内在涨落变化，其振幅取决于偏心率和其他元素的短周期变化。值得注意的是，方程（112）中的条件对应于最优稳定配置，因此可以为数值优化和轨道控制提供有用的指导。

IV Concluding Remarks
四、结束语

Detecting GWs with TianQin requires a stable three-satellite constellation, configured as closely to an equilateral triangle as possible. In high Earth orbits, gravitational perturbations, especially from lunar and solar influences, can distort this configuration. To quantify the impact, we have developed an analytical model delineating the effects of lunar and solar point masses on the TianQin constellation. This model provides expressions for three kinematic indicators, including arm-lengths, relative velocities, and breathing angles, derived from the first-order perturbation solution for satellite orbital elements.
用天琴探测 GW 需要一个稳定的三星星座，配置尽可能靠近等边三角形。在高地球轨道上，引力扰动，尤其是来自月球和太阳影响的引力扰动，可以扭曲这种配置。为了量化影响，我们开发了一个分析模型，描述了月球和太阳点质量对天琴星座的影响。该模型提供了三个运动学指标的表达式，包括臂长、相对速度和呼吸角，这些指标来自卫星轨道元件的一阶扰动解。

The analysis of these indicators has revealed that gravitational perturbations induce secular, long-period, and short-period variations in satellite orbital elements, leading to relative motion between satellites and distortions in the constellation. These distortions appear as three distinct types, i.e., linear drift, bias, and fluctuation. The drift, progressively increasing over time, is a primary destabilizing factor affecting arm-lengths and breathing angles but has almost no impact on relative velocities. To alleviate design constraints on onboard scientific payloads, these three distortions have been further optimized. It is demonstrated that both drift and bias can be eliminated, and fluctuation amplitude reduced, if the three orbits adhere to the following constraints on average: synchronized orbital periods, aligned orbital planes, equally spaced phases, and minimized eccentricities. The expressions for the optimized indicators are presented, revealing that the optimized TianQin constellation displays only fluctuation with amplitude dependent on eccentricities and short-period variations in other elements.
对这些指标的分析表明，引力扰动会引起卫星轨道元件的长期、长期和短期变化，从而导致卫星之间的相对运动和星座的扭曲。这些失真表现为三种不同的类型，即线性漂移、偏置和波动。漂移随着时间的推移逐渐增加，是影响臂长和呼吸角度的主要不稳定因素，但对相对速度几乎没有影响。为了减轻对机载科学有效载荷的设计限制，这三种失真得到了进一步优化。结果表明，如果三个轨道平均遵守以下约束条件，则可以消除漂移和偏置，并减小波动幅度：同步轨道周期、对齐的轨道平面、等距相位和最小化偏心率。给出了优化指标的表达式，揭示了优化后的天琴星座仅显示振幅取决于偏心率和其他元素的短周期变化的波动。

These results can provide valuable insights and guidelines for enhancing the stability of the GW observatory constellation, such as in numerical optimization and orbit control. For future works, this model will be extended to incorporate the influence of initial orbit errors [22, 48, 49, 50]. The perturbation solution can be applied to explore the dynamics of TianQin satellite eccentricity, closely linked to the constellation stability. Potential applications in celestial mechanics, especially for high-inclination TianQin-like orbits subject to the Kozai–Lidov effect [51, 52], may arise. Further discussions are deferred to future work.
这些结果可以为增强 GW 天文台星座的稳定性提供有价值的见解和指导，例如在数值优化和轨道控制方面。在未来的工作中，该模型将被扩展以纳入初始轨道误差的影响 [22， 48， 49， 50]。该微扰解可用于探究天琴卫星偏心率的动力学特性，与星座稳定性密切相关。在天体力学中的潜在应用，特别是受 Kozai-Lidov 效应影响的高倾角天琴样轨道 [51， 52]，可能会出现。进一步的讨论推迟到未来的工作中。

Acknowledgements. 确认。

The authors thank Jianwei Mei, Yunhe Meng, Defeng Gu, Liang-Cheng Tu, Cheng-Gang Shao, Yan Wang, and Jun Luo for helpful discussions and comments. Special thanks to the anonymous referee for valuable suggestions and comments. X. Z. is supported by the National Key R&D Program of China (Grant Nos. 2020YFC2201202 and 2022YFC2204600), NSFC (Grant No. 12373116), and Fundamental Research Funds for the Central Universities, Sun Yat-sen University (Grant No. 23xkjc001).
作者感谢 Jianwei Mei、Yunhe Meng、Defeng Gu、Liang-Cheng Tu、Cheng-Gang Shao、Yan Wang 和 Jun Luo 的有益讨论和评论。特别感谢匿名审稿人的宝贵建议和评论。X. Z. 得到了国家重点研发计划（批准号 2020YFC2201202 和 2022YFC2204600）、国家自然科学基金（批准号 12373116）和中央高校基本研究基金（批准号 23xkjc001）的支持。

Appendix A Table of Symbols
附录 A 符号表

Table LABEL:table:symbols below lists the main symbols used in the paper and their meanings for quick look-ups.
下面的表格 LABEL：table：symbols 列出了论文中使用的主要符号及其含义，以便快速查找。

Table 3: List of symbols and their meanings.
表 3：符号列表及其含义。

Symbols 符号	Meanings 意义
$t$	Time 时间
$t_{0}$	Reference epoch 参考纪元
UTC	Coordinated Universal Time 协调世界时
$a$	Semimajor axis 半长轴
$e$	Orbital eccentricity 轨道偏心率
$i$	Orbital inclination 轨道倾角
$\Omega$	Longitude of ascending node 升序节点的经度
$\omega$	Argument of perigee 近地点的论点
$\nu$	True anomaly 真正的异常
$M$	Mean anomaly 平均异常
$E$	Eccentric anomaly 偏心异常
$\text{u}\,(=\omega+\nu)$	Argument of latitude 纬度参数
$\lambda$	Defined as $(\omega+M)$ 定义为 $(\omega+M)$
$\xi\,(=e\cos\omega)$	Singularity-free variable for eccentricity 偏心率的无奇点变量
$\eta\,(=e\sin\omega)$	Singularity-free variable for eccentricity 偏心率的无奇点变量
$\sigma$	General representation of orbital elements 轨道元素的一般表示
$\sigma_{0}$	Initial orbital elements 初始轨道元件
Subscript e 下标 e	Refers to the Earth 指地球
Subscript s 下标 s	Refers to the Sun 指太阳
Subscript m 下标 m	Refers to the Moon 指月亮
Subscript $J_{2}$ 下标 $J_{2}$	Refers to the Earth’s oblateness $J_{2}$ 指地球的扁圆度 $J_{2}$
Subscript 0 下标 0	Denoting initial values or zeroth order 表示初始值或零阶
Subscript 1 下标 1	Denoting SC1 or first order 表示 SC1 或一阶
Subscript 2 下标 2	Denoting the Sun, SC2 or second order 表示太阳、SC2 或二阶
Subscript 3 下标 3	Denoting the Moon or SC3 表示月亮或 SC3
Subscript o 下标 o	Denoting nominal values 表示标称值
Notation $c$ 表示法 $c$	Denoting secular variation 表示长期变化
Notation $l$ 表示法 $l$	Denoting long-period variation 表示长周期变化
Notation $l[c]$ 表示法 $l[c]$	Denoting special long-period variation 表示特殊的长周期变化
Notation $l[l]$ 表示法 $l[l]$	Denoting general long-period variation 表示一般的长周期变化
Notation $s$ 表示法 $s$	Denoting short-period variation 表示短期变化
$i$ , $j$ , $k$ $i$ 、 $j$ 、 $k$	Represent satellites and take values of 1, 2, 3 表示卫星并取值 1、2、3
$\kappa$ , $p$ , $q$ $\kappa$ 、 $p$ 、 $q$	Represent components of the perturbation solution 表示扰动解的分量
$\dot{y}$	Time derivative of $y$ 的时间 $y$ 导数
$\overline{y}$	Mean orbital elements or average values 平均轨道元素或平均值
$\Delta y$	Difference (See Eq. (87)) 差异（见方程（87））
$\delta y$	Change relative to the nominal value $y_{\text{o}}$ (See Eqs. (77), (88), and (90)) 相对于标称值 $y_{\text{o}}$ 的变化（参见方程。（77）、（88）和（90））
$\delta_{\!\Delta}y$	See Eq. (96) 见方程（96）
$\mu\,(=GM_{\text{e}})$	Gravitational constant of the Earth 地球的引力常数
$\mu_{2}\,(=GM_{\text{s}})$	Gravitational constant of the Sun 太阳的引力常数
$\mu_{3}\,(=GM_{\text{m}})$	Gravitational constant of the Moon 月球的引力常数
$R_{\text{e}}$	Equatorial radius of the Earth 地球的赤道半径
$J_{2}$	Second zonal harmonic coefficient 二次区域谐波系数
$\epsilon$	Obliquity of the ecliptic 黄道的倾角
$U$	Gravitational potential 重力势
$U_{0}$	Central gravitational potential 中心引力势
$\mathcal{R}$	Perturbative potential 扰动电位
$P_{N}(x)$	Legendre polynomial of degree $N$ 度的勒让德多项式 $N$
$\mathcal{N}$	Truncation degree 截断度
$R_{x}(\gamma)$ , $R_{z}(\gamma)$ $R_{x}(\gamma)$ 、 $R_{z}(\gamma)$	Rotation matrices about the $x$ and $z$ axes by angle $\gamma$ 按角度 $\gamma$ 绕 $x$ 和 $z$ 轴的旋转矩阵
$f_{0}$ , $f_{1}$ $f_{0}$ 、 $f_{1}$	Zeroth and first order of the function $f$ 函数 $f$ 的零阶和一阶
$\delta_{\iota}$	Takes values 0 or 1 (see Eqs. (27) and (71)) 取值 0 或 1（参见 Eqs。（27）和（71））
$r$	Geocentric satellite distance 地心卫星距离
$\mathbf{r}$	Geocentric satellite position vector 地心卫星位置矢量
$\mathbf{\dot{r}}$	Geocentric satellite velocity vector 地心卫星速度矢量
$\mathbf{\hat{r}}$ , $\mathbf{\hat{r}_{2}}$ , $\mathbf{\hat{r}_{3}}$ $\mathbf{\hat{r}}$ 、 $\mathbf{\hat{r}_{2}}$ 、 $\mathbf{\hat{r}_{3}}$	Unit vectors for satellite, the Sun, and the Moon, respectively 分别用于卫星、太阳和月亮的单位向量
$\psi_{2},\psi_{3}$	Geocentric angles between satellite and the Sun, and the Moon 卫星与太阳和月球之间的地心角
$\varphi$	Geocentric latitude in Earth-fixed coordinate system 地球固定坐标系中的地心纬度
$r_{2}$	Sun-Earth distance 日地距离
$\overline{r}_{2}$	Mean Sun-Earth distance 平均日地距离
$\text{u}_{2}$	Ecliptic longitude of the Sun 太阳的黄道经度
$u_{2}$	Mean longitude of the Sun 太阳的平均经度
$U_{2}$	Defined as $(u_{2}-\overline{\Omega})$ 定义为 $(u_{2}-\overline{\Omega})$
$r_{3}$	Earth-Moon distance 地月距离
$\overline{r}_{3}$	Mean Earth-Moon distance 平均地月距离
$\overline{i}_{3}$	Mean inclination of the Moon’s orbit 月球轨道的平均倾角
$\overline{\Omega}_{3}$	Secular variation in the Moon’s longitude of ascending node 月球上升交点经度的长期变化
$\overline{M}_{3}$	Secular variation in the Moon’s mean anomaly 月球平均距平的长期变化
$\text{u}_{3}$	Latitude argument of the Moon 月球的纬度论点
$u_{3}$	Secular variation in the Moon’s latitude argument 月球纬度论点的长期变化
$\Delta u_{3}$	Defined as $(u_{3}-\overline{M}_{3})$ 定义为 $(u_{3}-\overline{M}_{3})$
$\theta_{3}$	Defined as $(\overline{\Omega}-\overline{\Omega}_{3})$ 定义为 $(\overline{\Omega}-\overline{\Omega}_{3})$
$n$	Mean motion of satellite with an orbit period of 3.64 days 轨道周期为3.64天的卫星平均运动
$n_{\Omega}^{\prime}$	Rate of change of $\Omega$ 变化 $\Omega$ 率
$n_{\lambda}^{\prime}$	Rate of change of $\lambda$ 变化 $\lambda$ 率
$n_{2}$	Mean motion of the Earth with a period of 365.2564 days 地球的平均运动周期为 365.2564 天
$n_{3}$	Rate of change of $u_{3}$ with a period of 27.21 days (draconic month) $u_{3}$ 周期为 27.21 天的变化率（draconic month）
$n_{\!M_{3}}$	Rate of change of $\overline{M}_{3}$ with a period of 27.55 days (anomalistic month) 周期为 27.55 天（异常月份）的变化 $\overline{M}_{3}$ 率
$n_{\Omega_{3}}$	Rate of change of $\overline{\Omega}_{3}$ with a period of 18.6 years 18.6 年的变化 $\overline{\Omega}_{3}$ 率
$\Delta n_{3}\,(=n_{3}-n_{\!M_{3}})$	Rate of change of $\Delta u_{3}$ with a period of 6.0 years 6.0 年期间的变化 $\Delta u_{3}$ 率
$n_{\theta_{3}}\,(=n_{\Omega}^{\prime}-n_{\Omega_{3}})$	Rate of change of $\theta_{3}$ 变化 $\theta_{3}$ 率
$n_{U_{2}}\,(=n_{2}-n_{\Omega}^{\prime})$	Rate of change of $U_{2}$ 变化 $U_{2}$ 率
$L_{ij}$	Arm-length formed by SC $i$ and SC $j$ 由 SC $i$ 和 SC $j$ 形成的臂长
$v_{ij}$	Relative line-of-sight velocity (rate of change of $L_{ij}$ ) 相对视距速度（变化率 $L_{ij}$ ）
$\alpha_{k}$	Breathing angle at SC $k$ SC $k$ 处的呼吸角度
$L_{\text{o}}$ , $v_{\text{o}}$ , $\alpha_{\text{o}}$ $L_{\text{o}}$ 、 $v_{\text{o}}$ 、 $\alpha_{\text{o}}$	Nominal values of arm-length, relative velocity, and breathing angle 臂长、相对速度和呼吸角的标称值
$f^{\sigma(N)}_{[q]}(i),f^{\sigma(N)}_{[p,\,q]}(i)$	Inclination functions in the lunar perturbation solution 月球扰动解中的倾角函数

Appendix B Model derivation and verification
附录 B 模型推导和验证

B.1 Derivation and motivation of Eq. (54)
B.1 方程（54）的推导和动机

For lunar point-mass perturbation, the perturbative potential $\mathcal{R}_{\text{m}}$ is represented as [35]
对于月球点质量扰动，扰动势 $\mathcal{R}_{\text{m}}$ 表示为 [35]

	$\displaystyle\mathcal{R}_{\text{m}}$	$\displaystyle=\mu_{3}\left(\frac{1}{\|\mathbf{r}-\mathbf{r_{3}}\|}-\frac{\mathbf% {r_{3}}}{r_{3}^{3}}\cdot\mathbf{r}\right)$		(116)
		$\displaystyle=\mu_{3}\left(\frac{1}{\sqrt{{r}^{2}-2\,r\,r_{3}\cos\psi_{3}+r_{3% }^{2}}}-\frac{r}{r_{3}^{2}}\cos\psi_{3}\right),$		(117)

where $r=\frac{a(1-e^{2})}{1+e\cos\nu}$ , $r_{3}=\frac{a_{3}(1-e_{3}^{2})}{1+e_{3}\cos\nu_{3}}$ , and $\cos\psi_{3}=\mathbf{\hat{r}_{3}}\cdot\mathbf{\hat{r}}$ . $\mathcal{R}_{\text{m}}$ , described by Eq. (117), exhibits a square root form, introducing challenges in solving the Lagrange perturbation equations. This complexity can be circumvented by expressing $\frac{1}{|\mathbf{r}-\mathbf{r}_{3}|}$ as an expansion of Legendre polynomials:
其中 $r=\frac{a(1-e^{2})}{1+e\cos\nu}$ ， $r_{3}=\frac{a_{3}(1-e_{3}^{2})}{1+e_{3}\cos\nu_{3}}$ ，和 $\cos\psi_{3}=\mathbf{\hat{r}_{3}}\cdot\mathbf{\hat{r}}$ . $\mathcal{R}_{\text{m}}$ ，由方程（117）描述，表现出平方根形式，给求解拉格朗日微扰方程带来了挑战。这种复杂性可以通过表示 $\frac{1}{|\mathbf{r}-\mathbf{r}_{3}|}$ 为勒让德多项式的扩展来规避：

	$\displaystyle\frac{1}{\|\mathbf{r}-\mathbf{r_{3}}\|}$	$\displaystyle=\frac{1}{r_{3}}\left[1-2\left(\frac{r}{r_{3}}\right)\cos\psi_{3}% +\left(\frac{r}{r_{3}}\right)^{2}\right]^{-1/2}$		(118)
		$\displaystyle=\frac{1}{r_{3}}\sum_{N=0}^{\infty}\left(\frac{r}{r_{3}}\right)^{% N}P_{N}(\cos\psi_{3}).$		(119)

Further substituting Eq. (119) into Eq. (117) and removing the first term $\frac{1}{r_{3}}$ (as $\frac{\partial}{\partial\sigma}(\frac{1}{r_{3}})=0$ after substitution into Lagrange’s equations) yields
进一步将方程（119）代入方程（117）并删除第一项 $\frac{1}{r_{3}}$ （如 $\frac{\partial}{\partial\sigma}(\frac{1}{r_{3}})=0$ 代入拉格朗日方程后）得到

\displaystyle\mathcal{R}_{\text{m}}=\frac{\mu_{3}}{r_{3}}\sum_{N=2}^{\infty}% \left(\frac{r}{r_{3}}\right)^{N}P_{N}(\cos\psi_{3}).

(120)

This formulation proves more advantageous for solving the Lagrange equations than Eq. (117). Based on estimated magnitudes and validation through numerical simulations, we set the maximum degree of Legendre polynomials in $\mathcal{R}_{\text{m}}$ at $N=6$ . Additionally, the solar potential $\mathcal{R}_{\text{s}}$ resembles Eq. (120), with the maximum degree set at $N=2$ .
事实证明，这个公式比方程（117）更有利于求解拉格朗日方程。根据估计的大小和通过数值模拟的验证，我们将勒让德多项式的最大次数设置为 $\mathcal{R}_{\text{m}}$ 。 $N=6$ 此外，太阳势 $\mathcal{R}_{\text{s}}$ 类似于方程（120），最大度数设置为 $N=2$ 。

B.2 Derivation of Eqs. (73)-(75)
乙.2 方程的推导。（73）-（75）

The analytical expressions for Eqs. (73)–(75) can be derived by applying perturbation methods to solve Eq. (70). To enhance the accuracy of the analytical solution, it is more advantageous to use the mean orbital elements $\overline{\sigma}(t)$ [34, 35], corresponding to a long-term precessing elliptical orbit, rather than the Keplerian orbit $\sigma^{(0)}(t)$ as the reference solution. Consequently, the perturbation solution’s form (73) is reformulated as
Eqs 的解析表达式。（73）–（75）可以通过应用微扰方法求解方程（70）来推导。为了提高解析解的准确性，使用对应于长期进动椭圆轨道的平均轨道单元 $\overline{\sigma}(t)$ [34， 35] 而不是开普勒轨道 $\sigma^{(0)}(t)$ 作为参考解更有利。因此，扰动解的形式（73）被重新表述为

\displaystyle\sigma(t)=\overline{\sigma}(t)+[\sigma_{l[l]}(t)+\sigma_{s}(t)]

(121)

with 跟

\displaystyle\overline{\sigma}(t)=\overline{\sigma}^{(0)}(t)+\Delta\sigma_{c}(% t)+\Delta\sigma_{l[c]}(t),

(122)

and 和

	$\displaystyle\overline{\sigma}^{(0)}(t)=\overline{\sigma}_{0}+\delta_{\iota}% \overline{n}(t-t_{0}),$		(123)
	$\displaystyle\overline{\sigma}_{0}=\sigma_{0}-[\sigma_{l[l]}(t_{0})+\sigma_{s}% (t_{0})],$		(124)

where $\overline{\sigma}^{(0)}$ represents the unperturbed secular variations, $\overline{\sigma}_{0}$ is the initial mean elements, $\Delta\sigma_{c}(t):=\sigma_{c}(t)-\sigma_{c}(t_{0})$ , and $\Delta\sigma_{l[c]}(t):=\sigma_{l[c]}(t)-\sigma_{l[c]}(t_{0})$ . Notably, $\sigma_{l[c]}(t)$ is incorporated into $\overline{\sigma}(t)$ , considering its short-term behavior akin to secular variation.
其中 $\overline{\sigma}^{(0)}$ 表示未受扰动的长期变化， $\overline{\sigma}_{0}$ 是初始均值元素 $\Delta\sigma_{c}(t):=\sigma_{c}(t)-\sigma_{c}(t_{0})$ ，和 $\Delta\sigma_{l[c]}(t):=\sigma_{l[c]}(t)-\sigma_{l[c]}(t_{0})$ 。值得注意的是， $\sigma_{l[c]}(t)$ 考虑到其类似于长期变化的短期行为，它被纳入 $\overline{\sigma}(t)$ 。

In relation to the left-side partitioning of Eq. (70) concerning $\sigma(t)$ , the function $f_{1}$ on the right is similarly decomposed into
关于方程（70）的左侧划分 $\sigma(t)$ ，右侧的函数 $f_{1}$ 同样被分解为

\displaystyle f_{1}=f_{1c}+f_{1l[c]}+f_{1l[l]}+f_{1s}.

(125)

$f_{1c}$ depends solely on $\overline{a}$ , $\overline{\xi}$ , $\overline{\eta}$ , and $\overline{i}$ . Both $f_{1l[c]}$ and $f_{1l[l]}$ involve trigonometric functions with arguments related to slow variables, such as $\Omega_{3}$ with an 18.6-year period and $u_{3}$ with a 27.21-day period, while $f_{1s}$ incorporates the fast variable $\lambda$ , which has a 3.64-day period, as the argument. The decomposition in Eq. (125) is achieved through averaging [53, 47, 16], where, for instance, $f_{1s}$ is obtained via
$f_{1c}$ 仅依赖于 $\overline{a}$ 、 $\overline{\xi}$ 、 $\overline{\eta}$ 和 $\overline{i}$ 。和 $f_{1l[l]}$ 都 $f_{1l[c]}$ 涉及三角函数，其参数与慢变量相关，例如 $\Omega_{3}$ 18.6 年周期和 $u_{3}$ 27.21 天周期，同时 $f_{1s}$ 包含快速变量 $\lambda$ ，其周期为 3.64 天，作为参数。方程（125）中的分解是通过平均[53， 47， 16]实现的，例如 $f_{1s}$ ，通过求法得到

\displaystyle f_{1s}=f_{1}-\left\langle f_{1}\right\rangle_{\lambda},\qquad% \left\langle f_{1}\right\rangle_{\lambda}:=\frac{1}{2\pi}\int_{0}^{2\pi}f_{1}% \,d\lambda.

(126)

A similar averaging over slow variables is applied to derive $f_{1c}$ , $f_{l[c]}$ , and $f_{l[l]}$ . Moreover, given that the TianQin orbits are nearly circular with $\overline{e}\simeq 0.0005$ [24], the terms on the right side of Eq. (125) consider only the leading-order effects of eccentricity for simplicity.
对慢速变量进行类似的平均应用于派生 $f_{1c}$ 、 $f_{l[c]}$ 和 $f_{l[l]}$ 。此外，鉴于天琴轨道与 [24] 几乎是圆形 $\overline{e}\simeq 0.0005$ 的，为简单起见，方程（125）右侧的项只考虑了偏心率的超前序效应。

By inserting the formal solution (121) into both sides of Eq. (70) and conducting a Taylor expansion around $\overline{\sigma}(t)$ , the comparison of coefficients for the same powers ( $\varepsilon^{0}$ , $\varepsilon^{1}$ , $\varepsilon^{2}$ , $\cdots$ ) yields [35]
通过将形式解（121）代入方程（70）的两侧并围绕 $\overline{\sigma}(t)$ 进行泰勒展开，相同幂（ $\varepsilon^{0}$ ， $\varepsilon^{1}$ ， $\varepsilon^{2}$ $\cdots$ ，）的系数比较得到 [35]

	$\displaystyle\overline{\sigma}^{(0)}(t)=\int^{t}_{t_{0}}[f_{0}]_{\overline{a}}% \,dt=\overline{\sigma}_{0}+\delta_{\iota}\overline{n}(t-t_{0}),$		(127)
	$\displaystyle\sigma_{c}^{(1)}(t)=\int^{t}[f_{1c}]_{\overline{\sigma}}\,dt,$		(128)
	$\displaystyle\sigma_{l[c]}^{(1)}(t)=\int^{t}[f_{1l[c]}]_{\overline{\sigma}}\,dt,$		(129)
	$\displaystyle\sigma_{l[l]}^{(1)}(t)=\int^{t}[f_{1l[l]}]_{\overline{\sigma}}\,dt,$		(130)
	$\displaystyle\sigma_{s}^{(1)}(t)=\int^{t}[\frac{\partial f_{0}}{\partial a}a_{% s}^{(1)}+f_{1s}]_{\overline{\sigma}}\,dt,$		(131)
	$\displaystyle\sigma_{c}^{(2)}(t)=\int^{t}\biggl{[}\frac{1}{2}\frac{\partial^{2% }\!f_{0}}{\partial a^{2}}[a_{s}^{(1)}]^{2}_{c}$
	$\displaystyle\qquad\quad\;\;\,\,+\biggl{(}\sum_{j=1}^{6}\frac{\partial f_{1}}{% \partial\sigma_{j}}[\sigma_{s}^{(1)}+\sigma_{l[l]}^{(1)}]_{j}\biggr{)}_{\!\!c}% \biggr{]}_{\!\overline{\sigma}}\,dt,$		(132)
	$\displaystyle\cdots.$

The superscript in parentheses denotes the order of the perturbation solution; in this paper, we focus on the first-order solution. Utilizing Eqs. (128)-(131), we derive explicit expressions for the four components of $\sigma^{(1)}(t)$ , presented in Appendix C. Particularly, for $\xi_{c}$ , $\xi_{l[c]}$ , $\eta_{c}$ , and $\eta_{l[c]}$ , it is more reasonable to directly solve the oscillation equations they satisfy (see Eqs. (191) and (192)) [35]; detailed derivations are provided in Appendix C.4.
括号中的上标表示扰动解的顺序;在本文中，我们重点介绍一阶解。利用方程。（128）-（131）中，我们推导出了附录 C 中所示的 $\sigma^{(1)}(t)$ 的四个组成部分的显式表达式。特别是，对于 $\xi_{c}$ 、、 $\xi_{l[c]}$ $\eta_{c}$ 和 $\eta_{l[c]}$ ，直接求解它们满足的振荡方程更合理（参见 Eqs。（191）和（192）） [35];附录 C.4 提供了详细的推导。

Note that in Eqs. (127)-(132), $\sigma$ on the right-hand side all take the form $\overline{\sigma}(t)$ defined in Eq. (122). For $i$ , $\Omega$ , and $\lambda$ , they are embedded in the trigonometric functions of $f_{1}$ . To enable integrable solutions, $\sigma_{l[c]}(t)$ in $\overline{\sigma}(t)$ is approximated as a linear term with a rate of change $\widetilde{n}_{\sigma}$ ,
请注意，在 Eqs.（127）-（132）， $\sigma$ 右侧都采用方程（122）中定义的形式 $\overline{\sigma}(t)$ 。对于 $i$ 、 $\Omega$ 和 $\lambda$ ，它们嵌入在的 $f_{1}$ 三角函数中。为了实现可积解， $\sigma_{l[c]}(t)$ in $\overline{\sigma}(t)$ 近似为具有变化 $\widetilde{n}_{\sigma}$ 率的线性项，

\displaystyle\widetilde{n}_{\sigma}=\frac{\left\langle\Delta\sigma_{l[c]}(t)% \right\rangle_{\tau}}{\tau},\qquad\left\langle f(t)\right\rangle_{\tau}:=\frac% {1}{\tau}\int_{t_{0}}^{t_{0}+\tau}f(t)\,dt,

(133)

where $\tau$ denotes the duration, implying
其中 $\tau$ 表示持续时间，表示

\overline{\sigma}(t)\simeq\left\{\begin{array}[]{ll}\overline{\sigma}_{0}+n_{% \sigma}^{\prime}(t-t_{0}),&\sigma\in\{\Omega,\lambda\},\vskip 3.0pt plus 1.0pt% minus 1.0pt\\ \left\langle i(t)\right\rangle_{\tau},&\sigma=i,\\ \end{array}\right.

(134)

with 跟

\displaystyle n_{\sigma}^{\prime}:=\delta_{\iota}\overline{n}+n_{\sigma}+% \widetilde{n}_{\sigma}.

(135)

Here, $n_{\sigma}$ represents the rate of change for $\sigma_{c}(t)$ , given by $n_{\sigma}=n_{\sigma\text{s}}+n_{\sigma\text{m}}+n_{\sigma\!_{J_{2}}}$ . Additionally, $\overline{i}(t)$ is approximated as the mean value $\left\langle i(t)\right\rangle_{\tau}$ , due to its small secular variation.
其中， $n_{\sigma}$ 表示 $\sigma_{c}(t)$ 的变化率，由给出。 $n_{\sigma}=n_{\sigma\text{s}}+n_{\sigma\text{m}}+n_{\sigma\!_{J_{2}}}$ 此外，由于其长期变化较小， $\overline{i}(t)$ 因此近似为平均值 $\left\langle i(t)\right\rangle_{\tau}$ 。

B.3 Verification of Eqs. (73)-(75)
乙.3 方程验证。（73）-（75）

To validate the derived analytical solution for satellite orbits, we conduct high-precision numerical orbit simulations using the NASA General Mission Analysis Tool (GMAT) [54]. The force models, consistent with those in Ref. [24], include the point-mass gravity fields of the Moon, Sun, and solar system planets (the ephemeris DE421), a $10\times 10$ spherical-harmonic model of the Earth’s gravity field (JGM-3), and the first-order relativistic correction. Non-gravitational perturbations, such as solar radiation pressure, are omitted as the satellites are drag-free controlled. Additionally, an adaptive step, ninth-order Runge-Kutta integrator with eighth-order error control (RungeKutta89) is employed, with the maximum integration step size set to 45 minutes. Initial orbital elements for the test orbits are detailed in Table 4. Orbit-1 corresponds to the nominal orbit of the TianQin satellite. In addition, three cases with different inclinations are considered to facilitate a more comprehensive validation, considering that the inclination is a crucial parameter in the analytical solution.
为了验证卫星轨道的推导解析解，我们使用 NASA 通用任务分析工具（GMAT） [54] 进行了高精度数值轨道模拟。力模型与参考文献[24]中的模型一致，包括月球、太阳和太阳系行星的点质量重力场（星历表DE421）、地球引力场的 $10\times 10$ 球谐模型（JGM-3）和一阶相对论校正。由于卫星是无阻力控制的，因此省略了非引力扰动，例如太阳辐射压力。此外，采用自适应步长、具有八阶误差控制的九阶 Runge-Kutta 积分器（RungeKutta89），最大积分步长设置为 45 分钟。表 4 中详细说明了测试轨道的初始轨道元件。Orbit-1 对应于天琴卫星的标称轨道。此外，三种具有不同倾向的情况是考虑到倾斜度是解析解中的关键参数，因此有助于进行更全面的验证。

Table 4: Initial orbital elements of the test orbits in the J2000-based Earth-centered ecliptic coordinate system at the epoch 1 Jan, 2034, 00:00:00 UTC.
表 4：在 2034 年 1 月 1 日 00：00：00 UTC 纪元的基于 J2000 的地心黄道坐标系中测试轨道的初始轨道元素。

Test orbits 测试轨道	$a_{0}$ (km) $a_{0}$ （公里）	$i_{0}$ (^∘) $i_{0}$ （^∘）	$\Omega_{0}$ (^∘) $\Omega_{0}$ （^∘）	$\nu_{0}$ (^∘) $\nu_{0}$ （^∘）
Orbit-1 轨道-1	100 000	94.7	210.4	60
Orbit-2 轨道 2	100 000	65.0	210.4	60
Orbit-3 轨道 3	100 000	35.0	210.4	60
Orbit-4 轨道 4	100 000	05.0	210.4	60

The comparison between analytical and numerical orbits reveals the errors $\Delta\sigma(t):=\sigma_{\text{ana}}(t)-\sigma_{\text{num}}(t)$ ( $\sigma\in\{a,e,i,\Omega,\omega,\lambda\}$ ). Statistical results, shown in Table 5, demonstrate that analytical expressions for $e$ , $i$ , $\Omega$ , and $\lambda$ are in good agreement with numerical simulations, with the relative mean deviation of $e$ being less than $6\%$ , and long-term deviations for $i$ , $\Omega$ , and $\lambda$ being small. In addition, there are relatively large errors in $a$ and $\omega$ , with the latter having a minor influence on the constellation stability (see, e.g., Eq. (89)). Table 5 also includes a comparison of satellite positions, denoted as $|\Delta\mathbf{r}(t)|:=|\mathbf{r}_{\text{ana}}(t)-\mathbf{r}_{\text{num}}(t)|$ . For the TianQin orbit, the average and maximum deviations over a 5-year period are approximately 87 km and 210 km, respectively.
解析轨道和数值轨道之间的比较揭示了误差 $\Delta\sigma(t):=\sigma_{\text{ana}}(t)-\sigma_{\text{num}}(t)$ （ $\sigma\in\{a,e,i,\Omega,\omega,\lambda\}$ ）。表 5 所示的统计结果表明，、、和 $i$ $\Omega$ $\lambda$ 的 $e$ 解析表达式与数值模拟非常一致，相对平均偏差 $e$ 小于 $6\%$ ，而 $i$ 、 $\Omega$ 和 $\lambda$ 的长期偏差很小。此外，在和 $\omega$ 中 $a$ 还存在相对较大的误差，后者对星座稳定性的影响很小（参见方程（89））。表 5 还包括卫星位置的比较，表示为 $|\Delta\mathbf{r}(t)|:=|\mathbf{r}_{\text{ana}}(t)-\mathbf{r}_{\text{num}}(t)|$ 。对于天琴轨道，5 年期间的平均和最大偏差分别约为 87 公里和 210 公里。

Table 5: Comparison of the analytical solution for satellite orbits with numerical simulations over 5 years, indicating mean errors (and maximum errors). To address the secular variation of

\Omega

, the second-order solution

\Omega_{c}^{(2)}(t)

obtained from Eq. (132) has been incorporated into the analytical solution. When propagating

\lambda

and, consequently,

\mathbf{r}

using Eqs. (73)-(75), the mean semimajor axes of the numerical orbits were employed.
表 5：卫星轨道的解析解与 5 年内的数值模拟的比较，表明平均误差（和最大误差）。为了解决

\Omega

的长期变化，从方程（132）得到的二阶解

\Omega_{c}^{(2)}(t)

已被纳入解析解中。当传播

\lambda

并因此

\mathbf{r}

使用 Eqs.（73）-（75），采用了数值轨道的平均半长轴。

Test orbits 测试轨道	$\Delta a$ (km) $\Delta a$ （公里）	$\Delta e/e$ (%) $\Delta e/e$ （%）	$\Delta i$ (^∘) $\Delta i$ （^∘）	$\Delta\Omega$ (^∘) $\Delta\Omega$ （^∘）	$\Delta\omega$ (^∘) $\Delta\omega$ （^∘）	$\Delta\lambda$ (^∘) $\Delta\lambda$ （^∘）	$\|\Delta\mathbf{r}\|$ (km) $\|\Delta\mathbf{r}\|$ （公里）
Orbit-1 轨道-1	$-$ 1.4 ( $-$ 16.1) $-$ 1.4 （ $-$ 16.1）	$+$ 5.7 ( $+$ 7.3) $+$ 5.7 （ $+$ 7.3）	$-$ 0.00 ( $-$ 0.02) $-$ 0.00 （ $-$ 0.02）	$-$ 0.00 ( $-$ 0.01) $-$ 0.00 （ $-$ 0.01）	$-$ 1.1 ( $+$ 29.9) $-$ 1.1 （ $+$ 29.9）	$-$ 0.05 ( $-$ 0.11) $-$ 0.05 （ $-$ 0.11）	087 (210) 087 （210）
Orbit-2 轨道 2	$+$ 0.4 ( $-$ 15.5) $+$ 0.4 （ $-$ 15.5）	$+$ 1.2 ( $+$ 3.0) $+$ 1.2 （ $+$ 3.0 ）	$+$ 0.01 ( $+$ 0.06) $+$ 0,01 （ $+$ 0,06）	$+$ 0.08 ( $+$ 0.12) $+$ 0,08 （ $+$ 0,12）	$-$ 1.1 ( $+$ 70.9) $-$ 1.1 （ $+$ 70.9）	$+$ 0.01 ( $-$ 0.10) $+$ 0,01 （ $-$ 0,10）	126 (240)
Orbit-3 轨道 3	$+$ 1.9 ( $+$ 16.2) $+$ 1.9 （ $+$ 16.2 英寸）	$-$ 4.4 ( $-$ 8.5) $-$ 4.4 （ $-$ 8.5）	$+$ 0.04 ( $+$ 0.14) $+$ 0.04 （ $+$ 0.14）	$+$ 0.08 ( $+$ 0.14) $+$ 0,08 （ $+$ 0,14）	$-$ 1.6 ( $-$ 44.4) $-$ 1.6 （ $-$ 44.4）	$+$ 0.08 ( $+$ 0.18) $+$ 0,08 （ $+$ 0,18）	291 (502)
Orbit-4 轨道 4	$+$ 2.2 ( $+$ 18.2) $+$ 2.2 （ $+$ 18.2 页）	$-$ 1.0 ( $+$ 4.9) $-$ 1.0 （ $+$ 4.9）	$-$ 0.02 ( $-$ 0.05) $-$ 0,02 （ $-$ 0,05）	$-$ 0.38 ( $-$ 0.96) $-$ 0.38 （ $-$ 0.96）	$+$ 7.3 ( $+$ 19.9) $+$ 7.3 （ $+$ 19.9）	$+$ 0.42 ( $+$ 1.02) $+$ 0.42 （ $+$ 1.02）	161 (426)

For future improvements, potential dominant sources causing the aforementioned errors are briefly outlined as follows. Firstly, simplified analytical coordinates for the Sun and Moon (see Sec. III.2) were utilized, instead of higher-precision ones with multiple trigonometric corrections [35, 45]. Secondly, smaller perturbative effects, including those from other planets in the solar system and the nonspherical gravitational field of the Sun and Moon, were omitted in Eq. (52). Lastly, the second-order solution was lacking, and the next-leading-order eccentricity effect was neglected in Eq. (125), etc.
为了将来的改进，导致上述错误的潜在主要来源简要概述如下。首先，使用了太阳和月亮的简化解析坐标（见第III.2节），而不是具有多次三角校正的高精度坐标[35,45]。其次，方程（52）中省略了较小的扰动效应，包括来自太阳系中其他行星的扰动效应以及太阳和月亮的非球形引力场。最后，在方程（125）中忽略了二阶解，忽略了次导阶偏心率效应。

B.4 Model verification for the three indicators
乙.4 3个指标的模型验证

Appendix B.3 verifies the analytical solution for satellite orbits, focusing on individual satellites. Additionally, this subsection presents the verification of the analytical expressions for the kinematic indicators of the three-satellite constellation: $L_{ij}(t)$ , $v_{ij}(t)$ , and $\alpha_{k}(t)$ , derived from Eq. (76) (or Eqs. (73)-(75)).
附录 B.3 验证了卫星轨道的分析解，重点关注单个卫星。此外，本小节还介绍了对三卫星星座运动学指标的分析表达式的验证： $L_{ij}(t)$ 、 $v_{ij}(t)$ 、和 $\alpha_{k}(t)$ ，源自方程（76）（或方程。（73）-（75））。

The time evolution of these three indicators in both analytical and numerical models is plotted in the left panel of Fig. 3 for a representative set of initial orbital elements provided in Table 6. In the numerical model, the considered perturbations, integrator, and step size align with those detailed in Appendix B.3. The right panel illustrates the time evolution of the deviations between the analytical and numerical models for these three quantities. Figure 3 indicates that the analytical model can effectively capture the long-term variations in the indicators, while noticeable periodic deviations exist. Numerical simulation results suggest that these deviations primarily arise from approximations in the Sun and Moon analytical dynamical model. By employing higher-precision models for solar and lunar motion [35, 45], incorporating numerous trigonometric correction terms in the Sun’s ecliptic longitude $\text{u}_{2}$ and the Moon’s latitude argument $\text{u}_{3}$ , these deviations can be effectively reduced. On the other hand, this enhanced complexity presents challenges in analytically solving the Lagrange equations, as $\text{u}_{2}$ and $\text{u}_{3}$ themselves involve trigonometric functions (see Eqs. (57) and (59)). Additionally, beyond the orbits specified in Table 6, the analytical model has been validated on two additional sets: nominal orbits (with SC1’s initial elements matching those of Orbit-1 in Table 4) and optimized orbits (refer to Table 3 in [24]), yielding consistent results.
这三个指标在分析模型和数值模型中的时间演变如图左侧面板所示。3 表示表6中提供的一组具有代表性的初始轨道元素。在数值模型中，考虑的扰动、积分器和步长与附录 B.3 中详述的一致。右图显示了这三个量的分析模型和数值模型之间偏差的时间演变。图3表明，分析模型可以有效地捕捉指标的长期变化，同时存在明显的周期性偏差。数值模拟结果表明，这些偏差主要源于太阳和月亮解析动力学模型中的近似值。通过采用更高精度的太阳和月球运动模型[35,45]，在太阳的黄道经度 $\text{u}_{2}$ 和月球的纬度参数 $\text{u}_{3}$ 中结合许多三角校正项，可以有效地减少这些偏差。另一方面，这种增加的复杂性给解析求解拉格朗日方程带来了挑战，因为 $\text{u}_{2}$ 拉 $\text{u}_{3}$ 格朗日方程本身涉及三角函数（参见方程。（57）和（59））。此外，除了表 6 中规定的轨道之外，分析模型还在另外两组上进行了验证：标称轨道（SC1 的初始元件与表 4 中 Orbit-1 的初始元件相匹配）和优化轨道（参见 [24] 中的表 3），产生了一致的结果。

Figure 3: Time evolution of three kinematic indicators for the analytical model compared to the numerical model. The left panel displays subplots illustrating variations in arm-lengths (

L_{ij}

), relative velocities (

v_{ij}

), and breathing angles (

\alpha_{k}

), respectively. Colors in each subplot represent the numerical model (blue for

L_{12}

v_{12}

, or

\alpha_{3}

; red for

L_{23}

v_{23}

, or

\alpha_{1}

; green for

L_{31}

v_{31}

, or

\alpha_{2}

), while black denotes the analytical model (only

v_{12}

is plotted for clarity in relative velocities). The right panel illustrates the temporal evolution of discrepancies between the analytical and numerical models for these three indicators.
图 3：与数值模型相比，解析模型的三个运动学指标的时间演变。左侧面板显示子图，分别说明了臂长（

L_{ij}

）、相对速度（

v_{ij}

）和呼吸角度（

\alpha_{k}

）的变化。每个子图中的颜色代表数值模型（蓝色代表、或 ;红色代表

L_{12}

L_{23}

、

v_{23}

或

\alpha_{1}

;绿色代表

L_{31}

、或

v_{31}

\alpha_{2}

），而黑色表示解析模型（仅

v_{12}

绘制为了相对速度的清晰度）。

\alpha_{3}

v_{12}

右图说明了这三个指标的解析模型和数值模型之间差异的时间演变。

Table 6: Initial orbital elements for simulated TianQin orbits in the J2000-based Earth-centered ecliptic coordinate system at the epoch of 22 May, 2034, 12:00:00 UTC. These initial elements deviate by approximately 1 km,

10^{-5}

, 0.2^∘, 0.2^∘, 0.1^∘, and 0.1^∘ in

a_{0}

e_{0}

i_{0}

\Omega_{0}

\omega_{0}

, and

\nu_{0}

, respectively, from the optimized orbits (cf. Table 3 in [24]). The subsequent orbital evolution is illustrated in Fig. 3.
表 6： 2034 年 5 月 22 日 12：00：00 UTC 纪元，基于 J2000 的地心黄道坐标系中模拟天琴轨道的初始轨道元素。这些初始元素在优化轨道上

a_{0}

e_{0}

i_{0}

\Omega_{0}

\omega_{0}

\nu_{0}

分别偏离了大约 1 km、

10^{-5}

、0.2^∘、0.2^∘、0.1^∘ 和 0.1^∘（参见 [24] 中的表 3）。随后的轨道演化如图 1 所示。3.

	$a_{0}$ (km) $a_{0}$ （公里）	$e_{0}$	$i_{0}$ (^∘) $i_{0}$ （^∘）
SC1	099 996.572 323 电话 099 996.572 323	0.000 440	94.897 997
SC2	100 010.400 095	0.000 010	94.904 363
SC3	099 992.041 899 电话 099 992.041 899	0.000 296	94.509 747
	$\Omega_{0}$ (^∘) $\Omega_{0}$ （^∘）	$\omega_{0}$ (^∘) $\omega_{0}$ （^∘）	$\nu_{0}$ (^∘) $\nu_{0}$ （^∘）
SC1	210.645 892	358.724 463	061.429 603 电话 061.429 603
SC2	210.240 199	359.900 000	180.130 706
SC3	210.644 582	359.901 624	299.812 164

Furthermore, the expressions (89), (97), and (101), derived from the series expansion of the three indicators, have been verified. The results suggest that to achieve a deviation magnitude similar to that before the series expansion, the second-order term $\delta L_{ij}(\delta\lambda(t)^{2})=\frac{\sqrt{3}}{4}a_{\text{o}}\delta% \lambda_{i}(t)\delta\lambda_{j}(t)-\frac{\sqrt{3}}{8}a_{\text{o}}[\delta% \lambda_{i}(t)^{2}+\delta\lambda_{j}(t)^{2}]$ in the arm-length, where $\delta\lambda(t)^{2}\sim\delta n^{2}\,t^{2}$ rapidly increases with time, needs to be taken into account. Notably, for the relative velocity and breathing angle, these second-order terms are both zero, $\delta v_{ij}(\delta\lambda(t)^{2})=\delta\alpha_{k}(\delta\lambda(t)^{2})=0$ .
此外，从三个指标的级数扩展中得出的表达式（89）、（97）和（101）也得到了验证。结果表明，为了获得与级数展开前相似的偏差幅度，需要考虑臂长中的二阶项，该 $\delta\lambda(t)^{2}\sim\delta n^{2}\,t^{2}$ 项 $\delta L_{ij}(\delta\lambda(t)^{2})=\frac{\sqrt{3}}{4}a_{\text{o}}\delta% \lambda_{i}(t)\delta\lambda_{j}(t)-\frac{\sqrt{3}}{8}a_{\text{o}}[\delta% \lambda_{i}(t)^{2}+\delta\lambda_{j}(t)^{2}]$ 随时间迅速增加。值得注意的是，对于相对速度和呼吸角，这些二阶项都是零。 $\delta v_{ij}(\delta\lambda(t)^{2})=\delta\alpha_{k}(\delta\lambda(t)^{2})=0$

Appendix C Explicit expressions for terms in Eq. (75)
附录 C 方程（75）中项的显式表达式

In this section, explicit expressions for each term of $\sigma^{(1)}(t)$ , as listed in Table 2, are presented. These expressions, categorized by the perturbations of the Sun, Moon, and Earth’s $J_{2}$ , are detailed in C.1, C.2, and C.3. Appendix C.4 provides the perturbation solutions for jointly solved $\xi$ and $\eta$ : $\xi_{c}$ , $\eta_{c}$ , $\xi_{l[c]}$ , and $\eta_{l[c]}$ , considering both solar and lunar perturbations. Note that in the subsequent expressions, the orbital elements $a$ , $\Omega$ , $\lambda$ , and $i$ , take the mean value $\overline{a}$ or the form defined by Eq. (134). Similarly, the Moon’s orbit inclination $i_{3}$ represents $\overline{i}_{3}$ .
在本节中，显示了表 2 中所列的每个术语 $\sigma^{(1)}(t)$ 的显式表达式。这些表达式按太阳、月亮和地球 $J_{2}$ 的扰动分类，在 C.1、C.2 和 C.3 中有详细说明。附录 C.4 提供了联合解决 $\xi$ 的 AND $\eta$ 的扰动解： $\xi_{c}$ 、、 $\eta_{c}$ $\xi_{l[c]}$ 和 $\eta_{l[c]}$ ，同时考虑了太阳和月球的扰动。请注意，在后续表达式中，轨道元素 $a$ 、 $\Omega$ 、 $\lambda$ 和 $i$ 采用平均值 $\overline{a}$ 或方程（134）定义的形式。同样，月球的轨道倾角 $i_{3}$ 表示 $\overline{i}_{3}$ 。

C.1 Solar perturbation solution $\sigma_{\text{s}}(t)$
C.1 太阳能扰动解决方案 $\sigma_{\text{s}}(t)$

The secular variation $\sigma_{\text{s}}^{c}(t)$ , long-period variation $\sigma_{\text{s}}^{l}(t)$ , and short-period variation $\sigma_{\text{s}}^{s}(t)$ within $\sigma_{\text{s}}(t)$ are presented as follows:
内部的长期变化 $\sigma_{\text{s}}^{c}(t)$ 、长期变化 $\sigma_{\text{s}}^{l}(t)$ 和短期变化 $\sigma_{\text{s}}^{s}(t)$ 表示 $\sigma_{\text{s}}(t)$ 如下：

(1) Secular terms with the form $\sigma_{\text{s}}^{c}(t)$ :
（1）形式 $\sigma_{\text{s}}^{c}(t)$ 为：

$\displaystyle a_{\text{s}}^{c}=$	$\displaystyle~{}i_{\text{s}}^{c}=0,$	(136)
$\displaystyle\Omega_{\text{s}}^{c}=$	$\displaystyle~{}n_{\Omega{\text{s}}}\,t,\qquad n_{\Omega{\text{s}}}:=-\frac{3}% {4}c_{\text{s}}\cos i,$	(137)
$\displaystyle\lambda_{\text{s}}^{c}=$	$\displaystyle~{}n_{\lambda\text{s}}\,t,\qquad n_{\lambda\text{s}}:=\frac{1}{8}% c_{\text{s}}(1-3\cos 2i),$	(138)

where $c_{\text{s}}:=\frac{\mu_{2}}{\overline{n}\,\overline{r}_{2}^{3}}$ . 其中 $c_{\text{s}}:=\frac{\mu_{2}}{\overline{n}\,\overline{r}_{2}^{3}}$ .

(2) Long-period terms with the form $\sigma_{\text{s}}^{l}(t)$ :
（2）格式为 $\sigma_{\text{s}}^{l}(t)$ ：

$\displaystyle a_{\text{s}}^{l[l]}=$	$\displaystyle~{}\xi_{\text{s}}^{l[l]}=\eta_{\text{s}}^{l[l]}=0,$	(139)
$\displaystyle i_{\text{s}}^{l[l]}=$	$\displaystyle~{}\frac{3}{8}\frac{c_{\text{s}}}{n_{U_{2}}}\sin i\cos 2U_{2},$	(140)
$\displaystyle\Omega_{\text{s}}^{l[l]}=$	$\displaystyle~{}\frac{3}{8}\frac{c_{\text{s}}}{n_{U_{2}}}\cos i\sin 2U_{2},$	(141)
$\displaystyle\lambda_{\text{s}}^{l[l]}=$	$\displaystyle~{}\frac{3}{16}\frac{c_{\text{s}}}{n_{U_{2}}}(\cos 2i-3)\sin 2U_{% 2},$	(142)

with 跟

	$\displaystyle U_{2}(t):=u_{2}(t)-\overline{\Omega}(t)=n_{U_{2}}(t-t_{0})+U_{2_% {0}},$		(143)
	$\displaystyle n_{U_{2}}=n_{2}-n_{\Omega}^{\prime},\qquad U_{2_{0}}=u_{2_{0}}-% \overline{\Omega}_{0},$		(144)

where $n_{\Omega}^{\prime}$ is the rate of change of $\overline{\Omega}(t)$ as defined in Eq. (134). These expressions reveal that solar perturbation induces general long-period variations in satellite orbital elements with an annual period tied to the solar apparent motion $n_{2}$ . The magnitudes of these variations are governed by $\overline{a}$ and $\overline{i}$ . $\xi_{\text{s}}^{l[l]}=\eta_{\text{s}}^{l[l]}=0$ arises from considering only the leading order $N=2$ within $P_{N}(\cos\psi_{2})$ (refer to Eq. (53)); when $N\geq 3$ , both $\xi$ and $\eta$ will exhibit periodic variations, as observed in the case of lunar perturbation (see Eqs. (168) and (169)).
其中 $n_{\Omega}^{\prime}$ 是方程（134）中定义的的变化 $\overline{\Omega}(t)$ 率。这些表达式表明，太阳扰动会引起卫星轨道元素的一般长周期变化，其年周期与太阳视运动 $n_{2}$ 有关。这些变化的大小由 $\overline{a}$ 和 $\overline{i}$ 控制。 $\xi_{\text{s}}^{l[l]}=\eta_{\text{s}}^{l[l]}=0$ 产生于仅考虑（ $P_{N}(\cos\psi_{2})$ 参见方程（53））中的前导顺序 $N=2$ ）;当 $N\geq 3$ 时，两者和 $\xi$ $\eta$ 都将表现出周期性变化，如在月球扰动的情况下观察到的那样（参见方程。（168）和（169））。

There are no special long-period variations in $a$ , $i$ , $\Omega$ , and $\lambda$ ,
$a$ 在、、 $i$ $\Omega$ 和 $\lambda$ 中没有特殊的长周期变化

\displaystyle\sigma_{\text{s}}^{l[c]}\equiv 0,\quad\text{for}\ \sigma\in\{a,i,% \Omega,\lambda\}.

(145)

For $\xi$ and $\eta$ , they exhibit special long-period variations coupled with lunar perturbation, as indicated in Eqs. (204) and (205).
对于 $\xi$ 和 $\eta$ ，它们表现出特殊的长周期变化和月球扰动，如方程所示。（204）和（205）。

(3) Short-period terms with the form $\sigma_{\text{s}}^{s}(t)$ :
（3）短期期限，表格 $\sigma_{\text{s}}^{s}(t)$ 如下：

$\displaystyle a_{\text{s}}^{s}$	$\displaystyle=\sum_{\begin{subarray}{c}p=-2\\ (\Delta p=2)\end{subarray}}^{2}a\times h^{a}_{(2,\,p)}\cos(2\,\lambda+p\,U_{2}),$	(146)
$\displaystyle\xi_{\text{s}}^{s}$	$\displaystyle=\sum_{\begin{subarray}{c}\kappa=1\\ (\Delta\kappa=2)\end{subarray}}^{3}\sum_{\begin{subarray}{c}p=-2\\ (\Delta p=2)\end{subarray}}^{2}h^{\xi}_{(\kappa,\,p)}\cos(\kappa\,\lambda+p\,U% _{2}),$	(147)
$\displaystyle\eta_{\text{s}}^{s}$	$\displaystyle=\sum_{\begin{subarray}{c}\kappa=1\\ (\Delta\kappa=2)\end{subarray}}^{3}\sum_{\begin{subarray}{c}p=-2\\ (\Delta p=2)\end{subarray}}^{2}h^{\eta}_{(\kappa,\,p)}\sin(\kappa\,\lambda+p\,% U_{2}),$	(148)
$\displaystyle i_{\text{s}}^{s}$	$\displaystyle=\sum_{\begin{subarray}{c}p=-2\\ (\Delta p=2)\end{subarray}}^{2}h^{i}_{(2,\,p)}\cos(2\,\lambda+p\,U_{2}),$	(149)
$\displaystyle\Omega_{\text{s}}^{s}$	$\displaystyle=\sum_{\begin{subarray}{c}p=-2\\ (\Delta p=2)\end{subarray}}^{2}h^{\Omega}_{(2,\,p)}\sin(2\,\lambda+p\,U_{2}),$	(150)
$\displaystyle\lambda_{\text{s}}^{s}$	$\displaystyle=\sum_{\begin{subarray}{c}p=-2\\ (\Delta p=2)\end{subarray}}^{2}h^{\lambda}_{(2,\,p)}\sin(2\,\lambda+p\,U_{2}),$	(151)

with 跟

\displaystyle h^{\sigma}_{(\kappa,\,p)}:=\frac{c_{\text{s}}}{\kappa\,n_{% \lambda}^{\prime}+p\,n_{U_{2}}}\times h^{\sigma}_{[\kappa,\,p]}(i)

(152)

for $\sigma\in\{a,\xi,\eta,i,\Omega\}$ , and for $\sigma\in\{a,\xi,\eta,i,\Omega\}$ 和

\displaystyle h^{\lambda}_{(\kappa,\,p)}:=\frac{3\,c_{\text{s}}\,n\times h^{% \lambda[a]}_{[\kappa,\,p]}(i)}{(\kappa\,n_{\lambda}^{\prime}+p\,n_{U_{2}})^{2}% }+\frac{c_{\text{s}}\times h^{\lambda[\lambda]}_{[\kappa,\,p]}(i)}{\kappa\,n_{% \lambda}^{\prime}+p\,n_{U_{2}}},

(153)

where $n_{\lambda}^{\prime}$ denotes the rate of change of $\overline{\lambda}(t)$ as defined in Eq. (134). In the specific case of solar perturbation alone, $n_{\lambda}^{\prime}=\overline{n}+n_{\lambda\text{s}}$ . The terms on the right side of Eq. (153) correspond to the integrals of the two terms in Eq. (131). The explicit forms of $h^{\sigma}_{[\kappa,\,p]}(i)$ are given by
其中 $n_{\lambda}^{\prime}$ 表示方程（134）中定义的的变化 $\overline{\lambda}(t)$ 率。在单独的太阳扰动的特定情况下， $n_{\lambda}^{\prime}=\overline{n}+n_{\lambda\text{s}}$ .方程（153）右侧的项对应于方程（131）中两项的积分。的 $h^{\sigma}_{[\kappa,\,p]}(i)$ 显式形式由下式给出

	$\displaystyle h^{a}_{[2,0]}=\frac{3}{2}\sin^{2}i,\qquad h^{a}_{[2,\pm 2]}=% \frac{3}{8}(3\mp 4\cos i+\cos 2i),$
	$\displaystyle h^{i}_{[2,0]}=\frac{3}{8}\sin 2i,\qquad h^{i}_{[2,\pm 2]}=\frac{% 3}{16}(\pm 2\sin i-\sin 2i),$
	$\displaystyle h^{\Omega}_{[2,0]}=\frac{3}{4}\cos i,\qquad h^{\Omega}_{[2,\pm 2% ]}=\frac{3}{8}(\pm 1-\cos i),$
	$\displaystyle h^{\lambda[a]}_{[2,0]}=-\frac{3}{4}\sin^{2}i,\qquad h^{\lambda[a% ]}_{[2,\pm 2]}=-\frac{3}{16}(3\mp 4\cos i+\cos 2i),$
	$\displaystyle h^{\lambda[\lambda]}_{[2,0]}=-\frac{3}{8}(3-\cos 2i),$
	$\displaystyle h^{\lambda[\lambda]}_{[2,\pm 2]}=-\frac{3}{16}(5\mp 6\cos i+\cos 2% i),$		(154)
	$\displaystyle h^{\xi}_{[1,0]}=\frac{1}{16}(7-15\cos 2i),$
	$\displaystyle h^{\xi}_{[1,\pm 2]}=\frac{3}{32}(7\mp 12\cos i+5\cos 2i),$
	$\displaystyle h^{\xi}_{[3,0]}=\frac{3}{8}\sin^{2}i,\qquad h^{\xi}_{[3,\pm 2]}=% \frac{3}{32}(3\mp 4\cos i+\cos 2i),$
	$\displaystyle h^{\eta}_{[1,0]}=-\frac{1}{16}(11-3\cos 2i),$
	$\displaystyle h^{\eta}_{[1,\pm 2]}=-\frac{3}{32}(11\mp 12\cos i+\cos 2i),$
	$\displaystyle h^{\eta}_{[3,0]}=\frac{3}{8}\sin^{2}i,\qquad h^{\eta}_{[3,\pm 2]% }=\frac{3}{32}(3\mp 4\cos i+\cos 2i).$

C.2 Lunar perturbation solution $\sigma_{\text{m}}(t)$
C.2 月球扰动解 $\sigma_{\text{m}}(t)$

The secular variation $\sigma_{\text{m}}^{c}(t)$ , special long-period variation $\sigma_{\text{m}}^{l[c]}(t)$ , general long-period variation $\sigma_{\text{m}}^{l[l]}(t)$ , and short-period variation $\sigma_{\text{m}}^{s}(t)$ within $\sigma_{\text{m}}(t)$ are shown as follows:
长期变化 $\sigma_{\text{m}}^{c}(t)$ 、特殊长期变化 $\sigma_{\text{m}}^{l[c]}(t)$ 、一般长期变化 $\sigma_{\text{m}}^{l[l]}(t)$ 和内部 $\sigma_{\text{m}}(t)$ 的短期变化 $\sigma_{\text{m}}^{s}(t)$ 如下所示：

(1) Secular terms with the form $\sigma_{\text{m}}^{c}(t)$ :
（1）形式 $\sigma_{\text{m}}^{c}(t)$ 为：

$\displaystyle a_{\text{m}}^{c}=$	$\displaystyle~{}i_{\text{m}}^{c}=0,$	(155)
$\displaystyle\Omega_{\text{m}}^{c}=$	$\displaystyle~{}n_{\Omega{\text{m}}}\,t,$	(156)
$\displaystyle\lambda_{\text{m}}^{c}=$	$\displaystyle~{}n_{\lambda\text{m}}\,t,$	(157)

where 哪里

$\displaystyle n_{\Omega{\text{m}}}:=$	$\displaystyle-\frac{3\,c_{\text{m}}^{(2)}}{16}\cos i\,(1+3\cos 2i_{3})-\frac{4% 5\,c_{\text{m}}^{(4)}}{32768}$
	$\displaystyle\times(9\cos i+7\cos 3i)(9+20\cos 2i_{3}+35\cos 4i_{3})$
	$\displaystyle-\frac{525\,c_{\text{m}}^{(6)}}{16777216}(50\cos i+45\cos 3i+33% \cos 5i)$
	$\displaystyle\times(50+105\cos 2i_{3}+126\cos 4i_{3}+231\cos 6i_{3}),$	(158)

$\displaystyle n_{\lambda\text{m}}:=$	$\displaystyle~{}\frac{1}{32}c_{\text{m}}^{(2)}(1-3\cos 2i)(1+3\cos 2i_{3})+% \frac{9\,c_{\text{m}}^{(4)}}{65536}$
	$\displaystyle\times(27+40\cos 2i-35\cos 4i)(9+20\cos 2i_{3}$
	$\displaystyle+35\cos 4i_{3})+\frac{75\,c_{\text{m}}^{(6)}}{33554432}(250+455% \cos 2i$
	$\displaystyle+294\cos 4i-231\cos 6i)(50+105\cos 2i_{3}$
	$\displaystyle+126\cos 4i_{3}+231\cos 6i_{3}),$	(159)

and $c_{\text{m}}^{(N)}:=\frac{\mu_{3}}{\overline{n}\,\overline{r}_{2}^{3}}(\frac{% \overline{a}}{\overline{r}_{2}})^{N-2}$ . 和 $c_{\text{m}}^{(N)}:=\frac{\mu_{3}}{\overline{n}\,\overline{r}_{2}^{3}}(\frac{% \overline{a}}{\overline{r}_{2}})^{N-2}$ .

(2) Special long-period terms with the form $\sigma_{\text{m}}^{l[c]}(t)$ :
（2）特殊长期条款，格式 $\sigma_{\text{m}}^{l[c]}(t)$ 为：

$\displaystyle a_{\text{m}}^{l[c]}$	$\displaystyle=0,$	(160)
$\displaystyle i_{\text{m}}^{l[c]}$	$\displaystyle=\sum_{\begin{subarray}{c}N=2\\ (\Delta N=2)\end{subarray}}^{\mathcal{N}}\sum_{q=1}^{N}f^{i(N)}_{(q)}\cos q\,% \theta_{3},$	(161)
$\displaystyle\Omega_{\text{m}}^{l[c]}$	$\displaystyle=\sum_{\begin{subarray}{c}N=2\\ (\Delta N=2)\end{subarray}}^{\mathcal{N}}\sum_{q=1}^{N}f^{\Omega(N)}_{(q)}\sin q% \,\theta_{3},$	(162)
$\displaystyle\lambda_{\text{m}}^{l[c]}$	$\displaystyle=\sum_{\begin{subarray}{c}N=2\\ (\Delta N=2)\end{subarray}}^{\mathcal{N}}\sum_{q=1}^{N}f^{\lambda(N)}_{(q)}% \sin q\,\theta_{3},$	(163)

with 跟

\displaystyle f^{\sigma(N)}_{(q)}:=\frac{c_{\text{m}}^{(N)}}{q\,n_{\theta_{3}}% }\times f^{\sigma(N)}_{[q]}(i),

(164)

and 和

	$\displaystyle\theta_{3}(t):=\overline{\Omega}(t)-\overline{\Omega}_{3}(t)=n_{% \theta_{3}}(t-t_{0})+\theta_{3_{0}},$		(165)
	$\displaystyle n_{\theta_{3}}=n_{\Omega}^{\prime}-n_{\Omega_{3}},\qquad\theta_{% 3_{0}}=\overline{\Omega}_{0}-\Omega_{3_{0}}.$		(166)

For explicit forms of $f^{\sigma(N)}_{[q]}(i)$ , see Appendix D.1.
有关的 $f^{\sigma(N)}_{[q]}(i)$ 明确形式，请参阅附录 D.1。

(3) General long-period terms with the form $\sigma_{\text{m}}^{l[l]}(t)$ :
（3）一般长期条款，格式 $\sigma_{\text{m}}^{l[l]}(t)$ 为：

\displaystyle a_{\text{m}}^{l[l]}

\displaystyle=0,

(167)

and 和

$\displaystyle\xi_{\text{m}}^{l[l]}$	$\displaystyle=\sum_{\begin{subarray}{c}N=3\\ (\Delta N=2)\end{subarray}}^{\mathcal{N}}\sum_{\begin{subarray}{c}p>0\\ (\Delta p=2)\end{subarray}}^{N}\sum_{q=-N}^{N}f^{\xi(N)}_{(p,\,q)}\cos(p\,u_{3% }+q\,\theta_{3}),$	(168)
$\displaystyle\eta_{\text{m}}^{l[l]}$	$\displaystyle=\sum_{\begin{subarray}{c}N=3\\ (\Delta N=2)\end{subarray}}^{\mathcal{N}}\sum_{\begin{subarray}{c}p>0\\ (\Delta p=2)\end{subarray}}^{N}\sum_{\begin{subarray}{c}q=-N\\ (q\neq 0)\end{subarray}}^{N}f^{\eta(N)}_{(p,\,q)}\sin(p\,u_{3}+q\,\theta_{3}),$	(169)
$\displaystyle i_{\text{m}}^{l[l]}$	$\displaystyle=\sum_{\begin{subarray}{c}N=2\\ (\Delta N=2)\end{subarray}}^{\mathcal{N}}\sum_{\begin{subarray}{c}p>0\\ (\Delta p=2)\end{subarray}}^{N}\sum_{\begin{subarray}{c}q=-N\\ (q\neq 0)\end{subarray}}^{N}f^{i(N)}_{(p,\,q)}\cos(p\,u_{3}+q\,\theta_{3}),$	(170)
$\displaystyle\Omega_{\text{m}}^{l[l]}$	$\displaystyle=\sum_{\begin{subarray}{c}N=2\\ (\Delta N=2)\end{subarray}}^{\mathcal{N}}\sum_{\begin{subarray}{c}p>0\\ (\Delta p=2)\end{subarray}}^{N}\sum_{q=-N}^{N}f^{\Omega(N)}_{(p,\,q)}\sin(p\,u% _{3}+q\,\theta_{3}),$	(171)
$\displaystyle\lambda_{\text{m}}^{l[l]}$	$\displaystyle=\sum_{\begin{subarray}{c}N=2\\ (\Delta N=2)\end{subarray}}^{\mathcal{N}}\sum_{\begin{subarray}{c}p>0\\ (\Delta p=2)\end{subarray}}^{N}\sum_{q=-N}^{N}f^{\lambda(N)}_{(p,\,q)}\sin(p\,% u_{3}+q\,\theta_{3}),$	(172)

with 跟

\displaystyle f^{\sigma(N)}_{(p,\,q)}:=\frac{c_{\text{m}}^{(N)}}{p\,n_{3}+q\,n% _{\theta_{3}}}\times f^{\sigma(N)}_{[p,\,q]}(i).

(173)

For explicit forms of $f^{\sigma(N)}_{[p,\,q]}(i)$ , see Appendix D.2.
有关的 $f^{\sigma(N)}_{[p,\,q]}(i)$ 明确形式，请参阅附录 D.2。

(4) Short-period terms with the form $\sigma_{\text{m}}^{s}(t)$ :
（4）短期条款，格式 $\sigma_{\text{m}}^{s}(t)$ 为：

$\displaystyle a_{\text{m}}^{s}$	$\displaystyle=~{}\sum_{N=2}^{\mathcal{N}}\sum_{\begin{subarray}{c}\kappa>0\\ (\Delta\kappa=2)\end{subarray}}^{N}\sum_{\begin{subarray}{c}p=-N\\ (\Delta p=2)\end{subarray}}^{N}\sum_{q=-N}^{N}a\,f^{a(N)}_{(\kappa,\,p,\,q)}% \cos\Lambda,$	(174)
$\displaystyle\xi_{\text{m}}^{s}$	$\displaystyle=\sum_{N=2}^{\mathcal{N}}\sum_{\begin{subarray}{c}\kappa>0\\ (\Delta\kappa=2)\end{subarray}}^{N+1}\sum_{\begin{subarray}{c}p=-N\\ (\Delta p=2)\end{subarray}}^{N}\sum_{q=-N}^{N}f^{\xi(N)}_{(\kappa,\,p,\,q)}% \cos\Lambda,$	(175)
$\displaystyle\eta_{\text{m}}^{s}$	$\displaystyle=\sum_{N=2}^{\mathcal{N}}\sum_{\begin{subarray}{c}\kappa>0\\ (\Delta\kappa=2)\end{subarray}}^{N+1}\sum_{\begin{subarray}{c}p=-N\\ (\Delta p=2)\end{subarray}}^{N}\sum_{q=-N}^{N}f^{\eta(N)}_{(\kappa,\,p,\,q)}% \sin\Lambda,$	(176)
$\displaystyle i_{\text{m}}^{s}$	$\displaystyle=\sum_{N=2}^{\mathcal{N}}\sum_{\begin{subarray}{c}\kappa>0\\ (\Delta\kappa=2)\end{subarray}}^{N}\sum_{\begin{subarray}{c}p=-N\\ (\Delta p=2)\end{subarray}}^{N}\sum_{q=-N}^{N}f^{i(N)}_{(\kappa,\,p,\,q)}\cos\Lambda,$	(177)
$\displaystyle\Omega_{\text{m}}^{s}$	$\displaystyle=\sum_{N=2}^{\mathcal{N}}\sum_{\begin{subarray}{c}\kappa>0\\ (\Delta\kappa=2)\end{subarray}}^{N}\sum_{\begin{subarray}{c}p=-N\\ (\Delta p=2)\end{subarray}}^{N}\sum_{q=-N}^{N}f^{\Omega(N)}_{(\kappa,\,p,\,q)}% \sin\Lambda,$	(178)
$\displaystyle\lambda_{\text{m}}^{s}$	$\displaystyle=\sum_{N=2}^{\mathcal{N}}\sum_{\begin{subarray}{c}\kappa>0\\ (\Delta\kappa=2)\end{subarray}}^{N}\sum_{\begin{subarray}{c}p=-N\\ (\Delta p=2)\end{subarray}}^{N}\sum_{q=-N}^{N}f^{\lambda(N)}_{(\kappa,\,p,\,q)% }\sin\Lambda,$	(179)

with $\Lambda:=\kappa\,\lambda+p\,u_{3}+q\,\theta_{3}$ , 其中 $\Lambda:=\kappa\,\lambda+p\,u_{3}+q\,\theta_{3}$ ，

\displaystyle f^{\sigma(N)}_{(\kappa,\,p,\,q)}:=\frac{c_{\text{m}}^{(N)}}{% \kappa\,n_{\lambda}^{\prime}+p\,n_{3}+q\,n_{\theta_{3}}}\times f^{\sigma(N)}_{% [\kappa,\,p,\,q]}(i)

(180)

for $\sigma\in\{a,\xi,\eta,i,\Omega\}$ , and for $\sigma\in\{a,\xi,\eta,i,\Omega\}$ 和

	$\displaystyle f^{\lambda(N)}_{(\kappa,\,p,\,q)}:=$	$\displaystyle~{}\frac{3\,c_{\text{m}}^{(N)}\,n}{(\kappa\,n_{\lambda}^{\prime}+% p\,n_{3}+q\,n_{\theta_{3}})^{2}}\times f^{\lambda(N)[a]}_{[\kappa,\,p,\,q]}(i)$
		$\displaystyle+\frac{c_{\text{m}}^{(N)}}{\kappa\,n_{\lambda}^{\prime}+p\,n_{3}+% q\,n_{\theta_{3}}}\times f^{\lambda(N)[\lambda]}_{[\kappa,\,p,\,q]}(i).$		(181)

For explicit forms of $f^{\sigma(N)}_{[\kappa,\,p,\,q]}(i)$ , see Appendix D.3.
有关 $f^{\sigma(N)}_{[\kappa,\,p,\,q]}(i)$ 的显式形式，请参阅附录 D.3。

C.3 Earth’s $J_{2}$ perturbation solution $\sigma_{\!J_{2}}(t)$
C.3 地球的 $J_{2}$ 扰动解 $\sigma_{\!J_{2}}(t)$

The Earth’s $J_{2}$ perturbation solution, $\sigma_{\!J_{2}\text{(eq)}}$ in the geocentric equatorial coordinate system, can be derived straightforwardly from Eqs. (55) and (128)-(131). For consistency with the solar and lunar perturbation solutions, the ecliptic representation $\sigma_{\!J_{2}}$ is essential. Given that the $J_{2}$ perturbation is significantly smaller than lunisolar perturbations, this study focuses on the secular components, including $i_{J2}^{c}$ , $\Omega_{J2}^{c}$ , and $\lambda_{J2}^{c}$ . Employing spherical trigonometry [35, 33] and variational method, we have
在地心赤道坐标系 $\sigma_{\!J_{2}\text{(eq)}}$ 中，地球的 $J_{2}$ 扰动解可以直接从 Eqs 推导出来。（55）和（128）-（131）。为了与太阳和月球扰动解保持一致，黄道表示 $\sigma_{\!J_{2}}$ 是必不可少的。 $J_{2}$ 鉴于扰动明显小于阴阳扰动，本研究侧重于长期分量，包括 $i_{J2}^{c}$ 、 $\Omega_{J2}^{c}$ 和 $\lambda_{J2}^{c}$ 。采用球面三角学 [35， 33] 和变分方法，我们有

$\displaystyle i_{J_{2}}^{c}$	$\displaystyle=\sin\Omega\sin\epsilon\,\Omega_{J_{2}\text{(eq)}}^{c}=n_{i_{J_{2% }}}t,$	(182)
$\displaystyle\Omega_{J_{2}}^{c}$	$\displaystyle=(\cos\epsilon+\cot i\cos\Omega\sin\epsilon)\,\Omega_{J_{2}\text{% (eq)}}^{c}=n_{\Omega_{J_{2}}}t,$	(183)
$\displaystyle\lambda_{J_{2}}^{c}$	$\displaystyle=\lambda_{J_{2}\text{(eq)}}^{c}-\frac{\cos\Omega}{\sin i}\sin% \epsilon\,\Omega_{J_{2}\text{(eq)}}^{c}=n_{\lambda_{J_{2}}}\,t,$	(184)

with 跟

$\displaystyle n_{i_{J_{2}}}$	$\displaystyle:=-\frac{3}{2}c_{J_{2}}\sin\Omega\sin\epsilon\cos i_{\text{eq}},$	(185)
$\displaystyle n_{\Omega_{J_{2}}}$	$\displaystyle:=-\frac{3}{2}c_{J_{2}}(\cos\epsilon+\cot i\cos\Omega\sin\epsilon% )\cos i_{\text{eq}},$	(186)
$\displaystyle n_{\lambda_{J_{2}}}$	$\displaystyle:=\frac{3}{2}c_{\small{J_{2}}}(1+2\cos 2i_{\text{eq}})+\frac{3}{2% }c_{J_{2}}\cos i_{\text{eq}}\frac{\cos\Omega}{\sin i}\sin\epsilon,$	(187)

where $\epsilon=23.439\,291^{\circ}$ is the obliquity of the ecliptic, $c_{\small{J_{2}}}:=\frac{\mu J_{2}R_{\text{e}}^{2}}{\overline{a}^{5}\,% \overline{n}}$ , and
其中 $\epsilon=23.439\,291^{\circ}$ 是黄道的倾角， $c_{\small{J_{2}}}:=\frac{\mu J_{2}R_{\text{e}}^{2}}{\overline{a}^{5}\,% \overline{n}}$ 和

\displaystyle\cos i_{\text{eq}}=\cos i\cos\epsilon-\sin i\cos\Omega\sin\epsilon.

(188)

C.4 Expressions for $\xi_{c}$ , $\eta_{c}$ , $\xi_{l[c]}$ , and $\eta_{l[c]}$
C.4 $\xi_{c}$ 、 $\eta_{c}$ 、 $\xi_{l[c]}$ 和 $\eta_{l[c]}$ 的表达式

As pointed out in Appendix B, deriving the “secular” perturbation solutions $\xi^{c}_{\text{new}}(t)$ and $\eta^{c}_{\text{new}}(t)$ ,
如附录 B 中所指出的，推导出“长期”扰动解 $\xi^{c}_{\text{new}}(t)$ 和 $\eta^{c}_{\text{new}}(t)$ ，

	$\displaystyle\xi^{c}_{\text{new}}(t)$	$\displaystyle:=\xi_{c}(t)+\xi_{l[c]}(t),$		(189)
	$\displaystyle\eta^{c}_{\text{new}}(t)$	$\displaystyle:=\eta_{c}(t)+\eta_{l[c]}(t),$		(190)

involves solving the oscillation equations:
涉及求解振荡方程：

	$\displaystyle\frac{\mathrm{d}\xi^{c}_{\text{new}}(t)}{\mathrm{d}t}$	$\displaystyle=h_{\xi}\,\eta^{c}_{\text{new}}(t)+f_{1\xi}^{l[c]}(t),$		(191)
	$\displaystyle\frac{\mathrm{d}\eta^{c}_{\text{new}}(t)}{\mathrm{d}t}$	$\displaystyle=h_{\eta}\,\xi^{c}_{\text{new}}(t)+f_{1\eta}^{l[c]}(t).$		(192)

Due to the minor eccentricity variations ( $\sim$ $10^{-6}$ ) induced by the $J_{2}$ perturbation, only the effects of the lunisolar perturbations are taken into account, resulting in expressions for $h_{\xi}$ and $h_{\eta}$ given by
由于扰动引起的微小偏心率变化（ $\sim$ $10^{-6}$ ），仅考虑阴阳扰动的影响，因此表达式 for $h_{\xi}$ 和 $h_{\eta}$ 由 $J_{2}$

\displaystyle h_{\xi}\simeq h_{\text{s}}^{\xi}+h_{\text{m}}^{\xi},\qquad h_{% \eta}\simeq h_{\text{s}}^{\eta}+h_{\text{m}}^{\eta},

(193)

with 跟

\displaystyle h_{\text{s}}^{\xi}=\frac{3}{8}c_{\text{s}}(1-5\cos 2i),\qquad h_% {\text{s}}^{\eta}=\frac{3}{2}c_{\text{s}},

(194)

and 和

$\displaystyle h_{\text{m}}^{\xi}=$	$\displaystyle~{}\frac{3\,c_{\text{m}}^{(2)}}{32}(1-5\cos 2i)(1+3\cos 2i_{3})-% \frac{45\,c_{\text{m}}^{(4)}}{65536}(3$
	$\displaystyle+12\cos 2i+49\cos 4i)(9+20\cos 2i_{3}+35\cos 4i_{3})$
	$\displaystyle-\frac{525\,c_{\text{m}}^{(6)}}{33554432}(30+71\cos 2i+114\cos 4i% +297\cos 6i)$
	$\displaystyle\times(50+105\cos 2i_{3}+126\cos 4i_{3}+231\cos 6i_{3}),$
$\displaystyle h_{\text{m}}^{\eta}=$	$\displaystyle~{}\frac{3\,c_{\text{m}}^{(2)}}{8}(1+3\cos 2i_{3})$
	$\displaystyle+\frac{45\,c_{\text{m}}^{(4)}}{8192}(3+5\cos 2i)(9+20\cos 2i_{3}+% 35\cos 4i_{3})$
	$\displaystyle+\frac{525\,c_{\text{m}}^{(6)}}{4194304}(15+28\cos 2i+21\cos 4i)$
	$\displaystyle\times(50+105\cos 2i_{3}+126\cos 4i_{3}+231\cos 6i_{3}).$	(195)

Moreover, the expressions for $f_{1\xi}^{l[c]}(t)$ and $f_{1\eta}^{l[c]}(t)$ are

	$\displaystyle f_{1\xi}^{l[c]}$	$\displaystyle=\sum_{q=-5}^{5}-f^{\xi}_{(q)}\sin(\Delta u_{3}+q\,\theta_{3}),$		(196)
	$\displaystyle f_{1\eta}^{l[c]}$	$\displaystyle=\sum_{\begin{subarray}{c}q=-5\\ (q\neq 0)\end{subarray}}^{5}f^{\eta}_{(q)}\cos(\Delta u_{3}+q\,\theta_{3}),$		(197)

with

	$\displaystyle f^{\xi}_{(q)}$	$\displaystyle:=c_{\text{m}}^{(3)}(4\,r_{3}^{\varepsilon}+u_{3}^{A})f^{\xi(3)}_% {[q]}(i)+c_{\text{m}}^{(5)}(6\,r_{3}^{\varepsilon}+u_{3}^{A})f^{\xi(5)}_{[q]}(% i),$		(198)
	$\displaystyle f^{\eta}_{(q)}$	$\displaystyle:=c_{\text{m}}^{(3)}(4\,r_{3}^{\varepsilon}+u_{3}^{A})f^{\eta(3)}% _{[q]}(i)+c_{\text{m}}^{(5)}(6\,r_{3}^{\varepsilon}+u_{3}^{A})f^{\eta(5)}_{[q]% }(i),$		(199)

and

	$\displaystyle\Delta u_{3}(t):=u_{3}(t)-\overline{M}_{3}(t)=\Delta n_{3}\,(t-t_% {0})+\Delta u_{3_{0}},$		(200)
	$\displaystyle\Delta n_{3}=n_{3}-n_{\!M_{3}},\qquad\Delta u_{3_{0}}=u_{3_{0}}-M% _{3_{0}},$		(201)

where $r_{3}^{\varepsilon}:=\frac{r_{3}^{A}}{\overline{r}_{3}}$ , and $f^{\eta}_{(0)}=0$ . Further, solving Eqs. (191) and (192) yields

	$\displaystyle\xi_{c}(t)=c_{1}\cosh(\sqrt{h_{\xi}}\sqrt{h_{\eta}}t)+c_{2}\frac{% \sqrt{h_{\xi}}}{\sqrt{h_{\eta}}}\sinh(\sqrt{h_{\xi}}\sqrt{h_{\eta}}t),$		(202)
	$\displaystyle\eta_{c}(t)=c_{2}\cosh(\sqrt{h_{\xi}}\sqrt{h_{\eta}}t)+c_{1}\frac% {\sqrt{h_{\eta}}}{\sqrt{h_{\xi}}}\sinh(\sqrt{h_{\xi}}\sqrt{h_{\eta}}t),$		(203)
	$\displaystyle\xi_{l[c]}(t)=\sum_{q=-5}^{5}\frac{f^{\xi}_{(q)}(\Delta n_{3}+q\,% n_{\theta_{3}})-h_{\xi}f^{\eta}_{(q)}}{(\Delta n_{3}+q\,n_{\theta_{3}})^{2}+h_% {\xi}h_{\eta}}\cos\Gamma,$		(204)
	$\displaystyle\eta_{l[c]}(t)=\sum_{q=-5}^{5}\frac{f^{\eta}_{(q)}(\Delta n_{3}+q% \,n_{\theta_{3}})+h_{\eta}f^{\xi}_{(q)}}{(\Delta n_{3}+q\,n_{\theta_{3}})^{2}+% h_{\xi}h_{\eta}}\sin\Gamma,$		(205)

where

\displaystyle c_{1}=\overline{\xi}_{\!0}-\xi_{l[c]}(t_{0}),\qquad c_{2}=

\displaystyle~{}\overline{\eta}_{\!0}-\eta_{l[c]}(t_{0}),

(206)

and $\Gamma:=\Delta u_{3}+q\,\theta_{3}$ . Equations (204) and (205) indicate that solar perturbation alone induces no special long-period variations in $\xi$ and $\eta$ , since the terms $f^{\xi}_{(q)}$ and $f^{\eta}_{(q)}$ are exclusively associated with lunar perturbation. Additionally, it is worth noting that introducing $\sigma_{l[c]}(t)$ in the reference solution (122), as well as in Eqs. (189) and (190), is crucial. Without this term, a significant increase in analytical solution errors would occur, leading to the disappearance of terms related to $\xi_{l[c]}(t_{0})$ and $\eta_{l[c]}(t_{0})$ in Eq. (206), as well as the terms associated with $h_{\xi}$ and $h_{\eta}$ in Eqs. (204) and (205).

Appendix D Explicit forms of inclination functions

The explicit forms of inclination functions within the lunar perturbation solution are shown below. These encompass $f^{\sigma(N)}_{[q]}(i)$ in Eqs. (164), (198), and (199) for the special long-period terms, $f^{\sigma(N)}_{[p,\,q]}(i)$ in Eq. (173) for the general long-period terms, and $f^{\sigma(N)}_{[\kappa,\,p,\,q]}(i)$ in Eqs. (180) and (181) for the short-period terms. For brevity, only the leading-order inclination functions with $N=2$ or $N=3$ are presented (cf. Table 7). Inclination functions for other orders can be derived using the methods outlined in Appendix B. Note that $i$ and $i_{3}$ in these functions represent mean values $\overline{i}$ and $\overline{i}_{3}$ , respectively.

Table 7: Relationship between the inclination functions and the Legendre polynomial degree

N

$\sigma$	$f^{\sigma(N)}_{[q]}(i)$	$f^{\sigma(N)}_{[p,\,q]}(i)$	$f^{\sigma(N)}_{[\kappa,\,p,\,q]}(i)$
$a$	$-$	$-$	$N=2,3,4,\cdots$
$\xi$ , $\eta$	$N=3,5,\cdots$	$N=3,5,\cdots$	$N=2,3,4,\cdots$
$i$ , $\Omega$ , $\lambda$	$N=2,4,\cdots$	$N=2,4,\cdots$	$N=2,3,4,\cdots$

D.1 Inclination functions for special long-period terms

The inclination functions $f^{\sigma(N)}_{[q]}(i)$ , associated with the special long-period terms $\sigma_{\text{m}}^{l[c]}(t)$ in Eqs. (204)-(205) and (161)-(163), are listed below in terms of the orbital elements $\xi$ , $\eta$ , $i$ , $\Omega$ , and $\lambda$ :

(1) Inclination functions with the form $f^{\xi(3)}_{[q]}(i)$ :

	$\displaystyle\displaystyle f^{\xi(3)}_{[\raisebox{0.0pt}{-}3]}=\frac{225}{3276% 8}\frac{(\cos i-\cos 3i)(3-4\cos 2i_{3}+\cos 4i_{3})}{1-\cos i_{3}},$
	$\displaystyle\displaystyle f^{\xi(3)}_{[0]}=\frac{45}{4096}(\sin i+5\sin 3i)(% \sin i_{3}+5\sin 3i_{3}),$
	$\displaystyle\displaystyle f^{\xi(3)}_{[\pm 1]}=-\frac{15}{16384}(\cos i+15% \cos 3i)(\pm 6-\cos i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\xi(3)}_{[\pm 1]}=}\pm 10\cos 2i_{3}-15% \cos 3i_{3}),$
	$\displaystyle\displaystyle f^{\xi(3)}_{[\pm 2]}=\frac{75}{8192}(\sin i-3\sin 3% i)(\sin i_{3}\pm 4\sin 2i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\xi(3)}_{[\pm 2]}=}-3\sin 3i_{3}),$
	$\displaystyle\displaystyle f^{\xi(3)}_{[3]}=-\frac{225}{16384}(\cos i-\cos 3i)% (2-\cos i_{3}-2\cos 2i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\xi(3)}_{[3]}=}+\cos 3i_{3}).$		(207)

(2) Inclination functions with the form $f^{\eta(3)}_{[q]}(i)$ :
（2）形式 $f^{\eta(3)}_{[q]}(i)$ 为的倾斜函数：

	$\displaystyle\displaystyle f^{\eta(3)}_{[\raisebox{0.0pt}{-}3]}=\frac{225}{163% 84}\frac{(1-\cos 2i)(3-4\cos 2i_{3}+\cos 4i_{3})}{1-\cos i_{3}},$
	$\displaystyle\displaystyle f^{\eta(3)}_{[\pm 1]}=\frac{15}{8192}(3+5\cos 2i)(6% \mp\cos i_{3}+10\cos 2i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\eta(3)}_{[\pm 1]}=}\mp 15\cos 3i_{3}),$		(208)
	$\displaystyle\displaystyle f^{\eta(3)}_{[\pm 2]}=\frac{75}{2048}\sin 2i(\pm% \sin i_{3}+4\sin 2i_{3}\mp 3\sin 3i_{3}),$
	$\displaystyle\displaystyle f^{\eta(3)}_{[3]}=\frac{225}{8192}(1-\cos 2i)(2-% \cos i_{3}-2\cos 2i_{3}+\cos 3i_{3}).$

(3) Inclination functions with the form $f^{i(2)}_{[q]}(i)$ :
（3）形式 $f^{i(2)}_{[q]}(i)$ 为的倾斜函数：

	$\displaystyle f^{i(2)}_{[1]}=-\frac{3}{8}\cos i\sin 2i_{3},$
	$\displaystyle f^{i(2)}_{[2]}=-\frac{3}{16}(1-\cos 2i_{3})\sin i.$		(209)

(4) Inclination functions with the form $f^{\Omega(2)}_{[q]}(i)$ :

	$\displaystyle f^{\Omega(2)}_{[1]}=\frac{3}{8}\cos 2i\csc i\sin 2i_{3},$
	$\displaystyle f^{\Omega(2)}_{[2]}=\frac{3}{16}\cos i(1-\cos 2i_{3}).$		(210)

(5) Inclination functions with the form $f^{\lambda(2)}_{[q]}(i)$ :

	$\displaystyle f^{\lambda(2)}_{[1]}=-\frac{3}{16}(3\cos i-\cos 3i)\csc i\sin 2i% _{3},$
	$\displaystyle f^{\lambda(2)}_{[2]}=-\frac{3}{32}(3-\cos 2i)(1-\cos 2i_{3}).$		(211)

D.2 Inclination functions for general long-period terms

The inclination functions $f^{\sigma(N)}_{[p,\,q]}(i)$ for the general long-period terms $\sigma_{\text{m}}^{l[l]}(t)$ in Eqs. (168)-(172) are listed in the order of $\xi$ , $\eta$ , $i$ , $\Omega$ , and $\lambda$ , as follows:

(1) Inclination functions with the form $f^{\xi(3)}_{[p,\,q]}(i)$ :

	$\displaystyle\displaystyle f^{\xi(3)}_{[1,\raisebox{0.0pt}{-}3]}=-\frac{225}{1% 6384}\frac{(\cos i-\cos 3i)(3-4\cos 2i_{3}+\cos 4i_{3})}{1-\cos i_{3}},$
	$\displaystyle\displaystyle f^{\xi(3)}_{[1,0]}=-\frac{45}{2048}(\sin i+5\sin 3i% )(\sin i_{3}+5\sin 3i_{3}),$
	$\displaystyle\displaystyle f^{\xi(3)}_{[1,\pm 1]}=\frac{15}{8192}(\cos i+15% \cos 3i)(\pm 6-\cos i_{3}\pm 10\cos 2i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\xi(3)}_{[1,\pm 1]}=}-15\cos 3i_{3}),$
	$\displaystyle\displaystyle f^{\xi(3)}_{[1,\pm 2]}=-\frac{75}{4096}(\sin i-3% \sin 3i)(\sin i_{3}\pm 4\sin 2i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\xi(3)}_{[1,\pm 2]}=}-3\sin 3i_{3}),$
	$\displaystyle\displaystyle f^{\xi(3)}_{[1,3]}=\frac{225}{8192}(\cos i-\cos 3i)% (2-\cos i_{3}-2\cos 2i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\xi(3)}_{[1,3]}=}+\cos 3i_{3}),$
	$\displaystyle\displaystyle f^{\xi(3)}_{[3,\raisebox{0.0pt}{-}2]}=\frac{75(3% \sin 3i-\sin i)(10\sin i_{3}-5\sin 3i_{3}+\sin 5i_{3})}{8192(3-4\cos i_{3}+% \cos 2i_{3})},$
	$\displaystyle\displaystyle f^{\xi(3)}_{[3,0]}=-\frac{75}{2048}(\sin i+5\sin 3i% )(3\sin i_{3}-\sin 3i_{3}),$
	$\displaystyle\displaystyle f^{\xi(3)}_{[3,\pm 1]}=\frac{75}{8192}(\cos i+15% \cos 3i)(\pm 2-\cos i_{3}\mp 2\cos 2i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\xi(3)}_{[3,\pm 1]}=}+\cos 3i_{3}),$
	$\displaystyle\displaystyle f^{\xi(3)}_{[3,2]}=\frac{75}{4096}(3\sin 3i-\sin i)% (5\sin i_{3}-4\sin 2i_{3}+\sin 3i_{3}),$
	$\displaystyle\displaystyle f^{\xi(3)}_{[3,\pm 3]}=\frac{75}{8192}(\cos i-\cos 3% i)(\pm 10-15\cos i_{3}\pm 6\cos 2i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\xi(3)}_{[3,\pm 3]}=}-\cos 3i_{3}).$		(212)

(2) Inclination functions with the form $f^{\eta(3)}_{[p,\,q]}(i)$ :

	$\displaystyle\displaystyle f^{\eta(3)}_{[1,\raisebox{0.0pt}{-}3]}=-\frac{225}{% 8192}\frac{(1-\cos 2i)(3-4\cos 2i_{3}+\cos 4i_{3})}{1-\cos i_{3}},$
	$\displaystyle\displaystyle f^{\eta(3)}_{[1,\pm 1]}=-\frac{15}{4096}(3+5\cos 2i% )(6\mp\cos i_{3}+10\cos 2i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\eta(3)}_{[1,\pm 1]}=}\mp 15\cos 3i_{3}),$
	$\displaystyle\displaystyle f^{\eta(3)}_{[1,\pm 2]}=-\frac{75}{1024}\sin 2i(\pm% \sin i_{3}+4\sin 2i_{3}\mp 3\sin 3i_{3}),$
	$\displaystyle\displaystyle f^{\eta(3)}_{[1,3]}=\frac{225}{4096}(\cos 2i-1)(2-% \cos i_{3}-2\cos 2i_{3}+\cos 3i_{3}),$
	$\displaystyle\displaystyle f^{\eta(3)}_{[3,\pm 1]}=-\frac{75}{4096}(3+5\cos 2i% )(2\mp\cos i_{3}-2\cos 2i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\eta(3)}_{[3,\pm 1]}=}\pm\cos 3i_{3}),$
	$\displaystyle\displaystyle f^{\eta(3)}_{[3,\pm 2]}=-\frac{75}{1024}\sin 2i(\pm 5% \sin i_{3}-4\sin 2i_{3}\pm\sin 3i_{3}),$
	$\displaystyle\displaystyle f^{\eta(3)}_{[3,\pm 3]}=-\frac{75}{4096}(1-\cos 2i)% (10\mp 15\cos i_{3}+6\cos 2i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\eta(3)}_{[3,\pm 3]}=}\mp\cos 3i_{3}).$		(213)

(3) Inclination functions with the form $f^{i(2)}_{[p,\,q]}(i)$ :

	$\displaystyle f^{i(2)}_{[2,\pm 1]}=-\frac{3}{16}\cos i(2\sin i_{3}\mp\sin 2i_{% 3}),$
	$\displaystyle f^{i(2)}_{[2,\pm 2]}=-\frac{3}{32}(\pm 3-4\cos i_{3}\pm\cos 2i_{% 3})\sin i.$		(214)

(4) Inclination functions with the form $f^{\Omega(2)}_{[p,\,q]}(i)$ :

	$\displaystyle f^{\Omega(2)}_{[2,0]}=-\frac{9}{16}\cos i(1-\cos 2i_{3}),$
	$\displaystyle f^{\Omega(2)}_{[2,\pm 1]}=\frac{3}{16}\cos 2i\csc i(\pm 2\sin i_% {3}-\sin 2i_{3}),$
	$\displaystyle f^{\Omega(2)}_{[2,\pm 2]}=\frac{3}{32}\cos i(3\mp 4\cos i_{3}+% \cos 2i_{3}).$		(215)

(5) Inclination functions with the form $f^{\lambda(2)}_{[p,\,q]}(i)$ :
（5）形式 $f^{\lambda(2)}_{[p,\,q]}(i)$ 为的倾斜函数：

	$\displaystyle f^{\lambda(2)}_{[2,0]}=\frac{3}{32}(1-3\cos 2i)(1-\cos 2i_{3}),$
	$\displaystyle f^{\lambda(2)}_{[2,\pm 1]}=-\frac{3}{32}(3\cos i-\cos 3i)\csc i(% \pm 2\sin i_{3}-\sin 2i_{3}),$
	$\displaystyle f^{\lambda(2)}_{[2,\pm 2]}=-\frac{3}{64}(3-\cos 2i)(3\mp 4\cos i% _{3}+\cos 2i_{3}).$		(216)

D.3 Inclination functions for short-period terms

The inclination functions $f^{\sigma(N)}_{[\kappa,\,p,\,q]}(i)$ for the short-period terms $\sigma_{\text{m}}^{s}(t)$ in Eqs. (174)-(179) are listed below in the order of $a$ , $\xi$ , $\eta$ , $i$ , $\Omega$ , and $\lambda$ :

(1) Inclination functions with the form $f^{a(2)}_{[\kappa,\,p,\,q]}(i)$ :

	$\displaystyle f^{a(2)}_{[2,\raisebox{0.0pt}{-}2,0]}=\frac{9}{32}(1-\cos 2i)(1-% \cos 2i_{3}),$
	$\displaystyle f^{a(2)}_{[2,\raisebox{0.0pt}{-}2,\pm 1]}=\frac{3}{16}(2\sin i% \pm\sin 2i)(2\sin i_{3}\pm\sin 2i_{3}),$
	$\displaystyle f^{a(2)}_{[2,\raisebox{0.0pt}{-}2,\pm 2]}=\frac{3}{64}(3\pm 4% \cos i+\cos 2i)(3\pm 4\cos i_{3}+\cos 2i_{3}),$
	$\displaystyle f^{a(2)}_{[2,0,0]}=\frac{3}{16}(1-\cos 2i)(1+3\cos 2i_{3}),$
	$\displaystyle f^{a(2)}_{[2,0,\pm 1]}=-\frac{3}{8}(\pm 2\sin i+\sin 2i)\sin 2i_% {3},$		(217)
	$\displaystyle f^{a(2)}_{[2,0,\pm 2]}=\frac{3}{32}(3\pm 4\cos i+\cos 2i)(1-\cos 2% i_{3}),$
	$\displaystyle f^{a(2)}_{[2,2,0]}=f^{a(2)}_{[2,\raisebox{0.0pt}{-}2,0]},$
	$\displaystyle f^{a(2)}_{[2,2,\pm 1]}=-\frac{3}{16}(2\sin i\pm\sin 2i)(2\sin i_% {3}\mp\sin 2i_{3}),$
	$\displaystyle f^{a(2)}_{[2,2,\pm 2]}=\frac{3}{64}(3\pm 4\cos i+\cos 2i)(3\mp 4% \cos i_{3}+\cos 2i_{3}).$

(2) Inclination functions with the form $f^{\xi(2)}_{[\kappa,\,p,\,q]}(i)$ :

	$\displaystyle\displaystyle f^{\xi(2)}_{[1,\raisebox{0.0pt}{-}2,0]}=\frac{3}{12% 8}(7-15\cos 2i)(1-\cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\xi(2)}_{[1,\raisebox{0.0pt}{-}2,\pm 1]}=\frac{3% }{64}(6\sin i\pm 5\sin 2i)(2\sin i_{3}\pm\sin 2i_{3}),$
	$\displaystyle\displaystyle f^{\xi(2)}_{[1,\raisebox{0.0pt}{-}2,\pm 2]}=\frac{3% }{256}(7\pm 12\cos i+5\cos 2i)(3\pm 4\cos i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\xi(2)}_{[1,\raisebox{0.0pt}{-}2,\pm 2]}% =}+\cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\xi(2)}_{[1,0,0]}=\frac{1}{64}(7-15\cos 2i)(1+3% \cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\xi(2)}_{[1,0,\pm 1]}=-\frac{3}{32}(\pm 6\sin i+% 5\sin 2i)\sin 2i_{3},$
	$\displaystyle\displaystyle f^{\xi(2)}_{[1,0,\pm 2]}=\frac{3}{128}(7\pm 12\cos i% +5\cos 2i)(1-\cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\xi(2)}_{[1,2,0]}=f^{\xi(2)}_{[1,\raisebox{0.0pt% }{-}2,0]},$
	$\displaystyle\displaystyle f^{\xi(2)}_{[1,2,\pm 1]}=-\frac{3}{64}(6\sin i\pm 5% \sin 2i)(2\sin i_{3}\mp\sin 2i_{3}),$
	$\displaystyle\displaystyle f^{\xi(2)}_{[1,2,\pm 2]}=\frac{3}{256}(7\pm 12\cos i% +5\cos 2i)(3\mp 4\cos i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\xi(2)}_{[1,2,\pm 2]}=}+\cos 2i_{3}),$		(218)
	$\displaystyle\displaystyle f^{\xi(2)}_{[3,\raisebox{0.0pt}{-}2,0]}=\frac{9}{12% 8}(1-\cos 2i)(1-\cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\xi(2)}_{[3,\raisebox{0.0pt}{-}2,\pm 1]}=\frac{3% }{64}(2\sin i\pm\sin 2i)(2\sin i_{3}\pm\sin 2i_{3}),$
	$\displaystyle\displaystyle f^{\xi(2)}_{[3,\raisebox{0.0pt}{-}2,\pm 2]}=\frac{3% }{256}(3\pm 4\cos i+\cos 2i)(3\pm 4\cos i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\xi(2)}_{[3,\raisebox{0.0pt}{-}2,\pm 2]}% =}+\cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\xi(2)}_{[3,0,0]}=\frac{3}{64}(1-\cos 2i)(1+3% \cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\xi(2)}_{[3,0,\pm 1]}=-\frac{3}{32}(\pm 2\sin i+% \sin 2i)\sin 2i_{3},$
	$\displaystyle\displaystyle f^{\xi(2)}_{[3,0,\pm 2]}=\frac{3}{128}(3\pm 4\cos i% +\cos 2i)(1-\cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\xi(2)}_{[3,2,0]}=f^{\xi(2)}_{[3,\raisebox{0.0pt% }{-}2,0]},$
	$\displaystyle\displaystyle f^{\xi(2)}_{[3,2,\pm 1]}=-\frac{3}{64}(2\sin i\pm% \sin 2i)(2\sin i_{3}\mp\sin 2i_{3}),$
	$\displaystyle\displaystyle f^{\xi(2)}_{[3,2,\pm 2]}=\frac{3}{256}(3\pm 4\cos i% +\cos 2i)(3\mp 4\cos i_{3}+\cos 2i_{3}).$

(3) Inclination functions with the form $f^{\eta(2)}_{[\kappa,\,p,\,q]}(i)$ :

	$\displaystyle\displaystyle f^{\eta(2)}_{[1,\raisebox{0.0pt}{-}2,0]}=-\frac{3}{% 128}(11-3\cos 2i)(1-\cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\eta(2)}_{[1,\raisebox{0.0pt}{-}2,\pm 1]}=-\frac% {3}{64}(6\sin i\pm\sin 2i)(2\sin i_{3}\pm\sin 2i_{3}),$
	$\displaystyle\displaystyle f^{\eta(2)}_{[1,\raisebox{0.0pt}{-}2,\pm 2]}=-\frac% {3}{256}(11\pm 12\cos i+\cos 2i)(3\pm 4\cos i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\eta(2)}_{[1,\raisebox{0.0pt}{-}2,\pm 2]% }=}+\cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\eta(2)}_{[1,0,0]}=-\frac{1}{64}(11-3\cos 2i)(1+% 3\cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\eta(2)}_{[1,0,\pm 1]}=\frac{3}{32}(\pm 6\sin i+% \sin 2i)\sin 2i_{3},$
	$\displaystyle\displaystyle f^{\eta(2)}_{[1,0,\pm 2]}=-\frac{3}{128}(11\pm 12% \cos i+\cos 2i)(1-\cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\eta(2)}_{[1,2,0]}=f^{\eta(2)}_{[1,\raisebox{0.0% pt}{-}2,0]},$
	$\displaystyle\displaystyle f^{\eta(2)}_{[1,2,\pm 1]}=\frac{3}{64}(6\sin i\pm% \sin 2i)(2\sin i_{3}\mp\sin 2i_{3}),$
	$\displaystyle\displaystyle f^{\eta(2)}_{[1,2,\pm 2]}=-\frac{3}{256}(11\pm 12% \cos i+\cos 2i)(3\mp 4\cos i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\eta(2)}_{[1,2,\pm 2]}=}+\cos 2i_{3}),$		(219)
	$\displaystyle\displaystyle f^{\eta(2)}_{[3,\raisebox{0.0pt}{-}2,0]}=\frac{9}{1% 28}(1-\cos 2i)(1-\cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\eta(2)}_{[3,\raisebox{0.0pt}{-}2,\pm 1]}=\frac{% 3}{64}(2\sin i\pm\sin 2i)(2\sin i_{3}\pm\sin 2i_{3}),$
	$\displaystyle\displaystyle f^{\eta(2)}_{[3,\raisebox{0.0pt}{-}2,\pm 2]}=\frac{% 3}{256}(3\pm 4\cos i+\cos 2i)(3\pm 4\cos i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\eta(2)}_{[3,\raisebox{0.0pt}{-}2,\pm 2]% }=}+\cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\eta(2)}_{[3,0,0]}=\frac{3}{64}(1-\cos 2i)(1+3% \cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\eta(2)}_{[3,0,\pm 1]}=-\frac{3}{32}(\pm 2\sin i% +\sin 2i)\sin 2i_{3},$
	$\displaystyle\displaystyle f^{\eta(2)}_{[3,0,\pm 2]}=\frac{3}{128}(3\pm 4\cos i% +\cos 2i)(1-\cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\eta(2)}_{[3,2,0]}=f^{\eta(2)}_{[3,\raisebox{0.0% pt}{-}2,0]},$
	$\displaystyle\displaystyle f^{\eta(2)}_{[3,2,\pm 1]}=-\frac{3}{64}(2\sin i\pm% \sin 2i)(2\sin i_{3}\mp\sin 2i_{3}),$
	$\displaystyle\displaystyle f^{\eta(2)}_{[3,2,\pm 2]}=\frac{3}{256}(3\pm 4\cos i% +\cos 2i)(3\mp 4\cos i_{3}+\cos 2i_{3}).$

(4) Inclination functions with the form $f^{i(2)}_{[\kappa,\,p,\,q]}(i)$ :

	$\displaystyle f^{i(2)}_{[2,\raisebox{0.0pt}{-}2,0]}=\frac{9}{64}(1-\cos 2i_{3}% )\sin 2i,$
	$\displaystyle f^{i(2)}_{[2,\raisebox{0.0pt}{-}2,\pm 1]}=\frac{3}{32}(\cos i\pm% \cos 2i)(2\sin i_{3}\pm\sin 2i_{3}),$
	$\displaystyle f^{i(2)}_{[2,\raisebox{0.0pt}{-}2,\pm 2]}=-\frac{3}{128}(\pm 3+4% \cos i_{3}\pm\cos 2i_{3})(2\sin i\pm\sin 2i),$
	$\displaystyle f^{i(2)}_{[2,0,0]}=\frac{3}{32}(1+3\cos 2i_{3})\sin 2i,$
	$\displaystyle f^{i(2)}_{[2,0,\pm 1]}=-\frac{3}{16}(\pm\cos i+\cos 2i)\sin 2i_{% 3},$		(220)
	$\displaystyle f^{i(2)}_{[2,0,\pm 2]}=-\frac{3}{64}(1-\cos 2i_{3})(\pm 2\sin i+% \sin 2i),$
	$\displaystyle f^{i(2)}_{[2,2,0]}=f^{i(2)}_{[2,\raisebox{0.0pt}{-}2,0]},$
	$\displaystyle f^{i(2)}_{[2,2,\pm 1]}=-\frac{3}{32}(\cos i\pm\cos 2i)(2\sin i_{% 3}\mp\sin 2i_{3}),$
	$\displaystyle f^{i(2)}_{[2,2,\pm 2]}=-\frac{3}{128}(\pm 3-4\cos i_{3}\pm\cos 2% i_{3})(2\sin i\pm\sin 2i).$

(5) Inclination functions with the form $f^{\Omega(2)}_{[\kappa,\,p,\,q]}(i)$ :
（5）形式 $f^{\Omega(2)}_{[\kappa,\,p,\,q]}(i)$ 为的倾斜函数：

	$\displaystyle f^{\Omega(2)}_{[2,\raisebox{0.0pt}{-}2,\raisebox{0.0pt}{-}1]}=% \frac{3}{32}\sec\frac{i}{2}\sin\frac{3i}{2}(2\sin i_{3}-\sin 2i_{3}),$
	$\displaystyle f^{\Omega(2)}_{[2,\raisebox{0.0pt}{-}2,0]}=\frac{9}{32}\cos i(1-% \cos 2i_{3}),$
	$\displaystyle f^{\Omega(2)}_{[2,\raisebox{0.0pt}{-}2,1]}=\frac{3}{32}\cos\frac% {3i}{2}\csc\frac{i}{2}(2\sin i_{3}+\sin 2i_{3}),$
	$\displaystyle f^{\Omega(2)}_{[2,\raisebox{0.0pt}{-}2,\pm 2]}=-\frac{3}{64}(1% \pm\cos i)(\pm 3+4\cos i_{3}\pm\cos 2i_{3}),$
	$\displaystyle f^{\Omega(2)}_{[2,0,\raisebox{0.0pt}{-}1]}=\frac{3}{16}\sec\frac% {i}{2}\sin\frac{3i}{2}\sin 2i_{3},$
	$\displaystyle f^{\Omega(2)}_{[2,0,0]}=\frac{3}{16}\cos i(1+3\cos 2i_{3}),$		(221)
	$\displaystyle f^{\Omega(2)}_{[2,0,1]}=-\frac{3}{16}\cos\frac{3i}{2}\csc\frac{i% }{2}\sin 2i_{3},$
	$\displaystyle f^{\Omega(2)}_{[2,0,\pm 2]}=-\frac{3}{32}(\pm 1+\cos i)(1-\cos 2% i_{3}),$
	$\displaystyle f^{\Omega(2)}_{[2,2,\raisebox{0.0pt}{-}1]}=-\frac{3}{32}\sec% \frac{i}{2}\sin\frac{3i}{2}(2\sin i_{3}+\sin 2i_{3}),$
	$\displaystyle f^{\Omega(2)}_{[2,2,0]}=f^{\Omega(2)}_{[2,\raisebox{0.0pt}{-}2,0% ]},$
	$\displaystyle f^{\Omega(2)}_{[2,2,1]}=-\frac{3}{32}\cos\frac{3i}{2}\csc\frac{i% }{2}(2\sin i_{3}-\sin 2i_{3}),$
	$\displaystyle f^{\Omega(2)}_{[2,2,\pm 2]}=-\frac{3}{64}(\pm 1+\cos i)(3\mp 4% \cos i_{3}+\cos 2i_{3}).$

(6) Inclination functions with the form $f^{\lambda(2)}_{[\kappa,\,p,\,q]}(i)$ :
（6）形式 $f^{\lambda(2)}_{[\kappa,\,p,\,q]}(i)$ 为的倾斜函数：

	$\displaystyle\displaystyle f^{\lambda(2)[a]}_{[2,\raisebox{0.0pt}{-}2,0]}=-% \frac{9}{64}(1-\cos 2i)(1-\cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\lambda(2)[a]}_{[2,\raisebox{0.0pt}{-}2,\pm 1]}=% \frac{-3}{32}(2\sin i\pm\sin 2i)(2\sin i_{3}\pm\sin 2i_{3}),$
	$\displaystyle\displaystyle f^{\lambda(2)[a]}_{[2,\raisebox{0.0pt}{-}2,\pm 2]}=% -\frac{3}{128}(3\pm 4\cos i+\cos 2i)(3\pm 4\cos i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\lambda(2)[a]}_{[2,\raisebox{0.0pt}{-}2,% \pm 2]}=}+\cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\lambda(2)[a]}_{[2,0,0]}=-\frac{3}{32}(1-\cos 2i% )(1+3\cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\lambda(2)[a]}_{[2,0,\pm 1]}=\frac{3}{16}(\pm 2% \sin i+\sin 2i)\sin 2i_{3},$
	$\displaystyle\displaystyle f^{\lambda(2)[a]}_{[2,0,\pm 2]}=-\frac{3}{64}(3\pm 4% \cos i+\cos 2i)(1-\cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\lambda(2)[a]}_{[2,2,0]}=f^{\lambda(2)[a]}_{[2,% \raisebox{0.0pt}{-}2,0]},$
	$\displaystyle\displaystyle f^{\lambda(2)[a]}_{[2,2,\pm 1]}=\frac{3}{32}(2\sin i% \pm\sin 2i)(2\sin i_{3}\mp\sin 2i_{3}),$
	$\displaystyle\displaystyle f^{\lambda(2)[a]}_{[2,2,\pm 2]}=-\frac{3}{128}(3\pm 4% \cos i+\cos 2i)(3\mp 4\cos i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\lambda(2)[a]}_{[2,2,\pm 2]}=}+\cos 2i_{% 3}),$
	$\displaystyle\displaystyle f^{\lambda(2)[\lambda]}_{[2,\raisebox{0.0pt}{-}2,% \raisebox{0.0pt}{-}1]}=-\frac{3}{64}\sec\frac{i}{2}(5\sin\frac{i}{2}+2\sin% \frac{3i}{2}-\sin\frac{5i}{2})(2\sin i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\lambda(2)[\lambda]}_{[2,\raisebox{0.0pt% }{-}2,\raisebox{0.0pt}{-}1]}=}-\sin 2i_{3}),$		(222)
	$\displaystyle\displaystyle f^{\lambda(2)[\lambda]}_{[2,\raisebox{0.0pt}{-}2,0]% }=-\frac{9}{64}(3-\cos 2i)(1-\cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\lambda(2)[\lambda]}_{[2,\raisebox{0.0pt}{-}2,1]% }=-\frac{3}{64}\csc\frac{i}{2}(5\cos\frac{i}{2}-2\cos\frac{3i}{2}-\cos\frac{5i% }{2})(2\sin i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\lambda(2)[\lambda]}_{[2,\raisebox{0.0pt% }{-}2,1]}=}+\sin 2i_{3}),$
	$\displaystyle\displaystyle f^{\lambda(2)[\lambda]}_{[2,\raisebox{0.0pt}{-}2,% \pm 2]}=-\frac{3}{128}(5\pm 6\cos i+\cos 2i)(3\pm 4\cos i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\lambda(2)[\lambda]}_{[2,\raisebox{0.0pt% }{-}2,\pm 2]}=}+\cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\lambda(2)[\lambda]}_{[2,0,\raisebox{0.0pt}{-}1]% }=-\frac{3}{32}\sec\frac{i}{2}(5\sin\frac{i}{2}+2\sin\frac{3i}{2}-\sin\frac{5i% }{2})\sin 2i_{3},$
	$\displaystyle\displaystyle f^{\lambda(2)[\lambda]}_{[2,0,0]}=-\frac{3}{32}(3-% \cos 2i)(1+3\cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\lambda(2)[\lambda]}_{[2,0,1]}=\frac{3}{32}\csc% \frac{i}{2}(5\cos\frac{i}{2}-2\cos\frac{3i}{2}-\cos\frac{5i}{2})\sin 2i_{3},$
	$\displaystyle\displaystyle f^{\lambda(2)[\lambda]}_{[2,0,\pm 2]}=-\frac{3}{64}% (5\pm 6\cos i+\cos 2i)(1-\cos 2i_{3}),$
	$\displaystyle\displaystyle f^{\lambda(2)[\lambda]}_{[2,2,\raisebox{0.0pt}{-}1]% }=\frac{3}{64}\sec\frac{i}{2}(5\sin\frac{i}{2}+2\sin\frac{3i}{2}-\sin\frac{5i}% {2})(2\sin i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\lambda(2)[\lambda]}_{[2,2,\raisebox{0.0% pt}{-}1]}=}+\sin 2i_{3}),$
	$\displaystyle\displaystyle f^{\lambda(2)[\lambda]}_{[2,2,0]}=f^{\lambda(2)[% \lambda]}_{[2,\raisebox{0.0pt}{-}2,0]},$
	$\displaystyle\displaystyle f^{\lambda(2)[\lambda]}_{[2,2,1]}=\frac{3}{64}\csc% \frac{i}{2}(5\cos\frac{i}{2}-2\cos\frac{3i}{2}-\cos\frac{5i}{2})(2\sin i_{3}$
	$\displaystyle\displaystyle\phantom{f^{\lambda(2)[\lambda]}_{[2,2,1]}=}-\sin 2i% _{3}),$
	$\displaystyle\displaystyle f^{\lambda(2)[\lambda]}_{[2,2,\pm 2]}=\frac{-3}{128% }(5\pm 6\cos i+\cos 2i)(3\mp 4\cos i_{3}+\cos 2i_{3}).$

References 引用

B. P. Abbott et al. [2016]
BP Abbott 等人 [2016] B. P. Abbott et al., Phys. Rev. Lett. 116, 061102 (2016).
BP Abbott 等人，Phys. Rev. Lett. 116， 061102 （2016）。
Ni [2016] 倪 [2016] W.-T. Ni, Int. J. Mod. Phys. D 25, 1630001 (2016).
W.-T. Ni，国际 J. Mod.Phys. D 25， 1630001 （2016）。
[3] NASA Gravitational-Wave Mission Concept Study Final Report (2012), https://www.cosmos.esa.int/documents/427239/442202/GW_Study_Rev3_Aug2012-Final-1.pdf.
NASA 引力波任务概念研究最终报告（2012 年），https://www.cosmos.esa.int/documents/427239/442202/GW_Study_Rev3_Aug2012-Final-1.pdf。
[4] P. Amaro-Seoane et al., arXiv:1702.00786.
P. Amaro-Seoane 等人，arXiv：1702.00786。
[5] M. Colpi et al., arXiv:2402.07571.
M. Colpi 等人，arXiv：2402.07571。
Seto et al. [2001]
Seto et al. [2001] N. Seto, S. Kawamura, and T. Nakamura, Phys. Rev. Lett. 87, 221103 (2001).
N. Seto、S. Kawamura 和 T. Nakamura，Phys. Rev. Lett. 87， 221103 （2001）。
J. Luo et al. [2016]
J. Luo 等 [2016] J. Luo et al., Class. Quantum Grav. 33, 035010 (2016).
J. Luo et al.，量子重力类.33， 035010 （2016）.
Hu and Wu [2017] 胡和吴 [2017] W.-R. Hu and Y.-L. Wu, Natl Sci. Rev. 4, 685 (2017).
W.-R. 胡和 Y.-L.Wu， Natl Sci. Rev. 4， 685 （2017）。
Hu et al. [2017]
胡 et al. [2017] Y.-M. Hu, J. Mei, and J. Luo, Natl Sci. Rev. 4, 683 (2017).
Y.-M. 胡， J. Mei，和 J. 罗，国家科学修订版 4， 683 （2017）.
J. Mei et al. [2020]
J. Mei 等人 [2020] J. Mei et al., Prog. Theor. Exp. Phys. 2020, 05A107 (2020).
J. Mei 等人，Prog。理论。物理实验 2020， 05A107 （2020）.
Gong et al. [2021]
Gong et al. [2021] Y. Gong, J. Luo, and B. Wang, Nature Astron. 5, 881 (2021).
Y. Gong、J. Luo 和 B. Wang，Nature Astron.5， 881 （2021）。
Tinto and Dhurandhar [2014]
Tinto 和 Dhurandhar [2014] M. Tinto and S. V. Dhurandhar, Living Reviews in Relativity 17, 1 (2014).
M. Tinto 和 S. V. Dhurandhar，《相对论生活评论》第 17 卷第 1 期（2014 年）。
Zhou et al. [2021]
周等 [2021] M.-Y. Zhou, X.-C. Hu, B. Ye, S. Hu, D.-D. Zhu, X. Zhang, W. Su, and Y. Wang, Phys. Rev. D 103, 103026 (2021).
周 M.-Y. X.-C.胡，叶斌， S. 胡， D.-D.Zhu， X. Zhang， W. Su，和 Y. Wang， Phys. Rev. D 103， 103026 （2021）。
Zheng et al. [2023]
Zheng et al. [2023] L. Zheng, S. Yang, and X. Zhang, Phys. Rev. D 108, 022001 (2023).
L. Zheng、S. Yang 和 X. Zhang，Phys. Rev. D 108， 022001 （2023）。
Folkner et al. [1997]
Folkner等人 [1997] W. M. Folkner, F. Hechler, T. H. Sweetser, M. A. Vincent, and P. L. Bender, Class. Quantum Grav. 14, 1405 (1997).
W. M. Folkner、F. Hechler、TH Sweetser、MA Vincent 和 PL Bender，量子重力。14， 1405 （1997）.
Roscoe et al. [2013]
Roscoe 等人 [2013] C. W. T. Roscoe, S. R. Vadali, and K. T. Alfriend, J. Astronaut. Sci. 60, 408 (2013).
CWT Roscoe、SR Vadali 和 KT Alfriend，J. 宇航员。科学 60， 408 （2013）.
Wu et al. [2019]
Wu et al. [2019] B. Wu, C.-G. Huang, and C.-F. Qiao, Phys. Rev. D 100, 122001 (2019).
B. Wu， C.-G.Huang 和 C.-F.乔，物理学修订版 D 100,122001 （2019）。
De Marchi et al. [2012]
De Marchi 等人 [2012] F. De Marchi, G. Pucacco, and M. Bassan, Class. Quantum Grav. 29, 035009 (2012).
F. De Marchi、G. Pucacco 和 M. Bassan，量子引力类。29， 035009 （2012）。
Qiao et al. [2023a]
Qiao et al. [2023a] D. Qiao, F. Jia, X. Li, and X. Zhou, Space Sci. Technol. 3, 0015 (2023a).
D. Qiao， F. Jia， X. Li，和 X. 周，空间科学技术3， 0015 （2023a）.
Dhurandhar et al. [2005]
Dhurandhar 等人 [2005] S. V. Dhurandhar, K. R. Nayak, S. Koshti, and J.-Y. Vinet, Class. Quantum Grav. 22, 481 (2005).
SV Dhurandhar、KR Nayak、S. Koshti 和 JYVinet，类 Quantum Grav.22， 481 （2005 年）。
Hu et al. [2018]
胡 et al. [2018] X.-C. Hu, X.-H. Li, Y. Wang, W.-F. Feng, M.-Y. Zhou, Y.-M. Hu, S.-C. Hu, J.-W. Mei, and C.-G. Shao, Class. Quantum Grav. 35, 095008 (2018).
X.-C. 胡， X.-H.李彦 Y. Wang， W.-F.冯 M.-Y.周 Y.-M.胡， S.-C.胡，J.-W.Mei 和 C.-G.Shao，类 Quantum Grav.35， 095008 （2018）。
Qiao et al. [2023b]
Qiao et al. [2023b] D. Qiao, X. Zhou, and X. Li, Astrodynamics 7, 271 (2023b).
D. Qiao， X. 周，和 X. Li，天体动力学 7， 271 （2023b）.
Jiao and Zhang [2023] 焦和张 [2023] L. Jiao and X. Zhang, Phys. Rev. D 107, 102004 (2023).
L. Jiao 和 X. Zhang，Phys. Rev. D 107， 102004 （2023）。
Ye et al. [2019]
Ye et al. [2019] B. Ye, X. Zhang, M. Zhou, Y. Wang, H. Yuan, D. Gu, Y. Ding, J. Zhang, J. Mei, and J. Luo, Int. J. Mod. Phys. D 28, 1950121 (2019).
叶斌，张晓明，周明，王妍，袁浩，顾大军，丁彦，张俊杰，梅俊杰，罗俊杰，国际J. Mod.Phys. D 28， 1950121 （2019 年）。
Tan et al. [2020]
Tan 等人 [2020] Z. Tan, B. Ye, and X. Zhang, Int. J. Mod. Phys. D 29, 2050056 (2020).
Z. Tan， B. Ye，和 X. Zhang， Int. J. Mod.Phys. D 29， 2050056 （2020）。
Martens and Joffre [2021]
Martens 和 Joffre [2021] W. Martens and E. Joffre, J. Astronaut. Sci. 68, 402 (2021).
W. Martens 和 E. Joffre，J. 宇航员。科学 68， 402 （2021）。
Nayak et al. [2006]
Nayak 等人 [2006] K. R. Nayak, S. Koshti, S. V. Dhurandhar, and J.-Y. Vinet, Class. Quantum Grav. 23, 1763 (2006).
KR Nayak、S. Koshti、SV Dhurandhar 和 JYVinet，类 Quantum Grav.23， 1763 （2006）。
Yi et al. [2008]
Yi et al. [2008] Z. Yi, G. Li, G. Heinzel, A. Ruediger, O. Jennrich, L. Wang, Y. Xia, F. Zeng, and H. Zhao, Int. J. Mod. Phys. D 17, 1005 (2008).
Z. Yi， G. Li， G. Heinzel， A. Ruediger， O. Jennrich， L. Wang， Y. Xia， F. Zeng， and H. Zhao， Int. J. Mod.Phys. D 17， 1005 （2008）。
Cerdonio et al. [2010]
Cerdonio 等人 [2010] M. Cerdonio, F. De Marchi, R. De Pietri, P. Jetzer, F. Marzari, G. Mazzolo, A. Ortolan, and M. Sereno, Class. Quantum Grav. 27, 165007 (2010).
M. Cerdonio， F. De Marchi， R. De Pietri， P. Jetzer， F. Marzari， G. Mazzolo， A. Ortolan，和 M. Sereno， Class.量子引力。27， 165007 （2010）.
Pucacco et al. [2010]
Pucacco 等人 [2010] G. Pucacco, M. Bassan, and M. Visco, Class. Quantum Grav. 27, 235001 (2010).
G. Pucacco、M. Bassan 和 M. Visco，类 Quantum Grav。27， 235001 （2010）.
Wu et al. [2020]
Wu et al. [2020] B. Wu, C.-G. Huang, and C.-F. Qiao, Phys. Rev. D 101, 064049 (2020).
B. Wu， C.-G.Huang 和 C.-F.乔，物理修订版 D 101,064049 （2020）。
Nie et al. [2019]
聂等人 [2019] T. Nie, P. Gurfil, and S. Zhang, Celest. Mech. Dyn. Astron. 131, 1 (2019).
T. Nie、P. Gurfil 和 S. Zhang，Celest。机械 Dyn.Astron. 131， 1 （2019）.
[33] M. Capderou, Handbook of Satellite Orbits: From Kepler to GPS (Springer Science & Business, 2014).
M. Capderou，《卫星轨道手册：从开普勒到GPS》（Springer Science & Business，2014年）。
Kozai [1973] 香菜 [1973] Y. Kozai, SAO Special Report 349 (1973).
Y. Kozai，SAO 特别报告 349 （1973）。
[35] L. Liu, Orbit Theory of Spacecraft (in Chinese) (National Defence Industry Press, Beijing, 2000).
刘丽，《航天器轨道理论》（国防工业出版社，北京，2000 年）。
Xu et al. [2011]
Xu et al. [2011] G. Xu, T. Xu, T.-K. Yeh, and W. Chen, Mon. Not. R. Astron. Soc. 410, 645 (2011).
G. Xu， T. Xu， T.-K. Yeh 和 W. Chen， Mon. Not. R. Astron. Soc. 410， 645 （2011）。
Ye et al. [2021]
Ye et al. [2021] B. Ye, X. Zhang, Y. Ding, and Y. Meng, Phys. Rev. D 103, 042007 (2021).
B. Ye、X. Zhang、Y. Ding 和 Y. Meng，Phys. Rev. D 103， 042007 （2021）。
Zhang et al. [2021]
Zhang et al. [2021] X. Zhang, C. Luo, L. Jiao, B. Ye, H. Yuan, L. Cai, D. Gu, J. Mei, and J. Luo, Phys. Rev. D 103, 062001 (2021).
Luo and Zhang [2022] 罗和张 [2022] C. Luo and X. Zhang, Phys. Rev. D 105, 102007 (2022).
C. Luo 和 X. Zhang，Phys. Rev. D 105， 102007 （2022）。
Yao et al. [2023]
Yao et al. [2023] Z. Yao, J. Lv, H. Chen, X. Zhang, and Y. Meng, Advances in Space Research 72, 997 (2023).
Z. Yao、J. Lv、H. Chen、X. Zhang 和 Y. Meng，空间研究进展 72， 997 （2023）。
Jia et al. [2023]
Jia et al. [2023] F. Jia, X. Li, D. Qiao, and X. Zhou, Acta Astronautica 202, 522 (2023).
F. Jia， X. Li， D. Qiao，和 X. 周，宇航学报 202， 522 （2023）.
Zhang et al. [2023]
Zhang et al. [2023] Z. Zhang, D. Gu, X. Zhang, B. Ye, Y. Ding, C. Wei, K. Shao, and J. Zhu, Class. Quantum Grav. 27, 235001 (2023).
Z. Zhang， D. Gu， X. Zhang， B. Ye， Y. Ding， C. Wei， K. Shao，和 J. Zhu，类量子引力。27， 235001 （2023）.
[43] G. Xu and J. Xu, Orbits: 2nd Order Singularity-free Solutions (Springer Science & Business Media, 2013).
G. Xu 和 J. Xu，轨道：二阶无奇点解决方案（Springer Science & Business Media， 2013）.
Folkner et al. [2014]
Folkner 等人 [2014] W. M. Folkner, J. G. Williams, D. H. Boggs, R. S. Park, and P. Kuchynka, IPN Progress Report 196, 42 (2014).
WM Folkner、JG Williams、DH Boggs、RS Park 和 P. Kuchynka，IPN 进展报告 196,42 （2014）。
[45] O. Montenbruck and E. Gill, Satellite Orbits: Models, Methods and Applications (Springer, Berlin, 2001).
O. Montenbruck 和 E. Gill，《卫星轨道：模型、方法和应用》（Springer，柏林，2001 年）。
[46] https://ssd.jpl.nasa.gov/horizons.cgi.
Prado [2003] 普拉多 [2003] A. F. B. A. Prado, J. Guid. Control. Dyn. 26, 33 (2003).
AF BAA 普拉多，J. Guid。控制。男人。26， 33 （2003）。
Zhou et al. [2022]
周等 [2022] X. Zhou, F. Jia, and X. Li, Aerospace 9, 519 (2022).
X. 周， F. Jia， and X. Li，航天 9， 519 （2022）.
[49] F. Jia, X. Li, X. Zhou, and D. Qiao, IEEE International Conference on Unmanned Systems (ICUS), Guangzhou, China, 2022, pp. 1233-1237.
F. Jia、X. Li、X. 周和 D. Qiao，IEEE 无人系统国际会议（ICUS），中国广州，2022 年，第 1233-1237 页。
An et al. [2022]
安等人 [2022] Z. An, K. Shao, D. Gu, J. Zhu, M. Li, L. Tong, and C. Wei, Class. Quantum Grav. 39, 245016 (2022).
Z. An， K. Shao， D. Gu， J. Zhu， M. Li， L. Tong， and C. Wei，量子重力类。39， 245016 （2022）.
Kozai [1962] 小杂 [1962] Y. Kozai, Astron. J. 67, 591 (1962).
Y. Kozai， Astron. J. 67， 591 （1962）。
Lidov [1962] 利多夫 [1962] M. L. Lidov, Planet. Space Sci. 9, 719 (1962).
M. L. Lidov，行星。空间科学 9， 719 （1962）。
Smith [1962] 史密斯 [1962] D. E. Smith, Planet. Space Sci. 9, 659 (1962).
D. E. 史密斯，行星。空间科学 9， 659 （1962）。
[54] https://gmat.atlassian.net/wiki/spaces/GW/overview.

Effects of lunisolar perturbations on TianQin constellation: An analytical model阴阳扰动对天琴星座的影响：一个解析模型

Abstract 抽象

I Introduction I 引言

II Fundamentals of stable TianQin constellation二、稳定天琴座的基本原理

II.1 Keplerian orbits of satellitesII.1 卫星的开普勒轨道

II.2 Orbit design of the TianQin constellation二.2 天琴座的轨道设计

III Effects of lunisolar perturbations on TianQin constellation三、阴阳扰动对天秦星座的影响

III.1 Dynamic model三.1 动态模型

III.2 Motion of the Sun and Moon三.2 太阳和月亮的运动

III.3 Lunisolar perturbations on the TianQin satellites三.3 天琴卫星上的阴阳扰动

III.4 Perturbed motion of the TianQin constellationIII.4 天琴星座的扰动运动

IV Concluding Remarks四、结束语

Acknowledgements. 确认。

Appendix A Table of Symbols附录 A 符号表

Appendix B Model derivation and verification附录 B 模型推导和验证

B.1 Derivation and motivation of Eq. (54)B.1 方程 （54） 的推导和动机

B.2 Derivation of Eqs. (73)-(75)乙.2 方程的推导。（73）-（75）

B.3 Verification of Eqs. (73)-(75)乙.3 方程验证。（73）-（75）

B.4 Model verification for the three indicators乙.4 3个指标的模型验证

Appendix C Explicit expressions for terms in Eq. (75)附录 C 方程 （75） 中项的显式表达式

C.1 Solar perturbation solution σs⁢(t)subscript𝜎s𝑡\sigma_{\text{s}}(t)italic_σ start_POSTSUBSCRIPT s end_POSTSUBSCRIPT ( italic_t )C.1 太阳能扰动解决方案 σs⁢(t)subscript𝜎s𝑡\sigma_{\text{s}}(t)italic_σ start_POSTSUBSCRIPT s end_POSTSUBSCRIPT ( italic_t )

C.2 Lunar perturbation solution σm⁢(t)subscript𝜎m𝑡\sigma_{\text{m}}(t)italic_σ start_POSTSUBSCRIPT m end_POSTSUBSCRIPT ( italic_t )C.2 月球扰动解 σm⁢(t)subscript𝜎m𝑡\sigma_{\text{m}}(t)italic_σ start_POSTSUBSCRIPT m end_POSTSUBSCRIPT ( italic_t )

Appendix D Explicit forms of inclination functions

D.1 Inclination functions for special long-period terms

D.2 Inclination functions for general long-period terms

D.3 Inclination functions for short-period terms

References 引用

Effects of lunisolar perturbations on TianQin constellation: An analytical model
阴阳扰动对天琴星座的影响：一个解析模型

II Fundamentals of stable TianQin constellation
二、稳定天琴座的基本原理

II.1 Keplerian orbits of satellites
II.1 卫星的开普勒轨道

II.2 Orbit design of the TianQin constellation
二.2 天琴座的轨道设计

III Effects of lunisolar perturbations on TianQin constellation
三、阴阳扰动对天秦星座的影响

III.1 Dynamic model
三.1 动态模型

III.2 Motion of the Sun and Moon
三.2 太阳和月亮的运动

III.3 Lunisolar perturbations on the TianQin satellites
三.3 天琴卫星上的阴阳扰动

III.4 Perturbed motion of the TianQin constellation
III.4 天琴星座的扰动运动

IV Concluding Remarks
四、结束语

Appendix A Table of Symbols
附录 A 符号表

Appendix B Model derivation and verification
附录 B 模型推导和验证

B.1 Derivation and motivation of Eq. (54)
B.1 方程（54）的推导和动机

B.2 Derivation of Eqs. (73)-(75)
乙.2 方程的推导。（73）-（75）

B.3 Verification of Eqs. (73)-(75)
乙.3 方程验证。（73）-（75）

B.4 Model verification for the three indicators
乙.4 3个指标的模型验证

Appendix C Explicit expressions for terms in Eq. (75)
附录 C 方程（75）中项的显式表达式

C.1 Solar perturbation solution $\sigma_{\text{s}}(t)$
C.1 太阳能扰动解决方案 $\sigma_{\text{s}}(t)$

C.2 Lunar perturbation solution $\sigma_{\text{m}}(t)$
C.2 月球扰动解 $\sigma_{\text{m}}(t)$