Formation-flying interferometry in geocentric orbits

Spacecraft formation flying serves as a method of astronomical instrumentation that enables the construction of large virtual structures in space. The formation-flying interferometry generally requires very-high control accuracy, and beyond-Earth orbits are typically selected. By contrast, this study proposes the use of geocentric orbits for formation-flying interferometry. A geocentric orbit is beneficial because of its economic accessibility and the availability of flight-proven technologies for formation-flying autonomy, safety, and management. Its feasibility depends on the existence of specific orbits that satisfy a small-disturbance environment with favorable observation conditions. This theory, developed based on celestial mechanics, indicates that small-perturbation regions tend to appear in higher-altitude and shorter-separation regions. Candidate orbits are identified in high Earth orbit for the triangular laser-interferometric gravitational-wave telescope, which is 100 km in size, and in medium Earth orbit 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. As shown in these examples, geocentric orbits are potentially applicable for various types of formation-flying interferometry.


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 space-based 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 Deci-hertz Interferometer Gravitational Wave Observatory (DECIGO; Kawamura et al. 2021) is a laser-interferometric, space-based gravitational-wave 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 near-circular heliocentric orbit.Fabry-Pérot cavities are employed in DECIGO, each comprising two mirrors to measure the arm-length changes.The Large Interferometer For Exoplanets (LIFE; Quanz et al. 2022) is a space-based 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 Sun-Earth Lagrange second (L2) point.The current design for LIFE assumes four collector spacecrafts and one beam-combiner 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.
Flight-proven 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 Satellite-VII (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 formation-flying 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 TerraSAR-X/TanDEM-X (TerraSAR-X 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 On-Board Autonomy-3 (PROBA-3) will demonstrate several formation activities, including coarse-to-fine control, autonomy, and safety management, to observe the solar coronagraph in a highly elliptical orbit.Its formation-flying precision is capable of submillimeter relative displacement and arcsecond-order pointing accuracy.This accuracy will be achieved using state-of-the-art Ito, T.: A&A, 682, A38 (2024) metrology sensors and millinewton-class thrusters (Peñin 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 formation-flying interferometry.
System-level, end-to-end, in-space verification for the pathfinder may be necessary to mature formation-flying interferometry before full-scale 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 small-disturbance 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 formation-flying autonomy, safety, and management.Such drawbacks may lead to formationflying interferometry becoming a larger and longer-term 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 flight-proven technologies for formation flying autonomy, safety, and management with GNSS.Furthermore, state-of-the-art 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 formation-flying interferometry.Representative and state-of-the-art research in this context conducted by Hansen & Ireland (2020) involved a linear astronomical interferometer in a near-circular 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 Earth-orbit 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 formation-flying interferometry as site characterization is for the deployment of ground-based telescopes.
This study investigates the feasibility of formation-flying 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 propellant-efficient control strategy for correcting the absolute and/or relative perturbing accelerations.Secular and long-period 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 formation-flying interferometry.The proposed approach as well as methods for satisfying the observation conditions will be useful for future applications.
This study considers only a near-circular Earth orbit because such is sufficient for various potential applications in formationflying interferometry.In terms of terminology, this study uses the chief-deputy expression to describe formation-flying satellites.A chief does not necessarily represent the actual chief satellite, but rather the virtual origin of formation-flying 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 relative-motion model for a near-circular 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 formation-flying interferometry.The concluding remarks are presented in Sect.6.

Unperturbed relative motion
The relative equation of motion of the deputy with respect to the chief in the Earth-centered inertial (ECI) frame is where µ e is the Earth gravitational constant, r I is the relative position vector of the deputy (r By assuming the unperturbed motion of both the chief and deputy ( f I p = 0) and the uncontrolled motion of the chief (u I c = 0), Eq. (1) in the ECI frame is transformed into the local-verticallocal-horizontal (LVLH) frame: 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 right-hand frame completes the setup.Assuming r ≪ R c (r = ||r||), the first term on the right-hand side of Eq. ( 2) can be linearized as follows: For a circular orbit, θ = n, θ = 0, 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 n = µ e /a 3 , the unperturbed and linear relative dynamics model for a circular orbit is obtained as follows: where A u1 and A u2 are matrices defined as follows: (5) Equation ( 4) is known as the Clohessy-Wiltshire (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 A u2 ṙ with the control accelerations.However, the control acceleration requirement tends to become large, which is inefficient and unrealistic for formation-flying 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 along-track drift) is 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 θ = n for the unperturbed circular orbit; therefore, r r is also a function of time.

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 un-modeled 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): where ω = [ω x , ω y , ω z ] T is the orbital angular velocity vector of the chief in the LVLH frame, f g2B is the linearized gravity gradient acceleration owing to the two-body gravitational field ( f g2B = −(µ e /R 3 c )[−2x, y, z] T ), f p 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) 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 İ and Ω from Gauss's variational equations into Eq.( 8) (Alfriend et al. 2010).
Next, Eq. ( 7) is transformed into 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 user-defined target angular velocity.In Eq. ( 9), f o takes the following form: 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 near-circular-orbit assumption, the argument of latitude can be approximated as follows: 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 φ = (n + ∆ φ).Moreover, the time derivatives ∆ φ, İ, and Ω are assumed sufficiently small with respect to n.By neglecting the second-and higherorder terms of the small variables, (A p1 − A t1 ) and (A p2 − A t2 ) in Eq. ( 10) can be expressed as follows: Relative fictitious perturbation owing to the perturbed motion of the reference formation. where The following relationships are used to obtain Eqs. ( 12) and ( 13): Then, the desired reference formation is selected as the natural solution of Eq. ( 6): where Θ x and Θ z are the new variables satisfying Θ x = (θ + α x ), Θ z = (θ + α z ), Θx = Θz = ω t , and Θx = Θz = ωt , and the parameters of (ρ x , ρ y , ρ z ) can be time-variant.The relationships between r r , ṙr , and rr can be written as follows: where 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: where ϵ = (r − r r ) is the relative positional error with respect to the reference formation.By setting the control command of the deputy to 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 (r = r r , ṙ = ṙr ), 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 F I c 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 control-effort trade-off 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.

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, long-period, and short-period 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.
A38, page 4 of 18 Ito, T.: A&A, 682, A38 (2024) Table 2. Orders of the control acceleration magnitudes to compensate for the relative fictitious perturbations ( f o and f f ) for the secular, long-period, and short-period terms, in comparison to those of the absolute physical perturbations (F c ).

Terms Secular
Long-period Short-period ( Notes.In the first column, (•) denotes the osculating orbital elements comprising the secular, long-period, and short-period terms as expressed in Eq. ( 20).The osculating orbital elements of (e x , e y ) are substituted to (•)n 2 r r , and those of (e x , e y , I, Ω, ϕ) are substituted to ( •)nr r and ( •)r r for the analysis.In the second to fourth columns, c (s) , c (lp) , and c (sp) indicate the orders of the control acceleration magnitudes to compensate for the absolute physical perturbations for the secular, long-period, and short-period terms, respectively.They are defined as 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 ( •)nr r and ( •)r r 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: ( •)na for (e x , e y , I, Ω, ϕ).At this point, the osculating orbital elements are modeled as follows: where (•) (s) is the secular term that is a linear function of time, δ(•) (lp) is the long-period term for which the angular frequency has an order of n l (≪ n), δ(•) (sp) is the short-period term for which the angular frequency has an order of n, (•) (s) 0 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, long-period, and short-period 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, long-period, 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 long-period 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 long-period 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.

Analytical perturbation models
The perturbation sources considered in this study are CW nonlinearity, Earth J 2 gravity potential, Earth J 3 gravity potential, lunisolar gravity, atmospheric drag, and SRP.Their analytical models are explained in the following.

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 right-hand side of Eq. ( 3) from the left-hand side.The CW nonlinearity can contribute to a part of f p .

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 low-to-medium Earth orbit.By considering only zonal effects, the J 2 -perturbed absolute acceleration vector acting on the chief (F J 2 , as a part of F c ) is given by (Alfriend et al. 2010) where R e is the mean equatorial radius of Earth.The linearized J 2 differential acceleration vector ( f J 2 , as a part of f p ) has the following form (Izzo et al. 2003): where A comparison of Eqs. ( 21) and ( 22) shows that || f J 2 || is always less than ||F J 2 || 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 A38, page 5 of 18 Ito, T.: A&A, 682, A38 (2024) elements owing to the J 2 potential for the near-circular 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.

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 long-period 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 long-period J 3 perturbation are as follows (Kozai 1959): where ω is the J 2 -perturbed secular drift rate of the argument of perigee given by ω = (3/4)J 2 (R e /a) 2 (5 cos 2 I − 1)n for the near-circular 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): 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 near-circular 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 .

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) where R I b is the position vector of the third body in the ECI frame, µ b denotes the gravity constant of the third body, C L I is the rotation matrix from the ECI frame to the rotational (LVLH) frame of the chief, and R b = ||R I b ||.The linearized gravity gradient acceleration vector owing to the third body ( f 3rd , as a part of f p ) is 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, long-periodic, and short-period motions in a near-circular 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 long-period satellite motion due to lunisolar perturbations up to the second degree of Legendre polynomials, but short-period motions were not analyzed.Subsequent work (Kozai 1973) analytically proposed the calculation of short-period satellite motion caused by lunisolar perturbations, but secular and long-period motions were calculated numerically.Recent work (Roscoe et al. 2013) developed secular and periodic motion models owing to third-body perturbations using a single-and double-averaged disturbing function.Analytical models were developed based on the assumption of a circular-restricted three-body problem; however, this may be problematic, particularly when long-term 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, near-circular satellite motion model based on to lunisolar gravity, and the results are given in Appendix C.
In addition, the long-term effects of the lunisolar gravitational field and J 2 term of the Earth's oblateness may be comparable for a medium-to-high 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) where ϵ e is the obliquity of the ecliptic (ϵ e = 23.4• ), Φ is the azimuth angle, and ω J 2 , ω m , and ω s are defined as 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 gravitational-wave detector, for obtaining the long-term stability on the orbital orientation with respect to the Sun, as discussed in Sect.5.1.

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 atmospheric-drag acceleration vector acting on the chief spacecraft in the LVLH frame (F D , as a part of F c ) is as follows: 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 cross-sectional area, C D is the drag coefficient, and m is the spacecraft mass.Equation ( 30) is obtained by assuming that the atmosphere co-rotates with the Earth (Montenbruck & Gill 2000) and V is approximately one order of magnitude higher than that of the Earth-fixed 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 Harris-Priester density model (Harris & Priester 1962) is relatively simple and is widely used as a standard atmospheric model.The model is given as follows: 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 lowinclination 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 atmospheric-density 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: 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 .

Solar radiation pressure
Solar radiation pressure (SRP) can be regarded as a function of the Sun-to-spacecraft 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): 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): 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 .

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 long-period effects on the eccentricity vector); the control strategy against the secular and long-period 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 long-period 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 long-period drifts using continuous control.For a near-circular 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 A38, page 7 of 18 Ito, T.: A&A, 682, A38 (2024) Fig. 1.Three methods for mitigating the secular and long-period eccentricity vector drift.The panels on the left ((a) e-drift compensation) show the approach for compensating for the eccentricity vector drift of the chief using continuous control.The center panels ((b) e-error 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) e-drift and e-error compensation) show the hybrid approach.The top and bottom panels for each approach show the time histories of the eccentricity and control accelerations, respectively.A check mark indicates that all marked perturbations are compensated for by the control, and a check mark with an asterisk indicates that either the absolute or relative perturbations are compensated for.The perturbation specified as "Lunisolar gravity" excludes its secular and long-period effects on the eccentricity vector drift.The e (s) and e (lp)  indicate the secular and long-period perturbations on the eccentricity vector, respectively.
variational equations for the eccentricity vector are as follows: where ėx and ėy denote the change rates of the eccentricity vector by the control.Terms with e and e 2 are neglected in these equations.Setting u x and u y as provides the necessary control to compensate for the eccentricity vector drift, where |ė| = ė2 x + ė2 y , cos ϕ 0 = −ė x /|ė|, and sin ϕ 0 = ėy /|ė|.In this context, ėx and ėy are considered as the time derivatives of the secular and long-period variations on the eccentricity vector without control.This control law does not yield the secular and long-period 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 long-period 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.

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 twostep analysis to first evaluate the factors for the perturbed orbital A38, page 8 of 18 Ito, T.: A&A, 682, A38 (2024) 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 long-period drift of the eccentricity vector owing to Earth's J 3 gravity, and atmospheric drag decrease with an increase in a, whereas those of the lunisolar gravity increase, and that of the SRP remains constant.The control magnitudes against the relative perturbations for Earth's J 2 gravity, long-period 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 long-term drift of the eccentricity vector) remain constant, and their secular and longperiod 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 long-period 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 long-term drift on the eccentricity vector, a more propellant-saving 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, small-perturbation regions tend to appear in the higher-altitude and shorter-separation regions.Thus, a reasonable approach for orbit selection is to search for the candidates of the semimajor axis that satisfy the small-disturbance conditions under the desired formation size.However, the orbital orientation parameters (I and Ω) are selected to satisfy the observation conditions, considering the long-term perturbations.Because the observation conditions are highly mission-oriented, the orbit selection approach is verified through case studies in the next section.

