Introduction

Astrodynamics and celestial mechanics have been since their birth focused on point dynamics, that is, how a single satellite (or planet) moves about a single attracting body. For unperturbed motion particularly, this point dynamics problem is well understood; not just the initial value problem of orbit propagation, but the boundary value, or Lambert, problem is as well. Here, the dynamics of unperturbed motion is examined in a new way: how a cloud of particles from a point source evolves as a density distribution under the same unperturbed two-body motion. This is significantly more difficult to compute than the point dynamics case, but the rewards are potentially very large.

There are two kinds of such clouds: real clouds, such as from a fragmentation, and virtual clouds, representing a probability distribution, as from uncertain orbits or forces. When an instantaneous fragmentation happens on orbit, a single object suddenly becomes multiple objects, and the orbits of each of these pieces will be different from the source orbit. In fact, the fragmentation adds a random impulsive delta-v to each piece, and analysis may proceed by treating each imparted delta-v as an equivalent maneuver.

Researchers have employed various techniques to study the evolution of the cloud of fragments. Jehn [12] divided the post-fragmentation time interval into four phases, labeled “A,” “B,” “C,” and “D.” In phase A, the cloud of fragments is a pulsating ellipsoid, being still confined relatively close to the source orbit. The pulsation is seen to occur every orbital period; the extent of the cloud decreases and hence its density increases as it returns to the fragmentation point. In phase B, the shape of the cloud is a torus with the fragmentation point a pinch point. This may occur as soon as a few hours after fragmentation, where the spreading of fragments into the toroidal ring is a velocity change component in the in-track direction that changes the orbital periods of the fragments. In phase C, perturbations, notably the J2 geopotential perturbation, cause the nodal regression which smears out the torus into a spherical shell truncated at the polar caps. Finally in phase D, if the orbit is low enough, the evenly distributed band loses population to atmospheric drag lowering and ultimately deorbiting the fragments. Simulation using the technique described here shows the phases A and B.

Hujsak [10], almost simultaneously with Jehn, derived a nonlinear dynamical model of the relative motion for the purpose of analysis of debris density from breakup through what Jehn would call phase B. By including select nonlinear terms and using a generalization of the Hill-Clohessy-Wiltshire equations for elliptical orbits, he was able to derive equations useful even far from the source orbit. He shows the inverse of the relative motion transformation is an approximate solution for the whole orbit Lambert (boundary value) problem. Including the geopotential J2 perturbation, he is able to analyze debris density differentially as a function of breakup velocity distribution. He found that particle density can vary by four orders of magnitude, and attributes increases in density away from the pinch point to perturbations. While we find the linearization of the relative motion approach too limiting and therefore take a different approach, we regard Hujsak’s study as a novel and pioneering idea for debris cloud evolution studies, and in a very significant sense inspired our present undertaking. His claim of four orders of magnitude variation is confirmed, but this variation is not due to perturbations exclusively.

In the following year, Housen [9] published an algorithm based on the Jacobian of the two-body initial value map and the Lambert solution, noted that it does not rely on linearization as had been done previously, and so produces an exact density based on the two-body force. The author, using initial velocity distributions that were described in the literature at the time, goes on to analyze the short-term effects of a point fragmentation by creating contour plots of the spatial density. We use the idea he described combining the Jacobian and Lambert solution, but we use different solution techniques for these, basing them on Battin [3].

Shortly after these publications, McInnes [19] introduced an analytical model using a hydrodynamic approximation that treats orbital debris in a continuum, fluid representation. From this, he is able to predict the effects of perturbations such as orbital drag, and the creation of new objects by launch and collisions to forecast the orbital debris environment many decades into the future. In the same year, Alan Jenkin presented at a conference and later published as [13] a paper in which he identified the Jacobian as the key part of the problem as a mapping from a debris cloud at a point to a spread in velocity at a later time of interest, and traces the history of the idea. His purpose was to characterize the dynamics of what he terms the non-stationary pinch point, and our own work confirms and extends his observations. A few years later, Frazzoli et al. [4] developed an analytic method based on integration over a fragmentation distribution of a single object density function. Though not noted in that paper, the transformation of variables technique that forms a central part of our approach is present in their equation (10). Their analysis produces general features of the cloud such as the pinch point and a line they call the “pinch line” opposite it, but in less detail than Jenkin.

More recently, Letizia [15] in her dissertation and papers coauthored with her advisers [16,17,18] continued and expanded McInnes’s approach using the continuity partial differential equation and other techniques from fluid mechanics. Like these investigators, we take a continuum approach and look at the density in space, but we are not exclusively considering an ensemble of orbital debris. The modeling here is for a much shorter time range and perturbations are not considered as in Letizia’s works; on the other hand, approximations are not made here.

Also in that year, Healy et al. [7] presented results based on the method described here that showed for the first time the exact spatial density for tens of orbits from a low-earth orbiting satellite over almost all of the source orbital plane. Additional features not previously understood or even reported, such as multiple bands, were observed and explained. The purpose of the present paper is to understand the patterns observed in that work that are due to orbital dynamics and independent of the initial velocity distribution. This was followed with Healy et al. [6], which investigated the persistence of patterns over time, and Healy et al. [8] which explored the effects of fragmentations in constellations and showed the three-dimensional pattern of a fragmentation.

Hansen [5] and Oltrogge and Vallado [20] presented independent work at the same conference addressing the evolution of debris shortly after a fragmentation. In the former, the author noted similar patterns to those we have observed, in which points in the interior of the initial velocity distribution become part of the exterior boundary of the spatial distribution at a later time. In the latter, the authors use a Lambert solver to simulate fragment evolution, and compare the results with the Iridium/Cosmos collision of 2009.

Often, fragmentation analyses proceed with explicit or implicit assumptions of the initial velocity distribution. Sometimes it is assumed, for example, to be isotropic in delta-v, or uniform in all directions, sometimes uniformly distributed in magnitude within some bound. The analysis may be based on some established fragmentation model, such as the NASA model Johnson et al. [14], but more often it is not. While the initial velocity distribution and orbital dynamics are both important in determining post-fragmentation density, there is much that can be concluded without a specific cloud creation model, and the goal of this paper is to explore observed structures that can be attributed to the dynamics and not the initial velocity distribution.

