Issue 
A&A
Volume 682, February 2024



Article Number  A38  
Number of page(s)  18  
Section  Astronomical instrumentation  
DOI  https://doi.org/10.1051/00046361/202348218  
Published online  02 February 2024 
Formationflying interferometry in geocentric orbits
Institute of Space and Astronautical Science,
311 Yoshinodai, Chuoward, Sagamihara,
Kanagawa
2520222,
Japan
email: ito.takahiro@jaxa.jp
Received:
10
October
2023
Accepted:
16
November
2023
Context. Spacecraft formation flying serves as a method of astronomical instrumentation that enables the construction of large virtual structures in space. The formationflying interferometry generally requires very high control accuracy, and extraterrestrial orbits are typically selected. To pave the way for comprehensive missions, proposals have been made for preliminary space missions, involving nano or small satellites, to demonstrate formationflying interferometry technologies, especially in low Earth orbits. From a theoretical perspective, however, it is unknown where and to what extent feasible regions for formationflying interferometry should exist in geocentric orbits.
Aims. This study aims to demonstrate the feasibility of formationflying interferometry in geocentric orbits in which various perturbation sources exist. Geocentric orbits offer the advantage of economic accessibility and the availability of proven formationflying technologies tailored for Earth orbits. Its feasibility depends on the existence of specific orbits that satisfy a smalldisturbance environment with good observation conditions.
Methods. Spacecraft motions in Earth orbits subjected to perturbations are analytically modeled based on celestial mechanics. The magnitudes of the accelerations required to counteract these perturbations are characterized by parameters such as the semimajor axis and the size of the formation.
Results. Smallperturbation regions tend to appear in higheraltitude and shorterseparation regions in geocentric orbits. Candidate orbits are identified in high Earth orbits for the triangular laserinterferometric gravitationalwave telescope, which is 100 km in size, and in medium Earth orbits for the linear astronomical interferometer, which is 0.5 km in size. A low Earth orbit with a separation of approximately 0.1 km may be suitable for experimental purposes.
Conclusions. Geocentric orbits are potentially applicable for various types of formationflying interferometry.
Key words: gravitation / gravitational waves / instrumentation: interferometers / space vehicles / telescopes / celestial mechanics
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1 Introduction
Spacecraft formation flying represents a crucial method in astronomical instrumentation, with the potential to revolutionize modern astronomy by creating large virtual structures in space. Some spacebased astronomical missions under conceptual study require rigid and continuous formation while maintaining a very high control accuracy of less than 1 μm and/or μrad, and interferometric techniques are typically used. The Decihertz Interferometer Gravitational Wave Observatory (DECIGO; Kawamura et al. 2021) is a laserinterferometric, spacebased gravitationalwave observatory used to measure coalescing binary neutron stars, the acceleration of the universe, and stochastic gravitational waves predicted by inflation. Each cluster comprises three spacecraft that are 1000 km apart and form an equilateral triangle, and four clusters are placed in a nearcircular heliocentric orbit. FabryPérot cavities are employed in DECIGO, each comprising two mirrors to measure the armlength changes. The Large Interferometer For Exoplanets (LIFE; Quanz et al. 2022) is a spacebased midinfrared nulling interferometer that detects and characterizes the thermal emissions of exoplanets. LIFE is based on earlier mission concepts, such as Darwin (Cockell et al. 2009) and the Terrestrial Planet Finder Interferometer (Lawson & Dooley 2005), and its primary orbital selection is an orbit around the SunEarth Lagrange second (L2) point. The current design for LIFE assumes four collector spacecrafts and one beamcombiner spacecraft, and their formation size is controlled to remain between tens and hundreds of meters. The common features of DECIGO and LIFE are that precise control is necessary to maintain the optical path length or difference between satellites and that their candidate orbits are beyond Earth’s orbit.
Flightproven examples in formation flying have been shown in geocentric orbits, and coarse formation control on the order of submeters to meters was typically sufficient to satisfy mission objectives. The Japanese Engineering Test SatelliteVII (Kawano et al. 2001) demonstrated autonomous formation flying, rendezvous, and docking in the late 1990s in a circular low Earth orbit (LEO). In the 2000s and 2010s, autonomous formationflying technology advanced for circular LEOs through leading space missions such as Prototype Research Instruments and Space Mission technology Advancement (PRISMA; Gill et al. 2007) and TerraSARX/TanDEMX (TerraSARX addon for Digital Elevation Measurement; Kahle et al. 2014) for technological demonstration and Earth observation. Their accuracy during formation flying reached the centimeter to meter level based on the navigation accuracy provided by the global navigation satellite system (GNSS). In the 2020s, Project for OnBoard Autonomy3 (PROBA3) will demonstrate several formation activities, including coarsetofine control, autonomy, and safety management, to observe the solar coronagraph in a highly elliptical orbit. Its formationflying precision is capable of submillimeter relative displacement and arcsecondorder pointing accuracy. This accuracy will be achieved using stateoftheart metrology sensors and millinewtonclass thrusters (Penin et al. 2020). As seen in these examples, formation flight technologies are gradually maturing; however, a technological gap still exists in the development of future formationflying interferometry.
Systemlevel, endtoend, inspace verification for the pathfinder may be necessary to mature formationflying interferometry before fullscale missions, as suggested in Monnier et al. (2019). However, the use of extraterrestrial orbits may be an inherent bottleneck for this purpose. The major benefits of using extraterrestrial orbits are the smalldisturbance environment and good observation conditions. However, its major limitations are the high transportation cost to the desired orbit and that GNSS cannot be used for formationflying autonomy, safety, and management. Such drawbacks may lead to formationflying interferometry becoming a larger and longerterm project.
An alternative is to use a geocentric orbit in which various perturbation sources exist. The expected advantages of using Earth orbits are the economic accessibility and availability of flightproven technologies for formation flying autonomy, safety, and management with GNSS. Furthermore, stateoftheart technologies are extending the use of GNSS from LEO to high Earth orbit (Ashman et al. 2018). However, few studies have applied Earth orbits to formationflying interferometry. Representative and stateoftheart research in this context conducted by Hansen & Ireland (2020) involved a linear astronomical interferometer in a nearcircular LEO, and Earth’s J_{2} gravity potential and atmospheric drag were considered to be major perturbations. The necessary ΔV (the time integral of the thrust acceleration) to maintain the strict formation requirements was at an acceptable level of a few tens of m s^{−1} per year, as calculated based on the numerical integration of the equations of motion. However, the numerical approach tends to rely on trial and error to determine feasible Earthorbit regions. From a theoretical perspective, understanding where and to what extent feasible regions should exist is difficult. The characterization of orbital perturbations is as crucial for formationflying interferometry as site characterization is for the deployment of groundbased telescopes.
This study investigates the feasibility of formationflying interferometry in geocentric orbits based on a celestialmechanics approach. This study makes four major contributions to the literature. The first contribution is the analytical modeling of the perturbed orbital motion to attain precise formation flying. We identify four types of absolute and relative perturbations that need to be mitigated. The second contribution is to offer a propellantefficient control strategy for correcting the absolute and/or relative perturbing accelerations. Secular and longperiod perturbations of the eccentricity vector play key roles in determining the control strategy. The third contribution enables a complete and explicit calculation of major perturbing accelerations. Using the proposed method, the disturbance environments for different formations can be analytically evaluated in various Earth orbits. By integrating these three contributions, the perturbation magnitudes can be characterized by the semimajor axis and size of the formation, and this information can be used as a guideline for finding candidate orbits. The final contribution is the provision of detailed examples of orbital selection and control approaches for formationflying interferometry. The proposed approach as well as methods for satisfying the observation conditions will be useful for future applications.
This study considers only a nearcircular Earth orbit because such is sufficient for various potential applications in formationflying interferometry. In terms of terminology, this study uses the chiefdeputy expression to describe formationflying satellites. A chief does not necessarily represent the actual chief satellite, but rather the virtual origin of formationflying satellites. By contrast, a deputy denotes a real spacecraft flying in proximity to the chief. These terms are used throughout the paper without any notation. Lastly, this study uses numerous symbols owing to the thorough analysis of perturbed motion. The major symbols, particularly those used in the theoretical development, are listed in Appendix A.
The remainder of this paper is organized as follows. Section 2 presents a review of unperturbed relative motion for a circular orbit. Section 3 presents the development of a perturbed relativemotion model for a nearcircular orbit against arbitrary perturbations and explains the four types of resultant perturbing accelerations. Section 4 first presents an analysis of the relationship between the perturbation periods and their contributions to perturbing acceleration magnitudes. Then, it presents analytical models of the perturbed orbital motions and perturbing accelerations as well as an effective approach for orbital selection and control. Section 5 presents numerical examples of the proposed approach to formationflying interferometry. The concluding remarks are presented in Sect. 6.
2 Unperturbed relative motion
The relative equation of motion of the deputy with respect to the chief in the Earthcentered inertial (ECI) frame is (1)
where μ_{e} is the Earth gravitational constant, r^{I} is the relative position vector of the deputy and are the position vectors of the chief and deputy, respectively, is the disturbance acceleration vector of the deputy relative to the chief and are the disturbance acceleration vectors for the chief and deputy, respectively, and are the control acceleration vectors applied to the chief and deputy, respectively, and . The superscript 1 indicates the position, velocity, or acceleration in the ECI frame.
By assuming the unperturbed motion of both the chief and deputy and the uncontrolled motion of the chief , Eq. (1) in the ECI frame is transformed into the localverticallocalhorizontal (LVLH) frame: (2)
where r = [x, y, z]^{T} denotes the position vector of the deputy relative to the chief in the LVLH frame, θ is the argument of latitude of the chief, and u_{d} = [u_{x}, u_{y}, u_{z}]^{T} is the control acceleration vector applied to the deputy in the LVLH frame. The position, velocity, or acceleration in the LVLH frame are written without superscript notations, unlike those in the ECI frame. The LVLH frame is defined such that the x direction is directed radially outward from the chief, the z direction is normal to the orbital plane and positive in the direction of the (instantaneous) angular momentum, and the y direction in the righthand frame completes the setup. Assuming r ≪ R_{c} (r = r), the first term on the righthand side of Eq. (2) can be linearized as follows: (3)
For a circular orbit, , and R_{c} = a, where n and a are the mean motion and semimajor axis of the chief, respectively. In addition, using Eqs. (2) and (3) and the relationship , the unperturbed and linear relative dynamics model for a circular orbit is obtained as follows: (4)
where A_{u1} and A_{u2} are matrices defined as follows: (5)
Equation (4) is known as the ClohessyWiltshire (CW) equation (Clohessy & Wiltshire 1960).
A possible approach for realizing formation flying under the orbital dynamics in Eq. (4) is to cancel the relative accelerations A_{u1} r and with the control accelerations. However, the control acceleration requirement tends to become large, which is inefficient and unrealistic for formationflying applications. Therefore, the reference formation orbit is typically selected using the natural solution of Eq. (4) when u_{d} = 0. The bounded CW solution (the solution without an alongtrack drift) is (6)
where ρ_{x}, ρ_{y}, ρ_{z}, α_{x}, and α_{z} denote the design parameters of the relative orbit (Clohessy & Wiltshire 1960). Clearly, θ is a linear function of time because for the unperturbed circular orbit; therefore, r_{r} is also a function of time.
3 Perturbed relative motion
In contrast to Keplerian dynamics, oscillating and drifting orbital motion can arise in the presence of orbital perturbations. The reference formation described in the previous section was obtained using unperturbed, circular, and linearized assumptions. Therefore, elliptical and nonlinear effects may disturb the actual formation. This section analyzes how orbital perturbations and unmodeled dynamics can affect the rigid and continuous formation control.
The perturbed relative equation of motion of the deputy in the LVLH frame is expressed as follows (Izzo et al. 2003): (7)
where ω = [ω_{x},ω_{y},ω_{z}]^{T} is the orbital angular velocity vector of the chief in the LVLH frame, f_{2B} is the linearized gravity gradient acceleration owing to the twobody gravitational field is the sum of the relative and physical perturbing accelerations applied to the deputy, and u_{c} and u_{d} are the control acceleration vectors applied to the chief and deputy in the LVLH frame, respectively. The elements of ω are (Alfriend et al. 2010) (8)
where Ω and I are the right ascension of the ascending node (RAAN) and inclination of the chief, respectively. We can verify ω_{y} = 0 by substituting I and from Gauss’s variational equations into Eq. (8) (Alfriend et al. 2010).
Next, Eq. (7) is transformed into (9)
We note that A_{t1} and A_{t2} are (3 × 3) matrices, where the elements of A_{u1} and A_{u2} in Eq. (5) are replaced by n → ω_{t}, where ω_{t} is the userdefined target angular velocity. In Eq. (9), f_{o} takes the following form: (10)
where A_{p1} and A_{p2} are (3 × 3) matrices corresponding to the first four terms in Eq. (7). It should be noted that f_{o} is the relative fictitious acceleration owing to the perturbed orbital motion of the chief, and it remains zero if the chief orbit is unperturbed and circular. Equation (10) suggests that the selection of the ω_{t} profile is important for minimizing f_{o}. Choosing ω_{t} close to ω_{z} is one of the best solutions, and Eq. (10) can be approximated further. Hereafter, ω_{t} is selected as ω_{z} throughout this study. First, under the nearcircularorbit assumption, the argument of latitude can be approximated as follows: (11)
where ϕ is the mean argument of latitude defined as ϕ = (l + ω), ω is the argument of perigee, l is the mean anomaly, e_{x} and e_{y} are the elements of the eccentricity vector defined as [e_{x}, e_{y}]^{T} = [e cos ω, e sin ω]^{T}, and e is the eccentricity of the chief. Then, let the time derivative of ϕ be expressed as . More over, the time derivatives , and are assumed sufficiently small with respect to n. By neglecting the second and higher order terms of the small variables, (A_{p1} – A_{t1}) and (A_{p2} – A_{t2}) in Eq. (10) can be expressed as follows: (12)
The following relationships are used to obtain Eqs. (12) and (13): (14)
Then, the desired reference formation is selected as the natural solution of Eq. (6): (15)
where Θx and Θz are the new variables satisfying Θx = (θ + α_{x}), , and , and the parameters of (ρ_{x}, ρ_{y}, ρ_{z}) can be timevariant. The relationships between , and can be written as follows: (16) (17)
We note that f_{f} is the relative fictitious acceleration owing to the perturbed motion of the reference formation.
Finally, by subtracting Eq. (16) from Eq. (9), we obtain the following: (18)
where ϵ = (r – r_{r}) is the relative positional error with respect to the reference formation. By setting the control command of the deputy to (19)
where u_{r} is the control acceleration to compensate for all relative disturbances, and by assuming that the current relative position and velocity are on the desired formation , the relative acceleration error becomes zero; thus, tracking on the reference formation is expected for the deputy.
Table 1 summarizes the four categories of perturbing accelerations that can be mitigated by the control. The simplest way to nullify the perturbation effects is to compensate for all absolute and relative perturbing accelerations (i.e., F_{c}, which is the expression of in the LVLH frame, f_{p} , f_{o} , and f_{f} ) using the control u_{d} . However, compensation for absolute perturbing accelerations typically leads to a significant waste of propellant. In such cases, another possible approach is to maintain most of the absolute perturbing accelerations, except for the small ones, and compensate for all the relative perturbations with respect to the chief. The feasibility of such a control strategy depends on whether the control acceleration u_{d} , and eventually F_{c} , f_{p} , f_{o} , and f_{f} are sufficiently small for the reference formation. Determining the differential disturbances of the deputy with respect to the chief is necessary to evaluate f_{p} . Moreover, determining the perturbed orbital motions of (e_{x}, e_{y}) as well as the first and second time derivatives of (e_{x}, e_{y}, I, Ω, ϕ) is necessary to evaluate f_{o}. The ω_{t}(= ω_{z}) profile affects f_{o} and f_{f}, and understanding the perturbed motions of ω_{z} and ω_{z}, which are functions of (e_{x}, e_{y} , I, Ω, ϕ), is essential. Furthermore, the controleffort tradeoff between compensating for or maintaining F_{c} remains unknown. Such concerns motivate the development of analytical models of perturbed orbital motions and perturbing accelerations, as described in the next section.
Categories of perturbing accelerations that can be mitigated by the control.
4 Perturbation mitigation methods
4.1 Perturbation periods and acceleration magnitudes
The physical perturbing accelerations (i.e., F_{c} and f_{p}) are determined simply from the positions and velocities of the chief and deputy, in addition to spacecraft and environmental parameters at a certain moment. In contrast, to calculate the fictitious perturbing accelerations (i.e., f_{o} and f_{f}), it is necessary to understand the perturbed orbital motion of the chief, as presented in the previous section. In addition, the control accelerations to compensate for the absolute physical perturbations (F_{c}) can be divided into the secular, longperiod, and shortperiod parts through Gauss’s variational equations. Whereas compensation for all parts of F_{c} could lead to a significant waste of propellant, compensation for small parts of F_{c} may be comparable to the relative fictitious perturbations. Therefore, before developing the perturbation models for each perturbing source, we show an analysis to understand how the perturbed orbital motion of the chief can affect the control accelerations to compensate for the relative fictitious perturbations, and how the relative fictitious perturbations are compared to the absolute physical ones with respect to the perturbation periods.
Equations (10), (12), and (13) suggest that the osculating orbital motion produces the relative fictitious acceleration of f_{o} , whose magnitude comprises the following terms: (·)n^{2}r_{r} for (e_{x}, e_{y}), and and for (e_{x}, e_{y}, I, Ω, ϕ), where (·) expresses the osculating orbital element and r_{r} = r_{r}. From Gauss’s variational equations, the control accelerations to compensate for the absolute perturbations comprise the following terms: for (e_{x}, e_{y}, I, Ω, ϕ). At this point, the osculating orbital elements are modeled as follows: (20)
where (·)^{(s)} is the secular term that is a linear function of time, δ(·)^{(lp)} is the longperiod term for which the angular frequency has an order of n_{l}(≪ n), δ(·)^{(sp)} is the shortperiod term for which the angular frequency has an order of n, is the initial value of the secular term, and is the drift rate of the secular term. By substituting Eq. (20) into relevant terms, the effects of the secular, longperiod, and shortperiod terms can be analyzed. The same discussion holds for the relative fictitious acceleration of f_{f}, provided that ρ_{x}, ρ_{y}, and ρ_{z} are derived from the natural solution in Eq. (6).
Table 2 summarizes the orders of the control acceleration magnitudes required to compensate for the relative fictitious perturbations (f_{o} and f_{f}) for the secular, longperiod, and shortperiod terms in comparison with those of the absolute physical perturbations (F_{c}). Table 2 suggests that most terms to compensate for the relative fictitious perturbations include the common terms of the absolute physical ones multiplied by (r_{r}/a) and/or (n_{l}/n), both of which are sufficiently smaller than 1. Therefore, the best approach is to compensate for the relative fictitious perturbations to gain less control acceleration. However, two exceptional cases appear in the secular and longperiod parts of the term with (·)n^{2}r_{r}. They have common terms of the absolute physical perturbations multiplied by not only (r_{r}/a) (≪ 1), but also (nt), which can be considerably larger than 1 when t is large, or (n_{l}/n)^{−1}, which is significantly larger than 1. This result implies that the secular and longperiod motions of the eccentricity vector (e_{x}, e_{y}) must be compared to determine whether to compensate for the absolute or relative perturbations to attain more efficient control. This insight is used to determine the orbit and formation control approach in Sect. 4.3.
Orders of the control acceleration magnitudes to compensate for the relative fictitious perturbations (f_{o} and f_{f}) for the secular, longperiod, and shortperiod terms, in comparison to those of the absolute physical perturbations (F_{c}).
4.2 Analytical perturbation models
The perturbation sources considered in this study are CW nonlinearity, Earth J2 gravity potential, Earth J_{3} gravity potential, lunisolar gravity, atmospheric drag, and SRP. Their analytical models are explained in the following.
4.2.1 CW nonlinearity
The relative acceleration owing to CW nonlinearity is produced by nonlinear terms that are neglected in the CW equation (Eq. (4)). This can be computed by subtracting the first term on the righthand side of Eq. (3) from the lefthand side. The CW nonlinearity can contribute to a part of f_{p} .
4.2.2 Earth J_{2} gravity potential
The J_{2} gravity potential of Earth’s gravity is known to be one of the major perturbation sources for satellites, particularly in lowtomedium Earth orbit. By considering only zonal effects, the J_{2} perturbed absolute acceleration vector acting on the chief ( , as a part of F_{c}) is given by (Alfriend et al. 2010) (21)
where R_{e} is the mean equatorial radius of Earth. The linearized J_{2} differential acceleration vector ( , as a part of f_{p}) has the following form (Izzo et al. 2003): (22)
A comparison of Eqs. (21) and (22) shows that is always less than under the assumption that r ≪ R_{c}.
The J_{2}perturbed motion of the chief can also yield the relative fictitious accelerations of f_{o} and f_{f} . The osculating orbital elements owing to the J_{2} potential for the nearcircular orbit are obtained by substituting e → 0 into the known analytical solutions (e.g., see Kozai 1959; Brouwer 1959), as given in Appendix B. Notably, the J_{2} perturbation produces a secular drift in the RAAN. This drift may be useful for the mission orbit to scan the visible sky area, as discussed in detail in Sect. 5.
4.2.3 Earth J_{3} gravity potential
The Earth J_{3} gravitational potential represents asymmetry with respect to the equatorial plane. Its coefficient is on the order of 10^{−6} , which is thee orders of magnitude smaller than that for the J_{2} term; therefore, most effects on the J_{3} term can be neglected in comparison to the J_{2}perturbed effects. The only effect to be considered is the longperiod J_{3} perturbation of the eccentricity vector in the presence of the J_{2} perturbation.
Assuming a small eccentricity and neglecting longitudedependent effects, Lagrange’s equations relating to the eccentricity vector for the longperiod J_{3} perturbation are as follows (Kozai 1959): (24)
where is the J_{2}perturbed secular drift rate of the argument of perigee given by for the nearcircular orbit. By applying the equilibrium condition of the eccentricity vector to an arbitrary inclination, the equilibrium eccentricity vector, known as a frozen eccentricity vector (e_{f}), is given as follows (Kozai 1959): (25)
Equation (25) implies that the mean eccentricity vector cannot be maintained at zero under the existence of Earth’s J_{2} and J_{3} gravity potentials, and its effect can be significant, particularly in LEO. For example, for a satellite in a nearcircular LEO at an altitude of 500 km and inclination of 100°, the frozen eccentricity becomes e_{f} = [0, 0.0011]^{T}. Such a nonzero eccentricity vector can contribute to the relative fictitious accelerations of f_{o} and f_{f}.
4.2.4 Lunisolar gravity
The absolute lunisolar gravity acceleration vector of the chief with respect to the Earth (F_{3rd}, as a part of F_{c}) is approximated as (Montenbruck & Gill 2000) (26)
where is the position vector of the third body in the ECI frame, μ_{b} denotes the gravity constant of the third body, is the rotation matrix from the ECI frame to the rotational (LVLH) frame of the chief, and .
The linearized gravity gradient acceleration vector owing to the third body (f_{3rd} , as a part of f_{p}) is (27)
Assume that r ≪ a; then, f_{3rd} is always less than F_{3rd}.
The perturbed motion of the chief by lunisolar gravity can also yield the relative fictitious accelerations of f_{o} and f_{f}; therefore, understanding its osculating motion is important. To the best of our knowledge, however, no appropriate analytical expressions describing the secular, longperiodic, and shortperiod motions in a nearcircular geocentric orbit due to lunisolar gravity have been proposed with reasonable accuracy for a relatively short timescale. Early work (Cook 1962) analytically calculated the secular and longperiod satellite motion due to lunisolar perturbations up to the second degree of Legendre polynomials, but shortperiod motions were not analyzed. Subsequent work (Kozai 1973) analytically proposed the calculation of shortperiod satellite motion caused by lunisolar perturbations, but secular and longperiod motions were calculated numerically. Recent work (Roscoe et al. 2013) developed secular and periodic motion models owing to thirdbody perturbations using a single and doubleaveraged disturbing function. Analytical models were developed based on the assumption of a circularrestricted threebody problem; however, this may be problematic, particularly when longterm orbital motion on the eccentricity vector is a major concern. It is known that the parallactic term (third harmonic) of the lunar disturbing function can change the eccentricity of a circular geocentric orbit (e.g., see the Introduction in Allan & Cook 1964) through the nonzero eccentricity of the lunar orbit. To capture this effect, the disturbing function of the third body is considered up to the third degree of the Legendre polynomials. Taking these points into account, this study analytically developed a relatively simple, nearcircular satellite motion model based on to lunisolar gravity, and the results are given in Appendix C.
In addition, the longterm effects of the lunisolar gravitational field and J_{2} term of the Earth’s oblateness may be comparable for a mediumtohigh Earth orbit. For this gravity field, equilibrium solutions exist in which the averaged chief orbits may be frozen. This plane is known as the classical Laplace plane, which is between the Earth’s equator and that of the ecliptic, passing the vernal equinox (Allan & Cook 1964; Tremaine et al. 2009; Rosengren et al. 2014). Assuming that the chief orbit is circular, the lunar orbit lies in the ecliptic, and the direction of the angular momentum for the chief lies on the principal plane defined by the Earth polar axis and normal direction of the Earth orbit (such that Ω = 0° or 180°), the frozen orbit condition is (Allan & Cook 1964) (28)
where ϵ_{e} is the obliquity of the ecliptic (ϵ_{e} = 23.4°), Φ is the azimuth angle, and , ω_{m}, and ω_{s} are defined as (29)
The subscripts m and s indicate the Moon or the Sun, respectively, and a and e with these subscripts represent the semimajor axis and eccentricity for the Moon or the Sun. For example, for a chief satellite at an altitude of 35 800 km (geosynchronous) with Ω = 0°, one can find that I = Φ = 7.4° is one of the solutions for the frozen orbits in Eq. (28). If the circular orbit does not lie in the classical Laplace plane, the orbital orientation will oscillate around the frozen point with a period of up to several tens of years. The frozen orbit lying in the classical Laplace plane may be useful particularly for the gravitationalwave detector, for obtaining the longterm stability on the orbital orientation with respect to the Sun, as discussed in Sect. 5.1.
4.2.5 Atmospheric drag
Atmospheric drag is a small force for most orbits, but it can be a major perturbation source for relative motion in LEOs. The atmosphericdrag acceleration vector acting on the chief spacecraft in the LVLH frame (FD, as a part of F_{c}) is as follows: (30)
where V is the velocity vector of the chief for the unperturbed and circular assumption in the LVLH frame (V = V= na), ρ is atmospheric density, B_{D}(= S C_{D}/m) is the coefficient for drag, S is the crosssectional area, C_{D} is the drag coefficient, and m is the spacecraft mass. Equation (30) is obtained by assuming that the atmosphere corotates with the Earth (Montenbruck & Gill 2000) and V is approximately one order of magnitude higher than that of the Earthfixed atmosphere; thus, atmospheric velocity is neglected.
In general, the atmospheric density has a large uncertainty, mainly owing to diurnal variation, periodic solar activity, and sudden solar activity (Montenbruck & Gill 2000). However, the most evident dependence of atmospheric density is that it decreases with increasing altitude. The HarrisPriester density model (Harris & Priester 1962) is relatively simple and is widely used as a standard atmospheric model. The model is given as follows: (31)
where h is the height above Earth’s reference ellipsoid, Ψ is the angle between the satellite position vector and the apex of the diurnal bulge, N is a numerical value set to 2 for low inclination orbits and 6 for polar orbits, and ρ_{M}(h) and ρ_{m}(h) are the apex and antapex densities, respectively. It should be noted that ρ(h) can be directly computed from the table of Harris Priester atmosphericdensity coefficients for mean solar activity (Long et al. 1989). When assuming a spherical Earth, height is given by h = (a − R_{e}); thus, ρ is written as a function of the semimajor axis.
Assuming that the two satellites are flying in proximity formation, such that V and ρ are the same for both spacecraft, the relative acceleration vector for the deputy caused by atmospheric drag ( f_{D} , as a part of f_{p}) can be modeled as follows: (32)
where δ B_{D} is the difference between the coefficients B_{D} of the two spacecraft for atmospheric drag.
In contrast to Earth J_{2} and lunisolar gravity, F_{D} and  f_{D} can be of the same order. In addition, drag is a small force in most Earth orbits. Therefore, we assumed that both F_{D} and f_{D} are compensated for by the control accelerations, and we neglected the effects of f_{o} and f_{f}.
4.2.6 Solar radiation pressure
Solar radiation pressure (SRP) can be regarded as a function of the Suntospacecraft distance, such that its magnitude is almost identical for spacecraft orbiting the Earth. A cannonball model was used, in which the satellite shape was assumed to be a uniform sphere. When the shade effect of the Earth and the small eccentricity of the Earth’s orbit are neglected, the acceleration vector by SRP acting on the chief spacecraft in the LVLH frame (F_{s}, as a part of F_{c}) is expressed as follows (Montenbruck & Gill 2000): (33)
where P_{s} is the SRP, C_{s} is the radiation pressure coefficient, B_{s} is the coefficient for SRP, and s is the unit vector pointing to the Sun in the LVLH frame. Assuming that two spacecraft are flying in proximity formation (P_{s} is the same for both spacecraft), the relative acceleration vector for the deputy caused by SRP (f_{s}, as a part of f_{p}) can be modeled as follows (Koenig et al. 2017): (34)
where δB_{s} is the difference between the coefficients B_{s} of the two spacecraft for SRP. Similar to atmospheric drag, F_{s} and f_{s} can be of the same order and their acceleration magnitudes are small. Therefore, we assumed that both F_{s} and f_{s} are compensated for by the control accelerations and neglected and the effects of f_{o} and f_{f}.
4.3 Control approach
The orbit and formation control approach used in this study is determined based on the characteristics of each perturbation. Table 3 summarizes the control approach used to mitigate each perturbation. The check mark in Table 3 indicates that all marked perturbations are compensated for by the control, and the check mark with an asterisk indicates that either the absolute or relative perturbations are compensated for. As summarized in Table 3, the CW nonlinearity must be nullified by the control to maintain the rigid formation; compensating for the relative perturbations for the Earth J_{2} gravity potential and lunisolar gravity is economical (except for the secular and longperiod effects on the eccentricity vector); the control strategy against the secular and longperiod changes on the eccentricity vector owing to the lunisolar and Earth J_{3} gravity potential are compared to decide whether to compensate for the absolute or relative perturbations; the absolute and relative accelerations for the atmospheric drag and SRP, both of which are small in most Earth orbits, are corrected by the control (although the compensation of the absolute accelerations may be unnecessary for certain applications).
We propose three possible methods for mitigating the secular and longperiod effects on the eccentricity vector, as shown in Fig. 1. The first approach is to maintain the eccentricity vector at zero by nullifying secular and/or longperiod drifts using continuous control. For a nearcircular orbit, the radial control force (u_{x}) is sensitive to the eccentricity vector and mean argument of latitude; the tangential force (u_{y}) is sensitive to the eccentricity vector and semimajor axis; and the normal force (u_{z}) is sensitive to the inclination and RAAN. Therefore, when u_{x} and u_{y} are used to compensate for the eccentricity vector drift, Gauss’s variational equations for the eccentricity vector are as follows: (35)
where and denote the change rates of the eccentricity vector by the control. Terms with e and e^{2} are neglected in these equations. Setting ux and uy as (36)
provides the necessary control to compensate for the eccentricity vector drift, where , and . In this context, and are considered as the time derivatives of the secular and longperiod variations on the eccentricity vector without control. This control law does not yield the secular and longperiod changes in a and ϕ. The second approach is to accept the small eccentricity and compensate for its relative perturbing acceleration using continuous control. As previously mentioned, the control acceleration can be calculated using Eqs. (10) and (17). The third approach combines the previous two methods: the eccentricity vector drift is left as long as it remains within an acceptable level, and the accumulated eccentricity is periodically nullified by impulsive or continuous control. A benefit of the third approach is that the propellant for the observation may be saved by using another propulsion system for the satellite bus system (if available) because the objective of this control is only to reconfigure the orbit of the chief, whereas the previous two methods must realize the formation maintenance for the observation and orbit control of the chief simultaneously; thus, the propellant for the observation cannot be saved.
For longperiod perturbations owing to J_{3} and lunisolar gravity, the eccentricity vector can remain around the frozen point. Thus, the first and second methods are compared. For secular perturbations due to lunisolar gravity, the eccentricity can increase linearly with time; thus, the secular drift of the eccentricity vector may not be acceptable. Thus, the first and third methods are compared.
Fig. 1 Three methods for mitigating the secular and longperiod eccentricity vector drift. The panels on the left ((a) edrift compensation) show the approach for compensating for the eccentricity vector drift of the chief using continuous control. The center panels ((b) eerror compensation) show the method that accepts the small eccentricity error and compensates for the relative perturbing acceleration of the deputy via continuous control. The panels on the right ((c) edrift and eerror compensation) show the hybrid approach. The top and bottom panels for each approach show the time histories of the eccentricity and control accelerations, respectively. 
Control approach for mitigating each perturbation in this study.
4.4 Orbit selection approach
By analyzing the factors that contribute the most to the control magnitudes, further insight is attained regarding the orbit selection approach. The developed perturbation models suggest that the semimajor axis (a) characterizes the absolute perturbations, whereas the semimajor axis and formation size (r_{r}) characterize the relative perturbations. Therefore, the orders of the control magnitudes for F_{c} and f_{p} can be evaluated by an explicit analysis of their accelerations as a function of a and/or r_{r} . Similarly, the control magnitudes for f_{o} and f_{f} can be evaluated using a two step analysis to first evaluate the factors for the perturbed orbital elements of the chief and then calculate the orders of their perturbing accelerations as a function of a and r_{r} . The three types of relative perturbations ( f_{p} , f_{o} , f_{f} , if available) have the same order against each perturbation source provided that ρ_{x}, ρ_{y}, and ρ_{z} are derived from the natural solution in Eq. (6).
Table 4 summarizes the factors that contribute the most to the control magnitudes for mitigating the absolute and relative perturbations, characterized by a, r_{r}, and other important parameters. The control magnitudes against the absolute perturbations relating to Earth’s J_{2} gravity, the longperiod drift of the eccentricity vector owing to Earth’s J3 gravity, and atmospheric drag decrease with an increase in a, whereas those of the luniso lar gravity increase, and that of the SRP remains constant. The control magnitudes against the relative perturbations for Earth’s J_{2} gravity, longperiod drift of the eccentricity vector owing to Earth’s J_{3} gravity, atmospheric drag, and SRP maintain the same trend as that for the absolute perturbations. Interestingly, those of the lunisolar gravity (except for the longterm drift of the eccentricity vector) remain constant, and their secular and long period effects on the eccentricity even decrease with an increase in a. Regarding the relationship with r_{r}, the control magnitudes against the relative perturbations increase linearly with r_{r} for Earth’s J_{2} gravity, the longperiod drift of the eccentricity vector owing to Earth’s J_{3} gravity, and lunisolar gravity, and are squared with r_{r} for the CW nonlinearity. In contrast, those for the atmospheric drag and SRP remain constant with respect to r_{r}.
To mitigate the longterm drift on the eccentricity vector, a more propellantsaving approach for Earth’s J_{3} gravity is the compensation for the absolute perturbations (with a factor of a^{−5}) in the larger semimajor axis and size of formation, and the relative perturbations (with a factor of a^{−4} r_{r}) for smaller a and r_{r} . However, a better approach for the lunisolar gravity is to compensate for the absolute perturbations in the smaller a and larger r_{r} as well as the relative perturbations in the larger a and smaller r_{r}.
Overall, smallperturbation regions tend to appear in the higheraltitude and shorterseparation regions. Thus, a reasonable approach for orbit selection is to search for the candidates of the semimajor axis that satisfy the smalldisturbance conditions under the desired formation size. However, the orbital orientation parameters (I and Ω) are selected to satisfy the observation conditions, considering the longterm perturbations. Because the observation conditions are highly missionoriented, the orbit selection approach is verified through case studies in the next section.
5 Results and discussion
Numerical examples of control acceleration magnitudes, observabilities, and visibilities are provided to verify the proposed approach to formationflying interferometry. The selected case studies were a laserinterferometric gravitationalwave telescope and an astronomical interferometer.
The analytical perturbation models developed in this study were used to analyze the control acceleration magnitudes and observation conditions. In addition, the control acceleration magnitudes were evaluated by averaging each perturbing acceleration per orbital revolution from the initial time and then calculating their rootsumsquare (Rss). Table 5 lists the common numerical conditions for all formation configurations. The areatomass ratio was assumed for a small class of satellites (e.g., M ≈ 500 kg and S ≈ 2.3 m^{2}), and its value was used to calculate atmospheric drag and SRP effects. To calculate the control acceleration by method (c) in Fig. 1, the secular drift of the eccentricity vector owing to lunisolar gravity was assumed to be periodically corrected to zero each year. The propellant consumption for this correction did not consider the control acceleration; the consumption to mitigate the relative acceleration due to the mean (nonzero) eccentricity vector accumulated over one year was considered only. The study did not account for propellant consumption related to formation reconfiguration maneuvers, such as altering observation directions or the size of the formation.
In the heliocentric orbit or halo orbit around the SunEarth L2 point, SRP is the dominant disturbance. Using the values in Table 5, the sum of the absolute and relative perturbing accelerations owing to SRP was calculated as approximately 4.5 × 10^{−8} m s^{−2} (1.4 m s^{−1} per year). Therefore, this study aims to find a small disturbance environment in geocentric orbits of less than 1.0 × 10^{−7} m s^{−2} (3.0 m s^{−1} per year) that maintains a magnitude only a few times larger than that in the heliocentric or SunEarth L2 orbits.
Largest contributing factors to the control magnitudes.
Typical numerical conditions for controlacceleration analysis.
Fig. 2 Reference formation orbit for the laserinterferometric gravitationalwave observatory. The distance between each spacecraft is given by (constant). 
5.1 Gravitationalwave telescope
The assumed laserinterferometric gravitationalwave telescope in this study is BDECIGO (Kawamura et al. 2021), which is a precursor to DECIGO. BDECIGO has three pairs of FabryPérot cavities, and one cluster of three satellites is placed equidistantly, 100 km apart, in a triangular formation.
The triangular reference formation for BDECIGO uses a general circular orbit (GCO; Alfriend et al. 2010), which is achieved by setting the parameters in Eq. (6) as ρ_{y} = 0, , and α_{z} = α_{x}. The GCO is always inclined by 60° with respect to the orbital plane of the chief, because where r_{r} = [x_{r}, y_{r}, z_{r}]^{T}. The reference formation can be configured by placing the three spacecraft in a GCO with orbit phase angles of 2π/3 from each other. The radius of the GCO (r_{r}) and the distance between each spacecraft (L_{12}, L_{23}, L_{31}) remain constant such that and . The control requirement for BDECIGO is imposed on the stability of each optical path length between the mirrors comprising the Fabry–Pérot cavity (e.g., ) within the order of nanometers. Based on the BDECIGO requirements, the sizing parameter was selected as . Spacecraft 1 (SC1), 2 (SC2), and 3 (SC3) were assigned α_{x} = 0, 2π/3, and 4π/3, respectively. Figure 2 illustrates the triangular formation of the GCO. For this reference formation, the first and second time derivatives of ρ_{x}, ρ_{y}, and ρ_{z} are zero, but is nonzero. Therefore, only the second term in Eq. (17) is nonzero, and the other terms are zero.
The most important observation requirement for the orbital orientation of BDECIGO is sunlight avoidance of the optical paths. The minimum avoidance angle was set to 15° in this study. This requirement can always be satisfied in a heliocentric orbit, whereas it can be violated in Earth orbits; therefore, a careful choice of orbital orientation parameters is necessary. To satisfy the sunlightavoidance requirement, the parameters were selected as in a high Earth orbit for the following reasons.
The target attitude frame is illustrated in Figs. 2 and 3 such that the two laser beams of each spacecraft point to [x_{t}, y_{t}, z_{t}]^{T} = [1, 0, 0]^{T} and and the z_{t}, axis is selected in the righthand frame. In addition, by approximating the inclination as the same as the obliquity of the ecliptic, the direction of the Sun is constrained within the xy plane in the LVLH frame. Then, the Sun vector in the target attitude frame (s^{t} = [s_{tx}, s_{ty}, s_{tz}]^{T}) is expressed as follows: (37)
and θ_{S} is the argument of latitude of the Sun, changing by 2π over one year. By defining the azimuth angle (ψ) and elevation angle (γ) with respect to the Sun vector as ψ = arctm2(s_{ty}, s_{tx}) and γ = sin^{−1} (s_{tz}), respectively (see Fig. 3), a pair of (ψ, γ) follows the locus similar to a “figure of eight” per orbit revolution of satellites passing the crosspoint (ψ, γ) = (Q + π/2, 0) twice, as shown in Fig. 4. In addition, angle Q drifts by 2π annually. Because the angle between the Sun vector and laser direction approaches zero when (ψ, γ) → (0, 0) or (π/3, 0) (see Fig. 3), each of the two laser axes in a single satellite is subject to interference by sunlight only once per year. Hence, excellent observability is expected. Furthermore, as the selected orbital plane approximately corresponds to the classical Laplace plane, longterm stability of the orbital orientation is expected, and this stability will maintain excellent observability.
To comprehensively understand the nature of geocentric orbit regions, the control acceleration magnitudes were computed over a wider range of spacecraft separations than the assumed value for BDECIGO. The left panel of Fig. 5 shows a contour map of the control acceleration magnitudes with respect to separation and altitude. The right panel shows the control acceleration magnitudes for each perturbing source at the desired separation for BDECIGO with respect to the semimajor axis.
The first feature is shown in Fig. 5 left, indicating that the smaller control acceleration magnitudes appeared in the shorterseparation and higheraltitude region. Therefore, a smalldisturbance region of less than 10^{−7} m s^{−2} was identified in high Earth orbit. The second feature is observed at the borderlines in Fig. 5 left to attain more propellantefficient control. In particular, the control approach to the secular drift of the eccentricity vector owing to lunar gravity switched from the mitigation of the absolute perturbing acceleration (in Fig. la with a factor of a^{2}) at an altitude lower than 51000 km to that of the relative perturbing acceleration (in Fig. 1c with a factor of a^{−1/2}rt). Consequently, the total control acceleration magnitudes were maintained as low for 10^{5} km along the semimajor axis, as shown in Fig. 5 right.
By selecting the candidate Earth orbit for BDECIGO at an altitude of 80000 km , the magnitude of the control acceleration was 7 × 10^{−8} m s^{−2}. This acceleration can be further broken down to compensate for the absolute and relative perturbations, which account for 17 % of the total acceleration (F_{c}), 20 % (f_{p}), 33 % (f_{o}), and 30 % (f_{f}). This result suggests that these four types of perturbed acceleration may contribute in the same order, and all of them must be considered in the control approach. The eccentricity (initially zero) in the selected orbital environment can increase to 0.0012 in a one year duration owing to the Moon’s gravity, based on the satellite motion model in Appendix C. Thus, the eccentricity error must be corrected periodically.
Then, the observation conditions were analyzed at an altitude of 80000 km. Figure 6 shows the rate of sunlight avoidance per orbit for one year. Each of the three laser cavities interfered twice per year, and the interference season lasted for approximately 40 days. The minimum satisfaction rate was 83% per orbit, which is approximately equivalent to 10 h of the 50 h (one orbital period) subject to interference by sunlight. One of the three cavities cannot be used during these periods. Nonetheless, sunlight interference might be avoided by taking appropriate measures, such as changing the target attitudes of the satellites while one of the three cavities suspends scientific observations.
It is important to note that an Earth eclipse can occur once per orbit. This is primarily because the orientation of the selected orbit is closely aligned with the ecliptic plane. As depicted in Fig. D.1, the duration of these eclipses can extend up to approximately 100 min during both the vernal and autumnal equinoxes. Such eclipse events may limit the observation time owing to potential thermal instability in the satellite or constraints on the electrical power system. While this factor should be considered in the design of satellites, a detailed exploration of its implications is beyond the scope of this study.
Fig. 3 Sun and laser directions in the target attitude frame. The azimuth and elevation angles (ψ, γ) of the Sun vector are also depicted. 
Fig. 4 Sunlightavoidance area and Sunvector locus per orbit with respect to the azimuth and elevation. The sunlightavoidance areas for the two laser axes of each spacecraft are depicted as gray circles. Equation (37) indicates that the Sunvector locus per orbit has azimuth and elevation amplitudes of nearly π/9 and π/3, respectively. The angle Q drifts by 2π over one year. 
Fig. 5 Control acceleration magnitudes for the laserinterferometric gravitationalwave observatory. The left panel shows the contour map of the control acceleration magnitudes (u_{d} = u_{d}) with respect to the separation (defined as r_{r}) and altitude. The gray lines show the borderlines that indicate when to change the mitigation methods for the secular and longperiod drifts of the eccentricity vector. For the secular and longperiod eccentricity drift owing to lunar gravity, the more propellantsaving approach appeared above the lines (specified by e^{(s)}(Moon) and e^{(lp)}(Moon)) when applying method (c) for the secular one and (b) for the longperiod one in Fig. 1, while it appeared below the lines when applying method (a) in Fig. 1. For the longperiod eccentricity drift owing to Earth’s J_{3} gravity, the more propellantsaving approach appeared above the line (specified by e^{(lp)}(J_{3})) when applying method (a) in Fig. 1, while it appeared below the line when applying method (b) in Fig. 1. The borderlines for the secular and longperiod drifts of the eccentricity vector owing to solar gravity were abbreviated because their magnitudes were sufficiently smaller than 10^{−9} m s^{−2} in the plotted region. The right panel shows the control acceleration magnitudes at the desired separation with respect to the semimajor axis. The atmospheric drag and secular and longperiod drift of the eccentricity vector owing to solar gravity were removed from the figure because their magnitudes were sufficiently smaller than 10^{−9} m s^{−2}. 
Fig. 6 Rate of sunlight avoidance per orbit during one year (January 1 – December 31, 2024). The six lines indicate the rates for the six lasertarget directions: from SC1 to SC2, from SC2 to SC3, from SC3 to SC1, and their opposite directions. Each laser direction is subject to interference by the Sun once per year; thus, each cavity is interfered twice per year. 
5.2 Astronomical interferometer
The astronomical interferometer assumed in this study uses a linear formation under orthogonal inertial pointing, similar to nanosatellite astronomical interferometers (Hansen & Ireland 2020; Matsuo et al. 2022). The target observation direction in the ECI frame is given by d^{I} = [cos δ cos α, cos δ sin α, sin δ]^{T}, where α and δ denote the right ascension and declination of the target observation direction, respectively. The reference formation can be configured by placing one spacecraft at the origin and the other two spacecraft symmetrically in the (y − z) plane in the LVLH frame (Hansen & Ireland 2020). A symmetrical spacecraft is placed by setting the parameters in Eq. (6) as ρ_{x} = 0, p_{z} = ρ_{y} tan p (−π/2 < p < π/2), and α_{z} = q (0 ≤ q ≤ 2π), and the other is placed symmetrically. Angles p and q are defined to satisfy the following relationships: (39)
where ∆Ω_{α} = (Ω − α). The distances between spacecraft 1 (origin) and spacecrafts 2 and 3 are not necessarily constant; however, their optical path difference must be sufficiently small to attain the interference signals from the target directions. Although specific missions were not assumed in this study, the selected sizing parameter was ρ_{y} = 0.25 km, similar to that in Hansen & Ireland (2020) as the typical mission. Figure 7 illustrates the linear formation under orthogonal inertial pointing in Earth’s orbit. For this reference, the first and second time derivatives of ρ_{x} and ρ_{y} are zero, but are nonzero; therefore, all four terms in Eq. (17) are nonzero. The first and second time derivatives of ρ_{z} are given by and , respectively, and ṗ and can be calculated from Eq. (39).
The most important observation condition for the orbital orientation of an astronomical interferometer is the coverage of the visible sky area. The key factors applied in this study are that the maximum angle between the observation target (star) and the Sun’s directions is less than π/3 (condition 1) to maintain the solar arrays of each spacecraft facing to the Sun and keep the sunlight from interfering with the optical paths of the science instruments; and the clearance to the target star occultation by the Earth and its atmospheric thickness during one orbit (condition 2). These factors are derived from the same ideas as those in previous research (Hansen & Ireland 2020); however, condition 2 is a stricter rule because temporal occultation during one orbit is allowed in Hansen & Ireland (2020), but not allowed in this study so that better observation conditions may be attained. In addition, this study imposes a further requirement that the angle p be less than π/3 to regulate the formation size up to (condition 3). The albedo effects of the Earth and Moon and the lunar occultation of the target stars were neglected in this study. Conditions 1–3 are given as follows: (40)
where s^{I} is the unit Sun vector in the ECI frame, and h_{t} is the Earth’s atmospheric thickness, assuming h_{t} = 100 km. A spherical Earth is assumed to derive condition 2. Finally, the secular drift of Ω owing to Earth’s J_{2} perturbation is actively used to increase the visible sky area during several years, so that the orbital inclination should be taken away from ±90° (see Eq. (B.2) for ). Considering these conditions, the orbital orientation parameters were selected as .
Lastly, the observation direction was set to α = 19^{h} 20^{m} and δ = 20° for the perturbing acceleration analysis. This selection was mainly aimed at understanding the typical acceleration magnitudes by taking the typical value of p as approximately 20^{°}, and was not intended at particular observations of astronomical objects.
The left panel of Fig. 8 shows a contour map of the control acceleration magnitudes with respect to the separation and altitude of the astronomical interferometer. The right panel shows the control acceleration magnitude at the desired separation (0.25 km) with respect to the semimajor axis. A trend similar to that of BDECIGO is shown in the left panel of Fig. 8: smaller control acceleration magnitudes appeared in the shorterseparation and higheraltitude region, and a smalldisturbance region of less than 10^{−7} m s^{−2} was identified in medium Earth orbit. The control approach to the longperiod drift of the eccentricity vector owing to Earth’s J_{3} gravity was to compensate for the relative perturbing acceleration (in Fig. 1b) for the lower altitude and shorter separation and the absolute perturbing acceleration (in Fig. 1a) for the higher altitude and longer separation.
By selecting the candidate Earth orbit for the astronomical interferometer at an altitude of 8000 km , the control acceleration magnitude was determined to be approximately 6 × 10^{−8} m s^{−2}. This acceleration can be further broken down to compensate for the absolute and relative perturbations, each accounting for 31% (F_{c}), 31% (f_{p}), 26% (f_{o}), and 12% (f_{f}).
Figure 9 shows the visible sky area for an astronomical interferometer at an altitude of 8000 km. The mission duration was assumed to be 5 years, where Ω drifted by 2π owing to the J_{2} perturbation. The visible area ranged for any possible α within approximately −75° < δ < 75° (depending on α), and the maximum visibility was approximately 25 %. In particular, good visibility appeared within the area of −30° < δ < 30°. Although Earth eclipses do occur in the selected orbit for the astronomical interferometer, they are not a constant phenomenon in each orbit. The maximum duration of these eclipses is approximately 40 min, as illustrated in Fig. D.l.
Fig. 7 Reference formation orbit for the astronomical interferometer. The two purple arrows are parallel and point in the same direction. 
Fig. 8 Control acceleration magnitudes for the linear astronomical interferometer. The left panel shows the contour map of the control acceleration magnitudes with respect to the separation (defined as ρ_{y}) and altitude. For the longperiod eccentricity drift owing to Earth’s J_{3} gravity, the more propellantsaving approach appeared above the line (specified by e^{(lp)} (J_{3})) when applying method (a) in Fig. 1, while it appeared below the line when applying method (b) in Fig. 1. The borderlines for the secular and longperiod drifts of the eccentricity vector owing to lunisolar gravity were abbreviated because their magnitudes were sufficiently smaller than 10^{−9} m s^{−2} in the plotted region. The right panel shows the control acceleration magnitudes at the desired separation (0.25 km) with respect to the semimajor axis. The atmospheric drag, lunisolar gravity, and secular and longperiod drifts of the eccentricity vector owing to the lunisolar gravity were removed from the figure because their magnitudes were sufficiently smaller than 10^{−9} m s^{−2}. 
Fig. 9 Visible sky area for the astronomical interferometer at an altitude of 8000 km and . The color map indicates the ratio of the visibility in all sky directions over 5 years. 
5.3 LEO environment
The use of LEOs has been considered in some formationflying interferometry missions, particularly for astronomical interferometers (Hansen & Ireland 2020; Matsuo et al. 2022). However, their disturbance characteristics tended to be analyzed within missionspecific baseline lengths and orbits. Therefore, it remains unclear where and to what extent a smalldisturbance environment can exist in LEO. Therefore, the LEO environment was assessed based on a perturbing acceleration analysis, assuming linear astronomical interferometers as potential applications. The Sun synchronous orbit was selected, where the RAAN drifts by 2π per year, and daylight is ensured in most seasons. The initial RAAN and inclination were set to and , respectively. The same observation direction in Sect. 5.2 was used for the perturbing acceleration analysis.
The left panel of Fig. 10 shows a contour map of the control acceleration magnitudes with respect to separation (up to 1 km) and altitude (up to 1000 km). The right panel shows the control acceleration magnitude at a specific separation (0.1 km) with respect to the semimajor axis. In contrast to the high and medium Earth orbits, the LEO has few smalldisturbance regions of less than 10^{−7} m s^{−2}. The largest perturbing acceleration sources were dominated by atmospheric drag (below approximately 600 km in altitude and 0.2 km in separation) and the Earth J_{2} gravity potential for the other regions in Fig. 10 left. However, a relatively small acceleration of 10^{−65} to 10^{−6}m s^{−2} (considering ∆V: 10–30 m s^{−1} per year) can be attained in LEO at an altitude > 500 km, with a separation of approximately 0.1 km. For example, the control acceleration magnitude at an altitude of 600 km was approximately 5 × 10^{−7} m s^{−2}. This acceleration can be further broken down to compensate for the absolute and relative perturbations, each accounting for 9% (F_{c}), 52% (f_{p},), 29% (f_{o}), and 10% (f_{f}). In the plotted region in Fig. 10 left, the better mitigation approach to the J_{3}perturbed longperiod drift of the eccentricity vector was to compensate for the relative perturbing acceleration (in Fig. 1b).
Figure 11 shows the visible sky area for an astronomical interferometer at an altitude of 600 km. The mission duration was assumed to be one year, where Ω drifted by 2π owing to the J_{2} perturbation. In comparison to the results in Fig. 9 and Hansen & Ireland (2020), the visible area reduced within approximately −14° < δ < 30° for any possible α, and the maximum visibility decreased to approximately 12 %. This reduction in visibility is attributed to condition 2 in Eq. (40). When condition 2 is omitted, the visibility returns to levels similar to those found in Fig. 9 and the work of Hansen & Ireland (2020). Earth eclipses also occur, but they are seasonally limited. Their maximum duration is approximately 20 min, as detailed in Fig. D.1.
A negative aspect of using LEO for formationflying interferometry is that the mission lifetime may be significantly shortened in comparison with those at higher geocentric orbits, owing to the limitation of the maximum load capability of the propellant. Another limitation is that precise formation control becomes more difficult with stronger perturbation effects. However, using LEO may outweigh these limitations from a benefittocost perspective, particularly for smallclass, experimental, and/or pathfinder missions. Moreover, using the latest lowthrust electrical propulsion with a high specific impulse can further improve the lifetime problem, as indicated in Hansen & Ireland (2020). The last point to discuss is safe autonomous formation flying in proximity. In particular, proximity formation flying at 0.1 km or less is considered highrisk for collisions (Monnier et al. 2019). In addition, the linear formation for the astronomical interferometer cannot apply the passively safe formation orbit known as the “eccentricity/inclination vector separation” technique (D’Amico & Montenbruck 2006). Therefore, active safety for collisionhazard detection and avoidance must be ensured using autonomous formationflying functions. However, if this technical bottleneck is resolved, the LEO environment will become an ideal platform for demonstrating various technologies for future space interferometry.
Fig. 10 Control acceleration magnitudes for the linear astronomical interferometer in LEO. The left panel shows the contour map of the control acceleration magnitudes with respect to the separation (defined as ρ_{y}) and altitude. The control magnitudes against the secular and longperiod drifts of the eccentricity vector owing to lunisolar gravity were sufficiently smaller than 10^{−9} m s^{−2} in the plotted region. The right panel shows the control acceleration magnitudes at the desired separation (0.1 km) with respect to the semimajor axis. The lunisolar gravity and the secular and longperiod drifts of the eccentricity vector owing to lunisolar gravity were removed from the figure because their magnitudes were sufficiently smaller than 10^{−9} m s^{−2}. 
Fig. 11 Visible sky area for the astronomical interferometer at an altitude of 600 km and . The color map indicates the ratio of the visibility in all sky directions over one year. 
6 Conclusions
This study demonstrates the feasibility of formationflying interferometry in geocentric orbits based on celestial mechanics. The perturbing accelerations to be mitigated by precise control can be categorized into four types: (1) absolute physical acceleration, (2) relative physical acceleration, (3) relative fictitious acceleration owing to the perturbed orbital motion of the chief, and (4) relative fictitious acceleration owing to the perturbed motion of the reference formation. In addition, the secular and/or longperiod perturbations on the eccentricity vector owing to Earth J_{3} and lunisolar gravity were unique because an efficient control approach could be changed to correct the absolute or relative perturbing accelerations. These results led to the characterization of the magnitude of each disturbance by the semimajor axis and the size of the formation. Consequently, smallperturbation conditions tended to appear in the higheraltitude and shorterseparation regions. Candidate orbits with a disturbance of less than 10^{−7} m s^{−2} were observed in the high Earth orbit for a laserinterferometric gravitationalwave telescope with a size of 100 km and in the medium Earth orbit for a linear astronomical interferometer with a size of 0.5 km. The LEO tended to have larger disturbance magnitudes; however, a LEO at a separation of approximately 0.1 km may be suitable for experimental purposes. These examples support the idea that the Earth’s orbit can be used for various types of formationflying interferometry.
Acknowledgements
This study was supported by JSPS KAKENHI Grant Number JP23K13501. We would like to thank Dr. Shuichi Sato of Hosei University and Dr. Taro Matsuo of Nagoya University for their valuable advice on the formation control requirements for the spacebased gravitationalwave observatory and astronomical interferometer.
Appendix A Nomenclature
Table A.1 lists the major symbols particularly used in the theoretical development of this study.
Nomenclature.
Appendix B J_{2}perturbed orbital elements
For a mean circular orbit, the J_{2}perturbed orbital motion of the chief up to the first order is given as follows: (B.1)
Appendix C Thirdbodyperturbed orbital elements
The disturbing function owing to the third body can be written as follows (Kozai 1973): (C.1)
where a_{b}, R_{b}, and n_{b} are the semimajor axis, radius, and mean motion of the disturbing body, respectively; P_{2}, P_{3}, and ⋯ are Legendre polynomials; cos S_{b} is the direction cosine between the radius to the spacecraft and that to the disturbing body; β = m_{b}/(m + m_{b}); and m and m_{b} are the masses of the Earth and disturbing body, respectively. The singleaveraged disturbing function 〈𝓡〉 was derived by obtaining the average with respect to the mean anomaly l of the chief. The doubleaveraged disturbing function 〈〈𝓡〉〉 is derived by averaging 〈𝓡〉 with respect to the mean anomaly l_{b} of the disturbing body. The longperiod and shortperiod terms of the disturbing function are calculated using 𝓡^{(lp)} = (〈𝓡〉 − 〈〈𝓡〉〉) and 𝓡^{(sp)} = (𝓡 − 〈𝓡〉), respectively. By replacing 𝓡 with 〈〈𝓡〉〉, 𝓡^{(lp)}, and 𝓡^{(sp)} in Lagrange’s equations, the orbital equations for the secular, longperiod, and shortperiod motions are obtained, and the orbital elements can be calculated analytically by integrating the orbital equations with respect to time.
If the terms with e^{2} and higher are neglected, the singleaveraged disturbing function up to P_{3}(cos S_{b}) of the Legendre polynomials is written as follows (Kozai 1973): (C.2)
The A and B are defined as follows: (C.3)
where α_{b} and δ_{b} denote the right ascension and declination of the disturbing body in the ECI frame, respectively. Assuming that the eccentricity of the third body (e_{b}) is small, the following approximations are applied: (a_{b}/R_{b})^{3} ≈ (1 + 3e_{b} cos l_{b}) and (a_{b}/R_{b})^{4} ≈ (1 + 4e_{b} cos l_{b}). This approximation is used to calculate 〈〈𝓡〉〉. By contrast, a further approximation of e_{b} = 0 (the circularorbit assumption for the third body) is used to calculate 𝓡^{(lp)}. The doubleaveraged disturbing function is calculated as follows: (C.4)
The constants K_{1} and K_{2} are introduced to obtain compact expressions: (C.6)
Subsequently, the orbital elements for the nearcircular satellite motion are provided. The secular parts are (C.7)
The longperiod parts are as follows: (C.9)
The shortperiod parts are as follows: (C.11)
To calculate the osculating orbital elements, the secular orbital elements are first computed. Then, the longperiod variation δ(·)^{(lp)} is calculated using the secular orbital elements, and the shortperiod variation δ(·)^{(sp)} is calculated using the secular and longperiod orbital elements.
Appendix D Eclipse effects of the Earth
Figure D.1 illustrates the Earth’s eclipse effects spanning from January 1, 2024, 00:00:00 UT to December 31, 2024, 24:00:00 UT for the orbits selected in Sect. 5. The three lines in Fig. D.1 represent the annual duration of eclipses per orbit for the orbits discussed in Sects. 5.1 (depicted in blue), 5.2 (in orange), and 5.3 (in yellow). These durations were computed under the assumption of circular orbits. For low and medium Earth orbits, (J_{2})perturbed secular motion was considered, whereas high Earth orbits were assumed to exhibit unperturbed motion. Furthermore, only the effects of the umbra were taken into account.
Fig. D.1 Eclipse effects of the Earth in various selected orbits. 
References
 Alfriend, K. T., Vadali, S. R., Gurfil, P., How, J., & Breger, L. S. 2010, Spacecraft Formation Flying: Dynamics, Control and Navigation, Elsevier Astrodynamics Series (Elsevier) [Google Scholar]
 Allan, R. R., & Cook, G. E. 1964, Proc. R. Soc. London Ser. A Math. Phys. Sci., 280, 97 [NASA ADS] [Google Scholar]
 Ashman, B., Bauer, F. H., Parker, J., & Donaldson, J. 2018, in 2018 SpaceOps Conference, https://arc.aiaa.org/doi/abs/10.2514/6.20182568 [Google Scholar]
 Brouwer, D. 1959, AJ, 64, 378 [NASA ADS] [CrossRef] [Google Scholar]
 Clohessy, W. H., & Wiltshire, R. S. 1960, J. Aerospace Sci., 27, 653 [CrossRef] [Google Scholar]
 Cockell, C. S., Herbst, T., & Léger, A. 2009, Exp. Astron., 23, 435 [NASA ADS] [CrossRef] [Google Scholar]
 Cook, G. E. 1962, Geophys. J. Int., 6, 271 [NASA ADS] [CrossRef] [Google Scholar]
 D’Amico, S., & Montenbruck, O. 2006, J. Guidance Control Dyn., 29, 554 [CrossRef] [Google Scholar]
 Gill, E., D’Amico, S., & Montenbruck, O. 2007, J. Spacecraft Rockets, 44, 671 [NASA ADS] [CrossRef] [Google Scholar]
 Hansen, J. T., & Ireland, M. J. 2020, PASA, 37, e019 [NASA ADS] [CrossRef] [Google Scholar]
 Harris, I., & Priester, W. 1962, J. Atmos. Sci., 19, 286 [NASA ADS] [CrossRef] [Google Scholar]
 Izzo, D., Sabatini, M., & Valente, C. 2003, in Proceedings of the 17th AIDAA National Congress, 1, 493 [Google Scholar]
 Kahle, R., Runge, H., Ardaens, J. S., et al. 2014, Acta Astron., 99, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Kawamura, S., Ando, M., Seto, N., et al. 2021, Progr. Theoret. Exp. Phys., 2021, 05A105 [CrossRef] [Google Scholar]
 Kawano, I., Mokuno, M., Kasai, T., & Suzuki, T. 2001, J. Spacecraft Rockets, 38, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Koenig, A. W., Guffanti, T., & D’Amico, S. 2017, J. Guidance Control Dyn., 40, 1749 [NASA ADS] [CrossRef] [Google Scholar]
 Kozai, Y. 1959, AJ, 64, 367 [NASA ADS] [CrossRef] [Google Scholar]
 Kozai, Y. 1973, A New Method to Compute Lunisolar Perturbations in Satellite Motions, Tech. rep., Smithsonian Astrophysical Observatory [Google Scholar]
 Lawson, P. R., & Dooley, J. A. 2005, Technology Plan for the Terrestrial Planet Finder Interferometer, Tech. rep., Jet Propulsion Laboratory [Google Scholar]
 Long, A., Cappellari, J., Velez, C., & Fuchs, A. 1989, Goddard Trajectory Determination System (GTDS) – Mathematical Theory – Revision 1, Tech. rep., Goddard Space Flight Center [Google Scholar]
 Matsuo, T., Ikari, S., Kondo, H., et al. 2022, J. Astron. Telesc. Instrum. Syst., 8, 015001 [NASA ADS] [CrossRef] [Google Scholar]
 Monnier, J., Aarnio, A., Absil, O., et al. 2019, A realistic roadmap to formation flying space interferometry, Tech. rep., Astro2020 APC White Paper [Google Scholar]
 Montenbruck, O., & Gill, E. 2000, Satellite Orbits: Models, Methods, and Applications (Springer) [CrossRef] [Google Scholar]
 Peñin, L. F., Scoarnec, Y., FernándezIbarz, J. M., et al. 2020, in 34th Annual Small Satellite Conference, Logan, Utah [Google Scholar]
 Quanz, S. P., Ottiger, M., Fontanet, E., et al. 2022, A&A, 664, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Roscoe, C., Vadali, S. R., & Alfriend, K. T. 2013, J. Astronaut. Sci., 60, 408 [NASA ADS] [CrossRef] [Google Scholar]
 Rosengren, A. J., Scheeres, D. J., & McMahon, J. W. 2014, Adv. Space Res., 53, 1219 [NASA ADS] [CrossRef] [Google Scholar]
 Tremaine, S., Touma, J., & Namouni, F. 2009, AJ, 137, 3706 [Google Scholar]
