Formation of the Lunar Dust Ejecta Cloud

The Lunar Dust Experiment (LDEX) on board the Lunar Atmosphere and Dust Environment Explorer (LADEE) mission orbited the Moon from 2014 September to 2015 April and observed a dynamic, permanently present dust cloud produced by continual meteoroid bombardment. For the latitudes observed by LDEX, the sporadic background contribution to the impacting dust flux is dominated by helion (HE), apex (AP), and antihelion (AH) sources oscillating with lunar phase. Using improved impact ejecta distributions, a three-dimensional model was implemented to estimate the inner and outer ejecta cone angles from LDEX plume measurements. Expanding upon this single-plume model and the derived ejecta cone angles, we implemented a global lunar model fitted to LDEX measurements of the sporadic background to constrain the product of meteoroid impactor fluxes and ejecta mass yield per source. We use the observed asymmetry between sunward and antisunward impactor fluxes to discuss the possible contributions of β-meteoroids, in addition to the HE and AH sporadic meteoroid sources. We find that if β-meteoroids are responsible for the day/night asymmetry, they must have an impact ejecta yield of at least 103.


Introduction
Impact bombardment is a primary source of space weathering for all airless bodies in the solar system, affecting both the topology and chemical distributions of their surfaces over time (Pieters & Noble 2016). Characterizing these meteoroid impactor sources is essential for understanding the origin and evolution of airless bodies within our solar system (DeMeo et al. 2009;Brunetto et al. 2015). The Moon, as an airless body, is of particular interest, as measurements of impact ejecta from the lunar surface serve to help characterize the meteoroid environment near Earth and guide model predictions for the dust environment and regolith evolution of asteroids.
Observations of the lunar ejecta environment by the Lunar Dust Experiment (LDEX) on board the Lunar Atmosphere and Dust Environment Explorer (LADEE) mission discovered a permanently present dust cloud produced by continual meteoroid bombardment . This dust cloud is highly sensitive to meteoroid showers and changes in impactor flux rates (Szalay et al. 2016;Szalay & Horányi 2016a). Of the sporadic meteoroid sources, impact ejecta was found to be produced by the helion (HE), apex (AP), and antihelion (AH) sources with a small contribution from the antiapex (AA) source and possible contributions from northern (NT) and southern (ST) toroidal sources, each oscillating with lunar phase , 2016aSzalay et al. 2019). The fluxes of these four sporadic sources can be constrained by LDEX observations through a forward-modeling approach. The β-meteoroids, dust grains on hyperbolic orbits moving away from the Sun, have been suggested to be an additional source of impact ejecta at the Moon due to the persistent dust ejecta density enhancement on the Moon's sunward side (Szalay et al. 2020a).
Constraints on mass yields, the ratio of the total mass of the produced dust ejecta over the mass of the initial impacting particle, are necessary for deriving impactor fluxes from modeled ejecta clouds. While the relationship between impactor mass, speed, and mass yield has been measured in a variety of experiments (e.g., Koschny & Grün 2001a, 2001b, these do not cover the relevant impactor compositions, sizes, speeds, or characteristics of the surface regolith for the lunar environment. It is the goal of this paper to constrain the product of impactor flux and mass yield per source via a forwardmodeling approach expanding a single-plume model (Bernardoni et al. 2019) to a global one. We also explore the possibility of β-meteoroids producing observable impact ejecta (Szalay et al. 2020a).

Particle Discretization
Ejecta for each impact plume were modeled discretely (Bernardoni et al. 2019), describing the initial ejecta velocity (v), mass (m), and angle from surface normal (f) distributions as Szalay et al. 2016;Szalay & Horányi 2016a;Bernardoni et al. 2019) 2.4 e km s −1 is the lunar escape speed, R m = 1737 km is the Moon's radius, and λ=200 km is the characteristic scale height derived from the dust density distribution measurements (Szalay & Horányi 2016a). The mass exponent (α=0.91±0.003) was measured by LDEX consistent across all local times , and the mass bounds ( =m 10 min 17 and =m 10 max 8 kg) were determined by the detection limits of LDEX at the spacecraft's apex velocity. The outer and inner ejecta cone angles, f 0 and f 1 (see Figure 1), were fixed at 8°and 0°, respectively (Bernardoni et al. 2019). While Sachse et al. (2015) suggested coupled speed and size ejecta distributions, the ejecta distributions used here are treated as separable, as there is a lack of altitude dependence in the detected impact charge distribution of LDEX (Bernardoni et al. 2019;Horányi et al. 2015). For simulating the contribution of a single plume's ejecta to LDEX's measured ejecta impact rate, we sampled 200,000 simulated ejecta particle trajectories based on f (f) and f (v). Positions and velocities for each simulated particle trajectory were calculated at fixed time steps of dt=1 s using elliptical orbit equations. As LDEX's lowest detectable mass threshold is velocity-dependent, f (m) was implemented as a weight for the simulated impact rate from ejecta into LDEX, I i , for impact i (Szalay et al. 2016), where v sc is the spacecraft speed with x chosen in the direction of the spacecraft velocity vector, ω is the angle between the detector's boresight and the ejecta's relative velocity vector v rel , dV is the simulated volume element, and A is the detector's effective area, which is maximal at ω=0 and zero for w   68 . Each plume in this model is considered to have these same initial ejecta distributions. The effects of impactor mass, speed, and obliquity of the impact are contained entirely within the mass yield (number of ejecta per plume) and the impactor flux (number of plumes per spatial-temporal bin). This is done for simplicity, as such effects are not sufficiently constrained over the relevant impact parameters. Formulating these additional effects, particularly asymmetric plumes produced from oblique impacts, could lead to a possible improvement of this model.