Results and discussion
Numerical examples of control acceleration magnitudes, observabilities, and visibilities are provided to verify the proposed approach to formation-flying interferometry.The selected case studies were a laser-interferometric gravitational-wave 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 root-sum-square (Rss).Table 5 lists the common numerical conditions for all formation configurations.The area-to-mass 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 Earth J 2 gravity a −4 a −5 r r Earth J 3 gravity (e (lp) ) a −5 a −4 r r Lunisolar gravity a r r Lunisolar gravity (e (s) ) a 2 a −1/2 r r t Lunisolar gravity (e (lp) ) a 2 a −1/2 r r Atmospheric drag ρa −1 ρa −1 Solar radiation pressure Const.Const.
Notes.For Earth J 3 gravity (e (lp) ) and lunisolar gravity (e (s) and e (lp) ), the eccentricity originally had a contributing factor of (en 2 r r ) ∝ (ea −3 r r ) to the relative fictitious accelerations, as shown in Sect.4.1.Subsequently, the contributing factors were calculated by substituting the secular and long-period parts of the eccentricity for Earth J 3 and lunisolar gravity into (ea −3 r r ).Notes.The initial orbital elements of a (s) 0 , I (s) 0 , and Ω (s) 0 were selected based on each mission application, as given in  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 Sun-Earth 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 Sun-Earth L2 orbits.