The purpose of this study is to introduce our approach to Housen’s general method and to use it to identify features of density dynamics from a point source that are independent of the initial velocity distribution. Our form of Housen’s method allows propagation over time intervals of any length. The features are studied by means of the dynamic admittance, which may be viewed as an extension of the analysis of Jenkin [13] who used a single Jacobian at a pinch point. In addition, the effects on the dynamic admittance for the presence of the earth (or any central body) is examined.

There are two component parts to the present form of Housen’s method: the transformation of variables method, and the solutions of the complete Lambert problem. The transformation of variables method tells us how to compute a density after mapping the domain space through a transformation. It needs two pieces of information: the inverse image points under the initial velocity to final position mapping, and the determinant of the Jacobian of the map at each of those points. In orbital terms, the initial velocity to final position mapping is simply the initial value problem for orbit dynamics. Here, unperturbed two-body (Kepler) motion is assumed, and the propagator uses the Lagrange coefficients f and g. The inverse map needed is obtained from the solution of the two-point boundary value problem, i.e., the Lambert problem. From every point r2 at any specified post-fragmentation time of interest and from the fixed initial location of the fragmentation r1, the task is to find all the initial velocities v1 that solve the two-body motion. Then, for each of those initial values, to find the determinant of the Jacobian of the orbit propagation problem.

Because our central purpose in this paper is to investigate the features of the spatial density flow due to the dynamics and independent of the source orbit and velocity distribution, we use the initial inertial velocity instead of the velocity relative to the source velocity. Our purpose here is specifically to step back from the application to look at the problem in a more straightfoward dynamical way; our publications referenced above provide some insight into the velocity-biased fragmentation behavior and so are only briefly repeated here. An added benefit of this analysis is, because the Lambert problem is axially symmetric, that only two dimensions need be considered, which makes the task easier than analyzing in three dimensions.

An Orbital Density Map

Although there are important differences, one view of the evolution of the orbital distribution is that of a fluid. In this analogy, the fragments are governed by the laws of orbital motion in the same way that the molecules of the fluid are governed by the forces of their environment. There are two fundamentally different ways to analyze fluid dynamics: the Lagrangian view and the Eulerian view [11]. In the former, each fluid element is described by giving its location as a function of time. In the latter, each location has a time-dependent function representing density. In orbital mechanics, the traditional, ephemeris, approach is the Lagrangian approach. We seek a shift from this approach to a new, Eulerian approach; at the same time, we make the shift to the continuum density from the discrete individual objects.

An orbital density map as we define it is the density of orbiting objects in Cartesian space around the earth. This density is a function of the position, and of time. A probability density or normalized number density function will integrate to one less the fraction lost to earth collision over the whole volume; a general density function will not have a specified integral, and might integrate to the total number of pieces, for example.

In order to find the density map, we will need to know the location of the source (e.g., fragmentation) r1 and the distribution of velocities. The distribution of locations r2(t) for an elapsed time t can be computed from orbit mechanics using the transformation of variables technique for density functions. The transformation (or change) of variables technique Sivia and Skilling [21] can be stated simply. Given a differentiable function F from a space X to a space Y, a density function GX in the space X can be transformed into one in the space Y by summing over all the points in the inverse image F− 1(y) (yY) divided by the absolute value of the determinant of the Jacobian matrix of F,

$$ {G}_{Y}(y) = \sum\limits_{x \in {F}^{-1}(y)} \frac{{G}_{X}(x)}{|\det J_{F}(x)|}, $$
(1)

where x and y can be vectors. The Jacobian matrix JF is the matrix of partial derivatives of the function F. See Fig. 1 for an illustration in one dimension: the points that transform into the y band come from three bands in x, and the amount contributed by each is inversely proportional to the absolute value of the derivative. When the inverse image is a continuum rather than a discrete set of points, the sum must be replaced with an integral.

Fig. 1
figure 1

Transformation of density in one dimension

In the present application, the space X is three dimensional, the velocity at the source point, and the space Y, also three dimensional, is the space in which we seek a density GY at a specific elapsed time later. We seek this spatial density function after an elapsed time given the velocity distribution GX at the time of the source event. For each value of the elapsed time after the source event, the function GY is computed for all points on a grid.

The spatial distribution of simulated fragments from a 900 km altitude circular orbit satellite using the NASA orbital debris model [14] is shown in Fig. 2. This plot, like all the density or dynamic admittance plots in this paper, is represented on a quasi-logarithmic scale, with black representing zero, and the darkest blue a small but positive density, ranging up to yellow and white for the highest density, progressing logarithmically. The well-known stationary “pinch point” where the fragmentation occurs is visible on the right, and high density bands are visible on the left, antipodal to the pinch point, i.e. on the opposite side of the earth from it. The point of fragmentation is in the middle vertically and to the right of the earth. This plot is shown in the source orbit plane; in three dimensions, the antipodal bands look something like bow ties [8].

Fig. 2
figure 2

Spatial density twenty-four hours after fragmentation

Dynamic Admittance

The orbital density map is an exact calculation of the resulting spatial density from the input velocity distribution at the single point source. It provides insight into the fragmentation dynamics, but it does not allow us to identify what in the observed pattern is due to the input velocity distribution, and what is due to the dynamics of orbit motion. To get an understanding of the latter, define the dynamic admittance

$$ D(y) = \sum\limits_{x \in {F}^{-1}(y)} \frac{1}{|\det J_{F}(x)|}. $$
(2)

While it is not possible to use this quantity directly to compute the resultant density, such a calculation will show in a general sense features that are exclusively due to the dynamics of the map F and not to the initial distribution. A higher dynamic admittance represents a place that is more accessible due to the dynamics alone than a lower one, and thus will be expected to have higher density, subject to the input velocity distribution. In the orbit application, the function F is the transformation from the source point velocity v1 to the final location r2, so the orbital dynamic admittance with a point source is

$$ D(\boldsymbol{v}_{1}) = \sum\limits_{\boldsymbol{v}_{1} \in {F}^{-1}(\boldsymbol{r}_{2})} \frac{1}{|\det J_{F}(\boldsymbol{v}_{1})|}, $$
(3)

with the Jacobian matrix being the 3 × 3 matrix of partial derivatives JF = r2/v1. Because of the cylindrical symmetry of the problem about the axis from the center of the earth to the source point, this need be computed only in a half plane from this axis.

