Abstract
Confidence regions for spacecraft state can be constructed in phase space which encapsulate some region where there is a likelihood for the state to reside. These regions can be treated as phase space distributions or structures. Structures, such as surfaces or volumes, are constrained to preserve specific properties as they evolve in phase space under Hamiltonian dynamics. Thus, spacecraft uncertainty is then constrained by Hamiltonian flow which can provide insight into state determination. This work examines the modified constraints in the presence of non-conservative forces which relate to both probabilistic and geometric properties of the evolving uncertainty structure. The modified constraints are then derived for a Two-Body and drag environment and are shown to be valid after comparison with alternative methods. Applying the modified constraints, the constrained evolution of the confidence region is then tied to a simple physical explanation for the changing knowledge in our spacecraft state, in the atmospheric drag environment and Poynting–Robertson drag environment.
Similar content being viewed by others
1 Introduction
As the Low Earth Orbit (LEO) space environment becomes increasingly occupied by commercial and government spacecraft (McDowell 2020), an in-depth understanding of orbit uncertainty becomes more valuable. In particular, constraints on the states of space objects allow operators to mitigate collision probability with intelligently designed guidance, navigation, and control. These objects may include active spacecraft and debris lacking propulsion capabilities. Several contributions highlight practical solutions to orbit uncertainty propagation for various orbital environments. Park and Scheeres (2006, 2007) outline how higher order expansions of the spacecraft solution flow, referred to as state transition tensors, can be utilized to propagate mean state perturbations and state covariance matrices. These tensors act as higher order corrections to the standard linear state transition matrix update. Moreover, Boone and McMahon (2021) showcase a computationally efficient method to compute such state transition tensors in highly non-linear environments.
Another approach for spacecraft uncertainty propagation was addressed in the context of heliocentric asteroid and near-Earth debris tracking. Only a substate of the detected object can be specified following the typical series of measurements in a given pass. Hence, some component of the state is uncertain which can be restricted by the dynamical constraints of the problem. These restrictions confine the unknown subcomponents of the state to a specific uncertainty surface or admissible region of phase space. From various passes of this object over our ground station, multiple measurements can be correlated to each other. By discretizing the uncertainty region into a grid of potential or “virtual” particles, the admissible regions for each measurement can be dynamically mapped to the same epoch. Based on the overlap of the admissible regions, virtual particles can be ranked and one can be assigned as the most probable orbit for the true particle (Milani et al. 2004; Tommei et al. 2007; Maruskin 2008).
The previously mentioned uncertainty propagation techniques are more focused around orbit determination where the uncertainty is treated either as a collection of “microscopic” particles or a matrix. While these treatments of uncertainty can provide insight into uncertainty evolution, they do not capture “macroscopic” aspects of uncertainty evolution itself which further constrain the potential orbits. To capture these aspects, more dynamical systems theory approaches have been taken to map uncertainty forward in time. Hamiltonian dynamical systems allow for useful constraints (Goldstein et al. 2002) to be applied to spacecraft state uncertainty distributions. Under Hamiltonian formalism, Liouville’s theorem and the other Poincaré-Cartan integral invariants along with Gromov’s Non-squeezing theorem can restrict various uncertainty structures in 2n-dimensional phase space. For instance, Scheeres et al. (2014) and Scheeres et al. (2012) demonstrate how these constraints and topological properties can provide an invariant lower bound on the covariance of a spacecraft state, and how such a bound is updated with measurements. Maruskin et al. (2007) showcases the power of this analysis with application to admissible regions and asteroid tracking in order to constrain uncertain asteroid orbits at a future time.
The dynamical systems theory approaches above assume the systems of consideration are Hamiltonian and all forces acting on particles embedded in the system are conservative. With relation to spacecraft applications, there are several orbital environments where this condition does not hold.
Over the past few decades, humans have drastically accelerated their presence in space ranging from enhanced communication capabilities with new communication satellites to mega-constellations injected into LEO (McDowell 2020). For the purposes of Space Situational Awareness (SSA) and tracking of these objects in addition to space debris, there is great interest surrounding spacecraft propagation in the LEO environment where atmospheric drag is a dominant non-conservative force. Drag leeches energy and angular momentum from the space object, driving orbit decay, and Hamilton’s equations no longer describe space object evolution. Such a system is non-Hamiltonian. Leaning away from spacecraft applications and towards solar system dynamics, the long term evolution of dust particles situated in heliocentric orbits is influenced by Poynting–Robertson drag effects. This aspect of solar system evolution is of particular interest because various dust structures are dispersed throughout the Solar System including the dust ring trapped in a 1 : 1 mean motion resonance with Earth (Reach 2010). Dust structures, formed from ejected particles from comets and collisional debris from asteroids, are composed of tiny particles on the order of micrometers to millimeters. These particles are highly sensitive under the Poynting–Robertson drag effect. The Poynting–Robertson drag effect is often considered in a relativistic context as it is fundamentally a consequence of general relativity.
It is then of interest to see how such constraints arising from integral invariants and canonical transformations are modified in the presence of non-conservative forces. If the modifications are well understood, we can constrain potential orbits for a larger range of orbital environments.
The structure of this paper is as follows. Section 1 will illustrate spacecraft dynamics and orbits in the context of Hamiltonian phase flow and the usefulness of this dynamical approach. Section 2 will introduce the classical Hamiltonian constraints upon uncertainty distributions and how they will evolve in a general, dynamical system. This section will also be used to illustrate how these constraints are preserved in Hamiltonian systems. Section 3 will present dynamical environments which are influenced by non-conservative forces and show how the modified Hamiltonian constraints yield strong agreement with conventional propagation techniques. Section 4 will discuss the implications of such results for constraining populations of orbiting objects.
2 Phase space structures and spacecraft applications
We often refer to phase space as a higher dimensional space spanned by the canonical coordinates of a given particle or a spacecraft. Canonical coordinates are intrinsic to particles subjected to a particular dynamical environment, and describe how such particles will evolve. As a simple example, one set of canonical coordinates is \(x,\;y,\;z,\;p_x,\;p_y,\;p_z\) indicating the Cartesian components of position and linear momentum respectively. The particle canonical coordinate set is then given by \({\textbf {x}}=({\textbf {q}},{\textbf {p}}) = (x,y,z,p_x,p_y,p_z)\). Given a full set of canonical coordinates at any given time, all instantaneous dynamical aspects of the particle can be computed including its position, velocity, energy, and momentum. The canonical coordinates of a particle obey differential equations set by the dynamical environment. Solving these differential equations will yield trajectories in phase space. At any time, the particle’s canonical coordinates will reside on this 1-dimensional trajectory in phase space. If the initial canonical coordinates for a spacecraft are chosen, a corresponding “phase” trajectory will intersect this initial point in phase space. This trajectory indicates all past and future canonical coordinates the spacecraft will attain. It is common to refer to this particle evolving or stepping along its trajectory as a phase point.
The advantage of viewing spacecraft trajectories in the context of canonical coordinate evolution is it allows one to develop a very intuitive, macroscopic nature for the dynamics. If we consider a collection of spacecraft that all start with slightly perturbed initial conditions around some nominal initial condition, the collection will appear as some “cloud” of phase points in phase space. Each phase point will evolve according to its own trajectory which may diverge or converge to each other. If we generate enough initial phase points around the nominal, we effectively paint in the “cloud” until it is a continuous distribution in phase space. This distribution will then evolve as a parcel of fluid under the dynamics. We can refer to this as a continuous, closed phase space distribution taken about some nominal state. Canonical formalism allows us to treat dynamical motion through phase space as 2n-dimensional fluid flow. Generally, propagating a large number of states within this distribution can provide insight into the global behaviour of the dynamical flow. However, it no longer becomes necessary to numerically propagate a large number of phase points discretizing this distribution because existing Hamiltonian constraints enforce certain macroscopic properties of the distribution.
We are interested in the dynamical evolution of macroscopic structures over individual phase points because we can only be confident to some degree that our particle will reside in a macroscopic region of phase space. Due to small perturbations to our dynamical models and errors in our measurements, spacecraft trajectories become stochastic in nature. A set of spacecraft canonical coordinates, often referred to as the spacecraft state, is described with a probability density function (pdf). As with any pdf, one can integrate over a region of the phase space to obtain a probability or likelihood of the state residing in this region. For Gaussian statistics, it is standard to construct n-\(\sigma \) uncertainty ellipses or ellipsoids in higher dimensions from the covariance matrix. The true state will lie within the ellipse or ellipsoid within some level of confidence. With canonical coordinates, a pdf can be defined for a spacecraft canonical state and a region in phase space can be traced out that returns a certain probability threshold when integrated over. This region then corresponds to a macroscopic phase space structure which will obey constraints alluded to earlier.
While these constraints and how they are applied in aerospace applications may be rather vague, it may be useful to overview their impact on a standard aerospace application. An uncertainty region associated with orbital debris state can be constructed from optical measurements. The constraints on this phase space structure can limit its evolution. This may then be able to inform operators where to point an optical telescope on the celestial sphere for the greatest chance of detection. This can also identify active satellites which cross through the uncertainty region at a future time and be at risk of collision.
3 Constraints on uncertainty distributions
We will examine the following constraints on uncertainty structures or distributions evolving in 2n-dimensional phase space: conservation of probability and the Poincaré-Cartan integral invariants.
3.1 Hamiltonian formalism
Our 2n-dimensional phase space is spanned by coordinates \((q_1,q_2, \ldots ,q_n,p_1,p_2, \ldots ,p_n)\) which denote the canonical position and momentum coordinates of a given particle evolving in phase space. We refer to a particle or phase point by the vector \({\textbf {x}} = ({\textbf {q}}, {\textbf {p}}) = (q_1,q_2, \ldots ,q_n,p_1,p_2, \ldots ,p_n)\); \({\textbf {x}}\in {\mathbb {R}}^{2n}\) where its evolution \({\textbf {x}}(t)\) will trace out a trajectory in phase space. Higher dimensional structures such as hyperplanes and hyperspheres are denoted as \(\beta _{2k}\) where 2k denotes the dimensionality of the structure. Every structure can be thought of as an ensemble of phase points, all glued together. One can then associate the dynamical mapping of an entire structure, \(\beta _{2k}(t)\), with dynamical mapping of the phase points that compose it.
Although canonical transformations and symplectomorphisms are quite general concepts which can refer to coordinate transformations between various canonical coordinate sets or to dynamical mapping, for the purposes of spacecraft applications, we will consider canonical transformation synonymous with dynamical mapping by the solution flow of the Hamiltonian system,
For a given dynamical system, the canonical coordinates of a phase particle evolve according to Hamilton’s equations. This is the trademark feature of canonical coordinates separating them from other types of coordinates. The general form of Hamilton’s equations (Tveter 1994) accounting for both generalized, conservative and non-conservative forces can be shown to be,
H is the system Hamiltonian constructed from conservative forces, \({\textbf {F}}\) is the sum of non-conservative forces acting on the particle, and \({\textbf {r}}\) is the particle’s position vector. For a system to be Hamiltonian, the canonical coordinate evolution must obey,
which implies Hamiltonian systems have no non-conservative forces. When a system is classified as Hamiltonian, an important result can be derived using infinitesimal perturbations to the identity canonical transformation and the property of superposition of canonical transformations (Goldstein et al. 2002). It can be shown a dynamical mapping of a particular region in 2n-dimensional phase space is equivalent to some canonical transformation of the region.
3.2 Conservation of probability
Spacecraft states are rather unknown due to uncertainties embedded in both measurements and dynamical models. These states are characterized by probability density functions that describe the likelihood of the spacecraft state residing in a certain region of phase space. Spacecraft probability density functions generally evolve according to the Fokker-Planck equation, a partial differential equation. The Fokker-Planck equation is defined as (Maybeck 1982),
where \(\rho ({\textbf {x}},t)\) is the probability density function, \({\textbf {f}}({\textbf {x}},t)\) refers to the system dynamics such that \((\dot{{\textbf {q}}},\dot{{\textbf {p}}}) = \dot{{\textbf {x}}} = {\textbf {f}}({\textbf {x}},t)\), and matrices G and Q characterize the diffusion in the stochastic component of the system dynamics. For spacecraft applications, deterministic effects dominate the dynamics even under the influence of non-conservative forces. Hence, it is often a valid approximation to discard the stochastic, diffusion terms in the Fokker-Planck equation. Then, the Fokker-Planck equation, Eq. (4), reduces to,
Macroscopic structures in phase space are of interest for spacecraft uncertainty applications. In particular, we will examine the time evolution of scalar quantities integrated over full dimensional phase structures \(\beta _{2n}\),
where \(M({\textbf {q}},{\textbf {p}},t)\) is some non-autonomous, continuous, scalar function of state. For \(I_M(t)\) to be dynamically invariant, \(M({\textbf {q}},{\textbf {p}},t)\) must obey the following partial differential equation (Greenwood 2000),
Note that no restrictions apart from differentiability have been placed on \({\textbf {f}}({\textbf {x}},t)\). Setting \(\rho ({\textbf {x}},t) = M({\textbf {q}},{\textbf {p}},t)\), \(\rho ({\textbf {x}},t)\) satisfies the invariant condition outlined above, Eq. (7), because it satisfies the Fokker-Planck differential equation, Eq. (5). Therefore,
is invariant to dynamical mapping both with and without non-conservative forces acting. However, this integral over a region of phase space \(\beta _{2n}(t)\) corresponds to the probability P the spacecraft state will reside in \(\beta _{2n}(t)\). The implication is that the probability a spacecraft resides in some initial confidence region in phase space is constant as this region is mapped forward in time. For example, one can outline an ellipsoid in phase space from an initial covariance matrix for spacecraft state assuming Gaussian statistics. The spacecraft state has some corresponding probability to reside in this ellipsoid. As the ellipsoid is dynamically mapped forward in time, the probability to remain inside the mapped ellipsoid is unchanged.
Regarding conservation of probability on lower dimensional structures \(\beta _{2k}(t), k < n\), one can always treat these structures as an infinitesimal limit of a full dimensional structure \(\beta _{2n}(t)\). Hence, probability is conserved on lower dimensional uncertainty structures as well.
3.3 First order Poincaré-Cartan integral invariant
The first order Poincaré-Cartan integral invariant constrains 2-dimensional structures or subvolumes which may extend into the full 2n-dimensional phase space. To continue, we must first introduce the concepts of exterior forms, differential forms, and exterior products (Arnold 1989).
Definition 1
An exterior form of degree 1, otherwise known as a 1-form, is a linear function \(\omega \): \({\mathbb {R}}^n \longrightarrow {\mathbb {R}}\) where
A 1-form must always take a vector as some input and return a corresponding scalar extracted from the features of the vector. A straightforward example of a 1-form is the coordinate 1-form which extracts a particular coordinate from an input vector. For example, consider the field \({\mathbb {R}}^3\) spanned by the coordinates (x, y, z) and let \(\zeta = (\zeta _x, \zeta _y, \zeta _z) \in {\mathbb {R}}^3\), we can then denote a 1-form as x() where \(x(\zeta ) = \zeta _x\). It is clear that this 1-form can be written as \(x(\zeta ) = \zeta \cdot {\hat{x}}\) for any arbitrary vector \(\zeta \). There are many other 1-forms which extract different properties from vectors but we will be focused primarily on 1-forms which extract coordinates from our higher dimensional phase space. 2-forms are defined in a similar manner with some additional properties.
Definition 2
An exterior form of degree 2, also referred to as a 2-form, is a function on pairs of vectors \(\omega ^2\): \({\mathbb {R}}^n \times {\mathbb {R}}^n \longrightarrow {\mathbb {R}}\), which is bilinear and skew symmetric:
Similarly, 2-forms take in two vectors and return a scalar dependent on the properties of the two vectors. We will examine the exterior product which is a 2-form that combines two 1-forms in a simple relation.
Definition 3
The value of the exterior product \(\omega _1 \wedge \omega _2\) on the pair of vectors \(\zeta _1, \zeta _2 \in {\mathbb {R}}^n\) is the oriented area of the image of the parallelogram on the \(\omega _1,\omega _2\)-plane where:
A common application of this 2-form exterior product is computing the signed area between two vectors in \({\mathbb {R}}^2\). Let \(\zeta _1, \zeta _2 \in {\mathbb {R}}^2\) where \(\zeta _1 = (\zeta _{1,x}, \zeta _{1,y})\) and \(\zeta _2 = (\zeta _{2,x}, \zeta _{2,y})\). Then \((x \wedge y)(\zeta _1, \zeta _2) = \zeta _{1,x}\zeta _{2,y} - \zeta _{1,y}\zeta _{2,x}\) which is the signed area of the parallelogram spanned by the two vectors. From Definition 3, we have a very useful property: \((\omega _1 \wedge \omega _1)(\zeta _1, \zeta _2) = 0\).
Consider a higher dimensional space with \(x_1, \ldots , x_n\) as coordinates. \(\zeta (t) = (x_1(t), \ldots ,x_n(t))\) denotes a curve in this higher dimensional space which is parameterized by a single parameter t. A general family of curves traces out a manifold. The tangent space to this manifold contains tangent vectors to the curves outlining such a manifold. At a given point in this higher dimensional space, the tangent vector to a given curve \(\zeta (t)\) is \(\gamma (t) = (\dot{x}_1(t), \ldots , \dot{x}_n(t))\).
Definition 4
A differential form of degree 1, or a differential 1-form, on a manifold M is a smooth map
of the tangent bundle of M to the line, linear on each tangent space \(TM_x\).
An example may clarify the definition. Going to \({\mathbb {R}}^2\), suppose there exists a plane curve defined as \(\zeta (t) = (x(t), y(t)) = (\cos (t), \sin (t))\) with \(\gamma (t) = (-\sin (t), \cos (t))\). The differential 1-forms dx(), dy() produce \(dx(\gamma (t)) = -\sin (t)\) and \(dy(\gamma (t)) = \cos (t)\) which extract coordinates from the tangent vector of \(\zeta (t)\).
We will now introduce some notation that will be utilized in the following derivations. We will often write quantities such as \({\textbf {p}} \wedge {\textbf {q}}\). This may be initially confusing as it appears to be the exterior product between two vectors which is not well-defined given previous definitions. While exterior products between vectors are well-defined, referred to as bivectors, the notation \({\textbf {p}} \wedge {\textbf {q}}\) is shorthand for a summation of 2-form exterior products composed of coordinate 1-forms across an even-dimensional phase space. We will define the following notation which will be used throughout the rest of this work,
Here p, q are not vectors but rather a collection of coordinate 1-forms which must extract scalar \(p_i\), \(q_i\) coordinates when provided input vectors as seen above. Using this notation, we can define two interesting exterior products. Let \(c \in {\mathbb {C}}\), then \(c p_i \wedge q_i\) is a valid exterior product as it is composed of two valid 1-forms obeying Definition 1. Then, we can define \(c {\textbf {p}} \wedge {\textbf {q}}\) as a valid shorthand. More importantly, let \(A \in {\mathbb {C}}^{n\times n}\) be a constant matrix where \(a_{ij}\) are its elements. Then, \(a_{ij}p_i \wedge q_i\) is a valid exterior product. Further, \((a_{i1}p_1 +\cdots +a_{in}p_n) \wedge q_i\) is a valid exterior product as it is composed of two valid 1-forms obeying Definition 1. Finally, \(\sum _{i=1}^n(a_{i1}p_1 + \cdots +a_{in}p_n) \wedge q_i\) is a valid sum of exterior products of coordinate 1-forms. Thus, in our shorthand notation, we can define
as a valid exterior product of 1-forms satisfying Definition 1.
Our 2n-dimensional phase space is generally taken to be \({\mathbb {R}}^{2n}\) with coordinates \(\{q_i,p_i\}_{i=1}^n\) denoting the canonical position and momentum. Extending from particles to structures in phase space, we can project such structures onto the n \(p_i\)-\(q_i\) or ‘symplectic’ planes of this phase space.
The first order Poincaré-Cartan integral invariant, with the use of multi-dimensional Stokes’ theorem, can be formulated as,
which states the sum of the projected signed area of a 2-dimensional phase space structure onto each symplectic plane is conserved under the dynamical flow of a Hamiltonian system. \(P_i()\) is the projection operator onto the ith symplectic plane. For an infinitesimal structure, we have,
where d() is a small displacement of the canonical coordinates in the tangent space, and \(\bar{C}\) is some constant. To study macroscopic behaviour of phase structures, we will examine the collective behaviour of the infinitesimal phase structure composing this macroscopic structure. For our purposes, we will then relabel \(\hbox {d}I(t)\longrightarrow I(t)\) where from now on it will be understood that I(t) refers to an infinitesimal structure.
It is rather unclear what \(\zeta _1(t),\zeta _2(t)\) are when evaluating the above expression. To determine this, it is useful to discuss the notion of the integral of a 1-form along a path. Let this 1-form be denoted as \(\omega ^1\). The integral of the form \(\omega ^1\) on the path \(\gamma \) is defined as a limit of Riemann sums. Every Riemann sum consists of the values of the form \(\omega ^1\) on some tangent vectors \(\zeta _i\). Let this path \(\gamma \) be parametrized by a single parameter \(s \in [0,1]\). Then one can increment this path into discretized segments \(\Delta _i: s_i \le s \le s_{i+1}\),
The tangent vectors \(\zeta _i\) can be constructed in the following way. \(\Delta _i\) can equivalently be considered as a tangent vector to the s-axis (pointing along it from \(s_i\) to \(s_{i+1}\)). \(\zeta _i\) is then defined as the image of the tangent vector \(\Delta _i\) in the tangent space at the point \(\gamma (s_i)\):
Hence the \(\zeta _i\) fed into the 1-form when integrating the 1-form over the parametrized path \(\gamma (s)\) are a consecutive, series of tangent vectors stemming from discrete points along the path. We can apply the same principles for the integration of a 2-form (exterior product) over some sub-dimensional surface. One can take our 2-dimensional surface \(\beta _2(t)\) and discretize it into a series of grid points. At a given node, there exists tangent vectors \(\zeta _{1}, \zeta _{2}\) which span the parallelogram connecting adjacent grid nodes in opposite corners. Therefore, the vectors input into the integrand of the first order Poincaré-Cartan integral invariant correspond to infinitesimal tangent vectors spanning the infinitesimal parallelograms across the entire surface \(\beta _2(t)\). For infinitesimal applications to a single parallelogram, we will denote these two tangent vectors to the surface as \(\zeta , \eta \in {\mathbb {R}}^{2n}\). Then, we have
For this derivation, we will apply the total time derivative operator to determine \(\dot{I}(t)\) and, hence, an “Euler” update for I(t) in the presence of non-conservative forces. To begin this derivation, we must first address the total time derivative operator applied to exterior and differential forms. In other words, is \(\frac{\hbox {d}}{\hbox {d}t}({\textbf {dp}}())\) a valid statement? If so, what is its interpretation?
This may seem rather challenging at first as the time derivative operator applied to an operator (an exterior or differential form) may not be easily defined. The key point to remember is that, when taking the time derivative of the 1st Order Poincaré-Cartan Integral Invariant, Eq. (12), we are not taking time derivatives of operators, we are taking time derivatives of operators acting on input vectors. In other words, the 1st Order Poincaré-Cartan Integral Invariant returns a scalar and it is this scalar we are taking the time derivative of. Since we are working with coordinate exterior products and differential 1-forms, we are simply taking the time derivatives of the coordinates and coordinate displacements of vectors which is a completely valid procedure. For example, if we have a particle tracing out a trajectory in 3D space, at a given moment in time it has a given position and velocity vector. The x-component of the velocity vector is \(\dot{x}=\frac{\hbox {d}}{\hbox {d}t}(x)=\frac{\hbox {d}}{\hbox {d}t}(x({\textbf {r}}))\) where x() is a coordinate 1-form which extracts the x-component of an input vector and \({\textbf {r}}\) is the particle’s position vector. When examining the time variation of this quantity, it is important to recognize that the right hand side is a scalar quantity and hence time differentiation can be applied to this exterior product. We can introduce the shorthand notation \(\dot{\omega }()\) to denote \(\frac{\hbox {d}}{\hbox {d}t}(\omega (\zeta ))\). It is understood that wherever \(\dot{\omega }()\) is written, it is implied that it is evaluated on a vector \(\zeta \). We can then write
following our previously defined notation. Lastly, we will assume the dynamics of this higher dimensional phase space obey the general Hamilton’s equations, also seen in Eq. (2),
When time differentiating I(t), we use the product rule for exterior forms above. Recall the following exterior and differential forms are evaluated on parallelogram spanning vectors \(\zeta , \eta \),
We introduce a rather trivial identity from calculus of variations
where \(\delta ()\) extracts the variation of \(\frac{dx(\zeta )}{dt}\) along \(\zeta \). Applying this identity, we arrive at
One can show through exterior product algebra that \(H_{qq}{} {\textbf {dq}} \wedge {\textbf {dq}}=0;\;H_{pp}{} {\textbf {dp}} \wedge {\textbf {dp}}=0\) and \(H_{qp}{} {\textbf {dp}} \wedge {\textbf {dq}}={\textbf {dp}} \wedge H_{pq}{} {\textbf {dq}}\) with the following two exterior product propositions.
Proposition 1
Let \({\textbf {dx}} = [\hbox {d}x_1, \hbox {d}x_2, \hbox {d}x_3]\) where \(dx_i\) are differential 1-forms and \(x_i\) are the coordinates of \({\mathbb {R}}^3\). Let \(A \in {\mathbb {R}}^{3\times 3}\) where A is a symmetric matrix and, thus, its elements are given by \(a_{ij} = a_{ji}\). Then, \(A{\textbf {dx}} \wedge {\textbf {dx}}=0\).
Proof
\(\square \)
Proposition 2
Let \({\textbf {dx}} = [dx_1, dx_2, dx_3], {\textbf {dy}} = [dy_1, dy_2, dy_3]\) where \(dx_i, dy_i\) are differential 1-forms and \(\{x_i,y_i\}\) are the coordinates of \({\mathbb {R}}^6\) and \(A \in {\mathbb {R}}^{3\times 3}\) where A is an arbitrary matrix. Then, \(A{\textbf {dx}} \wedge {\textbf {dy}}={\textbf {dx}} \wedge A^T{\textbf {dy}}\).
Proof
\(\square \)
\(H_{qq}\) and \(H_{pp}\) are symmetric matrices and \(H_{pq} = H_{qp}^T\). With these theorems, we have
Note that if our non-conservative terms \({\textbf {Q}}_q,{\textbf {Q}}_p\) are zero, we arrive at the Hamiltonian result of conservation: \(I(t+\Delta t) = I(t)\).
3.4 Second order Poincaré-Cartan integral invariant
Following the first order invariant, this second order integral invariant constraints 4-dimensional phase space structures. We define the second order invariant in a similar manner to the first order invariant,
for 2n-dimensional phase space. Here, we employ a 4-form exterior product, acting on four structure spanning vectors \(\zeta _1,\zeta _2,\zeta _3,\zeta _4\). This 4-form is composed of differential, coordinate 1-forms. \(S_{ij}\) refers to a 4-dimensional structure spanned by the two symplectic planes (\(p_i,q_i\) and \(p_j,q_j\)). In the presence of non-conservative forces, the second order integral invariant evolves according to the following an ordinary differential equation. However, we have yet to find a concise expression of this differential equation. So, we omit the study of this invariant for simplicity.
3.5 Third order Poincaré-Cartan integral invariant
In our standard 6-dimensional phase space, the third order invariant is also referred to as the maximum order invariant as it concerns phase space structures of the maximum dimension permitted by the phase space, i.e, the dimension of the phase space itself. One may recognize that these phase space structures are full dimensional phase volumes. This ties the integral invariant to Liouville’s Theorem.
Specifically, Liouville’s theorem, also known as the maximum order Poincaré-Cartan integral invariant, can be stated as follows: Consider a 2n-dimensional bounded structure in 2n-dimensional phase space, \(\beta _{2n}\). Then, under any canonical transformation or symplectomorphism \(\phi \), the 2n-dimensional volume of \(\beta _{2n}\) is conserved. Liouville’s theorem can then be recast as,
A 2n-form exterior product consisting of canonical coordinate 1-forms represents the 2n-dimensional volume of a 2n-dimensional parallelepiped. This is equivalent to an infinitesimal volume element. In other words, \(\hbox {d}p_1 \wedge \hbox {d}q_1 \wedge \cdots \wedge \hbox {d}p_n \wedge \hbox {d}q_n = \hbox {d}q_1 \hbox {d}q_2 \cdots \hbox {d}q_n \hbox {d}p_1 \hbox {d}p_2 \cdots \hbox {d}p_n= \hbox {d}V\).
This volume element is a 2n-form exterior product composed of 2n differential 1-forms that act on 2n vectors spanning a 2n-dimensional parallelepiped. As mentioned earlier, we are interested in infinitesimal phase space structures so, for clarity, we will relabel \(dV(t)\longrightarrow V(t)\). For \(V=dV\),
Substituting in the general form of Hamilton’s equations, Eq. (2), one can show all Hamiltonian terms cancel leaving us with
3.5.1 Higher dimensional fluid approach
Many parallels can be drawn between fluid mechanics and phase flow. Mathematically, bounded 2n-dimensional phase structures are equivalent to parcels of higher dimensional fluid evolving in time. A standard theorem in fluid mechanics, derived from Leibniz Integral Rule, is Reynold’s Transport Theorem which examines volume integrated quantities over fluids with deforming surfaces and changing boundaries. For a phase fluid moving along the dynamical flow with some 2n-dimensional velocity, Reynold’s Transport Theorem states for a scalar function f that:
where \(S(\beta _{2n}(t))\) is the surface enclosing the phase volume, \({\textbf {v}}\) is the velocity of the phase volume moving along the dynamical flow, and \(\hat{{\textbf {n}}}\) is a unit vector normal to each infinitesimal surface element of the phase volume. If we consider the case where \(f=1\), then we have:
where the last equality made use of the Divergence Theorem. The left-hand side of the above equation corresponds to the time rate of change of the volume. So, we arrive at:
This is consistent with our exterior product derivation above, Eq. (25). Substituting in our general Hamilton’s equations, Eq. (2), we then have,
This expression describes the general phase volume evolution under both non-conservative and conservative forces. Note that if our non-conservative forces \({\textbf {Q}}_q, {\textbf {Q}}_p\) are zero, the Hamiltonian result of volume preservation is recovered: \(\frac{\hbox {d}V}{\hbox {d}t}=0\).
4 Data for validation
When examining non-conservative environments to test against the theory developed above, we compute data for our invariant quantity evolution to probe the validity of the above analysis. This data comes from an alternative method to compute the integral invariant which itself has approximations. Hence, we will have the differential equations above and these alternative techniques to compute the same invariant and compare results. For phase volume propagation, we propagate a state covariance matrix alongside our analytical expressions using the state transition matrix (STM) of the dynamical system assuming a linear update. At every time step, the propagated state covariance matrix is used to compute an estimate of the 2n-dimensional phase volume contained within an instantaneous 2n-\(\sigma \) ellipsoid. This ellipsoid is defined by the entries of the covariance matrix. Under a Gaussian statistics assumption, this ellipsoid indicated a confidence region in phase space using the Mahalanobis distance. The volume contained within this ellipsoid is then provided using the determinant of the covariance matrix, \(V(t) = \sqrt{\text {det}(P)}\). One can show this by using generalized spherical coordinates to integrate over the region of the higher dimensional ellipsoid. The linear update for the covariance matrix is given by,
where \(P_k^-\) refers to the state covariance matrix at time \(t_k\) and \(\Phi _{(t_{k+1}, t_{k})}\) refers to the STM mapping a small perturbation to a nominal trajectory in a linear sense forward in time. The STM is integrated alongside the nominal trajectory of the spacecraft according to,
This alternative approach for computing volumes is valid because we are examining infinitesimal volume elements which span very small portions of phase space. With small deviations in phase space, a linear update of a covariance matrix assuming Gaussianity is a good approximation.
For the first order Poincaré-Cartan integral invariant, i.e, total signed area evolution, we consider an infinitesimal area “patch” which we will track in phase space. The patch consists of four phase points propagating together in the close vicinity of each other. The points are localized around some nominal phase trajectory. We treat the phase points at every single time instance as vertices of a higher dimensional polygon embedded in 6-dimensional phase space. We compute the signed areas of the images of this higher dimensional polygon on each of the symplectic planes. Each image itself is a 2D polygon on the \((q_i,p_i)\) symplectic plane.
Assuming each 2D polygon can be classified as a simple polygon where no portions of the polygon overlap or have varying signs of signed area, we determine the area of the projection onto each symplectic plane using the Shoelace formula. This formula takes in the x and y coordinates for the vertices of each polygon and computes the signed area spanned by the vertices of the 2D polygon,
For the four phase points propagating in close vicinity, we store the signed areas \(A_{q_1p_1}\), \(A_{q_2p_2}\), and \(A_{q_3p_3}\) and sum them together to produce a form of alternative data for I(t). Once again, the justification for this alternative method falls from our infinitesimal treatment. With vertices close enough together, the vectors joining the vertices can approximate tangent vectors associated with a phase point residing on the manifold of the flow.
5 Non-conservative environments
5.1 Earth two-body problem: atmospheric drag
The standard expression for the drag force acting on a given spacecraft passing through an atmosphere is given by,
where \(\rho ({\textbf {r}})\) now denotes the atmospheric density as a function of spacecraft altitude, \({\textbf {r}}\) denotes the spacecraft position vector, \(C_D\) denotes the spacecraft drag coefficient, A denotes the cross-sectional area of the spacecraft facing the ram direction, \({\textbf {v}}\) denotes the spacecraft velocity vector.
There are two common canonical coordinate choices for this environment: Cartesian canonical coordinates and Delaunay canonical coordinates. Cartesian canonical coordinates are used because the drag force is easily expressed in them, and the general Hamilton’s equations reduce quite nicely when the elements of the position vector are solely canonical position coordinates. Delaunay coordinates are used because of their simple relations with the Keplerian orbital elements, and, hence, rather simple dynamical evolution. In the unperturbed Two-Body problem, each of the Cartesian canonical coordinates will vary dramatically over a single orbit period. However, five of the six Delaunay canonical coordinates will remain constant while the 6th evolves in a rather straightforward manner following the mean anomaly M. This is because Delaunay canonical coordinates are action-angle coordinates. In a Hamiltonian system, the actions are integrals of motion and the angles evolve linearly in time according to frequencies defined by the actions. The unperturbed, Two-Body problem has two of the three frequencies zero and hence two angles and all three actions are constant. The coordinate sets are defined in the following way,
Cartesian
Delaunay
where relations between the Delaunay elements and the Keplerian orbital elements \((a, e, i, \Omega , \omega , M)\) are illustrated. The conservation of probability result will not be investigated in this environment as we have shown it to hold for both non-conservative and conservative environments, Eq. (8). Starting with Cartesian coordinates, our Hamilton’s equations reduce to,
where \(\bar{H} = -\frac{1}{2\,m}(p_x^2+p_y^2+p_z^2)-\frac{GMm}{\sqrt{x^2+y^2+z^2}}\), m is the mass of the spacecraft, and M is the mass of Earth. We will define \(\mu \equiv GM\). We can also recast the drag force in Cartesian canonical coordinates,
It is important to recall that atmospheric density \(\rho \) is generally a function of Cartesian position. So, we work with two atmospheric density profiles: (i) a simplified, constant density profile, (ii) and an exponential density profile. Using the Mass Spectrometer and Incoherent Scatter radar (MSIS) profile from 2005 sampled from 0 km to 999 km above the United States, a simple exponential profile was fit to model the atmospheric density,
where \(\rho _0\) is some reference atmospheric density at a reference altitude \(h_0\), \(H_s\) is a scale height characterizing the balance of gravitational and thermal effects in the atmosphere, r is the magnitude of the spacecraft’s geocentric position vector, and \(R_E\) is the radius of the Earth. Our constant density profile is just set to be \(\rho ({\textbf {r}}) = \rho _0\).
For both an exponential profile and constant profile, we have,
With Cartesian canonical coordinates, the position dependence of the exponential profile plays no role in the divergence of the phase flow velocity field. We consider the evolution of a microscopic or infinitesimal volume element. If all the infinitesimal volume elements evolve together in a similar way or exhibit the same behaviour, then conclusions can be made about the collective motion of these infinitesimal volumes and hence the macroscopic volume of \(\beta _{2n}\). Then, we have,
where the integral \(\int \rho ({\textbf {r}}) v dt\) can be evaluated in a straightforward manner via numerical integration given the atmospheric density profile and the spacecraft state evolution. With Cartesian canonical coordinates, it is the first order Poincaré-Cartan integral invariant where position dependence in the non-conservative force plays a significant role. For a constant density profile \(\rho ({\textbf {r}})=\rho _0\), the non-conservative forces become \({\textbf {Q}}_q = {\textbf {0}}; {\textbf {Q}}_p ={\textbf {Q}}_p({\textbf {p}})\).
Approximating the first order integral invariant evolution with an Euler update, we can simulate this invariant evolution as,
The signed area update is then, assuming \(\rho (\textbf{r}) = \rho _0\),
For an exponential profile, Eq. (38), we have \({\textbf {Q}}_q = {\textbf {0}}; {\textbf {Q}}_p = {\textbf {Q}}_p({\textbf {q}},{\textbf {p}})\) and our update becomes,
Switching over to Delaunay canonical coordinates, Barrio and Palacian (2003) feature the Gauss’ Planetary Equations for the Delaunay variables. Given the form of \({\textbf {F}}_D\) outlined for atmospheric drag above, we have,
for both constant and exponential density profiles. Note one can derive these equations of motion by using the general modified Hamilton’s equations of motion with the specific energy of the system provided by \(\varepsilon = -\frac{\mu ^2}{L^2}\), the spacecraft position \({\textbf {r}} = {\textbf {r}}(L,G,H,l,g,h)\), and the atmospheric drag \({\textbf {F}}_D\) expressed in Delaunay coordinates. Note that \(\dot{l}=\sqrt{\mu /a^3},\dot{g}=0,\dot{h}=0\) in the Two-Body problem. Looking at Eq. (50), \(\dot{l},\dot{g}\) are clearly different in the presence of atmospheric drag. The resulting terms correspond to \(\mathbf{Q_q}\) in the general Hamilton’s equations where \(\dot{l} = \sqrt{\mu /a^3} +Q_l, \dot{g} = Q_g\).
Following the same procedure as with Cartesian canonical coordinates for the Delaunay canonical coordinates, the volume and signed area updates will be much more tedious to derive analytically. For this reason, a Delaunay variable update for signed area will not be presented due to its algebraic length and complexity. In practice, one will always choose the most convenient canonical coordinate choice for their application which, for us, is the one that the non-conservative force is most simply expressed in. A Delaunay variable update for phase volume will only be presented in the case of a constant density profile for a circular orbit. Proceeding through the algebraic steps of computing the velocity field divergence, one can arrive at,
for a circular orbit. An approximate analytical expression can actually be derived from this differential equation which can describe the orbit-averaged phase volume evolution. This is only for a constant density atmospheric profile \(\rho ({\textbf {r}})=\rho _0\). Under this circular orbit assumption, true anomaly becomes synonymous with mean anomaly. One can show the time-averaged phase volume will obey a differential equation given by the non time-averaged differential equations with all time dependent quantities replaced by time-averaged versions of themselves. We then arrive at,
We rely on a common analytical technique used to solve ordinary differential equations: a series expansion. In particular, we Taylor expand our orbital elements about a reference solution via a small parameter. To first order, \(\bar{a}\) takes the form: \(\bar{a} = \bar{a}_0 + \dot{a}(t_0)(t - t_0)\). Imposing \(\rho ({\textbf {r}})=\rho _0\), we find,
where \(\bar{a}(t_0)\) can be obtained from the atmospheric drag Gauss’ Planetary Equations listed above. We refer to this V(t) expression above, Eq. (53), as the analytical, averaged expression. Similar infinitesimal volume evolution features can be seen between the two canonical coordinate sets. Both Delaunay and Cartesian volume updates indicate \(\frac{\hbox {d}V}{\hbox {d}t}\) will be non-increasing for all time. It is straightforward to see this with Cartesian coordinates. Since \(|2\sin ^2(f)| \le 2\), \(\frac{\hbox {d}V}{\hbox {d}t} \le 0\) for Delaunay coordinates. All volume elements will contract in some sense. It is less trivial to predict the outcome of the signed area evolution,
Our volumes are all initialized with some ellipsoidal region defined by the 1\(\sigma \) uncertainty ellipsoid of the following covariance matrix,
This \(P_0\) provides \(V_0\) via \(V = \sqrt{\text {det}(P)}\) and also supplies the validation data from linear propagation, Eq. (30). For comparison to the validation data, we must initialize some orbits. We consider a small object in the LEO environment with mass of m = 1 kg, A = 12 cm\(^2\), \(C_D\) = 2.2. This small object is initially placed into an orbit defined as (\(a_0,e_0,i_0,\Omega _0,\omega _0,M_0\)) = (6781 km, 0.0007, 0\(^\circ \), 0\(^\circ \), 0\(^\circ \), 0\(^\circ \)). Regarding atmospheric density, we take \(h_0\) = 405 km, \(\rho _0\) = 0.001069 kg/km\(^3\), and H = 47.946372 km. This spacecraft is then integrated for one nominal Two-body orbit period (5.5\(\times \)10\(^3\) s). To showcase evolution on equal footing, Fig. 1 illustrates the phase volume evolution for a constant density atmosphere.
Top: 2n-dimensional phase volume evolution of spacecraft uncertainty distribution subject to Earth Two-Body and constant density atmospheric drag effects (\(\rho (r) = \rho _0\)); Specifically, the change from the initial volume \(V_0\) is shown. This plot shows 30 orbit periods of volume evolution with our analytical and semi-analytical expressions. Bottom: This plot shows a zoom-in on one orbit period of volume evolution with the validation data. Oscillatory behaviour is seen in the Numerical Approx.: Delaunay curve likely due to our artificial hold of zero eccentricity
Since the numerically integrated orbit remains roughly circular, we can place our Delaunay results onto this plot as these results had an artificial zero eccentricity hold.
The Validation: Propagated Covariance curve captures the phase volume evolution using covariance matrix determinants. The Numerical: Cartesian curve shows the phase volume evolution from the numerical integration of \(V(t) = V_0e^{-2\frac{C_DA}{m}\int \rho ({\textbf {r}})v \hbox {d}t}\). The Numerical Approx.: Delaunay curve depicts the phase volume evolution from the numerical integration of \(V(t) = V_0e^{\int \left( \frac{\rho C_D A}{m} \right) \sqrt{\frac{\mu }{a}}(2\sin ^2f - 3) \hbox {d}t}\), which has an artificial \(e=0\) making it simplified approximation. Lastly, the Analytical-Averaged: Delaunay curve illustrates the phase volume evolution from
\(\bar{V}(t) = \bar{V}_0 \cdot \hbox {exp} \left( -\frac{1}{2}\sqrt{\mu }\left( \frac{\rho _0 C_D A}{m}\right) \left( \frac{8}{\dot{a}(t_0)}\sqrt{\bar{a}_0 + \dot{a}(t_0)(t - t_0)}\right) \right) \).
In general, one should still expect the Delaunay evolution and the Cartesian evolution to be consistent. Dissipative forces will destroy the canonical transformation property of dynamical mapping, but an instantaneous transformation between two canonical coordinate sets should still be canonical. Then, dynamical mapping in one set of canonical coordinates will be consistent with dynamical mapping in another set. From Fig. 1, there is strong agreement between the validation data from covariance propagation, the Cartesian results, and the Delaunay results. However, it is important to note that covariance propagation becomes numerically unstable due to varying orders of magnitude in the covariance matrix entries. This is why we do not plot past one nominal orbit period with the validation data. Figure 2 depicts the exponential atmospheric density profile volume evolution of both the validation data and \(V(t) = V_0e^{-2\frac{C_DA}{m}\int \rho ({\textbf {r}})v\hbox {d}t}\).
Top: 2n-dimensional phase volume evolution of spacecraft uncertainty distribution subject to Earth Two-Body and exponential density atmospheric drag effects (\(\rho (r) = \rho (r)\)); Specifically, the change from the initial volume \(V_0\) is shown. This plot shows 30 orbit periods of volume evolution with our semi-analytical expression. Bottom: This plot captures 1 orbit period of volume evolution with the validation data. Propagating past this timescale causes the propagated covariance to break down
Once again, strong agreement is seen until the validation data goes unstable.
For signed area analysis, we take the orbit (\(a_0,e_0,i_0,\Omega _0,\omega _0,M_0\)) from above and center our infinitesimal patch about this orbit. To construct this patch, we work with the Delaunay variables where the 2-dimensional patch is initially parallel to the L-l symplectic plane. That is, we use the phase point (\(a_0,e_0,i_0,\Omega _0,\omega _0,M_0\)) to generate (\(L_0,G_0,H_0,l_0,g_0,h_0\)), via Eq. (35). We then generate four additional phase points located at (\(L_0\pm \Delta L,G_0,H_0,l_0\pm \Delta l,g_0,h_0\)), where \(\Delta l = 0.005\) and \(\Delta L =0.5,2.5\) for exponential and constant drag profiles respectively. These four phase points form the nodes of the polygon approximating the infinitesimal patch, and are used to compute signed areas via the Shoelace formula and approximate the tangent vectors used to evaluate in Eq. (46). Figure 3 captures the signed area evolution for both the constant and exponential density profiles results alongside the validation data. This is shown in the top and bottom panels.
Top Panel: Sum of signed areas of projections of infinitesimal 2D patch projected on the symplectic planes; The patch was subject to Earth Two-Body and constant density atmospheric drag effects, for 5 nominal orbit periods; Bottom Panel: Sum of signed areas of projections of infinitesimal 2D patch projected on the symplectic planes; The patch was subject to Earth Two-Body and exponential density atmospheric drag effects, for 5 nominal orbit periods
Agreement is seen with these signed area evolutions as well. An interesting question is whether the changes to these integral invariants are physically ‘noticeable’ in the presence of these natural, non-conservative forces. For our constant density evolutions, the percent changes from initial values were as follows: Phase Volume Change − 0.0241\({\%}\) decrease; Signed Area Change − 3.05\({\%}\) decrease. These changes are rather insignificant. If we switch over to an exponential atmospheric density profile, then we find that much more drastic changes in the signed area evolution but not the volume: Phase Volume Change − 0.0220\({\%}\) decrease; Signed Area Change − 1138.8\({\%}\) decrease. Perhaps positional dependence in non-conservative forces generally manifest into large changes to lower order integral invariants rather than the higher order invariants.
5.2 Sun Two-Body problem: Poynting–Robertson drag
Under classical assumptions as seen in several heliocentric dust applications, the Poynting–Robertson drag effect can be approximated by a rather elegant expression,
where it is assumed the dust particle is spherical. W is the total power absorbed by the dust particle from incident sunlight on its cross-sectional area, c is the speed of light, \({\textbf {v}}\) is the dust grain orbital velocity, \(\alpha \) is the dust grain absorptivity assumed to be uniform and time-invariant over the grain, R is the spherical radius of the grain, \(L_{\odot }\) is the solar luminosity, m is the mass of the grain, r is the heliocentric separation between the Sun and the grain, and \({\textbf {p}}\) is the grain’s linear momentum.
For this application, our volumes start as an ellipsoidal region defined by the 1\(\sigma \) uncertainty ellipsoid of the following covariance matrix,
We take \(\alpha \) = 0.9 as dust particles can typically be modeled as near perfect blackbodies, \(L_\odot \) = \(3.828\times 10^{26}\) W, R = \(0.5 \upmu \)m, and m = \(1\times 10^{-16}\) kg. Our dust particle orbit is initialized for a nominal, heliocentric, circular orbit at 1 AU in an attempt to mimic a particle in the dust ring trapped in resonance with the Earth. The initial orbit was selected in heliocentric, Cartesian coordinates as \((x_0,y_0,z_0,\dot{x}_0,\dot{y}_0,\dot{z}_0)=(1\;\text {AU},0,0,0,\sqrt{\mu /x_0},0)\). Once again, due to the position dependence in this dissipative force, derivations involving Delaunay canonical coordinates are omitted. Repeating a similar procedure as with atmospheric drag, we study the divergence of the phase velocity field when subject to the Poynting–Robertson drag effect,
For the signed area update, \({\textbf {Q}}_q = {\textbf {0}}; {\textbf {Q}}_p = {\textbf {Q}}_p({\textbf {q}},{\textbf {p}})\) and our update becomes,
When simulating signed area, our infinitesimal patch is initially parallel to the \(p_x\)–x plane with \(\Delta x=5\) m and \(\Delta p_x = 5\times 10^{-13}\) kg m/s. Figures 4 and 5 capture both the infinitesimal volume and signed area evolution with respect to the validation data. The propagated covariance goes unstable rather quickly. However, it does appear to be consistent with the volume produced from our ODEs for at least the first 3 months of integration. The polygon propagation matches the signed area produced from ODEs for a much longer duration.
Top: 2n dimensional phase volume evolution of spacecraft uncertainty distribution subject to Sun Two-Body and Poynting–Robertson drag effects; Specifically, the change from the initial volume \(V_0\) is shown. The plot depicts 30 orbit periods of volume evolution with our semi-analytical expression. Bottom: The plot illustrates a zoom-in on the first 1 orbit period of volume evolution with the validation data. Propagating past this timescale causes the propagated covariance to break down
Sum of signed areas of projections of infinitesimal 2D patch projected on the symplectic planes; The patch was subject to Sun Two-Body and Poynting–Robertson drag effects
Given the \(1/r^4\) dependence in Eq. (63), we expect the magnitude of \(\frac{\hbox {d}I}{\hbox {d}t}\) to increase in time as the dust particle is stripped of its angular momentum and falls into the Sun. All of this is to show the applicability of an integral invariant perspective to a variety of dynamics problems.
6 Discussion
For the perturbations to phase volume evolution, we were able to argue that all microscopic volume elements contract inwards on themselves to some degree as the divergence of the phase velocity field under atmospheric drag and Poynting–Robertson drag is globally non-increasing in phase space. This implies all macroscopic phase volumes \(\beta _{2n}\) composed of microscopic volumes will collapse inwards on themselves. This macroscopic phase volume is often associated with a confidence or uncertainty region defined using some probability density function. Aspects of the spacecraft state will become more certain as it evolves in LEO will increase over time. This is because the probability of the spacecraft residing in \(\beta _{2n}\) is constant, but the volume occupied by \(\beta _{2n}\) falls, as seen in Fig. 2. Note this does not prevent other aspects of the spacecraft state from becoming more uncertain, as a volume collapse may represent a 6-dimensional ball collapsing to a spread out 5-dimensional pancake.
While both atmospheric drag and Poynting–Robertson drag will have collapsing uncertainty distributions, the structure of this collapse varies significantly. We can see this by inspecting the divergence of the velocity fields for both drag forces. Note that the divergence of the velocity field can be considered as an analog of the rate of collapse for the phase volume. With the constant density profile, the rate of collapse only relies on the particle momentum. With an exponential atmospheric density, the rate of collapse depends on the entirety of the particle state. With Poynting–Robertson drag, only particle position places a role in the rate of collapse. This observation is interesting when framed in the following way: sections of \(\beta _{2n}\) with larger momenta will collapse faster when constant atmospheric drag acts on the particle; sections of \(\beta _{2n}\) with smaller position and larger momenta will collapse faster when exponential drag acts; sections of \(\beta _{2n}\) with smaller position will collapse faster when Poynting–Robertson drag acts. Phase structures will collapse slower or faster based on the stretching of the structure along the momenta coordinate axes, position coordinate axes, or both.
While we only propagate volume evolution for a finite number of nominal orbit periods, it is clear that these volumes should asymptotically collapse to zero as time tends to infinity. Since probability is conserved, it may seem a finite probability can be assigned to a single phase point given an infinite amount of time. This would be a contradiction to the notion of a continuous probability density function. It is important to recall that a phase space structure does not need to collapse to a point to have zero 2n-dimensional volume. The structure must only collapse to a sub-volume for a zero volume. In the context of uncertainty, this corresponds to absolute knowledge of a particular aspect of the spacecraft motion whether it be a particular canonical coordinate or some more complicated relationship. The notion of our uncertainty region collapsing to zero volume can be best understood with a practical example. This example will ignore physical effects of atmospheric friction on the spacecraft such as increased temperatures on the spacecraft.
Consider a spacecraft passing through a dense atmosphere surrounding a planet and that we can ignore temperature effects. As angular momentum is stripped from the spacecraft, the spacecraft will asymptote to a rectilinear orbit. Gradually, the spacecraft will fall more radially inwards passing through an exponentially increasing density atmosphere. Independent of its incident velocity, the spacecraft will ultimately approach the terminal velocity dictated by the spacecraft altitude and atmosphere properties. This represents an asymptotic physical collapse of the phase space distribution. Although where in position space the spacecraft may achieve terminal velocity is uncertain, its velocity direction becomes more and more certain as time passes. Further, the terminal velocity magnitude becomes constrained by the position of the spacecraft. For Cartesian velocity components, this enforces a direction constraint on the components as well as a magnitude constraint. In other words, \(\dot{r},\dot{\theta },\dot{\phi }\) are all constrained. This would reduce a 6-dimensional uncertainty distribution to a 3-dimensional “surface” asymptotically. The distribution collapses to a subvolume and one would expect a zero volume despite imperfect knowledge about the full spacecraft state, particularly its position. Given infinite time, the spacecraft would lie on a rectilinear orbit and the collapse to a 3-dimensional volume would be complete.
For the first order Poincaré-Cartan integral invariant, the signed areas alone do not provide enough insight into the geometrical evolution of the uncertainty distribution. Rather, one requires both the geometry and scalar areas to understand the physical constraints imposed on spacecraft uncertainty evolution. We will present the discussion of the constant atmospheric drag scenario. It is for this type of analysis where the use of Delaunay canonical coordinates becomes quite practical. Figure 6 presents a ‘toy’ uncertainty distribution that has been evolved in a constant drag, Two-Body environment. A simple 2-dimensional ring initially aligned in the L–l plane was generated and modeled as a multi-vertex polygon. This ring then evolved and became deformed in both shape and orientation as evident by Fig. 6.
2-dimensional ring initially aligned with the L-l symplectic plane (Top, \(I(t)=78.54\)) and evolved under Earth Two-Body gravity and exponential density atmospheric drag (Bottom, \(I(t)\approx 38\))
Cartoon in each quadrant shows one of the four conditions. In all cases, we can see this integral invariant prevents perfect knowledge in more than one symplectic plane
The ring becomes heavily contorted. This implies a strong correlation between L and l. If we know \(l\in [l_0,l_1]\), then the possible values of L fall in a small range. This illustrates a constraint derived directly from the geometry of our uncertainty distribution.
6.1 Heisenberg’s uncertainty principle analogy
For Hamiltonian systems, the sum of signed areas does not constrain any area of a symplectic plane. Area can be traded between the symplectic planes as long as the sum is preserved. Each phase point on the ring above initially has the same G, g, H, h coordinates. Since \(\dot{h} = 0\), then the signed area of the H-h symplectic plane \(A_{Hh}\) will be zero. With \(A_{Hh}=0\) for all time, signed area is simply traded between the L-l and G-g planes. In the absence of drag, we would then have,
Let \(\sigma _{ii}\) denote the ‘spread’ along the ith coordinate axis in the i–j symplectic plane for the projection of the uncertainty distribution. This can also be thought of as the “uncertainty” in the ith coordinate. If \(\sigma _{ii} = 0\), then \(A_{ij}=0\). The following four conditions must hold for all time,
If absolute knowledge in one canonical coordinate is obtained, then the signed area of the projection of the uncertainty distribution into the symplectic plane containing the canonical coordinate is zero. This is because ith coordinate will fall along a line corresponding to this known coordinate value. In this case, the area on the remaining symplectic plane cannot collapse to zero to preserve the sum of signed areas. This represents an analog to the famous Heisenberg’s Uncertainty Principle in quantum mechanics where the four conditions describe how perfect knowledge in one canonical coordinate prevents perfect knowledge in any other canonical pair. Once again, these conclusions only hold for Hamiltonian systems. Figure 7 visually shows these four conditions and their implications on the geometry for the uncertainty region.
Transferring to a non-conservative system, we have \(A_{Ll}(t)+A_{Gg}(t) = I(t)\) where I(t) is a decreasing function under the influence of atmospheric drag. As long as I(t) is non-zero, the above four Hamiltonian conditions will hold for non-Hamiltonian systems. If we are certain (\(\sigma _{ii}=0\)) in a subcomponent of the spacecraft state, the remaining subcomponents are more tightly constrained with a decreasing I(t) by the previous arguments. If I(t) converges to zero, the uncertainty of two spacecraft canonical coordinates belonging to the two symplectic planes can simultaneously approach zero. Hence, our confidence could possibly increase. The other case could be that the two areas grow with opposite orientation to cancel each out leading to lower confidence in both symplectic coordinate pairs.
If I(t) increases, then absolute knowledge in one coordinate implies a more poorly constrained canonical pair in the other symplectic plane. In terms of aerospace applications, we seek sum of signed area evolutions that are decreasing in time.
7 Conclusion
In this work, we have addressed the modification of classical Hamiltonian phase space distribution constraints in non-Hamiltonian systems. We first outlined Hamiltonian phase flow in the context of stochastic spacecraft trajectories, and presented motivation for why a macroscopic, phase flow approach to propagating spacecraft uncertainty can lead to important constraints that are not available from standard uncertainty propagation techniques. Three Hamiltonian constraints were primarily examined: (i) Conservation of probability; (ii) the first order Poincaré-Cartan integral invariant; (iii) the third order Poincaré-Cartan integral invariant. Each constraint was re-derived in the presence of a dissipative force leading to differential equations describing the evolution of these integral invariants. Alongside validation data, the differential equations were integrated both numerically and analytically when possible to yield integral invariant evolutions. The modified constraints in terms of integral invariant evolutions were examined in the context of an atmospheric drag, Two-Body environment and a Poynting–Robertson Drag, Two-Body environment. Strong agreement was seen between the validation data and evolutions from the derived differential equations. The behaviour seen in the integral invariant evolutions was then connected to physical implications for our confidence in spacecraft state for both full and sub-dimensional uncertainty structures. For future work, we plan to design control laws to produce non-conservative forces with desired behaviour. The nature of volume evolution, for instance, is completely determined by the divergence of the non-conservative force in 2n-dimensional phase space. Hence, control laws will restrict geometric aspects of our confidence regions in phase spaces. Perhaps one may wish for the volume of the uncertainty region to expand or contract depending on what information is desired of a spacecraft at an instant in time.
Data availability
The software generated and data analyzed in this article will be shared on reasonable request to the corresponding author.
References
Arnold, V.I.: Mathematical Methods of Classical Mechanics. Springer, New York (1989)
Barrio, R., Palacian, J.: High-order averaging of eccentric artificial satellites perturbed by the Earth’potential and air-drag terms. In: Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, Vol. 459, (2003). https://doi.org/10.1098/rspa.2002.1089
Boone, S., McMahon, J.: Directional state transition tensors for capturing dominant nonlinear dynamical effects, vol. 08 (2021)
Goldstein, H., Poole, C., Safko, J.: Classical Mechanics. Addison Wesley, Boston (2002)
Greenwood, D.T.: Classical Dynamics. Dover Publications, New York (2000)
Maruskin, J.: On the dynamical propagation of subvolumes and on the geometry and variational principles of nonholonomic systems, vol. 01 (2008)
Maruskin, J., Scheeres, D., Bloch, A.: Dynamics of symplectic SubVolumes. SIAM J. Appl. Dyn. Syst. 8, 10 (2007). https://doi.org/10.1137/070697938
Maybeck, P.S.: Stochastic Models, Estimation and Control, vol. 2, pp. 159–271. Academic Press, New York, NY (1982)
McDowell, J.C.: The low earth orbit satellite population and impacts of the SpaceX Starlink constellation. Astrophys. J. Lett. 892, L36 (2020). https://doi.org/10.3847/2041-8213/ab8016
Milani, A., Gronchi, G.F., Vitturi, D., Michieli, M., Knezevic, Z.: Orbit determination with very short arcs. I - Admissible regions. Celest. Mech. Dyn. Astron. 90(07), 57–85 (2004). https://doi.org/10.1007/s10569-004-6593-5
Park, R., Scheeres, D.: Nonlinear mapping of Gaussian statistics: theory and applications to spacecraft trajectory design. J. Guid. Control Dyn. 29(11), 1367–1375 (2006). https://doi.org/10.2514/1.20177
Park, R., Scheeres, D.: Nonlinear semi-analytic methods for trajectory estimation. J. Guid. Control Dyn. 30(11), 1668–1676 (2007). https://doi.org/10.2514/1.29106
Reach, W.: Structure of the Earth’s circumsolar dust ring. Icarus 209(10), 848–850 (2010). https://doi.org/10.1016/j.icarus.2010.06.034
Scheeres, D., Gosson, M.D., Maruskin, J.: Applications of symplectic topology to orbit uncertainty and spacecraft navigation. J. Astronaut. Sci. 59(03), 63–83 (2014). https://doi.org/10.1007/s40295-013-0006-5
Scheeres, D., Hsiao, F.-Y., Park, R., Villac, B., Maruskin, J.: Fundamental limits on spacecraft orbit uncertainty and distribution propagation. J. Astronaut. Sci. 54, 12 (2012). https://doi.org/10.1007/BF03256503
Tommei, G., Milani, A., Rossi, A.: Orbit determination of space debris: admissible regions. Celest. Mech. Dyn. Astron. 97, 04 (2007). https://doi.org/10.1007/s10569-007-9065-x
Tveter, F.T.: Hamilton’s equations of motion for non-conservative systems. Celest. Mech. Dyn. Astron. 60, 409–419 (1994). https://doi.org/10.1007/BF00692025
Acknowledgements
Oliver Boodram and Daniel Scheeres are supported by AFOSR (Air Force Office of Scientific Research) Grant FA9550-21-1-0332.
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
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://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Boodram, O., Scheeres, D. Constrained evolution of Hamiltonian phase space distributions in the presence of natural, non-conservative forces. Celest Mech Dyn Astron 136, 1 (2024). https://doi.org/10.1007/s10569-023-10172-1
Received:
Revised:
Accepted:
Published:
DOI: https://doi.org/10.1007/s10569-023-10172-1