Issue 
A&A
Volume 594, October 2016



Article Number  A52  
Number of page(s)  25  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201628598  
Published online  12 October 2016 
Analytical formulation of lunar cratering asymmetries
School of Astronomy and Space Science and Key Laboratory of Modern Astronomy and Astrophysics in Ministry of Education, Nanjing University, 210046 Nanjing, PR China
email: zhoujl@nju.edu.cn
Received: 28 March 2016
Accepted: 3 July 2016
Context. The cratering asymmetry of a bombarded satellite is related to both its orbit and impactors. The inner solar system impactor populations, that is, the mainbelt asteroids (MBAs) and the nearEarth objects (NEOs), have dominated during the late heavy bombardment (LHB) and ever since, respectively.
Aims. We formulate the lunar cratering distribution and verify the cratering asymmetries generated by the MBAs as well as the NEOs.
Methods. Based on a planar model that excludes the terrestrial and lunar gravitations on the impactors and assuming the impactor encounter speed with Earth v_{enc} is higher than the lunar orbital speed v_{M}, we rigorously integrated the lunar cratering distribution, and derived its approximation to the first order of v_{M}/v_{enc}. Numerical simulations of lunar bombardment by the MBAs during the LHB were performed with an Earth–Moon distance a_{M} = 20−60 Earth radii in five cases.
Results. The analytical model directly proves the existence of a leading/trailing asymmetry and the absence of near/far asymmetry. The approximate form of the leading/trailing asymmetry is (1 + A_{1}cosβ), which decreases as the apex distance β increases. The numerical simulations show evidence of a pole/equator asymmetry as well as the leading/trailing asymmetry, and the former is empirically described as (1 + A_{2}cos2ϕ), which decreases as the latitude modulus  ϕ  increases. The amplitudes A_{1,2} are reliable measurements of asymmetries. Our analysis explicitly indicates the quantitative relations between cratering distribution and bombardment conditions (impactor properties and the lunar orbital status) like A_{1} ∝ v_{M}/v_{enc}, resulting in a method for reproducing the bombardment conditions through measuring the asymmetry. Mutual confirmation between analytical model and numerical simulations is found in terms of the cratering distribution and its variation with a_{M}. Estimates of A_{1} for crater density distributions generated by the MBAs and the NEOs are 0.101−0.159 and 0.117, respectively.
Key words: planets and satellites: surfaces / Moon / minor planets, asteroids: general / methods: analytical
© ESO, 2016
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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 socalled 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 lowlatitude impacts compared to highlatitude regions resulting from the concentration of lowinclination 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 mainbelt asteroids (MBAs), which dominated during the late heavy bombardment (LHB) ~ 3.9 Gya, and the nearEarth 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 semimajor 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 powerlaw 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 Subkm Main Belt Asteroid Survey (SMBAS), which each found 861 MBAs in Rband and 1001 MBAs in both B and Rbands, 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 Nbody 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 semianalytically 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 thermaldynamical 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 twomode 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 semiempirical description of the crater density (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 semianalytical 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 (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 spinorbit resonance can be estimated by (3)where τ_{tide} is the tidal circularization timescale, m_{⊕} is the Earth mass, n_{M}, a_{M}, R_{M}, m_{M}, and are the mean motion, semimajor axis, radius, mass, and effective tidal dissipation factor of the Moon. A measure of this timescale is τ_{syn} ≲ 10^{2} 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 preencounter stage, and the impactors are all massless particles with semimajor 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 where G is the gravitational constant, M_{⊙} is the solar mass, f_{p} and f_{⊕} are true anomalies, and the terrestrial semimajor 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, (6)Thus, the encounter velocity depends on a_{p} and e_{p} (or q_{p}): (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 (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 a_{p} = 2.0 AU and q_{p} = 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 v_{p} (black solid arrow) decreases as q_{p} increases, and thus v_{enc} (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 xaxis 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 yaxis 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: Its magnitude, that is, the impact speed, is (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 v_{enc} (km s^{1}) as a function of a_{p} and q_{p}. 
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 a_{M} = 60 R_{⊕}, with v_{M} = 1.0 km s^{1}. The impactors, distributed extensively enough to cover the lunar orbit (black circle), are in the common direction of v_{enc} (black dashed arrow). Where f_{M} = 0°,45°,...,315°, the lunar velocity v_{M} (black solid arrow) and the impact velocity v for v_{enc} = 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 yaxis points to the apex and the minus xaxis points to Earth. In this frame, 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 (16)and its cross section is (17)Denoting with e = (−cosλ,−sinλ) the normal vector, the incidence angle on the unit area can be derived with (18)and the normal speed v_{⊥}, the normal component of v is (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 (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 (21)The bombarded hemisphere is then bounded by the longitude interval [ λ_{l},λ_{u} ], where the lower and upper limits λ_{l,u} = λ_{⊥} ∓ 90°, and (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 (23)It turns out
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 v_{enc} = 10, 20, 30, and 40 km s^{1}) and their global averages (dashed lines in the same color for the same v_{enc}, except for those in panel b), which are all black to represent the invariant value) are calculated with v_{M} = 1.0 km s^{1}, d_{min} = 0.5 km, α_{p} = 1.75, ρ = 3.4 × 10^{25} km^{2}, t = 3.7 Gyr, and d_{c} = 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 (26)Integration of dN over the interval [ f_{l},f_{u} ] leads to (27)To simplify the expression, we define σ ∈ (0,π) by (28)which leads to Therefore, (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, , , and . The integration leads to where Functions F(φ  k^{2}) and F(k^{2}) are the incomplete and complete elliptic integrals of the first kind; E(φ  k^{2}) and E(k^{2}) 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 The constants c_{p,c}, γ_{p,c}, and given by different studies vary, but their values are not important to the analytical derivation here. Still, they are needed for illustration and postprocessing of numerical simulations, and we therefore adopt the forms given by Le Feuvre & Wieczorek (2011) for a nonporous regime (working for d_{c}> 20 km) with the Moon as the target: where the nonporous scaling parameters are K = 1.17, ν_{1} = 0.22, and ν_{2} = 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 , independent of its orbital distribution. Given a minimum of projectile size d_{min}, a normalized size distribution is (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 (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 , 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 (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 sizefrequency 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, (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 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 (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 on any lunar longitude, and thus the slopes of the size distributions of impactors and craters are related by (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_{M}t/ (2π): (58)This form describes the impact density after a bombardment duration t, but its relative variation does not differ at all from the oneperiodic 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°, Substituting λ with β, the integrated variables as functions of β are where 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 (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 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, (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_{M}t/ (2π): (80)
Apparently, the global averages of V, Θ, and V_{⊥} are the integrals ^{!}vd^{2}C, ^{!}sinθd^{2}C, and over the above λ and f_{M} intervals divided by C, respectively. The integration results in
where is defined by , 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 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 (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
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 v_{enc} = 10, 20, 30, and 40 km s^{1}) are calculated using exact formulations with v_{M} = 1.0 km s^{1}, d_{min} = 0.5 km, α_{p} = 1.75, ρ = 3.4 × 10^{25} km^{2}, t = 3.7 Gyr, and d_{c} = 25 km. The approximate variations (dashed curves in the same color as the solid curves for the same v_{enc}) 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 v_{enc}) 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 β, (98)where A_{0}> 0 and 0 <A_{1}< 1. This function is symmetric about the point (90°,A_{0}) and has one maximum A_{0}(1 + A_{1}) at the apex (β = 0°) and one minimum A_{0}(1−A_{1}) 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_{0} and A_{1} entirely determine the variation of Γ. It is true that for every Γ, A_{0} is exactly the approximate and also the value of Γ on the prime meridian (β = 90°). Additionally, A_{0} 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_{1}cosβ) is then equal to the normalized distribution ΔΓ. The second parameter A_{1} 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_{1} ∝ η, 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_{1}.
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 These are the results of the relations between the variables, and the first three equations involving A_{0}, 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_{1}, the asymmetry amplitudes, are degenerate relations between A_{0}: 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 is independent of this, and only involves α_{p}. That 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 on 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 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 N_{c} lead to (red, orange, yellow, green, cyan, and blue curves) are compared. Where v_{enc} = 10, 20, 30, and 40 km s^{1} on condition v_{M} = 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 leastsquares 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_{1}, η can be easily reproduced based on the relations between A_{1} 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, where the lefthand sides are fits of and , 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}, , (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 km, , , and (with uncertainties ≲0.1%), among which and 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 bestfit 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 coupledasymmetric 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^{5} 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 semimajor 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 semimajor 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^{12} 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 Nbody simulations. The total mass of the particles of 1 × 10^{18} 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 postprocessing 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 powerlaw breaks. In this work, we assumed the normalized size distribution (Eq. (52)) to be a single powerlaw 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.
Threedimensional Nbody simulations were run to integrate the orbits of the Sun, Earth, the Moon, and 3 × 10^{5} 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 RungeKuttaFehlberg 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 highlatitude 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^{7} 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 a_{M} = 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 semimajor 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 highlatitude 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, . This is how the size distribution of impactors is involved in postprocessing. 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 lowlatitude region Γ_{low} are always higher than those of highlatitude 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 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 threedimensional 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 lowinclination 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 ). 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 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.
In general, simulated shows a negative correlation with a_{M} as well. (We note that the relative variation of also represents that of C_{gl}.) The approximate expression of (Eq. (93)) indicates, however, that it is proportional to the constant v_{enc}, while the exact expression involving v_{M} (Eq. (80)) is consistent with this negative correlation. Of the five sets of data points, case 2 (a_{M} = 30 R_{⊕}) seems abnormal, where is oddly deficient and N_{near} is clearly smaller than N_{far}. Assuming that the reason is that the Earth’s shielding reduces the impacts on the near side and thus the global impact number, which is supported by the fact that N_{far} of this case is larger than case 3 (a_{M} = 40 R_{⊕}) as expected, case 1 with its smaller a_{M} should have suffered the same effect. However, in case 1, N_{near} is even larger than N_{far}. Therefore, we prefer to consider it as a possibility instead of confirming it as evidence of the near/far asymmetry. Because , it is natural to see that the variation of with a_{M} combines those of and . Despite of the wavy fluctuation, generally decreases as a_{M} increases. Additionally, is close to with excess of 8.6%−9.7% for all the cases, in agreement with Eq. (85).
In short, this preliminary statistics confirms the presence of the leading/trailing asymmetry in spatial distributions of all Γ and the pole/equator asymmetry in those of all but V. All but show negative correlations with a_{M} as indicated by the analytical model, which precisely predicts the value of all cases and mildly overestimates and thus and , while the correlations between global averages apply well.
3.3. Spatial distributions
The cratering distribution of case a_{M} = 20 R_{⊕} is shown in Fig. 10, and those of other cases have no qualitative differences. Every value of Γ on a given location is the average over a circular area centered on this location with a “smooth radius” of r_{s} = 35° (along a great circle), equivalent to 1061 km on the lunar surface. Smaller r_{s} can reveal more details, but because of the limited number of data points, we consider 35° as the proper choice for showing important features.
Fig. 10 Simulated spatial distribution of impact speed a); incidence angle b); normal speed c); crater diameter d); impact density e); and crater density f) for the case a_{M} = 20 R_{⊕}. The value at each location is the average over a circular area centered on this location with a radius of 1061 km, assuming km and d_{c} = 100 km. 
Fig. 11 Impact speed as a function of angular distance from the apex a); and the incidence angle b); normal speed c); crater diameter d); impact density e); and crater density f) as functions of longitude, all in terms of simulated global averages. a) Data from the simulated distribution of the impact speed (black squares), each of which is the average over the small circle where β is fixed, are fit with Eq. (111). The best fit (black solid curve) is generated with plotted . b)−f) Data from simulated distributions, each of which is the mean of the two values at location (λ, ± ϕ) with  ϕ  being 0°, 30°, 60°, and 90° (black, red, blue, and green squares), are fit with Eqs. (112)−(116). The best fits (solid curves in the same color as the squares for the same  ϕ ) are generated with plotted A_{0,1,2}. 
The distribution of V (Fig. 10a) is typical of a pure leading/trailing asymmetry: the maximum 10.0 km s^{1} is near the apex, while the minimum 7.8 km s^{1} is near the antapex; values on the longitudes of 0° and ± 180°, the dividing line between leading and trailing sides, are approximately equal to the global average 9.2 km s^{1}; there is an obvious trend for V to be monotonously decreasing from apex to antapex. The distribution of Θ (Fig. 10b), however, combines properties of both the leading/trailing and pole/equator asymmetries: the maximum 47.2° is still near the apex, while the minimum seems to slide from the antapex toward the poles along the longitude 90° and turns into a pair of valleys on both sides of the equator; the contours around the apex tend to stretch along the latitude lines. Distributions of V_{⊥}, D, N, and N_{c} have similar structures to Θ, but with relatively milder pole/equator asymmetries, which is shown by the closer proximity of their valleys. We describe in Sect. 3.5 how a spatial distribution evolves with increasing leading/trailing and pole/equator asymmetries and show that all the above distributions but that of V are productions of the coupled asymmetries. Again it is not found necessary to involve the near/far asymmetry to explain the variations, as these two hemispheres are roughly identical in all aspects.
As expected, it is found that although θ of all the impacts lie in an extended range covering 0−90° (Fig. 8), the range of Θ, the average of θ for every unit area, is quite limited (Fig. 10b). The analysis in Sect. 2.3 indicates that as the normal impact point moves on the lunar surface, on every unit area, impacts from normal to extremely oblique can all occur, and that the varying moving rate adds the normal impacts on the apex, resulting in the asymmetry. In other words, the enlargement near the apex is not caused by the general upward shift of θ, but is due to the increased relative fraction of great θ, which is now verified by simulations.
In addition to the nonuniform moving rate of the normal impact point, the leading/trailing asymmetry can be understood in this way: the lunar velocity adds to the impactors heading toward the Moon and subtracts from the impactors chasing it, so that N is higher near the apex; the subtraction of the lunar velocity from that of the impactors, how the impact velocities are derived, increases the velocity components along the apexantapex axis of the heading impactors, but decreases those of the chasing impactors, meaning that statistically, not only the magnitudes of the impact velocities are greater near the apex, but their directions are also closer to the apexantapex axis, so that V and Θ are also biased. The pole/equator asymmetry is triggered by the low relative inclinations between the lunar equator and projectile orbits in our simulations. This initialization limits the number of impacts arriving at the poles, especially the normal impacts, and thus decreases N and Θ there. Although near the poles there are fewer impacts and they are statistically more oblique, the magnitudes of the impact velocities are not influenced, so that the pole/equator asymmetry is absent in the V distribution. Since v_{⊥} and d_{c} depend on both v and θ, the distributions of V_{⊥} and D naturally introduce both asymmetries in themselves. Similarly, because N_{c} depends on V_{⊥} and N, whose increases provide both a greater probability of generating large craters and more craters in total, the N_{c} distribution also shows the coupled asymmetries.
To describe the pure leading/trailing asymmetry, we use the approximate formulations derived in Sect. 2.7, whose common form is Γ ∝ (1 + A_{1}cosβ), where in a threedimensional model the angular distance from the apex β is related to the longitude λ and latitude ϕ by (108)To describe the pure pole/equator asymmetry, we suggest an empirical function (109)with a monotonic decrease from the equator to the poles. A spatial distribution of the coupled asymmetries is thus assumed to be (110)It can be easily proved by integrating the righthand side over the spherical surface that still holds as Sect. 2.7 indicates. At the same time, A_{0} is the value of Γ at location (0°, ± 45°) or (180°, ± 45°). The parameters A_{1} and A_{2} are leading/trailing and pole/equator asymmetry amplitudes, respectively. They are measurements of the asymmetry degrees, and are assumed to be independent of each other. A detailed discussion of the coupled asymmetries can be found in Sect. 3.5. We note that because of r_{s}, the asymmetry amplitudes are reduced. The decrease factor of A_{1} is estimated to be sinr_{s}/r_{s}, meaning that the fits of A_{1} here are the products of this factor and A_{1} given in Sect. 2.7. For r_{s} = 35°, this decrease factor is 0.939.
We fit these formulations with the nonlinear leastsquares method to the simulated spatial distributions shown in Fig. 10: Figure 11 shows the agreement between the simulated and fit variations of the above variables, with the fits of A_{0,1,2} plotted. For a fixed latitude, the variables always decrease with increasing distance from longitude − 90° toward + 90°, with their variation amplitudes and vertical positions dependent on the latitude. Given fits of A_{0,1,2}, the above formulations can lead to fit spatial distributions whose configurations are determined by the relative magnitudes of A_{1} and A_{2} (Sect. 3.5).
3.4. Fit parameters and reproduced conditions
Fig. 12 Fit A_{0}a); A_{1}b); and A_{2}c) of spatial distributions of impact speed, incidence angle, normal speed, crater diameter, impact density, and crater density (squares connected by solid lines in red, orange, yellow, green, cyan, and blue) as functions of lunar semimajor axis. a) Every approximate global average A_{0} is in terms of the relevant simulated global average. , , and are compared to , , and (squares connected by dashed lines in yellow, green, and blue). b) and c) and are compared to and (squares connected by dashed lines in green and blue). 
Fig. 13 Speed ratios and lunar orbital speeds reproduced from fit leading/trailing asymmetry amplitudes of impact speed, incidence angle, normal speed, crater diameter, impact density, and crater density (squares connected by solid lines in red, orange, yellow, green, cyan, and blue) in comparison with predictions (black dashed curve). 
The qualitative properties of the cratering distribution of case a_{M} = 20 R_{⊕} shown above are common to any other cases. We fit the cratering distribution with Eqs. (111)−(116) for each simulation case. In this subsection, the fit parameters, their quantitative interrelations, and their variations with a_{M} are shown and compared with the analytical predictions. Then the bombardment conditions are reproduced.
Fits of A_{0} of each Γ are plotted in terms of simulated in Fig. 12a. Theoretically, regardless of the pole/equator asymmetry (Sect. 3.3). In Fig. 12a, the fit A_{0} are all slightly smaller than the simulated , but only by a few percents, which is partly due to the approximation process itself. The largest difference is no more than 7% between and , and the smallest is only 1% between and . Therefore, even though we did not rigorously prove that the pole/quator asymmetry can be formulated as Eqs. (109), (110) is clearly a good description of the coupledasymmetric distribution and the fit of A_{0} is an excellent reproduction of . The reproduction goodness was expected to be proportional to (inversely proportional to v_{M}) with a given impactor population, but this trend is not observed because it is obscured by the statistical fluctuation. The variation of with a_{M} has been discussed in Sect. 3.2. Additionally, the interrelations between A_{0} (Eqs. (99)−(101)) are verified regardless of the involved pole/equator asymmetry: for any simulation case, approximates to with a difference of ≲0.5%; is almost equivalent to with a difference of <3%; and is close to with a difference of ≲11%.
Our analytical model indicates that A_{1} ∝ v_{M}/v_{enc}, that is, with a given impactor population. As shown in Fig. 12b, the variation of fit with a_{M} reflects this relation very well: , 0.469 ± 0.008, 0.455 ± 0.013, 0.491 ± 0.004, and 0.493 ± 0.006 for case 1−5 in turn, which is almost a constant. Because the leading/trailing asymmetry in the Θ distribution is the faintest, the fit involves relatively great fluctuation, and its relation with a_{M} is not clearly seen. The distributions of V_{⊥} and D depend on those of both V and Θ, so that the fit and correlate mildly negatively with a_{M}. The wavy decrease of fit and the monotonic decrease of are also qualitatively consistent with the above theoretic relation. The interrelations between A_{1} (Eqs. (102) and (103)) also apply quite well. For all cases, the differences between and and those between and are ≲3%.
We have assumed that A_{2} is invariant with the lunar orbit (Sect. 3.3). As a result, no common dependence of fit A_{2} on a_{M} is seen in Fig. 12c. To determine whether the assumption is realistic or if the relation is too weak to show itself has to wait for further investigation in the future. However, referring to the interrelations between the parameters A_{0,1}, we find that similar rules are present for A_{2}: The comparisons between the left and right sides are also shown in Fig. 12c. We caution that Eqs. (117) and (118) are empirically derived and not guaranteed to be valid for other bombardment conditions, although they probably are.
Given the fit A_{1}, the speed ratio η and thus the lunar orbital speed v_{M} were easily reproduced according to Eqs. (87)−(92). Figure 13 shows the reproduced η and v_{M} in comparison with their predictions. The decrease factor was considered in the reproduction. Generally speaking, the reproductions are reasonably consistent with each other. Still, η reproduced from is obviously smaller than its prediction. Theoretically, , but the constant quotient of the fit and predicted η is about 0.53 (with the decrease factor considered), which is smaller than the coefficient π/ 4 derived in Sect. 2.7, so that the reproduced η is decreased. Although the reproduction from is meant to be lower (Sect. 2.8), it is not enough to explain the difference of 32%. The reason why the fit precisely follows but underestimates the predicted η and v_{M} may be that the leading/trailing asymmetry of the V distribution is decreased by the introduction of the pole/equator asymmetry in the threedimensional numerical model, so that the coefficient decreases from π/ 4 to somewhere lower. It may also be that the mean v_{enc} is greater than the predicted v_{enc} because of the gravitational acceleration, which is ignored in the analytical model, so that η is smaller than its prediction. Because the pole/equator asymmetry is not seen in the V distribution and substituting the predicted v_{enc} with the simulated only reduces the reproduction error to about 27%, other places where the numerical simulations deviate from the analytical model need to be checked. Although statistical fluctuation results in a great error of η reproduced from , fortunately the reproductions from asymmetry amplitudes of V_{⊥}, D, N, and N_{c} are much better, especially for N, which leads to reproduction errors of 1%−15% for all cases.
Again we followed the method of reproducing bombardment conditions described in Sect. 2.8 and present a simple estimation based on the fit of case 1, for only D and N_{c} are observables in surveys of crater records. Assuming v_{enc} = 8.32 km s^{1} and t = 10^{7} yr are known, the reproduced conditions are calculated to be v_{M} = 1.45 ± 0.02 km s^{1}, α_{p} = 3.23 ± 0.06, d_{min} = 4.71 ± 0.03 km, and F = (1.821 ± 0.164) × 10^{4} yr^{1}. Comparisons of the former three to the predictions, v_{M} = 1.78 km s^{1}, α_{p} = 3, and d_{min} = 5 km, lead to differences of 19%, 8%, and 6%, respectively. We note that the impactor flux F is derived with Eq. (105) where l_{gl} = 2πR_{M} is substituted by , as the numerical simulations are threedimensional. It is equivalent to a total number of impacts of 1821 ± 164 during an integration time of 10^{7} yr, which is very close to the actual result C = 1544. Taking the limited number of impacts in simulations into account, it is rather inspiring to find that this method of reproducing the bombardment conditions from the cratering distribution works well. The importance of this method is fully shown when we recall that very many lunar craters were generated by MBAs during the LHB, meaning that our current knowledge of the LHB and MBAs may provide clues to the early lunar dynamical history.
3.5. Coupling of leading/trailing and pole/equator asymmetries
Figure 10 illustrates the spatial distributions of Γ. All of the distributions involve the leading/trailing asymmetry, and all but the V distribution also show the pole/equator asymmetry. The coupling of the two asymmetries can be described using Eq. (110), whose normalized form is (119)The asymmetry amplitudes A_{1,2} must be in [ 0,1 ] to ensure that ΔΓ is nonnegative and its first and second factors are in negative correlations with β and  ϕ . Clearly, the coupled distribution is symmetric about both the equator and the apexantapex axis, with a maximum (1 + A_{1})(1 + A_{2}) always at the apex. If ϕ is constant, the function degenerates to , where and . This means that on every latitude line, ΔΓ is maximized and minimized at longitude − 90° and + 90° with the average and the variation amplitude dependent on  ϕ . If λ is constant, the function is a cubic polynomial of cosϕ, whose first factor monotonically decreases or increases as  ϕ  increases on the leading or trailing side, respectively, and the second factor always decreases with increasing  ϕ . Therefore, with a given λ ∈ (−180°,0°), ΔΓ is always maximized on the equator and minimized on the poles, while with a given λ ∈ (0°,180°), the conditions vary with varying A_{1,2}. The distributions are thus classified into three configurations according to the numbers and locations of the extrema.
Fig. 14 Normalized cratering distribution of the coupled asymmetries. a) Pure leading/trailing asymmetry with A_{2} = 0; b) pure pole/equator asymmetry with A_{1} = 0; c) configuration 1 with and ; d) transition between configurations 1 and 2 with and ; e) configuration 2 with and ; f) transition between configurations 1 and 3 with and ; g) configuration 3 with and ; h) transition between configurations 3 and 2 with and . 
We summarize the theory distribution of the coupled asymmetries as follows. The case A_{1} ≠ 0 and A_{2} = 0 corresponds to the pure leading/trailing asymmetry, when the contour lines on the lunar surface are all circles centered on the apex and the antapex, which are the locations of the only maximum and minimum (Fig. 14a). The case A_{1} = 0 and A_{2} ≠ 0 corresponds to the pure pole/equator asymmetry, leaving a distribution consisting of latitudeparallel strips with a peak at the whole equator and two minima at the poles (Fig. 14b). Excluding those two limiting conditions, there are two critical values of A_{2}, and , and two critical values of A_{1} with A_{2} given, and , which together define the configurations. We note that always holds and only when .