As an aid to understanding the origin of the structure of the observed patterns in the dynamic admittance, it is helpful to compute an energy-limited dynamic admittance which is the same expression shown in Eq. 3, but with the summation limited to those inverse images whose initial velocity represent a total energy that does not exceed some specified bound ε,

$$ D_{\varepsilon}(\boldsymbol{v}_{1}) = \sum\limits_{\boldsymbol{v}_{1} \in {F}^{-1}(\boldsymbol{r}_{2})} \frac{\displaystyle \sigma \left( \varepsilon - \frac{{v_{1}^{2}}}{2} + \frac{\mu}{{r}_{1}}\right)}{|\det J_{F}(\boldsymbol{v}_{1})|}, $$
(4)

where v1 = |v1|, r1 = |r1| and σ is a step function with σ(x) = 1 if x ≥ 0, and zero otherwise.

In the remainder of the paper, there are no units or scale given for any form of dynamic admittance because our only concern are the relative values.

There are clearly two pieces to this computation: finding the inverse images F− 1(r2) at each point of interest, and computing the forward Jacobian determinant. In the current context, the inverse image calculation amounts to the two-point boundary value, or Lambert, problem. The forward function F is the initial value problem for orbit mechanics; given an initial position r1 and velocity v1, find the position v2 and the velocity r2 after time t. We call this approach Housen’s method, after the description in [9]. Our computational technique, differing from Housen’s, is described in the next two sections. By obtaining all physical Lambert solutions, the method works for all time intervals.

Computation of Dynamic Admittance

The Complete Lambert Solution

The inverse image needed for the transformation of variables (3) is a matter of finding the correct initial velocities v1 such that together with the known source location r1 and the elapsed time t, the object reaches the final location r2 after the elapsed time. That is, t, r1 and r2 are given, and v1 must be found. This is the classic Lambert problem, or two-point boundary value problem for orbit mechanics. Over the centuries since it was first posed, numerous methods have been developed for solving this problem. Some methods do better or are only applicable for certain kinds of orbits, for example, short arcs, bound orbits, or those that go less than a whole orbit. For the present application, we cannot restrict the kinds of orbits, and so need a method that is good for as many kinds of orbits as possible.

Between any two position vectors, an orbit can go in one of two directions. If the geometric transfer angle is less than 180 degrees, the arc is considered short way, if between 180 and 360 degrees, it is considered long way. Note the terms refer to the angle; the amount of time is the same in both cases, as assumed in the problem. If the earth center and the two points are co-linear, there is no unique solution because the plane is undetermined, though some orbital elements can be determined. In this case, the inverse image will have a continuum solution.

If the trajectory passes the initial point r1 at least once before t has elapsed, the trajectory is a whole orbit trajectory (customarily called “multi-revolution”); it will pass r1N ≥ 1 times before the elapsed time is up. If it does not pass r1 before the elapsed time is up, it is a zero orbit case, or N = 0. All possible routes for non-negative integer N need to be examined in the Lambert algorithm.

A zero orbit trajectory can be of any conic section, but a whole orbit must be a bound orbit, i.e., circular or elliptical. While there is a mathematical solution for a zero orbit trajectory for any pair of points and time, assuming unlimited speed is possible, the whole orbit case has a minimum time. As N increases, this minimum time increases. Therefore, the solution search is bounded; if N is high enough that there are no solutions because the elapsed time is less than the minimum time, then no higher value of N need be searched.

Our application is unique among typical uses of Lambert solvers: we must obtain all solutions for a given set of boundary conditions and time; that is, all inverse images are necessary to form the sum in Eq. 1. This includes all conic sections circular, elliptical, parabolic, and hyperbolic, both directions, and all possible revolution counts N. Hence we call it the complete Lambert solution.

The preceding exposition refers to the mathematical solution of the Lambert problem. Not all mathematical solutions are physical solutions; if the geocentric distance becomes less than the radius of the earth at any point on the trajectory, that route is not physically possible. Therefore, after determining a solution mathematically, it must be checked for minimum geocentric distance. If it is a whole orbit trajectory, that is the perigee radius. If it is zero orbit, and the trajectory goes through perigee, then again the perigee radius is the minimum geocentric distance. If the trajectory doesn’t go through perigee, then the minimum geocentric distance is given by the smaller of the distances at the endpoints r1 or r2.

Regardless of the Lambert algorithm used, an iterative solver is always necessary, because all algorithms involve some sort of root solving of a function that is not invertible analytically. The solver is typically a standard algorithm such as a Newton method (if a derivative is available), or a bisection solver.

Figure 3 shows all possible routes between the two points indicated over the specified time; there are fourteen possible mathematically; three are not possible physically because the trajectories intersect the earth (indicated by a dashed line). The notation in the boxes attached to each route is the value of N, followed by “S” or “L” for short or long way, respectively. The designation “<” and “>” where N > 0 is explained below. The relation signs that follow if N > 0 are explained in the next section. Note that the velocities v2 (tangents to the orbits at the endpoint) are all very different; this is one way in which the “orbital fluid” is not like a normal fluid.

Fig. 3
figure 3

The complete Lambert solution for r1 = (7.278,0)Mm, r2 = (− 28,8.82)Mm, and an elapsed time of twenty-four hours

Lambert Algorithm

Battin’s hypergeometric method Battin (3, sec. 7.1) for solving the Lambert problem gives consistently complete sets of solutions for all types of orbits. The core Lambert algorithm computes a quantity x > − 1; the Lagrange coefficients are then used to solve for the initial velocity. Assuming the value of x is known, the Battin method computes a dimensionless elapsed time (Battin 3, (7.32))

$$ \sqrt{\frac{\mu}{a_{m}^{3}}} t = \frac{2 \pi N}{\left( 1-x^{2}\right)^{3/2}} + \frac{4}{3} {\eta}^{3}_{2} {F}_{1}\left( 3,1;\frac{5}{2};{S}_{1}\right) + 4 \lambda \eta, $$
(5)

where (see Battin equations (6.9), (6.10), (7.6) and (7.30)),

