EDP Sciences
Open Access
Issue
A&A
Volume 594, October 2016
Article Number A52
Number of page(s) 25
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201628598
Published online 12 October 2016

© ESO, 2016

Licence Creative Commons
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

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

The cratering asymmetry depends on not only the lunar orbit, but also on the impactor population. After analyzing the size distributions of craters, Strom et al. (2005, 2015) suggested that there are two impactor populations in the inner solar system: the main-belt asteroids (MBAs), which dominated during the late heavy bombardment (LHB) ~ 3.9 Gya, and the near-Earth objects (NEOs), which have dominated since about 3.83.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 ap = 2.0−3.5 AU were the same as today, their eccentricities and inclinations were greatly excited by migrating giant planets (Gomes et al. 2005). Instead, their size distribution has changed little after the first ~100 Myr (Bottke et al. 2005). Several surveys of the current main belt have reported power-law breaks of its size distribution. The Sloan Digital Sky Survey (SDSS; Ivezić et al. 2001) distinguished two types of asteroids with different albedos (0.04 and 0.16 for blue and red asteroids) and found that for the whole sample throughout the whole belt, the cumulative size distribution has a slope αp = 1.3 for diameters dp = 0.4−5 km and 3 for 540 km. Parker et al. (2008) also confirmed this conclusion based on the SDSS moving object catalog 4, including ~88 000 objects. The first (Yoshida et al. 2003) and second (Yoshida & Nakamura 2007) runs of the Sub-km Main Belt Asteroid Survey (SMBAS), which each found 861 MBAs in R-band and 1001 MBAs in both B and R-bands, reported αp = 1.19 ± 0.02 for dp = 0.5−1 km and αp = 1.29 ± 0.02 for dp = 0.6−1 km, claimed to be consistent with SDSS for MBAs smaller than 5 km. The size distribution of craters formed by MBAs also has a complex shape, whose cumulative slope is αc = 1.2 for crater diameters dc ≲ 50 km, 2 for 50100 km, and 3 for 100300 km (Strom et al. 2015).

The NEO population is relatively well understood and suggested to have been in steady state for the past ~3 Gyr, constantly resupplied mainly from the main belt (Bottke et al. 2002). Bottke et al. (2002) derived the debiased orbital and size distributions of NEOs by fitting a model population to the known NEOs: the orbits are constrained in a range ap = 0.5−2.8 AU, ep< 0.8, ip< 35°, the perihelion distance qp< 1.3 AU, and the aphelion distance Qp> 0.983 AU; the size distribution is characterized by a single slope αp = 1.75 ± 0.1 for dp = 0.2−4 km. The corresponding slope for the crater size distribution indicated by Strom et al. (2015) is αc = 2 for dc = 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 dc> 5 km on lunar highlands, formed by NEOs in the past ~1.1 Gyr.

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

Still, we list some efforts of ocean modelers here. Hansen (1982) was the first to introduce Laplace’s tidal equations in calculations of oceanic tidal torque. He modeled four cases, combinations of two idealized continentalities and two frictional resistance coefficients, and found the Earth–Moon distances at 4.5 Gya ranging from 38 to 53 R, leading to nearly the same values at 3.9 Gya when the LHB occurred. Webb (1980, 1982) developed an average ocean model and obtained the dates of the Gerstenkorn event as 3.9 and 5.3 Gya with and without solid Earth dissipation included, respectively. According to these two results, the Earth–Moon distance at 3.9 Gya should be no larger than 25 R or about 42 R. The author claimed, however, that the results should only be taken qualitatively. Ross & Schubert (1989) simulated a coupled thermal-dynamical evolution of the Earth–Moon system based on equilibrium ocean model, resulting in aM being 47 R and eM being 0.04 at 3.9 Gya when the evolution timescale and parameters such as final aM and eM 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 vorb = 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 (1)where β is the angular distance from the apex; using the debiased NEO model, Gallant et al. (2009) ran simulations with Earth–Moon distances of aM = 50, 38, 30, 20, and 10 R, and confirmed the negative correlation between aM and leading/trailing asymmetry degree. A semi-analytical method capable of producing cratering of the Moon caused by current asteroids and comets in the inner solar system was proposed by Le Feuvre & Wieczorek (2011), who suggested a fit relation (2)which is valid between 20 and 60 R.

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

2. Analytical deduction

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

2.1. Assumptions and precondition

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

  • 1.

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

  • 2.

    The orbits of Earth and the Moon are circular.

  • 3.

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

  • 4.

    Earth is treated as a particle.

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

A precondition for the cratering asymmetry is that the Moon must be synchronously rotating, because otherwise any longitudinal variation would vanish as the lunar hemisphere facing the Earth changes constantly. According to Zhou & Lin (2008), the timescale for the Moon to reach the 1:1 spin-orbit resonance can be estimated by (3)where τtide is the tidal circularization timescale, m is the Earth mass, nM, aM, RM, mM, and are the mean motion, semi-major axis, radius, mass, and effective tidal dissipation factor of the Moon. A measure of this timescale is τsyn ≲ 102 yr for aM = 10 R, several orders shorter than the timescale of the tidal evolution. This precondition is therefore expected to hold for the main lifetime of the Moon, including during the LHB.

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

2.2. Encounter geometry

We first consider the encounter condition before exhibiting the details of the impact in the following subsections. In the heliocentric frame, Earth orbits the Sun circularly at 1 AU, the Moon is ignored at this pre-encounter stage, and the impactors are all massless particles with semi-major axes ap> 1 AU and eccentricities 0 <ep< 1. If the orbit of an impactor intersects that of Earth on the ecliptic, meaning that its perihelion distance qp ≤ 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 venc = vpv, where vp and v are orbital velocities of the impactor and Earth. Without loss of generality, the impactor’s and Earth’s arguments of perihelion are both assumed to be , so that their velocities in heliocentric ecliptic coordinates are where G is the gravitational constant, M is the solar mass, fp and f are true anomalies, and the terrestrial semi-major axis is a = 1 AU. The mutual node is where fp = f = fenc, making the impactor’s heliocentric distance equal to a, that is, (6)Thus, the encounter velocity depends on ap and ep (or qp): (7)For a fixed ap, as shown in Fig. 1, the encounter speed venc is minimized when qp = 1 AU. As qp decreases (ep increases), even though the magnitudes of the two orbital velocities vp and v in the moment of encounter are invariant, the angle between them expands, and thus the encounter speed venc increases. On the other hand, if qp is fixed, larger ap can also lead to higher venc, because the impactor orbital speed in the moment of encounter (8)is greater. As shown in Fig. 2, venc increases rightward (ap increases) and downward (qp decreases). Given a population of impactors, its typical venc 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.

