Realistic gravitational focusing of meteoroid streams

The number density and flux of a meteoroid stream is enhanced near a massive body due to the phenomenon known as gravitational focusing. The greatest enhancement occurs directly opposite the massive body from the stream radiant: as an observer approaches this location, the degree of focusing is unbound for a perfectly collimated stream. However, real meteoroid streams exhibit some dispersion in radiant and speed that will act to eliminate this singularity. In this paper, we derive an analytic approximation for this smoothing that can be used in meteoroid environment models and is based on real measurements of meteor shower radiant dispersion.


INTRODUCTION
As meteoroids approach a massive body such as the Earth, they are accelerated by the body's gravitational field. The initial trajectory is bent and the meteoroid accelerates. The number density near the massive body increases due to the inward deflection of meteoroids, and the flux increases further still due to the increase in speed. This phenomenon is known as gravitational focusing, and when trajectories are blocked by intersection with the massive body, it is called planetary shielding. Both effects must be taken into account in order to model the meteoroid environment near a massive or sizeable body.
A global flux enhancement factor can be derived by invoking conservation of energy and angular momentum (Öpik 1951;Kessler 1972;Jones & Poole 2007). However, gravitational focusing is highly anisotropic (Divine 1992;Staubach et al. 1997;Matney 2002;Jones & Poole 2007) and the number density is at its highest directly opposite the meteoroid radiant. In fact, the anti-radiant line corresponds to a singularity in the gravitational focusing equations, in which one-dimensional rings of initial meteoroid trajectories are focused onto zero-dimensional points along the anti-radiant. As a result, any observer or spacecraft that wanders close to the anti-radiant is predicted to encounter a massively increased meteoroid flux.
Real meteoroid streams will not be perfectly collimated; some dispersion in velocity will always be present. Numerical simulations have demonstrated that a radiant with some "width" results in finite gravitational focusing at all loca-E-mail: althea.moorhead@nasa.gov tions (Jones & Poole 2007). However, an analytical approximation for this effect is needed for efficient modeling of the meteoroid environment. In this paper, we derive a modified analytical treatment of gravitational focusing and planetary shielding that incorporates smoothing terms and we tie these modifications to observed dispersions in meteor shower radiants.
There is a second singularity in the gravitational focusing equations: the flux or number density enhancement due to gravitational focusing tends to infinity as the meteoroids' initial speed relative to the massive body approaches zero. This behavior was noted byÖpik (1951) and is a result of the assumptions inherent in most treatments of gravitational focusing. These assumptions are that there are no other gravitational bodies present, that the stream of meteoroids is infinitely wide, and even that the universe is infinitely old. Under these conditions, a single massive body can indeed attract initially stationary meteoroids from arbitrarily large distances, resulting in an infinitely large flux. In this paper, we compute a velocity smoothing term that assumes gravitational focusing is limited by the massive body's Hill radius.
These two smoothing terms solve different meteoroid dynamics problems. Our angular smoothing term results in more realistic meteoroid fluxes near the anti-radiant; relevant scenarios include gravitational focusing of a meteor shower by the Earth at the Moon's location, by the Moon at the Earth's location, and by the Earth at the location of a high-orbiting satellite, such as one in geosynchronous orbit. Our velocity smoothing term results in more realistic fluxes at all angles for very slow meteoroids. This is more likely to be useful when modeling sporadic meteoroids or dust whose orbits have been circularized due to Poynting-Robertson drag and thus encounter the Earth at very low speeds.
One can, of course, directly simulate the motion of meteoroids in response to the Earth's gravity. This approach is accurate but computationally intensive. We thus presume that anyone applying these equations is doing so to reduce run time and so we prioritize numerical efficiency over fidelity.