$$\begin{array}{@{}rcl@{}} s &=& \frac{{r}_{1}+r_{2}+c}{2} \end{array} $$
(6a)
$$\begin{array}{@{}rcl@{}} a_{\mathrm{m}} &=& \frac{s}{2} \end{array} $$
(6b)
$$\begin{array}{@{}rcl@{}} \lambda &=& \pm \sqrt{1 - \frac{c}{s}} \end{array} $$
(6c)
$$\begin{array}{@{}rcl@{}} y &=& \sqrt{1 - \lambda^{2}\left( 1-x^{2}\right)} \end{array} $$
(6d)
$$\begin{array}{@{}rcl@{}} \eta &=& y - \lambda x \end{array} $$
(6e)
$$\begin{array}{@{}rcl@{}} S_{1} &=& \frac{1}{2} (1 - \lambda - x \eta); \end{array} $$
(6f)

am is the semimajor axis of the minimum energy orbit and c = |rv2rv1| is the length of the chord connecting the endpoints. The sign for λ is chosen by which way the orbit went, short (+ ) or long (−) way. The variable x > − 1 is the unknown that must be solved; since Eqs. 6a to 6c are independent of x, they may be pre-computed before starting the iteration. The function 2F1 is a hypergeometric function. The gravitational constant is given by μ; the factor multiplying the elapsed time t is the mean motion of the minimum energy orbit.

Many root-solving techniques (e.g. Newton’s method) can make use of the derivative. The derivative of the hypergeometric function is another hypergeometric function [1], so the derivative of the right hand side of Eq. 5 is

$$ \frac{6 \pi N x}{\left( 1-{x}^{2}\right)^{5/2}} -\frac{4}{5}(px+\eta){\eta}^{3}_{2} {F}_{1}\left( 4,2;\frac{7}{2};{S}_{1}\right) + 4 {\lambda} p + 4 p {\eta}^{2}_{2} {F}_{1}\left( 3,1;\frac{5}{2};{S}_{1}\right) $$
(7)

with

$$ p=\frac{x\lambda^{2}}{y}-\lambda. $$
(8)

For the zero orbit case, presuming there is enough velocity and no obstruction from the earth, there is no minimum time to get to another point (except the time it takes light to travel between the points), so a universal method, one that is capable of finding solutions for any conic section, is desirable. The Battin method is one such universal method. On the other hand, for the whole orbit case N ≥ 1, there is a minimum time tmin(N) such that if t < tmin(N), no solution exists. For each additional whole orbit, this minimum time increases tmin(N + 1) > tmin(N), so if it is found for a particular N that the minimum time is too high, the search with increasing N may be terminated for the input conditions. This means that there is a maximum orbit count, a function of the desired elapsed time, Nmax(t), beyond which no solutions are to be found. To be precise, there are two separate maximum counts, one for the short way and one for the long way. They usually have the same value but sometimes differ by one.

Figure 4 shows a curve of dimensionless time versus Battin’s Lambert quantity x by whole orbit N for r1 = (7.278,0)Mm, r2 = (− 10,3.75)Mm. For this case, each dimensionless time unit represents 1332.05s, so twenty-four hours is equivalent to 64.862 in dimensionless time, represented by the dashed line in the figure. The plot shows only short way curves; it is clear that N = 0 has one such solution and each of N = 1,…,9 has two solutions for a specified time of flight. The minimum time for N ≥ 10 is too large for there to be any solutions for twenty-four hours. The plot of long way curves looks similar. If x ≥ 1, the orbit is unbound, which is only possible if N = 0.

Fig. 4
figure 4

Typical plot of dimensionless time versus x by whole orbit N; short way only shown

The transformation of variables technique for density calculation requires that every Lambert solution be found. This means each route possible in the elapsed time; the whole orbit count N, the two directions (short and long way), and for N ≥ 1, x either less than or greater than the value of x when Eq. 5 is at its minimum (indicated by < and > in Fig. 3). For every set of boundary conditions r1, r2 and elapsed time t, there is one route per direction each for N = 0 and two of each for N ≥ 1.

The sequence of computation chosen is driven by efficiency of computation. In order to find the minimum time quickly and avoid its unnecessary recomputation, and to accelerate the determination of solutions at one time step by using the solutions obtained from a previous time step, a single processor is assigned to one r2 point, and the entire set of elapsed times for which a solution is desired is computed sequentially, starting with the shortest time. The N = 0 solutions are first found, then the whole orbit minimum times. The maximum N considered is Nmax(t) at any elapsed time. Then the solution for each of the possible routes is found. These constitute the complete Lambert solution for that time step.

The dimensionless time (5) is computed as a function of x, as is its derivative Eq. 7. Because the function is not defined for x ≤− 1, a Newton solver, which is an unconstrained optimization, is not ideal. It is easy to find cases where a guess of a valid x will produce a derivative such that the next iteration of the Newton method gives an invalid value x ≤− 1. Therefore, we use a hybrid bracketing and Newton solver method. In this approach for the first time step at which there is a solution, the solver moves toward the extreme values of x (depending on the derivative of the function at the current evaluation point), halving the distance to the boundaries x = ± 1 for N ≥ 1 until it has bracketed the desired value. The case of N = 0 uses linear extrapolation to x > 1 to find the lower bracketing value if x = 1 is not below. Then, a Newton solver is used to find the exact (within specified tolerance) value of x; these steps ensure it will converge. For subsequent time steps for that location r2, the Newton solver is applied directly to the solution obtained on the previous time step for that route (N, short/long way, and smaller/larger value of x). This requires the time steps be small enough that this always converges.

Solving for Orbits

Once the solver has found a value of x that solves (5), it is necessary to find the initial velocity v1 and Jacobian determinant J for that solution. This is done with a universal variable formulation of the Lagrange coefficients, ensuring that the solution is correct for both bound (circular, elliptical) and unbound (parabolic and hyperbolic) orbits. The universal variable α is the inverse of the semimajor axis and is thus always defined: negative for hyperbolic orbits, zero for parabolic, and positive for bound orbits. From Battin (3, (7.20)),

$$ \alpha = \frac{1-x^{2}}{a_{\mathrm{m}}}. $$
(9)

The other universal variable needed is the generalized anomaly χ. In order to define this, define σ as

$$ \sigma = \frac{\boldsymbol{r} \cdot \boldsymbol{v}}{\sqrt{\mu}}; $$
(10)

