Autonomous formation flight using solar radiation pressure

Autonomous formation flight enables new satellite missions for novel applications. The cost and limits of propulsion systems can be overcome if environmental resources are being benefitted of. Currently, atmospheric drag is used in low Earth orbit to this end. Solar radiation pressure, which is of similar order of magnitude as aerodynamic ram pressure, is, however, always neglected. We introduce this force and show that it can be exploited. We demonstrate through simulations that a formation geometry is established quicker if the solar radiation pressure is modeled.


Formation flight
The recent advent of small, low-cost and reduced-performance satellites, such as CubeSats, allows new, exciting mission concepts. Employing several of these satellites and combining their capacities promises entirely novel applications. Most prominently, Planet [1] and Spire [2] have launched constellations, i.e., hundred or more, small satellites for new, unparalleled commercial uses which were economically unachievable with larger satellites of the past. Others satellite swarms serve scientific purposes [3]. Recently, concepts of mega-constellations, i.e., constellations consisting of thousands or tens of thousands of satellites, have emerged to fill and develop new niche uses, such as for Internet of Things (IoT) applications [4]. Industry is ramping up production capacity to meet the expected high demand.
In recent years, the use of small satellites, among them CubeSats, for formation flight-that is, a small number of satellites flying in proximity-has received more attention. While theoretical foundations were established decades ago [5]- [7], formation flight has been further developed lately with the introduction of control laws [8]- [12] enabling autonomous operations. The source of the control force is typically a reaction control propulsion system. Only a few missions have been realized so far, such as Prisma [8], CanX-4/5 [13] and Hawkeye [14]. Very recently, the Net-Sats [15] formation flight mission has been launched; Yoon [16] demonstrates the use of drag for limited relative orbit control, i.e., the keeping distance between different satellites for collision avoidance. All these missions have been primarily experimental in nature.

Solar-aerodynamic controlled formations
A new method to generate control forces for formation flight is the use of aerodynamic forces such as drag and lift. While Planet Inc. [1] and others already use drag for distance maintenance, Ivanov [9] and Traub [17] showed that lift can also be used; lift is particularly well suited for formation flight as it enables out-of-orbital-plane forces.
The benefits of aerodynamic forces for orbital control are numerous. Avoiding the need for a propulsion system and its constituent fuel tanks reduces directly the cost of the space system. A secondary effect of this benefit is that the spacecraft size can shrink, which further reduces cost. Third, the amount of propellant limits the lifetime of classical satellite operations. Using environmental resources such as the atmosphere and solar radiation pressure prevents this constraint.
However, aerodynamic forces are small. In particular lift, which is a primary force for out-of-orbital-plane maneuvers, 1 3 is even very small [17]. Yet, these forces are available in limitless supply and, therefore, for unlimited duration. The integral of the small forces over time is appreciable and, therefore, exploitable, as we will show.
Until today, solar radiation pressure-which can create forces of similar magnitude to aerodynamic forces-has not been addressed for flight in low earth orbit (LEO), but only for higher altitudes such as geostationary orbits [18] or at the deep space Lagrange points [19]. This is surprising since it is known that they are of similar magnitude to aerodynamic forces in LEO [20]. Accounting for solar radiation pressure in a formation flight control algorithm is the main novelty of our research.
Solar radiation pressure-generated force can increase the availability of control forces in general and is in particular suitable for dusk/dawn sun-synchronous orbits where the solar light vector is near-perpendicular to the orbital plane and, therefore, efficient for out-of-plane maneuvers for which lift provides only a small force. Such a situation is shown in Fig. 1. Both forces act upon the surfaces of satellite and can be used for orbital and, hence, formation control. Dusk/dawn orbits experience only short durations of eclipse at the winter solstice. Solar radiation pressure is, therefore, almost permanently available to provide significant benefit to formation geometry control.
Besides being a resource for orbital control, the correct modeling of solar radiation pressure is also important for applications where satellite position needs to be maintained with high accuracy. This is the case, for instance, for formations with satellite-satellite VHF beam-forming techniques requiring close proximity of satellites [21].
Since the solar radiation pressure force is of appreciable amount, it also promises benefits for controlling formations around celestial bodies such as the Earth's moon, which is, however, not addressed here.