thumbnail Fig. 1

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

Open with DEXTER

2.3. Impact geometry

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

In the geocentric frame with the x-axis pointing to the lunar perigee, the direction of venc can be random. Since we average over the lunar period, it can be assumed to be in the positive y-axis direction without loss of generality. Figure 3a shows that the angle between the encounter velocity venc and the lunar orbital velocity vM is the Moon’s true anomaly fM, which varies uniformly at the rate of the mean motion nM. The impact velocity is the relative velocity between the impactors and the Moon: Its magnitude, that is, the impact speed, is (12)At any instant, the impact can only occur on one hemisphere, and the impactors in this hemisphere have an equal v since they share the same venc. 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.

thumbnail Fig. 2

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

Open with DEXTER

Given the encounter speed venc and lunar orbital speed vM, as long as venc>vM, the impact speed v is always minimized when the Moon is at the perigee (fM = 0°) and maximized at the apogee (fM = 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 venc 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 venc>vM is taken for granted hereafter in this work. This holds for the MBAs, whose minimum venc is 7.9 km s-1 (ap = 2.5 AU and qp = 1 AU), which is nearly equal to the maximum vM when aM = 1 R. It is also valid for those NEOs with ap> 1.2 AU, whose minimum encounter speed venc = 2.4 km s-1 (ap = 1.2 AU and qp = 1 AU) approximates to vM for aM = 11R, while aM during the dominant epoch of the NEOs is at least about 40 R.

thumbnail Fig. 3

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

Open with DEXTER

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 to 90°, the angle between v and the local horizon; S is the cross section. A unit area (length in planar model) at longitude λ on the bombarded hemisphere is (16)and its cross section is (17)Denoting with e = (−cosλ,−sinλ) the normal vector, the incidence angle on the unit area can be derived with (18)and the normal speed v, the normal component of v is (19)Denoting with ρ the uniform spatial density of the impactors, the impact flux, the number of impactors received per unit time by the unit area is (20)At any instant when the lunar true anomaly is fM, the longitude of the normal impact point is defined as λ, where its normal vector e is parallel to v, so that (21)The bombarded hemisphere is then bounded by the longitude interval [ λl,λu ], where the lower and upper limits λl,u = λ ∓ 90°, and (22)As the bombarded hemisphere moves westward along the lunar equator, the moment a certain position λ enters this hemisphere and the moment it leaves are when λ becomes the lower and upper limits λl,u, respectively, and when sinθ = 0 there. Therefore, the time interval while a certain position λ is in the bombarded hemisphere, characterized by [ fl,fu ], can be derived with (23)It turns out

thumbnail Fig. 4

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

Open with DEXTER

2.4. Exact formulations

thumbnail Fig. 5

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

Open with DEXTER

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

The impact density N is the number of impacts on a certain region divided by its area. For the unit area on longitude λ illustrated in Fig. 4, its impact density during a unit time is (26)Integration of dN over the interval [ fl,fu ] leads to (27)To simplify the expression, we define σ ∈ (0) by (28)which leads to Therefore, (31)The impact speed V, incidence angle Θ, and normal speed V, as functions of longitude, are averages of v, θ, and v over one lunar period, that is, , , and . The integration leads to where Functions F(φ | k2) and F(k2) are the incomplete and complete elliptic integrals of the first kind; E(φ | k2) and E(k2) 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 dc and projectile diameter dp with The constants cp,c, γp,c, and given by different studies vary, but their values are not important to the analytical derivation here. Still, they are needed for illustration and post-processing of numerical simulations, and we therefore adopt the forms given by Le Feuvre & Wieczorek (2011) for a non-porous regime (working for dc> 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 , independent of its orbital distribution. Given a minimum of projectile size dmin, a normalized size distribution is (52)It can be interpreted as the probability of an arbitrary impactor to be larger than dp. The mean size of the projectiles is (53)where the slope αp> 1 (otherwise, the mean projectile size would be infinite). This equation should hold everywhere on the Moon since the size distribution is independent of orbital distribution, and gravitations on impactors are ignored. Therefore, the spatial distribution of the periodic averaged crater diameter is , which is substituted with to avoid difficult analytical integration. The deviation of the substitute is no more than 1.5% for the ratio vM/venc as high as 0.25 (the case aM = 10 R and venc = 10 km s-1) and even lower for a lower speed ratio. Thus we derive (54)The crater density Nc( >dc) is the number of craters with diameters larger than dc per unit area. It is relevant to the age determination in the cratering chronology method. It depends on the projectile diameter dp, the projectile size-frequency distribution , and the impact density N. On the assumption that every impact leaves one and only one crater, meaning that saturation, erosion, and secondary craters are all ignored, (55)In every unit time dfM ∈ [ fl,fu ], the factor dN gives the total density of craters on the longitude λ; function dp gives the projectile size required to form a crater as large as dc there; function gives the fraction of impactors larger than dp, which is equivalent to the fraction of craters larger than dc; and finally the product of and dN determines the required part of dN. The substitute for the integral is , whose deviation is no more than 0.07% for a ratio vM/venc as high as 0.25 and even smaller for a lower speed ratio. Assuming dc is large enough to ensure dp(dc,V(λ)) ≥ dmin holds everywhere on the lunar surface, then (56)We note that Nc( >dc) describes not only the spatial distribution of the crater density, but also the craters’ size distribution. Equation (56) proves on any lunar longitude, and thus the slopes of the size distributions of impactors and craters are related by (57)Although the integration is made over one lunar period, all of the distributions but N(λ) and Nc( >dc) have no dependence on time, so that they are applicable to any time interval t (multiple of lunar period, theoretically) as long as vM and venc are constant. We rewrite N(λ) after multiplying it by nMt/ (2π): (58)This form describes the impact density after a bombardment duration t, but its relative variation does not differ at all from the one-periodic form, which is one of the reasons we suggest using the asymmetry amplitude to measure the relative variation (Sect. 2.7). Hereafter, this new expression of N(λ) takes the place of the previous one, while Nc( >dc) (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 dc [0, ], nevertheless, their periodic averages are what the figure shows. The most significant common feature is that the curves all peak at the apex (λ = −90°) and are lowest at the antapex (λ = + 90°), regardless of venc. 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 venc leads to an upward shift of all the curves, expect for Θ(λ). An increase in venc diminishes and enlarges the absolute amplitudes of D and Nc, while it almost does not influence those of N, V, and V. With vM fixed, the absolute amplitude of Θ that only depends on the speed ratio vM/venc also seems to decrease as venc increases. We provide the explanation in Sect. 2.7.

2.5. Near/far symmetry

Because N(λ) = N( ± 180°−λ), N(λ) is symmetric about λ = ± 90°, meaning that the distribution of the impact density on the lunar surface is symmetric about the line connecting the apex and the antapex. The same is easily found in terms of V, Θ, V, D, and Nc. 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 aM< 25 R; Gallant et al. (2009) claimed that there was very little asymmetry for aM in the range 1050 R; Le Feuvre & Wieczorek (2011) found the asymmetry negligible with about 0.1% more craters of the near side than the far side, in contrast with Le Feuvre & Wieczorek (2005), who claimed a factor of four enhancement on the near side using impactors coplanar with Earth. Our deduction shows that without the Earth’s effect, the near/far asymmetry cannot exist, while the leading/trailing asymmetry is the inherent and natural consequence of cratering.

We therefore use β, the angular distance from the apex, to describe the leading/trailing asymmetry and not the longitude λ. It ranges between 0° and 180°, where the apex and antapex are, respectively. For the near side, where − 90° ≤ λ ≤ + 90°, Substituting λ with β, the integrated variables as functions of β are where We note that the above equations are applicable to both the near and the far side.

Figure 6 illustrates the relative spatial variations of N, V, Θ, V, D, and Nc 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 venc (the lower the speed ratio vM/venc with given vM), the smaller the relative variation amplitudes (to be explained in Sect. 2.7).

2.6. Global averages

Here we derive the global averages of the impact density N, impact speed V, incidence angle Θ, normal speed V, crater diameter D, and crater density Nc through integrating. The impact number for a unit area on the longitude λ during a unit time is (78)where F is the impact flux (Eq. (20)), vx,y are elements of v (Eqs. (14) and (15)), and ρ is the uniform spatial density of impactors. The total number of impacts on the global lunar surface after one lunar period is the integral over the λ interval [ λl,λu ] and fM interval [ 0°,360° ] in turn, where λl,u(fM) are boundaries of the instantaneous bombarded hemisphere (Eq. (22)). Thus, (79)The global average of N for the bombardment duration t (multiple of lunar period) is the total impact number for one period, C, divided by the global area lgl = 2πRM (perimeter for the planar model) and multiplied by nMt/ (2π): (80)

Apparently, the global averages of V, Θ, and V are the integrals !vd2C, !sinθd2C, and over the above λ and fM intervals divided by C, respectively. The integration results in

where is defined by , and it is found . The global averages of D and Nc are and , but to avoid analytical integration, we again substitute them by and , functions (Eqs. (45) and (44)) of and . The substitutes are quite acceptable, for their deviations are <3% and <1%, respectively, for a speed ratio vM/venc as high as 0.25 (the case aM = 10 R and venc = 10 km s-1). A lower speed ratio leads to even smaller deviations. Thus, we derive The global averages are shown in Fig. 5 with horizontal dashed lines. They all exhibit positive relations with venc, except for . Equation (82) explains the exception by determining the global average of sinθ to be π/ 4 always, equivalent to , regardless of both the encounter speed venc and the lunar orbital speed vM.

2.7. Approximate formulations

We have formulated the spatial distributions N(β), V(β), Θ(β), V(β), D(β), and Nc( >dc) (Eqs. (59)(64)) and the global averages , , , , , and (Eqs. (80)(85)). Hereafter, we use Γ to denote any one (or every one) of the variables N, V, Θ, V, D, and Nc, and use to denote its global average. The normalized distribution is defined as the absolute distribution divided by the global average, which describes the relative variation (86)The above variables represent every aspect of cratering and their spatial distributions are what we call cratering distribution. Each of the formulations Γ(β) and except for the constant can be taken as the product of two factors involving the encounter speed venc and the speed ratio η = vM/venc respectively. Recalling the assumption vM<venc, 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 Nc(β) and , the errors of the approximate forms are 1% for η = 0.1 and 5% even for η = 0.25. In particular, the errors of V(β) are merely 0.04% and 0.25%, respectively. The approximations of Nc(β) and also depend on αp. Given η = 0.1, their errors are 1%5% and 0.7%1.6% with αp increasing from 1 to 3, which is acceptable. However, for extreme conditions when η and αp are both large, approximate forms of Nc(β) and should be used with caution. For reference, the current Earth–Moon distance aM = 60 R (vM = 1.0 km s-1) and the current impactor population of NEOs (venc ≃ 20 km s-1 according to Gallant et al. 2009) result in an η of only about 0.05.

thumbnail Fig. 6

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

Open with DEXTER

There are a few points about the approximate distributions to note. First, all the spatial distributions can be described by a simple function of β, (98)where A0> 0 and 0 <A1< 1. This function is symmetric about the point (90°,A0) and has one maximum A0(1 + A1) at the apex (β = 0°) and one minimum A0(1−A1) 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 A0 and A1 entirely determine the variation of Γ. It is true that for every Γ, A0 is exactly the approximate and also the value of Γ on the prime meridian (β = 90°). Additionally, A0 has a positive relation with venc (except for Θ), explaining the dependence of on venc shown in Fig. 5. (The effect of η on the exact is negligible.) The factor (1 + A1cosβ) is then equal to the normalized distribution ΔΓ. The second parameter A1 just determines the amplitude of the relative variation shown in Fig. 6 and therefore is called “asymmetry amplitude”. For all Γ, it always holds that A1η, that is, the higher the venc or the lower the vM (the larger aM), 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 → −dfM/ dt = −nM, 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 A1.

Third, some simple relations between A0,1 of different variables are found. Parameters A0,1 of N, V, Θ, V, D, and Nc 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 A0, the approximate , are identical to the connections between V, Θ, V, D, and Nc themselves (Eqs. (54) and (56)). Furthermore, the last two relations between A1, the asymmetry amplitudes, are degenerate relations between A0: the product of (A)γpαp and AN becomes their sum, and then the exponents of A, both γc and γpαp, become its coefficients. The absolute distributions of D and Nc both depend on the impactors’ size distribution, which is characterized by αp and dmin, but is independent of this, and only involves αp. That has no dependence on dc means that the normalized crater density distribution is not a function of dc: .

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

2.8. Reproducing bombardment conditions

thumbnail Fig. 7

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

Open with DEXTER

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 Nc, it is seen that the smaller the η, the smaller the difference between the fit and exact distributions, because smaller η leads to a better approximation at first. The errors of fit distributions are <0.3% when η = 0.1 and αp = 1.75, almost half the errors of the approximate ones, except that for Θ they are nearly equal. The best fits are for the V distribution, with errors <0.02%. Furthermore, given the fits of A1, η can be easily reproduced based on the relations between A1 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 Nc (given αp = 1.75) give relatively large errors, while V gives the smallest; Θ and D tend to overestimate, while N, V, and Nc tend to underestimate. We caution that the observational error and limited data will prevent reaching such good fits and reproductions in reality, but our estimation shows what the series of approximate formulations is capable of.

The problem is that of N, V, Θ, V, D, and Nc, 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 Nc. By fitting formulations D(β) and Nc(β) to the observed spacial distributions of D and Nc and then solving the system of equations, where the left-hand sides are fits of and , we may reproduce the complete set of bombardment conditions: the lunar orbital speed vM and the impactors’ encounter speed venc, slope of the size distribution αp, minimum diameter dmin, impact flux F, and dominant duration t. We note that F = 2RMρvenc is used instead of ρ only for convenience in observation, that (Eq. (53)), and that cc,p, γc,p, , (Eqs. (46)–(51)) lgl, and dc are all constants. First, even without any knowledge of the above conditions, αp and η = vM/venc can be directly derived from Eqs. (106) and (107). Second, with αp and η obtained, knowing any of vM, venc, dmin, and Ft results in acquirement of all of them, because vM and venc are related by η, while venc, dmin, 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 Nc are taken again as the observed ones under current conditions where the lunar impactor population are the NEOs (Strom et al. 2005): vM = 1.022 km s-1, venc = 20 km s-1 (Gallant et al. 2009), αp = 1.75 (Bottke et al. 2002), dmin = 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 dmin and dc = 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 ) results in km, , , and (with uncertainties 0.1%), among which and lead to η = (5.112 ± 0.003) × 10-2 and αp = 1.738 ± 0.003. If venc is known, then vM = 1.0224 ± 0.0006 km s-1, dmin = 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 vM, αp, dmin, and F, respectively. The goodness of reproduction sufficiently exhibits the viability of the reproducing method.

3. Numerical simulations

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

3.1. Numerical model

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

The initial orbit of Earth was assumed to be the same as today. In the ecliptic coordinate system, its semi-major axis was a = 1 AU, the eccentricity was e = 0.0167, the inclination was i = 0°, the longitude of the ascending node was Ω = 348.7°, the argument of perihelion was ω = 114.2°, and the mean anomaly was M = 0°. Because the lunar early history involves a great uncertainty as mentioned in Sect. 1, five cases were set with the lunar semi-major axis aM = 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 aM on the cratering. The lunar eccentricity eM was initially set to 0.04 (Ross & Schubert 1989), while the current value of the lunar inclination to the ecliptic iM = 5° was adopted. The Moon’s longitude of the ascending node ΩM, the argument of perihelion ωM, and the mean anomaly MM 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 ap 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 qp no larger than 1 AU. To maximize the impact probability (Morbidelli & Gladman 1998) and avoid the interference of the effect of qp, the initial qp was constrained to be 1 AU, resulting in 0.50 ≤ ep ≤ 0.71. The inclination was set to 0° ≤ ip ≤ 0.5° to examine the mechanism of the latitudinal cratering asymmetry. The orbital elements Ωp, ωp, and Mp 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 ep = 1−qp/ap. The mass of every particle was 4 × 1012 kg, equivalent to a diameter dp ~ 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 × 1018 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 dmin = 5 km referring to Ivezić et al. (2001), who determined αp = 3 for 5 km dp 40 km.

Three-dimensional N-body simulations were run to integrate the orbits of the Sun, Earth, the Moon, and 3 × 105 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 (aM = 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 venc = 8.32 km s-1 (estimated with Eq. (7) given qp = 1 AU and ap = 2.75 AU, the mean of the initial ap), αp = 3, dmin = 5 km ( km assuming that the size distribution of the impactors is globally invariant), t = 107 yr, and vM = 1.78, 1.45, 1.26, 1.12, and 1.03 km s-1 (with eM ignored) for cases 15. Hereafter the above modeled bombardment conditions and the results we derived from them using our analytical formulations are called “predicted”, while those directly derived form simulations are called “simulated”.

3.2. Preliminary statistics

thumbnail Fig. 8

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

Open with DEXTER

thumbnail Fig. 9

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

Open with DEXTER

Cases 15 with aM = 20, 30, 40, 50, and 60 R in turn result in the global impact number Cgl = 1544, 1413, 1482, 1388, and 1391, respectively. Figure 8 shows the impact distribution on the lunar surface for case 1, and those of the other cases are similar. Every impact location was recorded along with the impact condition, namely, the impact speed v and the incidence angle θ. After this, the normal speed v = vsinθ and the size of the formed crater centered on this impact location were calculated. In a specific area, the averages of v, θ, v, and dc 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 Nc is the number of craters larger than a given size dc divided by the area. We point out that this crater number was obtained by totalling the probability of generating a crater larger than dc for every impact in this area, that is, . This is how the size distribution of impactors is involved in post-processing. The given size dc was set to 100 km to ensure dp(dc,v) ≥ dmin, otherwise the descriptions of Nc distribution does not hold (Sect. 2.4). We note that the assumed values of and dc do not influence the relative spatial variations of D and Nc at all.

We summarize the regional Γ (denoting every/any one of V, Θ, V, D, N, and Nc following Sect. 2) for cases 15 in Fig. 9. The impact density of the leading side Nldg is much higher than that of the trailing side Ntrg for all the cases, with excesses of 25%49% (equivalent to the excesses of Cldg to Ctrg). 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 aM. It also shows unambiguously that except for V, the averages of the low-latitude region Γlow are always higher than those of high-latitude region Γhigh, which confirms the pole/equator asymmetry regardless of aM. 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, Nnear is 8% smaller than Nfar for case 2, but 1%5% greater for other cases (the same for Cnear and Cfar). Therefore, the near/far asymmetry is clearly negligible for the whole range 2060 R even if it is present, in accordance with Gallant et al. (2009) and Le Feuvre & Wieczorek (2011), who adopted the current impactors.

It is apparent that monotonically decreases from 9.15 to 8.80 km s-1 with aM increasing. The predicted (Eq. (81)) is 8.608.41 km s-1 for aM = 2060 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 venc = 8.32 km s-1 with km s-1, where the lunar escape speed is vesc = 2.38 km s-1. Then the predicted is recalculated to be 8.928.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 qp = 1 AU. For the whole population of the MBAs without this constraint as well as the NEOs, vesc is even more negligible, so that the exclusion of the acceleration in the analytical model is quite acceptable.

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

and show a mild dependence on aM, because their variations with aM are combinations of those of and ( and ). Owing to the fluctuation of the variation, their negative correlations with aM are slightly obscured. The dependence on also causes the simulated , 5.816.03 km s-1, to be lower by 9.3%12.2% than the predicted values of 6.616.75 km s-1 (Eq. (83)); this is also true for the simulated , 96.898.6 km, which are lower by 6.7%8.3% than the predictions of 105.4106.5 km (Eq. (84)). However, the difference between the simulated and the product of the simulated and is no more than 0.3% for all the cases, and the simulated is lower than for only 2.6% at most. This indicates that although the pole/equator asymmetry in the Θ distribution leads to an overestimation of by the analytical model, the relationships between the global averages derived in Sect. 2.6 are still well applicable.

In general, simulated shows a negative correlation with aM as well. (We note that the relative variation of also represents that of Cgl.) The approximate expression of (Eq. (93)) indicates, however, that it is proportional to the constant venc, while the exact expression involving vM (Eq. (80)) is consistent with this negative correlation. Of the five sets of data points, case 2 (aM = 30 R) seems abnormal, where is oddly deficient and Nnear is clearly smaller than Nfar. 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 Nfar of this case is larger than case 3 (aM = 40 R) as expected, case 1 with its smaller aM should have suffered the same effect. However, in case 1, Nnear is even larger than Nfar. Therefore, we prefer to consider it as a possibility instead of confirming it as evidence of the near/far asymmetry. Because , it is natural to see that the variation of with aM combines those of and . Despite of the wavy fluctuation, generally decreases as aM increases. Additionally, is close to with excess of 8.6%9.7% for all the cases, in agreement with Eq. (85).

In short, this preliminary statistics confirms the presence of the leading/trailing asymmetry in spatial distributions of all Γ and the pole/equator asymmetry in those of all but V. All but show negative correlations with aM as indicated by the analytical model, which precisely predicts the value of all cases and mildly overestimates and thus and , while the correlations between global averages apply well.

3.3. Spatial distributions

The cratering distribution of case aM = 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 rs = 35° (along a great circle), equivalent to 1061 km on the lunar surface. Smaller rs can reveal more details, but because of the limited number of data points, we consider 35° as the proper choice for showing important features.

thumbnail 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 aM = 20 R. The value at each location is the average over a circular area centered on this location with a radius of 1061 km, assuming km and dc = 100 km.

Open with DEXTER

thumbnail Fig. 11

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

Open with DEXTER

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 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 Nc 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 090° (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 dc depend on both v and θ, the distributions of V and D naturally introduce both asymmetries in themselves. Similarly, because Nc depends on V and N, whose increases provide both a greater probability of generating large craters and more craters in total, the Nc 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 + A1cosβ), where in a three-dimensional model the angular distance from the apex β is related to the longitude λ and latitude ϕ by (108)To describe the pure pole/equator asymmetry, we suggest an empirical function (109)with a monotonic decrease from the equator to the poles. A spatial distribution of the coupled asymmetries is thus assumed to be (110)It can be easily proved by integrating the right-hand side over the spherical surface that still holds as Sect. 2.7 indicates. At the same time, A0 is the value of Γ at location (0°, ± 45°) or (180°, ± 45°). The parameters A1 and A2 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 rs, the asymmetry amplitudes are reduced. The decrease factor of A1 is estimated to be sinrs/rs, meaning that the fits of A1 here are the products of this factor and A1 given in Sect. 2.7. For rs = 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: Figure 11 shows the agreement between the simulated and fit variations of the above variables, with the fits of A0,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 A0,1,2, the above formulations can lead to fit spatial distributions whose configurations are determined by the relative magnitudes of A1 and A2 (Sect. 3.5).

3.4. Fit parameters and reproduced conditions

thumbnail Fig. 12

Fit A0a); A1b); and A2c) 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 A0 is in terms of the relevant simulated global average. , , and are compared to , , and (squares connected by dashed lines in yellow, green, and blue). b) and c) and are compared to and (squares connected by dashed lines in green and blue).