see Battin problem 3–23, p. 128. At the terminal points, the values of σ are expressed in terms of x, λ (6c), η (6e), and am (6b) as

$$\begin{array}{@{}rcl@{}} \sigma_{1} &=& \frac{2\lambda a_{\text{m}} - r_{1}(\lambda + x\eta)}{\eta\sqrt{a_{\text{m}}}}, \end{array} $$
(11a)
$$\begin{array}{@{}rcl@{}} \sigma_{2} &=& -\frac{2\lambda a_{\text{m}} - r_{2}(\lambda + x\eta)}{\eta\sqrt{a_{\text{m}}}}; \end{array} $$
(11b)

see Battin eq. (7.35). The generalized anomaly χ is now computed from the known variables

$$ \chi = \alpha\sqrt{\mu} t + \sigma_{2} - \sigma_{1}; $$
(12)

see Battin eq. (4.86).

Lagrange (f and g) coefficients are used for propagation. Because an orbit lies in a plane, any point on the orbit and any velocity vector can be described by two vectors that span the plane. Because for realistic orbits they cannot be co-linear, the initial position r1 and velocity v1 may serve as those basis vectors, with suitable coefficients,

$$\begin{array}{@{}rcl@{}} \boldsymbol{r}_{2} &=& f(t)\boldsymbol{r}_{1} + g(t)\boldsymbol{v}_{1}, \end{array} $$
(14a)
$$\begin{array}{@{}rcl@{}} \boldsymbol{v}_{2} &=& \dot{f}(t)\boldsymbol{r}_{1} + \dot{g}(t)\boldsymbol{v}_{1}. \end{array} $$
(14b)

The initial conditions are must match these values when t = 0,

$$\begin{array}{@{}rcl@{}} f(0) &=& 1, \quad{} \dot{f}(0) = 0 \end{array} $$
(15a)
$$\begin{array}{@{}rcl@{}} g(0) &=& 0, \quad{} \dot{g}(0) = 1. \end{array} $$
(15b)

These form a set of vector differential equations which can be solved for the coefficients f and g. Although they are time dependent, we will drop the function arguments henceforth. However, it is important to remember that these “coefficients” are actually functions of r1 and v1, because we will need to take the partial derivatives for the Jacobian.

The Lagrange coefficients for two-body propagation are computed from the functions U1(χ;α) and U2(χ;α),

$$\begin{array}{@{}rcl@{}} f &=& 1 - \frac{U_{2}(\chi;\alpha)}{{r}_{1}}, \quad{} g = \frac{{r}_{1}U_{1}(\chi;\alpha) + \sigma_{1}U_{2}(\chi;\alpha)}{\sqrt{\mu}}, \end{array} $$
(16a)
$$\begin{array}{@{}rcl@{}} \dot{f} &=& -\frac{\sqrt{\mu}U_{1}(\chi;\alpha)}{{r}_{1}r_{2}} , \quad{} \dot{g} = 1 - \frac{U_{2}(\chi;\alpha)}{{r}_{2}}; \end{array} $$
(16b)

see Battin eq. (4.84). The U functions are defined alternately in terms of sines and cosines, either circular or hyperbolic (see Battin problem 4–21), or in terms of hypergeometric functions and continued fractions (Battin sec. 4.7). The former may have numerical computation problems near a parabolic (α = 0 or x = 1) orbit; however, for the present application, performance was found to be satisfactory.

The Lagrange coefficients (1) may be rearranged to solve for v1,

$$ \boldsymbol{v}_{1} = \frac{\boldsymbol{r}_{2} - f\boldsymbol{r}_{1}}{g}. $$
(17)

When computing a spatial density (1), this velocity v1 is needed for the look-up in the velocity distribution function GX at the source for the density transformation. For the dynamic admittance calculation, it is needed if an energy limitation Eq. 4 is applied. In either case, it is needed to compute the Jacobian determinant of the forward map (“Jacobian Determinant”).

This procedure determines the mathematical solutions to the Lambert problem. Not all of these are physical solutions, however; any trajectory that intersects the earth will not contribute to the spatial density or dynamic admittance after the collision. The three routes shown with dashed lines in Fig. 3 are mathematical but not physical solutions.

Collision determination depends on the whole orbit count N. If it is positive, the minimum radial distance over the orbit segment is the radius of perigee; if this is less than the earth radius, then the trajectory is not physical. If it is zero, the minimum distance over the orbit segment is computed from knowledge of whether the satellite is ascending or descending at the terminal points. If both points are on the ascending segment r1v1 > 0 and r2v2 > 0 (past perigee before apogee), and point 1 is higher than point 2 (r1 > r2), or if both points are on the descending segment r1v1 < 0 and r2v2 < 0 (past apogee before perigee), and point 1 is lower than point 2 (r1 < r2), then the orbit segment passes through perigee, so the minimum distance is the perigee radius. If none of these cases hold, then the trajectory does not pass through perigee and the minimum geocentric distance is the smaller of r1 and r2.

For the boundary conditions and time whose short way curves are plotted in Fig. 4, there are in total thirty-eight mathematical solutions (both short and long ways) which includes two physical solutions (one each short and long way) for each of N = 0, 7, 8, and 9, and no more, for a total of eight. The corresponding curves are marked in red on the plot; only one each of the whole orbit solutions is a physical solution, so there are four total short way solutions. There are also four long way solutions, ungraphed.

Jacobian Determinant

The Jacobian matrix needed in Eq. 3 is the 3 × 3 matrix

$$ \renewcommand{\arraystretch}{2} J = \frac{\partial \boldsymbol{r}_{2}}{\partial \boldsymbol{v}_{1}} =\left[ \begin{array}{r r r} \displaystyle \frac{\partial r_{2I}}{\partial {v}_{1I}} & \displaystyle \frac{\partial r_{2I}}{\partial {v}_{1J}} & \displaystyle \frac{\partial r_{2I}}{\partial {v}_{1K}} \\ \displaystyle \frac{\partial r_{2J}}{\partial {v}_{1I}} & \displaystyle \frac{\partial r_{2J}}{\partial {v}_{1J}} & \displaystyle \frac{\partial r_{2J}}{\partial {v}_{1K}} \\ \displaystyle \frac{\partial r_{2K}}{\partial {v}_{1I}} & \displaystyle \frac{\partial r_{2K}}{\partial {v}_{1J}} & \displaystyle \frac{\partial r_{2K}}{\partial {v}_{1K}} \end{array} \right], $$
(18)

