© ESO, 2016

1. Introduction

The nonuniformity of the cratering distribution on a satellite or a planet when it comes under bombardment is called cratering asymmetry. There are three types. The first is the leading/trailing asymmetry, which is due to the synchronous rotation of a satellite. When synchronously locked, a satellite’s velocity always points to its leading side, so that this hemisphere tends to gain a higher impact probability, higher impact speed and more normal impacts. Because the crater diameter depends on impact speed and incidence angle (e.g., Le Feuvre & Wieczorek 2011), this can also lead to an enhanced crater size on the leading side. This is the so-called apex/antapex effect (Zahnle et al. 2001; Le Feuvre & Wieczorek 2011). Zahnle et al. (1998, 2001) and Levison et al. (2000) investigated this effect on the giant planets and their satellites in detail. This work focuses on the Moon. Its leading/trailing asymmetry has been proposed by theoretical works (Zahnle et al. 1998; Le Feuvre & Wieczorek 2011), confirmed by numerical simulations (Gallant et al. 2009; Ito & Malhotra 2010), and directly verified by seismic observations (Kawamura et al. 2011; Oberst et al. 2012) and observations of rayed craters (Morota & Furumoto 2003). The second type, the pole/equator asymmetry (or latitudinal asymmetry), which is an enhancement of low-latitude impacts compared to high-latitude regions resulting from the concentration of low-inclination projectiles, has often been reported as well (Le Feuvre & Wieczorek 2008, 2011; Gallant et al. 2009; Ito & Malhotra 2010). The third type is the near/far asymmetry, which is most uncertain. It has been proposed that Earth focuses the projectiles onto the near side of the Moon as a gravitational lens and thus the craters there are enhanced (Wiesel 1971), or that Earth is an obstacle in the projectiles’ trajectories so that the near side is shielded from impacts (Bandermann & Singer 1973). Contradicting and limited conclusions have not led to consensus.

The cratering asymmetry depends on not only the lunar orbit, but also on the impactor population. After analyzing the size distributions of craters, Strom et al. (2005, 2015) suggested that there are two impactor populations in the inner solar system: the main-belt asteroids (MBAs), which dominated during the late heavy bombardment (LHB) ~ 3.9 Gya, and the near-Earth objects (NEOs), which have dominated since about 3.8−3.7 Gya. The two populations are different in their orbital and size distributions, fluxes, and origins, which means that there is no reason to expect the cratering asymmetries generated by them to be the same.

We note the MBAs referred to here are the asteroids that occupied the region where the current main belt is during the LHB. Although their semi-major axes a_p = 2.0−3.5 AU were the same as today, their eccentricities and inclinations were greatly excited by migrating giant planets (Gomes et al. 2005). Instead, their size distribution has changed little after the first ~100 Myr (Bottke et al. 2005). Several surveys of the current main belt have reported power-law breaks of its size distribution. The Sloan Digital Sky Survey (SDSS; Ivezić et al. 2001) distinguished two types of asteroids with different albedos (0.04 and 0.16 for blue and red asteroids) and found that for the whole sample throughout the whole belt, the cumulative size distribution has a slope α_p = 1.3 for diameters d_p = 0.4−5 km and 3 for 5−40 km. Parker et al. (2008) also confirmed this conclusion based on the SDSS moving object catalog 4, including ~88 000 objects. The first (Yoshida et al. 2003) and second (Yoshida & Nakamura 2007) runs of the Sub-km Main Belt Asteroid Survey (SMBAS), which each found 861 MBAs in R-band and 1001 MBAs in both B and R-bands, reported α_p = 1.19 ± 0.02 for d_p = 0.5−1 km and α_p = 1.29 ± 0.02 for d_p = 0.6−1 km, claimed to be consistent with SDSS for MBAs smaller than 5 km. The size distribution of craters formed by MBAs also has a complex shape, whose cumulative slope is α_c = 1.2 for crater diameters d_c ≲ 50 km, 2 for 50−100 km, and 3 for 100−300 km (Strom et al. 2015).

The NEO population is relatively well understood and suggested to have been in steady state for the past ~3 Gyr, constantly resupplied mainly from the main belt (Bottke et al. 2002). Bottke et al. (2002) derived the debiased orbital and size distributions of NEOs by fitting a model population to the known NEOs: the orbits are constrained in a range a_p = 0.5−2.8 AU, e_p< 0.8, i_p< 35°, the perihelion distance q_p< 1.3 AU, and the aphelion distance Q_p> 0.983 AU; the size distribution is characterized by a single slope α_p = 1.75 ± 0.1 for d_p = 0.2−4 km. The corresponding slope for the crater size distribution indicated by Strom et al. (2015) is α_c = 2 for d_c = 0.02−100 km. A few recent works have investigated the cratering asymmetry generated by NEOs based on the debiased NEO model described in Bottke et al. (2002). Gallant et al. (2009) ran N-body simulations and determined an apex/antapex ratio (ratio of the crater density at the apex to antapex) of 1.28 ± 0.01, a polar deficiency of ~10%, and the absence of a near/far asymmetry. In a similar work but with a different numerical model, Ito & Malhotra (2010) found a leading/trailing hemispherical ratio of 1.32 ± 0.01 and also a polar deficiency of ~10%. Le Feuvre & Wieczorek (2008) analytically predicted the lunar pole/equator ratio to be 0.90. Le Feuvre & Wieczorek (2011) followed and semi-analytically estimated that involving both longitudinal and latitudinal asymmetries, the lunar crater density was minimized at (90° E, ±65° N) and maximized at the apex with deviations of about 25% with respect to the average for the current Earth–Moon distance, resulting in an apex/antapex ratio of 1.37 and a pole/equator ratio of 0.80. From observations, Morota & Furumoto (2003) reported an apex/anntapex ratio of ~1.5 for 222 rayed craters with d_c> 5 km on lunar highlands, formed by NEOs in the past ~1.1 Gyr.

However, according to Strom et al. (2015), for craters with diameters larger than 10 km, those that the MBAs formed exceed the other population by more than an order of magnitude. This underlines the need for taking the craters and the cratering asymmetry that the MBAs contribute to into account, and thus requires the awareness of the lunar orbit during the LHB, namely, the dominance of MBAs. The Earth–Moon system evolved drastically in the early history. As a general trend, it is accepted that the Moon has been receding from Earth because of tidal dissipation, but a wide spectrum of opinions exists for the details of this picture. After the timescale problem was raised (Gerstenkorn 1955; MacDonald 1964; Goldreich 1966; Lambeck 1977; Touma & Wisdom 1994) and then solved through ocean models (Hansen 1982; Webb 1982; Ross & Schubert 1989; Kagan & Maslova 1994; Kagan 1997; Bills & Ray 1999), it is known today that the tidal dissipation factor of Earth must have been much larger in the distant past, or in other words, the tidal friction was much weaker than today (e.g., Bills & Ray 1999). Because tidal friction depends on not only the lunar orbit but also the terrestrial ocean shape, which is associated with continental drift, the exact history of Earth’s tidal dissipation factor and hence the early evolution of the lunar orbit cannot be confirmed.

Still, we list some efforts of ocean modelers here. Hansen (1982) was the first to introduce Laplace’s tidal equations in calculations of oceanic tidal torque. He modeled four cases, combinations of two idealized continentalities and two frictional resistance coefficients, and found the Earth–Moon distances at 4.5 Gya ranging from 38 to 53 R_⊕, leading to nearly the same values at 3.9 Gya when the LHB occurred. Webb (1980, 1982) developed an average ocean model and obtained the dates of the Gerstenkorn event as 3.9 and 5.3 Gya with and without solid Earth dissipation included, respectively. According to these two results, the Earth–Moon distance at 3.9 Gya should be no larger than 25 R_⊕ or about 42 R_⊕. The author claimed, however, that the results should only be taken qualitatively. Ross & Schubert (1989) simulated a coupled thermal-dynamical evolution of the Earth–Moon system based on equilibrium ocean model, resulting in a_M being 47 R_⊕ and e_M being 0.04 at 3.9 Gya when the evolution timescale and parameters such as final a_M and e_M were required to be realistic. Kagan & Maslova (1994) described a stochastic model that considered fluctuating effects of the continental drift. Their two-mode resonance approximations gave rise to evolutions of tidal energy dissipation that were consistent with global paleotide models. By adopting a reproduced tidal evolution with a timescale as close as possible to the realistic one, the Earth–Moon distance is estimated to be about 47 R_⊕ at 3.9 Gya. This timescale can vary by billions of years, of course, depending on the values of the resonance lifetime.

The uncertain lunar orbital status during the LHB means that a reliable relationship between the lunar orbit and the cratering asymmetry is needed. If the influence of the former on the latter were significant and an incorrect condition of the Earth–Moon system were assumed, the estimated cratering asymmetry would be also incorrect. Conversely, this also implies that the early history of the Earth–Moon system could be inferred from the observed crater record if the portion of craters that formed during the LHB were selected. The influence of the lunar orbit and of the impactor population on the cratering asymmetry have not been sufficiently studied so far. Numerical simulations can only offer an empirical estimation based on a limited number of cases: Zahnle et al. (2001), who simulated impacts of ecliptic comets on giant planet satellites, combined the results of case α_p = 2.0 and case α_p = 2.5 with the satellite orbital speed v_orb = 10.9 km s^-1 and the impactor speed at infinity in the rest frame of the planet fixed at v_∞ ≈ 5 km s^-1, and thus derived a semi-empirical description of the crater density ${\mathit{N}}_{\mathrm{c}}\mathrm{\propto }{\left(\mathrm{1}\mathrm{+}\frac{{\mathit{v}}_{\mathrm{orb}}}{\sqrt{\mathrm{2}{\mathit{v}}_{\mathrm{orb}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathrm{\infty }}^{\mathrm{2}}}}\cos \mathit{\beta }\right)}^{\mathrm{2.0}\mathrm{+}\mathrm{(}\mathrm{1.4}\mathit{/}\mathrm{3}\mathrm{)}{\mathit{\alpha }}_{\mathrm{p}}}\mathit{,}$(1)where β is the angular distance from the apex; using the debiased NEO model, Gallant et al. (2009) ran simulations with Earth–Moon distances of a_M = 50, 38, 30, 20, and 10 R_⊕, and confirmed the negative correlation between a_M and leading/trailing asymmetry degree. A semi-analytical method capable of producing cratering of the Moon caused by current asteroids and comets in the inner solar system was proposed by Le Feuvre & Wieczorek (2011), who suggested a fit relation ${\mathit{N}}_{\mathrm{c}}\mathrm{\left(}\mathrm{apex}\mathrm{\right)}\mathit{/}{\mathit{N}}_{\mathrm{c}}\mathrm{\left(}\mathrm{antapex}\mathrm{\right)}\mathrm{=}\mathrm{1.12}{\mathrm{e}}^{\mathrm{-}\mathrm{0.0529}\mathrm{(}{\mathit{a}}_{\mathrm{M}}\mathit{/}{\mathit{R}}_{\mathrm{\oplus }}\mathrm{)}}\mathrm{+}\mathrm{1.32}\mathit{,}$(2)which is valid between 20 and 60 R_⊕.

In this paper, we first derive the formulated lunar cratering distribution through rigorous integration based on a planar model in Sect. 2, where the leading/trailing asymmetry is unambiguously confirmed. Then we show the numerical simulations of the Moon under bombardment of MBAs during the LHB over various Earth–Moon distances in Sect. 3. The simulation results are compared to our analytical predictions, and the cratering distribution of coupled asymmetries is described. Last, Sect. 4 compares NEOs and MBAs and shows the consistence of our analytical model with related works.

2. Analytical deduction

This section shows how we analytically deduced formulation series describing the spatial distributions of impact speed, incidence angle, normal speed, crater diameter, impact density, and crater density on the lunar surface under bombardment, based on a planar model. These formulations directly prove the existence of the leading/trailing cratering asymmetry and the identity of the near and far sides. The amplitude of the leading/trailing asymmetry is proportional to the ratio of the impactor encounter speed to the lunar orbital speed, implying that it might be possible to derive the lunar orbit or the impactor population from crater record.

2.1. Assumptions and precondition

Our model includes the Sun, Earth, the Moon, and the impactors (asteroids and/or comets) in the inner solar system, which are assumed to be located farther away from the Sun than Earth. The main assumptions we adopt are listed below.

• 1.

  The orbits of Earth, the Moon, and the impactors as well as thelunar equator are coplanar.