One Bin
For simulating the contribution of multiple plumes to the detected impact rate, we first consider the simplified case of a single volume element, dV, with constant impactor flux, F imp . For a single simulated particle at time t i and azimuthal angle f i with respect to the plume's origin and spacecraft velocity (Figure 1), the contribution to the simulated impact rate, I i , is given by Equation (2). Thus, the measured average impact rate of LDEX, I LDEX , is where N p is the total number of plumes produced per simulated time step within the lunar surface element (dA), M p is the average total mass of one plume, and M s is the simulated mass of one plume. The simulated mass of one plume consists of the sum of all particles in the simulation multiplied by the average mass, based on the distribution from Equation (1), while we treat the total mass of an actual plume as an unknown quantity. Note that there are some additional complications to the total simulated mass, as the azimuthal component is treated as continuous while the initial ejecta distributions are sampled discretely. Since the impactor mass flux for characteristic impactor mass m imp can be written as F = m N dA p imp imp , we  where h is the altitude of the spacecraft, = Y M m p imp is the mass yield of a plume, and dz is the altitude bin size (Figure 1). Matching the simulation to LDEX data is then a matter of scaling F Y imp .

Multiple Bins
The previous case holds true for a constant impactor flux within the lunar surface patch considered. As the lunar impactor flux varies with local time, we expand this approach by considering local time bins, assuming that each bin encompasses a constant impactor flux. To standardize the notation of the following sections, i denotes a simulated impact sampled from ejecta particle trajectories as described in Section 2. A subscript of k denotes the observing surfacealtitude bin (i.e., the bin that contains the measured impact rate), while j denotes the surface bin being observed (i.e., surface patch possibly containing the plume origin, corresponding to simulated impact i). Thus, the subscript of dA ijk is read as "for simulated ejecta measurement i originating from surface bin j as measured within surface-altitude bin k." One subtlety that was glossed over was the location of dA relative to dV, as dA is only located directly beneath dV in the case where the angular distance of the simulated impact from the center of the plume q = 0 i . For a single simulated impact at time t i and position q f , i i ( ) relative to the center of the plume, dA is shifted from surface patch corresponding to dV. Thus, for a particular q f , i i ( ), the corresponding dA will overlap up to four local time bins, as shown in Figure 2.
For the illustrated case in Figure 2, dA projected from bin k overlaps the fixed surface bins j 1 , j 2 , j 3 , and j 4 with surface areas dA ij k 1 , dA ij k 2 , dA ij k 3 , and dA ij k 4 , respectively. This simulated impact contribution to the measured impact rate is Illustration of the possible plume impact location for a given simulated ejecta measurement i relative to the surface patch dA directly below the simulated volume element dV. This possible plume location is characterized by the projection of surface element dA given that the simulated measurement traveled a great circle distance of q i at an angle f i with respect to the spacecraft velocity. Dashed lines represent local time-latitude bins, each with a constant flux. For the illustrated case, the simulated impact i has contributions from bins j 1 , j 2 , j 3 , and j 4 weighted by the area ratio in Equation (5). Note that k is also binned by altitude with the corresponding bin determined by the height of simulated ejecta impact i.
where F j 1 , F j 2 , F j 3 , and F j 4 are the impactor fluxes for local time bins j 1 , j 2 , j 3 , and j 4 , respectively. Note that it is with the assumed constant flux that F = m N dA j ijk ijk imp for each i and k. Expanding this to all j bins and summing over all simulated impacts i gives the generalized form of Equation (4), Matching all surface and altitude bins of data becomes a minimization with YΦ j as parameters. In this formulation, the values of S jk contain all of the geometric information from the numerical simulation separated from the desired parameters to be estimated. Specifically, R i contains the geometry of the instrument, while the rest of S jk contains the geometry of the system.