1.
When is fixed, two configurations appear in turn as A_{1} increases. When , the contours begin twisting from the strips (Fig. 14c). The maximum can only be reached at the apex (− 90°,0°), while the minima slide symmetrically along the longitude 90° from the poles toward the antapex (90°,0°). The latitudes of the two minima are determined by (120)The antapex acts as a saddle point. We call this distribution configuration 1. The moment A_{1} increases to is when the minima join at the antapex and become one (Fig. 14d). After that when , there is only one maximum (1 + A_{1})(1 + A_{2}) at the apex and one minimum (1−A_{1})(1 + A_{2}) at the antapex (Fig. 14e). The contours near the apex and antapex resemble lying and standing ellipses. The higher the value of A_{1}, the closer to circles the contours. The distribution at this stage is called configuration 2.

2.
When is fixed, another configuration presents itself as well as the above two. When , configuration 1 is also observed, even though at the moment the two minima are still approaching each other. At this moment, the saddle at the antapex becomes two saddles that start approaching the minima, and the antapex itself becomes a third minimum (Fig. 14f). As increases, the saddle and minimum on each side of the equator come increasingly closer, while the antapex valley continuously deepens. This stage corresponds to configuration 3 (Fig. 14g), when the latitudes of the saddles are determined by (121)The moment is when each pair of approaching saddle and minimum combines and vanishes at the points (Fig. 14h). When , the distribution is that of configuration 2 based on our simple classification. However, it is necessary to point out that near the vanishing points, the contour line shapes are different from those with .

