© ESO 2022

1 Introduction

In the past, studies on the configuration of circumplanetary dust mainly focused on the grain particles around the giant planets and Mars (e.g., Tiscareno & Murray 2018; Spahn et al. 2019). The dynamics of dust particles in the Earth-Moon system has been less explored. Although the Munich Dust Counter did not detect the lunar ejecta dust cloud (Iglseder et al. 1996), a permanent and asymmetric dust cloud engulfing the Moon was indeed detected by the Lunar Dust Experiment (LDEX; Horányi et al. 2015). From the LDEX measurements, the size distribution of the dust grains was derived, and the density of particles was found to drop with altitude and to vary azimuthally around the Moon, reaching a peak at 5–7 h lunar local time. Six previously identified meteoroid populations, helion, antihelion, apex, antiapex, northern toroidal, and southern toroidal, were found to be plausible interplanetary projectiles that generate the lunar dust cloud via high-speed impacts (Szalay et al. 2019). The asymmetric nature of this dust cloud, the result of the asymmetric impactor flux, was analyzed, and a consistent mass production rate per unit surface as a function of impact direction was derived. Szalay & Horányi (2016) also analyzed particles that fall back to the lunar surface, which follow a cumulative size distribution with exponent 2.7. The dynamics of dust particles ejected from the Moon was previously studied by Colombo et al. (1966), and dust orbiting the near-Earth environment was studied by Peale (1966).

Investigating the lunar dust distributions at different times, heights, and positions was ever one of optional tasks of Chinese Lunar Exploration Mission, Chang’e 4 (Wang & Liu 2016). At the end of 2020, Chang’e 5 launched successfully. The dust distribution can be obtained from dust dynamics modeling; as such, the total dust flux received by the spacecraft can be calculated, and the dust hazard for a specific mission can be evaluated. Space activities in the Earth-Moon system are much more frequent compared with other planets, and thus the study of dust particles is of great significance for assessing the space environment and ensuring the security of explorations.

In this work, we focus on the dynamics and distribution of particles ejected from the surface of the Moon that escape the Moon’s gravity. The paper is structured as follows. In Sect. 2 the dynamical model for the motion of particles ejected from the lunar surface is presented, taking the solar radiation pressure and gravity of the Sun, the Earth, and the Moon into consideration. In Sect. 3 the mass production rate and the initial ejecta distribution on the surface of the Moon are deduced. In Sect. 4 the detailed simulation results are presented, with the distribution of sinks and the average life spans of particles as functions of size. The dust number density distribution in the Earth-centered inertial frame is also presented. Finally, the distributions and the evolution of osculating orbital elements of dust particles are given.

2 Dynamical model

Dust particles in interplanetary and circumplanetary space can be affected by various forces, including solar gravity, solar radiation pressure, Poynting-Robertson drag, the Lorentz force, and the gravity of the planets and other bodies in the Solar System. For particles ejected from the surface of the Moon, solar radiation pressure and the gravity of the Moon, Earth, and Sun are considered in this paper. We estimated the strength of the perturbation induced by the Lorentz force for typical strengths of the magnetic field and grain charges in the interplanetary magnetic field, the Earth’s plasma sphere, and the geomagnetic tail of the Earth. We find that outside 10 R_E the strength of the Lorentz force is small compared to the perturbations induced by solar radiation pressure, by at least two orders of magnitude for grains larger than 10 μm; even for grains of 1 μm, solar radiation pressure dominates by more than an order of magnitude. An Earth-centered J 2000 inertial frame Oxyz is used in our simulation. Here, the x axis points in the direction of vernal equinox at the J2000 epoch, the z axis is the normal of the Earth’s equatorial plane (north), and the y axis is determined by the right-handed rule.

Physical parameters for the Sun, Earth, and Moon used in this paper are shown in Table 1. All the ephemerides of the Sun, Earth, and Moon in our model are taken from NAIF SPICE toolkit¹.

The equation of motion of one particle ejected from the surface of the Moon reads $$\stackrel{\mathrm{..}}{r}={\stackrel{\mathrm{..}}{r}}_{{\text{G}}_{\text{E}}}+{\stackrel{\mathrm{..}}{r}}_{\text{RP}}+{\stackrel{\mathrm{..}}{r}}_{{\text{G}}_{\text{others}}},$$(1)