RADIANT DISPERSIONS
A meteor shower originates from a single body, usually a comet. However, even if the parent body's orbit remains unchanged over time, meteoroids are propelled away from the parent by sublimating gases; this process alone may produce velocity dispersions of nearly 1 km s −1 in some cases (Moorhead 2018). Meteoroid orbits are further modified by planetary perturbations and by size-dependent radiation-related effects such as solar radiation pressure, Poynting-Robertson drag, and the Yarkovsky effect. The resulting dispersions in radiant and speed are then obscured by differential deceleration in the atmosphere and by measurement error (Kresak 1992). The true dispersion in speed is particularly difficult to measure, as precise measurements of the initial speed cannot be obtained without modeling the deceleration of each meteor in the atmosphere (Vida et al. 2018b).
Meteor radiant measurements are less sensitive to atmospheric deceleration and can, if one assumes the velocity dispersion is isotropic, also be used to estimate the speed dispersion (Kresák & Porubčan 1970;. Radiant dispersions are not routinely measured, but some estimates do exist in the literature. Kresák & Porubčan (1970) presented meteor shower radiant dispersions for nine showers using high-quality double-station photographic data, and found that the observed radiant spreads were broadly consistent with a 1 km s −1 velocity dispersion. Betlem et al. (1997) presented meteor radiants from the 1995 Leonid outburst that appear to exhibit a dispersion much less than 1 • .  compared Perseid radiant dispersions in typical and outburst years and found the latter was significantly tighter. Several older studies obtained very tight radiant dispersions (e.g., 6.2' for the 1946 Draconids) based on the intersections of single-station meteor observations (Jacchia et al. 1950;Millman & McKinley 1963). We suspect this method, which averages many orbits, underestimates the true dispersion.
These historical measurements are sometimes limited by measurement error and, in many cases, do not reduce the observed dispersion to a single quantitative measurement. Therefore, in this section, we present real radiant dispersions for two showers observed using the Global Meteor Network and demonstrate that they can be described by a Rayleigh distribution. A more detailed analysis of these data will be the subject of a future paper and the distributions presented here are intended to be illustrative rather than definitive.  1 for their properties and stellar limiting magnitudes. A network of 7 stations with 16 mm lenses is located in western Croatia and all observe the same volume of the sky to ensure high quality meteor trajectory measurements. We will use this 16 mm data as a reference high-quality data set and investigate whether systems with 3.6 mm lenses observe the same shower radiant dispersions. All GMN stations run open source Raspberry Pi Meteor Station (RMS) software 2 . The details of the meteor detection algorithm and the software in general are given in Vida et al. (2016Vida et al. ( , 2018a. The initial astrometric calibration is manually performed for all stations using the GAIA DR2 star catalog (Prusti et al. 2016;Brown et al. 2018). The camera pointing (field of view center, rotation, and scale) is automatically recalibrated on a block of 256 averaged video frames around every meteor detection. This is necessary due to a camera pointing drift of up to 10' per night from thermal expansion and contraction of the camera bracket. Note that the magnitude of the expansion is larger than the plate scale (up to ten times larger for the 16 mm systems), making recalibration an absolute necessity to obtain precise measurements.
The trajectories are computed using the Monte Carlo method described in Vida et al. (2019b). Multi-station observations are paired based on temporal and spatial correlation, and only those solutions which have heights and velocities in the typical meteor range are taken. If a meteor is observed from three or more stations, observations from stations with angular fit residuals larger than 3σ from the mean are rejected and the trajectory is recomputed. From December 2018 until mid-November 2019, more than 35,000 orbits were collected using 77 GMN stations, and 2,500 orbits collected only using 16 mm systems. Trajectories computed using data from 3.6 mm systems have median trajectory fit residuals of about 60", while trajectories from 16 mm sys-1 GMN website: https://globalmeteornetwork.org/ 2 RMS library on GitHub: https://github.com/ CroatianMeteorNetwork/RMS tem have about 10". Vida et al. (2019c) show that in theory systems with 3.6 mm lenses should be able to achieve radiant measurement precision of about 0.1 • , which is much smaller than most modern multi-station estimates of radiant dispersions.
Meteors were associated with showers using the table of Sun-centered ecliptic shower radiant positions given in Jenniskens et al. (2018). This table, which is based on CAMS data, includes radiant drift and provides numerous, scattered sets of reference parameters per shower. For instance, within one degree of solar longitude near the Orionid peak, the table quotes 3894 reference radiants that span more than 15 • in longitude and 11 • in latitude. We assume that meteors within 1 • in solar longitude, 3 • in radiant, and 10% in geocentric velocity of a shower reference location are members of that shower. If there are multiple shower candidates in the vicinity, the shower with the minimum closeness score D sce is taken: where ∆λ is the difference in solar longitude between the meteor and the reference shower radiant, θ is the angular separation, and ∆w is the geocentric velocity (w) difference. Ideally, the 16 mm data would be used to estimate dispersions of all showers, but these stations produce a small fraction of the total number of meteors in the full data set. For instance, the 16 mm data set has only 51 Orionids, while the full data set has a total of 2384 Orionid meteors. To compensate, we apply aggressive quality cuts to the full data set: we keep only those meteors with (1) a convergence angle of at least 15 • , (2) a median fit error no more than 60", (3) an average velocity error no more than 2%, and (4) a total angular radiant error no more than 0.25 • for the Orionids or 0.35 • for the Perseids. After these cuts, 632 Orionids and 566 Perseids remain in the full data set. In section 2.2, we use this filtered data set to characterize the radiant dispersion; we also compare this dispersion with that of the 16 mm data and show that the two data sets produce consistent results.

Dispersion characterization
We present our geocentric meteor radiants in Sun-centered ecliptic longitude (λ g −λ ) and latitude (β g ) to minimize the movement of the shower radiant over time, but even in these coordinates some radiant drift is present. We must account for this drift in order to obtain the instantaneous radiant dispersion that will be relevant to our gravitational focusing calculations. We model the drift as a linear shift in both angular coordinates over time; for the Orionids, this drift is approximately: Subtracting equations 2 and 3 from each Orionid meteor radiant results in a more compact distribution. Fig. 1 presents the position of these radiants relative to the mean and includes their uncertainties as error bars. Note that the uncertainty in any individual meteor radiant is much smaller than the radiant dispersion, and thus measurement uncertainty does not make a significant contribution to the apparent dispersion.
−10 −5 0 5 10 We next computed the angular offset of each radiant from the average; see Fig. 2. We found that the entire distribution was not well-described by a Rayleigh distribution because of the long tail of radiant offsets larger than 1 − 2 • . We expect that this tail is due to contamination by sporadics or nearby showers. We eliminated these contaminants in a crude fashion by simply discarding those meteors that lay outside a cutoff value; we found that when we used a cutoff value of 1.4 • , the remaining Orionid meteors agreed well with a Rayleigh distribution. We therefore iterated our radiant drift fit to exclude meteors lying outside this boundary (equations 2 and 3 include this refinement). We also restricted our data to include only those meteors falling within 15 • of the median solar longitude for the shower.
In Fig. 2, we fit a Rayleigh distribution with a mode of 0.52 • to our 3.6 mm Orionid data. A Kolmogorov-Smirnov (K-S) test fails to reject the hypothesis that the data follow a Rayleigh distribution with a p-value of 0.6. A second K-S test fails to reject the hypothesis that the 16 mm data follow this same distribution (with a p-value of 0.5). We therefore conclude that a Rayleigh distribution is a reasonable description of these data, and that the 3.6 mm and 16 mm data exhibit the same Orionid radiant dispersion. These values are comparable to those of Kresák & Porubčan (1970), who give a median dispersion of 0.71 • for the Orionids and 1.07 • for the Perseids.
We applied a similar set of calculations to a set of Perseid meteors. The Perseids exhibit less of a radiant drift, with a slope of only -0.03 in longitude and 0.14 in latitude. They also display a higher dispersion in their radiant (see Fig. 3), and we obtained the closest resemblance to a Rayleigh distribution in this case by applying a radiant offset cutoff of 2.9 • . Here, the p-value for our Rayleigh distribution test between the 3.6 data and our best fit was 0.3, and the p-value between the 16 mm data and our best fit was 0.6. The mode of our best fit Rayleigh distribution was 1.2 • .
In this section we have demonstrated that, for both the Perseids and Orionids, the core component of the radiant dispersion can be described using a Rayleigh distribution (see appendix A for additional discussion of the speed and radiant distributions). In both cases, the "location" of the  Rayleigh distribution is fixed to 0, and so the only free parameter is the mode or scale of the distribution. Thus, we are now able to describe a meteor shower's radiant distribution with a single value. A more careful subtraction of the sporadic contamination from these radiant dispersions is slated for a future paper. In the meantime, we conclude that our simulated Rayleigh distribution with a mode of just over 1 • is a reasonable illustrative case.

NUMERICAL SIMULATIONS
Now that we have quantified several shower radiant distributions, we next investigate the impact of such a distribution on gravitational focusing. We do so by conducting a series of numerical simulations of meteoroid trajectories near a mas-  sive body. In these simulations, meteoroids approach a massive body on parallel or nearly parallel initial trajectories, and we sample their position at random times in order to determine the number density of meteoroids as a function of position.
We construct four scenarios to explore the effects of meteoroid speed and radiant dispersions on gravitational focusing and planetary shielding. In the first scenario, no dispersions are present: meteoroids are placed on perfectly parallel trajectories far from the massive body. In the last scenario, we include an idealized meteoroid velocity dispersion in which the distribution of each component of the meteoroids' velocity vector is described by a Gaussian with a standard deviation of 0.5 km s −1 . We construct two intermediate scenarios by dividing this velocity dispersion into two components: a speed dispersion and a radiant dispersion (see Fig. 4 and 5). In the former, we preserve only the variation in meteoroid speed and enforce a uniform direction of movement. In the latter case, we preserve only the variation in meteoroid directionality and enforce a constant speed. This approach assumes that the dispersion in meteoroid velocity is isotropic, and that the dispersion in speed can be determined from the dispersion in radiant and vice versa (see appendix A).
In all cases, we convert the initial trajectories of the meteoroids to orbital elements using the methodology of Murray & Dermott (1999), solve Kepler's equation (Gooding & Odell 1988) at random times, and convert the orbital elements back to positions at those times. We then bin these positions in x, y, and z and plot the z = 0 slice of the grid.
A set of simulations is shown in Fig. 6. In these simulations, the meteoroids have an initial speed of 25 km s −1 ; we have selected this speed, which is significantly slower than the geocentric speed of the Orionids or Perseids, in order to obtain a more compact illustration of gravitational focusing. The massive body is assumed to be the Earth and we include 100 km of meteoroid-blocking atmosphere in the Earth's shielding radius. The color scale provides the number density enhancement factor, or ratio of the local meteoroid number density to that at infinity. The massive body itself  and its "shadow" -those locations that are shielded by the massive body from the meteoroid stream -appear in black. When the meteoroid stream has no initial velocity or radiant dispersion (panel E1a of Fig. 6), the number density can be enhanced by two orders of magnitude near the anti-radiant line. The maximum value in our scale is limited only by our bin resolution (which is, in these simulations, 20 bins per R ⊕ ); a finer resolution would produce even higher enhancement values. Introducing a speed dispersion only (panel E3a) causes the "shadow" of the massive body to end less abruptly but does not eliminate the antiradiant singularity. Introducing a radiant dispersion (panel E2a), however, reduces the maximum enhancement factor to a factor of a few. We find that introducing a radiant dispersion only produces a result that is quite similar to that obtained with our full three-dimensional velocity dispersion (panel E4a). From these simulations, we conclude that, for the purpose of calculating gravitational focusing, it is more useful to measure a meteor shower's radiant dispersion than it is to measure its velocity dispersion.
We present a second set of simulations in Fig. 7. In this case, we set the initial speed of the meteoroid stream to 35 km s −1 and compute the number density over a distance that exceeds lunar perigee. We also include simulations in which we have used the Moon, rather than the Earth, as the massive body (panels M1b and M4b). For each massive body, we simulate the gravitational focusing of meteoroids on perfectly parallel paths (panels E1b and M1b) and those with a velocity dispersion of 0.5 km s −1 (and thus both a radiant and speed dispersion; panels E4b and M4b).
This simulation is designed to mimic the Geminid meteoroid stream, which has a geocentric speed of 35 km s −1 and lies near the ecliptic plane. Because of these properties, it has been suggested that the Moon could cause an enhancement of the Geminid flux visible at Earth via gravitational focusing (M. Matney, personal communication). Fig. 7(M1b) does indeed suggest that the Moon could produce a narrow but significant enhancement of the Geminid flux near the Earth. Based on this, one might expect short-lived, intense bursts of meteor shower activity whenever the Moon passes near the shower radiant. However, Fig. 7(M1b) indicates that the Moon and Earth would have to be very closely aligned with the Geminid radiant, while in fact the Moon never passes within five degrees of the radiant. Furthermore, Fig. 7(M4b) shows that a modest velocity dispersion of 0.5 km s −1 (which corresponds to a radiant dispersion of 0.8 • at 35 km s −1 ) obliterates any enhancement, and also substantially shortens the "shadow." Based on these results, we conclude that it is unlikely that the Moon could produce detectable flux enhancements in any meteor shower seen at Earth. However, Fig. 7(E4b) shows that the Earth could produce noticeable enhancements at the Moon's location if the geometry is favorable.
The simulations shown in Fig. 7 also demonstrate that for less massive bodies, we will need to consider the effect of radiant dispersion on planetary shielding as well as on gravitational focusing. A modest radiant dispersion can drastically reduce the maximum distance at which a spacecraft can employ shielding as a sheltering technique, for instance.

ANALYTICAL APPROXIMATIONS
The effects of gravitational focusing and shielding on a parallel stream of meteoroids can be computed analytically (Divine 1992;Staubach et al. 1997;Matney 2002;Jones & Poole 2007). This approach is far more computationally efficient than simulating a large number of individual meteoroid orbits. As a result, those meteoroid environment models that describe the flux near planets employ an analytical approach (Smith 1994;Staubach et al. 1997;Moorhead et al. 2019a). In this section, we concoct modified versions of the Staubach et al. (1997) equations that incorporate the effect of radiant and velocity dispersion while preserving numerical efficiency. The equations for computing gravitational focusing and shielding are scattered throughout this section; for the convenience of the reader, they are also presented in a more compact form in appendix C.

Gravitational focusing
If the velocity of a meteoroid relative to a massive body before entering its gravitational well is ì w, then the vis viva where G is the gravitational constant and M is the mass of the massive body. If we introduce the dimensionless variable f = GM/rw 2 , equation (4) can be re-written as: Note that our f is the inverse of the variable F used by Jones & Poole (2007). Let ξ represent the observer's angular offset from the anti-radiant; then, w r = w cos ξ is the component of the unfocused meteoroid velocity parallel to the line between the massive body and the observer. Using these parameters, we can re-write equation 7 of Staubach et al. (1997) to obtain the enhancement of the meteoroid number density due to gravitational focusing, η: Note that our b is equivalent to the unitless value B of Jones & Poole (2007), which, if multiplied by w, is the "auxiliary velocity" also named B by Staubach et al. (1997). The "±" sign is present in η because there are two paths by which meteoroids with a given radiant can reach a given point in space: one long, one short (see Fig. 8). The upper choice of sign corresponds to the short path, while the lower choice corresponds to the long path. Both must be included to determine the total enhancement.

Removal of the anti-radiant singularity
Putting the equation for gravitational focusing in this form reveals its singularities and opportunities for smoothing over those singularities. For instance, η tends toward infinity as ξ → 0; this results in an infinitely large enhancement along R sh or t pa th lon g pa th r observer w ξ Figure 8. Diagram of a massive body (gray circle), observer location, and the two paths by which meteoroids from a given radiant and with initial speed w can reach the observer.
the anti-radiant, as discussed in § 1. One can remove this singularity by introducing a minimum angle ξ min : η F = 1 2 1 ± sin 2 ξ 2 + sin 2 ξ min 2 + f b , where (8) b = sin 2 ξ 2 + sin 2 ξ min 2 sin 2 ξ 2 + sin 2 ξ min 2 + 2 f This limits the enhancement near the anti-radiant while ensuring that η F η when ξ is large. Including sin 2 ξ min /2 in the numerator of η as well as the denominator ensures that the total focusing factor tends to one when f and ξ are both small, replicating the behavior seen in Fig. 7(M4b).

Removal of the slow meteoroid singularity
The modifications made in equations (8) and (9) remove the anti-radiant singularity in the focusing factor, but a second singularity remains. As the speed w approaches zero, both f and the focusing factor eta approach infinity. This corresponds to an idealized scenario in which the massive body is the only gravitational body in existence. In such a case, the massive body is capable of attracting infinitely distant meteoroids so long as they have no initial velocity with respect to the planet. This is unrealistic for several reasons. First, meteoroid streams are not infinitely broad. Second, the gravitational influence of planets in a solar system is limited to their Hill sphere. In the case of the Earth, the Hill sphere is generally smaller than a meteoroid stream (although there can be exceptions; for instance Egal et al. 2019, predict Draconid filaments narrower than 0.01 au). The Earth's Hill radius corresponds to a "smoothing speed" of 48 m s −1 . This is much smaller than the speed of any meteor shower, but could be used to place an upper limit on the focusing of a highly circularized population of sporadic meteoroids or dust particles.

Planetary shielding
Equation (6) or (8) describe how a massive body's gravity modifies the local meteoroid number density, but, for many geometries, the massive body may block one or both trajectories by which a meteoroid could otherwise reach the observer. This is termed "planetary shielding" and must also be taken into account. A trajectory is blocked if the meteoroid encounters the observer post-periapse and the periapse distance lies inside the massive body's meteoroid-blocking radius. Staubach et al. (1997) provide equations for computing the periapse distance and true anomaly of the meteoroid at the observer's location that can be used to test these conditions. We provide simplified versions of these equations in our unitless terms.
First, rather than calculate the true anomaly, which requires us to perform several numerically expensive square roots and an arcsine, we instead calculate the radial component of the meteoroid's speed at the observer's location and determine its sign: where the upper choice of sign again corresponds to the short path. If equation (10) is true, then the meteoroid is moving away from the massive body at the observer's location and has thus passed through periapse. Second, the periapse distance q lies within the massive body's radius R when If both equations (10) and (11) are satisfied, the path is blocked and makes no contribution to the local number density. We were unable to concoct a simple modification of these equations that mimics the blurring effect that a radiant dispersion has on the planetary shielding pattern. In the case of the Earth, such a modification is probably unnecessary, as the shape of the "shadow" is very similar in the presence or absence of a radiant dispersion. In cases where planetary shielding plays a more prominent role than gravitational focusing, however, some alternative treatment of the shadow is clearly needed (see, e.g., panel M4b of Fig. 7). In this case, we use the following shielding requirement: sin ξ + cos ξ tan ξ min < R r When this equation is satisfied, the observer lies within an isosceles triangle whose two equal sides are tangent to the massive body and where the angle between them is 2ξ min 2ξmin ν0 Figure 9. Diagram of the nominally shielded region (red) and the "dispersion triangle" (blue). Note that the lines in red are curved and thus the nominally shielded region is not a triangle. In this case, the shielded region lies within the dispersion triangle.
(see Fig. 9). We assume that the "shadow" can be no larger than this "dispersion triangle." The repeated evaluation of this additional equation may be undesirable, especially for cases such as that depicted in Fig. 7(M4b), where the shielded region is clearly smaller than the triangle defined by ξ min . In such a case, it may be numerically advantageous to perform one calculation at the outset to determine whether (a) the nominally shielded region has a smaller radial extent than the dispersion triangle, or (b) the nominally shielded region lies entirely within the dispersion triangle.

The shielding-dominated case
Let us use r 0 to denote the largest distance from the center of the massive body at which the observer can experience classical planetary shielding. This distance is given by: where a = −GM/w 2 is the semi-major axis of the meteoroids. Condition (a) is therefore equivalent to: When the above inequality is true, the dispersion triangle lies entirely within the nominal shadowing region. If equation (12) is satisfied, we assume the short path is blocked. Contributions from the long path can be neglected, as η 0 when the lower choice of sign is taken and f is small compared to sin 2 ξ min /2.

The focusing-dominated case
The true anomaly of a meteoroid that grazes the massive body and intersects the observer along the anti-radiant line is given by cos ν 0 = a/(a − R). Thus, the nominally shadowed region lies entirely within the dispersion triangle when In this case, gravitational focusing determines the shape of the shadowed region and it is unnecessary to evaluate equation (12). Moon velocity 35 Table 2. Summary of simulations presented in this paper. Note that "intermediate" refers to a central body with no real-life analog whose mass is 0.11 M ⊕ and radius is 0.57 R ⊕ . A "velocity" dispersion indicates that each vectorial component of ì w follows a Gaussian distribution with σ = 0.5 km s −1 , while a "speed" dispersion indicates that only the magnitude, w, is allowed to vary.

Intermediate cases
Because meteoroids passing near the massive body do not follow straight paths, cases exist in which neither the nominally shielded region nor the dispersion triangle lie entirely within the other. This occurs when: In this case, one must evaluate both equations (10) and (12), as well as equation (11), to determine whether the observer lies within the shielded region. However, given that the radiant dispersion produces a fuzzy boundary, it is not likely to be of much value to determine whether the observer is in the small region of space that lies outside the dispersion triangle but within the nominally shaded region. We therefore suggest treating this case as equivalent to the focusing-dominated case.

Success of approximations
In this section, we compare our approximations (equations (8) and, where applicable, (12)) with our simulations of meteoroid streams with non-zero radiant and speed dispersion. We perform these tests on simulations E4b and M4b from Fig. 7, using the mode of the simulated radiant distribution, 0.8 • , as ξ min . In order to probe the "intermediate" case (section 4.2.3), we include an additional massive body with no real-world analog whose mass is 0.11M ⊕ and whose radius is 0.57R ⊕ . Fig. 10 presents the difference between our numerical simulations and our analytic approximations, expressed as a percentage. Table 2 lists the key differences between our various simulations.
When ξ is large (i.e., when the observer lies far from the anti-radiant line), the simulation and equation agree very well; the only visible differences are due to noise in the numerical simulation. For smaller values of ξ, our equation tends to underestimate the gravitational enhancement in the focusing-dominated and intermediate cases. In the region where focusing is the most intense, we underestimate the degree of focusing by about 25%. We have used the mode of the Rayleigh distribution -0.8 • -as our value of ξ min . We tested other values, such as the square root of the variance of the distribution, but none performed better than the mode. Thus, an approximation that uses only the mode of the (Rayleigh) radiant dispersion is able to describe the pattern of gravitational focusing produced by both this Rayleigh radiant dispersion and a corresponding speed dispersion.
In the regions where one or both trajectories are blocked according to §4.2, we underestimate the number density by up to 100%. This is because meteoroids can "bleed" into these regions in our simulations in a way that our analytic approximations do not capture. However, the physical extent of these regions is small and limited to the border of the shadow.
Overall, we deem our approximations to be in good agreement with our simulations: a 25% error in the number density enhancement is a significant improvement over the many-orders-of-magnitude errors that can occur when neglecting the radiant dispersion entirely.

TheÖpik test
In this paper, we are primarily concerned with reproducing the local gravitational focusing at a single observer's location. However, if one instead considers the total flux incident on the surface of a planet, a much simpler equation applies. Opik (1951) showed that, for a given value of w, gravity increases the flux incident on a planet's surface by a factor of: where H is the Heaviside function and f at the planet's surface is of course GM/Rw 2 . Equation (17) can thus be used to test whether a set of gravitational focusing equations behave as expected. Jones & Poole (2007) were not able to satisfy the "Öpik test" using equation (6); see Fig. 6 of Jones & Poole (2007). This appears to be due to a misunderstanding; equation (6) describes the factor by which the number density of meteoroids is increased due to gravitational focusing. It is not equivalent to the flux enhancement factor described byÖpik (1951); the flux is further enhanced due to the local increase in meteoroid speed. Moorhead et al. (2019a) demonstrate that equation (6) does, in fact, satisfy theÖpik test when applied as a number density enhancement factor.
We do not expect satisfaction of theÖpik test to be much compromised by our smoothing terms. The largest modification to the number density or flux enhancement occurs where ξ is small, and these locations will not experience any incoming flux due to planetary shielding. Nevertheless, as a test, we compute the total flux incident on the top of the Earth's atmosphere (i.e., we compute the sum of −η u r w where r = R ⊕ and u r < 0) for a variety of initial meteoroid speeds and dispersion angles and compare it with the expected value. Except where the initial speed is very low and the dispersion very high, the results satisfy theÖpik test to within 1% (see Fig. 11).

Conservation of flux
Perhaps a better test of the ramifications of our modified gravitational focusing treatment is whether the total flux in the absence of planetary shielding is zero. Fig. 12 displays the results of a series of such tests. In these tests, we place a sphere around the Earth and compute the total flux on the surface of this sphere: If we were to integrate η rather than η , we would obtain a total flux of zero. However, since the number density of material passing out of the sphere at low ξ tends to be artificially reduced by our smoothing factor, we expect our total to be slightly negative. To put the flux imbalance into perspective, we divide by the total incoming flux and express this fraction as a percentage. We perform these calculations only for distances lying outside the planet's shadow; inside the shadow, low-ξ locations are blocked and this test is not relevant. Fig. 12 shows that when the the incoming speed is low, up to 5% of the flux can be lost at certain distances. However, for typical shower speeds -those greater than 20 km s −1 -the discrepancy is smaller. We have assumed a radiant dispersion of 5 • in these simulations in order to exaggerate this effect; for a more typical dispersion of 1 • , the disagreement is at most 1% (see Fig. 13). Thus, at least in the case of the Earth, our analytic approximations are unlikely to introduce large errors in either the local number density or flux, the total flux incident on the Earth, or the total flux moving through a volume.

Apparent radiant
Our modified version of the Staubach et al. (1997) approach still has only two paths by which meteoroids may reach the observer. However, a diffuse radiant permits an observer to "see" meteoroids coming from a ring or arc of apparent radiants. Figure 14 shows the distribution of apparent radiants as seen at three observer locations in simulation E4b. Locations (a) and (b) lie in the region in which only one nominal trajectory is predicted to reach the observer, while location (c) lies very close to the anti-radiant line and is predicted to encounter meteoroids from both the "long" and "short" paths.
When the observer lies far from the anti-radiant, the "short" path is a reasonable description of the apparent meteoroid directionality. The vast majority of apparent radiants lie close to the predicted encounter angle. A few meteoroids appear to come from the "long" path, but this is a small fraction of the total. The agreement worsens as we approach the anti-radiant. Location (b) lies near the border of the region in which both trajectories are not blocked, but the long path is still predicted to be blocked. As a result, a substantial fraction of the meteoroids appear to lie near the long path. At locations even closer to the anti-radiant, such as (c), both trajectories are predicted to contribute, but the apparent radiant distribution resembles a ring rather than two points or arcs.
Thus, our approach is not capable of describing the extended directionality of meteoroids with non-zero radiant dispersion intersecting observer locations near the antiradiant line. A full description of the directionality can likely only be obtained by a more computationally intensive approach that considers the full distribution of meteoroid orbits.

APPLICATIONS
In this section, we discuss the ramifications of our study for several different scenarios.

Meteoroid environment models
We anticipate that the primary application of our numerical approximations will be in environment models. For instance, the NASA Meteoroid Environment Office (MEO) issues meteor shower forecasts (Moorhead et al. 2019b) that have been recently updated to incorporate the Staubach et al. (1997) treatment of gravitational focusing. We have found that spacecraft at high altitudes, such as those in geostationary orbit, often pass near the anti-radiant of active meteor showers and were thus predicted to encounter brief but intense spikes in the impactor flux. In this paper, we have shown that these spikes are not realistic, given our best estimates of meteor shower radiant dispersions. The modified algorithms we have derived here will be incorporated into future spacecraft-specific forecasts.
Our algorithms may also be useful for sporadic meteoroid models. For example, the MEO also generates a piece of software called the Meteoroid Engineering Model (MEM Moorhead et al. 2019a) that describes the sporadic meteoroid flux encountered along a spacecraft trajectory. MEM models the environment as a collection of meteoroid orbits that contribute to the local number density. If the velocity corresponding to one of these orbits is very closely parallel to the separation between the Earth and spacecraft, the anti-radiant singularity could produce what appears to be a "hot pixel" in the resulting flux map. Assigning some angular width to the apparent radiant would eliminate this possibility.
Neither MEM nor the MEO's shower forecasts model meteoroids whose speed relative to the massive body at infinity is less than 1 km s −1 . However, models of asteroidal material or dust particles may include particles on highly circularized orbits that encounter the Earth at extremely slow speeds. Several such studies (Wetherill & Cox 1985;Kortenkamp 2013;Pokorný et al. 2019) have noted that gravitational focusing enhancement of these particles can be extreme. For instance, Kortenkamp (2013) noted that low-inclination, quasi-satellite dust particles have enhance-ment factors of 3000. To avoid this, Pokorný et al. (2019) applied a softening parameter of 0.1 km s −1 toÖpik's relation. The minimum speed we derived from the Hill radius -0.048 km s −1 -is only a factor of two smaller and thus supports the decision made by Pokorný et al. (2019) to smooth over this velocity scale.

Gravitational enhancement of showers by the Earth or Moon
As discussed in section 3, it has been hypothesized that gravitational focusing by the Moon could produce enhanced meteor shower activity at the Earth (M. Matney, personal communication) or by the Earth at the Moon (see Fig. 7). We discussed the case of the Geminids, but, to give another example, in some years the Moon passes within a few degrees of the kappa and lambda Virginid (KVI and LVI) anti-radiants.
In the absence of any radiant dispersion, such an alignment would result in brief, intense enhancements of these meteor showers. Fig. 7(M4b) reveals that the Moon is highly unlikely to produce any observable modification of meteor rates at Earth. However, Fig. 7(E4b) indicates that gravitational focusing of a stream by the Earth could enhance the rate of meteoroids striking the Moon by a factor of a few. The meteoroid flux at the Moon can be measured via impact flashes (Suggs et al. 2014;Madiedo et al. 2014;Bonanos et al. 2016); determining whether and when such enhancements are visible could be worthy of further study.

Second-order focusing in the Earth-Moon system
Any meteoroid within the Moon's Hill sphere is also necessarily within the Earth's Hill sphere. To calculate the trajectory of such a particle with absolute accuracy, it is therefore necessary to compute the effects of both the Earth and the Moon. Similarly, a stream of meteoroids will experience second-order gravitational focusing in which the flux and number density are modified by both massive bodies. However, given the manner in which a modest radiant dispersion reduces the gravitational enhancement produced by the Earth at the Moon's location, and given also that scattering by the Earth will tend to increase dispersion in the meteoroid stream, we conclude that the effects of secondorder focusing can be ignored.

CONCLUSIONS
While the effect of a massive body on a meteoroid stream can be described analytically, these equations contains a singularity that results in an arbitrarily large enhancements along the anti-radiant line. An additional singularity causes similarly large enhancements as the initial speed of the meteoroids approaches zero. A real meteoroid stream, however, will exhibit a dispersion in its direction of motion that eliminates these singularities.
We found that meteoroid radiant dispersions can be described as a Rayleigh distribution and that the mode of this distribution can thus be used to describe the level of dispersion. We measured the dispersion of two showers and found that the Orionids had a radiant dispersion of 0.5 • and the Perseids 1.2 • . In the absence of measurements, we suggest assuming radiant dispersions of roughly 0.5-1 • . Additional radiant dispersion measurements are the subject of a future paper. Our approach assumes that the dispersion in the velocity vector of the stream is isotropic. If the dispersion in the transverse component of the meteoroids' speeds is larger or smaller than the dispersion in the speed in the direction of the stream's motion, or if the radiant distribution is not symmetric, this assumption is not valid.
We simulated the effect of modest radiant (and speed) dispersions on the behavior of a stream of meteoroids passing near a massive body. We assumed a 0.5 km s −1 dispersion in speed, corresponding to a Rayleigh radiant dispersion with a mode of 1.15 • for an average speed of 25 km s −1 or 0.8 • for 35 km s −1 . In all cases, we found that including a dispersion in radiant substantially reduced the maximum local number density enhancement. We also found that it limited the length of the region shielded from the stream by the massive body.
We derived modifications of the Staubach et al. (1997) treatment of gravitational focusing and planetary shielding that approximate the effects of a radiant dispersion. These approximations, which are summarized in appendix C, agree with our simulations to within 25%. We found that the mode of the Rayleigh radiant distribution worked well as a smoothing factor, which we call ξ min . An approximation using only the mode of the radiant distribution was able to describe the pattern of gravitational focusing produced by meteoroids with dispersions in both their speeds and radiants. Tests of our approximations showed that our modified gravitational focusing equations satisfy the so-calledÖpik test to within 1-2%, and conserve flux to within a few percent.
We anticipate that our algorithms could be useful in meteoroid environment modeling, where an efficient yet accurate description of the flux near a massive body is desired. Shower models would benefit from individual measurements of radiant dispersion; we plan to conduct a survey of meteor shower radiant dispersions in the near future.
Finally, the authors thank an anonymous reviewer for independently reproducing and checking our modeling results and for prompting us to include sections 4.3.3 and appendix A.
Step 13a. If the scenario is not shielding-dominated, determine whether the meteoroid trajectory in question is blocked by the massive body as follows: If both of the above statements are true, the trajectory is blocked. (Note that if r < R, the trajectory is also blocked.) Step 13b. If the scenario is shielding-dominated, determine whether the meteoroid trajectory in question is blocked by the massive body as follows: sin ξ + cos ξ tan ξ min < p (C16) If the above statement is true, the trajectory is blocked.
(Note that if r < R, the trajectory is also blocked.) Step 14. If the trajectory is blocked, η = 0. Otherwise, b = (s + s )(s + s + 2 f ) (C17) Sum both possible values of η to obtain the total number density enhancement factor.
Step 15. If desired, compute the velocity of the meteoroid relative to the observer.
r 2 w 1 − cos 2 ξ (C20) Note that the two paths will produce two separate apparent radiants at the observer's location.
Step 16. If desired, compute the meteoroid flux enhancement factor: As with the number density enhancement factor, sum the two possible values of ζ to obtain the total flux enhancement factor.
This paper has been typeset from a T E X/L A T E X file prepared by the author.