3.
When is fixed, the last stage of case belonging to configuration 2 will not appear. The reason is that only if can the critical value . When , A_{1} can never reach before increasing to 1.
Simply speaking, configurations 1, 2, and 3 are defined to have two, one, and three minima each and one maximum always. Our fit distributions of Θ, V_{⊥}, D, and N all belong to configuration 1 and 2 since none of their A_{2} is larger than in any simulation case. In particular, the Θ and D distributions of all the cases belong to configuration 1. The N_{c} distribution, with always larger than , belongs to configuration 2 for cases 1−3, configuration 1 for case 4 and configuration 3 for case 5. We note that we classify those distributions according to their fit A_{1} and A_{2}, but some simulated distributions are not quite symmetric about the equator, and therefore the pole/equator asymmetry amplitudes are sometimes underestimated by their fits. Nearly all the simulated distributions clearly exhibit the main properties of configuration 1 or 2 and are well fit. However, we cannot visually verify whether the N_{c} distribution of case 5 has a third minimum at the antapex, resulting form the small number statistics. More detailed studies in terms of both theory and observation are necessary to check the description of coupledasymmetric distribution.
4. Discussion
4.1. Comparison between NEOs and MBAs
Fig. 15 Contour map of the leading/trailing asymmetry amplitude of the crater density as a function of speed ratio η and sizedistribution slope α_{p}. The estimates of for NEOs and MBAs are indicated by dotted and dashed squares, respectively. 
Strom et al. (2005, 2015) suggested that the MBAs and NEOs that dominated during the LHB and since about 3.8−3.7 Gya in turn had formed craters on the Moon and even all the terrestrial planets. We review the current knowledge of two impactor populations in terms of orbital distribution, size distribution, etc. in Sect. 1. Now that we have formulated the cratering asymmetry based on our analytical model in Sect. 2, which has been verified in various aspects by numerical simulations in Sect. 3, we analytically estimate and compare the leading/trailing asymmetry amplitudes generated by the two impactor populations in this subsection.
Given a population of impactors, the typical encounter speed v_{enc} is determined by its orbital distribution, the minimum size d_{min} and slope α_{p} are determined by its size distribution, and the flux F is determined by both. The first three parameters influence the spatial distribution of D, and all of the four are involved in the N_{c} distribution. When the normalized distribution depending on just A_{1} is considered, then only v_{enc} and α_{p} matter.
Applying our analytical model, the inclination is not considered, and because we lack complete information about the orbital distribution of MBAs 3.9 Gya, we only take the difference between MBAs and NEOs in semimajor axis a_{p} into account, that is, the distributions of perihelion distance q_{p} are assumed to be the same. The semimajor axes of MBAs lie in the range of 2.0−3.5 AU, while those of NEOs lie in the range of 0.5−2.8 AU according to the “constrained target region” of the debiased NEO model built by Bottke et al. (2002). When we adopt the medians of a_{p}, 2.75 AU for MBAs and 1.65 AU for NEOs, and a median of q_{p}, 0.5 AU, the typical encounter speeds of MBAs and NEOs are calculated to be 25.4 km s^{1} and 22.1 km s^{1} (Eq. (7)), and thus their typical impact speeds are 25.5 km s^{1} and 22.3 km s^{1} (). This estimate of v_{imp} for NEOs is consistent with results of Gallant et al. (2009), Ito & Malhotra (2010) and Le Feuvre & Wieczorek (2011), who found the mean impact speed to be 20 km s^{1}, 22.4 km s^{1}, and 19.7 km s^{1}, respectively. The lunar orbital speeds v_{M} during the dominant epochs of the two impactor populations are different as well. For MBAs, v_{M} during the LHB may be from 1.12 to 1.26 km s^{1}, relevant to the uncertain a_{M} from 50 to 40 R_{⊕} (Sect. 1). For NEOs, the current value v_{M} = 1.03 km s^{1} is used since a_{M} is believed to have been beyond 40 R_{⊕} for most of the Moon lifetime. Here e_{M} is ignored in the calculation.
Therefore, the speed ratio η is 0.044−0.049 for MBAs and 0.046 for NEOs. A clear excess or deficiency of η of MBAs compared to NEOs is not seen. This is also true for the asymmetry amplitude because it is proportional to η (Eq. (106)): is 0.0205−0.0229 for MBAs and 0.0215 for NEOs. Conversely, the fit of an observed area on the lunar surface that is purely cratered by MBAs can lead to a speed ratio and thus an Earth–Moon distance during the LHB, helping to clarify the early lunar history. On the other hand, when of MBAs is precisely determined in theory, the fit of a cratered area can lead to the realistic ratio of crater numbers of the two impactor populations.
The asymmetry amplitude depends on not only η, but also on α_{p} (Eq. (107)). Assuming the size distribution slope of the NEOs is α_{p} = 1.75 (Bottke et al. 2002), Eq. (57) leads to the slope of craters generated by them of α_{c} = 2.07, which is well consistent with Strom et al. (2005, 2015), who found that the size distribution of craters on lunar young plains has a single slope of 2. Thus the of NEOs is calculated to be 0.117. The size distribution of the MBAs is not a single power law as mentioned in Sect. 1, and the variation of its slope is not slight, as Ivezić et al. (2001) suggested for d_{p} = 0.4−5 km and 5−40 km, α_{p} = 1.3 and 3, respectively. Still, all of our analysis is valid for the separated part of impactors in a given d_{p} interval with invariant α_{p} and the craters they generate. Therefore, the possible range of of MBAs shown in Fig. 15 is defined by the topright maximum 0.159, where η = 0.049 (a_{M} = 40 R_{⊕}) and α_{p} = 3, and the bottom left minimum 0.101, where η = 0.044 (a_{M} = 50 R_{⊕}) and α_{p} = 1.3. For the considered part of MBAs with d_{p} = 0.4−40 km, Fig. 15 shows that their is probably greater than that of the NEOs. The N_{c} distribution generated by impactors with a size distribution of multiple power laws will be formulated in a future work.
We caution that even the size distribution is not a single power law, and thus , the definition (Eq. (53)) still holds, so that the formulations of the D distribution and also apply. The approximate (Eq. (96)) only depends on and v_{enc}. Because v_{enc} of the MBAs is greater than that of NEOs and MBAs are generally larger than NEOs, this results in the separation into two crater populations in terms of size. A rough estimation leads to 33 km for MBAs and 4 km for NEOs, assuming 1 km and 0.1 km, respectively.
4.2. Apex/antapex and pole/equator ratios
Fig. 16 Reproduction of spatial distribution of the crater density in Le Feuvre & Wieczorek (2011). It is generated using Eq. (119) with asymmetry amplitudes A_{1} = 0.156 and A_{2} = 0.093. 
Fig. 17 Apex/antapex ratio of the crater density generated by the current impactors in related works compared with ours. Our relation between N_{c}(apex) /N_{c}(antapex) and a_{M} with v_{enc} = 22.1 km s^{1} and α_{p} = 1.75 (solid curve) is close to the empirical one in Le Feuvre & Wieczorek (2011), which works for a_{M} = 20−60 R_{⊕} (Eq. (2); dashed curve); the former with v_{enc} = 19.7 km s^{1} and α_{p} = 2.59 (dotted curve) is even nearly equal to the latter in its working range. Simulation results from Gallant et al. (2009) and Ito & Malhotra (2010) assuming the current lunar orbit (triangle and square) are also indicated. 
Earlier studies on cratering have often used the apex/antapex and pole/equator ratios to indicate the asymmetry degree of the spatial distribution. The former is the ratio of the value (typically the crater density N_{c}) at the apex to the antapex, which increases with increasing leading/trailing asymmetry, while the latter is the ratio of the value at the poles to the equator, which increases with decreasing pole/equator asymmetry. These ratios can be obtained easily and are especially convenient for simulations and observations, but their physical meanings were not understood before. Neither was it possible to completely define a spatial distribution through them. According to the formulation of coupled asymmetries that we established based on our analytical model and simulation results (Eq. (110)), the ratios are each related to the relevant asymmetry amplitude by (When A_{1} ≠ 0, Γ(equator) means the average over the equator.) Because a whole distribution can be directly reproduced given A_{1} and A_{2} using Eq. (119) and because the meaning of A_{1} has explicitly been determined in Sect. 2.7, the asymmetry amplitudes are the better measurements of the asymmetry degrees for theories.
We now reproduce the spatial distribution of N_{c} shown in Fig. 5 of Le Feuvre & Wieczorek (2011) through asymmetry amplitudes. They obtained it by assuming the current Earth–Moon distance and the current asteroids and comets in the inner solar system as impactors. The value of N_{c} varies from 0.80 at (90°, ± 65°) to 1.25 at the apex in terms of the global average, which leads to a maximumminimum ratio of 1.5, an apex/antapex ratio of 1.37, and a pole/equator ratio of 0.80. Using their apex/antapex ratio and minimum location, we reproduce the distribution shown in Fig. 16. It is derived from Eq. (119) with the asymmetry amplitudes = 0.156 and = 0.093, which are calculated with Eqs. (122) and (120) to ensure N_{c}(apex) /N_{c}(antapex) = 1.37 and the latitude of the minimum ϕ_{min} = ± 65°. Then its maximum and minima are 1.26 at apex and 0.88 at (90°, ± 65°), whose ratio is 1.44, well consistent with that of Le Feuvre & Wieczorek (2011), and our pole/equator ratio 0.83 is nearly equal to theirs. Not only the N_{c} distribution is well reproduced, but also the estimate α_{p} = 2.6 is derived from Eq. (107) with the reproduced , current v_{M}, and the mean impact speed on the Moon v_{enc} = 19.7 km s^{1} according to Le Feuvre & Wieczorek (2011). Therefore, a powerlaw size distribution with this slope is found to be a convenient substitute for the 10thorder polynomial that Le Feuvre & Wieczorek (2011) used for the empirical fit in their approach.
We have determined the leading/trailing asymmetry amplitude of the N_{c} distribution generated by NEOs to be = 0.117 (Sect. 4.1), which is equivalent to N_{c}(apex) /N_{c}(antapex) = 1.27. As reviewed in Sect. 1, Morota & Furumoto (2003) reported the apex/antapex ratio to be ~1.5 based on observations of young craters on the lunar highlands, Gallant et al. (2009) and Ito & Malhotra (2010) found the ratio to be 1.28 ± 0.01 and 1.32 ± 0.01 from their numerical simulations, and Le Feuvre & Wieczorek (2011) calculated it to be 1.37 using their semianalytical approach, all in good consistence with our estimate. Moreover, Le Feuvre & Wieczorek (2011) also empirically fit the dependence of N_{c}(apex) /N_{c}(antapex) on a_{M} and obtained Eq. (2) that is valid for a_{M} = 20−60 R_{⊕}. Treating the expression of (Eq. (107)) as a function of a_{M}, that is, where , we can also derive a relation (124)The two relations are compared in Fig. 17 with solid and dashed lines. We note that in calculating our relation, v_{enc} = 22.1 km s^{1} and α_{p} = 1.75 are the same as used in Sect. 4.1, but a_{M} is no longer fixed at 60 R_{⊕} since it is now considered a variable. Figure 17 shows that the two relations are quite close especially where that of Le Feuvre & Wieczorek (2011) is valid. When we substitute their mean impact speed on the Moon of 19.7 km s^{1} for our roughly estimated v_{enc} and adopt α_{p} = 2.6 as reproduced from their N_{c} distribution, our function is identical to their valid part with a difference of only <2%.
We briefly point out the difference between our work and that of Le Feuvre & Wieczorek (2011) in cratering modeling. Their model is able to produce the pole/equator asymmetry since the anisotropic impactor population was adopted, while our current analytical model alone can only formulate the leading/trailing asymmetry, but with our empirical description of the pole/equator asymmetry based on numerical simulations, we can also describe, reproduce, and measure the coupled asymmetries. More importantly, they integrated the cratering distribution numerically, while we give explicit analytical functions after rigorous deduction, which clearly shows the influence of every factor and makes the formulations very easy to apply.
4.3. Generalization of the analytical model
In addition to the Earth–Moon system, our analytical model (Sect. 2) can be applied to other planetsatellite systems. In the inner solar system, the NEOs and MBAs are the common impactor populations dominating different epochs. They encounter Mars with v_{enc} ≈ 19.3 and 22.9 km s^{1}, respectively, which is calculated with a_{⊕} replaced by the semimajor axis of Mars a_{Mars} = 1.5 AU in Eq. (7) and the typical a_{p} of NEOs being the median of the interval 0.76−2.8 AU instead of 0.5−2.8 AU because projectiles with a_{p}<a_{Mars}/ 2 on elliptic orbits cannot encounter Mars. The Martian moons, Phobos and Deimos, which lie on very circular orbits, have current orbital speeds of v_{orb} = 2.14 and 1.35 km s^{1} and negligible escape speeds v_{esc}. Thus, the NEOs with a typical speed ratio η of 0.111 for Phobos and 0.070 for Deimos are expected to generate a current leading/trailing asymmetry amplitude of the crater density of 0.28 and 0.18, respectively, both of which are much greater than for the Moon. Similarly, the orbital speeds of Phobos and Deimos and even their orbits during the LHB might be derived from the observed crater records if they had been tidallocked then.
The dominant impactors in the outer solar system are ecliptic comets, which are generally considered to originate in the Kuiper belt (Zahnle et al. 2001). Because their encounter speeds are generally lower than the orbital speeds of the satellites of the giants, our analytical model is not completely applicable. The divergence starts from the behavior of the normal impact point. With our assumption v_{enc}>v_{orb}, the normal impact point moving westward along the satellite’s equator (coplanar with its orbit and ecliptic) reaches the antapex and apex when the satellite is at its pericenter (where v_{enc} is in the same direction as v_{orb}) and apocenter, while the impact speed v is minimized and maximized (Sect. 2.3). However, if v_{enc}<v_{orb}, the normal impact point reaches the apex at both the pericenter and apocenter. Specifically, as the satellite’s true anomaly f increases from 0° (pericenter), the normal impact point leaves the apex, moves eastward until f = arccos(v_{enc}/v_{orb}) ∈ (0°,90°), turns around and returns to the apex when f = 180°; as f continues increasing, it continues moving westward until f = 360°−arccos(v_{enc}/v_{orb}) ∈ (270°,360°) and then returns to the apex when f = 360°. Its path is symmetric about the apex, bounded by the two return points at equal distances of β_{max} = arctan(v_{enc}/ from it. It is seen that β_{max}< 90°, so the satellite surface is divided into four areas: the area centered on the antapex with the longitude λ ∈ (β_{max},180°−β_{max}) never is hit by impactors; the area (β_{max}−180°,−β_{max}) centered on the apex is always under bombardment; and the two remaining symmetric areas on the near and far side each are periodically included in and excluded from the instant bombarded hemisphere. Apparently, the higher the ratio v_{orb}/v_{enc}, the smaller the distance β_{max}, and thus the greater the asymmetry. As a result, the leading/trailing asymmetry is much enhanced when v_{enc}<v_{orb} than otherwise. A generalized model including this case will be established in the future to examine the cratering of giants’ satellites.
Zahnle et al. (2001) suggested a semiempirical description of N_{c} distribution on satellites of giant planets shown as Eq. (1). We find qualitative consistence of our approximate formulation N_{c}(β) (Eq. (92)) with theirs since the asymmetry is amplified by speed ratio v_{orb}/v_{enc} and sizedistribution slope α_{p} in both forms. One difference is that α_{p} is a coefficient of cosβ in our form but acts as an exponent in theirs. Considering they derived their formulation based on simulation results, and referring to our analysis, there is a great possibility that α_{p} is still a coefficient (in approximate description) when v_{enc}<v_{orb} and the influence of α_{p} was overestimated by formulation of Zahnle et al. (2001), which will be examined by future model.
5. Conclusions
We analytically formulated the lunar cratering distribution and confirmed the derivation with numerical simulations. The formulations are quite easy to use. They are able to give results nearly identical to related works, avoiding the timeconsuming simulations and clarifying the physical meanings in quantitative relations.
Based on a planar model excluding the gravitations of Earth and the Moon on the impactors, we derived series of formulations of the cratering distribution on the synchronously rotating satellite through rigorous integration. The formulations directly and unambiguously proved the existence of a leading/trailing cratering asymmetry and the identity of the near and far sides in all aspects of cratering, that is, impact speed V, incidence angle Θ, normal speed V_{⊥}, crater diameter D, impact density N, and crater density N_{c}. Series expansion to the first order of v_{M}/v_{enc} resulted in the approximate Γ distributions in the common form , where is in positive correlation with v_{enc}, except that is the constant 51.8°, and A_{1} ∝ v_{M}/v_{enc}, that is, with a given v_{enc}. The relations between the cratering distribution and the bombardment conditions including properties of impactor population and the lunar orbit during the bombardment provides a viable method of reproducing the latter by measuring the former. The lower the speed ratio v_{M}/v_{enc}, the better the approximation to the exact Γ distributions and the better the reproduction of the bombardment conditions. In particular, the approximation of N_{c} can be improved by decreasing α_{p} as well.
We also numerically simulated the cratering on the Moon that is caused by the MBAs during the LHB with five simulation cases in which a_{M} = 20−60 R_{⊕}. Not only the leading/trailing asymmetry is present in the simulated cratering distribution in all aspects, but also the pole/equator asymmetry is seen, except for the V distribution, while the signs of the near/far asymmetry are not enough. The pole/equator asymmetry is empirically formulated with (1 + A_{2}cos2ϕ), leading to the description of coupled asymmetries , which is well fit to the simulated Γ distributions with errors of smaller than 6%. For each case, the analytical prediction of is only a few percent lower than the simulated , but the predicted is obviously larger than the simulated because the pole/equator asymmetry is involved. As predicted by the analytical model, except for all show negative correlations with a_{M} in general, and so do the fit A_{1} except for , whose variation with a_{M} involves greater statistical fluctuations, while the fit A_{2} seem independent of a_{M}. Additionally, the method of reproducing performs well, with reproduction errors being ~10% for bombardment conditions in case 1. This means that it is possible to speculate about the lunar orbital status during the LHB when the cratering distribution generated by MBAs is determined by observations.
Based on our analytical model, the leading/trailing asymmetry amplitude of the crater density that is generated by the MBAs is estimated to be , while that generated by the NEOs is , which is equivalent to N_{c}(apex) /N_{c}(antapex) = 1.27. We easily but precisely reproduced the N_{c} distribution and the variation of N_{c}(apex) /N_{c}(antapex) as a function of a_{M} derived by Le Feuvre & Wieczorek (2011), with α_{p} = 2.6, , and = 0.093. Our analytical model is applicable to other planetsatellite systems as long as v_{enc}>v_{orb} and will be generalized to the case v_{enc}<v_{orb} in future.
Acknowledgments
We thank the referee Bill Bottke for constructive and inspiring comments. This research has been supported by the Key Development Program of Basic Research of China (Nos. 2013CB834900), the National Natural Science Foundations of China (Nos. 11003010 and 11333002), the Strategic Priority Research Program “The Emergence of Cosmological Structures” of the Chinese Academy of Sciences (Grant No. XDB09000000), the Natural Science Foundation for the Youth of Jiangsu Province (No. BK20130547) and the 985 project of Nanjing University and Superiority Discipline Construction Project of Jiangsu Province.
References
 Bandermann, L. W., & Singer, S. F. 1973, Icarus, 19, 108 [NASA ADS] [CrossRef] [Google Scholar]
 Bills, B. G., & Ray, R. D. 1999, Geophys. Res. Lett., 26, 3045 [NASA ADS] [CrossRef] [Google Scholar]
 Bottke, W. F., Morbidelli, A., Jedicke, R., et al. 2002, Icarus, 156, 399 [NASA ADS] [CrossRef] [Google Scholar]
 Bottke, W. F., Durda, D. D., Nesvorný, D., et al. 2005, Icarus, 175, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Chambers, J. E. 1999, MNRAS, 304, 793 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Gallant, J., Gladman, B., & Ćuk, M. 2009, Icarus, 202, 371 [NASA ADS] [CrossRef] [Google Scholar]
 Gerstenkorn, H. 1955, Z. Astrophys., 36, 245 [NASA ADS] [Google Scholar]
 Goldreich, P. 1966, Rev. Geophys. Space Phys., 4, 411 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Gomes, R., Levison, H. F., Tsiganis, K., & Morbidelli, A. 2005, Nature, 435, 466 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Hansen, K. S. 1982, Rev. Geophys. Space Phys., 20, 457 [NASA ADS] [CrossRef] [Google Scholar]
 Ito, T., & Malhotra, R. 2010, A&A, 519, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ivezić, Ž., Tabachnik, S., Rafikov, R., et al. 2001, AJ, 122, 2749 [NASA ADS] [CrossRef] [Google Scholar]
 Kagan, B. 1997, Prog. Oceanogr., 40, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Kagan, B. A., & Maslova, N. B. 1994, Earth Moon and Planets, 66, 173 [NASA ADS] [CrossRef] [Google Scholar]
 Kawamura, T., Morota, T., Kobayashi, N., & Tanaka, S. 2011, Geophys. Res. Lett., 38, 15201 [NASA ADS] [CrossRef] [Google Scholar]
 Lambeck, K. 1977, Philosophical Transactions of the Royal Society of London A: Mathematical, Phys. Eng. Sci., 287, 545 [Google Scholar]
 Le Feuvre, M., & Wieczorek, M. A. 2005, in Lunar and Planetary Inst. Technical Report, eds. S. Mackwell, & E. Stansbery, 36 [Google Scholar]
 Le Feuvre, M., & Wieczorek, M. A. 2008, Icarus, 197, 291 [NASA ADS] [CrossRef] [Google Scholar]
 Le Feuvre, M., & Wieczorek, M. A. 2011, Icarus, 214, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Levison, H. F., Duncan, M. J., Zahnle, K., Holman, M., & Dones, L. 2000, Icarus, 143, 415 [NASA ADS] [CrossRef] [Google Scholar]
 MacDonald, G. J. F. 1964, Rev. Geophys. Space Phys., 2, 467 [Google Scholar]
 Morbidelli, A., & Gladman, B. 1998, Meteorit. Planet. Sci., 33, 999 [NASA ADS] [CrossRef] [Google Scholar]
 Morota, T., & Furumoto, M. 2003, Earth Planet. Sci. Lett., 206, 315 [NASA ADS] [CrossRef] [Google Scholar]
 Oberst, J., Christou, A., Suggs, R., et al. 2012, Planet. Space Sci., 74, 179 [NASA ADS] [CrossRef] [Google Scholar]
 Parker, A., Ivezić, Ž., Jurić, M., et al. 2008, Icarus, 198, 138 [NASA ADS] [CrossRef] [Google Scholar]
 Ross, M. N., & Schubert, G. 1989, J. Geophys. Res.: Solid Earth, 94, 9533 [NASA ADS] [CrossRef] [Google Scholar]
 Strom, R. G., Malhotra, R., Ito, T., Yoshida, F., & Kring, D. A. 2005, Science, 309, 1847 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Strom, R. G., Renu, M., Xiao, Z.Y., et al. 2015, RA&A, 15, 407 [Google Scholar]
 Touma, J., & Wisdom, J. 1994, AJ, 108, 1943 [NASA ADS] [CrossRef] [Google Scholar]
 Webb, D. J. 1980, Geophys. J., 61, 573 [NASA ADS] [CrossRef] [Google Scholar]
 Webb, D. J. 1982, Geophys. J., 70, 261 [NASA ADS] [CrossRef] [Google Scholar]
 Wiesel, W. 1971, Icarus, 15, 373 [NASA ADS] [CrossRef] [Google Scholar]
 Yoshida, F., & Nakamura, T. 2007, Planet. Space Sci., 55, 1113 [NASA ADS] [CrossRef] [Google Scholar]
 Yoshida, F., Nakamura, T., Watanabe, J.I., et al. 2003, PASJ, 55, 701 [NASA ADS] [CrossRef] [Google Scholar]
 Zahnle, K., Dones, L., & Levison, H. F. 1998, Icarus, 136, 202 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Zahnle, K., Schenk, P., Sobieszczyk, S., Dones, L., & Levison, H. F. 2001, Icarus, 153, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Zhou, J.L., & Lin, D. N. C. 2008, in IAU Symp. 249, eds. Y.S. Sun, S. FerrazMello, & J.L. Zhou, 285 [Google Scholar]
All Figures
Fig. 1 Encounter geometry for an impactor orbit with a_{p} = 2.0 AU and q_{p} = 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 v_{p} (black solid arrow) decreases as q_{p} increases, and thus v_{enc} (red solid arrow) shrinks. 

In the text 
Fig. 2 Contour map of the encounter speed v_{enc} (km s^{1}) as a function of a_{p} and q_{p}. 

In the text 
Fig. 3 Impact geometry seen in the rest frame of Earth a) and Moon b). a) The Moon is assumed to be where a_{M} = 60 R_{⊕}, with v_{M} = 1.0 km s^{1}. The impactors, distributed extensively enough to cover the lunar orbit (black circle), are in the common direction of v_{enc} (black dashed arrow). Where f_{M} = 0°,45°,...,315°, the lunar velocity v_{M} (black solid arrow) and the impact velocity v for v_{enc} = 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. 