where r is the Earth-centric radius vector of thedust particle and ${\stackrel{\mathrm{..}}{r}}_{{\text{G}}_{\text{E}}}$ is the acceleration caused by the gravity of the Earth, $${\stackrel{\mathrm{..}}{r}}_{{\text{G}}_{\text{E}}}=G{M}_{\text{E}}\nabla \left\{\frac{1}{r}\left[1-J_2{\left(\frac{R_{\text{E}}}{r}\right)}^2{P}_2(\cos \theta )\right]\right\}.$$(2)

Here, G is the gravitational constant, M_E is the mass of the Earth, R_E is the reference radius of Earth, θ is the colatitude in an Earth-centered body-fixed frame, and P₂ is the Legendre function of degree 2. In our simulation only the second-degree zonal harmonic, J₂ ≈ 1.082 × 10⁻³ (Pavlis et al. 2012), is considered. The variable ${\stackrel{\mathrm{..}}{r}}_{\text{RP}}$ denotes the acceleration due to the solar radiation pressure (Burns et al. 1979), $${\stackrel{\mathrm{..}}{r}}_{\text{RP}}=\frac{3Q_{\text{c}}{Q}_{\text{pr}}{\text{AU}}^2}{4{\left(r-r_{\text{S}}\right)}^2{\rho }_{\text{g}}{r}_{\text{g}}c}\left[1-\frac{\left(\stackrel{.}{r}-\stackrel{.}{r_{\text{S}}}\right)\cdot {\widehat{r}}_{\text{Sp}}}{c}\right]{\widehat{r}}_{Sp},$$(3)

where Q_c = 1.36 × 10³ W∕m² is the solar radiation energy flux at one astronomical unit (AU) distance and Q_pr is the solar radiation pressure efficiency factor, which depends on the size and material of the dust particle. The value of Q_pr was calculated from Mie theory for spherical particles using optical constants for silicates taken from Mukai (1989) (see our Fig. 1). The symbol ρ_g = 3500 kg m⁻³ denotes volumetric mass density for lunar dust (Solomon 1974), r_g is the radius of the particle, c is the speed of light, r_S is the vector from the Earth to the Sun, and $\widehat{r_{Sp}}$ is the unit vector from the Sun to the particle.

The acceleration by the gravity of Sun and Moon reads (e.g., Murray & Dermott 1999; Liu & Schmidt 2019) $${\stackrel{\mathrm{..}}{r}}_{{\text{G}}_{\text{others}}}=G{M}_{\text{S}}\left(\frac{r_{pS}}{r_{\text{pS}}^3}-\frac{r_S}{r_{\text{S}}^3}\right)+G{M}_{\text{M}}\left(\frac{r_{pM}}{r_{\text{pM}}^3}-\frac{r_M}{r_{\text{M}}^3}\right),$$(4)

where M_S and M_M are the masses of the Sun and the Moon. The symbol r_pS denotes the vector from the dust particleto the Sun, r_pM is the vector from the dust particle to the Moon, and r_M is the vectorfrom the Earth to the Moon.

In Sect. 4 we analyze the contribution of different perturbation forces to the evolution of eccentricity and inclination, d e∕d t and d i∕d t, which can be calculated from the perturbation equations as (Murray & Dermott 1999) $$\frac{\text{d}e}{\text{d}t}=\sqrt{a{\mu }^{-1}\left(1-e^2\right)}\left[\overline{R}\sin f+\overline{T}\left(\cos f+\cos E\right)\right]$$(5) $$\frac{\text{d}i}{\text{d}t}=\frac{r\overline{N}\cos \left(\omega +f\right)}{h}.$$(6)

Here, $\overline{R}$, $\overline{T}$, and $\overline{N}$ are the magnitudes of the radial, tangential, and normal components of the acceleration, respectively, which we determined from the equations of motion (Eq. (1)). The symbols a, e, i, ω, f, and E denote the osculating semimajor axis, eccentricity, inclination, argument of periapsis, true anomaly, and eccentric anomaly, respectively. We used the standard expressions from the two-body problem for radial distance, $$r=\frac{a\left(1-e^2\right)}{1+e\cos f},$$(7)

and angular momentum per unit mass, $$h=\sqrt{\mu a\left(1-e^2\right)}.$$(8)

The integration of a trajectory is terminated in our simulation if the grain hits the surface of the Earth or the Moon or if it escapes from the Earth-Moon system. We did not consider the effect of the Earth’s atmosphere. Collisions with the Earth or Moon could be easily missed due to the discrete time steps of the integrator. Therefore, cubic Hermite interpolation is used in our model to ensure the detection of collisions with the Earth or Moon if they occur between two consecutive time steps of the integrator (Liu et al. 2016).

Fig. 1 Solar radiation pressure efficiency, Qpr, for spherical silicate particles calculated from Mie theory.