• 2.

  The orbits of Earth and the Moon are circular.

• 3.

  The gravitation on the impactors due to Earth and the Moon are ignored.

• 4.

  Earth is treated as a particle.

The first assumption does not allow describing latitudinal cratering variation, which is investigated in our numerical simulations. It helps avoiding the coupling of two cratering asymmetries and concluding a pure leading/trailing asymmetry formulation. The Moon’s geometrical libration is neglected by the first two assumptions. The third partly means that the acceleration of the impactor velocity is not accounted for, which is valid on condition that the Earth–Moon distance is larger than ~17R_⊕, according to Le Feuvre & Wieczorek (2011). It also means the gravitational lensing by Earth is excluded, and moreover, the last assumption excludes the Earth’s shielding for the Moon. This causes the symmetry between the near and far sides.

A precondition for the cratering asymmetry is that the Moon must be synchronously rotating, because otherwise any longitudinal variation would vanish as the lunar hemisphere facing the Earth changes constantly. According to Zhou & Lin (2008), the timescale for the Moon to reach the 1:1 spin-orbit resonance can be estimated by ${\mathit{\tau }}_{\mathrm{syn}}\mathrm{~}{\left(\frac{{\mathit{R}}_{\mathrm{M}}}{{\mathit{a}}_{\mathrm{M}}}\right)}^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{tide}}\mathrm{=}\frac{\mathrm{4}{\mathit{Q}}_{\mathrm{M}}^{\mathrm{\prime }}}{\mathrm{63}{\mathit{n}}_{\mathrm{M}}}\frac{{\mathit{m}}_{\mathrm{M}}}{{\mathit{m}}_{\mathrm{\oplus }}}{\left(\frac{{\mathit{a}}_{\mathrm{M}}}{{\mathit{R}}_{\mathrm{M}}}\right)}^{\mathrm{3}}\mathit{,}$(3)where τ_tide is the tidal circularization timescale, m_⊕ is the Earth mass, n_M, a_M, R_M, m_M, and ${\mathit{Q}}_{\mathrm{M}}^{\mathrm{\prime }}$ are the mean motion, semi-major axis, radius, mass, and effective tidal dissipation factor of the Moon. A measure of this timescale is τ_syn ≲ 10² yr for a_M = 10 R_⊕, several orders shorter than the timescale of the tidal evolution. This precondition is therefore expected to hold for the main lifetime of the Moon, including during the LHB.

On condition that the Moon rotates synchronously, its hemispheres can be defined. The near side is defined as the hemisphere that always faces Earth, and the far side is its opposite. The leading side is defined as the hemisphere that the lunar orbital velocity points to, and the trailing side is its opposite. The centers of the leading and trailing sides are the apex and the antapex points.

2.2. Encounter geometry

We first consider the encounter condition before exhibiting the details of the impact in the following subsections. In the heliocentric frame, Earth orbits the Sun circularly at 1 AU, the Moon is ignored at this pre-encounter stage, and the impactors are all massless particles with semi-major axes a_p> 1 AU and eccentricities 0 <e_p< 1. If the orbit of an impactor intersects that of Earth on the ecliptic, meaning that its perihelion distance q_p ≤ 1AU, an encounter is considered to occur at the mutual node. The relative velocity between the impactor and Earth when they encounter is the encounter velocity v_enc = v_p−v_⊕, where v_p and v_⊕ are orbital velocities of the impactor and Earth. Without loss of generality, the impactor’s and Earth’s arguments of perihelion are both assumed to be 0°, so that their velocities in heliocentric ecliptic coordinates are $\begin{array}{ccc}{{v}}_{\mathrm{p}}& \mathrm{=}& \sqrt{\frac{\mathit{G}{\mathit{M}}_{\mathrm{\odot }}}{{\mathit{a}}_{\mathrm{p}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{p}}^{\mathrm{2}}\mathrm{)}}}\mathrm{(}\mathrm{-}\sin {\mathit{f}}_{\mathrm{p}}\mathit{,}{\mathit{e}}_{\mathrm{p}}\mathrm{+}\cos {\mathit{f}}_{\mathrm{p}}\mathrm{)}\mathit{,}\\ {{v}}_{\mathrm{\oplus }}& \mathrm{=}& \sqrt{\frac{\mathit{G}{\mathit{M}}_{\mathrm{\odot }}}{{\mathit{a}}_{\mathrm{\oplus }}}}\mathrm{(}\mathrm{-}\sin {\mathit{f}}_{\mathrm{\oplus }}\mathit{,}\cos {\mathit{f}}_{\mathrm{\oplus }}\mathrm{)}\mathit{,}\end{array}$where G is the gravitational constant, M_⊙ is the solar mass, f_p and f_⊕ are true anomalies, and the terrestrial semi-major axis is a_⊕ = 1 AU. The mutual node is where f_p = f_⊕ = f_enc, making the impactor’s heliocentric distance equal to a_⊕, that is, $\cos {\mathit{f}}_{\mathrm{enc}}\mathrm{=}\frac{{\mathit{a}}_{\mathrm{p}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{p}}^{\mathrm{2}}\mathrm{)}}{{\mathit{a}}_{\mathrm{\oplus }}{\mathit{e}}_{\mathrm{p}}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{e}}_{\mathrm{p}}}\mathrm{·}$(6)Thus, the encounter velocity depends on a_p and e_p (or q_p): ${{v}}_{\mathrm{enc}}\mathrm{\left(}{\mathit{a}}_{\mathrm{p}}\mathit{,}{\mathit{e}}_{\mathrm{p}}\mathrm{\right)}\mathrm{=}{{v}}_{\mathrm{p}}\mathrm{\left(}{\mathit{a}}_{\mathrm{p}}\mathit{,}{\mathit{e}}_{\mathrm{p}}\mathit{,}{\mathit{f}}_{\mathrm{enc}}\mathrm{\right(}{\mathit{a}}_{\mathrm{p}}\mathit{,}{\mathit{e}}_{\mathrm{p}}\mathrm{\left)}\mathrm{\right)}\mathrm{-}{{v}}_{\mathrm{\oplus }}\mathrm{\left(}{\mathit{f}}_{\mathrm{enc}}\mathrm{\right(}{\mathit{a}}_{\mathrm{p}}\mathit{,}{\mathit{e}}_{\mathrm{p}}\mathrm{\left)}\mathrm{\right)}\mathit{.}$(7)For a fixed a_p, as shown in Fig. 1, the encounter speed v_enc is minimized when q_p = 1 AU. As q_p decreases (e_p increases), even though the magnitudes of the two orbital velocities v_p and v_⊕ in the moment of encounter are invariant, the angle between them expands, and thus the encounter speed v_enc increases. On the other hand, if q_p is fixed, larger a_p can also lead to higher v_enc, because the impactor orbital speed in the moment of encounter ${\mathit{v}}_{\mathrm{p}}\mathrm{=}\sqrt{\mathit{G}{\mathit{M}}_{\mathrm{\odot }}\left(\frac{\mathrm{2}}{{\mathit{a}}_{\mathrm{\oplus }}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{a}}_{\mathrm{p}}}\right)}$(8)is greater. As shown in Fig. 2, v_enc increases rightward (a_p increases) and downward (q_p decreases). Given a population of impactors, its typical v_enc is determined by its orbital distribution and can affect the cratering distribution. Therefore, different impactor populations represent different cratering asymmetries with a given target, providing a method for deriving impactor properties form the observed crater record.

Fig. 1 Encounter geometry for an impactor orbit with ap = 2.0 AU and qp = 0.2−1.0 AU. Where the Earth orbit (blue dashed circle) intersects the impactor orbit (black solid ellipse), the angle between v⊕ (blue dashed arrow) and vp (black solid arrow) decreases as qp increases, and thus venc (red solid arrow) shrinks.

2.3. Impact geometry

When the encounter between one impactor orbit and the Earth orbit occurs, as seen in the Earth–Moon system, the impactors in the common orbit are treated as particles approaching along the parallel straight lines in the direction of their common encounter velocity (Fig. 3a). They are assumed to be uniformly distributed in space, and to be numerous enough to cover the Moon’s orbit, that is, the gravitational cross section of the Moon is always maximized (lunar diameter in the planar model).

In the geocentric frame with the x-axis pointing to the lunar perigee, the direction of v_enc can be random. Since we average over the lunar period, it can be assumed to be in the positive y-axis direction without loss of generality. Figure 3a shows that the angle between the encounter velocity v_enc and the lunar orbital velocity v_M is the Moon’s true anomaly f_M, which varies uniformly at the rate of the mean motion n_M. The impact velocity is the relative velocity between the impactors and the Moon: $\begin{array}{ccc}& & {{v}}_{\mathrm{enc}}\mathrm{=}\mathrm{\left(}\mathrm{0}\mathit{,}{\mathit{v}}_{\mathrm{enc}}\mathrm{\right)}\mathit{,}\\ & & {{v}}_{\mathrm{M}}\mathrm{=}{\mathit{v}}_{\mathrm{M}}\mathrm{(}\mathrm{-}\sin {\mathit{f}}_{\mathrm{M}}\mathit{,}\cos {\mathit{f}}_{\mathrm{M}}\mathrm{)}\mathit{,}\\ & & {v}\mathrm{=}\mathrm{(}{\mathit{v}}_{\mathrm{M}}\sin {\mathit{f}}_{\mathrm{M}}\mathit{,}{\mathit{v}}_{\mathrm{enc}}\mathrm{-}{\mathit{v}}_{\mathrm{M}}\cos {\mathit{f}}_{\mathrm{M}}\mathrm{)}\mathit{.}\end{array}$Its magnitude, that is, the impact speed, is $\mathit{v}\mathrm{=}\sqrt{{\mathit{v}}_{\mathrm{enc}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathrm{M}}^{\mathrm{2}}\mathrm{-}\mathrm{2}{\mathit{v}}_{\mathrm{enc}}{\mathit{v}}_{\mathrm{M}}\cos {\mathit{f}}_{\mathrm{M}}}\mathit{.}$(12)At any instant, the impact can only occur on one hemisphere, and the impactors in this hemisphere have an equal v since they share the same v_enc. A normal impact point is the center of this bombarded hemisphere, where the incidence angle and thus the normal impact speed are both largest at this instant. As the impactors are assumed to be uniformly distributed, this is also where the impact flux is highest.

Fig. 2 Contour map of the encounter speed venc (km s-1) as a function of ap and qp.

Given the encounter speed v_enc and lunar orbital speed v_M, as long as v_enc>v_M, the impact speed v is always minimized when the Moon is at the perigee (f_M = 0°) and maximized at the apogee (f_M = 180°), where the antapex and the apex become the normal impact points. Figure 3 shows that the normal impact point constantly moves westward on the lunar surface at a varying rate, and that lower v_enc brings not only the lower v, but also a longer time spent on the leading side for the normal impact point, which implies that the lower the encounter speed, the stronger the leading/trailing asymmetry. That v_enc>v_M is taken for granted hereafter in this work. This holds for the MBAs, whose minimum v_enc is 7.9 km s^-1 (a_p = 2.5 AU and q_p = 1 AU), which is nearly equal to the maximum v_M when a_M = 1 R_⊕. It is also valid for those NEOs with a_p> 1.2 AU, whose minimum encounter speed v_enc = 2.4 km s^-1 (a_p = 1.2 AU and q_p = 1 AU) approximates to v_M for a_M = 11R_⊕, while a_M during the dominant epoch of the NEOs is at least about 40 R_⊕.

Fig. 3 Impact geometry seen in the rest frame of Earth a) and Moon b). a) The Moon is assumed to be where aM = 60 R⊕, with vM = 1.0 km s-1. The impactors, distributed extensively enough to cover the lunar orbit (black circle), are in the common direction of venc (black dashed arrow). Where fM = 0°,45°,...,315°, the lunar velocity vM (black solid arrow) and the impact velocity v for venc = 5 or 10 km s-1 (red or blue solid arrows) are all plotted to the same scale. b) Under the same conditions, v is plotted pointing at the normal impact points on the lunar surface.