Gravitational-wave telescope
The assumed laser-interferometric gravitational-wave telescope in this study is B-DECIGO (Kawamura et al. 2021), which is a precursor to DECIGO.B-DECIGO has three pairs of Fabry-Pérot cavities, and one cluster of three satellites is placed equidistantly, 100 km apart, in a triangular formation.
The triangular reference formation for B-DECIGO uses a general circular orbit (GCO; Alfriend et al. 2010), which is achieved by setting the parameters in Eq. ( 6) as ρ y = 0, ρ z = √ 3ρ x , and α z = α x .The GCO is always inclined by 60 • with respect to the orbital plane of the chief, because tan −1 (z r / x 2 r + y 2 r ) = π/3, 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 r r = 2ρ x and L 12 = L 23 = L 31 = 2 √ 3ρ x .The control requirement for B-DECIGO is imposed on the stability of each optical path length between the mirrors comprising the Fabry-Pérot cavity (e.g., L12 (t)dt) within the order of nanometers.Based on the B-DECIGO requirements, the sizing parameter was selected as r r = 100/ √ 3 km.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 ωt 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 B-DECIGO 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 sunlight-avoidance requirement, the parameters were selected as I (s) 0 = 20 • and Ω (s) 0 = 0 • 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 [1/2, √ 3/2, 0] T and the z t axis is selected in the right-hand frame.In addition, by approximating the inclination as the same as the obliquity of the ecliptic, the direction 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. Sunlight-avoidance area and Sun-vector locus per orbit with respect to the azimuth and elevation.The sunlight-avoidance areas for the two laser axes of each spacecraft are depicted as gray circles.Equation (37) indicates that the Sun-vector locus per orbit has azimuth and elevation amplitudes of nearly π/9 and π/3, respectively.The angle Q drifts by 2π over one year.
of the Sun is constrained within the x-y 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: where 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 ψ = arctan2(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 cross-point (ψ, γ) = (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 A38, page 10 of 18 Ito, T.: A&A, 682, A38 (2024) Fig. 5.Control acceleration magnitudes for the laser-interferometric gravitational-wave 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 long-period drifts of the eccentricity vector.For the secular and long-period eccentricity drift owing to lunar gravity, the more propellant-saving 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 long-period one in Fig. 1, while it appeared below the lines when applying method (a) in Fig. 1.For the long-period eccentricity drift owing to Earth's J 3 gravity, the more propellant-saving 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 long-period 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 (100/ √ 3 km) with respect to the semimajor axis.The atmospheric drag and secular and long-period 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.
plane approximately corresponds to the classical Laplace plane, long-term 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 B-DECIGO.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 (100/ √ 3 km) for B-DECIGO 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 shorter-separation and higher-altitude 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 propellant-efficient 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. 1a with a factor of a 2 ) at an altitude lower than 51 000 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 B-DECIGO at an altitude of 80 000 km (a (s)  0 = 86 378 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 80 000 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 A38, page 11 of 18 Ito, T.: A&A, 682, A38 (2024) 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.

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, ρ 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: cos p = sin ∆Ω α sin I cos δ + cos I sin δ, sin p sin q = sin ∆Ω α cos I cos δ − sin I sin δ, sin p cos q = cos ∆Ω α cos δ, (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 ρz , ρz , and ωt are nonzero; therefore, all four terms in Eq. ( 17) are nonzero.The first and second time derivatives of ρ z are given by ρz = ( ṗρ y / cos 2 p) and ρz = [( p + 2 ṗ2 tan p)ρ y / cos 2 p], respectively, and ṗ and p 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 2ρ y (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: 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 I (s) 0 = 70 • and Ω (s) 0 = 0 • .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 B-DECIGO is shown in the left panel of Fig. 8: smaller control acceleration magnitudes appeared in the shorterseparation and higher-altitude region, and a small-disturbance region of less than 10 −7 m s −2 was identified in medium Earth orbit.The control approach to the long-period 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.
A38, page 12 of 18 Ito, T.: A&A, 682, A38 (2024) 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 long-period eccentricity drift owing to Earth's J 3 gravity, the more propellant-saving 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 long-period 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 .By selecting the candidate Earth orbit for the astronomical interferometer at an altitude of 8000 km (a (s) 0 = 14 378 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.1.

LEO environment
The use of LEOs has been considered in some formation-flying interferometry missions, particularly for astronomical interferometers (Hansen & Ireland 2020;Matsuo et al. 2022).However, their disturbance characteristics tended to be analyzed within mission-specific baseline lengths and orbits.Therefore, it remains unclear where and to what extent a small-disturbance 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 Ω (s) 0 = 10 • and I (s) 0 = 97.8• , 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 small-disturbance 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 −6.5 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 (a (s) 0 = 6 978 km) was approximately 5 × 10 −7 m s −2 .This acceleration can be further broken down to compensate for A38, page 13 of 18 Ito, T.: A&A, 682, A38 (2024) 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 long-period 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 long-period 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 plotted region in Fig. 10 left, the better mitigation approach to the J 3 -perturbed long-period 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 A negative aspect of using LEO for formation-flying 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 benefitto-cost perspective, particularly for small-class, experimental, and/or pathfinder missions.Moreover, using the latest low-thrust 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 high-risk 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 collision-hazard detection and avoidance must be ensured using autonomous formation-flying functions.However, if this technical bottleneck is resolved, the LEO environment will become an ideal platform for demonstrating various technologies for future space interferometry.

Conclusions
This study demonstrates the feasibility of formation-flying 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 A38, page 14 of 18 Ito, T.: A&A, 682, A38 (2024) and the size of the formation.Consequently, small-perturbation conditions tended to appear in the higher-altitude 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 laser-interferometric gravitational-wave 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 formation-flying interferometry.e x = J 16 14 sin 2 I (s) 0 cos 3ϕ (s) + 3 3 + 5 cos 2I (s) 0 cos ϕ (s) , e y = J 16 14 sin 2 I (s) 0 sin 3ϕ (s) + 3 1 + 7 cos 2I (s) 0 sin ϕ (s) , I = I (s) 0 + 3J 8 sin 2I (s) 0 cos 2ϕ (s) , Ω = Ω (s) 0 + Ωt +

Appendix C: Third-body-perturbed orbital elements
The disturbing function owing to the third body can be written as follows (Kozai 1973):  (lp) , and R (sp) in Lagrange's equations, the orbital equations for the secular, long-period, 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): R I c and R I d are the position vectors of the chief and deputy, respectively, f I p is the disturbance acceleration vector of the deputy relative to the chief ( f I p = F I d − F I c ), F I c and F I d are the disturbance acceleration vectors for the chief and deputy, respectively, u I c and u I d are the control acceleration vectors applied to the chief and deputy, respectively, R c = ||R I c ||, and R d = ||R I d ||.The superscript I indicates the position, velocity, or acceleration in the ECI frame.

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

Fig. 9 .
Fig. 9. Visible sky area for the astronomical interferometer at an altitude of 8000 km and I (s) 0 = 70 • .The color map indicates the ratio of the visibility in all sky directions over 5 years.

Fig. 11 .
Fig. 11.Visible sky area for the astronomical interferometer at an altitude of 600 km and I (s) 0 = 97.8• .The color map indicates the ratio of the visibility in all sky directions over one year.
cos S b ) + • • • , (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 single-averaged disturbing function ⟨R⟩ was derived by obtaining the average with respect to the mean anomaly l of the chief.The double-averaged disturbing function ⟨⟨R⟩⟩ is derived by averaging ⟨R⟩ with respect to the mean anomaly l b of the disturbing body.The long-period and short-period terms of the disturbing function are calculated using R (lp) = (⟨R⟩ − ⟨⟨R⟩⟩) and R (sp) = (R − ⟨R⟩), respectively.By replacing R with ⟨⟨R⟩⟩, R

Table 1 .
Categories of perturbing accelerations that can be mitigated by the control.

Table 3 .
Control approach for mitigating each perturbation in this study.

Table 4 .
Largest contributing factors to the control magnitudes.

Table 5 .
Typical numerical conditions for control-acceleration analysis.