Open with DEXTER

thumbnail 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).

Open with DEXTER

The qualitative properties of the cratering distribution of case aM = 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 aM are shown and compared with the analytical predictions. Then the bombardment conditions are reproduced.

Fits of A0 of each Γ are plotted in terms of simulated in Fig. 12a. Theoretically, regardless of the pole/equator asymmetry (Sect. 3.3). In Fig. 12a, the fit A0 are all slightly smaller than the simulated , but only by a few percents, which is partly due to the approximation process itself. The largest difference is no more than 7% between and , and the smallest is only 1% between and . Therefore, even though we did not rigorously prove that the pole/quator asymmetry can be formulated as Eqs. (109), (110) is clearly a good description of the coupled-asymmetric distribution and the fit of A0 is an excellent reproduction of . The reproduction goodness was expected to be proportional to (inversely proportional to vM) with a given impactor population, but this trend is not observed because it is obscured by the statistical fluctuation. The variation of with aM has been discussed in Sect. 3.2. Additionally, the interrelations between A0 (Eqs. (99)(101)) are verified regardless of the involved pole/equator asymmetry: for any simulation case, approximates to with a difference of 0.5%; is almost equivalent to with a difference of <3%; and is close to with a difference of 11%.

Our analytical model indicates that A1vM/venc, that is, with a given impactor population. As shown in Fig. 12b, the variation of fit with aM reflects this relation very well: , 0.469 ± 0.008, 0.455 ± 0.013, 0.491 ± 0.004, and 0.493 ± 0.006 for case 15 in turn, which is almost a constant. Because the leading/trailing asymmetry in the Θ distribution is the faintest, the fit involves relatively great fluctuation, and its relation with aM is not clearly seen. The distributions of V and D depend on those of both V and Θ, so that the fit and correlate mildly negatively with aM. The wavy decrease of fit and the monotonic decrease of are also qualitatively consistent with the above theoretic relation. The interrelations between A1 (Eqs. (102) and (103)) also apply quite well. For all cases, the differences between and and those between and are 3%.