Now we establish a rest frame of the Moon that is to be used in the integration. As shown in Fig. 4, the positive y-axis points to the apex and the minus x-axis points to Earth. In this frame, $\begin{array}{ccc}& & {v}\mathrm{=}\mathrm{\left(}{\mathit{v}}_{\mathit{x}}\mathit{,}{\mathit{v}}_{\mathit{y}}\mathrm{\right)}\mathit{,}\\ & & \\ & & \end{array}$while its magnitude v keeps its expression of Eq. (12). We introduce some variables: λ is the geometric longitude ranging from − 180° to + 180°, measured eastward from the center of near side; θ is the incidence angle ranging from 0° to 90°, the angle between v and the local horizon; S is the cross section. A unit area (length in planar model) at longitude λ on the bombarded hemisphere is $\mathrm{d}\mathit{l}\mathrm{=}{\mathit{R}}_{\mathrm{M}}\mathrm{d}\mathit{\lambda }$(16)and its cross section is $\mathrm{d}\mathit{S}\mathrm{=}{\mathit{R}}_{\mathrm{M}}\sin \mathit{\theta }\mathrm{d}\mathit{\lambda }\mathit{.}$(17)Denoting with e = (−cosλ,−sinλ) the normal vector, the incidence angle on the unit area can be derived with $\sin \mathit{\theta }\mathrm{=}\frac{{e}\mathrm{·}{v}}{\mathit{v}}\mathrm{=}\frac{{\mathit{v}}_{\mathit{x}}\cos \mathit{\lambda }\mathrm{+}{\mathit{v}}_{\mathit{y}}\sin \mathit{\lambda }}{\mathit{v}}\mathit{,}$(18)and the normal speed v_⊥, the normal component of v is ${\mathit{v}}_{\mathrm{\perp }}\mathrm{=}\mathit{v}\sin \mathit{\theta }\mathrm{=}{\mathit{v}}_{\mathit{x}}\cos \mathit{\lambda }\mathrm{+}{\mathit{v}}_{\mathit{y}}\sin \mathit{\lambda }\mathit{.}$(19)Denoting with ρ the uniform spatial density of the impactors, the impact flux, the number of impactors received per unit time by the unit area is $\mathrm{d}\mathit{F}\mathrm{=}\mathit{\rho v}\mathrm{d}\mathit{S}\mathrm{=}\mathit{\rho }{\mathit{R}}_{\mathrm{M}}{\mathit{v}}_{\mathrm{\perp }}\mathrm{d}\mathit{\lambda }\mathit{.}$(20)At any instant when the lunar true anomaly is f_M, the longitude of the normal impact point is defined as λ_⊥, where its normal vector e_⊥ is parallel to v, so that $\begin{array}{ccc}& \sin {\mathit{\lambda }}_{\mathrm{\perp }}\mathrm{=}\frac{{\mathit{v}}_{\mathit{y}}\mathrm{\left(}{\mathit{f}}_{\mathrm{M}}\mathrm{\right)}}{\mathit{v}\mathrm{\left(}{\mathit{f}}_{\mathrm{M}}\mathrm{\right)}}\mathit{,}& \cos {\mathit{\lambda }}_{\mathrm{\perp }}\mathrm{=}\frac{{\mathit{v}}_{\mathit{x}}\mathrm{\left(}{\mathit{f}}_{\mathrm{M}}\mathrm{\right)}}{\mathit{v}\mathrm{\left(}{\mathit{f}}_{\mathrm{M}}\mathrm{\right)}}\mathrm{·}\end{array}$(21)The bombarded hemisphere is then bounded by the longitude interval [ λ_l,λ_u ], where the lower and upper limits λ_l,u = λ_⊥ ∓ 90°, and $\begin{array}{ccc}& \sin {\mathit{\lambda }}_{\mathrm{l}\mathit{,}\mathrm{u}}\mathrm{=}\mathrm{\mp }\frac{{\mathit{v}}_{\mathit{x}}\mathrm{\left(}{\mathit{f}}_{\mathrm{M}}\mathrm{\right)}}{\mathit{v}\mathrm{\left(}{\mathit{f}}_{\mathrm{M}}\mathrm{\right)}}\mathit{,}& \cos {\mathit{\lambda }}_{\mathrm{l}\mathit{,}\mathrm{u}}\mathrm{=}\mathrm{\pm }\frac{{\mathit{v}}_{\mathit{y}}\mathrm{\left(}{\mathit{f}}_{\mathrm{M}}\mathrm{\right)}}{\mathit{v}\mathrm{\left(}{\mathit{f}}_{\mathrm{M}}\mathrm{\right)}}\mathrm{·}\end{array}$(22)As the bombarded hemisphere moves westward along the lunar equator, the moment a certain position λ enters this hemisphere and the moment it leaves are when λ becomes the lower and upper limits λ_l,u, respectively, and when sinθ = 0 there. Therefore, the time interval while a certain position λ is in the bombarded hemisphere, characterized by [ f_l,f_u ], can be derived with $\begin{array}{ccc}& \sin \mathit{\lambda }\mathrm{=}\mathrm{\mp }\frac{{\mathit{v}}_{\mathit{x}}\mathrm{\left(}{\mathit{f}}_{\mathrm{l}\mathit{,}\mathrm{u}}\mathrm{\right)}}{\mathit{v}\mathrm{\left(}{\mathit{f}}_{\mathrm{l}\mathit{,}\mathrm{u}}\mathrm{\right)}}\mathit{,}& \cos \mathit{\lambda }\mathrm{=}\mathrm{\pm }\frac{{\mathit{v}}_{\mathit{y}}\mathrm{\left(}{\mathit{f}}_{\mathrm{l}\mathit{,}\mathrm{u}}\mathrm{\right)}}{\mathit{v}\mathrm{\left(}{\mathit{f}}_{\mathrm{l}\mathit{,}\mathrm{u}}\mathrm{\right)}}\mathrm{·}\end{array}$(23)It turns out $\begin{array}{ccc}{\mathit{f}}_{\mathrm{l}}& \mathrm{=}& \arcsin \left(\frac{{\mathit{v}}_{\mathrm{M}}}{{\mathit{v}}_{\mathrm{enc}}}\sin \mathit{\lambda }\right)\mathrm{-}\mathit{\lambda ,}\\ {\mathit{f}}_{\mathrm{u}}& \mathrm{=}& \mathit{\pi }\mathrm{-}\arcsin \left(\frac{{\mathit{v}}_{\mathrm{M}}}{{\mathit{v}}_{\mathrm{enc}}}\sin \mathit{\lambda }\right)\mathrm{-}\mathit{\lambda }\mathit{.}\end{array}$

Fig. 4 Variables involved in integrating the cratering distribution. In the rest frame of the Moon, the common velocity of impactors is v (direction of the arrows). The hemisphere bounded by two positions where v (arrows on two sides) is parallel to the local horizons is the bombarded hemisphere, whose center, where v (arrow in the middle) is perpendicular to the lunar surface, is the normal impact point. Denotations are explained in the text.

2.4. Exact formulations

Fig. 5 Impact speed a); incidence angle b); normal speed c); crater diameter d); impact density e); and crater density f) as functions of longitude. Their variations with longitude (solid curves in red, orange, yellow, and green for venc = 10, 20, 30, and 40 km s-1) and their global averages (dashed lines in the same color for the same venc, except for those in panel b), which are all black to represent the invariant value) are calculated with vM = 1.0 km s-1, dmin = 0.5 km, αp = 1.75, ρ = 3.4 × 10-25 km-2, t = 3.7 Gyr, and dc = 25 km.

Here with v_enc and v_M given, the spatial distributions of impact density, impact speed, incidence angle, normal speed, crater diameter, and crater density are shown after integration. We note that the denotations of integrated variables are written uppercase to distinguish them from the instantaneous ones with lowercase letters.