where derivatives of each of the three Cartesian inertial components I, J, K are shown. This matrix is one 3 × 3 block of the 6 × 6 state transition matrix for orbit dynamics. It was computed by Battin (3, eq. (9.84)), and Battin’s presentation was consolidated by Arora et al. [2]. We compute starting with the known U1 and U2, χ and Lagrange coefficients f and g,

$$\begin{array}{@{}rcl@{}} U_{3} &=& \sqrt{\mu}(t-g) \end{array} $$
(19a)
$$\begin{array}{@{}rcl@{}} U_{4} &=& U_{1} U_{3} - \frac{1}{2} ({U_{2}^{2}} - \alpha {U_{3}^{2}}) \end{array} $$
(19b)
$$\begin{array}{@{}rcl@{}} U_{5} &=& \frac{1}{\alpha} \left( \frac{1}{6}\chi^{3}-U_{3} \right) \end{array} $$
(19c)
$$\begin{array}{@{}rcl@{}} \bar{C} &=& \frac{1}{\sqrt{\mu}} \left( 3U_{5} - \chi U_{4} -\sqrt{\mu}U_{2} t\right) \end{array} $$
(19d)
$$\begin{array}{@{}rcl@{}} \frac{\partial \boldsymbol{r}_{2}}{\partial \boldsymbol{v}_{1}} = \frac{{r}_{1}}{\mu}(1&-&f) \left[ (\boldsymbol{r}_{2}-\boldsymbol{r}_{1})\otimes \boldsymbol{v}_{1} - (\boldsymbol{v}_{2}-\boldsymbol{v}_{1})\otimes \boldsymbol{r}_{1} \right]\\ &&+ \frac{\bar{C}}{\mu} \boldsymbol{v}_{2} \otimes \boldsymbol{v}_{1} + gI, \end{array} $$
(19e)

where I is the 3 × 3 identity matrix and “⊗” represents the outer product, forming a 3 × 3 matrix from two 3-vectors. Note that U5 and thus \(\bar {C}\) and the Jacobian become unbounded as the orbit approaches parabolic (α = 0) from either the elliptical (α > 0) or hyperbolic α < 0 side. Arora et al. address its evaluation for large high elliptical eccentricity, but for truly parabolic orbits, the derivative is not defined. The transformation of variables computation uses the inverse absolute value of the determinant, which therefore goes to zero as an orbit approaches parabolic from either side. In other words, parabolic orbits do not contribute anything to orbital density. Hyperbolic orbits, being unbounded, only contribute transients; given a large enough elapsed time, the effects of hyperbolic orbits are seen only very far from the earth.

Simulation

Dynamic Admittance Structure from a Point Source

The dynamic admittance at a sequence of elapsed times from a source r1 = 7278 km from the center of the earth is shown in Figs. 567891011 and 12. In these plots, the source is on the right side of the earth. Video 1, in the electronic supplementary file accompanying this paper,Footnote 1 shows the evolution through eight hours.

Fig. 5
figure 5

Dynamic admittance at 20 minutes

Fig. 6
figure 6

Dynamic admittance at 40 minutes

Fig. 7
figure 7

Dynamic admittance at 55 minutes

Fig. 8
figure 8

Dynamic admittance at 1 hour 30 minutes

Fig. 9
figure 9

Dynamic admittance at 2 hours 20 minutes

Fig. 10
figure 10

Dynamic admittance at 3 hours 55 minutes

Fig. 11
figure 11

Dynamic admittance at 5 hours 30 minutes

Fig. 12
figure 12

Dynamic admittance at 24 hours

At twenty minutes elapsed time (Fig. 5), there is a crescent-shaped region of density with bound orbits which mates in the scalloped center of a V-shaped region that has unbound orbits. At forty minutes (Fig. 6), the bound orbit density region is almost closed on the antipodal side, and the unbound orbit region is almost completely covered in the field of view by the bound orbits. At fifty-five minutes (Fig. 7), there is just visible a high-density band at the inside edge of the unbound orbits density region (visible in the lower and upper right corners). At one hour thirty minutes, or approximately one orbital period of a circular orbit at the source radius, the well known stationary pinch point starts to form (Fig. 8). At about two hours and twenty minutes, the C-shaped region has closed on the antipodal side to form a band (Fig. 9). This is the first band of what will be several. At three hours fifty-five minutes (Fig. 10), the second band has formed, and at five and a half hours the third has formed (Fig. 11). At twenty-four hours, several more bands are visible (Fig. 12).

In the first hour and a half, the unbound orbits expand away from the earth on the source side; the spherical scallop gives way to a bubble which has a high-density ridge further away from the earth. See video 2.

Looking just at the twenty-four hour dynamic admittance, as the maximum total energy is lowered to − 0.25 circular orbit energy units (\(\sqrt {\mu /2r_{1}}\), or COEU), some antipodal infill between the bands disappear; see Fig. 13. As the maximum total energy is lowered further, Figs. 1415 and 16, more of the outer bands disappear. Video 3 shows the dynamic admittance as the maximum total energy increases from − 1 to 0.

Fig. 13
figure 13

Energy-limited dynamic admittance, limited to ε = − 0.25COEU

Fig. 14
figure 14

Energy-limited dynamic admittance, limited to ε = − 0.50COEU

Fig. 15
figure 15

Energy-limited dynamic admittance, limited to ε = − 0.75COEU

Fig. 16
figure 16

Energy-limited dynamic admittance, limited to ε = − 1.0COEU

As one might expect, the regions more distant from the earth are not accessible if the energy is not high enough. The observed sharp drop in dynamic admittance at the outer edges of each band coincides with a constant total energy. Figure 17 shows the dynamic admittance by location on the source axis; this plot is a cross-section of Fig. 12 from left to right, and in the middle vertically, except that the dynamic admittance (vertical axis) is plotted on a linear, rather than logarithmic scale. This axis has been truncated. The color coding represents the source energy of the dynamic admittance, the darker colors (purple) are lower energy and the lighter colors (green/yellow) higher energy. The band structure is clearly seen as spikes at the outer edge. The antipodal infill is observable as small light green (high energy) density.