All Tables
Orders of the control acceleration magnitudes to compensate for the relative fictitious perturbations (f_{o} and f_{f}) for the secular, longperiod, and shortperiod terms, in comparison to those of the absolute physical perturbations (F_{c}).
All Figures
Fig. 1 Three methods for mitigating the secular and longperiod eccentricity vector drift. The panels on the left ((a) edrift compensation) show the approach for compensating for the eccentricity vector drift of the chief using continuous control. The center panels ((b) eerror compensation) show the method that accepts the small eccentricity error and compensates for the relative perturbing acceleration of the deputy via continuous control. The panels on the right ((c) edrift and eerror compensation) show the hybrid approach. The top and bottom panels for each approach show the time histories of the eccentricity and control accelerations, respectively. 

In the text 
Fig. 2 Reference formation orbit for the laserinterferometric gravitationalwave observatory. The distance between each spacecraft is given by (constant). 

In the text 
Fig. 3 Sun and laser directions in the target attitude frame. The azimuth and elevation angles (ψ, γ) of the Sun vector are also depicted. 

In the text 
Fig. 4 Sunlightavoidance area and Sunvector locus per orbit with respect to the azimuth and elevation. The sunlightavoidance areas for the two laser axes of each spacecraft are depicted as gray circles. Equation (37) indicates that the Sunvector locus per orbit has azimuth and elevation amplitudes of nearly π/9 and π/3, respectively. The angle Q drifts by 2π over one year. 

