Mathieu Le Feuvre
email: mathieu.lefeuvre@univ-nantes.fr

Mark A Wieczorek

Nonuniform cratering of the Moon and a revised crater chronology of the inner solar system

Keywords: cratering, Moon, terrestrial planets, crater chronology, impact processes

inaccurate for periods prior to 3.5 Ga because of a different impactor population, or for craters smaller than a few kilometers on Mars and Mercury, due to the presence of subsurface ice and to the abundance of large secondaries, respectively. Standard parameter values allow for the first time to naturally reproduce both the size distribution and absolute number of lunar craters up to 3.5 Ga ago, and give self-consistent estimates of the planetary cratering rates relative to the Moon.

Email: mathieu.lefeuvre@univ-nantes.fr

We model the cratering of the Moon and terrestrial planets from the present knowledge of the orbital and size distribution of asteroids and comets in the inner solar system, in order to refine the crater chronology method. Impact occurrences, locations, velocities and incidence angles are calculated semianalytically, and scaling laws are used to convert impactor sizes into crater sizes. Our approach is generalizable to other moons or planets. The lunar cratering rate varies with both latitude and longitude: with respect to the global average, it is about 25 % lower at (±65 • N, 90 • E) and larger by the same amount at the apex of motion (0 • N, 90 • W) for the present Earth-Moon separation. The measured size-frequency distributions of lunar craters are reconciled with the observed population of near-Earth objects under the assumption that craters smaller than a few kilometers in diameter form in a porous megaregolith. Varying depths of this megaregolith between the mare and highlands is a plausible partial explanation for differences in previously reported measured size-frequency distributions. We give a revised analytical relationship between the number of craters and the age of a lunar surface. For the inner planets, expected size-frequency crater distributions are calculated that account for differences in impact conditions, and the age of a few key geologic units is given. We estimate the Orientale and Caloris basins to be 3.73 Ga old, and the surface of Venus to be 240 Ma old. The terrestrial cratering record is consistent with the revised chronology and a constant impact rate over the last 400 Ma. Better knowledge of the orbital dynamics, crater scaling laws and megaregoltih properties are needed to confidently assess the net uncertainty of the model ages that result from the combination of numerous steps, from the observation of asteroids to the formation of craters. Our model may be 1 Introduction

The counting of impact craters offers a simple method to estimate the ages of geologic units on planetary surfaces when in situ rock samples are lacking. The crater chronology method is based on the simple idea that old surfaces have accumulated more impact craters than more recent ones [START_REF] Baldwin | The face of the Moon[END_REF]. The relationship between geologic age and number of lunar craters, based on the radiometric dating of existing lunar rock samples, is found to be approximately linear from the present to about 3 Ga ago, and approximately exponential beyond that time (Neukum et al., 2001a;Stöffler and Ryder, 2001, and references therein). For statistical robustness, craters are generally counted within a number of consecutive diameter ranges, allowing one to recognize certain biases, such as erosion, resurfacing or crater saturation processes. Measurements over various geologic units has led to the postulate that the relative shape of the size-frequency distributions (SFD) of lunar impact craters was similar for all surfaces, but unfortunately, the exact shape of this production function in the 2-20 km kilometer range is still debated after decades of study.

The total predicted size-frequency distribution for any given time is obtained by multiplying the production function, assumed independent of age, by a time-variable constant. The age of a geologic unit is estimated by finding the best fit between the standardized and measured distributions.

The approach taken in this work models directly crater distributions from the current knowledge of the impactor population. This allows us to infer properties of the impact history of a planetary body that the sole observation of craters could not reveal. In particular, whereas the crater chronology method assumes that craters accumulate uniformly on the surface of the planetary body, this method allows to quantify possible spatial variations in the impact rate. Moreover, in the absence of dated samples with known geological context from Mercury, Venus and Mars, the generalization of the crater chronology to these planets requires an estimate of the relative cratering rates with respect to the Moon: in the same period of time, the number of craters of a given size that form on two planets differ according to both the different impact probabilities of the planet crossing objects and the impact conditions (e.g., impact velocity, surface gravity).

Using this bottom-up approach, Neukum et al. (2001b) proposed age estimates of mercurian geologic features based on a Mercury/Moon cratering rate ratio estimated from telescopic observations. Similarly, [START_REF] Hartmann | Cratering chronology and the evolution of Mars[END_REF] adapted the lunar cratering chronology to Mars, using the Mars/Moon cratering rate ratio of [START_REF] Ivanov | Mars/Moon cratering rate ratio estimates[END_REF]. The venusian surface, which contains a small number of craters that appear to be randomly distributed, has been the subject of several attempts of dating. The most recent estimates can be found in [START_REF] Korycansky | Modeling crater populations on Venus and Titan[END_REF] and [START_REF] Mckinnon | Cratering on Venus: models and observations[END_REF], and are both based on the population of Venus crossers as estimated by [START_REF] Shoemaker | Asteroid and comet flux in the neighborhood of Earth[END_REF]. In this study, we use improved estimates of the orbital characteristics and size-frequency distribution of the impactor population.

