A continuum model for the orbit evolution of self-propelled ‘smart dust’ swarms

A continuity equation is developed to model the evolution of a swarm of self-propelled ‘smart dust’ devices in heliocentric orbit driven by solar radiation pressure. These devices are assumed to be MEMs-scale (micro-electromechanical systems) with a large area-to-mass ratio. For large numbers of devices it will be assumed that a continuum approximation can be used to model their orbit evolution. The families of closed-form solutions to the resulting swarm continuity equation then represent the evolution of the number density of devices as a function of both position and time from a set of initial data. Forcing terms are also considered which model swarm sources and sinks (device deposition and device failure). The closed-form solutions presented for the swarm number density provide insights into the behaviour of swarms of self-propelled ‘smart dust’ devices an can form the basis of more complex mission design methodologies.


Introduction
Astrodynamics has traditionally been concerned with the orbital motion of artificial satellites described using sets of nonlinear ordinary differential equations (ODE). For example, the inverse square force due to the potential of a uniform, spherical mass leads to the well-known family of conic section orbits through the Kepler problem. Using the Kepler problem as the unperturbed solution to the artificial satellite problem, the influence of perturbations can then be investigated by constructing high order series solutions using a range of techniques. However, while nonlinear ODEs provide a set of powerful tools to investigate the rich behaviour of the artificial satellite problem, new tools are now required to understand the behaviour of B Colin R. McInnes colin.mcinnes@glasgow.ac.uk large ensembles of micro-satellites of diminishing length-scale. Rather than using the methods of ODEs to propagate the motion of individual satellites, a continuum approach using partial differential equations (PDE) can be used to provide an understanding of the behaviour of these new systems.
While an understanding of the behaviour of swarms of micro-satellites provides the contemporary motivation for a continuum approach to Astrodynamics, prior work on continuum dynamics for artificial satellites has been motivated by the need to understand the evolution of orbital debris clouds and early passive systems such as the West Ford needles [e.g. Heard (1976)]. In addition, prior literature on the orbit evolution of clouds of natural dust grains [e.g. Gor'kavyi et al. (1997)] and galactic dynamics [e.g. Binney and Tremaine (1987)] provides rich insights into the problem.
However, the need for new methodologies can be seen in the development of self-propelled MEMs-scale (micro-electromechanical systems) technologies which can enable the deployment of large swarms of self-propelled 'smart dust' devices (Kahn et al. 1999) driven by solar radiation pressure (Atchison and Peck 2008;Barker and Rodriguez-Salazar 2008;Atchison and Peck 2010). It is envisaged that swarms of such devices could be released from a dispenser (Manchester et al. 2013) and then utilise solar radiation pressure to slowly spiral inwards or outwards through the ecliptic plane (McInnes 1999). Distributed swarms could in principle deliver massively parallel sensing for space science applications by providing, for example, simultaneous multi-point measurements of magnetic field strength and direction (McInnes 1999;Barnhart et al. 2007;Atchison and Peck 2010).
In addition to the development of practical technologies for 'smart dust' devices (Manchester et al. 2013), recent work has led to an initial understanding of their orbit evolution (Colombo and McInnes 2011a, b;Lucking et al. 2012). This is a particularly intriguing problem given the impact of length-scaling (Atchison and Peck 2011); for a device of lengthscale L, its mass scales as L 3 while its surface area scales of L 2 , resulting in its area-to-mass ratio increasing with decreasing device size as L −1 . Smart dust devices are therefore strongly perturbed by surface forces such as aerodynamic drag and solar radiation pressure (Atchison and Peck 2011;McInnes et al. 2011).
In this paper a continuum model of a swarm of smart dust devices in heliocentric orbit is developed where the devices are self-propelled by solar radiation pressure. Again, rather than model the orbit evolution of a single device, a continuity equation is formed which represents the number density of devices. Families of closed-form solutions to the continuity equation are then found representing the orbit evolution of the swarm through its number density as function of both position and time. The work presented builds on prior analysis of the evolution of debris clouds (McInnes 1993;Letizia et al. 2015) and the evolution of geocentric micro-satellite swarms perturbed by air drag (McInnes 2000;Manchester and Peck 2011;McInnes and Colombo 2013). Importantly, since a continuum approach with a PDE is used, initial data at a given boundary is required to represent the spatial distribution of devices, with the characteristic curves of the PDE mapping the initial data onto the time evolution of the swarm. This is a quite different procedure to classical Astrodynamics where point initial conditions are required for a set of ODEs.
In order to develop solutions to the continuity equation, the dynamics of a single smart dust device are considered in Sect. 2 and approximate quasi-circular spiral trajectories obtained for a fixed device attitude. These trajectories later provide the characteristic curves of the continuity equation. Then, in Sect. 3 an azimuthally symmetric continuity equation is used to model the evolution of the number density of smart dust devices as a function of both orbit radius and time for a range of scenarios. These include consideration of a source term, representing the constant deposition of new devices, and a sink term representing the failure of devices during their lifetime, as presented in Sect. 4. Rarefaction and compression of the flow of devices is found, as expected, since the radial flow speed of the swarm due to solar radiation pressure is a strong function of orbit radius, while the swam either flows into a smaller (inward spiral) or larger (outward spiral) volume.
Extending the radial flow problem, azimuthally asymmetric flows are considered in Sect. 5 to model the orbit evolution of a given two-dimensional spatial pattern of smart dust devices. Here azimuthal shearing of the pattern is found, again as expected, since the azimuthal flow speed of the swarm is also a strong function of orbit radius distance through Kepler's third law. Lastly, some conclusions are drawn in Sect. 6 where it is proposed that the closed-form solutions obtained for both the radial flows and azimuthally asymmetric flows can provide insights into the orbit evolution of swarms and form the basis of future mission design tools.

Single device orbit evolution
In order to proceed, a simple model of a single smart dust device will be considered which is self-propelled by solar radiation pressure. It will be assumed that the device is a small, thin flat panel (Atchison and Peck 2010) which has non-perfect specular reflectivity η and that the front and rear emissivities are equal. The direct solar radiation pressure induced acceleration acting on such devices is significantly larger than the Poynting-Robertson drag, which can be neglected. Following McInnes (2014), the device acceleration a induced by solar radiation pressure can be written as with unit radial vector r and unit transverse vector t. In addition P is the solar radiation pressure, A is the device reference area, m is the device mass and α is the pitch angle of the device relative to the Sun-line in the range −π/2 ≤ α ≤ π/2, as shown in Fig. 1. It will be assumed that the device pitch angle is fixed, which can in principle be achieved by passive means through etching of the device surface into a sawtooth pattern (Atchison and Peck 2010). Again following McInnes (2014) the rate of change of (specific) orbit energy E of the device is given by The two-body orbit energy is defined by E = −μ (1 − β) /2r (t) for orbit radius r (t) and gravitational parameter μ. While the transverse component of acceleration a defined by Eq. (1) delivers work and so changes the device orbit energy, the radial component of acceleration can be incorporated into a reduced gravitational parameter (1 − β), as noted by Klačka (2014). For a slow inward or outward spiral, and so a quasi-circular orbit, the device velocity vector is approximated by the local Keplerian circular orbit velocity v(μ (1 − β) /r (t)) 1/2 t (Prussing and Conway 2012), where again the radial component of acceleration can be incorporated in a reduced gravitational parameter. Furthermore, the inverse square solar radiation pressure P is defined by P =P (r /r (t)) 2 , whereP is the solar radiation pressure at some distancer . The device lightness number is now defined as the ratio of the solar radiation pressure induced acceleration to the solar gravitational acceleration experienced by an ideal device with reflectivity η = 1. Since both solar gravity and solar radiation pressure scale as the inverse square of the device orbit radius r (t), the lightness number is found to be β = 2P Ar 2 /mμ. A more sophisticated analysis for dust grains is given by Klačka (2014). (2) it can be shown that the device orbit radius evolves as If the device has an initial orbit at heliocentric radiusr , such that ξ (0) = 1, the device orbit radius can then be obtained by integrating Eq. (3) so that Again assuming a slow quasi-circular spiral, the polar angle θ (τ ) can also be obtained from Kepler's third law such that Integrating Eq. (5) it can then be seen that where the initial device location is such that θ (τ ) = 0. From Eq.
(3) it can be seen that the rate of change of device orbit radius can be maximized by choosing the device pitch angle such that d dα to maximize the transverse acceleration experienced by the device, with α + ∼ +35 • and α − ∼ −35 • , representing the pitch attitude required for an outward and inward spiral respectively. Now that the orbit evolution of a single device has been considered, the evolution of a swarm of such devices can be investigated through a continuity equation, the family of solutions to which will provide the evolution of the swarm number density.

Swarm continuity equation
It will now be assumed that a swarm of identical devices is distributed in space within some bound volume and that the spatial distribution of the swarm can be represented by the swarm number density. At some initial time the swarm number density will therefore represent the so-called initial data of this continuum problem. Solutions to the continuity equation will propagate the initial data forward to describe the evolution of the swarm number density as a function of both position and time.
First, an azimuthally symmetric swarm will be considered whose number density is a function of orbit radius and time only. In order for this assumption to hold it can be noted that if the timescale for the inward (or outward) flow of the devices is slow relative to the local circular orbit period, the swarm can be approximated as being an azimuthally symmetric disk. An asymmetric swarm whose number density is a function of orbit radius, polar angle and time will be considered later in Sect. 5. In order to bound this approximation it can be seen from Eq. (6) that an individual self-propelled device will drift azimuthally along a quasi-circular spiral such that θ (τ ) = log (1 + λτ ) /λ, while a dispenser on a circular Keplerian orbit at unit radius will move uniformly such that θ (τ ) = τ . The difference in the azimuthal position of the device and the dispenser can therefore be approximated by δθ = τ − log (1 + λτ ) /λ ≈ λτ 2 /2. In order for the device to complete one revolution relative to the dispenser, such that δθ = 2π requires a duration τ ≈ √ 4π/λ. For a device with lightness number β = 0.01 (Atchison and Peck 2010) and α = α + this corresponds to a duration equivalent to 6.5 orbits of the dispenser. As will be seen later in Sect. 4, for β 1 the evolution of the swarm over much longer timescales will be considered, so the initial approximation of azimuthal symmetry is reasonable.
An individual device will now be located using (non-dimensional) cylindrical polar coordinates(ξ, θ), as shown in Fig. 2, with velocity components v ξ , v θ defined through Eq. (3) and Eq. (5). Due to the action of solar radiation pressure, its radial velocity v ξ will induce a slow inward (or outward) drift. Following McInnes (2000) a control volume will be defined as an annulus in the radial interval (ξ, ξ + ξ ) which will contain a total number of devices 2πξ ξn (ξ, τ ), for a given swarm number density n (ξ, τ ). The rate of change of the swarm number density can therefore be obtained by considering the flow of devices into and out of the control volume, along with any source and sink terms representing the deposition or removal of devices. To conserve the total number of devices it is required that whereṅ − (ξ, τ ) is a sink term, representing the rate of removal of devices whileṅ + (ξ, τ ) is a source term, representing the rate of deposition new devices. The radial velocity term v ξ (ξ + ξ ) and number density term n (ξ + ξ, τ ) in Eq. (8) can be Taylor expanded to first order. Taking the limit ξ → 0 it can be shown that  (2000), differentiating out the second term in Eq. (9) and noting that the radial velocity is a function of orbit radius only, the continuity equation may be written as ∂n (ξ, τ ) ∂τ where ( ) represents a differentiation with respect to ξ . It should be noted that the total number of physical devices N in an annulus of radius ξ and width ξ is 2πξ ξn (ξ, τ ), so that the swarm density n (ξ, τ ) can be interpreted as N /2πξ ξ, which is an areal density since a disk is being considered. In addition, for a fixed number of physical devices N and annulus width ξ the density will therefore scale as ξ −1 , as will be discussed in Sect. 4.3. It should also be noted that the choice of an annulus as control volume is driven by the symmetry of the problem and the representation of the continuity equation in polar coordinate form. However, a global conservation law can be written in integral form which can then be represented locally in any suitable coordinate system, as is common for other continuum problems (Anderson 2001).

Solutions to the swarm continuity equation
The continuity equation representing the evolution of the number density of devices as a function of orbit radius and time can now be solved using the method of characteristics (Murphy 1993). Along the characteristic curve, the partial differential equation becomes an ordinary differential equation amenable to solution by standard methods. This yields a solution to the original partial differential equation when combined with suitable initial data at a given boundary, for example the initial swarm density function n (ξ, 0). The method therefore transforms the single partial differential equation defined by Eq. (10) into two ordinary differential equations defined by Firstly, integrating Eq. (11a) yields the characteristic curves of the partial differential equation and substituting from Eq. (3) yields the family of characteristics G (ξ, τ ) defined by where C 1 = λC. Then, assuming for the present that the sink termṅ − (ξ, τ ) = 0 and the source termṅ + (ξ, τ ) = 0 in Eq. (11b), the device number density can be obtained from ( 1 4 ) which can be solved directly to obtain To obtain the general solution to the partial differential equation defined by Eq. (10), an arbitrary function of the characteristic equation ψ (G (ξ, τ )) is added in Eq. (15) which is determined from the initial data of the problem n (ξ, 0) using standard methods. The general solution can therefore be written as where (G (ξ, τ )) = ex p (ψ (G (ξ, τ ))). Then, at τ = 0 it can be seen from Eq. (13) that G (ξ, 0) = ξ 3/2 so defining z = G (ξ, 0) and so ξ (z) = z 2/3 the functional form of (G (ξ, τ )) can be obtained from which again is determined from the initial data of the problem n (ξ, 0). Now that the general solution to the continuity equation has been obtained a number of specific cases will be considered with suitable initial data representing each scenario. While Eq. (16) represents the solution to the continuity equation withṅ − (ξ, τ ) = 0 andṅ + (ξ, τ ) = 0 the addition of source and sink terms will be considered later.

Evolution of an infinite sheet
First, for illustration, an infinite sheet of devices will be considered with initial data n (ξ, 0) = 1 in the range 0 < ξ < ∞. This represents an ensemble of smart dust devices uniformly distributed throughout the entire domain of the continuity equation. Using this initial data In order to illustrate the form of the solution, the lightness number of each smart dust device is now defined as β = 0.01 (Atchison and Peck 2010) and the device attitude is assumed to be α = α − , so that the devices flow inwards. First, the family of characteristic curves defined by Eq. (13) are shown in Fig. 3. Again, the characteristics will map the initial data n (ξ, 0) = 1 to form the solution surface of the continuity equation n (ξ, τ ). The solution surface itself is shown in Fig. 4 for 0.1 ≤ ξ ≤ 1, where dimensional units are used for illustration with unit radius corresponding to 1 AU. It can be seen that the disk interior to the Earth's orbit is uniformly filled, but as the ensemble of smart dust devices flow inwards from outside the Earth's orbit the initially density grows, particularly at small orbit radii. The physical interpretation of this scaling of the swarm density with orbit radius is again discussed later in Sect. 4.3.

Evolution of a finite disk
A uniform disk of devices will now be considered with initial data n (ξ, 0) = 1, but in the range 0 < ξ ≤ 1. This represents an ensemble of smart dust devices uniformly distributed interior to a circle of unit radius. Using this initial data the arbitrary function of the characteristic equation in Eq. (16) can again be determined. Then, using the methods detailed above, it can be shown that the resulting solution to the swarm continuity equation is given by where H is the Heaviside step function. The solution surface defined by Eq. (20) is shown in Fig. 5 where it can be seen that the disk interior to the Earth's orbit is again initially uniform, but as the ensemble of smart dust devices flow inwards the disk contracts in radius while its density grows. The contraction of the outer edge of the disk is captured by the step function in Eq. (20), while the scaling of the density follows the same functional form as Eq. (19).

Constant device density at one boundary
Rather than a uniform disk, an ensemble of smart dust devices will now be considered with the device density held constant at unit radius. The initial data for this modified problem is therefore n (1, τ ) = 1 in the range 0 < ξ ≤ 1. Again, the arbitrary function of the characteristic equation in Eq. (16) can be determined from this initial data and it can be shown that the solution to the swarm continuity equation is now given simply by Here, the solution surface is time independent since the initial data on the boundary is such that n (1, τ ) = 1 for all time, while it can be seen that the growth in density with decreasing orbit radius scales as ξ −1/2 . To understand this scaling law, if the swarm density is time independent Eq. (9) yields d dξ ξv ξ (ξ ) n (ξ ) = 0 ( 2 2 )

Fig. 5
Solution surface n (ξ, τ ) for a uniform disk with initial data n (ξ, 0) = 1 in the range 0 < ξ ≤ 1, and with lightness number β = 0.01, η = 1 and device attitude α = α − so that ξv ξ (ξ ) n (ξ ) =C, for some constantC. However, it is then clear thaṫ is also constant, whereṄ (ξ ) is the total flux of smart dust devices crossing an annulus of radius ξ . Therefore, due to the requirement for continuity to conserve the total number of devices ξv ξ (ξ ) n (ξ ) = v ξ (1) n (1) and so from Eq. (24) and Eq. (3) it can be seen that which is equivalent to the result obtained by the method of characteristics presented in Eq.
(20). The resulting solution surface defined by Eq. (21) is shown in Fig. 6 which displays growing density with decreasing orbit radius as expected.
The continuity argument resulting in Eq. (25) also arises due to the geometry of the problem, since the circumference 2πξ of any annular boundary of radius ξ shrinks faster than the inflow speed v ξ (ξ ) increases with diminishing orbit radius. Indeed for the general solution defined by Eq. (16) the radial inflow of devices is faster at smaller orbit radii which acts to reduce the swarm density since the density n (ξ, τ ) scales as v ξ (ξ ) −1 . However, v ξ (ξ ) scales as ξ −1/2 while the compression due to the geometric term in Eq. (16) discussed in Sect. 3 scales as ξ −1 so that there is an overall density scaling of ξ −1/2 .

Constant device deposition rate at one boundary
While the uniform disk provides useful insights into the solution to the continuity equation it is clearly a rather unrealistic scenario. A more realistic scenario will now be considered with a dispenser in orbit at unit radius ξ = 1 which forms a source term in Eq. (11b). Although a single dispenser will form a point source, an azimuthally symmetric swarm will again be considered such that the timescale for the inward (or outward) flow of the devices is slow relative to the local circular orbit period, as discussed in Sect. 3. Therefore, adding a source termṅ + (ξ, τ ) =Cδ (ξ − 1) in Eq. (9) the swarm continuity equation is now given by ∂n (ξ, τ ) ∂τ whereC is the magnitude of the source term, to be determined later, and δ (ξ ) is the Dirac delta function. In order to obtain a solution to the continuity equation with a source term, Eq. (11b) must firstly be solved with the non-homogeneous forcing term such that Since the equation is non-homogeneous, and integrating factor I is required such that and so Eq. (27) reduces to which then yields Then, with initial data n (ξ, 0) = 0 in the range 0 < ξ ≤ 1 the arbitrary function of the characteristic equation in Eq. (30) can be determined so that where again H is the Heaviside step function. In order to understand the physical significance of the magnitude of the source termC the total deposition rate of devicesṄ can be determined froṁ so thatC =Ṅ /2π. Alternatively, from Eq. (31) the swarm density can be defined at the boundary ξ = 1 such that n (1, τ ) = 1 for τ = 0 and soC = 2λ/3. The solution surface corresponding to Eq. (31) with n (1, τ ) = 1 is shown in Fig. 7. It can be seen that n (ξ, 0) = 0 so that the swarm density initially vanishes everywhere while the deposition term introduces devices at ξ = 1 which then flow inwards eventually filling the space interior to the unit disk 0 < ξ ≤ 1. From Eq. (4) it can be seen that the disk is filled after a time τ = 1/λ. Again, the density of the swarm scales spatially as ξ −1/2 as expected. It can be seen from Eq. (31) that forC = 2λ/3 the evolution of the swarm density has an asymptotic form which represents the steady-state condition with equilibrium between the source term at ξ = 1 and an effective sink at ξ = 0.

Swarm evolution with on-orbit failures
In addition to a source term representing the deposition of new devices, a sink terṁ n − (ξ, τ ) = −εn (ξ, τ ) can be added to the swarm continuity equation to represent failure of devices within the swarm (McInnes 2000). From Eq. (26) the new continuity equation is therefore given by where 1/ε represents the mean device failure rate. With initial data n (ξ, 0) = 0 in the range 0 < ξ ≤ 1 it can then be shown that swarm density is given by

Fig. 7
Solution surface n (ξ, τ ) for a swarm with deposition rate fixed at one boundary with initial data n (ξ, 0) = 0 in the range 0 < ξ ≤ 1, and with lightness number β = 0.01, η = 1 and device attitude α = α − SelectingC such that n (1, τ ) = 1 and selecting ε such that the mean device failure rate is 10 years, the solution surface corresponding to Eq. (35) is shown in Fig. 8. It can be seen that the swarm density falls due to device failures, but that this is offset by the geometric scaling (the devices are flowing into annuli of diminishing circumference) leading to a growth in density at small orbit radii.

Extension to two-dimensional system
The analysis of an azimuthally symmetric swarm presented in Sect. 4 can now be extended to a two-dimensional planar system using polar coordinates. Here the swarm of devices flows in the radial direction, but the transverse motion of the swarm is also captured. Given that the azimuthal flow speed of the swarm is a function of orbit radius through Kepler's third law, this implies that a swarm of devices will be sheared in the transverse direction, as will be seen later. It can be shown by extending the continuity argument of Sect. 3 to include flow in the azimuthal direction that the new continuity equation is given by where n (ξ, θ, τ ) is again the device density and v θ (ξ ) is the local transverse speed at radius ξ . Again, the method of characteristics transforms the single partial differential equation defined by Eq. (36) into coupled ordinary differential equations defined by The new continuity equation therefore possess an additional characteristic curve which can now be obtained from and so substituting from Eq. (4) and Eq. (5) yields the family of characteristics H (θ, τ ) defined by The additional family of characteristic curves are shown in Fig. 9, where again the lightness number of each device is defined as β = 0.01 and the device attitude is α = α − , so that the devices flow inwards. Neglecting the source and sink terms it can then be shown that the general solution to the new continuity equation is given by n (ξ, θ, τ ) = 1 ξv ξ (ξ ) (G (ξ, τ ) , H (θ, τ )) (40) where (G (ξ, τ ) , H (θ, τ )) is a function of the two characteristics defined by Eq. (13) and Eq. (39) which can again be determined from the initial data. In order to proceed, for illustration the initial data n (ξ, θ, 0) will now be represented by the unit step function where the step functionU (z 1 , z 2 ) is defined as The initial data defined by Eq. (41) therefore represents a region of unit density centred on (ξ o , θ o ) and of side ( ξ, θ ) which can be used to determine the function form of (G (ξ, τ ) , H (θ, τ )). Again for illustration, a swarm of unit density centred on ξ o = 1 and θ o = 0 with width ξ = 1/4 and θ = π/4 will now be propagated using the solution to the two-dimensional continuity equation. Figure 10a shows inward flow (α = α − ) while Fig. 10b shows outward flow (α = α + ). For ease of illustration a lightness number β = 0.1 is assumed in order to clearly visualise the evolution of the swarm. For inward flow it can be seen that the initial region of unit density shears in the azimuthal direction while its radial extent is stretched resulting in a thin filament which wraps around itself, since the azimuthal and radial flow speeds are functions of orbit radius. It can also be seen that the density increases at the inner edge of the filament since the density of the swarm scales radially as ξ −1/2 . For outward flow it can be seen that the azimuthal shearing is less pronounced while the radial extent of the swarm compresses. Since the radial flow speed is an inverse function of orbit radius the inner edge of the swarm begins to catch up with the outer edge, resulting in increased number density.

Conclusions
Using a continuum approximation, a continuity equation has been developed to represent the evolution of a swarm of self-propelled 'smart dust' devices in heliocentric orbit under the action of solar radiation pressure. The family of solutions to the continuity equation describe the evolution of the swarm from a set of initial data at some boundary. By considering both radial and transverse motion, the spatial evolution of the swarm can be obtained and azimuthal shearing observed since the local transverse speed is a function of orbit radius. The addition of source and sinks terms to the continuity equation allows the deposition of new devices and on-orbit failures to be considered. It is proposed that the closed-form solutions presented for both the radial flows and azimuthally asymmetric flows of devices can provide useful insights into the orbit evolution of swarms of self-propelled 'smart dust' devices which can form the basis of future mission design tools. reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.