The impact density N is the number of impacts on a certain region divided by its area. For the unit area on longitude λ illustrated in Fig. 4, its impact density during a unit time is $\mathrm{d}\mathit{N}\mathrm{=}\frac{\mathrm{d}\mathit{F}\mathrm{d}\mathit{t}}{\mathrm{d}\mathit{l}}\mathrm{=}\frac{\mathit{\rho }}{{\mathit{n}}_{\mathrm{M}}}\mathrm{(}{\mathit{v}}_{\mathit{x}}\cos \mathit{\lambda }\mathrm{+}{\mathit{v}}_{\mathit{y}}\sin \mathit{\lambda }\mathrm{)}\mathrm{d}{\mathit{f}}_{\mathrm{M}}\mathit{.}$(26)Integration of dN over the interval [ f_l,f_u ] leads to $\mathit{N}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\mathit{\rho }}{{\mathit{n}}_{\mathrm{M}}}\left[\sqrt{{\mathit{v}}_{\mathrm{enc}}^{\mathrm{2}}\mathrm{-}{\mathit{v}}_{\mathrm{M}}^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\lambda }}\mathrm{-}{\mathit{v}}_{\mathrm{M}}\sin \mathit{\lambda }\arccos \left(\frac{{\mathit{v}}_{\mathrm{M}}}{{\mathit{v}}_{\mathrm{enc}}}\sin \mathit{\lambda }\right)\right]\mathit{.}$(27)To simplify the expression, we define σ ∈ (0,π) by $\mathit{\sigma }\mathrm{=}\arccos \left(\frac{{\mathit{v}}_{\mathrm{M}}}{{\mathit{v}}_{\mathrm{enc}}}\sin \mathit{\lambda }\right)\mathit{,}$(28)which leads to $\begin{array}{ccc}\sin \mathit{\sigma }& \mathrm{=}& \mathrm{\left|}\cos \mathrm{\right(}{\mathit{f}}_{\mathrm{l}\mathit{,}\mathrm{u}}\mathrm{+}\mathit{\lambda }\mathrm{\left)}\mathrm{\right|}\mathrm{=}\sqrt{\mathrm{1}\mathrm{-}{\left(\frac{{\mathit{v}}_{\mathrm{M}}}{{\mathit{v}}_{\mathrm{enc}}}\sin \mathit{\lambda }\right)}^{\mathrm{2}}}\mathit{,}\\ \cos \mathit{\sigma }& \mathrm{=}& \sin \mathrm{(}{\mathit{f}}_{\mathrm{l}\mathit{,}\mathrm{u}}\mathrm{+}\mathit{\lambda }\mathrm{)}\mathrm{=}\frac{{\mathit{v}}_{\mathrm{M}}}{{\mathit{v}}_{\mathrm{enc}}}\sin \mathit{\lambda }\mathit{.}\end{array}$Therefore, $\mathit{N}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\mathit{\rho }{\mathit{v}}_{\mathrm{enc}}}{{\mathit{n}}_{\mathrm{M}}}\mathrm{(}\sin \mathit{\sigma }\mathrm{-}\mathit{\sigma }\cos \mathit{\sigma }\mathrm{)}\mathit{.}$(31)The impact speed V, incidence angle Θ, and normal speed V_⊥, as functions of longitude, are averages of v, θ, and v_⊥ over one lunar period, that is, $\mathit{V}\mathrm{=}\mathrm{\left(}{{}^{\mathrm{\int }}\mathit{f}_{\mathrm{l}}}_{{\mathit{f}}_{\mathrm{u}}}^{}\mathit{v}\mathrm{d}\mathit{N}\mathrm{\right)}\mathit{/}\mathit{N}$, $\sin \mathrm{\Theta }\mathrm{=}\mathrm{\left(}{{}^{\mathrm{\int }}\mathit{f}_{\mathrm{l}}}_{{\mathit{f}}_{\mathrm{u}}}^{}\sin \mathit{\theta }\mathrm{d}\mathit{N}\mathrm{\right)}\mathit{/}\mathit{N}$, and ${\mathit{V}}_{\mathrm{\perp }}\mathrm{=}\mathrm{\left(}{{}^{\mathrm{\int }}\mathit{f}_{\mathrm{l}}}_{{\mathit{f}}_{\mathrm{u}}}^{}{\mathit{v}}_{\mathrm{\perp }}\mathrm{d}\mathit{N}\mathrm{\right)}\mathit{/}\mathit{N}$. The integration leads to $\begin{array}{ccc}& & \\ & & \\ & & \end{array}$where $\begin{array}{ccc}{\mathit{I}}_{\mathrm{0}}& \mathrm{=}& \\ {\mathit{I}}_{\mathrm{1}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{3}{\mathit{v}}_{\mathrm{M}}}\mathrm{\left\{}\mathrm{\right[}\mathrm{(}{\mathit{v}}_{\mathrm{enc}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathrm{M}}^{\mathrm{2}}\mathrm{)}\cos \mathit{\lambda }\mathrm{+}\mathrm{2}{\mathit{v}}_{\mathrm{M}}{\mathit{v}}_{\mathrm{enc}}\sin \mathit{\sigma }\mathrm{]}{\mathit{I}}_{\mathrm{+}}\\ & & \mathrm{-}\mathrm{\left[}\mathrm{\right(}{\mathit{v}}_{\mathrm{enc}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathrm{M}}^{\mathrm{2}}\mathrm{)}\cos \mathit{\lambda }\mathrm{-}\mathrm{2}{\mathit{v}}_{\mathrm{M}}{\mathit{v}}_{\mathrm{enc}}\sin \mathit{\sigma }\mathrm{]}{\mathit{I}}_{\mathrm{-}}\\ & & \mathrm{+}\mathrm{(}{\mathit{v}}_{\mathrm{enc}}\mathrm{+}{\mathit{v}}_{\mathrm{M}}\mathrm{)}\mathrm{(}{\mathit{v}}_{\mathrm{enc}}^{\mathrm{2}}\mathrm{+}\mathrm{7}{\mathit{v}}_{\mathrm{M}}^{\mathrm{2}}\mathrm{)}\mathrm{\left(}\sin \mathit{\lambda }\mathrm{\right)}\mathrm{\Delta }\mathit{E}\\ & & \mathrm{-}\mathrm{(}{\mathit{v}}_{\mathrm{enc}}\mathrm{+}{\mathit{v}}_{\mathrm{M}}\mathrm{)}\mathrm{(}{\mathit{v}}_{\mathrm{enc}}\mathrm{-}{\mathit{v}}_{\mathrm{M}}{\mathrm{)}}^{\mathrm{2}}\mathrm{(}\sin \mathit{\lambda }\mathrm{\left)}\mathrm{\Delta }\mathit{F}\mathrm{\right\}}\mathit{,}\\ {\mathit{I}}_{\mathrm{2}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{3}{\mathit{v}}_{\mathrm{M}}^{\mathrm{2}}}\mathrm{\left\{}\mathrm{\right[}\mathrm{(}\mathrm{2}{\mathit{v}}_{\mathrm{enc}}^{\mathrm{2}}\mathrm{-}\mathrm{3}{\mathit{v}}_{\mathrm{M}}^{\mathrm{2}}\mathrm{)}\cos \mathit{\lambda }\mathrm{-}{\mathit{v}}_{\mathrm{M}}{\mathit{v}}_{\mathrm{enc}}\sin \mathit{\sigma }\mathrm{\left]}\mathrm{\right(}\sin \mathit{\lambda }\mathrm{)}{\mathit{I}}_{\mathrm{+}}\\ & & \mathrm{-}\mathrm{\left[}\mathrm{\right(}\mathrm{2}{\mathit{v}}_{\mathrm{enc}}^{\mathrm{2}}\mathrm{-}\mathrm{3}{\mathit{v}}_{\mathrm{M}}^{\mathrm{2}}\mathrm{)}\cos \mathit{\lambda }\mathrm{+}{\mathit{v}}_{\mathrm{M}}{\mathit{v}}_{\mathrm{enc}}\sin \mathit{\sigma }\mathrm{]}\mathrm{\left(}\sin \mathit{\lambda }\mathrm{\right)}{\mathit{I}}_{\mathrm{-}}\\ & & \mathrm{-}\mathrm{(}{\mathit{v}}_{\mathrm{enc}}\mathrm{+}{\mathit{v}}_{\mathrm{M}}\mathrm{)}\mathrm{\left[}\mathrm{\right(}{\mathit{v}}_{\mathrm{enc}}^{\mathrm{2}}\mathrm{-}\mathrm{2}{\mathit{v}}_{\mathrm{M}}^{\mathrm{2}}\mathrm{\left)}\cos \mathrm{\right(}\mathrm{2}\mathit{\lambda }\mathrm{)}\mathrm{+}\mathrm{3}{\mathit{v}}_{\mathrm{M}}^{\mathrm{2}}\mathrm{]}\mathrm{\Delta }\mathit{E}\\ & & \mathrm{+}\mathrm{(}{\mathit{v}}_{\mathrm{enc}}\mathrm{+}{\mathit{v}}_{\mathrm{M}}\mathrm{)}\mathrm{(}{\mathit{v}}_{\mathrm{enc}}\mathrm{-}{\mathit{v}}_{\mathrm{M}}{\mathrm{)}}^{\mathrm{2}}\cos \mathrm{(}\mathrm{2}\mathit{\lambda }\mathrm{\left)}\mathrm{\Delta }\mathit{F}\mathrm{\right\}}\mathit{,}\\ {\mathit{I}}_{\mathrm{3}}& \mathrm{=}& \\ {\mathit{I}}_{\mathrm{\pm }}& \mathrm{=}& \sqrt{{\mathit{v}}_{\mathrm{enc}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathrm{M}}^{\mathrm{2}}\cos \mathrm{\left(}\mathrm{2}\mathit{\lambda }\mathrm{\right)}\mathrm{\pm }\mathrm{2}{\mathit{v}}_{\mathrm{M}}{\mathit{v}}_{\mathrm{enc}}\cos \mathit{\lambda }\sin \mathit{\sigma }}\mathit{,}\\ \mathrm{\Delta }\mathit{E}& \mathrm{=}& \mathit{E}\left(\frac{\mathit{\lambda }\mathrm{-}\mathit{\sigma }}{\mathrm{2}}\mathrm{+}\frac{\mathit{\pi }}{\mathrm{4}}\mathrm{|}{\mathit{k}}^{\mathrm{2}}\right)\mathrm{-}\mathit{E}\left(\frac{\mathit{\lambda }\mathrm{+}\mathit{\sigma }}{\mathrm{2}}\mathrm{+}\frac{\mathit{\pi }}{\mathrm{4}}\mathrm{-}\mathit{\delta \pi }\mathrm{|}{\mathit{k}}^{\mathrm{2}}\right)\mathrm{-}\mathrm{2}\mathit{\delta E}\mathrm{(}{\mathit{k}}^{\mathrm{2}}\mathrm{)}\mathit{,}\\ \mathrm{\Delta }\mathit{F}& \mathrm{=}& \mathit{F}\left(\frac{\mathit{\lambda }\mathrm{-}\mathit{\sigma }}{\mathrm{2}}\mathrm{+}\frac{\mathit{\pi }}{\mathrm{4}}\mathrm{|}{\mathit{k}}^{\mathrm{2}}\right)\mathrm{-}\mathit{F}\left(\frac{\mathit{\lambda }\mathrm{+}\mathit{\sigma }}{\mathrm{2}}\mathrm{+}\frac{\mathit{\pi }}{\mathrm{4}}\mathrm{-}\mathit{\delta \pi }\mathrm{|}{\mathit{k}}^{\mathrm{2}}\right)\mathrm{-}\mathrm{2}\mathit{\delta F}\mathrm{(}{\mathit{k}}^{\mathrm{2}}\mathrm{)}\mathit{,}\\ {\mathit{k}}^{\mathrm{2}}& \mathrm{=}& \\ \mathit{\delta }& \mathrm{=}& \{\begin{array}{c}\\ & \\ & \end{array}\end{array}$Functions F(φ | k²) and F(k²) are the incomplete and complete elliptic integrals of the first kind; E(φ | k²) and E(k²) are those of the second kind.

The crater size generated by an impact depends on the impact speed, incidence angle, the projectile size, the target surface gravity, the densities of the two objects, and so on. Assuming the cratering efficiency of an oblique impact depends on the normal speed, the crater scaling law allows converting between crater diameter d_c and projectile diameter d_p with $\begin{array}{ccc}{\mathit{d}}_{\mathrm{p}}\mathrm{\left(}{\mathit{d}}_{\mathrm{c}}\mathit{,}{\mathit{v}}_{\mathrm{\perp }}\mathrm{\right)}& \mathrm{=}& \\ {\mathit{d}}_{\mathrm{c}}\mathrm{\left(}{\mathit{d}}_{\mathrm{p}}\mathit{,}{\mathit{v}}_{\mathrm{\perp }}\mathrm{\right)}& \mathrm{=}& \end{array}$The constants c_p,c, γ_p,c, and ${\mathit{\gamma }}_{\mathrm{p}\mathit{,}\mathrm{c}}^{\mathrm{\prime }}$ given by different studies vary, but their values are not important to the analytical derivation here. Still, they are needed for illustration and post-processing of numerical simulations, and we therefore adopt the forms given by Le Feuvre & Wieczorek (2011) for a non-porous regime (working for d_c> 20 km) with the Moon as the target: $\begin{array}{ccc}{\mathit{c}}_{\mathrm{p}}& \mathrm{=}& \\ {\mathit{\gamma }}_{\mathrm{p}}& \mathrm{=}& \\ {\mathit{\gamma }}_{\mathrm{p}}^{\mathrm{\prime }}& \mathrm{=}& \\ {\mathit{c}}_{\mathrm{c}}& \mathrm{=}& \\ {\mathit{\gamma }}_{\mathrm{c}}& \mathrm{=}& \\ {\mathit{\gamma }}_{\mathrm{c}}^{\mathrm{\prime }}& \mathrm{=}& \end{array}$where the non-porous scaling parameters are K = 1.17, ν₁ = 0.22, and ν₂ = 0.31, the lunar surface gravity is g = 1.6 m s^-2, the target density is ρ_t = 2.8 g cm^-3, the projectile density is ρ_p = 3 g cm^-3, and the transition diameter is d_∗ = 8.5 km for the Moon. Again we point out that the following derivation is based on Eqs. (44) and (45), but not on Eqs. (46)−(51).

The cumulative size distribution of an impactor population can be commonly assumed as ${\mathit{N}}_{\mathrm{p}}\mathrm{\propto }{\mathit{d}}_{\mathrm{p}}^{\mathrm{-}{\mathit{\alpha }}_{\mathrm{p}}}$, independent of its orbital distribution. Given a minimum of projectile size d_min, a normalized size distribution is $\mathit{N̅}\mathrm{p}\mathrm{(}\mathit{>}{\mathit{d}}_{\mathrm{p}}\mathrm{)}\mathrm{=}{\left(\frac{{\mathit{d}}_{\mathrm{p}}}{{\mathit{d}}_{\min }}\right)}^{\mathrm{-}{\mathit{\alpha }}_{\mathrm{p}}}\mathit{,} \mathrm{(}{\mathit{d}}_{\mathrm{p}}\mathrm{\ge }{\mathit{d}}_{\min }\mathrm{)}\mathit{.}$(52)It can be interpreted as the probability of an arbitrary impactor to be larger than d_p. The mean size of the projectiles is $\mathit{d̅}\mathrm{p}\mathrm{=}{\mathrm{\int }}_{\mathrm{+}\mathrm{\infty }}^{{\mathit{d}}_{\min }}{\mathit{d}}_{\mathrm{p}}\mathrm{d}\mathit{N̅}\mathrm{p}\mathrm{=}\frac{{\mathit{\alpha }}_{\mathrm{p}}}{{\mathit{\alpha }}_{\mathrm{p}}\mathrm{-}\mathrm{1}}{\mathit{d}}_{\min }\mathit{,}$(53)where the slope α_p> 1 (otherwise, the mean projectile size would be infinite). This equation should hold everywhere on the Moon since the size distribution is independent of orbital distribution, and gravitations on impactors are ignored. Therefore, the spatial distribution of the periodic averaged crater diameter is $\mathit{D}\mathrm{=}\mathrm{\left[}{{}^{\mathrm{\int }}\mathit{f}_{\mathrm{l}}}_{{\mathit{f}}_{\mathrm{u}}}^{}{\mathit{d}}_{\mathrm{c}}\mathrm{\right(}\mathit{d̅}{}_{\mathrm{p}}\mathit{,}{\mathit{v}}_{\mathrm{\perp }}\mathrm{\left)}\mathrm{d}\mathit{N}\mathrm{\right]}\mathit{/}\mathit{N}$, which is substituted with to avoid difficult analytical integration. The deviation of the substitute is no more than 1.5% for the ratio v_M/v_enc as high as 0.25 (the case a_M = 10 R_⊕ and v_enc = 10 km s^-1) and even lower for a lower speed ratio. Thus we derive $\mathit{D}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}{\mathit{c}}_{\mathrm{c}}\mathit{d̅}\begin{array}{c}{\mathit{\gamma }}_{\mathrm{c}}^{\mathrm{\prime }}\\ \mathrm{p}\end{array}\mathrm{\left[}{\mathit{V}}_{\mathrm{\perp }}\mathrm{\right(}\mathit{\lambda }\mathrm{)}{\mathrm{]}}^{{\mathit{\gamma }}_{\mathrm{c}}}\mathit{.}$(54)The crater density N_c( >d_c) is the number of craters with diameters larger than d_c per unit area. It is relevant to the age determination in the cratering chronology method. It depends on the projectile diameter d_p, the projectile size-frequency distribution , and the impact density N. On the assumption that every impact leaves one and only one crater, meaning that saturation, erosion, and secondary craters are all ignored, $\mathrm{d}{\mathit{N}}_{\mathrm{c}}\mathrm{(}\mathit{>}{\mathit{d}}_{\mathrm{c}}\mathrm{)}\mathrm{=}\mathit{N̅}\mathrm{p}\mathrm{(}\mathit{>}{\mathit{d}}_{\mathrm{p}}\mathrm{(}{\mathit{d}}_{\mathrm{c}}\mathit{,}{\mathit{v}}_{\mathrm{\perp }}\mathrm{\left)}\mathrm{\right)}\mathrm{d}\mathit{N}\mathit{.}$(55)In every unit time df_M ∈ [ f_l,f_u ], the factor dN gives the total density of craters on the longitude λ; function d_p gives the projectile size required to form a crater as large as d_c there; function gives the fraction of impactors larger than d_p, which is equivalent to the fraction of craters larger than d_c; and finally the product of and dN determines the required part of dN. The substitute for the integral ${{}^{\mathrm{\int }}\mathit{f}_{\mathrm{l}}}_{{\mathit{f}}_{\mathrm{u}}}^{}\mathit{N̅}{}_{\mathrm{p}}\mathrm{(}\mathit{>}{\mathit{d}}_{\mathrm{p}}\mathrm{(}{\mathit{d}}_{\mathrm{c}}\mathit{,}{\mathit{v}}_{\mathrm{\perp }}\mathrm{\left)}\mathrm{\right)}\mathrm{d}\mathit{N}$ is , whose deviation is no more than 0.07% for a ratio v_M/v_enc as high as 0.25 and even smaller for a lower speed ratio. Assuming d_c is large enough to ensure d_p(d_c,V_⊥(λ)) ≥ d_min holds everywhere on the lunar surface, then ${\mathit{N}}_{\mathrm{c}}\mathrm{(}\mathit{>}{\mathit{d}}_{\mathrm{c}}\mathit{,\lambda }\mathrm{)}\mathrm{=}{\left(\frac{{\mathit{d}}_{\min }}{{\mathit{c}}_{\mathrm{p}}{\mathit{d}}_{\mathrm{c}}^{{\mathit{\gamma }}_{\mathrm{p}}^{\mathrm{\prime }}}}\right)}^{{\mathit{\alpha }}_{\mathrm{p}}}\mathrm{\left[}{\mathit{V}}_{\mathrm{\perp }}\mathrm{\right(}\mathit{\lambda }\mathrm{\left)}{\mathrm{]}}^{{\mathit{\gamma }}_{\mathrm{p}}{\mathit{\alpha }}_{\mathrm{p}}}\mathit{N}\mathrm{\right(}\mathit{\lambda }\mathrm{)}\mathit{.}$(56)We note that N_c( >d_c,λ) describes not only the spatial distribution of the crater density, but also the craters’ size distribution. Equation (56) proves ${\mathit{N}}_{\mathrm{c}}\mathrm{\propto }{\mathit{d}}_{\mathrm{c}}^{\mathrm{-}{\mathit{\gamma }}_{\mathrm{p}}^{\mathrm{\prime }}{\mathit{\alpha }}_{\mathrm{p}}}$ on any lunar longitude, and thus the slopes of the size distributions of impactors and craters are related by ${\mathit{\alpha }}_{\mathrm{c}}\mathrm{=}{\mathit{\gamma }}_{\mathrm{p}}^{\mathrm{\prime }}{\mathit{\alpha }}_{\mathrm{p}}\mathit{.}$(57)Although the integration is made over one lunar period, all of the distributions but N(λ) and N_c( >d_c,λ) have no dependence on time, so that they are applicable to any time interval t (multiple of lunar period, theoretically) as long as v_M and v_enc are constant. We rewrite N(λ) after multiplying it by n_Mt/ (2π): $\mathit{N}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}\frac{\mathit{\rho t}{\mathit{v}}_{\mathrm{enc}}}{\mathit{\pi }}\mathrm{(}\sin \mathit{\sigma }\mathrm{-}\mathit{\sigma }\cos \mathit{\sigma }\mathrm{)}\mathit{.}$(58)This form describes the impact density after a bombardment duration t, but its relative variation does not differ at all from the one-periodic form, which is one of the reasons we suggest using the asymmetry amplitude to measure the relative variation (Sect. 2.7). Hereafter, this new expression of N(λ) takes the place of the previous one, while N_c( >d_c,λ) (Eq. (56)) is still valid.