Orbital environment
CubeSats are exclusively used in LEO due to conceptual and technological constraints such as the utilization of Commercial-Off-The-Shelf (COTS) components, power needed for space-to-Earth communication, availability of a planetary magnetic field for basic attitude control or reliable and free availability of navigation information such as two-lineelements (TLE). Notable exceptions are the MarCO satellites [22], which performed a Mars fly-by. However, these exceptional nanosatellites are not typical CubeSats: neither in terms of technology nor financial budget.
CubeSat are usually found on circular Earth orbits for which the orbital speed v and the orbital period T, which depend solely on the radius r of the orbit, are defined as the following: The symbol, μ, is the standard gravitational parameter of Earth.

Residual atmosphere
In LEO, a highly rarefied yet appreciable atmosphere exists, which is called thermosphere due to its relatively high gas temperature. The high temperature of this upper atmosphere layer is the consequence of the absorption of solar ultraviolet (UV) radiation causing heating. The UV radiation is highly dependent on the solar cycle and, therefore, the heating of the thermosphere, its temperature, its chemical composition and its neutral density.
Neutral density causes aerodynamic drag which, over long periods of time, leads to a decay of orbital altitude of CubeSats and eventually their demise upon re-entering denser atmospheric layers.
Orbital lifetime calculations of Qiao [23] enable the formulation of a rule of thumb for CubeSats: • at an altitude of 400 km, the orbital lifetime is around one year. This is typically the minimum lifetime for meaningful satellite operations • at an altitude of 600 km, the orbital life time is around 25 years. This is the maximum permissible lifetime according to the IADC Debris Mitigation Guidelines [24] Therefore, a suitable altitude range for CubeSats is 400-600 km. In the following, the analyses focus on altitudes from within that range.
Several atmospheric models exist. Here, we selected the NRLMSISE-00 model [25] as it is publicly available and provides the necessary information such as temperature, chemical composition, neutral density and variation with diurnal, annual and the repeating 11-year solar cycle. The  Fig. 2 [26]. Moreover, the model resolves variations with Earth's longitude and latitude.
To derive atmospheric conditions for our subsequent aerodynamic analyses, we investigate the variability of the atmosphere during solar cycle 23, which occurred between the years 1996 and 2007 (minimal to minimal solar activity). The atmospheric conditions at the beginning of this cycle will serve as reference for the simulations in the following sections. Figure 3 shows the chemical composition of the atmosphere in the desired 400-600 km range as well as the temperature according to the chosen model. It can be seen that within the considered altitude range, the atmosphere mostly consists of atomic oxygen, which decreases approximately exponentially with altitude. The exact chemical processes leading to such models are not fully understood and are thought to be responsible for inaccurate neutral density and thus drag predictions [27]. The temperature is constant over altitude at around 740 K. The pressure exerted onto a perpendicular surface is the change of momentum due to impact of the atmospheric particles per unit surface. If we assume the particles stick to the surface, the aerodynamic pressure p aero is: where ρ ∞ is the atmospheric density and v ∞ is the orbital velocity as per computed with Eq. (1).
For fully elastic impacts, the change of momentum is doubled and, therefore, also the pressure. The actual impact will be a blend of fully elastic and fully sticking gas-surface interactions (GSI). It will be addressed in the following sections. Figure 4 shows the aerodynamic ram pressure as a function of altitude and selected instants during solar cycle 23 on the night side of the Earth at the equator.
According to the atmospheric model, the aerodynamic pressure can change by almost two orders of magnitude over the solar cycle. Thus, also the available control forces change significantly over the solar cycle. Figure 5 shows the density's latitude dependence for the year 1996 low-solar-activity reference at night time and for the high solar activity 2004 during the day.
From the figure, it can be inferred that the latitude dependence is moderate-particularly for lower altitudes and the night time. Comparing the 2002 night time neutral densities shown in Fig. 5 with the 2002 night time neutral densities shown in Fig. 4, it can be seen that daytime densities are higher than night time densities as is expected for ideal gasses heated by solar radiation during the day.
During the solar cycle minima and night time, the atmosphere is less dense and, therefore, provides smaller aerodynamic forces for satellite orbit control. Thus, the low solar activity and night time atmospheric conditions serve as a conservative values and are used for our analyses. Control will be better for all other conditions. Ultimately, a control algorithm shall account for exact atmospheric conditions. Surface physics models used in the remainder of this study depend on the product atomic oxygen partial density and atmospheric temperature, n o T. Its variation over the solar cycle, latitude and altitude is shown in Fig. 6