In the text 
Fig. 5 Control acceleration magnitudes for the laserinterferometric gravitationalwave observatory. The left panel shows the contour map of the control acceleration magnitudes (u_{d} = u_{d}) with respect to the separation (defined as r_{r}) and altitude. The gray lines show the borderlines that indicate when to change the mitigation methods for the secular and longperiod drifts of the eccentricity vector. For the secular and longperiod eccentricity drift owing to lunar gravity, the more propellantsaving approach appeared above the lines (specified by e^{(s)}(Moon) and e^{(lp)}(Moon)) when applying method (c) for the secular one and (b) for the longperiod one in Fig. 1, while it appeared below the lines when applying method (a) in Fig. 1. For the longperiod eccentricity drift owing to Earth’s J_{3} gravity, the more propellantsaving approach appeared above the line (specified by e^{(lp)}(J_{3})) when applying method (a) in Fig. 1, while it appeared below the line when applying method (b) in Fig. 1. The borderlines for the secular and longperiod drifts of the eccentricity vector owing to solar gravity were abbreviated because their magnitudes were sufficiently smaller than 10^{−9} m s^{−2} in the plotted region. The right panel shows the control acceleration magnitudes at the desired separation with respect to the semimajor axis. The atmospheric drag and secular and longperiod drift of the eccentricity vector owing to solar gravity were removed from the figure because their magnitudes were sufficiently smaller than 10^{−9} m s^{−2}. 