In the text 
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. 

In the text 
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 v_{enc} = 10, 20, 30, and 40 km s^{1}) and their global averages (dashed lines in the same color for the same v_{enc}, except for those in panel b), which are all black to represent the invariant value) are calculated with v_{M} = 1.0 km s^{1}, d_{min} = 0.5 km, α_{p} = 1.75, ρ = 3.4 × 10^{25} km^{2}, t = 3.7 Gyr, and d_{c} = 25 km. 

In the text 
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 v_{enc} = 10, 20, 30, and 40 km s^{1}) are calculated using exact formulations with v_{M} = 1.0 km s^{1}, d_{min} = 0.5 km, α_{p} = 1.75, ρ = 3.4 × 10^{25} km^{2}, t = 3.7 Gyr, and d_{c} = 25 km. The approximate variations (dashed curves in the same color as the solid curves for the same v_{enc}) 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 v_{enc}) 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). 

In the text 
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 N_{c} lead to (red, orange, yellow, green, cyan, and blue curves) are compared. Where v_{enc} = 10, 20, 30, and 40 km s^{1} on condition v_{M} = 1.0 km s^{1} are indicated (vertical dot lines). 

In the text 
Fig. 8 Impact distribution on the lunar surface for the case a_{M} = 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. 

In the text 
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 semimajor 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 highlatitude regions. 