Grid Motivation
While the latitude coverage of the LDEX data set is limited (Figure 4), expanding this model to a full three-dimensional sphere may prove useful for future data sets, as well as to minimize errors introduced by assuming a purely equatorial orbit. In the case of a grid covering the entire surface of a sphere, Equation (6) still holds true with different values for dA ijk . However, complications arise with the increasing number of parameters and computation time. For example, surface grid cells defined by equal latitude and longitude steps are no longer of equal area, as was the case in Figure 2, and whose overlap with arbitrary q i and f i proves computationally expensive to perform for each simulated impact.
Polygonal grids also introduce an additional parameter in the relative orientation of the projected cell to the considered cell. For example, in the case of Figure 2, dA ij k 1 would have a different value from the one illustrated if the projected square was tilted by an arbitrary angle.
Here we use overlapping spherical caps as surface grid cells with grid points defined as the center of these spherical caps ( Figure 3). As this will introduce error in the form of surface patches not covered or double-counted, we construct the grid cell size such that the sum of all surface grid cell areas equals the surface area of the sphere. This represents our best estimate, as the double-counted regions should come close to accounting for the regions not covered by the cells. Additionally, surface grid cells are positioned to meet the following criteria: 1. maximize the smallest distance between any two grid points and 2. minimize the number of unique grid point-to-grid point distances within an effective range.
The first criterion is imposed to achieve a close-to-uniform covering of the sphere, while the second criterion is imposed to reduce the computation time of the simulation. As the distance between grid points will become a parameter in determining S jk (Appendix A, Equations (A2)-(A3)), reducing the number of unique grid point-to-grid point distances significantly increases the performance of the simulation. The "effective range" mentioned in the criteria (and Appendix B, Figure 9) refers to surface grid cells within the maximum q i of all simulated impacts plus twice the radius of the surface grid cells, as any cell beyond this range does not contribute to the simulated impact rate. We used the HEALPix framework (Gorski et al. 2005) to generate the locations of these grid cells, as it has a reasonable agreement with our criteria and the flexibility to scale the resolution as needed. To balance accuracy with computation time, we chose a resolution of = N 2 side 5 , resulting in 12,288 surface grid cells.

Minimization
With the S jk values simulated for Equation (6) over all relevant surface-altitude bins, we can fit the measured LDEX data I k LDEX, using the parameters YΦ j . Starting with a leastsquares approach, we minimize With no additional constraints, n is the number of surface grid cells. To directly connect to physical quantities, we will use the four dominant impactor sources: HE, AH, AP, and AA. These are equatorial sources at local times of 10.3, 1.7, 6, and 18 LT, respectively . Each bin YΦ j is broken up by source α, In this form, the minimization is now a linear regression without a constant bias, F a Y are the regression coefficients, and x αk are the independent variables. Since the dependence of Y on the obliqueness of an impact is described by x αk , Y now represents the mass yield for impacts with normal incidence to the surface.