All of the absolute spatial distributions are shown in Fig. 5 with solid curves. We caution that on each unit area, there must be impacts with all the (instantaneous) incidence angles θ ∈ [0°, 90°] theoretically, and thus all the normal speeds v_⊥ ∈ [0, v] and crater diameters d_c ∈ [0, ], nevertheless, their periodic averages are what the figure shows. The most significant common feature is that the curves all peak at the apex (λ = −90°) and are lowest at the antapex (λ = + 90°), regardless of v_enc. The leading/trailing asymmetry is unambiguously detected in distributions of all the investigated variables, which is to be confirmed again through approximate expressions (Sect. 2.7). Additionally, it is apparent that higher v_enc leads to an upward shift of all the curves, expect for Θ(λ). An increase in v_enc diminishes and enlarges the absolute amplitudes of D and N_c, while it almost does not influence those of N, V, and V_⊥. With v_M fixed, the absolute amplitude of Θ that only depends on the speed ratio v_M/v_enc also seems to decrease as v_enc increases. We provide the explanation in Sect. 2.7.

2.5. Near/far symmetry

Because N(λ) = N( ± 180°−λ), N(λ) is symmetric about λ = ± 90°, meaning that the distribution of the impact density on the lunar surface is symmetric about the line connecting the apex and the antapex. The same is easily found in terms of V, Θ, V_⊥, D, and N_c. Therefore, the Moon’s near and far sides are exact mirror images of each other, with no sign of a near/far asymmetry.

This conclusion is based on the assumption that the Earth’s gravitation on the impactors and its volume are not considered, meaning that Earth is treated as if it were transparent to the impactors. In reality, it may block impactors like a shield, preventing them from impacting the Moon, or gravitationally focusing them like a lens, increasing the flux that the near side receives. Bandermann & Singer (1973) confirmed the former effect for the condition that a_M< 25 R_⊕; Gallant et al. (2009) claimed that there was very little asymmetry for a_M in the range 10−50 R_⊕; Le Feuvre & Wieczorek (2011) found the asymmetry negligible with about 0.1% more craters of the near side than the far side, in contrast with Le Feuvre & Wieczorek (2005), who claimed a factor of four enhancement on the near side using impactors coplanar with Earth. Our deduction shows that without the Earth’s effect, the near/far asymmetry cannot exist, while the leading/trailing asymmetry is the inherent and natural consequence of cratering.