Solar radiation pressure
Whereas the solar cycle causes high variability of UV radiation, which is responsible for the variability of the upper atmosphere, the total solar irradiance is counterintuitively largely constant [28]. Consequently, the solar radiation pressure is also invariable over the solar cycle. Its magnitude is shown for comparison purposes in Figs. 4 and 5. In analogy to the aerodynamic forces, the solar radiation pressure force depends on change of momentum during impact of the particle, which here is a photon. If it is fully absorbed, the pressure reads: The variable G is the solar constant and c is the speed of light. If the surface is fully reflective, the momentum change is doubled and, hence, also the solar radiation pressure.
In low earth orbit, a spacecraft experiences periods of eclipse during which solar radiation pressure is unavailable. For sun-synchronous dusk-dawn orbits (SSO, 6 am/pm local time of ascending/descending node [LTAN/LTDN]), these periods are minimal. They are computed and illustrated in Fig. 7 for the two extreme altitudes of those considered using Analytical Graphics Inc.'s (AGI) Systems Tool Kit version 12. The 23.5° day/night inclined terminator line can be seen. The orbital plane is inclined by 97° for the 400 km altitude  Figure 8 illustrates the duration of the eclipse over the year and for the considered altitude range. Eclipses are longer for higher altitudes. Accounting for the longer orbital periods in higher altitude as per Eq. (2), one can compute the percentage of eclipse. It is highest at the winter solstice for the lower altitude and amounts to 27%. In the following astrodynamics analyses, we use this conservative value for our further analyses. For flight software, the actual eclipse time should be used for optimal formation flight control.
The second source of radiation from the Sun, i.e., solar wind consisting of particles, is negligible and is, therefore, not accounted for here [29].
Within this research, we further neglect gravitational forces such as those by the J2 and J4 of the Earth and disturbances caused by the Moon and the Sun. As these forces are additional forces that can be exploited, our algorithm is, hence, conservative. If shown feasible, our algorithms and underlying assumptions suffice to solve the engineering problem. Future versions of our algorithm will make use of these forces for optimal formation flight control.

Surface physics
The forces exerted onto the satellites by the residual atmosphere and the solar light depend on the surface's properties. For both phenomena, the interactions are analogous: momentum is transferred from the impinging particle-either oxygen atom or a photon-to the spacecraft. In either case, the momentum transferred and, thus, the resultant pressure is dependent on whether the impact is plastic, i.e., the particle remains on the surface, or elastic, i.e., the particle is bounced back. In the latter case, the momentum transfer is doubled that of the former. The detailed surface physics differ for both phenomena, which is outlined below.

Gas-surface interaction
Particles impinging on a surface may be absorbed or may immediately bounce back. In the case of absorption, they may eventually desorb or react and leave the surface as a molecule with released bond energy [30]. However, this process is unlikely for the highly rarefied flow regime in LEO because of the scarcity of the reactants.
The actual gas-surface interaction will be a blend of the immediate back-bouncing-also known as specular reflection-and adsorption-desorption process. This was first modeled by Maxwell [31]. The process depends on the surface material-gas particle bonding energy and the temperature of the material.
The amount of energy accommodated, i.e., transferred, on the surface can be modeled as in [32] using the energy accommodation coefficient α: where T, T r and T W are the temperatures of the impinging atoms, the reflected atom and the surface wall, respectively. α becomes zero if reflected atoms have the same temperature as impinging atoms and are, therefore, not leaving while accommodating any energy on the surface. Conversely, α becomes unity if reflected atoms depart at wall temperature and, hence, accommodate almost all energy at the surface. As we will see later, in the first case, the aerodynamic forces are higher than in the second. Because higher forces are desirable for satellite control, materials with low energy accommodation are advantageous. Pilinski [33] proposes the following semi-empirical model for α: where n o is number density of atmospheric atomic oxygen and T is the temperature of the incident atom. For the relevant range of the parameter n o T as outlined in Sect. 2.1 and shown in Fig. 6, α is shown in Fig. 9.
For the solar minimum and night time, when the atmosphere is thinner, the parameter α is very small or negligible. The decreased aerodynamic forces due to the thin atmosphere are partially offset by better aerodynamic performance due to a lower energy accommodation coefficient, α.

Optical surface properties
The interaction of light, i.e., photons, with the surface of the spacecraft depend on its optical properties. There are four main mechanisms: absorption a, diffusive and specular reflection r D , r S and transmission t. As energy is conserved, their sum equals one, as in Eq. (7): Solar panel covers are glassy and allow transmission of light. However, the underlying material is nontransparent. Hence, effectively, the transmissivity of solar cells is zero.
List [34] provides typical surface property values for common solar cells, which are tabulated in Table 1 differentiated for beginning-of-life (BOL) and end-of-life properties.
Solar radiation pressure forces are higher for higher reflectivities. Materials with higher reflectivities are, therefore, advantageous. However, satellite surfaces are typically entirely covered with solar cells to maximize power generation. To this end, absorption of the solar light is necessary which is a directly conflicting surface property requirement.