Results
Using the LDEX impact rate measurements with the coverage shown in Figure 4, fitted values for F Y were calculated for the three equatorial source radiants 10.3, 1.7, and 6 LT, corresponding to HE, AH, and AP, respectively. These sources are labeled by their radiants in Figure 5, as these values represent all possible sources originating from that direction, not just HE, AH, and AP. We excluded data within 5 days of the Geminid and Quadrantid meteor showers to isolate the contribution from sporadic background sources (Szalay & Horányi 2016b). The segment of data coverage in Figure 4 from 12 to 22 LT was excluded from the fit, as the Sun was within LDEX's field of view generating elevated noise levels in the instrument. As this excluded range covers the majority of the 18 LT (AA) source's contribution to the lunar ejecta environment, the AA source was not included in this fit. Similarly, nonequatorial sources, such as the NT/ST, that may also contribute to ejecta production (Szalay et al. 2019) were not included, as the latitude coverage of Figure 4 is not sufficient to constrain them. Using the ejecta mass production weights from this latitude coverage for AP and NT/ST sources, where 18°is the average latitude of all LDEX measurements taken near 6 LT. While we assume that the entirety of the 6 LT source is from AP meteoroids for the following discussion, note that up to 22% may be from toroidal sources. Note that the fitted values for F Y are highly dependent upon m max , partially from the weight in Equation (2) but primarily from M s in Equations (3)-(6). For this reason, we include estimates for =m 10 max 8 kg (opaque points in Figure 5) corresponding to 100 μm ejecta taken as the peak in the sporadic background mass distribution and =m 10 max 11 kg (transparent points in Figure 5) corresponding to 10 μm ejecta and the largest ejecta particle detected by LDEX. The error bars for Figure 5 represent the standard error of that parameter for the linear regression. As shown in Figures 5 and 6, LT (AP) is the dominant contributing source to the lunar dust environment for both cases with a mission-averaged impact ejecta mass production (mass yield times impactor mass flux) of ( ´-  Whether the HE/AH asymmetry is an enhancement on the dayside of the moon or a deficit on the nightside is unclear. However, as we only have possible physical explanations for an enhancement on the dayside, we will proceed under that assumption. For comparison, the following section introduces the empirical formula fit from Pokorný et al. (2019) into our model. It should be noted that while a thermal dependence on the soil's mass yield has been suggested as a potential cause for this enhancement, no physical mechanism has been proposed to link the temperature of the lunar surface to its impact ejecta yield, nor has there been any experimental evidence to suggest that this is the case. As a possible physical explanation, we discuss the potential contribution from β-meteoroids as an additional impactor source in Section 4.2 (Szalay et al. 2020a).

Yield Variance
Our understanding of how the lunar regolith responds to meteoroid bombardment is limited. Effects such as surface temperature, UV radiation, or solar wind have been suggested as possible dependencies that may enhance ejecta mass production on the lunar dayside (Janches et al. 2018; Figure 5. Values for the mass yield times impactor mass flux averaged over the entire lunar surface (total lunar ejecta mass flux) per sporadic background source radiant fitted to LDEX measurements binned in synodic day increments. These values are produced from scaling F a Y ( ) in Equation (9) to the Moon's cross section using =m 10 max 8 kg (opaque points) and =m 10 max 11 kg (transparent points). The opaque and transparent points are staggered for clarity. Measurements only include those with quality flags of 3 or greater, excluding data within 5 days of the identified Geminid and Quadrantid meteor showers and values between 12 and 22 LT. Error bars for each fit represent the standard error for each parameter in the linear regression. Figure 6. Values for the mass yield times impactor mass flux averaged over the entire lunar surface (total lunar ejecta mass flux) per sporadic background source derived from Figure 5 with the additional contribution of Equation (11). Note that the yield in these fitted values should be treated as the lunar nightside yield, as the enhancement on the dayside is pulled into a S j in Equation (9). Pokorný et al. 2019). An empirical relation for the "excess" dayside lunar ejecta was found by comparing to expectations from dynamical models . To incorporate this effect within our model, we introduce the following term into a S j of Equation (9) With this additional local time dependence on the mass yield, we reproduce Figure 5 for =m 10 max 8 kg as Figure 6. Note that the fitted values plotted are of the Y on the lunar nightside, as the enhancement on the dayside is absorbed into a S j . Under this empirical formula, HE and AH sources have a much closer agreement with mission-averaged values of 5.9±0.1 t day −1 for the HE source and 5.0±0.1 t day −1 for the AH source. It should be noted that this formula is derived from the observed HE/AH asymmetry in lunar ejecta and not from surface response relations.

β-meteoroids
As a possible physical explanation for the excess dayside impact ejecta, we consider the contribution of an additional impactor source, β-meteoroids (Szalay et al. 2020a). Since β-meteoroids are of a similar radiant to the HE source at 11 LT (Szalay et al. 2020a), we would expect an enhancement of the total lunar ejecta mass flux on that side. As such, we label the HE/AH contribution as the lower of the two in Figure 7. We cannot take the difference between the total ejecta mass flux of 10.3 and 1.7 LT fitted values as the β-meteoroid contribution, as the mass distribution of ejecta may not be the same as that of Equation (1). However, when considering the ejecta size distributions per local time bin, the values for α remain consistent, with no significant deviations from the fit. For this reason, we use the same ejecta mass distribution for the β-meteoroid source as the sporadic background sources, with only the value for m max changing, as we do not expect an impactor to produce ejecta larger than itself. To provide an estimate range, we consider the cases of  Figure 7, values for the β-meteoroid total lunar ejecta mass flux are similar for the =m 10 kg, max 8´case 3 10 6 ( difference). Based on a variety of spacecraft measurements, the total number flux of β-meteoroids is expected to range from 10 to 600 km −2 s −1 (Berg & Grün 1973;Wehry & Mann 1999;Zaslavsky et al. 2012;Szalay et al. 2020b). For the size range of -0.5 0.2 μm with a bulk density of 2 g cm −3 , the expected impactor mass flux for the β-meteoroids is Φ β =2-300 g day -1 . Using the mission-averaged values from Figure 7, this corresponds to mass yields at a normal incident of =b Y 2.3 -10 3.4 10 3 5 using =m 10 max 14 kg for 0.5 and 0.2 μm β-meteoroids and =´-b Y 1.6 10 2.3 10 3 5 using = m max -10 15 kg. Extrapolating from experimentally measured yields to impactor speeds of 100 km s −1 , we expect high yields for β-meteoroids of -10 10 3 4 , consistent with our results (Koschny & Grün 2001a;Szalay et al. 2020a). Note that while β-meteoroids appear as a possible solution to the HE/AH ejecta production asymmetry, it is by no means certain under this cursory comparison. Further work is necessary to narrow down what effects are contributing to the observed asymmetry.

Conclusion
Expanding upon the description of a single plume (Bernardoni et al. 2019), we constructed a global multiplume model and constrained the product of the mass yield and impactor flux for the three dominant equatorial sources of the sporadic meteoroids bombarding the Moon: AP, HE, and AH ( Figure 5). Further latitude coverage or constraints relative to these equatorial sources are necessary to introduce additional sources into the model, such as NT/ST. Asymmetry between HE and AH sources persists despite considerations for spacecraft trajectory. Considering the cases of surface variation in ejecta mass yield ) and a β-meteoroid source (Szalay et al. 2020a), we constrained the values for the mass yield and impactor mass flux kg (transparent points) for the β-meteoroid source. The opaque and transparent points are staggered for clarity. As β-meteoroids are of a similar radiant as the HE source, the F Y for the β-meteoroids was estimated by requiring that the total lunar ejecta mass flux for the HE source equal that of the AH source derived from Figure 5. required (Figures 6 and 7). Restrictions on the required mass yield extrapolated from the experiment for β-meteoroid ejecta are consistent with our results for the expected impactor mass flux with ejecta mass yields on the order of -10 10 3 5 . This represents a significant increase in mass yield compared with values on the order of 1-10 for the sporadic background sources, which have average impactor sizes much larger than β-meteoroids. The larger inferred yields for β-meteoroids may be a combination of their higher impact speeds and the fact that, for smaller impactors, the surface will appear more solid and less regolith in nature. A similar mission to LADEE carrying an LDEX-type dust instrument on a near-polar orbit could provide the missing observations to precisely describe the ejecta production on the lunar surface. This could be critical to gauge the loss and accumulation of volatiles in the polar permanently shadowed regions, as well as to assess the effects of interplanetary dust bombardment on the evolution of the surface regolith of other airless bodies near 1 au. This work was supported by NASA Headquarters under the NASA Earth and Space Science Fellowship Program, grant 80NSSC17K0471, and NASA's Lunar Data Analysis Program (LDAP), grant 80NSSC17K0702. Support is also acknowledged from the Institute for Modeling Plasma, Atmospheres, and Cosmic Dust of NASA's Solar System Exploration Research Virtual Institute. The LDEX data are available through NASA's Planetary Data System, https://sbn.psi.edu/pds/resource/ldex.html.

Appendix A Grid Derivation
To achieve the area described in Section 3.3.1, the angular radius (r c ) of the spherical caps is given as grid is the number of grid cells on the surface. To calculate dA ijk in Equation (6), we use the following equation for the overlap between two spherical caps of radius r c separated by a great circle distance between their centers of d ijk (Tovchigrechko & Vakser 2001), Note that for this section, k denotes the surface bin directly below the surface-altitude bin under consideration. Additionally, the great circle distance between the projected area of a simulated impact q f , i i ( ) and surface grid cell j is given by the