3 Mass production rate and initial ejecta mass distribution

For the starting velocity we adopted a simple model in which the ejection direction is perpendicular to the lunar surface and all the particles have the same velocity distribution regardless of size. We assumed that the starting velocity follows a power law with exponent q. By normalizing ${\displaystyle {\int }_{v_0}^{v_{\text{max}}}p(v)}\text{d}v$, we obtain the form of the initial velocity distribution, $$p(v)=\frac{1-q}{v_{\text{max}}^{1-q}-v_0^{1-q}}{v}^{-q},$$(9)

where q = 2.2 was used, suitable for a regolith-covered surface (Krivov et al. 2003), v₀ = 130 m s⁻¹ is the minimum speed for lunar ejecta (Horányi et al. 2015), and v_max = 2v_esc ≈ 4.68 km s⁻¹ is adopted for the maximum ejection speed.

Via normalization to the total mass production rate, ${\displaystyle {\int }_{m_{\text{min}}}^{m_{\text{max}}}mp\left(m\right)\text{d}m}=M_{\text{total}}^{+}$, the distribution of mass can be expressed as $$p(m)=M_{\text{total}}^{+}\frac{1-\alpha }{m_{\text{max}}^{1-\alpha }-m_{\text{min}}^{1-\alpha }}{m}^{-\left(1+\alpha \right)},$$(10)

where α = 0.9 was inferred from the LDEX data (Horányi et al. 2015), m_min is the minimum ejecta mass, and m_max is the maximum ejecta mass. The ejecta size-distribution is obtained from the transformation $p\left(m\right)\text{d}m=p(r_{\text{g}})\text{d}{r}_{\text{g}}$ (e.g., Liu & Schmidt 2018), $$p(r_{\text{g}})=\frac{3}{r_{\text{max}}}\frac{M_{\text{total}}^{+}}{m_{\text{max}}}\frac{1-\alpha }{1-{\left(\frac{m_{\text{min}}}{m_{\text{max}}}\right)}^{1-\alpha }}{\left(\frac{r_{\text{g}}}{r_{\text{max}}}\right)}^{-1-3\alpha },$$(11)

where $r_{\text{max}}=\sqrt[3]{\frac{3m_{\text{max}}}{4\pi {\rho }_{\text{g}}}}$.

To apply these distributions of mass or radius, we must specify the total mass production rate, $M_{\text{total}}^{+}$. For the first step, the mass production rate per unit surface is $$\overline{M^{+}}=m_{\text{p}}{F}_{\text{p}}\cos \varphi Y_{\text{p}},$$(12)

where F_p is the number flux of impactors with characteristic mass m_p, treated here as a parallel beam that hits the surface with an angle ϕ from the surface normal, and cosϕ is the projection area factor. The yield, Y_p, is defined as the ratio of the ejecta particles’ mass to the impactors’ mass, which is also a function of the material of the target surface as well as the mass and velocity of the impactors (Koschny & Grün 2001a).

Considering the effects of oblique impacts from experiments by Gault (1973) and the dependence on projectile mass and velocity, we have $$Y_{\text{p}}\propto m_{\text{p}}^{0.23}{v}_{\text{p}}^{2.46}{\cos }^2\varphi .$$(13)

The mass production on the surface of the Moon varies with local time because the angle ϕ varies for projectile populations that approach the Earth-Moon system from different directions. Six meteoroid populations, helion (HE), antihelion (AH), apex (AP), antiapex (AA), northern toroidal (NT), and southern toroidal (ST), are found to be plausible sources of high-speed impacts (Szalay et al. 2019). Following Szalay & Horányi (2015), and taking Eqs. (13) and (14) into account, the mass production rate per unit surface reads $$\overline{M^{+}}=C{\sum }_s\underset{w_s}{\underbrace{F_s{m}_s^{1.23}{v}_s^{2.46}}}{\cos }^3{\varphi }_s{\Theta }_{\text{H}}\left(\cos {\varphi }_s\right).$$(14)

Here, C is a normalization constant, s labels the different projectile populations, F_s, m_s, and v_s are the number flux, characteristic mass, and impact velocity for each source, and Θ_H denotes the Heaviside function, which accounts for the fact that a given projectile population can only reach one hemisphere of the Moon. We followed Szalay et al. (2019) and used relative contributions w_s = (0.198, 0.198, 0.303, 0.025, 0.138, 0.138) for the HE, AH, AP, AA, NT, and ST sources of meteoroids. The symbol ϕ_s is the impact angle for each source from the surface normal, $$\cos {\varphi }_s=\sin \theta \sin {\theta }_s\cos (\phi -{\phi }_s)+\cos \theta \cos {\theta }_s.$$(15)

