Optimizing orbits for TianQin
Bo-Bing Ye, Xuefeng Zhang, Ming-Yue Zhou, Yan Wang, Hui-Min Yuan, Defeng Gu, Yanwei Ding, Jinxiu Zhang, Jianwei Mei, Jun Luo
OOptimizing orbits for TianQin
Bo-Bing Ye , , Xuefeng Zhang , , Ming-Yue Zhou , Yan Wang ,Hui-Min Yuan , , Defeng Gu , , Yanwei Ding , , Jinxiu Zhang , , JianweiMei , , Jun Luo , TianQin Research Center for Gravitational Physics, Sun Yat-sen University, Zhuhai 519082,China School of Physics and Astronomy, Sun Yat-sen University, Zhuhai 519082, China MOE Key Laboratory of Fundamental Physical Quantities Measurements,Hubei Key Laboratory of Gravitation and Quantum Physics,School of Physics, Huazhong University of Science and Technology, Wuhan 430074, ChinaE-mail: [email protected]
Abstract.
TianQin is a geocentric space-based gravitational-wave observatory missionconsisting of three drag-free controlled satellites in an equilateral triangle with an orbital radius of10 km. The constellation faces the white-dwarf binary RX J0806.3 + Keywords : space gravitational-wave detection, TianQin, geocentric orbit, orbit design,orbit optimization
1. Introduction
TianQin is a proposed space-borne science mission to detect gravitational waves (GW)in the mHz frequency band [1]. The mission concept relies on a constellation of threeidentical drag-free controlled spacecraft in high Earth orbits at an altitude about 10 km (figure 1). The constellation forms a nearly equilateral triangle, and the nominalorbital plane stands almost perpendicular to the ecliptic, facing the white-dwarf binaryRX J0806.3 + a r X i v : . [ g r- q c ] D ec ptimizing orbits for TianQin Figure 1.
An illustration of the TianQin constellation consisting of three spacecraft SC1, SC2,and SC3 (figure reproduced from [1]). The direction to J0806 is shown. careful orbit analysis and optimal design are necessary to alleviate demands on on-board instruments.At the geocentric distance of 10 km, the gravitational perturbations to the TianQinspacecraft are primarily caused by, in descending order of magnitudes, the Moon, theSun, and the Earth’s J oblateness. Their relative magnitudes are, respectively, of4 × − , 2 × − , and 6 × − , compared with the central force from the Earth [4].Hence, the Moon constitutes the largest factor affecting the formation stability.In regard to orbit stability and optimization for other space GW detection missions,the well-known heliocentric LISA design [5–7] has been extensively studied with analyticand numerical methods [8–17]. Particularly, the cost-function method based on carefullychosen performance measures has proved effective in numerical optimization [12–16].For the ASTROD-GW mission [18], optimized orbits near the Sun-Earth Lagrangepoints (L3, L4, and L5) have been acquired through tuning the average periods andeccentricities of the orbits [19–22]. Our optimization scheme benefits from these studies.An example of simulated TianQin orbits has been given in [1] (also [23]), however,without providing much details and further discussion. In this study, we intend to fill inthe blanks through an independent verification, and extend the current understandingin TianQin’s geocentric orbit design. The paper is organized as follows. In section 2, weintroduce the nominal circular orbits of TianQin and preliminary stability requirementson constellation geometry (estimates on eccentricities given in Appendix A). In section3, we describe the simulation tool and initial setup for orbit integration. In section 4,the methods and three steps of optimization are presented (detailed derivation deferredto Appendix B and Appendix C). The optimized TianQin orbits are shown in section 5,and in section 6, we analyze observed secular pointing drift from the reference source.In section 7, we present the optimization results for other orbital pointings. The paperconcludes in section 8. ptimizing orbits for TianQin
2. TianQin mission requirements
The nominal orbits of the TianQin constellation can be readily given in terms of Kepler’scircular orbits. For a design baseline, the three spacecraft have the same orbit radius andform an equilateral formation revolving around the Earth. They fly in the same orbitalplane oriented constantly towards the designated reference source J0806. Using theJ2000-based Earth-centered ecliptic coordinate system (EarthMJ2000Ec), we prescribethe nominal orbits as follows ‡ [1]: a = 10 km ,e = 0 ,i = 94 .
704 035 ◦ , Ω = 210 .
443 557 ◦ ,ω = 0 ◦ ,ν k = ν + ( k − · ◦ , k = 1 , , . (1)The orbital elements listed above include the semi-major axis a , eccentricity e ,inclination i , longitude of ascending node Ω, argument of periapsis ω , and true anomaly ν k ( ν for the spacecraft SC1, etc.). The geocentric formation travels in prograde motionwith a period of ∼ . i , the nominal orbital plane standsalmost upright to the ecliptic plane, which helps reduce direct sunlight into the opticalassemblies. Moreover, the normal of the triangle is aiming at J0806. The fixed detectorpointing of the nominal constellation differs greatly from the LISA orbit design whichfeatures a yearly pointing variation.Due to various perturbing forces in space, real-world orbits deviate from the simpleKeplerian approximation. The arm lengths, relative line-of-sight velocities (range rates),and breathing angles (between arms) of the spacecraft formation, as well as the orbitalplanes, undergo continuous changes. To accommodate instrumentation and scienceoperations such as Doppler measurement, telescope steering, and drag-free controlon board the spacecraft [6], the orbits should be designed to meet certain stabilityrequirements. For this study alone, we have assumed the following permissible variationsof constellation geometry in table 1 [1]. Note that no requirement is imposed on pointingstability at J0806 since we expect it to be non-critical to the mission (in fact, largerpointing variations benefit sky localization of GW sources).
3. Orbit propagation
The orbit propagation is implemented by the NASA General Mission Analysis Tool(GMAT) [24], which is an open-source, flight qualified software extensively used forspace mission design [25]. The force models we have adopted include a 10 ×
10 spherical-harmonic model of the Earth’s gravity field (JGM-3 [26]), the point-mass gravity field ‡ In the geocentric equatorial coordinate system, one has i = 74 .
541 611 ◦ and Ω = 211 .
596 667 ◦ . ptimizing orbits for TianQin Table 1. L ij ± × ( √ × ) kmRelative velocity v ij ±
10 m / s for 5 years ± / s for the first 2 yearsBreathing angle α i ± . ◦ for 5 years ± . ◦ for the first 2 years from the Moon, Sun and solar system planets (the ephemeris DE421 [27]), and thefirst-order relativistic correction § . As the spacecraft are drag-free controlled, we onlyconsider, for a design baseline, purely gravitational trajectories and assume no orbitcorrection maneuvers performed during the formation flight. A ninth-order Runge-Kutta integrator with eighth-order error control (RungeKutta89) is used. GMAT alsoprovides an optimization solver (fmincon) through an interface with MATLAB.Under planetary perturbation, orbit propagation from the initial elements providedby the nominal orbits, unless under fortuitous circumstances, fails to satisfy the stabilityrequirements. Generally in these cases, a long-term linear drift in arm lengths andbreathing angles can be observed, indicating one spacecraft chasing another withinthe constellation (see, e.g., [19, 28]). Therefore, one needs to adjust the initial orbitalelements (positions and velocities), as free variables, to stabilize the ensuing relativeorbital motion. The nominal orbits provide a suitable initial guess for such anoptimization procedure. Hence we have assumed the initial elements and epoch (22May, 2034 12:00:00 UTC, for testing optimization only) in table 2. Table 2.
The initial elements from the nominal orbits of the TianQin constellation in theJ2000-based Earth-centered ecliptic coordinate system (EarthMJ2000Ec) at the epoch 22 May,2034 12:00:00 UTC. a (km) e i ( ◦ ) Ω ( ◦ ) ω ( ◦ ) ν ( ◦ )SC1, 2, 3 10
4. Optimization method
The optimization starts with orbit propagation from a set of initial elements σ ≡ ( a , e , i , Ω , ω , ν ) (or equivalently, initial positions and velocities) at an epoch t ,which is usually derived from nominal orbits. Based on the resulting orbital behavior,one can make adjustment to σ and test new orbital elements σ (cid:48) from the same t , andthen repeat the process to ensure that the stability requirements can be met. § Other small effects, such as higher-order ( >
10) Earth gravity, the Earth tides, the non-sphericalgravity of the Moon and the Sun, have been tested or estimated to be negligible to the optimizationresults (Table 4), hence not included for the sake of computational efficiency. ptimizing orbits for TianQin
Step 1 : For each spacecraft, we use the iterative relation (see Appendix B forderivation) r (cid:48) = (cid:18) ε ε ¯ a (cid:48) − ¯ a ¯ a (cid:19) r , v (cid:48) = (cid:18) − ε ε ¯ a (cid:48) − ¯ a a (cid:19) v , (2)to have the average semi-major axis approaching the desired value ¯ a (cid:48) = 10 km. Here( r , v ) are the initial position and velocity before adjustment, and ( r (cid:48) , v (cid:48) ) the newones. The value of ¯ a is determined by averaging the semi-major axis over the entireorbit generated from ( r , v ). Additionally, we define the parameter ε = (¯ a − a ) /a .By applying the relation (2) repeatedly, one can eliminate the long-term linear drift inarm lengths and breathing angles.Furthermore, we use the iterations (see Appendix C for derivation) i (cid:48) = (cid:18) (cid:15) (cid:15) ¯ i (cid:48) − ¯ i ¯ i (cid:19) i , Ω (cid:48) = Ω + (cid:0) ¯Ω (cid:48) − ¯Ω (cid:1) , (3)to set the three spacecraft on the same average orbital plane. The primed and un-primednotation above is interpreted similarly as in equation (2), and (cid:15) = (¯ i − i ) /i . Step 2 : To further improve the result from the previous step, we use numericaloptimization to minimize the following cost function CF [12]: CF = k CF + k CF , (4) CF ≡ c (cid:90) t f t ( | v | + | v | + | v | ) d t,CF ≡ c (cid:90) t f t (cid:0) ( α − ◦ ) + ( α − ◦ ) + ( α − ◦ ) (cid:1) d t, (5)with weights k and k (typically, 0 .
5) and normalization constants c and c . The formof CF is to bring down the relative velocities | v ij | and confine the breathing angles α k close to 60 ◦ . The imposed constraints are directly taken from the stability requirements,i.e., | v ij | (cid:54) | α k | (cid:54) . ◦ for the first two years, and loosened to | v ij | (cid:54)
10 m/sand | α k | (cid:54) . ◦ for the following three years. The independent variables one may varyinclude the initial eccentricities, arguments of periapsis, and true anomalies of the threespacecraft.Minimizing the cost function CF can help reduce the average eccentricities, which,as our simulations have shown, strongly affect the formation stability (cf. Appendix A)in the long run. To achieve better results in numerical search, the step 2 may berepeated. Step 3 : Redo the iteration (2) from the first step if the semi-major axes stray from10 km after the step 2. ptimizing orbits for TianQin
5. Optimized orbits for TianQin
The initial orbital elements obtained from optimization are listed in table 3. We present,in figure 2, the 5-year evolutions of the arm lengths L (black), L (blue), L (red),and the relative velocities v (black), v (blue), v (red), and the breathing angles α (black), α (blue), α (red), as well as the pointing deviation φ from J0806 (the anglebetween the direction to J0806 and the normal direction of the triangle). The result hasbeen independently verified by other orbit simulators. For more details, we summarizethe stability performance in table 4. One can see that the optimized TianQin orbitsfulfill the stability requirements provided in table 1. Table 3.
The initial elements of the optimized TianQin orbits in the EarthMJ2000Ec (Keplerian,ecliptic) and EarthMJ2000Eq (Cartesian, equatorial) coordinates at the epoch 22 May, 203412:00:00 UTC. The subsequent orbital evolution is shown in figure 2. a (km) e i ( ◦ ) Ω ( ◦ ) ω ( ◦ ) ν ( ◦ )SC1 99 995.572 323 0.000 430 94.697 997 210.445 892 358.624 463 61.329 603SC2 100 011.400 095 0.000 000 94.704 363 210.440 199 0.000 000 179.930 706SC3 99 993.041 899 0.000 306 94.709 747 210.444 582 0.001 624 299.912 164 x (km) y (km) z (km) v x (km/s) v y (km/s) v z (km/s)SC1 −
46 746.087 307 −
51 973.844 583 71 473.835 818 1.448 401 0.471 646 1.291 321SC2 86 220.582 041 46 448.360 669 20 269.217 366 0.085 035 0.663 048 − −
39 378.654 985 5547.379 475 −
91 728.424 823 − − Table 4.
A summary of the optimization results regarding their average orbital planes over 5years and stability performance of the constellation.Orbital plane Stability for 5 years (and the first 2 years)Result ¯Ω ( ◦ ) ¯ i ( ◦ ) | ∆ l ij | max (%) | v ij | max (m/s) | ∆ α k | max ( ◦ ) ¯ φ +( φ max − ¯ φ ) − ( ¯ φ − φ min ) ( ◦ )TianQin a . +1 . − . (0 . +0 . − . )P1 210.18 91.63 0.156 (0.098) 4.993 (4.130) 0.150 (0.090) 3 . +0 . − . (3 . +0 . − . )P2 240.00 88.98 0.151 (0.125) 5.260 (4.626) 0.139 (0.098) 30 . +0 . − . (30 . +0 . − . )P3 270.00 86.58 0.164 (0.126) 6.005 (4.793) 0.160 (0.102) 60 . +0 . − . (60 . +0 . − . )P4 120.00 100.00 0.148 (0.131) 5.423 (4.319) 0.132 (0.102) 88 . +1 . − . (88 . +1 . − . )P5 b . +1 . − . (61 . +0 . − . )P6 c . +0 . − . (30 . +0 . − . ) a Initial epoch 22 May, 2034 12:00:00 UTC; for the others, 01 January, 2034 00:00:00 UTC. b Retrograde orbits. c The detector pointing deviates from the Galactic Center by 6 . +0 . − . (6 . +0 . − . ) degree. ptimizing orbits for TianQin Figure 2.
Evolution of constellation geometry of the optimized TianQin orbits generated fromthe initial elements of table 3. The plots show 5-year variations of the arm lengths L (black), L (blue), L (red), and the relative velocities v (black), v (blue), v (red), and the breathingangles α (black), α (blue), α (red), and the pointing deviation from J0806.
6. Secular pointing drift
In figure 2 (the lower right panel), one observes an overall increase in the deviationangle φ away from J0806. In what follows, we show that the pointing shift is mainlycaused by orbital precession under the combined influence of lunisolar and the Earth’s J oblateness perturbations.In terms of the elements i and Ω, the unit normal vector ˆ n of an orbital plane canbe expressed asˆ n = sin i sin Ω − sin i cos Ωcos i . (6)Thus one can calculate the shift angle from ˆ n ( i , Ω ) to ˆ n t ( i ( t ) , Ω( t )) according tocos φ ( t ) = ˆ n t · ˆ n with ˆ n pointing at J0806. Let Ω( t ) = Ω +∆Ω( t ) and i ( t ) = i +∆ i ( t ).Taking into account that both ∆ i and ∆Ω are small, we have an approximation φ ≈ ∆ i + ∆Ω sin i . (7)As figure 3 shows, ∆Ω outgrows ∆ i and sin i ∼
1. Therefore the deviation angle φ isdominated by ∆Ω in the long run ( | ∆Ω | max = 2 . ◦ , | ∆ i | max = 0 . ◦ ), which manifests as ptimizing orbits for TianQin
8a common trend in both figure 2 (lower right) and figure 3 (right). One can understandthe behavior of ∆Ω and ∆ i as a generic property that Ω possesses secular change, whilenot so much for i [29]. Our numerical tests over 85 ◦ ≤ i ≤ ◦ also confirm thisobservation. Figure 3.
Evolution of the inclinations and longitudes of the ascending nodes of the threespacecraft (marked black, blue, and red), for the optimized TianQin orbits of table 3. Bothcurves appear blue due to extensive overlap.
From perturbative analysis, the mean rate of secular change of Ω is given by [29](also [30], see p 614)˙Ω m/s = − GM m/s (1 − e ) (1 + 5 e )4 a m/s n (1 − e m/s ) / cos i mo (8)for the third-body perturbation from the Moon ( m ) or the Sun ( s ). Here M m/s denotesthe mass of the third body, and n the mean angular motion of the spacecraft. All theorbital elements, e.g., e and a m/s , take on mean values over certain periods (e.g., onelunar month for the Moon). Note that i mo represents the mean inclination to the thirdbody’s orbital plane. Likewise for the Earth’s J perturbation, one has [31] (see p 139)˙Ω J = − J R ⊕ n a (1 − e )] cos i eq , (9)where i eq denotes the mean inclination to the Earth’s equator, and R ⊕ the Earth’sequatorial radius and J = 1 . × − . The equation is commonly found in thecontext of Sun-synchronous orbits.The formulas (8) and (9) give rise to a total secular change of 2 . ◦ in Ω, whichcombines the contributions from the Moon (+1 . ◦ ), the Sun (+1 . ◦ ), and the J oblateness ( − . ◦ ). The value agrees with the numerical result | ∆Ω | max = 2 . ◦ within8%. It confirms that the long-term pointing shift φ in figure 2 is primarily driven bylunisolar gravitational perturbations.One can see from (8) and (9) that there will be no orbital precession if i mo/eq = 90 ◦ .However, the equality cannot be achieved simultaneously for the three perturbingbodies. Because the Moon’s effect account for the largest one, we can arrange themean orbital plane perpendicular to the mean lunar orbital plane to reduce | ∆Ω | max , ptimizing orbits for TianQin | ∆ φ | max as well. This is demonstrated in the next section (see thelast plot of figure 4).
7. Optimized orbits for other pointings
From the previous analysis, we have found that the orientation of the TianQin orbitalplane varies about 2 . ◦ during a 5-year mission lifetime, indicating a rather stabledetector pointing at J0806. Considering the future possibilities of new reference sources,we present the optimization results for other pointing directions, in this section.As mentioned before, the long-term stability of the constellation depends mainlyon the magnitude of the average eccentricities of the orbits. From our numerical tests,we have noticed that the average eccentricity can attain small values if one sets theorbital plane roughly perpendicular to the mean orbital plane of the Moon. Thereforeto ease our search in optimization, we consider six detector pointing directions P1-6, allapproximately aligned with the lunar orbital plane, and 30 ◦ apart to spread over a halfcircle. Regarding the mean lunar orbital plane (from 1 January, 2034 to 1 January, 2039)in the EarthMJ2000Ec coordinates, it has the inclination 5 .
161 139 ◦ and the longitudeof ascending node 138 .
590 584 ◦ [32]. Our results are shown in table 5 and figures 4-9.Their stability performance is summarized in table 4.For P1, we point out that the orbital configuration resembles TianQin’s (see table4). The main difference lies in that the long-term growth of the deviation angle (figure2) is suppressed in P1, and taken over by semi-annual fluctuation of 0 . ◦ due to theinclination (figure 4).In the cases of P3 and P4 (figures 6, 7), the time evolution of the breathing anglescan wander over ± . ◦ in the first two years, but only at a few peaks and by a smallamount (0 . ◦ , see table 4). Hence we still include them for future consideration.Unfortunately in these two cases, further improvement to suppress the angle excursionwithin ± . ◦ appears difficult and time-consuming.Regarding P6 (figure 9), it is worth noting that the orbital pointing differs fromthe direction to the Galactic Center (the compact astronomical radio source SagittariusA*) for only about 6 . ◦ (see table 4), the smallest among all the cases.
8. Conclusion
In order to achieve a successful mission, it is vital for the TianQin constellation toattain high stability in orbit design and operation. In this study, we apply optimizationmethods and manage to stabilize the TianQin constellation down to the level of ± . ± ± . ◦ in breathing angles, for aperiod of 2 years along free-fall orbits. The optimized orbit configuration fulfills, withmargins, the assumed 5-year stability requirements from instrumentation. Though nodirect constraint is imposed, the detector pointing can be made quite stable, but asmall amount of pointing variation ( < . ◦ in 5 years) is unavoidable due to third-body ptimizing orbits for TianQin Figure 4.
Constellation evolution of the optimized orbits P1 (see figure 2 for the colorassignment).
Figure 5.
Constellation evolution of the optimized orbits P2 (see figure 2 for the colorassignment). ptimizing orbits for TianQin Table 5.
The initial elements of the optimized orbits with the pointings P1-6 in theEarthMJ2000Ec coordinates and at the epoch 01 January 2034 00:00:00 UTC. Their timeevolutions correspond to figures 4-9, respectively.Ptn a (km) e i ( ◦ ) Ω ( ◦ ) ω ( ◦ ) ν ( ◦ )P1 SC1 99 988.451 891 0.000 000 91.447 831 210.436 268 0.000 000 60.019 828SC2 100 047.700 990 0.000 804 91.444 569 210.444 886 180.451 960 359.620 218SC3 99 985.159 485 0.000 829 91.445 774 210.436 648 84.711 013 215.305 422P2 SC1 99 985.313 256 0.000 694 88.977 641 240.599 044 319.646 851 100.370 180SC2 100 063.708 616 0.000 544 88.974 726 240.598 772 179.927 490 0.012 642SC3 99 975.723 887 0.000 928 88.984 674 240.600 376 57.832 469 242.028 312P3 SC1 99 992.403 653 0.000 644 86.746 616 270.774 164 314.121 261 103.509 222SC2 100 033.833 610 0.000 000 86.749 894 270.770 211 0.312 080 177.221 271SC3 99 983.589 342 0.000 703 86.746 990 270.776 583 63.576 742 233.967 848P4 SC1 99 984.187 480 0.000 607 100.086 542 117.446 322 289.515 853 174.493 735SC2 100 007.078 264 0.000 232 100.082 574 117.432 353 226.666 677 357.274 707SC3 100 008.968 474 0.000 208 100.084 516 117.444 713 0.056 315 343.906 747P5 SC1 99 993.147 430 0.000 091 89.995 647 328.718 323 0.024 628 60.172 294SC2 100 011.119 344 0.000 274 89.989 041 328.724 788 234.731 987 305.458 669SC3 99 995.665 243 0.000 000 89.984 838 328.717 370 36.310 671 263.883 611P6 SC1 100 000.269 197 0.000 003 89.751 362 181.247 381 131.469 664 180.050 283SC2 99 994.745 911 0.000 027 89.759 135 181.250 998 359.603 358 71.845 840SC3 100 011.728 446 0.000 769 89.754 793 181.247 289 190.317 712 1.173 907 perturbations. We also consider six other detector pointings (P1-6) spreading over themean lunar orbital plane, and present the corresponding optimized orbits, which onemay considered as backups or alternatives for the current design. Particularly, it allowsa possibility of adding a second constellation of three spacecraft (e.g., P4) with theorbital plane roughly perpendicular to the first constellation. In this way, year-roundGW observation can be arranged without interruption due to sunlight.In future studies, employment of orbital corrections that can restrain relativemotion between satellites for longer periods will be investigated. The maneuvers cantake advantage of the semi-annual transition periods (3 months or less) when sunlightis roughly aligned with the orbital plane and interferes with GW observation [1].Furthermore, the requirement on delivery accuracy into the target orbits is being workedout, and the preliminary estimation has shown promise. Combined schemes, such asinter-satellite laser ranging, Chinese Deep Space Network, BeiDou/GPS, satellite laserranging, etc., will be considered to assess the future orbit determination capability. Acknowledgements
The authors thank Gang Wang, Shoucun Hu, Yi-Ming Hu, and Hsien-Chi Yeh for helpfuldiscussion. Our gratitude extends to the developers of GMAT. The work is supported ptimizing orbits for TianQin Figure 6.
Constellation evolution of the optimized orbits P3 (see figure 2 for the colorassignment). by NSFC 11805287, 91636111, 11690022, 11475064, 11690021, 11503007, and 41274041.
Appendix A. Estimation on eccentricities
To provide some intuition, one can use a two-body model (see also [33]) to roughlyestimate the stability requirements on the eccentricities of the spacecraft orbits. Tobegin with, we assume that the Keplerian orbits in the same orbital plane are given by X k = a (cos ψ k − e ) , Y k = a √ − e sin ψ k , k = 1 , , . (A.1)The eccentric anomaly ψ k is determined from the mean anomaly M k and the eccentricity e : ψ k − e sin ψ k = M k , M k ≡ n t − M k , (A.2)with the mean angular motion n , the time t , and a constant M k . Thereby the eccentricanomaly can be obtained as ψ k = M k + e sin M k + e cos M k sin M k + O ( e ) . (A.3)Moreover, we haveM k = n t + 2 π k − , (A.4) ptimizing orbits for TianQin Figure 7.
Constellation evolution of the optimized orbits P4 (see figure 2 for the colorassignment). We request the deviation angle ≤ ◦ , hence the broken line in the lower rightplot. which follows from assuming equal arm lengths L = L = L when e = 0. Forthese elliptic orbits, we can infer that the formation stability is closely related to theeccentricity by the following relations: | ∆ l ij ( t ) | max ≡ (cid:12)(cid:12)(cid:12)(cid:12) L ij ( t ) − L L (cid:12)(cid:12)(cid:12)(cid:12) max ≈ e ≈ e . × , | v ij ( t ) | max ≡ (cid:12)(cid:12)(cid:12) ˙ L ij ( t ) (cid:12)(cid:12)(cid:12) max ≈ √ an e ≈ (cid:114) a e . ×
10 (m/s) , | ∆ α k ( t ) | max ≡ (cid:12)(cid:12)(cid:12) α k ( t ) − π (cid:12)(cid:12)(cid:12) max ≈ √ e ≈ e . × . ◦ ) , (A.5)with L = √ a . The estimate e ∼ − , as an upper bound, is consistent with ouroptimized orbits from previous sections. Appendix B. Derivation of equations (2)
To remove long-term drift in arm lengths and breathing angles, one needs to make surethat the mean angular velocities, or equivalently, the mean semi-major axes by Kepler’slaw, are the same for the three spacecraft (see, e.g., [19, 28]).For real-world, perturbed orbits at a given moment, the Keplerian description canstill apply. A perturbed orbit can be described by a set of mean elements ¯ σ , together ptimizing orbits for TianQin Figure 8.
Constellation evolution of the optimized orbits P5 (see figure 2 for the colorassignment).
Figure 9.
Constellation evolution of the optimized orbits P6 (see figure 2 for the colorassignment). ptimizing orbits for TianQin σ c and periodic terms σ ls (including both long-periodic and short-periodic terms) [30, 31]: σ ( t ) = ¯ σ ( t ) + σ ls ( t ) , ¯ σ ( t ) = ¯ σ + η ¯ n ( t − t ) + σ c ( t ) , ¯ σ = σ − σ ls ( t ) , (B.1)where σ ∈ { a, e, i, Ω , ω, M } , and η = 1 if σ = M, and η = 0 if σ (cid:54) = M, and ¯ σ denotesthe initial mean elements.On the form of σ ( t ), the secular terms are composed of linear functions orpolynomials of ( t − t ), and the long-periodic terms of trigonometric functions of (Ω, ω ), and the short-periodic terms of trigonometric functions of M. In addition, theircoefficients are functions of ( a , e , i ) [31] (Liu 2000, pp 96, 126, 268). For example, wehave a ls ( t ) ∼ F ( a, e, i ) · F (Ω , ω, M) , (B.2)Ω c ( t ) ∼ F ( a, e, i ) · ( t − t ) , (B.3)with some functions F , F , and F (also F ( e, i, Ω , ω, M) in equation (B.6)).In the case of conservative perturbations, one has a c ( t ) = 0 for the semi-major axis a ( t ) [29–31] (Liu 2000, p 117; Vallado 1997, p 588). From equation (B.1), we obtain themean element ¯ a ( t ) as¯ a ( t ) = a − a ls ( t ) . (B.4)Taking an variation of (B.4), we arrive at δ ¯ a = δa − δa ls ( t ) (B.5)with δa = a (cid:48) − a . From perturbative analysis [31] (Liu2000, pp 126, 268), one has a ls ( t ) ∼ a γ F ( e, i, Ω , ω, M) , (B.6)where γ = − J perturbation and γ = 4 for the lunisolar perturbations.Because the lunisolar effects are far greater in magnitude than the J perturbation, weadopt the relation a ls ( t ) ≈ a F ( e, i, Ω , ω, M). Therefore, we have δa ls ( t ) = ∂a ls ∂a δa ≈ a ls ( t ) a δa , (B.7)assuming that all other elements are fixed. One can regard ¯ a as a perturbation of a :¯ a = a + εa , (B.8)with | ε | ∼ − for a geocentric spacecraft orbiting at an altitude of ∼ km.Substituting equations (B.4) and (B.8) into equation (B.7), we obtain δa ls ( t ) ≈ − εδa , (B.9)and by equation (B.5), δ ¯ a ≈ (1 + 4 ε ) δa . (B.10)Combining (B.8) and (B.10), we find δa a ≈ ε ε δ ¯ a ¯ a , a (cid:48) ≈ (cid:18) ε ε ¯ a (cid:48) − ¯ a ¯ a (cid:19) a , (B.11) ptimizing orbits for TianQin ε = (¯ a − a ) /a from (B.8).In practice, we would prefer using the equation (B.11) in terms of the spacecraft’sposition and velocity, as in [19,20]. As mentioned before, geometric relations in ellipticalmotion can still apply to a perturbed Keplerian orbit, e.g., r = a (1 − e cos ψ ) , (B.12)where each quantity now varies with time. Considering the case with δe = 0 and δ M = 0,we obtain from (B.12) δrr = δaa . (B.13)For perturbed Keplerian orbits, the kinetic and potential energy is related to the semi-major axis by mv − GM ⊕ mr = − GM ⊕ m a (B.14)at each given time. Taking the variation of (B.14) gives vδv = GM ⊕ a δa − GM ⊕ r δr. (B.15)Combining equations (B.13), (B.14) and (B.15), we obtain δvv = − δa a . (B.16)By equations (B.11), (B.13) and (B.16), we find δ r r = δr r = δa a ≈ ε ε δ ¯ a ¯ a , δ v v ≈ − ε ε δ ¯ a a , (B.17)or, r (cid:48) ≈ (cid:18) ε ε ¯ a (cid:48) − ¯ a ¯ a (cid:19) r , v (cid:48) ≈ (cid:18) − ε ε ¯ a (cid:48) − ¯ a a (cid:19) v , (B.18)which completes the derivation of equations (2). Appendix C. Derivation of equations (3)
We also adjust the orbital planes of the three spacecraft so that they can stay in thesame average plane. The inclination i and the longitude of ascending node Ω determinethe orientation of the orbital plane. For the inclination, it has little secular change.Hence, similar to the treatment of a , we apply the same form of the iteration (B.11) on i , i.e., i (cid:48) ≈ (cid:18) (cid:15) (cid:15) ¯ i (cid:48) − ¯ i ¯ i (cid:19) i , (C.1)with the average inclination ¯ i , the desired inclination ¯ i (cid:48) , and (cid:15) = (¯ i − i ) /i . For thelongitude of ascending node, we know from (B.1) that¯Ω( t ) = Ω − Ω ls ( t ) + Ω c ( t ) . (C.2) ptimizing orbits for TianQin T , one has¯Ω T ≡ T (cid:90) t + Tt ¯Ω( t ) dt ≈ Ω + 1 T (cid:90) t + Tt Ω c ( t ) dt, (C.3)where the small contribution from Ω ls ( t ) is neglected. By the relation (B.3) itfollows that δ ¯Ω T ≈ δ Ω . Using ¯Ω( t ) = Ω( t ) − Ω ls ( t ) from (B.1), we arrive at¯Ω ≡ T (cid:82) t + Tt Ω( t ) dt = ¯Ω T andΩ (cid:48) ≈ Ω + (cid:16) ¯Ω (cid:48) − ¯Ω (cid:17) , (C.4)where ¯Ω (cid:48) is the desired mean longitude of ascending node. References [1] Luo J et al
Class. Quantum Grav. Class. Quantum Grav. S809–17[3] Hu Y-M, Mei J-W and Luo J 2017 Science prospects for space-borne gravitational-wave missions
National Science Review Satellite Orbits: Models, Methods, and Applications (New York:Springer) p 114[5] Vincent M A and Bender P L 1987
Proc. Astrodynamics Specialist Conf. (Kalispell) vol 1 (SanDiego, CA: Univelt) p 1346[6] Folkner W M, Hechler F, Sweetser T H, Vincent M A and Bender P L 1997 LISA orbit selectionand stability
Class. Quantum Grav. Class. QuantumGrav. S429–35[8] Dhurandhar S V, Nayak K R, Koshti S and Vinet J-Y 2005 Fundamentals of the LISA stable flightformation
Class. Quantum Grav. Class. Quantum Grav. Class.Quantum Grav. Int. J. Mod. Phys. D [13] Povoleri A and Kemble S 2006 LISA orbits AIP Conf. Proc.
Int. J. Mod. Phys. D Sci. China Phys. Mech. Astron. J. Phys.: Conf. Ser.
Class. Quantum Grav. Int. J. Mod. Phys. D Chin. Astron. Astrophys. ptimizing orbits for TianQin [20] Wang G and Ni W-T 2013 Orbit optimization for ASTROD-GW and its time delay interferometrywith two arms using CGC ephemeris Chin. Phys. B Chin. Phys. B Int. J. Mod. Phys. D AIAA/AAS Astrodynamics Specialist Conference,AIAA SPACE Forum , (AIAA 2014-4151)[26] Tapley B D et al
J. Geophys. Res.
Chinese Space Science and Technology Planet. SpaceSci. Fundamentals of Astrodynamics and Applications (New York: McGraw Hill)[31] Liu L 2000
Orbit Theory of Spacecraft (in Chinese) (Beijing: National Defence Industry Press)[32] https://ssd.jpl.nasa.gov/horizons.cgi[33] Hu X-C, Li X-H, Wang Y, Feng W-F, Zhou M-Y, Hu Y-M, Hu S-C, Mei J-W and Shao C-G 2018Fundamentals of the orbit and response for TianQin
Class. Quantum Grav.35