Fig. 17
figure 17

Dynamic admittance along the source axis by energy (color) at twenty-four hours

Table 1 shows the location of each band’s outer edge. The second column gives the total energy (in circular orbit energy units) at which the outer band edge appears in its final location, the third column is the computed distance of the bend edge from the center of the earth along the antipodal line as will be explained in “The Band Structure”, and the last column is the distance estimated from the dynamic admittance computation, as given in Fig. 17, and therefore only available to the nearest 40 km, the spacing between pixels.

Table 1 Total energy and geocentric distance at each of the outer band edges on the antipodal line

The dynamic admittance is significantly affected by the presence of the earth; approximately half the routes at twenty-four hours for the simulation shown above are lost to earth collision. In effect, the earth casts a dynamic “shadow” on the rest of space; this shadow does not follow the rules of optics, however. Some of the features in the dynamic admittance plots are due to this dynamic occlusion, and one way to understand their origin and nature is to change the earth radius in the earth intersection computation (“Solving for Orbits”) from zero to the radial distance of the point source, r1, for a fixed elapsed time.

The dynamic admittance effect on a source 7278 km at an elapsed time of twenty-four hours from a planet with earth’s gravity, but no radius, is shown in Fig. 18. The gravitation center is at the right end of the high-density line (lighter color) in the middle. Figures 1920 and 21 show the effect on dynamic admittance as the planet radius is increased to 0.25, 0.50 and 0.75 of the source distance. In each figure, the upper plot shows the dynamic admittance on the same quasi-logarithmic scale as Fig. 18, and the lower plot show the ratio of the dynamic admittance to the unoccluded dynamic admittance in that figure. The earth is depicted as a green disk in these figures. The gray-scale ranges from black (ratio 0.0, or dynamic umbra), through gray values representing dynamic penumbra, to white (ratio 1.0) with no dynamic shadow. Finally, Fig. 22 shows the occlusion effect if the source is at the surface of the planet; note that there is not complete absorption by the planet, but some density at the source point and on the antipodal line. Videos 4 and 5 show, respectively, the dynamic admittance and ratio to the unoccluded dynamic admittance as the occlusion radius grows from 0 to r1.

Fig. 18
figure 18

Dynamic admittance from source 7278 km from gravitation center with earth μ but no radius

Fig. 19
figure 19

Dynamic admittance and ratio to unoccluded admittance for planet radius r1/4

Fig. 20
figure 20

Dynamic admittance and ratio to unoccluded admittance for planet radius r1/2

Fig. 21
figure 21

Dynamic admittance and ratio to unoccluded admittance for planet radius 3r1/4

Fig. 22
figure 22

Dynamic admittance from source 7278 km from gravitation center with earth μ and radius equal to the source

The Band Structure

To understand why the bands exist, consider the antipodal line only. To simplify the problem initially, assume that the velocity is perpendicular to the radial vector at the source. Therefore, each trajectory starts at perigee or apogee and arrives at the antipodal line, 180 degrees from the source point, at apogee or perigee respectively. The time of flight is thus an odd multiple of one half the orbital period of this orbit. In other words, the trajectory makes \(N+\frac {1}{2}\) orbits for integer N ≥ 0 between the source and the arrival at the antipodal point.

This means that when looking at the antipodal line at an elapsed time t, the semimajor axis has a discrete set of possible values,

$$ a = \sqrt[3]{\mu \left[ \frac{t}{(2N + 1)\pi} \right]^{2}}. $$
(20)

As a consequence, the antipodal geocentric distance is

$$ r_{2} = 2a-r_{1} = 2\sqrt[3]{\mu \left[ \frac{t}{(2N + 1)\pi} \right]^{2}} -r_{1} $$
(21)

because the apse (major axis) lies on the r1r2 line, i.e., there is no radial component of velocity at those points. It is easy to see the bands in Eq. 21; as the sequence N = 0,1,… produces successively smaller terms in the cube root, which means r2 is smaller. These values are computed for r1 = 7278 km and twenty-four hours in the third column of Table 1. So the outermost band is the half-orbit case, the next one closer to the earth is the one-and-a-half orbit case, and so on. There is a maximum N for a given time t because otherwise the semimajor axis becomes low enough that perigee will be inside the earth (radius r), so

$$ N < \frac{1}{2} \left[ \frac{t}{\pi} \sqrt{\mu \left( \frac{2}{{r}_{\boldsymbol{\oplus}}+r_{1}} \right)^{3}} -1 \right]. $$
(22)

This limit will increase as the elapsed time increases. After sufficient time t, low values of N are limited by the maximum energy available. To travel only half an orbit over a long elapsed time t requires that semimajor axis, and thus apogee radius, be very high, and if the energy is limited, that will not be possible.

When there is a radial component to the velocity vr≠ 0, finding the location of the bands in geocentric distance on the antipodal axis is more complicated, but reduces to a simple formula. The presence of a radial component means that the true anomaly of the source point (and antipodal point) on the new orbit can be any value (and not only 0 or π), so the antipodal distance must be computed generally.

The velocity components at the source point may be computed from the velocity in the perifocal frame by rotation by the source point true anomaly 𝜃,

$$ \left[ \begin{array}{r} v_{\mathrm{r}} \\ {v}_{\perp} \end{array} \right] = \sqrt{\frac{\mu}{p}} \left[ \begin{array}{rr} \cos \theta & \sin \theta \\ -\sin \theta &\cos \theta \end{array} \right] \left[ \begin{array}{r} -\sin \theta \\ e + \cos \theta \end{array} \right] = \sqrt{\frac{\mu}{p}} \left[ \begin{array}{r} e \sin\theta \\ 1+e\cos\theta \end{array} \right] $$
(23)

with eccentricity e, semilatus rectum p, and gravitational constant μ. The antipodal distance (at true anomaly 𝜃 + π) is then seen to be a function of the tangential velocity only

$$ r_{2} = \frac{p}{1-e\cos\theta} = \frac{\mu p}{2\mu-hv_{\perp}} = \frac{{{r}_{1}^{2}}{v}_{\perp}^{2}}{2\mu-r_{1}{v}_{\perp}^{2}}, $$
(24)