We have assumed that A2 is invariant with the lunar orbit (Sect. 3.3). As a result, no common dependence of fit A2 on aM 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 A0,1, we find that similar rules are present for A2: 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 A1, the speed ratio η and thus the lunar orbital speed vM were easily reproduced according to Eqs. (87)(92). Figure 13 shows the reproduced η and vM in comparison with their predictions. The decrease factor was considered in the reproduction. Generally speaking, the reproductions are reasonably consistent with each other. Still, η reproduced from is obviously smaller than its prediction. Theoretically, , but the constant quotient of the fit and predicted η is about 0.53 (with the decrease factor considered), which is smaller than the coefficient π/ 4 derived in Sect. 2.7, so that the reproduced η is decreased. Although the reproduction from is meant to be lower (Sect. 2.8), it is not enough to explain the difference of 32%. The reason why the fit precisely follows but underestimates the predicted η and vM 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 venc is greater than the predicted venc 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 venc with the simulated only reduces the reproduction error to about 27%, other places where the numerical simulations deviate from the analytical model need to be checked. Although statistical fluctuation results in a great error of η reproduced from , fortunately the reproductions from asymmetry amplitudes of V, D, N, and Nc are much better, especially for N, which leads to reproduction errors of 1%15% for all cases.

Again we followed the method of reproducing bombardment conditions described in Sect. 2.8 and present a simple estimation based on the fit of case 1, for only D and Nc are observables in surveys of crater records. Assuming venc = 8.32 km s-1 and t = 107 yr are known, the reproduced conditions are calculated to be vM = 1.45 ± 0.02 km s-1, αp = 3.23 ± 0.06, dmin = 4.71 ± 0.03 km, and F = (1.821 ± 0.164) × 10-4 yr-1. Comparisons of the former three to the predictions, vM = 1.78 km s-1, αp = 3, and dmin = 5 km, lead to differences of 19%, 8%, and 6%, respectively. We note that the impactor flux F is derived with Eq. (105) where lgl = 2πRM is substituted by , as the numerical simulations are three-dimensional. It is equivalent to a total number of impacts of 1821 ± 164 during an integration time of 107 yr, which is very close to the actual result C = 1544. Taking the limited number of impacts in simulations into account, it is rather inspiring to find that this method of reproducing the bombardment conditions from the cratering distribution works well. The importance of this method is fully shown when we recall that very many lunar craters were generated by MBAs during the LHB, meaning that our current knowledge of the LHB and MBAs may provide clues to the early lunar dynamical history.