We therefore use β, the angular distance from the apex, to describe the leading/trailing asymmetry and not the longitude λ. It ranges between 0° and 180°, where the apex and antapex are, respectively. For the near side, where − 90° ≤ λ ≤ + 90°, $\begin{array}{ccc}\mathit{\lambda }\mathrm{=}\mathit{\beta }\mathrm{-}\mathrm{90}°\mathit{.}& & \end{array}$Substituting λ with β, the integrated variables as functions of β are $\begin{array}{ccc}& & \mathit{N}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}\mathrm{=}\frac{\mathit{\rho t}{\mathit{v}}_{\mathrm{enc}}}{\mathit{\pi }}\mathrm{(}\sin \mathit{\sigma }\mathrm{-}\mathit{\sigma }\cos \mathit{\sigma }\mathrm{)}\mathit{,}\\ & & \mathit{V}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}\mathrm{=}{\mathit{I}}_{\mathrm{1}}^{\mathrm{\prime }}\mathit{/}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\prime }}\mathit{,}\\ & & \sin \mathrm{\Theta }\mathrm{\left(}\mathit{\beta }\mathrm{\right)}\mathrm{=}{\mathit{I}}_{\mathrm{2}}^{\mathrm{\prime }}\mathit{/}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\prime }}\mathit{,}\\ & & {\mathit{V}}_{\mathrm{\perp }}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}\mathrm{=}{\mathit{I}}_{\mathrm{3}}^{\mathrm{\prime }}\mathit{/}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\prime }}\mathit{,}\\ & & \mathit{D}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}\mathrm{=}{\mathit{c}}_{\mathrm{c}}\mathit{d̅}\begin{array}{c}{\mathit{\gamma }}_{\mathrm{c}}^{\mathrm{\prime }}\\ \mathrm{p}\end{array}\mathrm{\left[}{\mathit{V}}_{\mathrm{\perp }}\mathrm{\right(}\mathit{\beta }\mathrm{)}{\mathrm{]}}^{{\mathit{\gamma }}_{\mathrm{c}}}\mathit{,}\\ & & {\mathit{N}}_{\mathrm{c}}\mathrm{(}\mathit{>}{\mathit{d}}_{\mathrm{c}}\mathit{,\beta }\mathrm{)}\mathrm{=}{\left(\frac{{\mathit{d}}_{\min }}{{\mathit{c}}_{\mathrm{p}}{\mathit{d}}_{\mathrm{c}}^{{\mathit{\gamma }}_{\mathrm{p}}^{\mathrm{\prime }}}}\right)}^{{\mathit{\alpha }}_{\mathrm{p}}}\mathrm{\left[}{\mathit{V}}_{\mathrm{\perp }}\mathrm{\right(}\mathit{\beta }\mathrm{\left)}{\mathrm{]}}^{{\mathit{\gamma }}_{\mathrm{p}}{\mathit{\alpha }}_{\mathrm{p}}}\mathit{N}\mathrm{\right(}\mathit{\beta }\mathrm{)}\mathit{,}\end{array}$where $\begin{array}{ccc}{\mathit{I}}_{\mathrm{0}}^{\mathrm{\prime }}& \mathrm{=}& \mathrm{2}\mathrm{(}\sin \mathit{\sigma }\mathrm{-}\mathit{\sigma }\cos \mathit{\sigma }\mathrm{)}\mathit{,}\\ {\mathit{I}}_{\mathrm{1}}^{\mathrm{\prime }}& \mathrm{=}& \frac{{\mathit{v}}_{\mathrm{enc}}}{\mathrm{3}\mathit{\eta }}\{\mathrm{\left[}\mathrm{\right(}\mathrm{1}\mathrm{+}{\mathit{\eta }}^{\mathrm{2}}\mathrm{)}\sin \mathit{\beta }\mathrm{+}\mathrm{2}\mathit{\eta }\sin \mathit{\sigma }\mathrm{]}{\mathit{I}}_{\mathrm{+}}^{\mathrm{\prime }}\\ & & \mathrm{-}\mathrm{\left[}\mathrm{\right(}\mathrm{1}\mathrm{+}{\mathit{\eta }}^{\mathrm{2}}\mathrm{)}\sin \mathit{\beta }\mathrm{-}\mathrm{2}\mathit{\eta }\sin \mathit{\sigma }\mathrm{]}{\mathit{I}}_{\mathrm{-}}^{\mathrm{\prime }}\\ & & \mathrm{-}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\eta }\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{7}{\mathit{\eta }}^{\mathrm{2}}\mathrm{)}\mathrm{\left(}\cos \mathit{\beta }\mathrm{\right)}\mathrm{\Delta }\mathit{E}\\ & & \mathrm{+}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\eta }\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\eta }{\mathrm{)}}^{\mathrm{2}}\mathrm{(}\cos \mathit{\beta }\mathrm{)}\mathrm{\Delta }\mathit{F}\}\mathit{,}\\ {\mathit{I}}_{\mathrm{2}}^{\mathrm{\prime }}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{3}{\mathit{\eta }}^{\mathrm{2}}}\{\mathrm{-}\mathrm{\left[}\mathrm{\right(}\mathrm{2}\mathrm{-}\mathrm{3}{\mathit{\eta }}^{\mathrm{2}}\mathrm{)}\sin \mathit{\beta }\mathrm{-}\mathit{\eta }\sin \mathit{\sigma }\mathrm{]}\mathrm{\left(}\cos \mathit{\beta }\mathrm{\right)}{\mathit{I}}_{\mathrm{+}}^{\mathrm{\prime }}\\ & & \mathrm{+}\mathrm{\left[}\mathrm{\right(}\mathrm{2}\mathrm{-}\mathrm{3}{\mathit{\eta }}^{\mathrm{2}}\mathrm{)}\sin \mathit{\beta }\mathrm{+}\mathit{\eta }\sin \mathit{\sigma }\mathrm{]}\mathrm{\left(}\cos \mathit{\beta }\mathrm{\right)}{\mathit{I}}_{\mathrm{-}}^{\mathrm{\prime }}\\ & & \mathrm{+}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\eta }\mathrm{)}\mathrm{\left[}\mathrm{\right(}\mathrm{1}\mathrm{-}\mathrm{2}{\mathit{\eta }}^{\mathrm{2}}\mathrm{\left)}\cos \mathrm{\right(}\mathrm{2}\mathit{\beta }\mathrm{)}\mathrm{-}\mathrm{3}{\mathit{\eta }}^{\mathrm{2}}\mathrm{]}\mathrm{\Delta }\mathit{E}\\ & & \mathrm{-}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\eta }\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\eta }{\mathrm{)}}^{\mathrm{2}}\cos \mathrm{(}\mathrm{2}\mathit{\beta }\mathrm{)}\mathrm{\Delta }\mathit{F}\}\mathit{,}\end{array}$$\begin{array}{ccc}& & {\mathit{I}}_{\mathrm{3}}^{\mathrm{\prime }}\mathrm{=}{\mathit{v}}_{\mathrm{enc}}\mathrm{\left[}\mathrm{\right(}\mathrm{1}\mathrm{+}\mathrm{2}{\cos }^{\mathrm{2}}\mathit{\sigma }\mathrm{)}\mathit{\sigma }\mathrm{-}\mathrm{3}\sin \mathit{\sigma }\cos \mathit{\sigma }\mathrm{]}\mathit{,}\\ & & {\mathit{I}}_{\mathrm{\pm }}^{\mathrm{\prime }}\mathrm{=}\sqrt{\mathrm{1}\mathrm{-}{\mathit{\eta }}^{\mathrm{2}}\cos \mathrm{\left(}\mathrm{2}\mathit{\beta }\mathrm{\right)}\mathrm{\pm }\mathrm{2}\mathit{\eta }\sin \mathit{\beta }\sin \mathit{\sigma }}\mathit{,}\\ & & \mathrm{\Delta }\mathit{E}\mathrm{=}\mathit{E}\left(\frac{\mathit{\beta }\mathrm{-}\mathit{\sigma }}{\mathrm{2}}\mathrm{|}{\mathit{k}}^{\mathrm{2}}\right)\mathrm{-}\mathit{E}\left(\frac{\mathit{\beta }\mathrm{+}\mathit{\sigma }}{\mathrm{2}}\mathrm{-}\mathit{\delta \pi }\mathrm{|}{\mathit{k}}^{\mathrm{2}}\right)\mathrm{-}\mathrm{2}\mathit{\delta E}\mathrm{(}{\mathit{k}}^{\mathrm{2}}\mathrm{)}\mathit{,}\\ & & \mathrm{\Delta }\mathit{F}\mathrm{=}\mathit{F}\left(\frac{\mathit{\beta }\mathrm{-}\mathit{\sigma }}{\mathrm{2}}\mathrm{|}{\mathit{k}}^{\mathrm{2}}\right)\mathrm{-}\mathit{F}\left(\frac{\mathit{\beta }\mathrm{+}\mathit{\sigma }}{\mathrm{2}}\mathrm{-}\mathit{\delta \pi }\mathrm{|}{\mathit{k}}^{\mathrm{2}}\right)\mathrm{-}\mathrm{2}\mathit{\delta F}\mathrm{(}{\mathit{k}}^{\mathrm{2}}\mathrm{)}\mathit{,}\\ & & {\mathit{k}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{4}\mathit{\eta }}{{\left(\mathrm{1}\mathrm{+}\mathit{\eta }\right)}^{\mathrm{2}}}\mathit{,}\\ & & \mathit{\delta }\mathrm{=}\{\begin{array}{c}\\ & \\ & \end{array}\\ & & \mathit{\sigma }\mathrm{=}\arccos \left(\mathrm{-}\mathit{\eta }\cos \mathit{\beta }\right)\mathit{,}\\ & & \sin \mathit{\sigma }\mathrm{=}\sqrt{\mathrm{1}\mathrm{-}{\left(\mathit{\eta }\cos \mathit{\beta }\right)}^{\mathrm{2}}}\mathit{,}\\ & & \cos \mathit{\sigma }\mathrm{=}\mathrm{-}\mathit{\eta }\cos \mathit{\beta ,}\\ & & \mathit{\eta }\mathrm{=}\frac{{\mathit{v}}_{\mathrm{M}}}{{\mathit{v}}_{\mathrm{enc}}}\mathrm{·}\end{array}$We note that the above equations are applicable to both the near and the far side.

Figure 6 illustrates the relative spatial variations of N, V, Θ, V_⊥, D, and N_c as functions of β with solid curves. The relative variation is the absolute variation divided by its global average (to be derived in Sect. 2.6). All variables decrease monotonically with increasing β, as we show in Fig. 5. Moreover, the greater the encounter speed v_enc (the lower the speed ratio v_M/v_enc with given v_M), the smaller the relative variation amplitudes (to be explained in Sect. 2.7).

2.6. Global averages

Here we derive the global averages of the impact density N, impact speed V, incidence angle Θ, normal speed V_⊥, crater diameter D, and crater density N_c through integrating. The impact number for a unit area on the longitude λ during a unit time is ${\mathrm{d}}^{\mathrm{2}}\mathit{C}\mathrm{=}\mathrm{d}\mathit{F}\mathrm{d}\mathit{t}\mathrm{=}\frac{\mathit{\rho }{\mathit{R}}_{\mathrm{M}}}{{\mathit{n}}_{\mathrm{M}}}\mathrm{(}{\mathit{v}}_{\mathit{x}}\cos \mathit{\lambda }\mathrm{+}{\mathit{v}}_{\mathit{y}}\sin \mathit{\lambda }\mathrm{)}\mathrm{d}\mathit{\lambda }\mathrm{d}{\mathit{f}}_{\mathrm{M}}\mathit{,}$(78)where F is the impact flux (Eq. (20)), v_x,y are elements of v (Eqs. (14) and (15)), and ρ is the uniform spatial density of impactors. The total number of impacts on the global lunar surface after one lunar period is the integral ${}^{!}\mathrm{d}^{\mathrm{2}}\mathit{C}$ over the λ interval [ λ_l,λ_u ] and f_M interval [ 0°,360° ] in turn, where λ_l,u(f_M) are boundaries of the instantaneous bombarded hemisphere (Eq. (22)). Thus, $\mathit{C}\mathrm{=}\frac{\mathrm{8}\mathit{\rho }{\mathit{R}}_{\mathrm{M}}}{{\mathit{n}}_{\mathrm{M}}}\mathrm{(}{\mathit{v}}_{\mathrm{enc}}\mathrm{+}{\mathit{v}}_{\mathrm{M}}\mathrm{)}\mathit{E}\mathrm{\left(}{\mathit{k}}^{\mathrm{2}}\mathrm{\right)}\mathit{.}$(79)The global average of N for the bombardment duration t (multiple of lunar period) is the total impact number for one period, C, divided by the global area l_gl = 2πR_M (perimeter for the planar model) and multiplied by n_Mt/ (2π): $\mathit{N̅}\mathrm{=}\frac{\mathrm{2}\mathit{\rho t}}{{\mathit{\pi }}^{\mathrm{2}}}\mathrm{(}{\mathit{v}}_{\mathrm{enc}}\mathrm{+}{\mathit{v}}_{\mathrm{M}}\mathrm{)}\mathit{E}\mathrm{\left(}{\mathit{k}}^{\mathrm{2}}\mathrm{\right)}\mathit{.}$(80)

Apparently, the global averages of V, Θ, and V_⊥ are the integrals ^!vd²C, ^!sinθd²C, and ${}^{!}\mathit{v}_{\mathrm{\perp }}{\mathrm{d}}^{\mathrm{2}}\mathit{C}$ over the above λ and f_M intervals divided by C, respectively. The integration results in

$\begin{array}{ccc}& & \\ & & \\ & & \end{array}$where is defined by $\sin \mathrm{\Theta ̅}\mathrm{=}\overline{\sin \mathit{\theta }}$, and it is found . The global averages of D and N_c are and , but to avoid analytical integration, we again substitute them by and , functions (Eqs. (45) and (44)) of and . The substitutes are quite acceptable, for their deviations are <3% and <1%, respectively, for a speed ratio v_M/v_enc as high as 0.25 (the case a_M = 10 R_⊕ and v_enc = 10 km s^-1). A lower speed ratio leads to even smaller deviations. Thus, we derive $\begin{array}{ccc}\mathit{D̅}& \mathrm{=}& \\ \mathit{N̅}\mathrm{c}& \mathrm{=}& \end{array}$The global averages are shown in Fig. 5 with horizontal dashed lines. They all exhibit positive relations with v_enc, except for . Equation (82) explains the exception by determining the global average of sinθ to be π/ 4 always, equivalent to , regardless of both the encounter speed v_enc and the lunar orbital speed v_M.

2.7. Approximate formulations

We have formulated the spatial distributions N(β), V(β), Θ(β), V_⊥(β), D(β), and N_c( >d_c,β) (Eqs. (59)−(64)) and the global averages , , , , , and (Eqs. (80)−(85)). Hereafter, we use Γ to denote any one (or every one) of the variables N, V, Θ, V_⊥, D, and N_c, and use to denote its global average. The normalized distribution is defined as the absolute distribution divided by the global average, which describes the relative variation $\mathrm{\Delta \Gamma }\mathrm{\left(}\mathit{\beta }\mathrm{\right)}\mathrm{=}\mathrm{\Gamma }\mathrm{\left(}\mathit{\beta }\mathrm{\right)}\mathit{/}\mathrm{\Gamma ̅}\mathit{.}$(86)The above variables represent every aspect of cratering and their spatial distributions are what we call cratering distribution. Each of the formulations Γ(β) and except for the constant can be taken as the product of two factors involving the encounter speed v_enc and the speed ratio η = v_M/v_enc respectively. Recalling the assumption v_M<v_enc, we can derive their series expansions around η = 0. To first order, the formulations are simplified to be

$\begin{array}{ccc}& & \\ & & \end{array}$$\begin{array}{ccc}& & \\ & & \\ & & \\ & & {\mathit{N}}_{\mathrm{c}}\left(\mathit{\beta }\right)\mathrm{=}\frac{\mathit{\rho t}{\mathit{v}}_{\mathrm{enc}}}{\mathit{\pi }}{\left[\frac{{\mathit{d}}_{\min }}{{\mathit{c}}_{\mathrm{p}}{\mathit{d}}_{\mathrm{c}}^{{\mathit{\gamma }}_{\mathrm{p}}^{\mathrm{\prime }}}}{\left(\frac{\mathit{\pi }{\mathit{v}}_{\mathrm{enc}}}{\mathrm{4}}\right)}^{{\mathit{\gamma }}_{\mathrm{p}}}\right]}^{{\mathit{\alpha }}_{\mathrm{p}}}\\ & & \\ & & \\ & & \\ & & \\ & & \\ & & \end{array}$All of the approximate distributions are shown in Fig. 6 with dashed curves for comparison with the exact ones. Smaller η leads to better approximation. Except for N_c(β) and , the errors of the approximate forms are ≲1% for η = 0.1 and ≲5% even for η = 0.25. In particular, the errors of V_⊥(β) are merely 0.04% and 0.25%, respectively. The approximations of N_c(β) and also depend on α_p. Given η = 0.1, their errors are 1%−5% and 0.7%−1.6% with α_p increasing from 1 to 3, which is acceptable. However, for extreme conditions when η and α_p are both large, approximate forms of N_c(β) and should be used with caution. For reference, the current Earth–Moon distance a_M = 60 R_⊕ (v_M = 1.0 km s^-1) and the current impactor population of NEOs (v_enc ≃ 20 km s^-1 according to Gallant et al. 2009) result in an η of only about 0.05.

Fig. 6 Normalized impact speed a); incidence angle b); normal speed c); crater diameter d); impact density e); and crater density f) as functions of β. The exact variations with β (solid curves in red, orange, yellow, and green for venc = 10, 20, 30, and 40 km s-1) are calculated using exact formulations with vM = 1.0 km s-1, dmin = 0.5 km, αp = 1.75, ρ = 3.4 × 10-25 km-2, t = 3.7 Gyr, and dc = 25 km. The approximate variations (dashed curves in the same color as the solid curves for the same venc) are calculated with the same parameters but using the approximate series of formulations. The fit variations (dotted curves in the same color as the solid curves for the same venc) are best fits of exact variations to the approximate formulations. All the variations are in terms of the relevant exact global averages (horizontal dotted lines).