Here, θ is the colatitude of the surface element, θ_s is the colatitude for sources, θ_NT = 30°, θ_ST = 150° with φ_NT = φ_ST = 0° (Szalay et al. 2019), φ is defined as the longitude from the apex of the motion of the Earth-Moon system for the surface element, φ_s is the source-specific longitude of approach in the ecliptic plane (HE, AH, AP, and AA), φ_HE = 65°, φ_AH = 295°, φ_AP = 0°, and φ_AA = 180° (Szalay & Horányi 2016). Therefore, Eq. (14) can be written as a function of the longitude, φ, and colatitude, θ, on the Moon. The total mass production rate, $M_{\text{total}}^{+}$, for the Moon was estimated to be approximately 0.2 kg s⁻¹ (Szalay et al. 2019; Pokornỳ et al. 2019). By normalization, ${\displaystyle {\int }_0^{\pi }{R}_{\text{M}}^2\sin \theta \text{d}\theta }{\displaystyle {\int }_0^{2\pi }\text{d}\phi \overline{M^{+}}(\phi ,\theta )}=M_{\text{total}}^{+}$, the normalization constant, C, in Eq. (14) can be fixed.

In our simulation, 200 starting positions on the surface of the Moon, including 20 longitudes and 10 colatitudes, are selected, which divides the whole surface of the Moon into 20 × 10 elements. Longitudes, φ, are equidistant from 0 to 2π, and colatitudes, θ, vary from 0 to π. Accordingly, the mass production rate for a given surface element (i, j) reads $$M^{+}(i,j)={\displaystyle {\int }_{{\theta }_j^1}^{{\theta }_j^2}{R}_{\text{M}}^2\sin \theta }\text{d}\theta {\displaystyle {\int }_{{\phi }_i^1}^{{\phi }_i^2}\text{d}\phi }\overline{M^{+}}(\phi ,\theta ),$$(16)

where i is the longitude index from 1 to 20, j is the colatitude index from 1 to 10, R_M is the radius of the Moon, ${\phi }_i^1$ and ${\phi }_i^2$ are boundary longitudes of surface element (i, j), and ${\theta }_j^1$ and ${\theta }_j^2$ are boundary colatitudes of surface element (i, j). The distribution of mass production rates is shown in Fig. 2.

Correspondingly, the size distribution for surface element (i, j) reads $$p(r_{\text{g}},i,j)=\frac{3}{r_{\text{max}}}\frac{M_{\text{+}}(i,j)}{m_{\text{max}}}\frac{1-\alpha }{1-{\left(\frac{m_{\text{min}}}{m_{\text{max}}}\right)}^{1-\alpha }}{\left(\frac{r_{\text{g}}}{r_{\text{max}}}\right)}^{-1-3\alpha }.$$(17)

4 Simulation results

To cover the size distribution of ejecta, we selected nine dust particle radii: 0.2, 0.5, 1, 2, 5, 10, 20, 50, and 100 μm. Here, r_min = 0.2 μm and r_max = 100 μm, so $m_{\text{min}}=\frac{4}{3}{\rho }_{\text{g}}\pi r_{\text{min}}^3$ and$m_{\text{max}}=\frac{4}{3}{\rho }_{\text{g}}\pi r_{\text{max}}^3$. For the initial velocity, ten velocities in the range [0.95 v_esc, 2 v_esc] were used, where v_esc is the escape velocity of the Moon. We note that the velocity distribution is normalized to unity in the range [v₀, v_max = 2v_esc], where v₀ = 130 m s⁻¹ is the minimal starting velocity of ejecta. Integrating over the velocity distribution in the interval [0.95v_esc, 2v_esc] then removes in our normalization the grains in [v₀, 0.95v_esc], which all fall back to the surface. A small fraction of grains in the range [0.95v_esc, v_esc] escape the Moon due to three-body effects. They are considered in our evaluation of the simulations.

To average over seasonal effects, we assumed that particles start at 12 different times, covering equidistantly one period of the ascending node precession of the Moon. We verified that this choice also covers the different phases of the Moon fairly uniformly. The number of grains necessary to cover the ranges of radius, velocity, time, colatitude, and longitude are 9, 10, 12, 10, and 20, respectively. Thus, the total number of particles is 216 000, which requires a huge amount of CPU time. The long-term simulations were carried out on the supercomputer located at the Finnish CSC-IT Center for Science.