Solar-aerodynamic forces
Sentman [35] provides equations to determine the aerodynamic force coefficients for the rarefied flight regime using the energy accommodation coefficient as input. They are repeated in Table 2 for convenience. Here, θ is the angle between the surface normal and the incoming particle direction, v ∞ is the speed of the impinging particle, α is the energy accommodation coefficient as per Eq. (6), R is the ideal gas constant and T i is atmospheric temperature.
Similarly, List [34] provides with the equations forces exerted on a surface due to solar radiation pressure, which are given in Table 3. p SRP solar radiation pressure as per Eq. (4), r D and r S are optical surface properties as per Eq. (7) and Table 1; A is the surface area.
The outcome of these equations, i.e., the solar and aerodynamic force coefficients, is given in Figs. 10 and 11.

Solar-Aerodynamic formation flight
For autonomous formation flight, we implemented the Hill-Clohessy-Wiltshire (HCW) equations in the notations of Ivanov [37]. They model the local movement of the satellites relative to each other.
The HCW equations are a set of ordinary differential equations with a right-hand side accounting for external forces, which are in our case those of aerodynamics and of solar radiation pressure.
In contrast to the approaches of others, we implemented rotations around all three axes using Euler angles, i.e., roll, pitch and yaw. We apply an often employed convention using two extrinsic and one intrinsic rotation. The rotation axes are shown in Fig. 1. For this research, permissible angles are multiples of 45°. In total, four units of aerodynamic surfaces are considered for a 1U CubeSat. Such surfaces consisting of a frame, solar cells and a spring-enabled deployment mechanisms are commercially available and frequently used for CubeSats. A conceptual drawing is given in Fig. 1.

Control law
We employ a classical Linear-Quadratic-Regulator (LQR) control algorithm following the example of Ivanov [37]. It controls at each time step all three axes. The regulator determines the control forces required for optimal control. Similarly to Ivanov, we employ a distributed formation flight control method: • The maximum x-component of the three computed control forces (per satellite) is determined. The amount is subtracted from each satellite's x-component control force. Thus, the updated x-component control force is negative or zero. This shift is needed because drag can Table 2 Aerodynamic coefficients for rarefied orbital flight regime Table 3 Force induced by solar radiation pressure only provide forces opposite of the direction of flight (negative forces). • For the directions perpendicular to the x-axis, i.e., y and z, the average control force is computed and subtracted from the corresponding control force for each satellite. This allows minimizing the required control force for each satellite. It benefits the formation flight control as the available forces are typically smaller than the required control force.
The shift and averaging of required control forces require the exchange of this information from one satellite to another through an inter-satellite communication link. To this end, the normally available telemetry and telecommand communication system can be used. By design, its range is at least 3000 km as needed for LEO-ground communication. Thus, the range is vastly sufficient for the purposes considered here.
A routine computes the available solar-aerodynamic forces for the given 45° granularity of the Euler angles (see Fig. 1). Restricting the permissible Euler angles prevents the controller from performing optimally and reduces stability because the required direction of the required control cannot usually not be achieved. However, this restriction also reduces the computational load, which is a scarce resource on a nanosatellite. This engineering choice has been seen favorable for our simulations. A detailed analysis will be made in our future research.
We implemented the methods in MATLAB 2019 and executed the code on a standard office laptop with Intel I7 CPU in one computational thread. Run time was about 3 h per flight case (altitude, modeling) as presented below.