In the text 
Fig. 10 Simulated spatial distribution of impact speed a); incidence angle b); normal speed c); crater diameter d); impact density e); and crater density f) for the case a_{M} = 20 R_{⊕}. The value at each location is the average over a circular area centered on this location with a radius of 1061 km, assuming km and d_{c} = 100 km. 

In the text 
Fig. 11 Impact speed as a function of angular distance from the apex a); and the incidence angle b); normal speed c); crater diameter d); impact density e); and crater density f) as functions of longitude, all in terms of simulated global averages. a) Data from the simulated distribution of the impact speed (black squares), each of which is the average over the small circle where β is fixed, are fit with Eq. (111). The best fit (black solid curve) is generated with plotted . b)−f) Data from simulated distributions, each of which is the mean of the two values at location (λ, ± ϕ) with  ϕ  being 0°, 30°, 60°, and 90° (black, red, blue, and green squares), are fit with Eqs. (112)−(116). The best fits (solid curves in the same color as the squares for the same  ϕ ) are generated with plotted A_{0,1,2}. 

In the text 
Fig. 12 Fit A_{0}a); A_{1}b); and A_{2}c) of spatial distributions of impact speed, incidence angle, normal speed, crater diameter, impact density, and crater density (squares connected by solid lines in red, orange, yellow, green, cyan, and blue) as functions of lunar semimajor axis. a) Every approximate global average A_{0} is in terms of the relevant simulated global average. , , and are compared to , , and (squares connected by dashed lines in yellow, green, and blue). b) and c) and are compared to and (squares connected by dashed lines in green and blue). 