Fig. 2 Mass production rate per surface element on the lunar surface. Zero longitude is defined by the direction of the apex of the motion of the Earth-Moon system.

4.1 Particle sinks and particle lifetimes

The final state of particles ejected from the surface of the Moon after 100 yr is shown in Table 2. Most dust particles (>88%) ultimately leave the Earth-Moon system after 100 Earth years,regardless of size.

Only a very small fraction of grains (>2 μm) remain in orbit after 100 yr. We do not integrate these trajectories further, to save CPU time. The fraction of grains that rapidly re-impact the Moon within one orbital period is approximately 22%, which only weakly depends on grain size. We note that this fraction depends on our choice of v_min = 0.95v_esc, so the fraction of grains falling back to the Moon within one orbital period is not listed in Table 2. The grains that re-impact the Moon after a longer time are listed in Table 2 in the column “Hit Moon”. Small particles are more sensitive to solar radiation pressure. The fraction of grains that hit the Moon increases with grain size. The fraction that hit the Earth grows from 1.2 to 10% at a grain size of 2 μm and then begins to drop to a value of 3% for 100 μm grains.

Figure 3a shows the average life spans as a function of the size of the dust particles. Small particles, including those sized 0.2, 0.5, and 1 μm, have short lifetimes, less than 0.1 yr. For particles between 1 and 10 μm, the life span increases rapidly with grain size. Large particles with sizes bigger than 10 μm have longer lifetimes, approximately 0.7 Earth years. The parameter β denotes the ratio of solar radiation pressure to gravity (Burns et al. 1979): $$\beta =\frac{3Q_{\text{c}}{Q}_{\text{pr}}{\text{AU}}^2}{4G{M}_{\text{s}}{\rho }_{\text{g}}{r}_{\text{g}}c}.$$(18)

As shown in Fig. 3b, the effect of solar radiation pressure is large for small particles, which induces high eccentricities, such that small grains are rapidly expelled from the Earth-Moon system and the life spans of those particles are short. With increasing grain size, β decreases rapidly and becomes negligible for grains of 10 μm or larger. Thus, the lifetimes of these particles are much longer.

Fig. 3 Average lifetimes and β for grains of different sizes. Panel a: average lifetimes. Panel b: β.

4.2 Steady-state configuration

A cylindrical grid $\left({\rho }_{\text{c}},{\varphi }_{\text{c}},z_{\text{c}}\right)$ is utilized in the simulation process, where ${\rho }_{\text{c}}=\sqrt{x^2+y^2}$, ϕ_c = atan2(y, x) and z_c = z. Indexes i_cell, j_cell, and k_cell are used to label grid cells. Hence, the number density of a given grid cell (i_cell, j_cell, k_cell) reads $$n(i_{\text{cell}},j_{\text{cell}},k_{\text{cell}})={\displaystyle \sum _i{\displaystyle \sum _j{\displaystyle {\int }_{r_{\text{min}}}^{r_{\text{max}}}p}}}(r_{\text{g}},i,j)\text{d}{r}_{\text{g}}{\displaystyle {\int }_{v_{\text{min}}}^{v_{\text{max}}}p}(v)\text{d}v\text{nn }\frac{n(i_{\text{cell}},j_{\text{cell}},k_{\text{cell}},r_{\text{g}},v,i,j)\Delta t}{V(i_{\text{cell}},j_{\text{cell}},k_{\text{cell}})N_{\text{start}}(r_{\text{g}},v,i,j)}.$$(19)

For all trajectories, the grain positions are stored at equal time intervals, Δt. The symbol n(i_cell, j_cell, k_cell, r_g, v, i, j) is the plain number of dust particles with radius r_g, starting velocity v, and starting surface element (i, j) recorded in the cell (i_cell, j_cell, k_cell), V (i_cell, j_cell, k_cell) is the volume of the cell (i_cell, j_cell, k_cell), and N_start(r_g, v, i, j) is the number of grains with a given starting radius r_g, velocity v, and starting element (i, j).

In a similar manner we can obtain the geometric optical depth covered by the grains in a cell, by integrating over $\pi r_g^2$ in Eq. (21) and dividing by, instead of volume, the surface of the cell (along a given line of sight). The normal optical depth and number density distribution are illustrated in Fig. 4. Figure 4a shows the normal optical depth of the dust in a Sun-oriented frame, where the positive x axis always points in the direction of the Sun. Particles occupy a torus between the Earth and the Moon with an outer edge that is separated by roughly one lunar Hill radius from the orbit of the Moon. The torus is asymmetric and slightly offset toward the Sun. From Fig. 4b, we can see that dust particles are distributed nearly symmetrically about the equator of the Earth, covering a wide vertical range of up to tens of Earth radii. We note, however, that the configuration shown in Fig. 4b was evaluated in an inertial frame, averaging over grains ejected from the Moon at 12 times equidistantly over one period of the ascending node precession of the Moon, and over an evolution time of more than 100 yr. Thus, the effects of the inclination of the Sun and the Moon are averaged out. The instantaneous dust configuration is expected to differ somewhat from the one shown in Fig. 4b because the instantaneous Laplace plane of the dynamical problem does not coincidewith the equatorial plane of the Earth.