using angular momentum \(h = \sqrt {\mu p}\) and rv = h.

The antipodal distance can be expressed in terms of the semimajor axis a, radial component of velocity vr, and source distance r1. Dividing numerator and denominator of the right-hand side of Eq. 24 by \(r_{1}{v}_{\perp }^{2}\) gives

$$ r_{2} = r_{1} \left( \frac{2\mu}{{r}_{1}{v}_{\perp}^{2}} -1 \right)^{-1}. $$
(25)

The transverse component of velocity is dependent on the semimajor axis and radial component of velocity,

$$ {v}_{\perp}^{2} = \mu\left( \frac{2}{{r}_{1}} - \frac{1}{a} \right) - {v}_{\mathrm{r}}^{2}. $$
(26)

Substituting Eq. 26 into Eq. 25, and rearranging gives the solution of r2 for any value of vr,

$$ r_{2} = 2\left( \frac{1}{a}+\frac{v_{\mathrm{r}}^{2}}{\mu}\right)^{-1}-r_{1}, $$
(27)

which clearly reduces to r2 = 2ar1 when vr = 0.

The semimajor axis is related to the elapsed time by generalization of Eq. 20 for any radial velocity, using the mean anomaly M1 at r1, and the antipodal mean anomaly M2,

$$ a = \sqrt[3]{\mu \left[ \frac{t}{2\pi N + M_{2} - M_{1}} \right]^{2}}. $$
(28)

While M1 and M2 are dependent on the radial velocity and their computation involves Kepler’s equation, the difference is always at least half a revolution M2M1π, with equality holding only when vr = 0. By substituting Eq. 28 into Eq. 27, it is thus seen that r2 is at a maximum when vr = 0, holding everything else equal. Therefore, the sharp outer edge of each band is associated with zero radial velocity, and width of the band toward the earth is due to positive radial velocity.

Validation

The calculations described in “An Orbital Density Map” and “Dynamic Admittance” may be validated in two different ways, by propagation to check the results of the Lambert solution, and sampling of initial velocities to check the density calculation. The purpose of the first check is to verify that the algorithms have been implemented correctly and to adjust the tolerance on the solver determining x so that the correct r2 point is hit at the elapsed time of interest. With an absolute tolerance of 1.0 × 10− 6, all points after an elapsed time of twenty-four hours are within 2 km, and about 2.4 × 10− 5 points, or twenty-four out of a million, are between 1 km and 2 km off. This was judged acceptable for the purposes of this study.

The results of the density map calculation can be validated by sampling the orbits from the initial velocity distribution, propagating to the desired time, counting numbers at each location and dividing by the total number. Here, the sample is selected from a low discrepancy sequence.

Since typical densities are found to be in the range 10− 13 to 10− 9km− 3, to obtain the expectation of one particle in a cube 24 km on a side (13824 km3), a minimum of 107 particles are needed in the simulation, with more providing an enhanced ability to check fidelity. We have chosen instead about 108 particles, to give more than the minimum and get better resolution, so this should give at least one significant figure of resolution.

The particle initial velocities are determined from a Sobol sequence in three dimensions, which is a low-discrepancy sequence designed to fill the space with minimum variation. The initial velocity change distribution is uniform and isotropic up to a maximum of 2 kms− 1difference from a circular orbit velocity at r1. The points that fall outside the velocity sphere are not propagated; the remaining particles are propagated for the desired elapsed time using the method of Lagrange coefficients. The final location of each particle is then determined, and from that the nearest node, which determines its bin. A bin is determined from the grid by finding the nearest node. For the third dimension, which is computed by rotating the Lambert solution cylindrically along the source axis, the nearest node is computed on the approximated straight line, not an arc. The length of this box depends on the cylindrical radius and the rotation angle; for small angles, and locations close to the axis, the fidelity of this approximation is high.

Figure 23 shows a comparison of computed and sampled density in the source orbital plane using this technique. It can be seen that the agreement is quite good.

Fig. 23
figure 23

Density approximation from sampling and exact density from Housen’s method

Conclusions

The spatial density resulting from Keplerian orbital motion from a point source (for example, a fragmentation) may be determined exactly using Housen’s method. Dynamic admittance, defined as the sum over all inverse images of the absolute Jacobian determinant, is independent of the initial velocity distribution, and thus shows the characteristics of the density arising from the orbital dynamics. The inverse images are found by solving the Lambert problem using Battin’s hypergeometric method, and for each the Jacobian determinant is found from the partial derivatives r2/v1 using universal variable Lagrange coefficients.

The structure observed in the dynamic admittance is a series of bands on the opposite side of the earth from the source point, the spacing between the bands increasing as distance from the earth increases. Near the line opposite the source point, there is high density relative to other areas. Within each band, moving radial outward, the density increases gradually. At the outer edge of the band, it drops suddenly to the next band. These sudden drops define the division between the band edges. The bands blend into each other moving away from antipodal line. This continuum high density region narrows to a pinch point of very high density at the source point.

Each of the band edges moving outward is found to originate from increasing total energy, with the separation in total energy about 0.05 circular orbit energy units \(\left (\sqrt {\mu /2r_{1}}\right )\) for a point r1 = 7278 km from the earth center (corresponding to 900 km altitude) for the inner rings, rising to about 0.10 circular orbit energy units at 66760 km. There is a small gap of zero density between the earth surface and the first band, the first and second band, and the second and third band. These gaps partially fill in near the antipodal line at an energy of − 0.171 circular orbit energy units. The bands are explained as arising from odd half integer orbits, and the predicted locations agree closely with those determined by examination of the dynamic admittance.

The number of Lambert solutions and thus the dynamic admittance is reduced by the presence of the earth, through which many trajectories would pass were it not present. Simulating all the trajectories with a zero-radius planet, the densities are found to be much higher, particularly on the antipodal line, and around the pinch point, where it is a ball of high density, rather than a point. The planet creates an teardrop-shaped zero-density region. As the planet size is increased, the pinch point forms. An examination of the ratio of the occluded to unoccluded densities, shows that, in contrast to light optics, the dynamic “shadow” occurs on the same side of the occlusion as the source, rather than the opposite side. In the dynamic case, a conical region on the opposite side is completely free of shadow. On the same side as the source, there are penumbral bands.