3.5. Coupling of leading/trailing and pole/equator asymmetries

Figure 10 illustrates the spatial distributions of Γ. All of the distributions involve the leading/trailing asymmetry, and all but the V distribution also show the pole/equator asymmetry. The coupling of the two asymmetries can be described using Eq. (110), whose normalized form is (119)The asymmetry amplitudes A1,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 + A1)(1 + A2) always at the apex. If ϕ is constant, the function degenerates to , where and . This means that on every latitude line, ΔΓ is maximized and minimized at longitude − 90° and + 90° with the average and the variation amplitude dependent on | ϕ |. If λ is constant, the function is a cubic polynomial of cosϕ, whose first factor monotonically decreases or increases as | ϕ | increases on the leading or trailing side, respectively, and the second factor always decreases with increasing | ϕ |. Therefore, with a given λ ∈ (−180°,0°), ΔΓ is always maximized on the equator and minimized on the poles, while with a given λ ∈ (0°,180°), the conditions vary with varying A1,2. The distributions are thus classified into three configurations according to the numbers and locations of the extrema.

thumbnail Fig. 14

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

Open with DEXTER

We summarize the theory distribution of the coupled asymmetries as follows. The case A1 ≠ 0 and A2 = 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 A1 = 0 and A2 ≠ 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 A2, and , and two critical values of A1 with A2 given, and , which together define the configurations. We note that always holds and only when .

  • 1.

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

  • 2.

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

  • 3.

    When is fixed, the last stage of case belonging to configuration 2 will not appear. The reason is that only if can the critical value . When , A1 can never reach before increasing to 1.