The differential size distribution in the torus between the Earth and the Moon is shown in Fig. 5. The size distribution in the torus after long-term simulations remains close to the initial ejecta distribution but has flattened in the size range of 1–10 μm. This can be attributed to the rapidly changing life span for grains in this size range (see Fig. 3a).

Fig. 4 Optical depth and number density of particles ejected from the surface of the Moon. Panel a: normal geometric optical depth. Panel b: azimuth-averaged number density in the ρc − zc plane. The red circle is the Earth, and the red line denotes the orbit of the Moon.

Fig. 5 Steady-state differential size distribution for the torus from simulations. The red line denotes the initial ejecta size distribution from $p(r_{\text{g}})\propto r_{\text{g}}^{-q}$ normalized to unity, and the blue line denotes the steady-state size distribution in the torus from Fig. 4.

4.3 Analysis of orbital elements

Figure 6 shows the distribution of eccentricity for 0.2 μm particles. The eccentricities are dispersed from 0 to 15, that is, a large number (~80%) of particles are hyperbolic, which explains their short life spans. These grains leave the Earth-Moon system in a very shorttime and contribute little to the number density. The evolution of eccentricity for a typical 0.2 μm particle with a lifetime of 0.12 yr is shown in Fig. 7a. The initial eccentricity is near 0.5, that is, this particle starts from an elliptical orbit. Initially, the eccentricity fluctuates slightly within a few days and then grows continually, representing a typical evolution from an ellipse to a hyperbola, which is a consequenceof the strong effect of solar radiation pressure. Figure 7b depicts the contribution of the perturbations induced by lunar gravity, solar gravity, and radiation pressure on the growth rate d e∕dt, calculated from Eq. (5). Solar gravity has a negligible influence on de∕dt. Initially, thede∕dt due to solar radiation pressure is a small negative or positive value. After 0.08 yr, solar radiation pressure dominates and d e∕dt remains positive. As a result, the eccentricity increases monotonically, and finally, the orbit becomes hyperbolic.

Large particles remain bound to the Earth-Moon system for a much longer time, and they make a significant contributionto the dust population. The distributions of the semimajor axis, eccentricity, and inclination of all 100 μm particles are shown in Fig. 8. The semimajor axis peaks around 20 R_E, mainly ranging from 10 R_E to 50 R_E. The eccentricity has a centered distribution below unity, peaking near 0.7. The inclination exhibits two peaks, at values of 20° and 140°, that is, a significant fraction of grains are on retrograde orbits.

We also calculated the distributions of the solar angle for 1 and 100 μm grains (Fig. 9). The solar angle is defined as the angle between the grain’s orbital pericenter and the Sun as seen from Earth (Hamilton 1993): $${\varphi }_{\odot }=\Omega +\omega -{\psi }_{\odot }.$$(20)

Here, Ω is the longitude of the ascending node, ω is the argument of periapsis, and ψ_⊙ is the solar longitude in the inertial frame.

From Fig. 9, the solar angle distribution is approximately symmetric about ϕ_⊙ = 180°. The solar angle of 1 μm particles peaks at 180°, that is, the population of these grains is offset toward the Sun. For 100 μm particles, the curve is smoother but shows two peaks, at 0° and 180°, which indicates that the offset of the torus formed by 100 μm particles is much smaller compared to that of the 1 μm grains. For particles of different sizes, the distribution exhibits different kinds of asymmetry, with offsets toward or away from the Sun. However, averaging over grain size, the torus is slightly offset toward the Sun (see Fig. 4a).

Figures 10 and 11 show the evolution of radial distance and inclination for a single 100 μm particle. This is one of the rare particles in our simulations that still remained in orbit after 100 yr. The grain covers a range between 10 R_E and 50 R_E, which is consistent with the distribution in Fig. 4a. The inclination oscillates from 0° to 60°. This particle is always moving on a prograde orbit within its life span. Such particles make a large contribution to the distribution of inclination (see Fig. 8c).