There are a few points about the approximate distributions to note. First, all the spatial distributions can be described by a simple function of β, $\mathrm{\Gamma }\mathrm{\left(}\mathit{\beta }\mathrm{\right)}\mathrm{=}{\mathit{A}}_{\mathrm{0}}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{A}}_{\mathrm{1}}\cos \mathit{\beta }\mathrm{)}\mathit{,}$(98)where A₀> 0 and 0 <A₁< 1. This function is symmetric about the point (90°,A₀) and has one maximum A₀(1 + A₁) at the apex (β = 0°) and one minimum A₀(1−A₁) at the antapex (β = 180°), with a monotonic decrease from the former to the latter. Thus, the leading/trailing asymmetry is again proved here through analysis.

Second, the two parameters A₀ and A₁ entirely determine the variation of Γ. It is true that for every Γ, A₀ is exactly the approximate and also the value of Γ on the prime meridian (β = 90°). Additionally, A₀ has a positive relation with v_enc (except for Θ), explaining the dependence of on v_enc shown in Fig. 5. (The effect of η on the exact is negligible.) The factor (1 + A₁cosβ) is then equal to the normalized distribution ΔΓ. The second parameter A₁ just determines the amplitude of the relative variation shown in Fig. 6 and therefore is called “asymmetry amplitude”. For all Γ, it always holds that A₁ ∝ η, that is, the higher the v_enc or the lower the v_M (the larger a_M), the fainter the leading/trailing asymmetry. This relation can be understood by imaging the movement of the normal impact point along the lunar equator, whose rate is dλ_⊥/ dt. If it were moving uniformly, that is, if η → 0 and thus dλ_⊥/ dt → −df_M/ dt = −n_M, then the conditions (impact number, size distribution of craters, etc.) in every unit area would be the same, leading to a uniform spatial distribution. The varying moving rate adds to the number of normal impacts near the apex, and thus results in the biased cratering in all aspects. In this way, η determines A₁.

Third, some simple relations between A_0,1 of different variables are found. Parameters A_0,1 of N, V, Θ, V_⊥, D, and N_c are denoted with superscripts N, V, Θ, ⊥, D, and c, respectively. It is seen that $\begin{array}{ccc}{\mathit{A}}_{\mathrm{0}}^{\mathrm{\perp }}& \mathrm{=}& \\ {\mathit{A}}_{\mathrm{0}}^{\mathit{D}}& \mathrm{=}& \\ {\mathit{A}}_{\mathrm{0}}^{\mathrm{c}}& \mathrm{=}& \\ {\mathit{A}}_{\mathrm{1}}^{\mathit{D}}& \mathrm{=}& \\ {\mathit{A}}_{\mathrm{1}}^{\mathrm{c}}& \mathrm{=}& \end{array}$These are the results of the relations between the variables, and the first three equations involving A₀, the approximate , are identical to the connections between V, Θ, V_⊥, D, and N_c themselves (Eqs. (54) and (56)). Furthermore, the last two relations between A₁, the asymmetry amplitudes, are degenerate relations between A₀: the product of (A^⊥)^γ_pα_p and A^N becomes their sum, and then the exponents of A^⊥, both γ_c and γ_pα_p, become its coefficients. The absolute distributions of D and N_c both depend on the impactors’ size distribution, which is characterized by α_p and d_min, but ${\mathit{A}}_{\mathrm{1}}^{\mathit{D}}$ is independent of this, and ${\mathit{A}}_{\mathrm{1}}^{\mathrm{c}}$ only involves α_p. That ${\mathit{A}}_{\mathrm{1}}^{\mathrm{c}}$ has no dependence on d_c means that the normalized crater density distribution is not a function of d_c: .

Last but not the least, the above facts inspire the observation. Except for D and N_c, the other four, N, V, Θ, and V_⊥, are unobservable in surveys of crater record. Now the dependence of ${\mathit{A}}_{\mathrm{1}}^{\mathit{D,}\mathrm{c}}$ on ${\mathit{A}}_{\mathrm{1}}^{\mathit{N,}\mathrm{\perp }}$ provides a way to directly obtain the relative variations of N and V_⊥. Furthermore, when η is derived from the normalized distributions of D and N_c, those of all the unobservable variables can be obtained. Since ${\mathit{A}}_{\mathrm{1}}^{\mathrm{c}}$ has no dependence on d_c, when observed craters are counted on the Moon, it is better to choose a small d_c for accuracy based on more data, as long as the exclusion of saturation and secondaries is ensured and d_p(d_c,V_⊥(β)) >d_min holds (Sect. 2.4). Conversely, we can also determine whether secondary craters are included in a sample or whether a sample is complete by verifying if is different from when d_c is greater. It is even more meaningful that formulations D(β) and N_c( >d_c,β) enable us to derive from the crater record the bombardment conditions, which includes not only the impactor properties and bombardment duration, but also the lunar orbit in the dominant epoch of the impactors (Sect. 2.8).

2.8. Reproducing bombardment conditions

Fig. 7 Reproduction error of the speed ratio η. For a given exact η, after fitting the approximate formulations to the exact cratering distribution (assuming αp = 1.75) and using the theoretical relations between asymmetry amplitude and speed ratio to reproduce the latter with the fits of the former, the reproduction error is calculated as the relative difference between the reproduced and the exact η. The reproduction errors that the spatial distributions of V, Θ, V⊥, D, N, and Nc lead to (red, orange, yellow, green, cyan, and blue curves) are compared. Where venc = 10, 20, 30, and 40 km s-1 on condition vM = 1.0 km s-1 are indicated (vertical dot lines).

To verify the validity of the approximate formulations (Eqs. (87)−(92)), we fit them using the least-squares method to the exact distributions (Eqs. (59)−(64)) that we considered as observed data. The fit cratering distributions are shown in Fig. 6 with dotted curves. Given α_p, which only influences the relative variation of N_c, it is seen that the smaller the η, the smaller the difference between the fit and exact distributions, because smaller η leads to a better approximation at first. The errors of fit distributions are <0.3% when η = 0.1 and α_p = 1.75, almost half the errors of the approximate ones, except that for Θ they are nearly equal. The best fits are for the V_⊥ distribution, with errors <0.02%. Furthermore, given the fits of A₁, η can be easily reproduced based on the relations between A₁ and η. Figure 7 shows how the reproduced η varies with the exact η. When η = 0.1 and α_p = 1.75, the errors of the reproduced η are lower than 1% for all the variables. As expected, smaller η leads to a better reproduction of itself. Moreover, fitting distributions of different variables results in different reproductions: Θ and N_c (given α_p = 1.75) give relatively large errors, while V_⊥ gives the smallest; Θ and D tend to overestimate, while N, V, and N_c tend to underestimate. We caution that the observational error and limited data will prevent reaching such good fits and reproductions in reality, but our estimation shows what the series of approximate formulations is capable of.

The problem is that of N, V, Θ, V_⊥, D, and N_c, the first four are not observables of formed craters (except elliptic rim of a crater may imply extreme oblique impact). We therefore propose a method for reproducing the bombardment conditions (the lunar orbit and its impactor population) when a given crater record was formed through combined observations of D and N_c. By fitting formulations D(β) and N_c(β) to the observed spacial distributions of D and N_c and then solving the system of equations, $\begin{array}{ccc}{\mathit{A}}_{\mathrm{0}}^{\mathit{D}}& \mathrm{=}& \\ {\mathit{A}}_{\mathrm{0}}^{\mathrm{c}}& \mathrm{=}& \\ {\mathit{A}}_{\mathrm{1}}^{\mathit{D}}& \mathrm{=}& \\ {\mathit{A}}_{\mathrm{1}}^{\mathrm{c}}& \mathrm{=}& \end{array}$where the left-hand sides are fits of ${\mathit{A}}_{\mathrm{0}\mathit{,}\mathrm{1}}^{\mathit{D}}$ and ${\mathit{A}}_{\mathrm{0}\mathit{,}\mathrm{1}}^{\mathrm{c}}$, we may reproduce the complete set of bombardment conditions: the lunar orbital speed v_M and the impactors’ encounter speed v_enc, slope of the size distribution α_p, minimum diameter d_min, impact flux F, and dominant duration t. We note that F = 2R_Mρv_enc is used instead of ρ only for convenience in observation, that (Eq. (53)), and that c_c,p, γ_c,p, ${\mathit{\gamma }}_{\mathrm{c}\mathit{,}\mathrm{p}}^{\mathrm{\prime }}$, (Eqs. (46)–(51)) l_gl, and d_c are all constants. First, even without any knowledge of the above conditions, α_p and η = v_M/v_enc can be directly derived from Eqs. (106) and (107). Second, with α_p and η obtained, knowing any of v_M, v_enc, d_min, and Ft results in acquirement of all of them, because v_M and v_enc are related by η, while v_enc, d_min, and Ft are related through Eqs. (104) and (105). Third, with Ft obtained, knowing either t or F leads to the reproduction of the other.

We now show a sample reproduction following the above procedure. The exact distributions of D and N_c are taken again as the observed ones under current conditions where the lunar impactor population are the NEOs (Strom et al. 2005): v_M = 1.022 km s^-1, v_enc = 20 km s^-1 (Gallant et al. 2009), α_p = 1.75 (Bottke et al. 2002), d_min = 0.5 km, F = 7.5 × 10^-13 yr^-1 (for craters larger than 1 km; Le Feuvre & Wieczorek 2011), and t = 3.7 Gyr (age of the Orientale basin; Le Feuvre & Wieczorek 2011). The value of d_min and d_c = 25 km are chosen for consistency with the working range of Eqs. (45) and (44). We note that our synthetic distributions do not represent the real observation because the real minimum size of NEOs is smaller than what we assume here, and the impact flux is defined to take all the craters into account and not only craters larger than 1 km. However, the goodness of reproduction are not affected. Fitting the synthetic distributions that each involve 181 data points (uniform β interval of 1°) results in ${\mathit{A}}_{\mathrm{0}}^{\mathit{D}}\mathrm{=}\mathrm{33.0}$ km, ${\mathit{A}}_{\mathrm{1}}^{\mathit{D}}\mathrm{=}\mathrm{0.0238}$, ${\mathit{A}}_{\mathrm{0}}^{\mathrm{c}}\mathrm{=}\mathrm{1.63}\mathrm{\times }{\mathrm{10}}^{-7}{\mathrm{km}}^{-1}$, and ${\mathit{A}}_{\mathrm{1}}^{\mathrm{c}}\mathrm{=}\mathrm{0.129}$ (with uncertainties ≲0.1%), among which ${\mathit{A}}_{\mathrm{1}}^{\mathit{D}}$ and ${\mathit{A}}_{\mathrm{1}}^{\mathrm{c}}$ lead to η = (5.112 ± 0.003) × 10^-2 and α_p = 1.738 ± 0.003. If v_enc is known, then v_M = 1.0224 ± 0.0006 km s^-1, d_min = 0.495 ± 0.004 km, and Ft = (2.82 ± 0.04) × 10^-3. Furthermore, if t is given, then F = (7.6 ± 0.1) × 10^-13 yr^-1. The reproduction errors, differences between the best-fit and the exact values, are 0.05%, 0.16%, 0.82%, and 1.44% for v_M, α_p, d_min, and F, respectively. The goodness of reproduction sufficiently exhibits the viability of the reproducing method.

3. Numerical simulations

As indicated by Strom et al. (2015), the contribution of MBAs to the crater density exceeds that of NEOs by more than an order of magnitude for craters larger than 10 km in diameter, therefore it is quite important to determine the consequence of the bombardment of MBAs. We numerically simulated this process in the background of the LHB, dominance of the MBAs, with the Earth–Moon distance a_M varying in different cases. We found both the leading/trailing and the pole/equator cratering asymmetries. We established the formulation of the coupled-asymmetric distribution. In addition, the analytical model and numerical simulations confirm each other in various aspects.

3.1. Numerical model