Simply speaking, configurations 1, 2, and 3 are defined to have two, one, and three minima each and one maximum always. Our fit distributions of Θ, V, D, and N all belong to configuration 1 and 2 since none of their A2 is larger than in any simulation case. In particular, the Θ and D distributions of all the cases belong to configuration 1. The Nc distribution, with always larger than , belongs to configuration 2 for cases 13, configuration 1 for case 4 and configuration 3 for case 5. We note that we classify those distributions according to their fit A1 and A2, 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 Nc 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.

4. Discussion

4.1. Comparison between NEOs and MBAs

thumbnail Fig. 15

Contour map of the leading/trailing asymmetry amplitude of the crater density as a function of speed ratio η and size-distribution slope αp. The estimates of for NEOs and MBAs are indicated by dotted and dashed squares, respectively.

Open with DEXTER

Strom et al. (2005, 2015) suggested that the MBAs and NEOs that dominated during the LHB and since about 3.83.7 Gya in turn had formed craters on the Moon and even all the terrestrial planets. We review the current knowledge of two impactor populations in terms of orbital distribution, size distribution, etc. in Sect. 1. Now that we have formulated the cratering asymmetry based on our analytical model in Sect. 2, which has been verified in various aspects by numerical simulations in Sect. 3, we analytically estimate and compare the leading/trailing asymmetry amplitudes generated by the two impactor populations in this subsection.