In the text 
Fig. 6 Rate of sunlight avoidance per orbit during one year (January 1 – December 31, 2024). The six lines indicate the rates for the six lasertarget directions: from SC1 to SC2, from SC2 to SC3, from SC3 to SC1, and their opposite directions. Each laser direction is subject to interference by the Sun once per year; thus, each cavity is interfered twice per year. 

In the text 
Fig. 7 Reference formation orbit for the astronomical interferometer. The two purple arrows are parallel and point in the same direction. 

In the text 
Fig. 8 Control acceleration magnitudes for the linear astronomical interferometer. The left panel shows the contour map of the control acceleration magnitudes with respect to the separation (defined as ρ_{y}) and altitude. For the longperiod eccentricity drift owing to Earth’s J_{3} gravity, the more propellantsaving approach appeared above the line (specified by e^{(lp)} (J_{3})) when applying method (a) in Fig. 1, while it appeared below the line when applying method (b) in Fig. 1. The borderlines for the secular and longperiod drifts of the eccentricity vector owing to lunisolar gravity were abbreviated because their magnitudes were sufficiently smaller than 10^{−9} m s^{−2} in the plotted region. The right panel shows the control acceleration magnitudes at the desired separation (0.25 km) with respect to the semimajor axis. The atmospheric drag, lunisolar gravity, and secular and longperiod drifts of the eccentricity vector owing to the lunisolar gravity were removed from the figure because their magnitudes were sufficiently smaller than 10^{−9} m s^{−2}. 

In the text 
Fig. 9 Visible sky area for the astronomical interferometer at an altitude of 8000 km and . The color map indicates the ratio of the visibility in all sky directions over 5 years. 

In the text 
Fig. 10 Control acceleration magnitudes for the linear astronomical interferometer in LEO. The left panel shows the contour map of the control acceleration magnitudes with respect to the separation (defined as ρ_{y}) and altitude. The control magnitudes against the secular and longperiod drifts of the eccentricity vector owing to lunisolar gravity were sufficiently smaller than 10^{−9} m s^{−2} in the plotted region. The right panel shows the control acceleration magnitudes at the desired separation (0.1 km) with respect to the semimajor axis. The lunisolar gravity and the secular and longperiod drifts of the eccentricity vector owing to lunisolar gravity were removed from the figure because their magnitudes were sufficiently smaller than 10^{−9} m s^{−2}. 

In the text 
Fig. 11 Visible sky area for the astronomical interferometer at an altitude of 600 km and . The color map indicates the ratio of the visibility in all sky directions over one year. 

In the text 
Fig. D.1 Eclipse effects of the Earth in various selected orbits. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.