Design test case
Orbital control for formation flight is required for four situations: 1. Deployment, i.e., the establishment of the default formation geometry from the initial configuration after launch. Typically, the satellites are in the same orbital plane and close to each other in the beginning. 2. Reconfiguration, i.e., establishment of a formation geometry from an existing formation geometry.
3. Re-establishment of the formation geometry from an arbitrary configuration; for instance, after a formation flight control anomaly 4. Maintenance, i.e., the keeping of the formation geometry in view of disturbances, control sensor uncertainty and actuator imperfections. This situation is the nominal one during which routine operations for accomplishing the mission objectives, for instance sensor use, is carried out.
Within this research, we focus on the first and demonstrate our algorithm with a mission scenario we coined Cubesat of Luxembourg's University for Space Technology and Earth Research (CLUSTER). It consists of three 1U CubeSats. The choice was made because, conventionally, a CubeSat deployer allows the ejection of three CubeSat units. Our three satellites would be simultaneously ejected with a very small differential speed dictated by the required separation springs [38] between the individual CubeSats. The algorithm is capable to simulate an arbitrary number of satellites. The main parameters for this mission are listed in Table 4. The demonstration formation geometry consists of, in Local Horizontal Local Vertical (LHLV) coordinates, two satellites flying at constant distance and a third satellite circling an intermediate position. The target formation geometry is illustrated in Fig. 12.
In Earth-centered coordinates, satellites 1 and 3 share the same Kepler elements except for the true anomaly. Satellite 2′s eccentricity and right ascension node is slightly different from those of the other two satellites leading to the circulating motion in LHLV. The formation features particular small satellite distances. Applications of such small formations are numerous. For instance, distributed sensors architectures requires such geometries. The small-distance flight requires the accounting of small disturbances such as aerodynamic drag and solar radiation pressure. In this article, we show not only how to integrate the disturbance but also how to make use of it for formation flight control. Figure 13 shows the trajectories of the three satellites relative to each other over a period of 120 orbital periods as computed by our algorithm that we coined CosmosAlpha in the local 3D LHLV coordinate system. All three satellites start at the origin. Satellite 3 moves toward the -x direction, overshoots its target location and then returns to it at x = − 230 m. The overshoot is due to the imperfect control algorithm; however, its robustness is demonstrated by the satellite's eventual return. Satellite 2 starts to move out of the z-x plane approaching the target relative circular movement with a period of one orbital period (compare the results of the solution obtained the control algorithm shown in Fig. 13 to the analytical solution in Fig. 12).

Results: effect of solar radiation pressure
In the following, we focus on analyzing the trajectory of satellite 2 only. The satellite moves in the 3D space and is, therefore, particularly suitable to show the features of our algorithm. Figure 14, shows the coordinates x, y and z (from top to bottom) over time of the position of the 2nd satellite for different altitudes in between the 400 km and 600 km limits. Comparing the computational cases in which only aerodynamic forces are modeled (dashed lines), it clearly can be seen that satellite 2 moves faster to its target trajectory at lower altitudes. This is expected since the denser atmosphere lends to higher aerodynamic forces available to the formation flight controller. Figure 14 shows 120 orbits of the flight for the x-coordinate. For the y and z coordinates, only a zoom to the first 20 orbits is shown for clarity of the display.
The solid lines show the trajectories of satellite 2 with modeled solar radiation pressure. Direction of the solar radiation and eclipse time correspond to a 6 am LTAN orbit and is, therefore, approximately in the direction of the negative y-axis. It can be seen that exploiting solar radiation pressure in addition to aerodynamic forces leads to a faster establishment of the formation geometry.
It can also be seen that at the upper limit of the 400-600 km range, the satellite does not reach the target trajectory. It is a consequence of very small available aerodynamic forces. Solar-aerodynamic formation flight would require larger control surfaces for formation flight.

Verification of the effect of eclipse
The simulations for which results are presented in the previous sections are conservative as they model the loss of solar radiation pressure due to eclipse for a fraction of the orbital duration. For the considered 6:00 am LTAN orbit, eclipse only takes place around the winter solstice as illustrated in Sect. 2.2, Fig. 8. The effect of eclipse is investigated exemplarily here for one altitude, i.e., 525 km. Figure 15 shows the trajectory of the second satellite with modeled eclipse and without. As before, 120 orbits are shown for the x-coordinate, but only 20 orbits for the y-and z-coordinate for clarity. For comparison purposes, we also show the trajectory not accounting for solar radiation pressure. It can be seen that the accounting for eclipse results in a slightly shortened duration of the formation deployment. The results are in-line with the expectations as the eclipse period is relatively short. Fig. 13 Trajectories of satellite 2 (beige) and 3 (magenta) relative to satellite 1 computed with our algorithm using solar-aerodynamic forces at an altitude of 525 km

Summary and conclusions
We report on our progress in developing the formation flight algorithm coined Cosmos in its version Alpha. We accounted for the first time for the presence and use of solar radiation pressure and show that the force is of significance compared to other disturbances in low earth orbit. We demonstrate that solar radiation pressure can be used to control the flight of satellites to establish a CubeSat formation in conjunction with aerodynamic forces. Emphasis was put on the use of conservative assumptions for these forces and their duration. Additional forces available for formation flight control were also neglected. Hence, our algorithm is realistic. A future version accounting for more effects will be higher performing.
Our simulations show that the formation geometry is reached faster with modeling of the solar radiation pressure. Consequently, it can also be concluded that neglecting of the solar radiation pressure leads to suboptimal formation flight control.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.