Given a population of impactors, the typical encounter speed venc is determined by its orbital distribution, the minimum size dmin 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 Nc distribution. When the normalized distribution depending on just A1 is considered, then only venc 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 ap into account, that is, the distributions of perihelion distance qp are assumed to be the same. The semi-major axes of MBAs lie in the range of 2.03.5 AU, while those of NEOs lie in the range of 0.52.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 ap, 2.75 AU for MBAs and 1.65 AU for NEOs, and a median of qp, 0.5 AU, the typical encounter speeds of MBAs and NEOs are calculated to be 25.4 km s-1 and 22.1 km s-1 (Eq. (7)), and thus their typical impact speeds are 25.5 km s-1 and 22.3 km s-1 (). This estimate of vimp 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 vM during the dominant epochs of the two impactor populations are different as well. For MBAs, vM during the LHB may be from 1.12 to 1.26 km s-1, relevant to the uncertain aM from 50 to 40 R (Sect. 1). For NEOs, the current value vM = 1.03 km s-1 is used since aM is believed to have been beyond 40 R for most of the Moon lifetime. Here eM is ignored in the calculation.

Therefore, the speed ratio η is 0.0440.049 for MBAs and 0.046 for NEOs. A clear excess or deficiency of η of MBAs compared to NEOs is not seen. This is also true for the asymmetry amplitude because it is proportional to η (Eq. (106)): is 0.02050.0229 for MBAs and 0.0215 for NEOs. Conversely, the fit of an observed area on the lunar surface that is purely cratered by MBAs can lead to a speed ratio and thus an Earth–Moon distance during the LHB, helping to clarify the early lunar history. On the other hand, when of MBAs is precisely determined in theory, the fit of a cratered area can lead to the realistic ratio of crater numbers of the two impactor populations.

The asymmetry amplitude depends on not only η, but also on αp (Eq. (107)). Assuming the size distribution slope of the NEOs is αp = 1.75 (Bottke et al. 2002), Eq. (57) leads to the slope of craters generated by them of αc = 2.07, which is well consistent with Strom et al. (2005, 2015), who found that the size distribution of craters on lunar young plains has a single slope of 2. Thus the of NEOs is calculated to be 0.117. The size distribution of the MBAs is not a single power law as mentioned in Sect. 1, and the variation of its slope is not slight, as Ivezić et al. (2001) suggested for dp = 0.4−5 km and 540 km, αp = 1.3 and 3, respectively. Still, all of our analysis is valid for the separated part of impactors in a given dp interval with invariant αp and the craters they generate. Therefore, the possible range of of MBAs shown in Fig. 15 is defined by the top-right maximum 0.159, where η = 0.049 (aM = 40 R) and αp = 3, and the bottom left minimum 0.101, where η = 0.044 (aM = 50 R) and αp = 1.3. For the considered part of MBAs with dp = 0.4−40 km, Fig. 15 shows that their is probably greater than that of the NEOs. The Nc distribution generated by impactors with a size distribution of multiple power laws will be formulated in a future work.

We caution that even the size distribution is not a single power law, and thus , the definition (Eq. (53)) still holds, so that the formulations of the D distribution and also apply. The approximate (Eq. (96)) only depends on and venc. Because venc of the MBAs is greater than that of NEOs and MBAs are generally larger than NEOs, this results in the separation into two crater populations in terms of size. A rough estimation leads to 33 km for MBAs and 4 km for NEOs, assuming 1 km and 0.1 km, respectively.

4.2. Apex/antapex and pole/equator ratios

thumbnail Fig. 16

Reproduction of spatial distribution of the crater density in Le Feuvre & Wieczorek (2011). It is generated using Eq. (119) with asymmetry amplitudes A1 = 0.156 and A2 = 0.093.

Open with DEXTER

thumbnail Fig. 17

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

