Analytical formulation of lunar cratering asymmetries

We formulate the lunar cratering distribution and verify the cratering asymmetries generated by the main-belt asteroids (MBAs) as well as the near-Earth objects (NEOs). Based on a planar model that excludes the terrestrial and lunar gravitations on the impactors and assuming the impactor encounter speed with Earth $v_{\rm{enc}}$ is higher than the lunar orbital speed $v_{\rm{M}}$, we rigorously integrated the lunar cratering distribution, and derived its approximation to the first order of $v_{\rm{M}}/v_{\rm{enc}}$. Numerical simulations of lunar bombardment by the MBAs during the late heavy bombardment were performed with an Earth-Moon distance $a_{\rm{M}}$ = 20--60 Earth radii in five cases. 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\beta)$, which decreases as the apex distance $\beta$ 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\varphi)$, which decreases as the latitude modulus $|\varphi|$ 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 \propto v_{\rm{M}}/v_{\rm{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_{\rm{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.


Introduction
The nonuniformity of the cratering distribution on a satellite or a planet when it comes under bombardment is called cratering asymmetry. There are three types. The first is the leading/trailing asymmetry, which is due to the synchronous rotation of a satellite. When synchronously locked, a satellite's velocity always points to its leading side, so that this hemisphere tends to gain a higher impact probability, higher impact speed and more normal impacts. Because the crater diameter depends on impact speed and incidence angle (e.g., Le Feuvre & Wieczorek 2011), this can also lead to an enhanced crater size on the leading side. This is the so-called apex/antapex effect (Zahnle et al. 2001;Le Feuvre & Wieczorek 2011). Zahnle et al. (1998Zahnle et al. ( , 2001 and Levison et al. (2000) investigated this effect on the giant planets and their satellites in detail. This work focuses on the Moon. Its leading/trailing asymmetry has been proposed by theoretical works (Zahnle et al. 1998;Le Feuvre & Wieczorek 2011), confirmed by numerical simulations (Gallant et al. 2009;Ito & Malhotra 2010), and directly verified by seismic observations (Kawamura et al. 2011;Oberst et al. 2012) and observations of rayed craters (Morota & Furumoto 2003). The second type, the pole/equator asymmetry (or latitudinal asymmetry), which is an enhancement of low-latitude impacts compared to high-latitude regions resulting from the concentration of low-inclination projectiles, has often been reported as well (Le Feuvre & Wieczorek 2008Gallant 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. (2005Strom et al. ( , 2015 suggested that there are two impactor populations in the inner solar system: the main-belt asteroids (MBAs), which dominated during the late heavy bombardment (LHB) ∼ 3.9 Gya, and the near-Earth objects (NEOs), which have dominated since about 3.8−3.7 Gya. The two populations are different in their orbital and size distributions, fluxes, and origins, which means that there is no reason to expect the cratering asymmetries generated by them to be the same.
We note the MBAs referred to here are the asteroids that occupied the region where the current main belt is during the LHB. Although their semi-major axes a p = 2.0−3.5 AU were the same as today, their eccentricities and inclinations were greatly excited by migrating giant planets (Gomes et al. 2005). Instead, their size distribution has changed little after the first ∼100 Myr (Bottke et al. 2005). Several surveys of the current main belt have reported power-law breaks of its size distribution. The Sloan Digital Sky Survey (SDSS; Ivezić et al. 2001) The NEO population is relatively well understood and suggested to have been in steady state for the past ∼3 Gyr, constantly resupplied mainly from the main belt (Bottke et al. 2002). Bottke et al. (2002) derived the debiased orbital and size distributions of NEOs by fitting a model population to the known NEOs: the orbits are constrained in a range a p = 0.5−2.8 AU, e p < 0.8, i p < 35 • , the perihelion distance q p < 1.3 AU, and the aphelion distance Q p > 0.983 AU; the size distribution is characterized by a single slope α p = 1.75 ± 0.1 for d p = 0.2−4 km. The corresponding slope for the crater size distribution indicated by Strom et al. (2015) is α c = 2 for d c = 0.02−100 km. A few recent works have investigated the cratering asymmetry generated by NEOs based on the debiased NEO model described in Bottke et al. (2002). Gallant et al. (2009) ran N-body simulations and determined an apex/antapex ratio (ratio of the crater density at the apex to antapex) of 1.28 ± 0.01, a polar deficiency of ∼10%, and the absence of a near/far asymmetry. In a similar work but with a different numerical model, Ito & Malhotra (2010) found a leading/trailing hemispherical ratio of 1.32 ± 0.01 and also a polar deficiency of ∼10%. Le Feuvre & Wieczorek (2008) analytically predicted the lunar pole/equator ratio to be 0.90. Le Feuvre & Wieczorek (2011) followed and semi-analytically estimated that involving both longitudinal and latitudinal asymmetries, the lunar crater density was minimized at (90 • E, ±65 • N) and maximized at the apex with deviations of about 25% with respect to the average for the current Earth-Moon distance, resulting in an apex/antapex ratio of 1.37 and a pole/equator ratio of 0.80. From observations, Morota & Furumoto (2003) reported an apex/anntapex ratio of ∼1.5 for 222 rayed craters with d c > 5 km on lunar highlands, formed by NEOs in the past ∼1.1 Gyr.
However, according to Strom et al. (2015), for craters with diameters larger than 10 km, those that the MBAs formed exceed the other population by more than an order of magnitude. This underlines the need for taking the craters and the cratering asymmetry that the MBAs contribute to into account, and thus requires the awareness of the lunar orbit during the LHB, namely, the dominance of MBAs. The Earth-Moon system evolved drastically in the early history. As a general trend, it is accepted that the Moon has been receding from Earth because of tidal dissipation, but a wide spectrum of opinions exists for the details of this picture. After the timescale problem was raised (Gerstenkorn 1955;MacDonald 1964;Goldreich 1966;Lambeck 1977;Touma & Wisdom 1994) and then solved through ocean models (Hansen 1982;Webb 1982;Ross & Schubert 1989;Kagan & Maslova 1994;Kagan 1997;Bills & Ray 1999), it is known today that the tidal dissipation factor of Earth must have been much larger in the distant past, or in other words, the tidal friction was much weaker than today (e.g., Bills & Ray 1999). Because tidal friction depends on not only the lunar orbit but also the terrestrial ocean shape, which is associated with continental drift, the exact history of Earth's tidal dissipation factor and hence the early evolution of the lunar orbit cannot be confirmed.
Still, we list some efforts of ocean modelers here. Hansen (1982) was the first to introduce Laplace's tidal equations in calculations of oceanic tidal torque. He modeled four cases, combinations of two idealized continentalities and two frictional resistance coefficients, and found the Earth-Moon distances at 4.5 Gya ranging from 38 to 53 R ⊕ , leading to nearly the same values at 3.9 Gya when the LHB occurred. Webb (1980Webb ( , 1982 developed an average ocean model and obtained the dates of the Gerstenkorn event as 3.9 and 5.3 Gya with and without solid Earth dissipation included, respectively. According to these two results, the Earth-Moon distance at 3.9 Gya should be no larger than 25 R ⊕ or about 42 R ⊕ . The author claimed, however, that the results should only be taken qualitatively. Ross & Schubert (1989) simulated a coupled thermal-dynamical evolution of the Earth-Moon system based on equilibrium ocean model, resulting in a M being 47 R ⊕ and e M being 0.04 at 3.9 Gya when the evolution timescale and parameters such as final a M and e M were required to be realistic. Kagan & Maslova (1994) described a stochastic model that considered fluctuating effects of the continental drift. Their two-mode resonance approximations gave rise to evolutions of tidal energy dissipation that were consistent with global paleotide models. By adopting a reproduced tidal evolution with a timescale as close as possible to the realistic one, the Earth-Moon distance is estimated to be about 47 R ⊕ at 3.9 Gya. This timescale can vary by billions of years, of course, depending on the values of the resonance lifetime.
The uncertain lunar orbital status during the LHB means that a reliable relationship between the lunar orbit and the cratering asymmetry is needed. If the influence of the former on the latter were significant and an incorrect condition of the Earth-Moon system were assumed, the estimated cratering asymmetry would be also incorrect. Conversely, this also implies that the early history of the Earth-Moon system could be inferred from the observed crater record if the portion of craters that formed during the LHB were selected. The influence of the lunar orbit and of the impactor population on the cratering asymmetry have not been sufficiently studied so far. Numerical simulations can only offer an empirical estimation based on a limited number of cases: Zahnle et al. (2001), who simulated impacts of ecliptic comets on giant planet satellites, combined the results of case α p = 2.0 and case α p = 2.5 with the satellite orbital speed v orb = 10.9 km s −1 and the impactor speed at infinity in the rest frame of the planet fixed at v ∞ ≈ 5 km s −1 , and thus derived a semi-empirical description of the crater density where β is the angular distance from the apex; using the debiased NEO model, Gallant et al. (2009) 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.

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.

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 the lunar 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 ∼17 R ⊕ , according to Le Feuvre & Wieczorek (2011). It also means the gravitational lensing by Earth is excluded, and moreover, the last assumption excludes the Earth's shielding for the Moon. This causes the symmetry between the near and far sides.
A precondition for the cratering asymmetry is that the Moon must be synchronously rotating, because otherwise any longitudinal variation would vanish as the lunar hemisphere facing the Earth changes constantly. According to Zhou & Lin (2008), the timescale for the Moon to reach the 1:1 spin-orbit resonance can be estimated by where τ tide is the tidal circularization timescale, m ⊕ is the Earth mass, n M , a M , R M , m M , and Q M 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.

Encounter geometry
We first consider the encounter condition before exhibiting the details of the impact in the following subsections. In the heliocentric frame, Earth orbits the Sun circularly at 1 AU, the Moon is ignored at this pre-encounter stage, and the impactors are all massless particles with semi-major axes a p > 1 AU and eccentricities 0 < e p < 1. If the orbit of an impactor intersects that of Earth on the ecliptic, meaning that its perihelion distance q p ≤ 1AU, an encounter is considered to occur at the mutual node. The relative velocity between the impactor and Earth when they encounter is the encounter velocity u enc = u p − u ⊕ , where u p and u ⊕ 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 semi-major axis is a ⊕ = 1 AU. The mutual node is where f p = f ⊕ = f enc , making the impactor's heliocentric distance equal to a ⊕ , that is, Thus, the encounter velocity depends on a p and e p (or q p ): u enc (a p , e p ) = u p (a p , e p , f enc (a p , e p )) − u ⊕ ( f enc (a p , e p )).
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 u p and u ⊕ 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 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. 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 u ⊕ (blue dashed arrow) and u p (black solid arrow) decreases as q p increases, and thus u enc (red solid arrow) shrinks.

Impact geometry
When the encounter between one impactor orbit and the Earth orbit occurs, as seen in the Earth-Moon system, the impactors in the common orbit are treated as particles approaching along the parallel straight lines in the direction of their common encounter velocity (Fig. 3a). They are assumed to be uniformly distributed in space, and to be numerous enough to cover the Moon's orbit, that is, the gravitational cross section of the Moon is always maximized (lunar diameter in the planar model).
In the geocentric frame with the x-axis pointing to the lunar perigee, the direction of u 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 u enc and the lunar orbital velocity u 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 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. 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 ⊕ . Now we establish a rest frame of the Moon that is to be used in the integration. As shown in Fig. 4, the positive y-axis points to the apex and the minus x-axis points to Earth. In this frame, 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 u and the local horizon; S is the cross section. A unit area (length in planar model) at longitude λ on the bombarded hemisphere is and its cross section is Denoting with e = (− cos λ, − sin λ) the normal vector, the incidence angle on the unit area can be derived with and the normal speed v ⊥ , the normal component of u is 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 A52, page 4 of 25 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 u, so that The bombarded hemisphere is then bounded by the longitude interval [λ l , λ u ], where the lower and upper limits λ l,u = λ ⊥ ∓90 • , and 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 It turns out

Exact formulations
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 Integration of dN over the interval [ f l , f u ] leads to A52, page 6 of 25 To simplify the expression, we define σ ∈ (0, π) by which leads to Therefore, 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, 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 γ p,c given by different studies vary, but their values are not important to the analytical derivation here. Still, they are needed for illustration and post-processing of numerical simulations, and we therefore adopt the forms given by Le Feuvre & Wieczorek (2011) for a non-porous regime (working for d c > 20 km) with the Moon as the target: where the non-porous 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 N p ∝ d −α p p , independent of its orbital distribution. Given a minimum of projectile size d min , a normalized size distribution is It can be interpreted as the probability of an arbitrary impactor to be larger than d p . The mean size of the projectiles is 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 di- ) 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 The crater density N c (>d c ) is the number of craters with diameters larger than d c per unit area. It is relevant to the age determination in the cratering chronology method. It depends on the projectile diameter d p , the projectile size-frequency distribu-tionN p , 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, In every unit time d f 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; functionN p 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 N p and dN determines the required part of dN. The substitute for the integral , 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 We note that N c (>d c , λ) describes not only the spatial distribution of the crater density, but also the craters' size distribution.
on any lunar longitude, and thus the slopes of the size distributions of impactors and craters are related by 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π): This form describes the impact density after a bombardment duration t, but its relative variation does not differ at all from the one-periodic form, which is one of the reasons we suggest using the asymmetry amplitude to measure the relative variation (Sect. 2.7). Hereafter, this new expression of N(λ) takes the place of the previous one, while N c (>d c , λ) (Eq. (56)) is still valid. All of the absolute spatial distributions are shown in Fig. 5 with solid curves. We caution that on each unit area, there must be impacts with all the (instantaneous) incidence angles θ ∈ [0 • , 90 • ] theoretically, and thus all the normal speeds v ⊥ ∈ [0, v] and crater diameters d c ∈ [0, d c (d p , v)], 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.

Near/far symmetry
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 Substituting λ with β, the integrated variables as functions of β are where A52, page 8 of 25 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).

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 where F is the impact flux (Eq. (20)), v x,y are elements of u (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 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π): Apparently, the global averages of V, Θ, and V ⊥ are the integrals vd 2 C, sin θd 2 C, and v ⊥ d 2 C over the above λ and f M intervals divided by C, respectively. The integration results in whereΘ is defined by sinΘ = sin θ, and it is foundV (45) and (44)) ofV ⊥ andN. 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 derivē 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Θ = 51.8 • , regardless of both the encounter speed v enc and the lunar orbital speed v M .