A considerable part of the particles in our simulations are on a retrograde orbit. The fractions of retrograde particles for different sizes are shown in Fig. 12. Some grains are ejected from the Moon directly into retrograde orbits (Colombo et al. 1966), but a larger number of grains evolve from a prograde to a retrograde orbit in the long-term simulation. To show this, a 1 μm particle is used here as an example. The evolution of inclination of this particle and contributions of specific forces to d i∕dt calculated from Eq. (6) are shown in Figs. 13 and 14. From Fig. 13, we see that the initial inclination of the grain is near 60°, and the final one is 170°. This particle evolves from a prograde orbit to a retrograde orbit multiple times. As can be seen in Fig. 14, the effects of gravity of the Moon and Sun on the evolution of d i∕dt are much smaller than the one induced by the solar radiation pressure, although there are some peaks from encounters with the Moon. The approximate change in inclination can be evaluated via the integration of the contribution to d i∕dt due to solar radiation pressure, which is about 100°.

Fig. 6 Distribution of eccentricity for 0.2 μm particles.

Fig. 7 Evolution of eccentricity and de∕dt for one 0.2 μm particle. Panel a: evolution of eccentricity. Panel b: evolution of de∕dt.

Fig. 8 Distribution of orbital elements for 100 μm particles ejected from the surface of the Moon. Panel a: distribution of the semimajor axis. Panel b: distribution of eccentricity. Panel c: distribution of inclination.

Fig. 9 Distribution of solar angle for 1 and 100 μm particles ejected from the surface of the Moon.

Fig. 10 Evolution of r for a 100 μm particle.

Fig. 11 Evolution of inclination for a 100 μm particle.

Fig. 12 Fraction of retrograde particles.

Fig. 13 Evolution of inclination for a 1 μm particle.

5 Conclusions

In this work we have derived the steady-state configuration of particles in the Earth-Moon system that have been ejected from the lunar surface in hypervelocity impacts of micrometeoroids (Horányi et al. 2015). A variety of forces, including the solar radiation pressure and the gravity of the Sun, the Earth, and the Moon, are considered in a numerical exploration of the trajectories of the ejected particles in the system. Expanding on previous work in the literature (Szalay & Horányi 2015, 2016; Szalay et al. 2019), we calculated the initial ejecta distribution and mass production rate on the lunar surface, which we used as initial conditions for the integrations. The final states (sinks), the average life spans, and the fraction of retrograde grains have been derived as functions of grain size. We also show examples for the evolution and distribution of orbital elements of particles. Small particles tend to rapidly evolve into hyperbolic orbits (on timescales of weeks), while most large particles remain in elliptical orbits for a longer time (up to a year). Manyof the grains develop high inclinations, and we find that a substantial fraction evolve into retrograde orbits.

From our long-term integrations we find that about 3.6 × 10⁻³ kg s⁻¹ (1.8%) particles escape from the lunar gravity. These particles form a tenuous broad torus, the densest part of which roughly spans a distance of 40 R_E from the Earth. We also estimated the mass rate of particles coming into the Earth’s atmosphere, $M_{\text{E}}^{+}=2.3\times {10}^{-4}\text{kg}{\text{s}}^{-1}$, which is a non-negligible value. As the space activities in the Earth-Moon system are much more frequent compared with other planets, dust around the Earth seriously threatens not only the mechanical structure of spacecraft but potentially also the health of astronauts. In Apollo missions, astronauts were exposed to the lunar dust environment, and respiratory, dermal, and ocular irritations from lunar dust were reported (Turci et al. 2015). NASA established the Lunar Airborne Dust Toxicity Advisory Group (LADTAG) to assess the risk degree for spacecraft and astronauts in lunar missions (Khan-Mayberry 2008). Thus, analysis of the characters for particles in the torus around the Earth is of great significance for assessing the space environment and ensuring thesecurity of explorations.

We characterize the torus in terms of the normal geometric optical depth and number density. Owing to the effect of solar radiation pressure, the torus is mildly offset toward the solar direction. The peak normal optical depth we obtain is on the order of 10⁻¹², which is small compared to the optical depth of the Thebe extension of Jupiter’s gossamer rings, ~ 10⁻⁹ (De Pater et al. 2018), and the limits inferred for the yet undetected Phobos ring of Mars, ~ 10⁻⁸ (Showalter et al. 2006). The peak number density we find for the lunar torus is on the order of 10⁻⁹ m⁻³, which is about four orders of magnitude smaller than the density of the Deimos ring inferred for radii of 15 μm by Juhász & Horányi (1995) and four and five orders of magnitude smaller than the number for the Deimos and Phobos rings (≥ 0.5 μm) determined by Liu & Schmidt (2021), respectively. For an in situ detector on a circular orbit at a distance of 10 R_E pointing in the prograde orbit direction, this translates into a flux of about 50 grains per m² and per year (taking into account that a fraction of the grains are on retrograde orbits). From the Interplanetary Meteoroid Engineering Model (IMEM; Dikarev et al. 2005). we find that this is a factor of several smaller than the flux of micron-sized interplanetary micrometeoroids at 1 AU. The directional distributions of the lunar grains and the interplanetary particles, however, would be different.