Open with DEXTER

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 Nc) at the apex to the antapex, which increases with increasing leading/trailing asymmetry, while the latter is the ratio of the value at the poles to the equator, which increases with decreasing pole/equator asymmetry. These ratios can be obtained easily and are especially convenient for simulations and observations, but their physical meanings were not understood before. Neither was it possible to completely define a spatial distribution through them. According to the formulation of coupled asymmetries that we established based on our analytical model and simulation results (Eq. (110)), the ratios are each related to the relevant asymmetry amplitude by (When A1 ≠ 0, Γ(equator) means the average over the equator.) Because a whole distribution can be directly reproduced given A1 and A2 using Eq. (119) and because the meaning of A1 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 Nc 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 Nc 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 = 0.156 and = 0.093, which are calculated with Eqs. (122) and (120) to ensure Nc(apex) /Nc(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 Nc distribution is well reproduced, but also the estimate αp = 2.6 is derived from Eq. (107) with the reproduced , current vM, and the mean impact speed on the Moon venc = 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 Nc distribution generated by NEOs to be = 0.117 (Sect. 4.1), which is equivalent to Nc(apex) /Nc(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 Nc(apex) /Nc(antapex) on aM and obtained Eq. (2) that is valid for aM = 20−60 R. Treating the expression of (Eq. (107)) as a function of aM, that is, where , we can also derive a relation (124)The two relations are compared in Fig. 17 with solid and dashed lines. We note that in calculating our relation, venc = 22.1 km s-1 and αp = 1.75 are the same as used in Sect. 4.1, but aM 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 venc and adopt αp = 2.6 as reproduced from their Nc distribution, our function is identical to their valid part with a difference of only <2%.

We briefly point out the difference between our work and that of Le Feuvre & Wieczorek (2011) in cratering modeling. Their model is able to produce the pole/equator asymmetry since the anisotropic impactor population was adopted, while our current analytical model alone can only formulate the leading/trailing asymmetry, but with our empirical description of the pole/equator asymmetry based on numerical simulations, we can also describe, reproduce, and measure the coupled asymmetries. More importantly, they integrated the cratering distribution numerically, while we give explicit analytical functions after rigorous deduction, which clearly shows the influence of every factor and makes the formulations very easy to apply.

4.3. Generalization of the analytical model

In addition to the Earth–Moon system, our analytical model (Sect. 2) can be applied to other planet-satellite systems. In the inner solar system, the NEOs and MBAs are the common impactor populations dominating different epochs. They encounter Mars with venc 19.3 and 22.9 km s-1, respectively, which is calculated with a replaced by the semi-major axis of Mars aMars = 1.5 AU in Eq. (7) and the typical ap of NEOs being the median of the interval 0.762.8 AU instead of 0.52.8 AU because projectiles with ap<aMars/ 2 on elliptic orbits cannot encounter Mars. The Martian moons, Phobos and Deimos, which lie on very circular orbits, have current orbital speeds of vorb = 2.14 and 1.35 km s-1 and negligible escape speeds vesc. 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 venc>vorb, 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 venc is in the same direction as vorb) and apocenter, while the impact speed v is minimized and maximized (Sect. 2.3). However, if venc<vorb, 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(venc/vorb) ∈ (0°,90°), turns around and returns to the apex when f = 180°; as f continues increasing, it continues moving westward until f = 360°−arccos(venc/vorb) ∈ (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(venc/ 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 vorb/venc, the smaller the distance βmax, and thus the greater the asymmetry. As a result, the leading/trailing asymmetry is much enhanced when venc<vorb 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 Nc distribution on satellites of giant planets shown as Eq. (1). We find qualitative consistence of our approximate formulation Nc(β) (Eq. (92)) with theirs since the asymmetry is amplified by speed ratio vorb/venc 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 venc<vorb and the influence of αp was overestimated by formulation of Zahnle et al. (2001), which will be examined by future model.

5. Conclusions

We analytically formulated the lunar cratering distribution and confirmed the derivation with numerical simulations. The formulations are quite easy to use. They are able to give results nearly identical to related works, avoiding the 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 Nc. Series expansion to the first order of vM/venc resulted in the approximate Γ distributions in the common form , where is in positive correlation with venc, except that is the constant 51.8°, and A1vM/venc, that is, with a given venc. 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 vM/venc, the better the approximation to the exact Γ distributions and the better the reproduction of the bombardment conditions. In particular, the approximation of Nc 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 aM = 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 + A2cos2ϕ), leading to the description of coupled asymmetries , which is well fit to the simulated Γ distributions with errors of smaller than 6%. For each case, the analytical prediction of is only a few percent lower than the simulated , but the predicted is obviously larger than the simulated because the pole/equator asymmetry is involved. As predicted by the analytical model, except for all show negative correlations with aM in general, and so do the fit A1 except for , whose variation with aM involves greater statistical fluctuations, while the fit A2 seem independent of aM. Additionally, the method of reproducing performs well, with reproduction errors being ~10% for bombardment conditions in case 1. This means that it is possible to speculate about the lunar orbital status during the LHB when the cratering distribution generated by MBAs is determined by observations.

Based on our analytical model, the leading/trailing asymmetry amplitude of the crater density that is generated by the MBAs is estimated to be , while that generated by the NEOs is , which is equivalent to Nc(apex) /Nc(antapex) = 1.27. We easily but precisely reproduced the Nc distribution and the variation of Nc(apex) /Nc(antapex) as a function of aM derived by Le Feuvre & Wieczorek (2011), with αp = 2.6, , and = 0.093. Our analytical model is applicable to other planet-satellite systems as long as venc>vorb and will be generalized to the case venc<vorb in future.

Acknowledgments

We thank the referee Bill Bottke for constructive and inspiring comments. This research has been supported by the Key Development Program of Basic Research of China (Nos. 2013CB834900), the National Natural Science Foundations of China (Nos. 11003010 and 11333002), the Strategic Priority Research Program “The Emergence of Cosmological Structures” of the Chinese Academy of Sciences (Grant No. XDB09000000), the Natural Science Foundation for the Youth of Jiangsu Province (No. BK20130547) and the 985 project of Nanjing University and Superiority Discipline Construction Project of Jiangsu Province.

References

All Figures

thumbnail Fig. 1

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

Open with DEXTER
In the text
thumbnail Fig. 2

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

Open with DEXTER
In the text
thumbnail Fig. 3

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

Open with DEXTER
In the text
thumbnail Fig. 4

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

Open with DEXTER
In the text
thumbnail Fig. 5

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

Open with DEXTER
In the text
thumbnail Fig. 6

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

Open with DEXTER
In the text
thumbnail Fig. 7

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

Open with DEXTER
In the text
thumbnail Fig. 8

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

Open with DEXTER
In the text
thumbnail Fig. 9

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

Open with DEXTER
In the text
thumbnail 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 aM = 20 R. The value at each location is the average over a circular area centered on this location with a radius of 1061 km, assuming km and dc = 100 km.

Open with DEXTER
In the text
thumbnail Fig. 11

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

Open with DEXTER
In the text
thumbnail Fig. 12

Fit A0a); A1b); and A2c) 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 A0 is in terms of the relevant simulated global average. , , and are compared to , , and (squares connected by dashed lines in yellow, green, and blue). b) and c) and are compared to and (squares connected by dashed lines in green and blue).

Open with DEXTER
In the text
thumbnail 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).

Open with DEXTER
In the text
thumbnail Fig. 14

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

Open with DEXTER
In the text
thumbnail Fig. 15

Contour map of the leading/trailing asymmetry amplitude of the crater density as a function of speed ratio η and size-distribution slope αp. The estimates of for NEOs and MBAs are indicated by dotted and dashed squares, respectively.

Open with DEXTER
In the text
thumbnail Fig. 16

Reproduction of spatial distribution of the crater density in Le Feuvre & Wieczorek (2011). It is generated using Eq. (119) with asymmetry amplitudes A1 = 0.156 and A2 = 0.093.

Open with DEXTER
In the text
thumbnail Fig. 17

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

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.