Our numerical model includes the Sun, Earth, the Moon, and 3 × 10⁵ small particles in each case, which represent the primitive MBAs. All the large bodies were treated as rigid spheres. Only gravitation was involved and the effects between particles were ignored. Efforts were made to produce the circumstances of the LHB and a relatively realistic orbit of the Moon whose eccentricity and inclination are considered, so that the simplification of our analytical model can be assessed. The initialization was set at the moment after the MBAs had been disturbed, so that the giants can be excluded. Still, we caution that the orbits of the MBAs are artificially biased and are not modified by the giant planets.

The initial orbit of Earth was assumed to be the same as today. In the ecliptic coordinate system, its semi-major axis was a_⊕ = 1 AU, the eccentricity was e_⊕ = 0.0167, the inclination was i_⊕ = 0°, the longitude of the ascending node was Ω_⊕ = 348.7°, the argument of perihelion was ω_⊕ = 114.2°, and the mean anomaly was M_⊕ = 0°. Because the lunar early history involves a great uncertainty as mentioned in Sect. 1, five cases were set with the lunar semi-major axis a_M = 20, 30, 40, 50, and 60 R_⊕ in turn, probably covering most durations of the lunar history. This also helps to examine the influence of a_M on the cratering. The lunar eccentricity e_M was initially set to 0.04 (Ross & Schubert 1989), while the current value of the lunar inclination to the ecliptic i_M = 5° was adopted. The Moon’s longitude of the ascending node Ω_M, the argument of perihelion ω_M, and the mean anomaly M_M were all 0°, because the lunar orbital period and the period of its precession along the ecliptic are extremely short relative to the timescale of the tidal evolution.

Primitive MBAs mean the asteroids that occupied the region where the contemporary main belt is during the LHB, that is, 2.0 AU ≤ a_p ≤ 3.5 AU. Their orbits were disturbed by the migration of giants (Gomes et al. 2005), while their size distribution can be taken as invariant (Bottke et al. 2005). We only needed those MBAs that have the potential for encountering the Earth–Moon system, namely, those with perihelion distances q_p no larger than 1 AU. To maximize the impact probability (Morbidelli & Gladman 1998) and avoid the interference of the effect of q_p, the initial q_p was constrained to be 1 AU, resulting in 0.50 ≤ e_p ≤ 0.71. The inclination was set to 0° ≤ i_p ≤ 0.5° to examine the mechanism of the latitudinal cratering asymmetry. The orbital elements Ω_p, ω_p, and M_p all ranged from 0 to 360°. Each particle’s initial orbital elements were randomly generated in corresponding ranges following uniform distributions, except that its eccentricity was derived from e_p = 1−q_p/a_p. The mass of every particle was 4 × 10¹² kg, equivalent to a diameter d_p ~ 1 km given a density of 3 g cm^-3. This means that the size distribution is not considered at the stage of N-body simulations. The total mass of the particles of 1 × 10¹⁸ kg is negligible to that of the Moon, therefore generating a size distribution in the simulations will not lead to statistically different results, but will take many more simulations to derive meaningful statistics. Therefore, we only considered this in the post-processing when the crater diameter and crater density are calculated. Although the size distribution of current MBAs has been studied by several surveys (Sect. 1), it cannot be easily modeled because of its power-law breaks. In this work, we assumed the normalized size distribution (Eq. (52)) to be a single power-law characterized by α_p = 3 and d_min = 5 km referring to Ivezić et al. (2001), who determined α_p = 3 for 5 km ≲ d_p ≲ 40 km.

Three-dimensional N-body simulations were run to integrate the orbits of the Sun, Earth, the Moon, and 3 × 10⁵ particles in each case. The integration time was 10 Myr as a compromise between computing efficiency and the fact that the LHB lasted for ~0.1 Gyr (Gomes et al. 2005). In addition, this is long enough to allow nearly all the potential impacts to occur, and it is short enough to omit the tidal evolution. The integration consists of two steps. We first integrated Earth, the Moon, and all the particles in the heliocentric ecliptic reference frame using the hybrid symplectic integrator in the Mercury software package (Chambers 1999), which was altered to adapt it to our purpose. When a particle had a close encounter with the Moon, the information of Earth, the Moon, and this particle was recorded. Then in the next step, we used the Runge-Kutta-Fehlberg method to integrate the recorded particles and the Earth–Moon system until all the particles had impacted the Moon. In this step, the perturbation of the Sun was ignored since the integration time it took is only ~10 min, while the lunar orbital period is five days (a_M = 20 R_⊕) at least.

Given that the Moon rotated synchronously during the LHB (Sect. 2.1), every impact in the simulations was precisely located on the lunar surface. Because the Moon has a longitudinal geometrical libration due to its eccentricity (the lunar spin axis was assumed to be perpendicular to the lunar orbital plane so that there is no libration in latitude), the near side is defined as the hemisphere facing Earth only when the Moon is at its perigee and the far side is its opposite. The longitude line that the center of the near side lies on is the prime meridian, with the east longitude values being positive. The apex point, the center of the leading side, is at (− 90°,0°), while the antapex point, the center of the trailing side, is at (90°,0°). The low- and high-latitude regions are also defined here as the area with a latitude between ± 45° and the remaining area on the lunar surface, respectively.

Based on this numerical model, the bombardment conditions are v_enc = 8.32 km s^-1 (estimated with Eq. (7) given q_p = 1 AU and a_p = 2.75 AU, the mean of the initial a_p), α_p = 3, d_min = 5 km ( km assuming that the size distribution of the impactors is globally invariant), t = 10⁷ yr, and v_M = 1.78, 1.45, 1.26, 1.12, and 1.03 km s^-1 (with e_M ignored) for cases 1−5. Hereafter the above modeled bombardment conditions and the results we derived from them using our analytical formulations are called “predicted”, while those directly derived form simulations are called “simulated”.

3.2. Preliminary statistics

Fig. 8 Impact distribution on the lunar surface for the case aM = 20 R⊕ with 1544 impacts in total. Every impact in the simulations is shown at the location it occurs with a point, whose size and color represent the magnitudes of impact speed and incidence angle for this impact.

Fig. 9 Regional impact speed a); incidence angle b); normal speed c); crater diameter d); impact density e); and crater density f) as functions of the lunar semi-major axis. The black squares connected by black solid lines are the global averages. The red squares connected by the red dashed and dotted lines are the averages on the near and far sides. The blue squares connected by the blue dashed and dotted lines are the averages on the leading and trailing sides. The green squares connected by the green dashed and dotted lines are the averages on the low- and high-latitude regions.

Cases 1−5 with a_M = 20, 30, 40, 50, and 60 R_⊕ in turn result in the global impact number C_gl = 1544, 1413, 1482, 1388, and 1391, respectively. Figure 8 shows the impact distribution on the lunar surface for case 1, and those of the other cases are similar. Every impact location was recorded along with the impact condition, namely, the impact speed v and the incidence angle θ. After this, the normal speed v_⊥ = vsinθ and the size of the formed crater centered on this impact location were calculated. In a specific area, the averages of v, θ, v_⊥, and d_c are denoted by V, Θ, V_⊥, and D, following Sect. 2; the impact density N is the number of impacts C divided by the area, while the crater density N_c is the number of craters larger than a given size d_c divided by the area. We point out that this crater number was obtained by totalling the probability of generating a crater larger than d_c for every impact in this area, that is, ${\sum }_{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{C}}\mathit{N̅}{}_{\mathrm{p}}\mathrm{(}\mathit{>}{\mathit{d}}_{\mathrm{p}}\mathrm{(}{\mathit{d}}_{\mathrm{c}}\mathit{,}{\mathit{v}}_{\mathrm{\perp }\mathit{,j}}\mathrm{\left)}\mathrm{\right)}$. This is how the size distribution of impactors is involved in post-processing. The given size d_c was set to 100 km to ensure d_p(d_c,v_⊥) ≥ d_min, otherwise the descriptions of N_c distribution does not hold (Sect. 2.4). We note that the assumed values of and d_c do not influence the relative spatial variations of D and N_c at all.

We summarize the regional Γ (denoting every/any one of V, Θ, V_⊥, D, N, and N_c following Sect. 2) for cases 1−5 in Fig. 9. The impact density of the leading side N_ldg is much higher than that of the trailing side N_trg for all the cases, with excesses of 25%−49% (equivalent to the excesses of C_ldg to C_trg). Other Γ_ldg also remain larger than Γ_trg with notable excesses. These are clear signs of the leading/trailing asymmetry, indicating its existence in all the investigated aspects of cratering at any a_M. It also shows unambiguously that except for V, the averages of the low-latitude region Γ_low are always higher than those of high-latitude region Γ_high, which confirms the pole/equator asymmetry regardless of a_M. We note that of the investigated variables, V is the most sensitive to the leading/trailing asymmetry, but does not involve a pole/equator asymmetry at all, while Θ has the greatest dominance in the pole/equtor asymmetry over the leading/trailing asymmetry. No one pair of Γ_near and Γ_far has a distinct and lasting discrepancy. For example, N_near is 8% smaller than N_far for case 2, but 1%−5% greater for other cases (the same for C_near and C_far). Therefore, the near/far asymmetry is clearly negligible for the whole range 20−60 R_⊕ even if it is present, in accordance with Gallant et al. (2009) and Le Feuvre & Wieczorek (2011), who adopted the current impactors.

It is apparent that monotonically decreases from 9.15 to 8.80 km s^-1 with a_M increasing. The predicted (Eq. (81)) is 8.60−8.41 km s^-1 for a_M = 20−60 R_⊕, only 6%−4% smaller than simulated values. Although we consider this already an excellent agreement, we continue to take the gravitational acceleration of the Moon into account, which is ignored in our analytical model, by substituting v_enc = 8.32 km s^-1 with $\sqrt{{\mathit{v}}_{\mathrm{enc}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathrm{esc}}^{\mathrm{2}}}\mathrm{=}\mathrm{8.65}$ km s^-1, where the lunar escape speed is v_esc = 2.38 km s^-1. Then the predicted is recalculated to be 8.92−8.74 km s^-1, only 2.5%−0.7% smaller than the simulated values. The remaining differences may be due to the ignored acceleration from Earth in addition to the statistical errors. It should be pointed out that the unusually low impact speeds result from the fixed initial perihelion distances of projectiles q_p = 1 AU. For the whole population of the MBAs without this constraint as well as the NEOs, v_esc is even more negligible, so that the exclusion of the acceleration in the analytical model is quite acceptable.

Unlike the precise prediction of , the simulated are all obviously smaller than the analytical prediction of 51.8° (Eq. (82)). The difference between and π/ 4 is 13.7%−16.4% for the five cases. The reason probably is that our analytical model is planar, describing the variation of Θ on the equator, while in the three-dimensional simulations, the impacts elsewhere on the lunar surface are all statistically more oblique without isotropic impactors, and thus the simulated are diminished. Since theoretically is constant regardless of how severe the leading/trailing asymmetry is (Sect. 2.6), and the pole/equtor asymmetry resulting from the concentration of low-inclination impactors should have no dependence on v_M, the simulated were expected to be invariant with a_M. In fact, we do find of all the cases lying in a narrow range 41.1−42.6° with a fluctuation of only 2%.

and show a mild dependence on a_M, because their variations with a_M are combinations of those of and ( and $\mathit{D̅}\mathrm{\propto }\mathit{V̅}\begin{array}{c}{\mathit{\gamma }}_{\mathrm{c}}\\ \mathrm{\perp }\end{array}$). Owing to the fluctuation of the variation, their negative correlations with a_M are slightly obscured. The dependence on also causes the simulated , 5.81−6.03 km s^-1, to be lower by 9.3%−12.2% than the predicted values of 6.61−6.75 km s^-1 (Eq. (83)); this is also true for the simulated , 96.8−98.6 km, which are lower by 6.7%−8.3% than the predictions of 105.4−106.5 km (Eq. (84)). However, the difference between the simulated and the product of the simulated and is no more than 0.3% for all the cases, and the simulated is lower than ${\mathit{c}}_{\mathrm{c}}\mathit{d̅}\begin{array}{c}{\mathit{\gamma }}_{\mathrm{c}}^{\mathrm{\prime }}\\ \mathrm{p}\end{array}\mathit{V̅}\begin{array}{c}{\mathit{\gamma }}_{\mathrm{c}}\\ \mathrm{\perp }\end{array}$ for only 2.6% at most. This indicates that although the pole/equator asymmetry in the Θ distribution leads to an overestimation of by the analytical model, the relationships between the global averages derived in Sect. 2.6 are still well applicable.