The early lunar cratering record indicates that impactors were hundred of times more numerous than today, possibly due to the massive injection in the inner system of main-belt and/or Kuiper-belt objects about 700 Ma after the Moon formed, an event known as the Late Heavy Bombardment (LHB) (see [START_REF] Gomes | Origin of the cataclysmic late heavy bombardment period of the terrestrial planets[END_REF][START_REF] Tera | Isotopic evidence for a terminal lunar cataclysm[END_REF][START_REF] Hartmann | The time-dependent intense bombardment of the primordial Earth/Moon system[END_REF]. More than 3 Ga ago, the impactor population appears to have reached a relative state of equilibrium, being replenished both in size and orbit by, respectively, collisions inside the main asteroid belt and the ejection by resonances with the giant planets [START_REF] Bottke | Debiased orbital and absolute magnitude distribution of the Near-Earth Objects[END_REF]. Under the assumption of a steady state distribution of impactors, the distribution of craters on ∼ 3 Ga old surfaces should be consistent with the present astronomically inferred cratering rates.

The lunar cratering record, in the form of the standardized production functions, agrees reasonably well with telescopic observations of planet-crossing objects down to a few kilometers in diameter [START_REF] Stuart | Observational constraints on the number, albedos, sizes, and impact hazards of the near-Earth asteroids[END_REF][START_REF] Werner | NOTE: The Near-Earth Asteroid Size-Frequency Distribution: A Snapshot of the Lunar Impactor Size-Frequency Distribution[END_REF].

Smaller craters appear to be not numerous enough to have been formed by an impacting population similar to the present one. But, in [START_REF] Ivanov | Exogenic dynamics, cratering and surface ages[END_REF], it was suggested that the presence of a porous lunar megaregolith may reduce the predicted size of small craters, accounting for this observation. In [START_REF] Strom | The origin of planetary impactors in the inner solar system[END_REF], the distinction was made between pre-and post-LHB lunar crater distributions. Older distributions, depleted in small craters, would reflect the SFD of Main Belt asteroids -massively provided by the LHBrather than the SFD of Near-Earth asteroids. [START_REF] Marchi | A new chronology for the moon and mercury[END_REF] revised the crater chronology method by deriving a new production function from impact modeling that differs significantly from those based on measurements. In this study, we attempt to fully reconcile the measured lunar production functions with the telescopic observations of planet-crossing objects.

Modeling the impact bombardment also allows us to quantify spatial cratering asymmetries that are not accounted for in the traditional crater chronology method. The presence of such asymmetries would bias the ages based on crater densities according to the location of the geologic unit. The Moon is subject to nonuniform cratering, since it is not massive enough to gravitationally homogenize encounter trajectories. Latitudinal asymmetries are produced by the anisotropy of encounter inclinations [START_REF] Le Feuvre | Nonuniform cratering of the terrestrial planets[END_REF], whereas longitudinal asymmetries result from synchronous rotation since the relative encounter velocity, hence impact rate, are maximized at the apex of motion (e.g., [START_REF] Morota | Influence of the asymmetrical cratering rate on the lunar cratering chronology[END_REF]. In addition, the Earth may focus low inclination and velocity projectiles onto the nearside of the satellite or act as a shield at small separation distances.

Several studies have attempted to give estimates of the lunar cratering asymmetries. [START_REF] Wiesel | The meteorite flux at the lunar surface[END_REF] used a simplified asteroid population, and [START_REF] Bandermann | Calculation of meteoroid impacts on moon and earth[END_REF] used analytical formulations based on strongly simplifying assumptions in order to calculate impact locations on a planet. These formulations did not allow to investigate any latitudinal effects. Wood (1973) numerically integrated the trajectories of ecliptic projectiles; [START_REF] Pinet | Lunar impact flux distribution and global asymmetry revisited[END_REF] numerically studied the asymmetries caused by geocentric projectiles that were potentially present early in the lunar history. [START_REF] Horedt | Cratering rate over the surface of a synchronous satellite[END_REF], [START_REF] Shoemaker | Cratering time scales for the Galilean satellites[END_REF], [START_REF] Zahnle | Cratering rates on the Galilean satellites[END_REF] and [START_REF] Zahnle | Differential cratering of synchronously rotating satellites by ecliptic comets[END_REF] all proposed analytical formulations for the apex/antapex effect, but the range of predicted amplitudes is very large (the first, in particular, claimed that this effect is negligible for the Moon). Moreover, these four studies based their results on isotropic encounter inclinations. [START_REF] Gallant | Current bombardment of the earthmoon system: Emphasis on cratering asymmetries[END_REF] numerically modeled projectile trajectories in the Sun-Earth-Moon system, and reported a significant apex/antapex effect. In our study, we derive a semi-analytic approach for calculating cratering rates of synchronously locked satellites, and apply it to the Moon. This method yields results that are nearly identical to full numerical simulations, is computationally very rapid, and easily generalizable to other satellites.

The "asteroid to crater" modeling requires several steps to create synthetic crater distributions. First, one needs to know the impactor population as a function of orbital elements, size and time. Second, impact probabilities are calculated over precession and revolution cycles of both the projectile and target, based on geometrical considerations [START_REF] Öpik | Collision probability with the planets and the distribution of planetary matter[END_REF][START_REF] Wetherill | Collisions in the asteroid belt[END_REF][START_REF] Greenberg | Orbital interactions: a new geometrical formalism[END_REF][START_REF] Bottke | Asteroidal collision probabilities[END_REF]. As these impact probabilities predict typically that a given planetary crosser will collide with a given planet over timescales of about 10 Ga, whereas the typical lifetime of these objects is thought to be a thousand times lower [START_REF] Michel | Origin and dynamics of near Earth objects[END_REF] as a result of ejection from the solar system or collision with another body, the calculated bombardment must be seen as the product of a steady-state population, where vacant orbital niches are continuously reoccupied [START_REF] Ivanov | Exogenic dynamics, cratering and surface ages[END_REF]. Third, scaling laws derived from laboratory experiments and dimensional analyses are used to calculate the final crater size produced for given a impact condition (e.g., impactor size, velocity and cohesion of the target material), for which several different scaling laws have been proposed [START_REF] Schmidt | Some recent advances in the scaling of impact and explosion cratering[END_REF][START_REF] Holsapple | Point source solutions and coupling parameters in cratering mechanics[END_REF][START_REF] Gault | Impact cratering[END_REF][START_REF] Holsapple | The scaling of impact processes in planetary sciences[END_REF][START_REF] Holsapple | A crater and its ejecta: An interpretation of Deep Impact[END_REF]. Subsequent gravitational modification of the transient cavity that gives rise to the final crater size have been deduced from crater morphological studies [START_REF] Pike | Formation of complex impact craters -Evidence from Mars and other planets[END_REF][START_REF] Croft | The scaling of complex craters[END_REF].

The rest of this paper is organized as follows. In the following section, we describe the employed orbital and size distribution of the impactor population.

In section three, the necessary equations used in generating synthetic cratering rates are summarized; the full derivations are given in the appendix. In section four, we compare our synthetic lunar crater distribution with observations, we describe the lunar spatial asymmetries, and give the predicted planet/Moon cratering ratios. We also provide simple analytic equations that reproduce the predicted cratering asymmetries on the Moon, as well as the latitudinal variations predicted for the terrestrial planets in Le Feuvre and [START_REF] Le Feuvre | Nonuniform cratering of the terrestrial planets[END_REF]. Section five is dedicated to revising the crater chronology method.

New age estimates are proposed for geologic units on the Moon, Earth, Venus and Mercury.

Impactor population

Let us denote a, e, i and d respectively the semi-major axis, eccentricity, inclination and diameter of those objects whose orbits can intersect the inner planets. The entire population can be written

n(> d, a, e, i) = n(> 1) × o(a, e, i) × s(> d) , (1)
where n can be expressed as the product of n(d > 1), the total number of objects with a diameter greater than 1 km; o(a, e, i), the relative number of objects with a given set of orbital elements, normalized so that o(a, e, i) da de di = 1; and s(> d), the normalized cumulative number of objects larger than a given diameter, such that s(> 1 km) = 1. This formulation assumes that no correlations exist between the size of the object and its orbit, which is consistent with the observations of [START_REF] Stuart | Bias-corrected population, size distribution, and impact hazard for the near-Earth objects[END_REF] for diameters ranging from ∼10 m to ∼10 km.

Orbital distribution

The orbital distribution of near Earth-objects (NEOs) is taken from the model of [START_REF] Bottke | Debiased orbital and absolute magnitude distribution of the Near-Earth Objects[END_REF], which provides a debiased estimate of the orbital distribution of bolides that can potentially encounter the Earth-Moon system.

ously replenished by the influx coming from source regions associated with the main asteroid belt or the transneptunian disk. The model was determined through extensive numerical integrations of test particles from these sources, and calibrated with the real population observed by the Spacewatch survey.

The relative number of objects is discretized in orbital cells that spans 0.25 AU in semi-major axis, 0.1 in eccentricity, and 5 • in inclination.

In order to model the martian impact flux, we have amended this model by including the known asteroids that cross the orbit of Mars, but which are not part of the NEO population. These are taken from the database provided by E.

Bowell (Lowell Observatory), that gives the orbit and absolute magnitude H of telescopic discoveries. We consider the population of H < 15 objects to be the best compromise between a sufficient number of objects (in order to avoid sparseness in the orbital space), and completeness. Their distribution as a function of perihelion is very similar in shape to brighter (hence larger) H < 13 objects (see Le Feuvre and [START_REF] Le Feuvre | Nonuniform cratering of the terrestrial planets[END_REF], the latter being large enough to ensure they do not suffer observational bias in the martian neighborhood.

The H < 15 population is therefore considered as complete, and is scaled to match the modeled NEO population. From the relationship between absolute magnitude and geometric albedo of [START_REF] Bowell | Application of photometric models to asteroids[END_REF], and a mean albedo of 0.13 [START_REF] Stuart | Bias-corrected population, size distribution, and impact hazard for the near-Earth objects[END_REF], H < 15 corresponds to d >∼ 4 km.

Size distribution

The size distribution of impactors is taken from a compilation of various estimates of the size-dependent impact rate on Earth, as shown in figure 1. This compilation gathers atmospheric recordings of meteoroids and impact probabilities calculated with Öpik equations from debiased telescopic observations.

Estimates from telescopic observations are those of [START_REF] Rabinowitz | A reduced estimate of the number of kilometre-sized near-Earth asteroids[END_REF], [START_REF] Morbidelli | From magnitudes to diameters: the albedo distribution of Near Earth Objects and the Earth collision hazard[END_REF], [START_REF] Harris | A New Estimate of the Population of Small NEAs[END_REF] and [START_REF] Stuart | Bias-corrected population, size distribution, and impact hazard for the near-Earth objects[END_REF].

Estimates based on atmospheric recordings include [START_REF] Halliday | Detailed data for 259 fireballs from the Canadian camera network and inferences concerning the influx of large meteoroids[END_REF], ReVelle (2001) and [START_REF] Brown | The flux of small near-Earth objects colliding with the Earth[END_REF]. Concerning the LINEAR survey, estimates of [START_REF] Harris | A New Estimate of the Population of Small NEAs[END_REF] have been scaled at large sizes to the estimates of [START_REF] Stuart | Bias-corrected population, size distribution, and impact hazard for the near-Earth objects[END_REF]. In contrast to [START_REF] Morbidelli | From magnitudes to diameters: the albedo distribution of Near Earth Objects and the Earth collision hazard[END_REF] and [START_REF] Stuart | Bias-corrected population, size distribution, and impact hazard for the near-Earth objects[END_REF], [START_REF] Rabinowitz | A reduced estimate of the number of kilometre-sized near-Earth asteroids[END_REF] did not use a debiased albedo distribution, but rather a constant mean geometric albedo of 0.1 in order to convert their telescopic observations from magnitude to diameter at all sizes.

Though small objects are expected to possess a larger albedo as they are generally younger than large objects, we did not attempt to correct for this effect, since estimates of [START_REF] Rabinowitz | A reduced estimate of the number of kilometre-sized near-Earth asteroids[END_REF] show a general agreement with other studies. In atmospheric flash estimates, masses have been deduced from kinetic energy, the latter being estimated from luminous energy. We convert kinetic energies of [START_REF] Brown | The flux of small near-Earth objects colliding with the Earth[END_REF] into diameters using a mean density of 2700 kg m -3 and a mean impact velocity of 20 km s -1 [START_REF] Stuart | Bias-corrected population, size distribution, and impact hazard for the near-Earth objects[END_REF]. In order to increase the range of sizes in our compilation, and to reduce the statistical uncertainties associated with the larger objects, we have also included the size-frequency distribution of Mars-crossing objects with sizes greater than 4 km, and scaled these to the terrestrial impact rates of [START_REF] Stuart | Bias-corrected population, size distribution, and impact hazard for the near-Earth objects[END_REF].

Various assumptions have led to all these estimates. Among them, the assumed impact velocity and bolide density are only of moderate influence. As an example, varying the density from 2700 to 2000 kg m -3 , or the mean impact velocity from 20 to 17 km s -1 , changes the estimates of [START_REF] Brown | The flux of small near-Earth objects colliding with the Earth[END_REF] only by about 10 %. Of major influence are the luminous efficiency used to obtain kinetic energy from flashes (see [START_REF] Ortiz | Detection of sporadic impact flashes on the Moon: Implications for the luminous efficiency of hypervelocity impacts and derived terrestrial impact rates[END_REF], and the debiasing process in the case of telescopic observations, which are difficult to assign a statistical uncertainty to. Consequently, we simply fit a 10th-order polynomial to the entire dataset, assuming each data is error free, and that the average combination of all estimates gives a good picture of the impactor population.

[Table 1 about here.]

We express the resulting analytic size-frequency distribution as the product of two terms: the normalized size distribution log s(> d) = 10 i=0 s i (log d km ) i , whose coefficients are listed in table 1, and the Earth's impact rate for objects larger than 1 km, φe (d > 1) = 3.1 × 10 -6 Ga -1 km -2 . An overbar is appended to the Earth's impact rate symbol, denoting that this quantity is spatially averaged over the planet's surface. The size-frequency distribution of impactors is here assumed to be the same for all bodies in the inner solar system. The relative impact rates for these bodies with respect to Earth are calculated in Section 3 using the orbital distribution of the planet crossing objects.

[Fig. 1 about here.]

3 From asteroids to craters

Here we describe how is calculated the number of craters that form on the Moon and planets, per unit area and unit time, as a function of the crater diameter and location on the surface. We first need to calculate the encounter conditions generated by the impactor population, then the corresponding im-pact rate and conditions (impact velocity and incidence angle), in order to finally obtain a cratering rate by the use of impact crater scaling laws. For the reader's convenience, derivations are given in the appendix.

Encounter probability

Following Öpik (1951), the assumptions under which an encounter is considered to occur can be summarized as follows:

• An encounter between the target (Moon or planet) and impactor occurs at the geometrical point of crossing of the two orbits (the mutual node).

The geometry of encounter is given by the relative velocity vector U at this point, which is expressed here in a frame where the X-axis points towards the central body (planet or Sun), (XY) defines the target's orbital plane, and the Z-axis points upward.

• The relative encounter velocity does not account for the acceleration generated by the mass of the target.

• The impactor, as seen by the target, is treated as if it were approaching from an infinite distance, under only the gravitational influence of the target.

The encounter trajectory is therefore hyperbolic in the reference frame of the target.

• For a given U, there are an infinite number of hyperbolic trajectories that can actually strike the target, that are distributed uniformly on a circle perpendicular to U with a surface equal to the target's gravitational cross section.

These approximations hold as long as the radius of the target's Hill sphere is large enough with respect to the target size. We performed three-body numerical simulations that show that a factor of ten between the Hill sphere and target radii suffices to ensure the validity of the above approximations.

For the terrestrial planets, this condition is largely verified. For the Moon, it corresponds to a minimum Earth-Moon separation of ∼ 17 Earth's radii.

In the case of planets in nearly circular orbits, the encounter geometry U and probability P (providing the two orbits intersect) are simply given by the well known Öpik equations ( Öpik, 1951). In order to account for the eccentricity of the target (which is important for Mercury and Mars, but not for the Earth), we use the improved formulation of [START_REF] Greenberg | Orbital interactions: a new geometrical formalism[END_REF] and [START_REF] Bottke | Asteroidal collision probabilities[END_REF]. The probability is largest for low inclination encounters, and for encounters occurring near the projectile's pericenter and apocenter.

Singularities of the encounter probability are avoided following [START_REF] Dones | Dynamical lifetimes and final fates of small bodies: orbit integrations vs Öpik calculations[END_REF]. For a given orbital geometry, the encounter probability is proportional to the gravitational cross section, whose radius is

τ = R √ 1 + Γ , ( 2
)
where R is the target's radius and

Γ = 2 GM RU 2 (3)
is the Safronov parameter, with G the gravitational constant and M the target's mass.

Note that the calculated impact probabilities are long-term averages over precession cycles of both the projectile and target (i.e., the longitude of node and argument of pericenter can take any value between 0 and 2π). We account for secular variations of the planetary orbital elements using the probability distributions given as a function of time in [START_REF] Laskar | Chaotic diffusion in the solar system[END_REF]. Our results are only sensitive to secular variations for Mars, and in the following, two values are quoted for this planet, that correspond to the present day value and to a long-term average (1 Ga and 4 Ga averages yield nearly identical results).

By calculating the encounter probability P and velocity U associated with a given orbital element set (a, e, i), by weighting this probability with the relative number of objects o(a, e, i), and by summing over the entire planet-crossing population for each terrestrial planet, we build the probability distribution of the encounter conditions, p(U), using bins of 1 km s -1 for each component of the encounter velocity. Analytically, we have

p(U) = p (U) U p (U) dU , (4)
with

p (U) = D P (D) o(a, e, i) δ(U(D) -U ) dD , (5)
where the integration is performed over the 8-dimensional domain D = (U X , U Y , U Z , a, e, i, w 0 , ∆Ω), with w 0 the target's argument of perihelion, ∆Ω the difference between the target and projectile's longitudes of the ascending node (see [START_REF] Greenberg | Orbital interactions: a new geometrical formalism[END_REF], and δ the Kronecker function which equals 1 when U = U and 0 otherwise. The impact rate relative to the Earth is

r = U p (U) dU U p e (U) dU . ( 6
)
where p e is calculated from equation 5 with the Earth being the target.

The lunar case requires a specific treatment, which is detailed in appendix A.2.

For simplicity and without altering the results, it is assumed that the lunar orbit is circular about the Earth and possess a zero inclination with respect to the ecliptic. We first calculate encounters probabilities P and velocities

V with the entire Earth-Moon system, whose expression is the same as for the Earth, except that the gravitational cross section radius is replaced in the Öpik equation by what we call here the lunar orbit cross section, defined as [START_REF] Zahnle | Cratering rates on the Galilean satellites[END_REF])

τ = a m 1 + 2 v 2 m V 2 , ( 7
)
where a m and v m are respectively the lunar semi-major axis and velocity. On the cross sectional disk, the distance from the center is denoted by the impact parameter b. Only when b ≤ τ is it possible for a hyperbolic orbit to impact the Moon (note that the gravitational cross section of the Moon itself is accounted for later in the calculation, see appendix A.2). The probability distribution of encounter conditions p(V) is first calculated according to this new definition of the encounter probability. Then, the relative lunar encounter velocity U and probability P m are derived analytically for each hyperbolic orbit crossing the Earth-Moon system (equations A.19 -A.18, A.24 -A.25, A.30 -A.33).

The probability distribution of lunar encounter conditions, p(U), is then determined by integrating numerically P m for all possible hyperbolic orbits of each encounter V, and for all encounter velocities. Similarly, the impact probability with Earth P e is determined for each hyperbola, and the Moon/Earth impact ratio r is calculated over all possible encounters. Mathematically, p(U)

and r are given by equations (A.34-A.36).

Note that there is a dependence of the lunar impact rate on the Earth-Moon separation distance, a m . This distance has evolved outward with time, and we test various separation distances in the simulations. A major difference between our approach and previous investigations [START_REF] Shoemaker | Cratering time scales for the Galilean satellites[END_REF][START_REF] Zahnle | Cratering rates on the Galilean satellites[END_REF][START_REF] Zahnle | Differential cratering of synchronously rotating satellites by ecliptic comets[END_REF] is that the argument of pericenter of the hyperbolic orbits is not assumed to precess uniformly within the Earth-Moon system, but is explicitly given by the encounter geometry. Our formulation allows to calculate explicitly lateral asymmetries in the lunar cratering rate.

Impact rate

Let us express the cumulative impact flux, that is the number of objects with diameters greater than d that hit the planet per unit time and area, as

φ(> d, λ, ϕ) = φ(> d) × ∆φ(λ, ϕ) , (8)
where λ and ϕ are respectively the latitude and longitude, φ(> d) is the spatially-averaged impact rate for projectiles larger than d, and ∆φ(λ, ϕ) is the relative impact rate as a function of position, normalized to the global average.

Using our normalized impactor SFD, the average impact rate expresses as

φ(> d) = φ(> 1) s(> d) , (9)
and φ(d > 1) is obtained from the impact ratio between the target and Earth, r, and the terrestrial impact rate φe (d > 1) as φ(> 1) = r φe (> 1) .

The net spatial asymmetry ∆φ(λ, ϕ) is found by integrating the spatial assymetries δφ(λ, ϕ, U) associated with each encounter geometry:

∆φ(λ, ϕ) = U δφ(λ, ϕ, U) p(U) dU , ( 11
)
where δφ is given in appendix A.3 as a function of the Safronov parameter Γ and obliquity of the target (equations A.47 -A.49). The impact flux is homogeneous for Γ = ∞, that is, for encounter velocities negligible with respect to the target's surface gravitational potential (equation 3). On the other hand, for Γ = 0, encounter trajectories are straight lines, and the impact flux is a simple geometrical projection of the spatially uniform encounter flux (on a plane perpendicular to the radiant) onto the target's spherical surface. The associated impact velocities and incidence angles are calculated from equations (A.51) and (A.54). The impact velocity u is only dependent on U and Γ, while the incidence angle θ further depends on position.

Cratering rate

To obtain cratering rates from impact rates, we need to convert the impactor diameters into crater diameters. For this purpose, we use equations that have been derived in the framework of π-scaling dimensional analysis [START_REF] Holsapple | Point source solutions and coupling parameters in cratering mechanics[END_REF][START_REF] Holsapple | The scaling of impact processes in planetary sciences[END_REF], where the transient crater is given as a function of the projectile diameter, impact velocity, surface gravity and projectile/target density ratio (see appendix A.4). It is assumed that only the vertical component of the impact velocity, whose value is obtained from the impact angle, contributes to the crater size [START_REF] Pierazzo | A reevaluation of impact melt production[END_REF], though other relations could be easily incorporated into this analysis. The scaling equation and parameters are taken from the summary of [START_REF] Holsapple | A crater and its ejecta: An interpretation of Deep Impact[END_REF] for the case of porous and non-porous scaling. It will be seen that both formation regimes are necessary to reconcile the impactor and crater SFDs.

We only consider craters that form in the gravity regime, where the tensile strength of rock is negligible, that is, craters larger than a few hundred meters in competent rock, and larger than a few meters in consolidated soils. An increase of the transient crater diameter by wall slumping and rim formation is accounted for as given in Melosh (1989, 253 pp.). Finally, large craters collapse due to gravity, becoming complex craters, and the relationship between simple and complex crater diameters is taken from [START_REF] Holsapple | The scaling of impact processes in planetary sciences[END_REF] The cratering rate, that is the number of craters larger than D that form at (λ, ϕ) per unit time and area, is

C(> D, λ, ϕ) = U φ(> d, λ, φ) p(U) dU , (12)
where

d = d(D, λ, φ, U) . ( 13
)
For convenience, we separate the cratering rate into

C(> D, λ, ϕ) = C(> D) × ∆C(> D, λ, ϕ) . ( 14
)
where C(> D) is the spatially averaged rate and ∆C(> D, λ, ϕ) is the relative spatial variation. Note that ∆C depends on D, though in practice, this dependence is moderate (see next section). 20

Crater size-frequency distributions

We first present our synthetic size-frequency distribution of lunar craters. Following the terminology of [START_REF] Marchi | A new chronology for the moon and mercury[END_REF], we refer to this as a model production function. We compare our model production function with the two standard measured production functions of Neukum (Neukum, 1983, 186 pp.; [START_REF] Neukum | Crater size distributions and impact probabilities on Earth from lunar, terrestrial-planet, and asteroid cratering data[END_REF] and Hartmann (Basaltic Volcanism Study Project, 1981;[START_REF] Hartmann | Martian cratering VI. Crater count isochrons and evidence for recent volcanism from Mars Global Surveyor[END_REF]. We note that the two are in good agreement over the crater diameter range from 300 m to 100 km, but differ between 2 and 20 km, with a maximum discrepancy of a factor 3 at 5 km.

Using the traditional non-porous scaling relations and a standard target density of 2800 kg m -3 , we calculate that 2.88 × 10 -11 craters larger than 1 km would be created each year on the lunar surface by the present impactor population. Using the time-dependence established by Neukum (1983, 186 pp.) that predicts a quasi-constant impact flux over the last ∼ 3 Ga, the Hartmann and Neukum production functions return respective values of 7.0 × 10 -13 and 8.2 × 10 -13 , which are about forty times lower, implying that the present flux must be considerably larger than the time averaged value.

However, by using the porous scaling law instead, in order to account for the presence of megaregolith on the lunar surface, our calculated spatially averaged lunar cratering rate is

Cm (D > 1) = 7.5 × 10 -13 yr -1 km -2 , (15)
a value in excellent agreement with the two measured production functions under the assumption of a constant impact flux over the last ∼ 3 Ga.

Let us now reconcile the entire shape of the measured production functions with the observed impactor population. As shown in figure 2, the two measured distributions are very well fitted by using the porous regime for small craters (D < 2 km), and the non-porous regime for larger craters (D > 20 km). We model a simple smooth transition between the two regimes by considering that the impactor size d required to create a crater of diameter D is a linear combination of the sizes required from the porous and non-porous scaling relations, the influence of each regime depending on the depth of material excavated by the crater. The depth of excavation z T is about 1/10 of the transient crater diameter D T , and does not seem to depend on target properties [START_REF] Melosh | Impact crater collapse[END_REF]. The impactor size is averaged over the depth of excavation:

d = 1 z T z T 0 d z (z) dz , (16)
where d z is the impactor size required by the material at depth z, given by the porous regime at the surface, by the non-porous regime at depths larger than a given "megaregolith thickness" T , and by a linear combination between z = 0 and z = T :

d z (z) = 1 T ((T -z) d p + z d np ) for z ≤ T , d z (z) = d np for z ≥ T , (17)
with d p and d np the impactor diameters respectively required from the porous and non-porous regimes. In the calculation of d np , the target density is set to 2800 kg m -3 (solid rocks), whereas we assume in calculating d p that the density of the porous material is 2500 kg m -3 , based on [START_REF] Bondarenko | Thickness of lunar regolith and Soderblom's crater parameter Dl[END_REF]Shkuratov (1999, abstract no. 1196) who inferred an upper regolith density comprised between 2300 and 2600 kg m -3 from correlations between the surface regolith thickness and the Soderblom's crater parameter [START_REF] Soderblom | Technique for rapid determination of relative ages of lunar areas from orbital photography[END_REF]. We note that given the simplicity of our crater-scaling procedure in the transition zone, the correspondance between T and the actual megaregolith thickness should not be expected to be exact.

As shown in figure 2 The use of porous scaling was first suggested by [START_REF] Ivanov | Earth/Moon impact rate comparison: searching constraints for lunar secondary/primary cratering proportion[END_REF] (see also [START_REF] Ivanov | Size-frequency distribution of asteroids and impact craters: estimates of impact rate[END_REF][START_REF] Ivanov | Exogenic dynamics, cratering and surface ages[END_REF], and is a natural consequence of a highly fractured megaregolith on airless bodies. Also natural is to expect that the thickness of the megaregolith will depend upon both age and local geology.

We note that the need for a transition regime falls within the diameter range where the measured production functions differ the most (excluding very large craters). We suspect that this is partially a result of the Hartmann production function being based on crater counts performed solely over mare units, whereas the Neukum production function also includes older highlands ter-rains (Neukum et al., 2001a).

Estimated megaregolith thicknesses are roughly consistent with seismic models of the lunar crust (e.g., [START_REF] Warren | Structure in the upper lunar crust[END_REF][START_REF] Lognonné | A new seismic model of the moon: implications for structure, thermal evolution and formation of the moon[END_REF] that generally predict reduced seismic velocities for the upper km, which is attributed to an increased porosity and fractures. Furthermore, seismic data at the Apollo 17 landing site, overlaid by mare basalt, indicates that the upper 250 / 400 m show a very low P-wave velocity with respect to the deeper basalt [START_REF] Kovach | Apollo 17 seismic profiling: Probing the lunar crust[END_REF][START_REF] Cooper | Lunar near-surface structure[END_REF], the lower estimate being in agreement with our calculated megaregolith depth of 250 m for the Hartmann production function. Finally, [START_REF] Thompson | Blocky craters: Implications about the lunar megaregolith[END_REF] show by analysis of radar and infrared data (which are dependent on the amount of near surface rocks) that craters overlying highlands show different signatures for craters greater and less than 12 km, and that mare craters down to 4 km in diameter possess a similar signature to that of highlands craters greater than 12 km. They attribute this difference to the presence of a pulverized megaregolith layer that is thicker in the older highlands than the younger mare.

By the use of a porous regime dictated by the properties of a megaregolith, our model production function reproduces the measured crater distributions in shape and in the absolute number of craters formed over the past 3 Ga, under the assumption of a constant impact flux. We caution that our simple formulation of the porous / non-porous transition does not account for the temporal evolution of the megaregolith and that the inferred megaregolith thicknesses are only qualitative estimates.

The present Earth-Moon distance has been used in the above calculation of the lunar cratering rate, and temporal variations in the lunar semimajor axis could, in principle, modify the Earth/Moon impact ratio and encounter velocity distribution with the Moon. Nevertheless, it is found in our simulations that, for a lunar semi-major axis as low as 20 Earth radii, the average lunar cratering rate is changed by less than 3%. This implies that both the shielding and gravitational focusing of projectiles by the Earth are of very moderate effects, especially since the Moon is believed to have spent the vast majority of its history beyond 40 Earth radii.

[Table 2 [START_REF] Le Feuvre | Nonuniform cratering of the terrestrial planets[END_REF][START_REF] Ivanov | Mars/Moon cratering rate ratio estimates[END_REF]. In addition to calculating the present day martian cratering rate, we also used the probability distribution of the martian eccentricity provided by [START_REF] Laskar | Chaotic diffusion in the solar system[END_REF] to calculate an average over the past 1 Ga (note that this value is nearly insensitive for longer averages).

Planetary size-frequency distributions are generally expressed with respect to the lunar one. This is done by defining the size-dependent quantity R c , which is the cratering ratio with respect to the Moon,

R c (> D) = C(> D) Cm (> D) . ( 18
)
For illustrative purpose, R c is shown for the inner planets in figure 3 by assuming that craters with diameters less than 10 km form in a porous soil on both the planet and Moon, while craters with greater sizes form in solid rocks (except for the Earth and Venus where only the non-porous regime is used). Note that R c can be easily calculated from equation 18 and table 2 for a different (and more realistic) transition between porous and non-porous regimes.

[Fig. 3 The mean impact velocity on the Moon is calculated to be ūm = 19.7 km s -1 . The full probability distribution of impact velocities for each planet is

shown in figure 4. The quantities φ/ φm , ū/ū m and g/g m , that give the relative impact flux, impact velocity and surface gravity with respect to the Moon, are given in table 3 for the inner planets. Mars experiences a high impact rate with respect to the Moon (about 3) due to its proximity to the main asteroid belt. In comparison, the martian cratering ratio is reduced (between about 0.5 and 2.5) because the impact velocity on Mars is significantly lower than on the Moon, requiring larger (and hence less numerous) impactors to create a crater of a given size. Mercury exhibits also a high impact rate, and the impact velocity is about twice as large as on the Moon, resulting in a high value of the cratering ratio R c , comprised between 2 and 4. The impact rate is larger on the Earth and Venus than on the Moon, as these planets possess a higher gravitational attraction. Their higher surface gravities compensate the differences in impact velocities with the Moon, and R c is comprised between 0.5 and 1.5 for the Earth and between 1 and 2 for Venus.

[Table 3 about here.]

Spatial variations

The relative spatial cratering variations on the Moon, ∆C(> D, λ, ϕ), are shown in figure 5 for the present Earth-Moon distance of about 60 Earth radii, and for crater diameters larger than 1 km. The cratering rate varies from approximately -20% to +25% with respect to the global average. It is minimized at about (±65 • N, 90 • E), whereas the maximum, which is a factor 1.5 higher, is located at the apex of motion (0 • N, 90 • W ).

[Fig. 5 about here.]

Two effects conjugate to give such a distribution. First, a latitudinal effect, detailed in Le Feuvre and [START_REF] Le Feuvre | Nonuniform cratering of the terrestrial planets[END_REF], comes from the higher proportion of low inclination asteroids associated with the higher probabilities of low inclination encounters. The pole/equator ratio is 0.80. Second, a longitudinal apex/antapex effect comes from the synchronous rotation of the Moon and the higher relative encounter velocities at the apex. The lunar orbital velocity is added to the projectile velocity for impacts at the apex, whereas it is subtracted at the antapex. The apex/antapex ratio is 1.37. We note that there is a negligible nearside/farside effect: the nearside experiences the formation of about 0.1% more craters that the farside. The Earth does indeed concentrate very low inclination (and moderate velocity) projectiles onto the lunar nearside, but these are not numerous enough to influence the global distribution.

The lateral cratering variations depend on the crater size since the sizefrequency distribution of impactors s(> d) is not a simple power law, and the impact conditions are not everywhere the same. Nevertheless, the maximum/minimum cratering ratio varies only by about 5% for D ranging from 30 m to 300 km. Consequently, for the following discussion, we shall consider that ∆C(> D, λ, ϕ) ∆C(> 1, λ, ϕ).

For smaller Earth-Moon distances a m , the apex/antapex effect increases as the lunar velocity increases. Between 20 and 60 Earth radii, this dependency is fit by the simple equation

C(apex)/C(antapex) = 1.12 e -0.0529 am Re + 1.32 , (19)
where R e is the Earth radius. Over this range of Earth-Moon separations, latitudinal variations and the nearside/farside effect are found to vary by less than 1%. These calculations assume that the lunar obliquity stayed equal to its present value in the past.

Figure 6 plots the relative cratering rate ∆C as a function of the angular distance from the apex of motion for the present Earth-Moon distance, and compares this with the counts of rayed craters with diameters greater than 5 km given in [START_REF] Morota | Influence of the asymmetrical cratering rate on the lunar cratering chronology[END_REF]. Rayed craters are younger than about 1 Ga [START_REF] Wilhelms | Bombardment as a cause of the lunar asymmetry[END_REF], which should corresponds to an Earth-Moon separation distance very close to the present one [START_REF] Sonett | Neoproterozoic Earth-Moon dynamics: Rework of the 900 Ma Big Cottonwood Canyon tidal laminae[END_REF][START_REF] Eriksson | Quantifying the oldest tidal record: The 3.2 Ga Moodies Group, Barberton Greenstone Belt, South Africa[END_REF]. As is seen, the model compares favorably to the data.

[Fig. 6 about here.]

We note that the impact rate exhibits nearly the same behavior as the cratering rate, but with a reduced amplitude. The pole/equator ratio is 0.90, whereas the apex/antapex ratio is 1.29. The latitudinal cratering variations are enhanced with respect to the impact rate variations, as the mean impact angle and impact velocity are smaller at the poles than at the equator (respectively by 2.5 degrees and 500 m/s), requiring a larger projectile to create the same crater size. As large projectiles are less numerous than small ones, the impact rate at the poles is smaller than the cratering rate (Le Feuvre and [START_REF] Le Feuvre | Nonuniform cratering of the terrestrial planets[END_REF]. This is also true for the apex/antapex asymmetry, as the average impact velocity is 500 m/s higher at the apex than at the antapex.

However, the increase is moderate, as the mean impact angle is only about 1.5 degree smaller.

It is seen in figure 7 that our predicted apex/antapex effect differs from that of [START_REF] Zahnle | Differential cratering of synchronously rotating satellites by ecliptic comets[END_REF]. These authors describe their variations of the impact rate as a function of γ, the angular distance to the apex, as

∆φ(γ) =   1 + v m 2v 2 m + V 2 cos γ   2 , (20)
where V 19 km s -1 is the mean encounter velocity with the Earth-Moon system. We are able to reproduce their analytical solution, but only under the condition where we force the encounter inclinations with respect to the lunar orbit plane to be isotropic in space. These authors used Öpik equations [START_REF] Shoemaker | Cratering time scales for the Galilean satellites[END_REF] for hyperbolic orbits that were assumed to precess uniformly inside the planet-moon system. We nevertheless point out that approximation might be valid.

[Fig. 7 about here.]

We next provide analytical solutions for the relative variations of both the impact and cratering rates on the Moon. We also give solutions for the latitudinal variations presented on the terrestrial planets in Le Feuvre and [START_REF] Le Feuvre | Nonuniform cratering of the terrestrial planets[END_REF]. Two values are quoted for Mars, one that corresponds to its present obliquity and eccentricity and the other to averaged results using variations over 3 Ga as given in [START_REF] Laskar | Long term evolution and chaotic diffusion of the insolation quantities of Mars[END_REF] (see Le Feuvre and [START_REF] Le Feuvre | Nonuniform cratering of the terrestrial planets[END_REF].

Spatial variations in the impact flux and cratering rate are parameterized by a sum a spherical harmonic functions

∆C(λ, ϕ) = ∞ l=0 l m=-l C lm Y lm (λ, ϕ) , (21)
where Y lm is the spherical harmonic function of degree l and order m, C lm is the corresponding expansion coefficient, and (λ, ϕ) represents position on the sphere in terms of latitude λ and longitude ϕ, respectively. The real spherical harmonics are defined as

Y lm (λ, ϕ) =                P lm (sin λ) cos mϕ if m ≥ 0 P l|m| (sin λ) sin |m|ϕ if m < 0, (22)
and the corresponding unnormalized Legendre functions are listed in Table 4.

Many of the expansion coefficients are nearly zero since the cratering rate is symmetric about both the equator and the axis connecting the apex and antapex of motion. Only the most significant coefficients are listed in Tables 5 and6, which for most cases reproduce the data to better than 0.2%. The coefficients for the latitudinal variations on the terrestrial planets are listed in Tables 7 and8.

[ Neukum, 1983, 186 pp.;Strom and Neukum, 1988;[START_REF] Neukum | Crater size distributions and impact probabilities on Earth from lunar, terrestrial-planet, and asteroid cratering data[END_REF] established empirically the following relationship between the number of craters with diameters greater than 1 km and the age of the geologic unit:

N (D > 1, t) = a e bt -1 + c t , (23)
where N (D > 1, t) is given per 10 6 km 2 , t is the age expressed in Ga, and a = 5.44 × 10 -14 , b = 6.93 and c = 8.38 × 10 -4 . This relationship is essentially linear over the last 3.3 Ga (constant cratering rate in time) and approximately exponential beyond. The data used to construct this empirical curve are obtained from radiometric ages of the Apollo and Luna rock samples, compared to the crater density covering the associated geologic unit. We emphasize that no agreed upon calibration data exist between 1 and 3 Ga and beyond 3.9

Ga [START_REF] Stöffler | Stratigraphy and isotope ages of lunar geologic units: Chronological standard for the inner Solar System[END_REF]. We also note that equation ( 23) was originally obtained using age estimates of the highland crust and Nectaris impact basin, both of which are disputed and considerably older than 3.9 Ga.

Accounting for our calculated spatial variations, we first convert the measured crater density at a given site, N (D > 1, λ, ϕ), into the corresponding spatially averaged quantity:

N (D > 1) = N (D > 1, λ, ϕ)/∆C(D > 1, λ, ϕ) . ( 24
)
The cratering asymmetry ∆C is given in table 2, and we use the function that corresponds to the present Earth-Moon separation for cratered surfaces that are less than 1 Ga (consistent with the tidal deposit data of [START_REF] Sonett | Neoproterozoic Earth-Moon dynamics: Rework of the 900 Ma Big Cottonwood Canyon tidal laminae[END_REF]). We choose a Earth-Moon separation of 40 Earth radii for units that are older than 3 Ga (based on the tidal deposits data of [START_REF] Eriksson | Quantifying the oldest tidal record: The 3.2 Ga Moodies Group, Barberton Greenstone Belt, South Africa[END_REF]). This lunar semimajor axis value corresponds to a lunar orbital velocity twice as large as the present one. We further assume that the lunar obliquity was equal to its present value (nearly zero) for the entire time between 3.9 Ga and the present.

The data used to calibrate the crater density versus age relationship are listed in table 9, along with their corrections accounting for spatial variations in the cratering rate. We use the crater density and ages values quoted in [START_REF] Stöffler | Stratigraphy and isotope ages of lunar geologic units: Chronological standard for the inner Solar System[END_REF]. We did not attempt to fit crater distributions with our model production function to re-estimate N (D > 1) values, since our model already reproduces very well the Neukum production function that was used to estimate this quantity. For the case of very young calibration surfaces, we suspect that a thinner megaregolith might change the crater distribution with respect to the Neukum production function at sizes larger than one kilometer, but, for small exposure times, the largest observed craters are below this diameter value. Based on the data and interpretations of [START_REF] Norman | The lunar cataclysm: Reality or 'mythconception' ?[END_REF], we exclude the Nectaris basin from consideration, as it is possible that samples assigned to this basin have instead an Imbrium provenance. We thus assume that the Descartes formation is not Nectaris ejecta, but rather Imbrium ejecta with an age of 3.85 Ga. We also exclude the Crisium basin from the fit, due to the uncertain provenance of samples dated at the Luna 20 site. Finally, we use the recent crater counts performed on Copernicus deposits by Hiesinger et al.

(2010, abstract no 1533) with high-resolution Lunar Reconnaissance Orbiter images, which agree with a constant impact flux during the last 800 Ma, in contrast to previous studies.

As seen in table 9, the spatial correction is in general moderate, since the calibration terrains are located in the central portion of the nearside hemisphere, far from the extrema of the spatially-dependent cratering rate. Nevertheless, we point out that after correction, the Apennines, Fra Mauro and Descartes formations, which are all Imbrian in age, exhibit nearly the same globally averaged crater density, which is consistent with our assignment of an Imbrium origin to the Descartes formation, as recently suggested by [START_REF] Norman | The lunar cataclysm: Reality or 'mythconception' ?[END_REF] (see also [START_REF] Haskin | The case for an imbrium origin of the apollo th-rich impact-melt breccias[END_REF].

We perform a new fit of the calibration point, using the same functional form as in equation 23. Our resulting best parameters are

a = 1.893 × 10 -26 b = 14.44 c = 7.960 × 10 -4 , (25)
and our proposed curve is shown in figure 8, along with the curves proposed by Neukum et al. (2001a) and [START_REF] Marchi | A new chronology for the moon and mercury[END_REF]. It is seen that our model favors a longer period of nearly constant impact flux (up to about 3.5 Ga), and that the agreement with the calibration points is improved. Crisium basin is also plotted (but was not included in the fit) for the tentative age of 3.895 given by [START_REF] Swindle | Searching for crisium basin ejecta -chemistry and ages of luna 20 impact melts[END_REF]Swindle et al. ( , abstract no 1229)), showing an excellent agreement with our model.

We note that the only calibration point that is not well fit to our relationship is the "young" group of basalts at the Apollo 11 landing site (point 9 in figure 8). As discussed in [START_REF] Stöffler | Stratigraphy and isotope ages of lunar geologic units: Chronological standard for the inner Solar System[END_REF], four different flow units having distinct compositions and ages were sampled at this site, with ages of 3.58, 3.7, 3.8, and 3.85 Ga. Previously measured crater densities were assigned to the 3.58 and 3.8 Ga age groups, and we acknowledge the possibility that one, or perhaps both, of these assignments could possibly be in error. As an example, we note that if the 3.7 Ga age were assigned to the younger crater density that it would lie directly on our best fit line. In contrast, if the older crater density were assigned the 3.7 Ga age, this would require a different form to the fit between 3.41 and 3.75 Ga. Given the discrepancy of the 3.58 Ga crater density, ages derived from our crater chronology between 3.41 and 3.75 Ga should be used with caution. [Table 9 about here.]

The number of craters larger than D that accumulate on a given planet, at the location (λ, ϕ), and over a period of t Ga, is

N (> D, λ, ϕ, t) = C(> D, λ, ϕ) × T (t) , (26)
where the temporal dependency is

T (t) = a e b t -1 + c t a (e b -1) + c , (27)
and a , b and c are given in equation 25. The dating of a planetary geologic unit is performed by finding the value of t that allows the best fit between N (D, λ, ϕ, t) and the data.

At a given location (λ, ϕ) (note that ϕ plays a role only in the lunar case and that, apart from Mars where the cratering rate is expected to be 30 % less at the equator than at the pole, the latitudinal effect on the terrestrial planets is less than 10%), the best age is found by solving the linear least-square problem as a function of T (t):

N obs (> D k ) = C(> D k , λ, ϕ) × T (t),
T (t) = K k=1 C(> D k , λ, ϕ) N obs (> D k )/σ 2 k K k=1 C 2 (> D k , λ, ϕ)/σ 2 k (28)
The a posteriori uncertainty is given by

δT (t) = K k=1 C 2 (> D k , λ, ϕ) σ 2 k -1 2 . ( 29
)
The best age t (in Ga) and its uncertainty are finally determined from T (t) ± δT (t).

We next show an application of our procedure by estimating new ages of key geologic units on the Moon, Earth, Venus and Mercury. Figure 9a shows our [START_REF] Hughes | A new approach to the calculation of the cratering rate of the earth over the last 125±20myr[END_REF] is for an age of 138 ± 6 Ma (figure 9b), in excellent agreement with a constant impact flux over this period. Craters smaller than 15 km in diameter were not included in the fit, being depleted as a consequence of erosional processes (note that the megaregolith thickness was obviously set to 0 in this calculation). Our result implies that craters larger than about 10-20 km have not been erased by erosional processes over the last 100 Ma. We also fit our model to the crater distribution of a number of North American and European Phanerozoic cratons that have been independently estimated by [START_REF] Grieve | The terrestrial cratering record. II -the crater production rate[END_REF] to be 375 Ma old on average. [START_REF] Neukum | Crater size distributions and impact probabilities on Earth from lunar, terrestrial-planet, and asteroid cratering data[END_REF] noted that the crater chronology, once applied to this crater population, returns an age of about 700 Ma. Our best fit corresponds to an age of 390 ± 75

Ma, in agreement with [START_REF] Grieve | The terrestrial cratering record. II -the crater production rate[END_REF].

Using crater counts derived from the Magellan mission, the average age of the Venusian surface has generally been estimated to be between 650 and 750 Ma [START_REF] Neukum | Crater size distributions and impact probabilities on Earth from lunar, terrestrial-planet, and asteroid cratering data[END_REF][START_REF] Mckinnon | Cratering on Venus: models and observations[END_REF][START_REF] Korycansky | Modeling crater populations on Venus and Titan[END_REF]. We fit our model distribution to the data, excluding craters smaller than 25 km in diameter: their depletion with respect to the model is attributed to atmospheric shielding. Our estimate, based on the use of the non-porous scaling alone, is shown in figure 9c, and corresponds to an age of 240 ± 13 Ma (figure 9c), significantly younger than other estimates. The difference lies principally in the impactor population used in the calculations.

In particular, both [START_REF] Korycansky | Modeling crater populations on Venus and Titan[END_REF] and [START_REF] Mckinnon | Cratering on Venus: models and observations[END_REF] used the Venus-crossing population estimated by [START_REF] Shoemaker | Asteroid and comet flux in the neighborhood of Earth[END_REF], which produces an impact rate that is significantly higher than our value calculated from the NEO model of [START_REF] Bottke | Debiased orbital and absolute magnitude distribution of the Near-Earth Objects[END_REF]. This high impact rate is compensated in their study by a very efficient atmospheric shielding model that has the effect of decreasing the proportion of small craters. In particular, the number of 30 km craters is reduced by about a factor 5 between their airless and atmosphere-shielded simulations in order to fit the data, whereas our airless model fits the crater distribution at that size, suggesting that the atmosphere is of negligible effect. [START_REF] Neukum | Crater size distributions and impact probabilities on Earth from lunar, terrestrial-planet, and asteroid cratering data[END_REF] assumed that the mean impact velocity on Venus is 19 km s -1 (we calculate 25 km s -1 ), and that the projectile flux at infinity was the same for Venus and the Moon. Moreover, our venusian size-frequency distribution is not of the same shape as in [START_REF] Neukum | Crater size distributions and impact probabilities on Earth from lunar, terrestrial-planet, and asteroid cratering data[END_REF], because it was not constructed from the lunar production function, where porous megaregolith decreases the proportion of small craters, but rather directly from the impactor size distribution. We finally note that our inferred age is more similar to the estimate of [START_REF] Strom | The cratering record on mercury and the origin of impacting objects[END_REF], who

calculated the global surface of Venus to be 290 Ma old.

[Fig. 9 about here.]

Figure 9d shows crater counts performed on Mercury, both on plains interior to the Caloris basin and on the lineated Caloris basin sculpture [START_REF] Fassett | Caloris impact basin: Exterior geomorphology, stratigraphy, morphometry, radial sculpture, and smooth plains deposits[END_REF]. Our best fits are respectively for megaregolith thicknesses of 250 m and 900 m, and for ages of 3.30 ± 0.3 and 3.73 ± 0.2 Ga. According to [START_REF] Fassett | Caloris impact basin: Exterior geomorphology, stratigraphy, morphometry, radial sculpture, and smooth plains deposits[END_REF], the size-frequency distribution of craters on the lineated sculpture (resulting from the basin formation) should be more representative of the Caloris basin formation time than the distribution measured on the Caloris rim, that exhibits a loss of small craters due to erosional processes and bad lighting conditions for crater counts, both being the consequences of the rim slope. Therefore, we attribute the age of 3.69 Ga to the Caloris basin. Our preferred values are consistent with the interpretation that older surfaces present a thicker layer of megaregolith than younger ones. Moreover, the megaregolith appears to be thicker on Mercury than on the Moon, which is consistent with the fact that Mercury has both impact rate and mean impact velocity that are about twice as great as for the Moon. Finally, we suggest the possibility that differences in crater size-frequency distribution between pre-and post-LHB surfaces may partially be explained by megaregolith thickness variations rather than changes in the impactor size distribution, the latter hypothesis being supported by [START_REF] Strom | The origin of planetary impactors in the inner solar system[END_REF].

Discussion

Here we compare our calculations with published results, and discuss the sensitivity of our model to various parameters. Our present Earth/Moon impact ratio and mean lunar impact velocity are calculated to be, respectively, 1.58 and 19.7 km s -1 . This is comparable to the values of [START_REF] Stuart | Bias-corrected population, size distribution, and impact hazard for the near-Earth objects[END_REF] (respectively 1.61 and 19.3 km s -1 ), who used Öpik calculations for an isolated

Moon possessing an Earth-like orbit about the Sun. These authors used the NEO population of [START_REF] Stuart | A Near-Earth Asteroid population estimate from the LINEAR survey[END_REF], estimated from the debiased observations of the LINEAR survey, and calculated the cratering rate on the Earth and Moon using various scaling laws. They did not account for a porous megaregolith. [START_REF] Gallant | Current bombardment of the earthmoon system: Emphasis on cratering asymmetries[END_REF] used numerical simulations to estimate lunar asymmetries. From the NEO population of [START_REF] Bottke | Debiased orbital and absolute magnitude distribution of the Near-Earth Objects[END_REF], they calculated the probability of encounter with the Earth-Moon system with Öpik probabilities, and used a 4-body (projectile, Moon, Earth and Sun) numerical approach to calculate the encounter trajectory, starting at 0.02 AU from the Earth. They reported a mean lunar impact velocity of 20 km s -1 . Our analytical approach yields consistent values with their predicted apex/antapex cratering ratio as a function of the Earth-Moon distance: our leading/trailing hemispheric ratio is 1.17 for the present lunar semi-major axis a m = 60 Earth radii, and 1.20

for a m = 40 Earth radii, while [START_REF] Gallant | Current bombardment of the earthmoon system: Emphasis on cratering asymmetries[END_REF] gave respectively 1.15 and 1.18. Their latitudinal effect is found to be 10%, whereas ours is 20%.

The difference is due to the crater scaling employed: [START_REF] Gallant | Current bombardment of the earthmoon system: Emphasis on cratering asymmetries[END_REF] did not include the effect of the impact angle on the crater diameter, while we consider that only the vertical component of the impact velocity contributes to the crater size. By ignoring the effect of impact angle, we indeed obtain a similar latitudinal effect of ∼ 10% (see Le Feuvre and [START_REF] Le Feuvre | Nonuniform cratering of the terrestrial planets[END_REF].

The "vertical component" scaling appears to be the safest assumption for a single target body, though the impact angle dependence of the average crater efficiency may vary from planet to planet (see Collins et al. 2009, Elbeshausen 2009). Nevertheless, this angle dependence is of negligible influence in the calculation of planetary cratering rate ratios with respect to the Moon, and only affects the amplitude of the cratering asymmetries on a given body. [START_REF] Gallant | Current bombardment of the earthmoon system: Emphasis on cratering asymmetries[END_REF] also estimated the terrestrial impact rate as a function of local time (see figure 9 of their paper). Although it is not the purpose of this paper, we note here that our semi-analytical approach reproduces their results.

Recently, [START_REF] Marchi | A new chronology for the moon and mercury[END_REF] proposed a revised crater chronology. The main differences with our approach (excluding the assumption of spatially uniform cratering rates in the latter) are the following:

• We use the orbital distribution of near-Earth objects of [START_REF] Bottke | Debiased orbital and absolute magnitude distribution of the Near-Earth Objects[END_REF], modified for Mars, which is assumed to be in steady state and independent of bolide size. [START_REF] Marchi | A new chronology for the moon and mercury[END_REF], in contrast, use a subset of this model based on orbital integrations of test particles coming from only the 3:1 and ν 6 resonances. Furthermore, the relative importance of these two sources is

• We assume that the size frequency distribution of objects impacting the planets is the same for all planets and that the probability of an object impacting a planet is independent of size. In contrast, in [START_REF] Marchi | A new chronology for the moon and mercury[END_REF], the size frequency distribution differs for each planet and is dependent on how much time the impactors spend in the asteroid belt. Furthermore, the impact probabilities are assumed to depend upon the projectile size.

• In the impact crater scaling relations, we use K=1.17 for impacts in water for the gravity regime, whereas [START_REF] Marchi | A new chronology for the moon and mercury[END_REF] use K=0.93 for impacts in wet soils and rock in the strength regime (though employed in the general scaling relationship that contains both strength and gravity terms).

• When converting transient crater diameters to final crater diameters, we use a multiplicative factor of 1.56 as suggested by Melosh (1989, 253 pp.) and [START_REF] Melosh | Crater scaling program[END_REF], whereas [START_REF] Marchi | A new chronology for the moon and mercury[END_REF] assume that the transient crater diameter is equivalent to the final simple crater diameter for their preferred impact scaling law that is based on the equations in [START_REF] Holsapple | A crater and its ejecta: An interpretation of Deep Impact[END_REF].

• Both studies treat the case of impact crater scaling in the porous megaregolith differently. We use a linear transition from porous (i.e., "sand or cohesive soils" from [START_REF] Holsapple | A crater and its ejecta: An interpretation of Deep Impact[END_REF]) to nonporous scaling to fit either the Neukum or Hartmann production functions. In contrast, [START_REF] Marchi | A new chronology for the moon and mercury[END_REF] place the transition between porous and nonporous scaling at projectile diameters of 0.5 km (corresponding to a crater diameter of about 5 km), and additionally use an average strength parameter in their impact crater scaling relationship that depends on crater size.

It is difficult to quantify how each of these differences affect the final crater size-frequency distribution on a planetary object, and hence the derived ages of a surface. Nonetheless, we note that the different bolide size-frequency distributions and the different crater scaling laws could be significant. As an example, if [START_REF] Marchi | A new chronology for the moon and mercury[END_REF] were to have used the same crater scaling constants and to have multiplied their transient crater diameters by the same factor as us, this would have increased their predicted crater densities by about a factor of 5. Excluding these factors from our analysis would have decreased our crater densities by the same factor.

Regardless of these differences, we note that our model production function is in closer agreement with the Neukum production function than that proposed by [START_REF] Marchi | A new chronology for the moon and mercury[END_REF]. In particular, our crater densities differ by a factor of 1.1 and 0.8 at crater diameters of 1 and 50 km, respectively (1.3 and 0.7 with respect to the Hartmann production function), whereas for Marchi et al.

(2009) the respective factors are about 1.5 and 0.3 (1.7 and 0.2 with respect to the Hartmann production function). Though a factor of 5 increase in crater densities would improve their fit to the Neukum production function at large crater diameters, the misfit would be significantly worse for small craters.

It should be noted that, according to [START_REF] Stewart | Martian subsurface properties and crater formation processes inferred from fresh impact crater geometries[END_REF], morphological analysis of young Martian craters combined with theoretical crater scaling relationships implies a factor of about 1.4 for modification of the transient cavity by wall slumping and rim uplift (in addition to gravity collapse of complex craters), whereas Melosh (1989, 253 pp.) suggested a factor of 1.56. Using the [START_REF] Stewart | Martian subsurface properties and crater formation processes inferred from fresh impact crater geometries[END_REF] factor would make our crater densities differ at most by 40 % from the Neukum production function for crater sizes between 100 m and 300 km. Using an impactor density of 2000 instead of 2700 kg m -3 would produce the same effect. Given the uncertainties in both the transient estimate, that is in closer agreement with the observed production functions.

The density and porosity of the lunar megaregolith are not known with any certainty. It is likely that porosity reduces the density of the first few hundred meters of the lunar crust. Nevertheless, the subsurface study of the Ries crater on Earth [START_REF] Pohl | The Ries impact crater[END_REF] indicates that the density reduction is moderate.

We used a value of 2500 kg m -3 for the upper density of the megaregolith in the calculation of the impact crater sizes, based on [START_REF] Bondarenko | Thickness of lunar regolith and Soderblom's crater parameter Dl[END_REF]Shkuratov (1999, abstract no. 1196). Regardless, we note that taking an extreme value of 2000 kg m -3 would still produce a broad agreement with the Neukum production function, with a maximum discrepancy of 50 % for crater sizes between 500 m and 300 km.

We test the sensitivity of our calculated lunar cratering asymmetries by using the NEO orbital distribution given by telescopic observations only, either for diameters larger than 1 km or 4 km, instead of the orbital model provided by [START_REF] Bottke | Debiased orbital and absolute magnitude distribution of the Near-Earth Objects[END_REF]. Using the d > 4 km population leads to a strong enhancement of lunar latitudinal asymmetries. The pole/equator is in this case equal to 0.56, due to the high proportion of these objects in the ecliptic plane. For the d > 1 km population, latitudinal variations are less than 1 % with respect to the average. These observations suggest that the observation of low inclination objects with d < 1 km is incomplete, which is not unexpected since this type of detection is more difficult than for high inclination objects.

Regardless, in both cases, the apex/antapex ratio does not change significantly.

The porous regime is required so that the calculated lunar production rate for craters smaller than about 10 km matches the Hartmann and Neukum production functions. As reported by [START_REF] Strom | The origin of planetary impactors in the inner solar system[END_REF], counts of post-LHB lunar craters seem to be adequately fit by the observed NEO size-frequency distribution and the sole use of the non-porous cratering regime. An alternative possibility is that young geologic units that are less impacted, possess a thinner porous megaregolith layer than older units, and that this has acted to change the apparent size-frequency distribution of lunar craters with time.

The evolution of megaregolith thickness with time remains to be determined, as well as the equivalent thickness on Mercury and Mars, for which the distribution of craters smaller than a few kilometers is still to be better understood before our model can be confidently used: Mercury may exhibit an abundance of large secondaries, whereas the presence of subsurface ice may significantly reduce the upper porosity of the martian crust.

Finally, we caution that our model may be inaccurate for periods prior to 3.5

Ga, both because of the Late Heavy Bombardment and a possible reorientation of the Moon [START_REF] Wieczorek | Did a large impact reorient the moon?[END_REF]. Our study is generalizable to different impactor populations or planetary reorientations, but the affects of these events are difficult to assess without a better knowledge. Waiting in particular for a suitable dynamical model for small bodies in the early solar system, we have applied our cratering model to times prior to 3.5 Ga, based on the following considerations. First, the size-frequency distributions of craters on the oldest lunar terranes are potentially explained by varying the megaregoltih thickness, without the need for a different impactor population.

Second, even if the cratering asymmetries would have been modified by these events, the influence of spatial variations on the age estimates is very moderate for times prior to 3.5 Ga (which correspond to the exponential part of the chronology curve).

We have presented a complete set of equations that allows to calculate the cratering rate as a function of crater size and location, on the Moon and inner planets. Spatial variations of the cratering rate are calculated semianalytically, hence considerably faster than by the use of a fully numerical method. The approach is generalizable to any other planet or moon, providing the orbital and size distributions of impactors are known. For our purpose, the NEO orbital distribution model of [START_REF] Bottke | Debiased orbital and absolute magnitude distribution of the Near-Earth Objects[END_REF] has been used, and the size distribution of bolides has been fitted to a compilation of various estimates [START_REF] Rabinowitz | A reduced estimate of the number of kilometre-sized near-Earth asteroids[END_REF][START_REF] Morbidelli | From magnitudes to diameters: the albedo distribution of Near Earth Objects and the Earth collision hazard[END_REF][START_REF] Harris | A New Estimate of the Population of Small NEAs[END_REF][START_REF] Stuart | Bias-corrected population, size distribution, and impact hazard for the near-Earth objects[END_REF][START_REF] Halliday | Detailed data for 259 fireballs from the Canadian camera network and inferences concerning the influx of large meteoroids[END_REF]ReVelle, 2001;[START_REF] Brown | The flux of small near-Earth objects colliding with the Earth[END_REF].

Significant cratering asymmetries, both latitudinal and longitudinal, are found on the lunar surface. These are the result of both the high proportion of low inclination encounters and the synchronous rotation of the satellite. If the resulting bias in the estimate of absolute ages is in general moderate with respect to other uncertainties associated with the crater chronology method, and therefore should not invalidate most of the published ages, this systematic effect could be as large as 500 Ma for geologic units of 2 Ga close to (0 • N, 90 • W )

or (±65 • N, 90 • E), and should be accounted for when determining relative ages between different surfaces.

Our synthetic size-frequency distribution of lunar craters matches the standard production functions (see Neukum et al., 2001a), providing that a porous megaregolith is accounted for at crater diameters smaller than a few kilometers. The absolute number of craters formed during the last 3 billion years is consistent with a constant impact flux and the present day cratering rate.

The shape of the size-frequency distribution is predicted to be time-dependent, since older surfaces would possess a thicker layer of fractured megaregolith.

We have provided an updated crater chronology curve that excludes some debated calibration points based on geologic considerations, and accounts for the spatially-dependent cratering rate at the calibration units. The proposed curve favors a longer period of constant impact flux, extending up to 3.5 Ga ago. We have also given polynomial fits for the crater size-frequency distributions calculated on Mercury, Venus, the Earth and Mars, allowing the easy calculation of age estimates on these planets.

Our model is consistent with the terrestrial cratering records [START_REF] Hughes | A new approach to the calculation of the cratering rate of the earth over the last 125±20myr[END_REF][START_REF] Grieve | The terrestrial cratering record. II -the crater production rate[END_REF] for the last 400 Ma. For the Moon, we find that the age of the Orientale impact basin is 3.73 Ga, which is 30 Ma older than previous estimates. This corresponds to a shorter Late Heavy Bombardment, since it is generally considered that the Caloris basin marks the end of this event. On Mercury, we attribute the same age to the Caloris basin. We finally estimate the venusian surface to be only 240 Ma old on average, which is 2-3 times younger than most previous estimates. Such an age is similar to the oldest oceanic crust on Earth and implies that volcanic resurfacing occurs on Venus at only a slightly reduced rate as on Earth.

Let us call V the encounter conditions with the Earth-Moon system, calculated as explained in section 3.1, and expressed in a right-handed coordinate frame where the origin is at the center of the Earth, the x-axis point towards the Sun, (xy) is the Earth orbital plane, and the z-axis points upward. The symmetry induced by the lunar rotation about the Earth makes the problem only dependent on V and V z . Let us take V y = 0 and V x ≤ 0. Consider also the frame (x yz ), where the x -axis and V are antiparallel. This frame is useful to express the specific angular momentum and eccentricity vector, since it avoids complications due to the initial position of the projectile at infinity.

The rotation matrix from (x yz ) to (xyz) is therefore

M =               √ V 2 -V 2 z V 0 Vz V 0 1 0 -Vz V 0 √ V 2 -V 2 z V               . (A.16)
The encounter velocity and initial position of the projectile on the cross section disk, perpendicular to x , expressed in (x yz ) as

V =               -V 0 0               (x yz ) , r =               ∞ b cos α b sin α               (x yz ) , (A.17)
where b is the impact parameter and α defines the angular position on the cross section disk with respect to the y-axis.

50

From these initial conditions, the specific angular momentum of the projectile and eccentricity vector are easily obtained by equations (A.4) and (A.5), and expressed in the (xyz) frame by the use of the rotation matrix M (equation A.16). The norm of the eccentricity vector is

E = 1 + b a m 2 V v m 4 , (A.18)
where a m and v m are the lunar semi-major axis and orbital velocity. The inclination with respect to the system orbit plane (xy) is (equation A.7)

I = cos -1   1 - V 2 z V 2 cos α   . (A.19)
Note that the eccentricity E and inclination I of the hyperbolic trajectory have been capitalized in order to distinguish them from the orbital elements e and i for elliptic orbits about the Sun. The pericenter distance q is (equation A.2)

q = a m (E -1)/ V v m 2 , (A.20)
the nodal vector is calculated with equation (A.8), and the argument of pericenter is finally expressed by equation (A.11) as

cos ω = Vz V √ E 2 -1 -1 -V 2 z V 2 sin α E sin I . (A.21)
For an intersection to occur between the projectile and Moon, the projectile must cross the lunar orbit plane at a distance a m from the Earth. Since V z is indifferently positive or negative ( Öpik, 1951), the symmetry with respect to the system orbit plane allows to restrict ourselves to one of the nodes, here chosen to be the descending node. In this case, Θ = π implies (equation A.15) cos f = -cos ω , (A.22) and r = a m gives (equation A.1)

cos f = q/a m (1 + E) -1 E . (A.23)
By using the two above equations, with E, q and cos ω expressed as a function of V, b and α from equations (A.18), (A.20) and (A.21), it is found that the value of the impact parameter that allows intersection is .24) and that there is only one value of b that allows intersection on the descending node for a given angular position on the lunar orbit cross section α.

b m = a m 2   V z V sin I 2 + 4 1 + tan α tan I v m V 2 - V z V sin I   . (A
For hyperbolic orbits, the relative encounter velocity with the Moon U is given in [START_REF] Shoemaker | Cratering time scales for the Galilean satellites[END_REF] as a function of q, e and I. For the appropriate value of the impact parameter that allows collision, we rewrite these equations as

                                       U 2 = V 2 + v 2 m -2v m U Y , U 2 X = V 2 1 -bm am 2 + 2v 2 m , U Y = V bm am cos I -v m , U 2 Z = V 2 bm am 2 sin 2 I . (A.25)
where U is expressed in the frame (XY Z) attached to the center of the Moon, where the X-axis point towards the Earth, the Y-axis points to the antapex of motion, (XY) defines the Lunar orbital plane, and the Z-axis points upward.

Since we have restricted ourselves to the descending node, U Z is negative. But symmetry requires that U Z is positive with the same probability. The sign of U X is to be discussed in the following.

between 0 and 2π with an equal probability. While this is appropriate for elliptic orbits that precess about the Sun, the ω of our hyperbolic trajectory is constrained by the geometry of encounter with the Earth-Moon system. The encounter probability with the Moon can be expressed as follows

P m = ∆(b 2 ) τ 2 × ∆t T m , (A.26)
where ∆(b 2 ) is the interval around b 2 m that allows the minimum distance between the projectile's and target's orbits to be smaller than the gravitational cross section radius of the Moon, τ m . We are interested in the square of b here, whose probability distribution on the lunar orbit cross section is uniform. The ratio ∆(b 2 ) τ 2 corresponds to the fraction in space allowing the collision, since b 2 can take values between 0 and τ 2 . Similarly, the term ∆t Tm is the fraction of time that allows the collision, T m being the orbital period of the Moon. This latter term is given by (see Öpik, 1951): .27) At this point we need to express ∆(b 2 ) in terms of the relevant variables of our problem. At the node, Θ = π; denoting L max the distance between r Θ=π and a m such as the minimum distance between the two orbits equals τ m , we have

∆t T m = τ m 4a m U U 2 -U 2 Y . ( A
2L max ∆(b 2 ) ∂r Θ=π ∂b 2 b=bm , (A.28) with L max = τ m 1 + U 2 X /U 2 Z . (
P m = τ 2 m a 2 m U |U Z | v 2 m 2v 2 m + V 2 1 sin I sin I + 1 -Vz V 2 sin α -Vz|U Z | v 2 m sin I 2 2 sin I + 2 1 -Vz V 2 sin α -Vz|U Z | v 2 m sin I . (A.30)
If U X ≤ 0 and q ≤ R e simultaneously, R e being the radius of the Earth, a collision with the planet occurs, and P m is set to 0 (the Moon cannot be impacted, as the sign of U X indicates that the projectile has already cross the Earth). Denoting b e the impact parameter corresponding to q = R e we have from equation A.20: A.31) and q ≤ R e corresponds to b m ≤ b e . The sign of U X is known from

b e = 2R e a m v m V 2 + R 2 e , (
U X = -U 2 X for √ E 2 -1 sin α/ tan I ≥ 1 U X = U 2 X for √ E 2 -1 sin α/ tan I < 1 , (A.32)
where the first and second cases correspond respectively to an argument of pericenter comprised between 0 and π and between π and 2π (equation A.12).

For instance, in the first case, the descending node is reached after the minimum distance to Earth, and the projectile is moving away from the Earth.

Finally, the impact probability with Earth is given by A.33) where the factor 1/2 comes from the restriction to one of the two possible taken into account at the end of this section). For a zero obliquity, we have cos β = (cos λ cos ϕ , cos λ sin ϕ , sin λ) × U ||U|| .

P e = 1 2 b 2 e τ 2 -P m , (
(A.49)

We finally need the impact velocity u and incidence angle with respect to the surface, θ. The conservation of energy implies

U 2 2 = u 2 2 - GM R , (A.50) giving u = U √ 1 + Γ , (A.51)
The specific angular momentum at the time of impact is (equation A.4)

h = R u sin(θ + π/2) . (A.52)
The conservation of h allows identification with equation (A.42), giving cos θ = b/τ , (A.53) and the impact angle is finally given as a function of β as

θ = cos -1 1 + µ 2 √ 1 + Γ 1 -cos 2 β . (A.54)
Let us finally account for the obliquity of the target. The obliquity reduces spatial variations of the impact rate, since the target's sphere is moved under the projectile flux by precession of the spin axis. Let us denote ϕ the longitude of the spin axis in the (XY Z) frame, and W the matrix 2004)). Our compilation at large sizes is augmented by including the observed size-frequency distribution of Mars-crossing objects with sizes greater than 4 km, scaled to the terrestrial impact rates of [START_REF] Stuart | Bias-corrected population, size distribution, and impact hazard for the near-Earth objects[END_REF] (grey diamonds). The red curve is the best fit polynomial.

W =               cos 2 ϕ + sin 2 ϕ cos cos ϕ sin ϕ (1 -cos ) -sin ϕ sin cos ϕ sin ϕ (1 -cos ) sin 2 ϕ + cos 2 ϕ cos cos ϕ sin sin ϕ S -cos ϕ sin cos               -1 , (
C 0,0 1 1 1
C 1,-1 -1.7779020 × 10 -1 -1.4591830 × 10 -1 -1.2670400 × 10 -1 C 2,0 -6.3209891 × 10 -2 -6.2923420 × 10 -2 -6.2755592 × 10 -2 C 2,2 -1.7937283 × 10 -3 -1.2524090 × 10 -3 -9.8998262 × 10 -4

C 4,0 -9.9735381 × 10 -3 -1.0141450 × 10 -2 -1.0222370 × 10 -2 C 2,0 -1.3954110 × 10 -1 -1.3874170 × 10 -1 -1.3831914 × 10 -1 C 2,2 -3.3499412 × 10 -3 -2.2729362 × 10 -3 -1.7822560 × 10 -3 C 3,-1 2.9040180 × 10 -3 2.3263713 × 10 -3 1.9948010 × 10 -3 C 4,0 2.7412530 × 10 -3 2.7863080 × 10 -3 2.8083000 × 10 -3

Table 7 Spherical harmonic coefficients of the relative impact flux ∆φ(λ) for the terrestrial planets.

Planet

  . Putting this altogether, we obtain the relation d(D, u, θ) that gives the impactor diameter d as a function of the crater diameter D, impact velocity u and incidence angle θ (appendix A.4, equations (A.63 -A.59)). The impactor diameter d required to create a crater of size D is ultimately a function of λ, φ and U (see equations A.51 and A.54).

  , our model reproduces both the Hartmann and Neukum production functions within the 100 m -300 km diameter range, for respective values of T equal to 250 and 700 m, respectively. For these diameter ranges, the maximum discrepancy between our model and the Neukum production function is only 30 % at 200 m, and always less than 20 % for craters larger than 500 m. Below 100 m, we note that our model is in reasonable agreement with the Neukum production function, and we leave the implications for the contribution of secondary craters to further investigations. The model production function proposed by Marchi et al. (2009) is also shown in figure 2. A detailed comparison with this latter study is given in the discussion section. [Fig. 2 about here.]

  about here.] [Fig. 4 about here.]

Figure 8

 8 Figure 8 also shows the range of bias that could occur if spatial variations are not taken into account in the crater chronology method. Radiometric ages obtained from future missions, or crater count data from images, that are far from the Apollo and Luna sampling sites might need to have their associated local crater densities corrected by a factor up to ∼25% to obtain the corresponding globally averaged value. Since the age/density relationship is approximatively linear for the last 3.5 Ga, ages could be biased by the same

  where N obs (> D k ) are the K measured cumulative crater densities with k ∈ [1, K]. Each data possesses an uncertainty σ k = N obs (> D k )/A, where A is the area where craters have been counted. The solution (e.g., Tarantola 2001) is given by

  synthetic lunar SFD fitted to crater counts performed over the Orientale basin (data fromNeukum et al. (2001a)). The crater SFD is shown in R-plot form,R(> D) = D 3 dN (>D)dD , where changes in slope are emphasized with respect to a reference cumulative distribution having a power-law slope of -2. We use a best-fitting megaregolith thickness of 500 m, and exclude in the fit craters smaller than 1.5 km in diameter that have reached saturation level. We find an age of 3.73 ± 0.01 Ga. The difference with Neukum et al. (2001a), who quoted an age of 3.70, lies principally in the updated coefficient values of the calibration curve. Without accounting for spatial cratering variations, the age would have been estimated to be 3.75 Ga. The bias in age from spatial cratering variations is here moderate, as we are in the exponential part of the calibration curve. On Earth, Hughes (2000) used a nearest neighbors technique and known crater ages to calculate the size-dependent formation rate of craters over the last 125 ± 20 Ma. The best fit between our model production function and the estimates of

Fig. 1 .

 1 Fig. 1. Impact probability per year on Earth, for impactors larger than a given diameter. Estimates come from atmospheric records or Öpik probabilities derived from telescopic observations (black triangles: Halliday et al. (1996); black diamonds: ReVelle (2001); grey triangles: Brown et al. (2002); grey squares and white triangles: Rabinowitz et al. (2000); grey circles: Harris, A. W. (2002); white circles: Morbidelli et al. (2002); black squares: Stuart and Binzel (2004)). Our compilation at large sizes is augmented by including the observed size-frequency distribution of Mars-crossing objects with sizes greater than 4 km, scaled to the terrestrial impact rates of[START_REF] Stuart | Bias-corrected population, size distribution, and impact hazard for the near-Earth objects[END_REF] (grey diamonds). The red curve is the best fit polynomial.

Fig. 2 .Fig. 4 .

 24 Fig. 2. Model production function of lunar craters, for one year, in comparison with the Hartmann and Neukum measured production functions, and the model production function of Marchi et al. (2009). Respective megaregolith thicknesses of 700 and 250 m allow to fit either the Neukum or Hartmann production functions in the diameter range 2-20 km. The thin dotted red curve is obtained by using only the non-porous scaling relation.

Fig. 5 .Fig. 6 .Fig. 7 .

 567 Fig. 5. Relative cratering rate on the Moon for the current Earth-Moon distance and for crater diameters larger than 1 km. This image is plotted in a Mollweide projection centered on the sub-Earth point.

Fig. 8 .

 8 Fig. 8. Lunar crater chronology calibration curve. Calibration points, listed in table 9, are shown before (grey circles) and after (red circles) accounting for the spatially-dependent cratering rate. Our proposed analytic fit (red line) is compared to the curves of Neukum et al. (2001a) (solid black) and Marchi et al. (2009) (dashed grey). Crisium basin, that was not included in the fit, is shown in blue. Also shown is the range of bias in age that could be produced by spatial variations in the cratering rate (purple area).

Fig. 9 .

 9 Fig. 9. Age estimates. Distributions plotted either in cumulative form or R-plot. Data represented by white circles are not included in the fit, either because of saturation (Moon), erosion (Earth) or atmospheric shielding (Venus). (a) Orientale basin: 3.73 ± 0.01 Ga with a megaregolith thickness of 500 m. Data from Neukum et al. (2001a) (b) Terrestrial cratering record: 138±6 Ma for the estimates of Hughes (2000) (circles), 390 ±75 Ma for the Phanerozoic cratons investigated by Grieve and Dence (1979) (squares) (c) Venusian surface: 240 ± 13 Ma. Data from the Magellan mission (d) Caloris basin lineated sculpture (grey): 3.73±0.2 Ga with a megaregolith thickness of 900 m; Caloris interior plains (black): 3.30±0.3 Ga with a megaregolith thickness of 250 m. Data from Fassett et al. (2009).

Table 4

 4

	about here.]
	[Table 5 about here.]
	[Table 6 about here.]
	[Table 7 about here.]
	[Table 8 about here.]
	5 Crater chronology
	Neukum et al. (2001a) (see also

Table 1

 1 Impactor size distribution: log s(> d) = 10 i=0 s i (log d km ) i . Planetary crater size-frequency distributions for 1 year (km -2 ): log C(> D) = 10 i=0 C i (log D km ) i , D ∈ [0.1 -1000] km.

	s 0	s 1	s 2	s 3	s 4	s 5
	1.0	3.1656E-01 1.0393E-01 5.7091E-02 -8.1475E-02 -2.9864E-02
	s 6	s 7	s 8	s 9	s 10
	1.3977E-02 5.8676E-03 -4.6476E-04 -3.8428E-04 -3.7825E-05

Table 3

 3 Impact rate, mean impact velocity and surface gravity for the inner planets, normalized to the Moon's.

		Moon Mercury Venus Earth	Mars	Mars
						(long term) (today)
	impact rate ratio	1	1.82	1.75	1.58	2.76	3.20
	mean velocity ratio	1	2.16	1.28	1.04	0.53	0.54
	surface gravity ratio	1	2.2	5.3	5.9	2.2	2.2

Table 5

 5 Spherical harmonic coefficients of the lunar relative impact flux ∆φ(λ, ϕ) for Earth-Moon separations of 30, 45, and 60 Earth radii.

	C lm	30 Earth radii	45 Earth radii	60 Earth radii

Table 6

 6 Spherical harmonic coefficients of the lunar relative cratering rate ∆C(D > 1, λ, ϕ) for Earth-Moon separations of 30, 45, and 60 Earth radii.

	C lm	30 Earth radii	45 Earth radii	60 Earth radii
	C 0,0	1	1	1

C 1,-1 -2.2715950 × 10 -1 -1.8571530 × 10 -1 -1.6092760 × 10 -1

Table 8

 8 Spherical harmonic coefficients of the relative cratering rate ∆C(D > 1, λ) for the terrestrial planets.

	Planet	C 00	C 20	C 40
	Mercury	1	4.5850560 × 10 -2	2.4749320 × 10 -3
	Venus	1 -1.8545722 × 10 -3	3.5865970 × 10 -4
	Earth	1 -7.6586370 × 10 -2	2.4353234 × 10 -4
	Mars (today)	1	3.3900970 × 10 -1	2.2655340 × 10 -3
	Mars (long-term average)	1	1.7986312 × 10 -1 -4.3484250 × 10 -4

Table 9

 9 Lunar crater chronology calibration points used in this study, and their spatial correction.

	Geologic unit	Location	Age	Measured Globally averaged
				crater density	crater density
			(Ga)	(10 -4 km -2 )	(10 -4 km -2 )
	1 Cone crater	-3.7 • N, -17.5 • E 0.025	0.21	0.169
	2 North Ray crater	-9 • N, -15.6 • E 0.053	0.44	0.390
	3 Tycho crater	-43.4 • N, -11.1 • E 0.109	0.90	0.824
	4 Copernicus	9.7 • N, -20.1 • E	0.80	7.15	5.77
	5 Ocean Procellarum (A12)	18.4 • N, -57.4 • E	3.15	36	29.7
	6 Mare Crisium (L24)	17.0 • N, 59.1 • E	3.22	30	30.2
	7 Mare Imbrium (A15)	32.8 • N,-15.6 • E	3.30	32	27.7
	8 Mare Fecunditatis (L16)	-7.8 • N,51.3 • E	3.41	33	32.4
	9 Mare Tranquilitatis (young) (A11)	8.5 • N,31.4 • E	3.58	64	60.1
	10 Mare Serenitatis (A17)	28.0 • N,17.5 • E	3.75	100	93.6
	11 Mare Tranquilitatis (old) (A11)	8.5 • N,31.4 • E	3.80	90	84.5
	12 Imbrium/Apennines	18.9 • N,-3.7 • E	3.85	350	301
	13 Fra Mauro (A14)	-3.7 • N,-17.5 • E	3.85	370	298
	14 Descartes (A16)	-9.0 • N,15.6 • E	3.85	340	306

Acknowledgements

We thank B. Ivanov and S. Werner for constructive reviews and comments that improved this manuscript. We thank W. Bottke and A. Morbidelli for

their orbital distribution model. MLF has been supported by a grant from the Centre National d'Etudes Spatiales (CNES) and acknowledges the Laboratory of Planetology and Geodynamics of Nantes, France, for hosting and fruitful discussions. Figures were created using the Generic Mapping Tools of [START_REF] Wessel | Free software helps map and display data[END_REF].

A Derivations

Here we describe in detail how we compute the lunar and planetary cratering rates. We start by recalling how to obtain orbital elements in the 2-body problem, which is required for calculating hyperbolic trajectories of objects as they encounter a planet or Moon. Second, we adapt Öpik impact probabilities, originally intended for planets, to the case of a moon in synchronous rotation.

Third, we derive equations to calculate the impact rate as a function of position on the target surface (planet or moon). Finally, we recall the equations that allow us to calculate crater diameters from impactor diameters and impact conditions.

A.1 Orbital elements

Consider a massless particle (our projectile) orbiting around a massive body M (our target) isolated in space (two-body problem), described in a Cartesian coordinate system (oxyz) whose origin is the center of the target, with x, ŷ and ẑ the associated unit vectors. The plane (oxy) is chosen as the orbital plane of the target body. All the following definitions can be found in standard celestial mechanics textbooks (e.g., Murray and Dermott, 2000, 606 pp.). The distance of the projectile to the center of the target is

where p is the semilatus rectum, e the eccentricity and f the true anomaly.

The pericentre distance (f = π) is given by

The semilatus rectum is

where M is the mass of the target body, and the specific angular momentum, constant and perpendicular to the projectile orbital plane, is

The eccentricity vector is

and the eccentricity is simply e = ||e||. Alternatively, the eccentricity can be calculated from the particle's specific energy ξ as

The inclination of the orbit with respect to the (oxy) plane is calculated from

where i ∈ [0, π]. The nodal vector, that points towards the ascending node, is

and the longitude of the ascending node, measured from (ox), is calculated

The argument of pericentre ω is given by

(A.12)

The true anomaly is calculated from

Finally, defining Θ as the angle that the projectile makes with respect to the line of nodes, we have

and the projectile crosses the target's orbital plane at Θ = 0 and Θ = π, respectively on an ascending and descending trajectory.

nodes. For a given encounter V with the Earth-Moon system, P m , U and P e are calculated for α ∈ [0, 2π]. The probability distribution of the lunar encounter conditions is

where p(V) is the probability distribution of the encounter distribution with the Earth-Moon system and δ is the Kronecker function. Finally, the Earth/Moon impact ratio is The impact rate at β is proportional to the area of the annulus 2π b db divided by the elementary target surface 2πR 2 sin β dβ. The total impact rate is proportional to the gravitational cross sectional area, πτ 2 , divided by the total area of the target surface 4πR 2 . The relative impact rate, normalized to the global average, is therefore

where f i is the projectile's true anomaly at the time of impact, and f ∞ is the initial true anomaly at infinity, we have

(A.38)

An impact occurs when r = R. The corresponding cosine of the true anomaly is (equation A.1) .39) Since the distance between the projectile and target is initially infinite, the cosine of the intial true anomaly is (equation A.1)

The true anomalie at infinity is negative (since r ∞ • U ∞ < 0). This is also true for f i , as the collision occurs necessarily before the projectile is at pericenter

The specific angular momentum is given by equation (A.

The angular distance β is then simply expressed as a function of latitude λ and longitude ϕ by calculating the scalar product between a given position on the target and the radiant (equation A.44). Note that the dependency in longitude ϕ corresponds to a target body in synchronous rotation. When it is not the case, the flux is longitudinally uniform (the obliquity of the target is then the latitude and longitude (λ , ϕ ) in the (XY Z) frame express as a function of the geographic coordinates (λ, ϕ) as sin λ = W 31 cos λ cos ϕ + W 32 cos λ sin ϕ + W 33 sin λ , (A.56)

and

The impact rate is finally given as a function of the geographic position as

where the precession of the spin axis is simulated by ϕ taking uniform values between 0 and 2π.

A.4 Crater scaling

For our purpose, we express the projectile diameter d required to create a crater with a given final diameter D, under given impact conditions. Let us call D s the final diameter of a simple crater. Large craters collapse due to gravity, such that D > D s . According to [START_REF] Holsapple | The scaling of impact processes in planetary sciences[END_REF],

where D * is a transition diameter believed to be approximately inversely proportional to the surface gravity, and is about D * m = 8.5 km for the Moon,

59

that is

where g m is the lunar surface gravity.

Before wall slumping and rim formation, the diameter of the transient cavity D T is smaller than the final diameter of the simple crater D s by a factor 1.2 × 1.3 = 1.56 (Melosh, 1989, 253 pp.):

Denoting D ⊥ the transient cavity for vertical impacts (θ = π/2) we assume

that is, only the vertical component of the impact velocity has to be accounted for. Other assumed functional dependencies on the incidence angle can easily be used in place of equation (A.62).

In the gravity regime, where the tensile strength of rock is negligible, the projectile diameter is finally given by [START_REF] Holsapple | A crater and its ejecta: An interpretation of Deep Impact[END_REF] as

with ρ and ρ 0 respectively the density of the projectile and target, and g the surface gravity. The scaling parameter are taken from [START_REF] Holsapple | A crater and its ejecta: An interpretation of Deep Impact[END_REF] as K = 1.17, ν 1 = 0.22 and ν 2 = 0.31 for the non-porous gravity scaling regime (impacts in water), and K = 1.03, ν 1 = 0.17 and ν 2 = 0.332 for the porous scaling. For a given crater diameter D, the impactor size d is known from the successive use of the above equations, from equation (A.63) to equation (A.59).

List of Figures