In the text 
Fig. 13 Speed ratios and lunar orbital speeds reproduced from fit leading/trailing asymmetry amplitudes of impact speed, incidence angle, normal speed, crater diameter, impact density, and crater density (squares connected by solid lines in red, orange, yellow, green, cyan, and blue) in comparison with predictions (black dashed curve). 

In the text 
Fig. 14 Normalized cratering distribution of the coupled asymmetries. a) Pure leading/trailing asymmetry with A_{2} = 0; b) pure pole/equator asymmetry with A_{1} = 0; c) configuration 1 with and ; d) transition between configurations 1 and 2 with and ; e) configuration 2 with and ; f) transition between configurations 1 and 3 with and ; g) configuration 3 with and ; h) transition between configurations 3 and 2 with and . 

In the text 
Fig. 15 Contour map of the leading/trailing asymmetry amplitude of the crater density as a function of speed ratio η and sizedistribution slope α_{p}. The estimates of for NEOs and MBAs are indicated by dotted and dashed squares, respectively. 

In the text 
Fig. 16 Reproduction of spatial distribution of the crater density in Le Feuvre & Wieczorek (2011). It is generated using Eq. (119) with asymmetry amplitudes A_{1} = 0.156 and A_{2} = 0.093. 

In the text 
Fig. 17 Apex/antapex ratio of the crater density generated by the current impactors in related works compared with ours. Our relation between N_{c}(apex) /N_{c}(antapex) and a_{M} with v_{enc} = 22.1 km s^{1} and α_{p} = 1.75 (solid curve) is close to the empirical one in Le Feuvre & Wieczorek (2011), which works for a_{M} = 20−60 R_{⊕} (Eq. (2); dashed curve); the former with v_{enc} = 19.7 km s^{1} and α_{p} = 2.59 (dotted curve) is even nearly equal to the latter in its working range. Simulation results from Gallant et al. (2009) and Ito & Malhotra (2010) assuming the current lunar orbit (triangle and square) are also indicated. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.