Approximate formulations
We have formulated the spatial distributions N(β), V(β), Θ(β), V ⊥ (β), D(β), and N c (> d c , β) (Eqs. (59)−(64)) and the global averagesN,V,V ⊥ ,Θ,D, andN c (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 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 (β) andN c , 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 (β) andN c 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 (β) andN c 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.
There are a few points about the approximate distributions to note. First, all the spatial distributions can be described by a simple function of β, 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 (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 → −d f 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 A D 1 is independent of this, and A c 1 only involves α p . That A c 1 has no dependence on d c means that the normalized crater density distribution is not a function of d c : N c (>d c , β)/N c (>d c ) = ∆N 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 A D,c 1 on A N,⊥ 1 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 A c 1 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 N c (>d c , β)/N c (>d c ) 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).

Reproducing bombardment conditions
To verify the validity of the approximate formulations (Eqs. (87)−(92)), we fit them using the least-squares method to the exact distributions (Eqs. (59)−(64)) that we considered as observed data. The fit cratering distributions are shown in Fig. 6 with dotted curves. Given α p , which only influences the relative variation of N c , it is seen that the smaller the η, the smaller the difference between the fit and exact distributions, because smaller η leads to a better approximation at first. The errors of A52, page 10 of 25 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 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). 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 left-hand sides are fits of A D 0,1 and A c 0,1 , 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, thatd p = d p (α p , d min ) (Eq. (53)), and that c c,p , γ 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 A D 0 = 33.0 km, A D 1 = 0.0238, A c 0 = 1.63 × 10 −7 km −1 , and A c 1 = 0.129 (with uncertainties 0.1%), among which A D 1 and A c 1 lead to η = (5.112 ± 0.003) × 10 −2 and α p = 1.738 ± 0.003. If v enc is known, then v M = 1.0224 ± 0.0006 km s −1 , d min = 0.495 ± 0.004 km, and Ft = (2.82 ± 0.04) × 10 −3 . Furthermore, if t is given, then F = (7.6 ± 0.1) × 10 −13 yr −1 . The reproduction errors, differences between the best-fit and the exact values, are 0.05%, 0.16%, 0.82%, and 1.44% for v M , α p , d min , and F, respectively. The goodness of reproduction sufficiently exhibits the viability of the reproducing method.

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

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 semi-major axis was a ⊕ = 1 AU, the eccentricity was e ⊕ = 0.0167, the inclination was i ⊕ = 0 • , the longitude of the ascending node was Ω ⊕ = 348.7 • , the argument of perihelion was ω ⊕ = 114.2 • , and the mean anomaly was M ⊕ = 0 • . Because the lunar early history involves a A52, page 12 of 25 great uncertainty as mentioned in Sect. 1, five cases were set with the lunar semi-major axis a M = 20, 30, 40, 50, and 60 R ⊕ in turn, probably covering most durations of the lunar history. This also helps to examine the influence of a M on the cratering. The lunar eccentricity e M was initially set to 0.04 (Ross & Schubert 1989), while the current value of the lunar inclination to the ecliptic i M = 5 • was adopted. The Moon's longitude of the ascending node Ω M , the argument of perihelion ω M , and the mean anomaly M M were all 0 • , because the lunar orbital period and the period of its precession along the ecliptic are extremely short relative to the timescale of the tidal evolution.
Primitive MBAs mean the asteroids that occupied the region where the contemporary main belt is during the LHB, that is, 2.0 AU ≤ a p ≤ 3.5 AU. Their orbits were disturbed by the migration of giants (Gomes et al. 2005), while their size distribution can be taken as invariant (Bottke et al. 2005). We only needed those MBAs that have the potential for encountering the Earth-Moon system, namely, those with perihelion distances q p no larger than 1 AU. To maximize the impact probability (Morbidelli & Gladman 1998) and avoid the interference of the effect of q p , the initial q p was constrained to be 1 AU, resulting in 0.50 ≤ e p ≤ 0.71. The inclination was set to 0 • ≤ i p ≤ 0.5 • to examine the mechanism of the latitudinal cratering asymmetry. The orbital elements Ω p , ω p , and M p all ranged from 0 to 360 • . Each particle's initial orbital elements were randomly generated in corresponding ranges following uniform distributions, except that its eccentricity was derived from e p = 1 − q p /a p . The mass of every particle was 4 × 10 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 N-body 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 post-processing when the crater diameter and crater density are calculated. Although the size distribution of current MBAs has been studied by several surveys (Sect. 1), it cannot be easily modeled because of its power-law breaks. In this work, we assumed the normalized size distribution (Eq. (52)) to be a single power-law characterized by α p = 3 and d min = 5 km referring to Ivezić et al. (2001), who determined α p = 3 for 5 km d p 40 km.
Three-dimensional N-body simulations were run to integrate the orbits of the Sun, Earth, the Moon, and 3 × 10 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 Runge-Kutta-Fehlberg method to integrate the recorded particles and the Earth-Moon system until all the particles had impacted the Moon. In this step, the perturbation of the Sun was ignored since the integration time it took is only ∼10 min, while the lunar orbital period is five days (a M = 20 R ⊕ ) at least.
Given that the Moon rotated synchronously during the LHB (Sect. 2.1), every impact in the simulations was precisely located on the lunar surface. Because the Moon has a longitudinal geometrical libration due to its eccentricity (the lunar spin axis was assumed to be perpendicular to the lunar orbital plane so that there is no libration in latitude), the near side is defined as the hemisphere facing Earth only when the Moon is at its perigee and the far side is its opposite. The longitude line that the center of the near side lies on is the prime meridian, with the east longitude values being positive. The apex point, the center of the leading side, is at (−90 • , 0 • ), while the antapex point, the center of the trailing side, is at (90 • , 0 • ). The low-and high-latitude regions are also defined here as the area with a latitude between ±45 • and the remaining area on the lunar surface, respectively.
Based on this numerical model, the bombardment conditions are v enc = 8.32 km s −1 (estimated with Eq. (7) given q p = 1 AU and a p = 2.75 AU, the mean of the initial a p ), α p = 3, d min = 5 km (d p = 7.5 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".

Preliminary statistics
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 ⊥ = v sin θ and the size of the formed crater centered on this impact location d c = d c (d p , v ⊥ ) 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, C j=1N p (> d p (d c , v ⊥, j )). This is how the size distribution of impactors is involved in post-processing. The given size d c was set to 100 km to ensure d p (d c , v ⊥ ) ≥ d min , otherwise the descriptions of N c distribution does not hold (Sect. 2.4). We note that the assumed values ofd p and d c do not influence the relative spatial variations of D and N c at all.
We summarize the regional Γ (denoting every/any one of V, Θ, V ⊥ , D, N, and N c following Sect. 2) for cases 1−5 in Fig. 9. The impact density of the leading side N ldg is much higher than that of the trailing side N trg for all the cases, with excesses of 25%−49% (equivalent to the excesses of C ldg to C trg ). Other Γ ldg also remain larger than Γ trg with notable excesses. These are clear signs of the leading/trailing asymmetry, indicating its existence in all the investigated aspects of cratering at any a M . It also shows unambiguously that except for V, the averages of the low-latitude region Γ low are always higher than those of 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 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. 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) andWieczorek (2011), who adopted the current impactors.
It is apparent thatV monotonically decreases from 9.15 to 8.80 km s −1 with a M increasing. The predictedV (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 v 2 enc + v 2 esc = 8.65 km s −1 , where the lunar escape speed is v esc = 2.38 km s −1 . Then the predictedV 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 ofV, the simulatedΘ are all obviously smaller than the analytical prediction of 51.8 • (Eq. (82)). The difference between sinΘ and π/4 is 13.7%−16.4% for the five cases. The reason probably is that our analytical model is planar, describing the variation of Θ on the equator, while in the three-dimensional simulations, the impacts elsewhere on the lunar surface are all statistically more oblique without isotropic impactors, and thus the simulatedΘ are diminished. Since theoreticallyΘ is constant regardless of how severe the leading/trailing asymmetry is (Sect. 2.6), and the pole/equtor asymmetry resulting from the concentration of low-inclination impactors should have no dependence on v M , the simulatedΘ were expected to be invariant with a M . In fact, we do findΘ of all the cases lying in a narrow range 41.1−42.6 • with a fluctuation of only 2%.V ⊥ andD show a mild dependence on a M , because their variations with a M are combinations of those ofV andΘ (V ⊥ = V sinΘ andD ∝V γ c ⊥ ). Owing to the fluctuation of theΘ variation, their negative correlations with a M are slightly obscured. The dependence onΘ also causes the simulatedV ⊥ , 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 simulatedD, 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 simulatedV ⊥ and the product of the sim-ulatedV and sinΘ is no more than 0.3% for all the cases, and the simulatedD is lower than c cd γ c pV γ c ⊥ 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, simulatedN shows a negative correlation with a M as well. (We note that the relative variation ofN also represents that of C gl .) The approximate expression ofN (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, whereN 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. Becausē N c ∝V γ p α p ⊥N , it is natural to see that the variation ofN c with a M combines those ofV ⊥ andN. Despite of the wavy fluctuation, N c generally decreases as a M increases. Additionally,N c is close to [d min /(c p d γ p c )] α pV γ p α p ⊥N 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 theV value of all cases and mildly overestimatesΘ and thusV ⊥ andD, while the correlations between global averages apply well.

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.
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 apex-antapex 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 apex-antapex 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 three-dimensional model the angular distance from the apex β is related to the longitude λ and latitude ϕ by cos β = − sin λ cos ϕ. (108) To describe the pure pole/equator asymmetry, we suggest an empirical function with a monotonic decrease from the equator to the poles. A spatial distribution of the coupled asymmetries is thus assumed to be It can be easily proved by integrating the right-hand side over the spherical surface that A 0 =Γ 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 sin r 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 least-squares method to the simulated spatial distributions shown in Fig. 10: N(λ, ϕ) = A N 0 (1 + A N 1 cos β)(1 + A N 2 cos 2ϕ), (115) N c (λ, ϕ) = A c 0 (1 + A c 1 cos β)(1 + A c 2 cos 2ϕ).
(116) 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).

Fit parameters and reproduced conditions
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, A 0 ≈Γ 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 A c 0 andN c , and the smallest is only 1% between A V 0 andV. 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 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, assumingd p = 7.5 km and d c = 100 km. the coupled-asymmetric distribution and the fit of A 0 is an excellent reproduction ofΓ. The reproduction goodness was expected to be proportional to √ a M (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, Our analytical model indicates that A 1 ∝ v M /v enc , that is, A 1 ∝ a −1/2 M with a given impactor population. As shown in Fig. 12b, the variation of fit A V 1 with a M reflects this relation very well: A V 1 (a M /R ⊕ ) 1/2 = 0.490 ± 0.010, 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 A Θ 1 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 A ⊥ 1 and A D 1 correlate mildly negatively with a M . The wavy decrease of fit A N 1 and the monotonic decrease of A c 1 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 γ c A ⊥ 1 and A D 1 and those between γ p α p A ⊥ 1 + A N 1 and A c 1 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 : 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 A V 0,1 . 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 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 A V 1 is obviously smaller than its prediction. Theoretically, η = A V 1 /(π/4), but the constant quotient of the fit A V 1 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 A52, page 18 of 25 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 semi-major axis. a) Every approximate global average A 0 is in terms of the relevant simulated global average.
squares connected by dashed lines in yellow, green, and blue). b) and c) A D 1,2 and A c 1,2 are compared to γ c A ⊥ 1,2 and γ p α p A ⊥ 1,2 + A N 1,2 (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 reproduced η is decreased. Although the reproduction from A V 1 is meant to be lower (Sect. 2.8), it is not enough to explain the difference of 32%. The reason why the fit A V 1 precisely follows A V 1 ∝ a −1/2 M 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 three-dimensional 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 simulatedV 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 A Θ 1 , 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 A D,c 0,1 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 S gl = 4πR 2 M , as the numerical simulations are three-dimensional. 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.
The asymmetry amplitudes A 1,2 must be in [0, 1] to ensure that ∆Γ is non-negative and its first and second factors are in negative correlations with β and |ϕ|. Clearly, the coupled distribution is symmetric about both the equator and the apex-antapex axis, with a maximum (1 + A 1 )(1 + A 2 ) always at the apex. If ϕ is constant, the function degenerates to ∆Γ = A 0 (1 − A 1 sin λ), where A 0 = 1 + A 2 cos 2ϕ and A 1 = A 1 cos ϕ. This means that on every latitude line, ∆Γ is maximized and minimized at longitude −90 • and +90 • with the average A 0 and the variation amplitude A 1 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. 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 latitude-parallel 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 , 1 7 and 3 5 , and two critical values of A 1 with A 2 given, 4A 2 5A 2 +1 and 2A 2 3(1−A 2 ) , which together define the configurations. We note that 4A 2 5A 2 +1 ≤ 2A 2 3(1−A 2 ) always holds and 4A 2 5A 2 +1 = 2A 2 3(1−A 2 ) = 1 3 only when A 2 = 1 7 .
1. When A 2 ∈ (0, 1 7 ] is fixed, two configurations appear in turn as A 1 increases. When A 1 ∈ (0, 4A 2 5A 2 +1 ], 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 The antapex acts as a saddle point. We call this distribution configuration 1. The moment A 1 increases to 4A 2 5A 2 +1 is when the minima join at the antapex and become one (Fig. 14d). After that when A 1 ∈ ( 4A 2 5A 2 +1 , 1], 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 A 2 ∈ ( 1 7 , 3 5 ] is fixed, another configuration presents itself as well as the above two. When A 1 ∈ (0, 4A 2 5A 2 +1 ], configuration 1 is also observed, even though at the moment A 1 = 4A 2 5A 2 +1 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 3(1−A 2 ) ] 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 The moment A 1 = 2A 2 3(1−A 2 ) is when each pair of approaching saddle and minimum combines and vanishes at the points (90 • , ± arccos( 1 3A 1 )) ( Fig. 14h). When A 1 ∈ ( 2A 2 3(1−A 2 ) , 1], 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 A 2 ∈ (0, 1 7 ]. 3. When A 2 ∈ ( 3 5 , 1] is fixed, the last stage of case A 2 ∈ ( 1 7 , 3 5 ] belonging to configuration 2 will not appear. The reason is that only if A 2 ≤ 3 5 can the critical value 2A 2 3(1−A 2 ) ≤ 1. When A 2 ∈ ( 3 5 , 1], A 1 can never reach 2A 2 3(1−A 2 ) 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 1 7 in any simulation case. In particular, the Θ and D distributions of all the cases belong to configuration 1. The N c distribution, with A c 2 always larger than 1 7 , 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 coupled-asymmetric distribution. Strom et al. (2005Strom et al. ( , 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. 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 A 2 ≤ 1 7 and A 1 ≤ 4A 2 5A 2 +1 ; d) transition between configurations 1 and 2 with A 2 ≤ 1 7 and A 1 = 4A 2 5A 2 +1 ; e) configuration 2 with A 2 ≤ 1 7 and A 1 > 4A 2 5A 2 +1 ; f) transition between configurations 1 and 3 with A 2 ∈ ( 1 7 , 3 5 ] and A 1 = 4A 2 5A 2 +1 ; g) configuration 3 with A 2 ∈ ( 1 7 , 3 5 ] and A 1 ∈ ( 4A 2 5A 2 +1 , 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 semi-major axis a p into account, that is, the distributions of perihelion distance q p are assumed to be the same. The semi-major 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 (v 2 imp = v 2 enc + v 2 esc ). 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.

Comparison between NEOs and MBAs
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 A D 1 because it is proportional to η (Eq. (106)): A D 1 is 0.0205−0.0229 for MBAs and 0.0215 for NEOs. Conversely, the fit A D 1 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 A D 1 of MBAs is precisely determined in theory, the fit A D 1 of a cratered area can lead to the realistic ratio of crater numbers of the two impactor populations.
The asymmetry amplitude A c 1 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. (2005Strom et al. ( , 2015, who found that the size distribution of craters on lunar young plains has a single slope of 2. Thus the A c 1 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 A c 1 of MBAs shown in Fig. 15 is defined by the top-right 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 A c 1 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 thusd p α p α p −1 d min , the definitiond p = d min +∞ d p dN p (Eq. (53)) still holds, so that the formulations of the D distribution andD also apply. The approximateD (Eq. (96)) only depends ond p 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 toD ∼ 33 km for MBAs and D ∼ 4 km for NEOs, assumingd p ∼ 1 km andd p ∼ 0.1 km, respectively.

Apex/antapex and pole/equator ratios
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 A52, page 22 of 25 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 maximum-minimum 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 A c 1 = 0.156 and A c 2 = 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 A c 1 , 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 power-law size distribution with this slope is found to be a convenient substitute for the 10th-order 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 A c 1 = 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 semi-analytical 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 A c 1 (Eq. (107)) as a function of a M , that is, 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.

Generalization of the analytical model
In addition to the Earth-Moon system, our analytical model (Sect. 2) can be applied to other planet-satellite 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 semi-major 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 tidal-locked 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 u enc is in the same direction as u 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 / v 2 orb − v 2 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 semi-empirical 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 size-distribution 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.

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 time-consuming 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 Γ(β) =Γ(1 + A 1 cos β), whereΓ is in positive correlation with v enc , except thatΘ is the constant 51.8 • , and A 1 ∝ v M /v enc , that is, A 1 ∝ a −1/2 M 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 cos 2ϕ), leading to the description of coupled asymmetries Γ(λ, ϕ) =Γ(1 + A 1 cos β)(1 + A 2 cos 2ϕ), which is well fit to the simulated Γ distributions with errors of Γ smaller than 6%. For each case, the analytical prediction ofV is only a few percent lower than the simulatedV, but the pre-dictedΘ 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 A Θ 1 , 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 A c 1 = 0.101−0.159, while that generated by the NEOs is A c 1 = 0.117, 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, A c 1 = 0.156, and A c 2 = 0.093. Our analytical model is applicable to other planet-satellite systems as long as v enc > v orb and will be generalized to the case v enc < v orb in future.