Our results on the grain lifetimes, their sinks, and their typical orbital evolution are robust because the relevant perturbation forces are taken into account and the respective parameters (masses of the gravitational perturbers, the solar constant, and the solar radiation pressure efficiency factor) are known up to small uncertainties. Large uncertainties, however, remain in our calibration of number density and optical depth. On one hand, our simulation results are based on the mass production rate of $M_{\text{total}}^{+}=0.2\text{kg}{\text{s}}^{-1}$ reported by Szalay et al. (2019). If we follow the method used in the literature for an analysis of ejecta from the Galilean moons (Krivov et al. 2003; Sremčević et al. 2003), employ the mass flux at 1AU from the IMEM model (Dikarev et al. 2005), and use the empirical expression for the yield suggested by Koschny & Grün (2001a), we obtain a mass production rate that is about one order of magnitude higher. Although in our analysis we prefer to use the value from the direct measurement, this suggests that the uncertainty in the calibration of our model is at least one order of magnitude. On the other hand, we assume a power law distribution for the ejection velocity, whose high velocity tail determines the fraction of ejecta that escape the Moon in our model. In the literature, a different functional form was fitted successfully to the bound particles detected by LDEX (Szalay & Horányi 2016). If the tail of the velocity distribution of the lunar ejecta turns out to be different from the one assumed in our paper, then the dust densities derived from our model will also change.

Fig. 14 Variation in di∕dt due to various forces for a 1 μm particle. Panel a: variation in di∕dt due to lunar gravity. Panel b: variation in di∕dt due to solar gravity. Panel c: variation in di∕dt due to solar radiation pressure.

Acknowledgements

This work was supported by National Natural Science Foundation of China (No. 12002397) and by the China Scholarship Council (CSC, 201906220134). We acknowledge CSC – IT Center for Science for the allocation of computational resources.

References

All Tables

All Figures

	Fig. 1 Solar radiation pressure efficiency, Qpr, for spherical silicate particles calculated from Mie theory.
In the text

	Fig. 2 Mass production rate per surface element on the lunar surface. Zero longitude is defined by the direction of the apex of the motion of the Earth-Moon system.
In the text

	Fig. 3 Average lifetimes and β for grains of different sizes. Panel a: average lifetimes. Panel b: β.
In the text

	Fig. 4 Optical depth and number density of particles ejected from the surface of the Moon. Panel a: normal geometric optical depth. Panel b: azimuth-averaged number density in the ρc − zc plane. The red circle is the Earth, and the red line denotes the orbit of the Moon.
In the text

	Fig. 5 Steady-state differential size distribution for the torus from simulations. The red line denotes the initial ejecta size distribution from $p(r_{\text{g}})\propto r_{\text{g}}^{-q}$ normalized to unity, and the blue line denotes the steady-state size distribution in the torus from Fig. 4.
In the text

	Fig. 6 Distribution of eccentricity for 0.2 μm particles.
In the text

	Fig. 7 Evolution of eccentricity and de∕dt for one 0.2 μm particle. Panel a: evolution of eccentricity. Panel b: evolution of de∕dt.
In the text

	Fig. 8 Distribution of orbital elements for 100 μm particles ejected from the surface of the Moon. Panel a: distribution of the semimajor axis. Panel b: distribution of eccentricity. Panel c: distribution of inclination.
In the text

	Fig. 9 Distribution of solar angle for 1 and 100 μm particles ejected from the surface of the Moon.
In the text

	Fig. 10 Evolution of r for a 100 μm particle.
In the text

	Fig. 11 Evolution of inclination for a 100 μm particle.
In the text

	Fig. 12 Fraction of retrograde particles.
In the text

	Fig. 13 Evolution of inclination for a 1 μm particle.
In the text

	Fig. 14 Variation in di∕dt due to various forces for a 1 μm particle. Panel a: variation in di∕dt due to lunar gravity. Panel b: variation in di∕dt due to solar gravity. Panel c: variation in di∕dt due to solar radiation pressure.
In the text