11institutetext: University of Catania, Catania, CT 95123, ITALY,
11email<EMAIL_ADDRESS>
11email<EMAIL_ADDRESS>
# Unleashing the Power of User Reviews: Exploring Airline Choices at Catania
Airport, Italy

Vincenzo Miracula 11 Antonio Picone 11

###### Abstract

This study aims to investigate the possible relationship between the
mechanisms of social influence and the choice of airline, through the use of
new tools, with the aim of understanding whether they can contribute to a
better understanding of the factors influencing the decisions of consumers in
the aviation sector. We have chosen to extract user reviews from well-known
platforms: Trustpilot, Google, and Twitter. By combining web scraping
techniques, we have been able to collect a comprehensive dataset comprising a
wide range of user opinions, feedback, and ratings. We then refined the BERT
model to focus on insightful sentiment in the context of airline reviews.
Through our analysis, we observed an intriguing trend of average negative
sentiment scores across various airlines, giving us deeper insight into the
dynamics between airlines and helping us identify key partnerships, popular
routes, and airlines that play a central role in the aeronautical ecosystem of
Catania airport during the specified period. Our investigation led us to find
that, despite an airline having received prestigious awards as a low-cost
leader in Europe for two consecutive years 2021 and 2022, the ”Catanese” user
tends to suffer the dominant position of other companies. Understanding the
impact of positive reviews and leveraging sentiment analysis can help airlines
improve their reputation, attract more customers, and ultimately gain a
competitive edge in the marketplace.

###### keywords:

computational social science, network flight, sentiment analysis, artificial
intelligence

##  1 Introduction

In many social and biological systems, individuals rely on the observation of
others to adapt their behaviors, revise their judgments, or make decisions
[1],[2],[3]. In human populations, access to social information has been
greatly facilitated by the continued growth of communication technology.
Indeed, people are constantly exposed to a steady stream of others’ opinions,
advice, and judgments about political ideas, new technologies, or commercial
products [4]. When facing peer opinions on a given issue, people tend to
filter and integrate the social information they receive and adjust their
beliefs accordingly [5]. At the scale of a group, repeated local influences
among group members can give rise to complex patterns of opinion dynamics such
as consensus formation, polarization, or fragmentation [6],[7],[8]. For
example, it has been found that people who share similar extreme views, such
as racial prejudice, tend to strengthen their judgment and trust after
interacting with each other [9].

Similar mechanisms of opinion dynamics can take place in a variety of social
contexts, such as within a group of friends exchanging opinions about their
willingness to be vaccinated against the flu [10],[11]. On an even larger
scale, local influences among friends, family members, or colleagues, often
combined with global mass media effects, constitute an effective mechanism
that drives opinion formation during elections, shaping cultural markets [12],
producing amplification or attenuation of risk perception [13], and shaping
public opinion on social issues, such as atomic energy or climate change [14].

Given the remarkably broad scope of social phenomena that are shaped by social
influence and opinion dynamics, it is surprising that the behavioral
mechanisms underlying these processes remain poorly understood. Major
questions remain unanswered: how do people adjust their judgments during
social interactions? What are the heuristics underlying opinion adaptation?
And how do these local influences ultimately generate global patterns of
opinion change? Much of the existing modeling work on opinion dynamics has
been approached from a physics-based point of view, where the basic mechanisms
of social influence derive from analogies with physical systems, particularly
the spin systems [15], [16].

The wide variety of existing models assumes that individuals hold different
rules for adapting their opinions, such as imitation [17], averaging with
people with similar opinions [18], or following the majority [19]. Although
informative about the complex dynamics that can potentially emerge in a
collective context, these simulation-based contributions share a common
drawback: the absence of empirical verification of the models’ hypotheses.

Indeed, it’s quite a challenge to trace and measure how opinions change under
experimental conditions, as these changes depend on numerous social and
psychological factors, such as individuals’ personalities, confidence levels,
credibility, social status, or persuasive power. In other fields like
cognitive sciences and social sciences, laboratory experiments have been
conducted to examine how people integrate feedback from others to revise their
initial responses to factual questions [5], [19], [20]. However, the outcomes
of local rules for opinion adaptation have yet to be employed in studying the
collective dynamics of the system, and it remains unclear how social influence
unfolds over time in large-scale social contexts [21]

This study aims to investigate the relationship between social influence
mechanisms and airline choice while seeking to address some fundamental
questions. What are the factors that sway individuals towards one airline over
another? Does the user have real freedom to choose the airline carrier? Is
there a hierarchy among airports? We specifically focus on flights departing
from and arriving at Catania Airport, given the significance of this airport
hub in the heart of the Mediterranean.

##  2 Materials and Method

To gather data for our analysis, we have chosen to extract user reviews from
two well-known platforms: Trustpilot and Google. Trustpilot offers a platform
for users to share their feedback and reviews about various businesses and
services, while Google allows users to leave reviews for apps, products, and
services.

By collecting reviews from Trustpilot and Google, our objective is to capture
a wide spectrum of user opinions and experiences. These platforms attract a
substantial user base that shares feedback on various aspects of the apps,
products, or services they have utilized. Analyzing reviews from both
platforms enables us to gain valuable insights into the sentiments and
perspectives of users within our specific context.

In order to gather the necessary data for our analysis, we utilized Python in
conjunction with the BeautifulSoup library. The process involved web scraping
to extract information from Trustpilot, while an official API was employed to
extract data from Google. For the scraping part, we leveraged BeautifulSoup’s
powerful parsing capabilities to navigate the HTML structure of Trustpilot’s
web pages. With BeautifulSoup, we were able to locate and extract relevant
information such as user reviews, ratings, and reviewers’ names. By
programmatically traversing the HTML elements, we systematically collected the
necessary data for further analysis. On the other hand, when it came to
extracting Google reviews, we relied on an official API. They provide a
structured and documented way to interact with a particular platform or
service. In this case, the Google API allowed us to access and retrieve
specific data related to reviews, ratings, and other relevant details directly
from the Google platform. By following the API’s guidelines and utilizing the
appropriate endpoints, we could efficiently retrieve the desired data in a
standardized format.

It’s important to note that while scraping involves extracting data from the
presentation layer of a web page, APIs offer a more streamlined and controlled
approach to accessing data from the source directly. APIs often provide
specific methods and endpoints tailored to retrieve data in a structured
manner, ensuring consistent and reliable results. This eliminates the need for
parsing and navigating HTML structures, simplifying the data extraction
process and ensuring compatibility with the platform’s terms of service.
Scraping certainly requires higher skills than the use of an API and a further
cleaning of the data, sometimes soiled by the presence of HTML tags that must
be removed, in order to obtain text in its classic form.

There is a common phenomenon known as “negativity bias”[23] that suggests
people are more inclined to share negative experiences and write negative
reviews on the internet. Several factors contribute to this trend.

Firstly, negative experiences often have a more significant emotional impact
on individuals compared to positive experiences. As a result, people may feel
a stronger urge to express their dissatisfaction or frustration by writing
negative reviews as a means of catharsis. Secondly, negative reviews can serve
as a warning to others, potentially preventing them from encountering similar
negative experiences. People may feel a sense of responsibility to share their
negative encounters to protect or inform others.

Thirdly, there is a perception that negative reviews have more weight and
credibility. Some individuals may believe that expressing criticism or
highlighting flaws will lead to improvements in the product, service, or
overall customer experience. This belief can incentivize people to voice their
negative opinions.

Twitter - on the other side - can indeed serve as a platform to mitigate the
negativity bias issue when it comes to online reviews and feedback. The real-
time nature of Twitter provides an advantage in mitigating the negativity bias
issue. Users can share their thoughts and experiences on the platform as they
happen, allowing for a more balanced representation of both positive and
negative sentiments. This immediacy enables individuals to express their
satisfaction, appreciation, or dissatisfaction with airlines or their services
at the moment, capturing a diverse range of opinions and experiences. Data
from this source were also extracted by using the official Twitter API.

By combining web scraping and API usage, we were able to collect a
comprehensive dataset from Trustpilot, Google, and Twitter, encompassing a
wide range of user opinions, feedback, and ratings. This rich dataset allows
for in-depth analysis and insights into the experiences and perspectives of
users within our specific research context. The dataset, consisting of the
cleaned textual reviews and the metadata mentioned, was saved in serialized
JSON (JavaScript Object Notation) format, for easy reading by both a human
being and a computer, as further analysis of sentiment followed, obtained
through the use of machine learning techniques applied to Natural Language
Processing (NLP).

To ensure accurate results, we initially performed text pre-processing, which
involved cleaning and preparing the airline reviews for analysis. This step
included removing unnecessary characters, punctuation, and stopwords, as well
as tokenizing the text into meaningful units. Next, we employed BERT
(Bidirectional Encoder Representations from Transformers)[22] for sentiment
analysis. BERT is a state-of-the-art NLP model that has demonstrated
exceptional performance in understanding language context. Instead of training
BERT from scratch, we utilized a pre-trained BERT model that had already
learned from a vast amount of text data. To tailor the pre-trained BERT model
specifically for sentiment analysis of airline reviews, we further fine-tuned
it. Fine-tuning involved training the model on a labeled dataset that included
positive and negative sentiment labels. By leveraging the pre-existing
knowledge and linguistic understanding of the BERT model, we refined it to
focus on discerning sentiment in the context of airline reviews.

This approach allowed us to benefit from BERT’s powerful language
representation capabilities and effectively adapt it to the task of sentiment
analysis for airline reviews. The fine-tuned BERT model was then able to
provide more accurate and reliable sentiment predictions, facilitating a
comprehensive analysis of the sentiments expressed in the reviews.

The data gathered thus far, while valuable, do not encompass the entirety of
the information required to conduct our analysis comprehensively. In addition
to the existing data, we need to include data from airline schedules,
including scheduled departures and arrivals; real departures and arrivals
times along with data about the airline of the aircraft. By integrating these
additional data points, we can form a more holistic view of the airline
industry and the dynamics of flights to and from Catania Airport.

To retrieve and save flight data for flights to and from Catania airport, we
developed Python code leveraging the ADSBexchange.com website.
ADSBexchange.com provides online access to ADS-B (Automatic Dependent
Surveillance-Broadcast) data, enabling precise tracking of aircraft movements,
including takeoff and landing times. We utilized Python to interact with the
APIs to get the required data from the website, extracting information about
specific flights associated with Catania airport. The retrieved data,
including flight details and timestamps, were then stored in a MongoDB
database. MongoDB, a popular NoSQL database, offered a flexible and scalable
solution to store and organize the flight data efficiently. By combining web
scraping techniques, the ADSBexchange.com APIs, and MongoDB, we were able to
gather and store accurate flight information for comprehensive analysis and
further insights.

##  3 Results

For our analysis, we have chosen a specific period from February 18, 2023, to
May 18, 2023, spanning three months. This timeframe allows us to capture a
substantial amount of data for in-depth analysis. During this period, a total
of 16,590 flights were recorded, encompassing both arrivals and departures at
the airport. To gain insights into the prominent airlines operating at the
airport, we extracted data on the number of flights operated by each airline
to and from the airport. By evaluating the flight volume, we were able to
identify the top 10 airlines (see Table 1) with the highest number of flights,
providing us with valuable information on the major players in the aviation
industry at Catania airport during the specified timeframe.

Table 1: Flights per airline # | Airline | Flights
---|---|---
1 | Ryanair | 6257
2 | ITA Airways | 2271
3 | Wizz Air | 2021
4 | easyJet | 1796
5 | Eurowings | 332
6 | Air Malta | 274
7 | Vueling | 261
8 | Poste Air Cargo | 257
9 | Volotea | 237
10 | Lufthansa | 234

It should be noted that among the top 10 airlines, ”Poste Air Cargo” appears,
which will not be examined because it is a flight carrying cargo and not a
commercial passenger one.

In addition to this analysis, we conducted a network analysis (Figure 1) to
obtain a more comprehensive understanding of the airline ecosystem at Catania
airport during the specified period. Using flight route data and the involved
airlines, we constructed a network model that represents the connections
between the airlines based on the operated flights.

Figure 1: Flight route network Catania airport

Through this network analysis, we were able to identify the relationships and
interactions among the airlines, including both direct and indirect
connections. We evaluated the centrality of airlines in the network,
considering metrics such as the degree of connectivity, the importance of
connections, and the strategic position within the network.

This network analysis provided us with a deeper insight into the dynamics
among the airlines and helped us identify key partnerships, the most frequent
routes, and the airlines playing a central role in the aviation ecosystem at
Catania airport during the specified period. This information can be valuable
in understanding the competitive structure of the aviation industry and the
collaboration strategies among the involved airlines.

###  3.1 Sentiment analysis for top 10 company

The collection of text data starts from the date of 26 April 2022, chosen as
the starting point for the collection of reviews because it corresponds to the
date of the first available review for ITA Airways - which is the newest
airline in the list. By starting data collection from this date, we ensure
that our data set remains balanced and includes revisions for the same period.
This approach allows us to maintain fairness and completeness in our analysis
by including an equal representation of customer sentiment across different
airlines over time (see Table 2).

Table 2: Sentiment Analysis on airline’s review

Airline | Positive | Negative
---|---|---
Ryanair | 10.56% | 89.44%
ITA Airways | 8.94% | 91.06%
Wizz Air | 3.6% | 96.4%
easyJet | 10.94% | 89.06%
Eurowings | 5.06% | 94.94%
Air Malta | 6.56% | 93.44%
Vueling | 3.88% | 96.12%
Volotea | 56.85% | 43.15%
Lufthansa | 3.29% | 96.71%

Throughout our analysis, we have observed an intriguing trend of average
negative sentiment scores across various airlines. Despite variations in
individual airline performance, the sentiment analysis consistently reveals a
predominant negative sentiment among customer reviews. This finding prompts us
to delve deeper into understanding the underlying factors contributing to this
overall negative sentiment. By conducting further analysis and examining
specific aspects of customer experiences, such as flight delays, and
conducting further analysis on the most common keywords in the reviews - we
aim in the next future - to gain insights into the common pain points and
areas for improvement within the airline industry as a whole.

At the same time, only one airline represents an exception to this negative
trend. We have been intrigued by the consistently high positive sentiment
scores for Volotea and sought to understand the underlying reasons. Our
investigation led us to discover that Volotea has been the recipient of
prestigious accolades such as the Europe’s Leading Low-Cost Airline for two
consecutive years: 2021 and 2022.

These awards recognize excellence in the travel industry and are bestowed upon
airlines based on criteria such as customer satisfaction, service quality, and
overall performance. It appears that Volotea’s commitment to delivering
exceptional travel experiences has resonated positively with passengers,
reflected in their favorable sentiments and high sentiment analysis scores. By
consistently meeting or exceeding customer expectations, Volotea has
established a reputation for providing outstanding service, leading to its
remarkable sentiment scores and recognition in the form of esteemed industry
awards.

The importance of using sentiment analysis for outcome prediction lies in its
ability to provide data-driven insights that guide decision-making and
strategic planning. By leveraging sentiment analysis, airlines and industry
stakeholders can gauge customer sentiments and identify areas for improvement.

This information can then be utilized to enhance service quality, address
customer concerns, and align business strategies with customer expectations.
Predicting the outcome of events such as the Europe’s Leading Low-Cost Airline
through sentiment analysis allows airlines to benchmark themselves against
competitors, set performance goals, and make informed decisions to maintain
their competitive edge.

##  4 Conclusion

The decision-making process when choosing airlines is multifaceted, and there
are several reasons why people may opt for airlines with low sentiment
analysis scores.

One factor is the limited availability of alternatives in certain regions or
for specific routes, which may leave passengers with no choice but to use an
airline with lower sentiment scores.

Price considerations also play a significant role, as passengers may
prioritize affordability over sentiment analysis. If a low-rated airline
offers significantly lower fares compared to its competitors, budgetary
constraints may drive passengers to choose it. Loyalty programs can also
influence decisions, as passengers who are loyal to a particular airline’s
frequent flyer program may continue to use it despite lower sentiment analysis
scores. The benefits and rewards offered through loyalty programs can outweigh
negative sentiments expressed by others.

Additionally, individual experiences can vary greatly, even within the same
airline. Passengers may have had positive experiences with certain aspects of
an airline, such as customer service or in-flight amenities, despite lower
sentiment scores overall.

Lastly, convenience and accessibility factors, such as flight schedules,
direct routes, and airport proximity, can heavily impact airline choices.
Passengers may prioritize convenience and choose an airline that offers the
most suitable options for their travel plans, even if sentiment analysis
suggests lower satisfaction levels. Overall, sentiment analysis provides a
general overview, but individual preferences, constraints, and past
experiences play significant roles in airline selection.

While it is true that users may have limited options when it comes to choosing
an airline based on available routes, an interesting question arises when two
airlines share the same route.

In such a scenario, where both airlines offer flights to the desired
destination, users are faced with a decision-making dilemma. In this case, we
have observed a noticeable trend where user preference tends to lean towards
the airline with higher percentages of positive reviews. While factors like
price and convenience play significant roles in decision-making, the sentiment
expressed by previous customers holds substantial influence. Positive reviews
serve as a testament to the airline’s service quality, customer satisfaction,
and overall experience. Users often rely on the experiences and opinions of
others when making their choices, seeking reassurance and validation from
positive reviews. These reviews create a sense of trust and confidence in the
airline’s ability to deliver a satisfactory travel experience. As a result,
airlines that consistently receive positive reviews and higher sentiment
scores have a higher likelihood of being selected by users when faced with
competing options on the same route. Understanding the impact of positive
reviews and leveraging sentiment analysis can help airlines enhance their
reputation, attract more customers, and ultimately gain a competitive edge in
the market.

## References

  * [1] Couzin, I. D., Ioannou, C. C., Demirel, G., Gross, T., Torney, C. J., Hartnett, A., … & Leonard, N. E. (2011). Uninformed individuals promote democratic consensus in animal groups. science, 334(6062), 1578-1580.
  * [2] Bikhchandani, S., Hirshleifer, D., & Welch, I. (1992). A theory of fads, fashion, custom, and cultural change as informational cascades. Journal of political Economy, 100(5), 992-1026.
  * [3] Moussaid, M., Garnier, S., Theraulaz, G., & Helbing, D. (2009). Collective information processing and pattern formation in swarms, flocks, and crowds. Topics in Cognitive Science, 1(3), 469-497.
  * [4] Wu, F., & Huberman, B. A. (2007). Novelty and collective attention. Proceedings of the National Academy of Sciences, 104(45), 17599-17601.
  * [5] Yaniv, I. (2004). Receiving other people’s advice: Influence and benefit. Organizational behavior and human decision processes, 93(1), 1-13.
  * [6] Mäs, M., Flache, A., & Helbing, D. (2010). Individualization as driving force of clustering phenomena in humans. PLoS computational biology, 6(10), e1000959.
  * [7] Isenberg, D. J. (1986). Group polarization: A critical review and meta-analysis. Journal of personality and social psychology, 50(6), 1141.
  * [8] Galam, S., & Moscovici, S. (1991). Towards a theory of collective phenomena: Consensus and attitude changes in groups. European Journal of Social Psychology, 21(1), 49-74.
  * [9] Myers, D. G., & Bishop, G. D. (1970). Discussion effects on racial attitudes. Science, 169(3947), 778-779.
  * [10] Funk, S., Salathé, M., & Jansen, V. A. (2010). Modelling the influence of human behaviour on the spread of infectious diseases: a review. Journal of the Royal Society Interface, 7(50), 1247-1256.
  * [11] Funk, S., Gilad, E., Watkins, C., & Jansen, V. A. (2009). The spread of awareness and its impact on epidemic outbreaks. Proceedings of the National Academy of Sciences, 106(16), 6872-6877.
  * [12] Salganik, M. J., Dodds, P. S., & Watts, D. J. (2006). Experimental study of inequality and unpredictability in an artificial cultural market. science, 311(5762), 854-856.
  * [13] Kasperson, R. E., Renn, O., Slovic, P., Brown, H. S., Emel, J., Goble, R., … & Ratick, S. (1988). The social amplification of risk: A conceptual framework. Risk analysis, 8(2), 177-187..
  * [14] Latané, B. (1981). The psychology of social impact. American psychologist, 36(4), 343.
  * [15] Castellano, C., Fortunato, S., & Loreto, V. (2009). Statistical physics of social dynamics. Reviews of modern physics, 81(2), 591.
  * [16] Lorenz, J. (2007). Continuous opinion dynamics under bounded confidence: A survey. International Journal of Modern Physics C, 18(12), 1819-1838.
  * [17] Deffuant, G., Neau, D., Amblard, F., & Weisbuch, G. (2001). Mixing beliefs among interacting agents. Advances in Complex Systems, (3), 11.
  * [18] Rainer, H., & Krause, U. (2002). Opinion dynamics and bounded confidence: models, analysis and simulation.
  * [19] Lorenz, J., Rauhut, H., Schweitzer, F., & Helbing, D. (2011). How social influence can undermine the wisdom of crowd effect. Proceedings of the national academy of sciences, 108(22), 9020-9025.
  * [20] Soll, J. B., & Larrick, R. P. (2009). Strategies for revising judgment: How (and how well) people use others’ opinions. Journal of experimental psychology: Learning, memory, and cognition, 35(3), 780.
  * [21] Mason, W. A., Conrey, F. R., & Smith, E. R. (2007). Situating social influence processes: Dynamic, multidirectional flows of influence within social networks. Personality and social psychology review, 11(3), 279-300.
  * [22] Devlin, J., Chang, M., Lee, K. & Toutanova, K. BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding. (2018), https://arxiv.org/abs/1810.04805
  * [23] Vaish, A., Grossmann, T., & Woodward, A. (2008). Not all emotions are created equal: The negativity bias in social-emotional development. Psychological Bulletin, 134(3), 383–403.


Vol.0 (20xx) No.0, 000–000

11institutetext:  Key Laboratory of Orogenic Belts and Crustal Evolution,
School of Earth and Space Sciences, Peking University, 100871 Beijing, China;
<EMAIL_ADDRESS>
22institutetext: Earth Dynamics Research Group, School of Earth and Planetary
Sciences, Curtin University, 6102 WA, Australia;
33institutetext: State Key Laboratory of Remote Sensing Science, Aerospace
Information Research Institute, Chinese Academy of Sciences, 100101 Beijing,
China<EMAIL_ADDRESS>
Received 20xx month day; accepted 20xx month day

# Lunar Cratering Asymmetries with High Orbital Obliquity and Inclination of
the Moon

Huacheng Li 11 Nan Zhang 1122 Zongyu Yue 33 Yizhuo Zhang 11

###### Abstract

Accurate estimation of cratering asymmetry on the Moon is crucial for
understanding Moon evolution history. Early studies of cratering asymmetry
have omitted the contributions of high lunar obliquity and inclination. Here,
we include lunar obliquity and inclination as new controlling variables to
derive the cratering rate spatial variation as a function of longitude and
latitude. With examining the influence of lunar obliquity and inclination on
the asteroids population encountered by the Moon, we then have derived general
formulas of the cratering rate spatial variation based on the crater scaling
law. Our formulas with addition of lunar obliquity and inclination can
reproduce the lunar cratering rate asymmetry at the current Earth-Moon
distance and predict the apex/ant-apex ratio and the pole/equator ratio of
this lunar cratering rate to be 1.36 and 0.87, respectively. The apex/ant-apex
ratio is decreasing as the obliquity and inclination increasing. Combining
with the evolution of lunar obliquity and inclination, our model shows that
the apex/ant-apex ratio does not monotonically decrease with Earth-Moon
distance and hence the influences of obliquity and inclination are not
negligible on evolution of apex/ant-apex ratio. This model is generalizable to
other planets and moons, especially for different spin-orbit resonances.

###### keywords:

Moon, meteorites, meteors, meteoroids, planets and satellites: surfaces

##  1 Introduction

Cratering asymmetry on the lunar surface has been recognized in many studies
(Le Feuvre & Wieczorek 2011; Wang & Zhou 2016). Understanding of such
asymmetry alters the basis of lunar cratering chronology (Hiesinger et al.
2000; Fassett et al. 2012), because it has assumed cratering rate is spatially
uniform on the whole Moon (McGill 1977), which eventually influences the
fundamental understanding of lunar evolution. Quantifying the asymmetry can
rectify the deviation in counting the lunar craters sampled by Apollo and Luna
missions (Hartmann 1970; Neukum et al. 1975; Neukum 1984). Cratering asymmetry
has been also generalized to the surface datings of other planets or moons
(Horedt & Neukum 1984; Neukum et al. 2001a, b; Hartmann & Neukum 2001; Zahnle
et al. 2001; Korycansky & Zahnle 2005). Various factors affecting the
cratering asymmetry on the Moon have been intensively investigated (Hartmann
1970; Neukum et al. 1975; Neukum 1984; Le Feuvre & Wieczorek 2011; Wang & Zhou
2016), and the key factors affecting the cratering asymmetry include (1) the
speed and inclination of asteroids encountering the Moon (Le Feuvre &
Wieczorek 2011) and (2) the distance between the Earth and the Moon (Zahnle et
al. 2001; Le Feuvre & Wieczorek 2011; Wang & Zhou 2016).

Three types of cratering asymmetries, i.e., the leading/trailing asymmetry,
pole/equator asymmetry, and near/far-side asymmetry have been recognized
(e.g., Le Feuvre & Wieczorek 2011; Wang & Zhou 2016). The leading/trailing
asymmetry has been explained by both theoretic derivations (Horedt & Neukum
1984; Le Feuvre & Wieczorek 2011; Wang & Zhou 2016) and numerical simulations
(Gallant et al. 2009; Wang & Zhou 2016). It has been confirmed that the
leading surface receives more impactor fluxes and higher impact speed than the
trailing surface due to the synchronous rotating, while this difference
declines with the Earth-Moon distance increased (Le Feuvre & Wieczorek 2011;
Wang & Zhou 2016). The pole/equator asymmetry has also been numerically
modelled (Gallant et al. 2009; Wang & Zhou 2016), which suggested that low
latitude of the Moon receives more impactor fluxes for the gathering of low
inclination asteroids (Le Feuvre & Wieczorek 2008, 2011). In addition, the
pole/equator asymmetry is found to vary by less than %1 when the Earth-Moon
distance is between 20 and 60 Earth radii (Le Feuvre & Wieczorek 2011). The
mechanism of near/far-side asymmetry has not reached a consensus (Wiesel 1971;
Bandermann & Singer 1973). In previous studies, two factors affecting impact
asymmetry, i.e., orbital obliquity and inclination of the Moon (relative to
the ecliptic), have been usually neglected (Le Feuvre & Wieczorek 2011; Wang &
Zhou 2016). However, these two factors might be important within the first 35
Earth radii of Earth-Moon distance when the Moon quickly left Earth (Ćuk et
al. 2016; Ward 1975). Therefore, it is necessary to investigate the influences
of these two factors on the lunar cratering asymmetries.

In this study, we derive the impact asymmetry reliance on the orbital
obliquity and inclination of the Moon by improving previous empirical models
of leading/trailing and pole/equator asymmetries (Le Feuvre & Wieczorek 2011)
and extending two-dimensional analytic formulas (Wang & Zhou 2016) to the
complete formulas based on three-dimensional geometry. Le Feuvre & Wieczorek
(2011) assumed the orbital obliquity of the Moon was constant when the Earth-
Moon distance is larger than 20 Earth radii. Wang & Zhou (2016) calculated the
cratering asymmetries in a planar model which excludes the influences of the
orbital obliquity and inclination of the Moon. Our analytical formulation
including obliquity and inclination can reveal more features of lunar
leading/trailing asymmetry (Le Feuvre & Wieczorek 2011) and add explicit term
for the pole/equator asymmetry (Wang & Zhou 2016). In Section 2, we derived
the formulas for the distribution of impact flux, normal speed, and cratering
rate on the Moon using the concentration of asteroids encountering with the
Moon and scaling laws that convert asteroids velocities and diameters to the
diameters of craters (Holsapple & Housen 2007). Section 3 shows the resultant
distributions of impact flux, normal speed, and cratering rate based on
formulas in Section 2. This result section also estimates the evolution of the
apex/ant-apex ratio of cratering rate according to the evolution of orbital
obliquity and inclination with different Earth-Moon distances (Ćuk et al.
2016). In Section 4, we verify formulas in Section 2 by comparing with
previous results and explain how the orbital obliquity and inclination
influence the lunar cratering rate asymmetry. Additionally, the influences of
orbital obliquity and inclination of the Moon on the concentration of
asteroids encountering with the Moon are detailed in the appendix.

##  2 Method

This section shows how we calculate the distribution of asteroids impact flux,
impact speed and cratering rate using variables in Table 1. Section 2.1
introduces assumptions and coordinate systems with which we derive the
expression of asteroid’s velocity $\vec{v_{p}}$ and the normal vector
$\vec{n}$ at the impact site. $\vec{v_{p}}$ and $\vec{n}$ will be used in the
following calculations. Section 2.2 uses equations from Wang & Zhou (2016) and
Le Feuvre & Wieczorek (2011) to estimate the impact flux at different impact
sites. These equations are rewritten as functions of $\vec{v_{p}}$ and
$\vec{n}$. Section 2.3 calculates the cratering rate variation using scaling
law from Holsapple & Housen (2007). Obtaining the cratering rate variation
requires the impact flux variation and impact normal speed variation. The
former has been calculated in section 2.2 and the later can be calculated with
minor changes in calculation of impact flux.

Table 1:  Variables or Parameters Used in the Method. Description | Notation | Range
---|---|---
inclination of asteroids’ encounter velocity | $\phi_{p}$ | $[-\frac{\pi}{2},\frac{\pi}{2}]$
azimuth of asteroids’ encounter velocity | $\lambda_{p}$ | $[0,2\pi]$
encounter speed of asteroids | $v_{p}$ | $[19km/s,\sim 20km/s]$
lunar orbit inclination | $i_{1}$ | $[0,\frac{\pi}{2}]$
lunar obliquity relative to the ecliptic | $i_{2}$ | $[0,\frac{\pi}{2}]$
azimuth of lunar orbit normal | $\omega_{1}$ | $[0,2\pi]$
azimuth of lunar spin axis | $\omega_{2}$ | $[0,2\pi]$
lunar true anomaly | $f_{m}$ | $[0,2\pi]$
lunar eccentric anomaly | $E$ | $[0,2\pi]$
lunar mean anomaly | $M$ | $[0,2\pi]$
lunar argument of perihelion | $\omega_{3}$ | $[0,2\pi]$
longitude of impact sites | $\lambda$ | $[0,2\pi]$
latitude of impact sites | $\phi$ | $[-\frac{\pi}{2},\frac{\pi}{2}]$
semi-major axis of lunar orbit | $a_{m}$ | $[25R_{e},60R_{e}]$
eccentricity of lunar orbit | $e$ | [0,1)

###  2.1 Asteroids Velocity and Normal Vector at Impact Site

This model assumes the orbit of the Moon is an ellipse with the Earth as a
focus. Then in the geocentric ecliptic coordinate system (Z-axis is parallel
to the ecliptic normal and X-axis is towards mean equinox of the J2000 epoch),
the position and velocity of the Moon are $\vec{r_{m}}$ and $\vec{v_{m}}$. We
note that the influence of variation of $i_{1},\omega_{1}$ or $\omega_{3}$ on
the lunar velocity can be estimated using Eqs. (2-5) of Ćuk & Burns (2004) and
is $<1\%$ compared to the influence of variation of $f_{m}$. We hence ignored
the variation of $i_{1},\omega_{1}$ or $\omega_{3}$ in deriving the lunar
velocity. In Eq. (2), $G$ and $M_{e}$ are the gravitational constant and the
mass of the Earth respectively.

$\displaystyle\vec{r_{m}}$
$\displaystyle=R_{z}(\frac{\pi}{2}+\omega_{1})R_{x}(i_{1})R_{z}(\omega_{3})\vec{r}$
(1) $\displaystyle\vec{r}$
$\displaystyle=\left[\cfrac{a_{m}(1-e^{2})}{1+e\cos{f_{m}}}\cos{f_{m}},\cfrac{a_{m}(1-e^{2})}{1+e\cos{f_{m}}}\sin{f_{m}},0\right]^{T}$
$\displaystyle\vec{v_{m}}$
$\displaystyle=R_{z}(\frac{\pi}{2}+\omega_{1})R_{x}(i_{1})R_{z}(\omega_{3})\vec{v}$
(2) $\displaystyle\vec{v}$
$\displaystyle=\left[-\sqrt{\frac{GM_{e}}{a_{m}(1-e^{2})}}\sin{f_{m}},\sqrt{\frac{GM_{e}}{a_{m}(1-e^{2})}}(e+\cos{f_{m}}),0\right]^{T}$
$\displaystyle R_{x}(\theta)=\left[\begin{array}[]{ccc}1&0&0\\\
0&\cos{\theta}&-\sin{\theta}\\\ 0&\sin{\theta}&\cos{\theta}\\\
\end{array}\right],\
R_{z}(\theta)=\left[\begin{array}[]{ccc}\cos{\theta}&-\sin{\theta}&0\\\
\sin{\theta}&\cos{\theta}&0\\\ 0&0&1\end{array}\right]$ (9)

The population concentration $C_{0}$ of asteroids encountering with the Moon
is defined as the distribution of the relative number of asteroids that
encounter with the Moon within a unit time and it can be determined by their
velocities (e.g., Figure 5 of Le Feuvre & Wieczorek 2008). In Eq. (3) and (4),
$\vec{v_{p}}$ and $\vec{e_{z}}$ are the asteroids’ encounter velocity in the
geocentric ecliptic coordinate system and an unit vector parallel to the
positive Z-axis respectively. $v_{p}$ is the average encounter velocity of the
asteroids related to the Earth.

$\displaystyle C_{0}$ $\displaystyle=p(\vec{v_{p}})=p(\lambda_{p},\phi_{p})$
(10) $\displaystyle\vec{v_{p}}$
$\displaystyle=v_{p}R_{z}(\frac{\pi}{2}+\lambda_{p})R_{x}(\frac{\pi}{2}-\phi_{p})\vec{e_{z}}$
(11)

In this model, $\vec{v}_{p}$ is determined by $\lambda_{p}$, $\phi_{p}$, and
$v_{p}$. The concentration of asteroids encountering with the Moon is assumed
unaffected by the orbital obliquity and inclination of the Moon (see the
appendix A). Eq. (A.20) indicates the concentration of asteroids encountering
with the Moon can be estimated by the concentration encountering with the
Earth. Then $C_{0}$ should be function of $(v_{p},\lambda_{p},\phi_{p})$.
Spectrum of $v_{p}$ is not considered in this study and it is set as the
average encounter speed (Horedt & Neukum 1984; Zahnle et al. 2001). Then
$C_{0}$ is independent of $v_{p}$ and a function of $(\lambda_{p},\phi_{p})$.
Since the precession of lunar orbit, the asteroids’ azimuth distribution will
not affect the cratering asymmetries. Therefore, only the marginal
distribution of Eq. (3) $\int_{-\pi}^{\pi}C_{0}d\lambda_{p}$ is required in
the calculation of cratering asymmetries. This marginal distribution is taken
from Le Feuvre & Wieczorek (2008) and it has been shown in Figure 6 of Le
Feuvre & Wieczorek (2008).

Figure 1: The coordinate systems used in calculation. The geocentric ecliptic
coordinate system is $OX_{1}Y_{1}Z_{1}$. The lunar fixed coordinate system is
$OX_{2}Y_{2}Z_{2}$ and its origin is translated to the Earth. The gray plane
$C_{1}$ is the ecliptic plane. The plane $C_{2}$ is the lunar equatorial
plane. $\overline{OA}$ is the intersection of $C_{1}$ and $C_{2}$.

The Moon is assumed to be synchronously rotating with a constant angular
velocity and the prime meridian is determined by the mean sub-Earth point
(GSFC 2008). In the lunar fixed coordinate system whose X-axis is the
intersection of the lunar equator plane and prime meridian plane, and Z-axis
is the lunar spin axis, the normal at lunar surface is $\vec{n}(\lambda,\phi)$
and the transformation matrix $T$ from this lunar fixed coordinate system to
the geocentric ecliptic coordinate system are determined by $\omega_{2}$,
$i_{2}$ and $M$. The relationship between coordinate systems used in this
section is illustrated in Figure 1.

$\displaystyle\vec{n}(\lambda,\phi)$
$\displaystyle=R_{z}(\frac{\pi}{2}+\lambda)R_{x}(\frac{\pi}{2}-\phi)\vec{e_{z}}$
(12) $\displaystyle T$
$\displaystyle=R_{z}(\frac{\pi}{2}+\omega_{2})R_{x}(i_{2})R_{z}(M+M_{0})$ (13)

In Eq. (6), $M_{0}$ is a parameter related to the position of the mean sub-
Earth point and it will be determined by Eqs. (7-9). When the Moon at the
perigee, the center of Earth passes through the lunar prime meridian plane.

$\displaystyle\frac{\vec{r_{m}}}{|\vec{r_{m}}|}\Big{|}_{f_{m}=0}+T\cdot\vec{n}(0,\phi)=0$
(14)

Solving Eq. (14), this gives

$\displaystyle\cos{M_{0}}$
$\displaystyle=\frac{\cos{i_{1}}\sin{(\omega_{1}-\omega_{2})}\sin{\omega_{3}}-\cos{(\omega_{1}-\omega_{2})}\cos{\omega_{3}}}{\sqrt{1-(\sin{i_{2}}\sin{(\omega_{1}-\omega_{2})}\cos{\omega_{3}}-\cos{i_{1}}\sin{i_{2}}\cos{(\omega_{1}-\omega_{2})}\sin{\omega_{3}}-\sin{i_{1}}\cos{i_{2}}\sin{\omega_{3}})^{2}}}$
(15) $\displaystyle\sin{M_{0}}$
$\displaystyle=\frac{-\cos{i_{1}}\cos{i_{2}}\cos{(\omega_{1}-\omega_{2})}\sin{\omega_{3}}-\cos{i_{2}}\sin{(\omega_{1}-\omega_{2})}\cos{\omega_{3}}+\sin{i_{1}}\sin{i_{2}}\sin{\omega_{3}}}{\sqrt{1-(\sin{i_{2}}\sin{(\omega_{1}-\omega_{2})}\cos{\omega_{3}}-\cos{i_{1}}\sin{i_{2}}\cos{(\omega_{1}-\omega_{2})}\sin{\omega_{3}}-\sin{i_{1}}\cos{i_{2}}\sin{\omega_{3}})^{2}}}$
(16)

###  2.2 Distribution of Impact Flux

For given $\vec{v_{p}}$, $\vec{r_{m}}$ and $\vec{v_{m}}$, the velocity of
asteroids relative to the Moon is $\vec{v_{p}}-\vec{v_{m}}$. Define $\hat{x}$
as
$\\-\cfrac{\vec{v_{p}}-\vec{v_{m}}}{|\vec{v_{p}}-\vec{v_{m}}|}\cdot(T\vec{n}(\lambda,\phi))$.
The impact flux $\delta F$ is defined as the distribution of the number of
asteroids that impact on the lunar surface within a unit area and a unit time.
According to Eq. (26) of Wang & Zhou (2016) and Eq. (A.47) of Le Feuvre &
Wieczorek (2011), the impact flux $\delta F$ is a function of $\hat{x}$. In
Eq. (13), $M_{m}$ and $R_{m}$ are the mass and radius of the Moon
respectively.

$\displaystyle\delta F$
$\displaystyle=C_{0}|\vec{v_{p}}-\vec{v_{m}}|f(\hat{x})$ (17) $\displaystyle
f_{1}(\hat{x})$
$\displaystyle=\left\\{\begin{array}[]{lcl}\hat{x}&,&\hat{x}\geq 0\\\
0&,&\hat{x}<0\end{array}\right.$ (20) $\displaystyle f_{2}(x)$
$\displaystyle=\left\\{\begin{array}[]{lcl}\frac{1}{4}(1+\Gamma)^{-1}(1+\mu^{-1})(\Gamma+(1+\mu)\hat{x})&,&\hat{x}\geq\frac{-\Gamma}{2+\Gamma}\\\
0&,&\hat{x}<\frac{-\Gamma}{2+\Gamma}\end{array}\right.$ (23)
$\displaystyle\mu=$ $\displaystyle\sqrt{1+\cfrac{2\Gamma}{1+\hat{x}}}\
,\Gamma=\cfrac{2GM_{m}}{|\vec{v_{p}}-\vec{v_{m}}|^{2}R_{m}}$ (24)

where $f_{1}(\hat{x})$ and $f_{2}(\hat{x})$ are two forms of $f(\hat{x})$.
$f_{1}(\hat{x})$ is from Wang & Zhou (2016) which assumes the trajectories of
asteroids are straight lines in the direction of their common encounter
velocity. While $f_{2}(\hat{x})$ is from Le Feuvre & Wieczorek (2011) in which
trajectories of asteroids are treated as hyperbolic curve with a focus at the
center of the Moon. Because $\Gamma<0.02$, we can expand Eq. (23) around 0
with Taylor series.

$\displaystyle
f_{2}(\hat{x})=\hat{x}+(\frac{1}{2}-2\hat{x})\Gamma+\frac{4\hat{x}^{3}+6\hat{x}^{2}-3}{4(1+\hat{x})^{2}}\Gamma^{2}+o(\Gamma^{3})\
\ \hat{x}>\frac{-\Gamma}{2+\Gamma}$ (25)

Obviously, $f_{1}(\hat{x})$ is the first order approximation of
$f_{2}(\hat{x})$. The absolute relative difference between $f_{1}$ and $f_{2}$
is less than $\%3.5$. For simplicity, we use $f(\hat{x})=f_{1}(\hat{x})$ in
following calculations. Then the average flux within a period of lunar orbit
is

$\displaystyle F=\frac{1}{2\pi}\int_{0}^{2\pi}\delta
FdM=\frac{1}{2\pi}\int_{0}^{2\pi}C_{0}|\vec{v_{p}}-\vec{v_{m}}|f(\hat{x})dM$
(26)

It is known today that $\omega_{3}$ changes with a period of 8.85 years and
$\omega_{2}$ changes with a period of 18.61 years. The secular average flux
independent of them is

$\displaystyle\overline{F}(\lambda,\phi;a_{m},e,i_{1},i_{2},v_{p})=\int_{0}^{2\pi}\frac{d\omega_{3}}{2\pi}\int_{0}^{2\pi}\frac{d\omega_{2}}{2\pi}\int_{0}^{2\pi}\frac{d\lambda_{p}}{2\pi}\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}F\frac{d\phi_{p}}{\pi}$
(27)

###  2.3 Distribution of Normal Impact Speed and Cratering Rate

In this model, the impact angle of asteroids at lunar surface is $\theta$ (Eq.
(A.54) of Le Feuvre & Wieczorek 2011) and the normal impact speed is
$V_{\perp}$ (Eqs. (A.50-A.51) of Le Feuvre & Wieczorek 2011).

$\displaystyle V_{\perp}$
$\displaystyle=|\vec{v_{p}}-\vec{v_{m}}|\sqrt{1+\Gamma}\sin{\theta}=|\vec{v_{p}}-\vec{v_{m}}|g(\hat{x})$
(28)

$g(\hat{x})$ can also be written as two different forms $g_{1}(\hat{x})$ and
$g_{2}(\hat{x})$. $g_{1}(\hat{x})$ is from Wang & Zhou (2016) and
$g_{2}(\hat{x})$ is from Le Feuvre & Wieczorek (2011).

$\displaystyle g_{1}(\hat{x})$ $\displaystyle=\hat{x}/\sqrt{1+\Gamma}$ (29)
$\displaystyle g_{2}(\hat{x})$
$\displaystyle=\sqrt{1+\Gamma-(\frac{1+\mu}{2})^{2}(1-\hat{x}^{2})}$
$\displaystyle=|\hat{x}|+\frac{1}{2}sgn(\hat{x})\Gamma-\frac{sgn(\hat{x})}{4(1+\hat{x})}\Gamma^{2}+o(\Gamma^{3})$
(30)

Similar to $f(\hat{x})$, we expand $g_{2}(\hat{x})$ around $\Gamma=0$.
$g_{1}(\hat{x})$ is the first order approximation of $g_{2}(\hat{x})$.
Substituting $g_{1}(\hat{x})$ in to Eq. (17). We obtain the average normal
speed

$\displaystyle\overline{V}_{\perp}(\lambda,\phi;a_{m},e,i_{1},i_{2},v_{p})=\frac{1}{\overline{F}(\lambda,\phi;a_{m},e,i_{1},i_{2},v_{p})}$
$\displaystyle\cdot\int_{0}^{2\pi}\frac{d\omega_{3}}{2\pi}\int_{0}^{2\pi}\frac{d\omega_{2}}{2\pi}\int_{0}^{2\pi}\frac{d\lambda_{p}}{2\pi}\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\frac{d\phi_{p}}{\pi}\int_{0}^{2\pi}\delta
FV_{\perp}\frac{dM}{2\pi}$ (31)

Combing Eq. (16) and Eq. (20), we finally obtain the cratering rate
expression. Similar to Eq. (56) of Wang & Zhou (2016) and applying the scaling
law of crater diameters(e.g., Holsapple & Housen (2007)), the cratering rate
in our model is

$\displaystyle
N_{c}(\lambda,\phi;a_{m},e,i_{1},i_{2},v_{p})\propto(\overline{V}_{\perp}(\lambda,\phi;a_{m},e,i_{1},i_{2},v_{p}))^{\gamma_{p}\alpha_{p}}\overline{F}(\lambda,\phi;a_{m},e,i_{1},i_{2},v_{p})$
(32)

Here the cratering rate calculation only takes account the the near-Earth
objects: $\gamma_{p}\alpha_{p}=0.987$ (Bottke et al. 2002; Holsapple & Housen
2007). $\alpha_{p}$ is an exponent in the cumulative size distribution of
near-Earth objects diameter (Bottke et al. 2002). $\gamma_{p}$ is a parameter
in the scaling law (Holsapple & Housen 2007; Le Feuvre & Wieczorek 2011; Wang
& Zhou 2016).

##  3 Result

In this section, we describe the cratering rate asymmetry produced by Eq.
(21). Section 3.1 demonstrates the spatial variation of impact flux, normal
speed, and cratering rate. Section 3.2 reveals the influences of orbital
obliquity and inclination of the Moon on lunar cratering rate asymmetry.
Section 3.3 provides the evolution of apex/ant-apex ratio with orbital
obliquity and inclination of the Moon.

###  3.1 Spatial Variations of Impact Flux, Normal Speed, and Cratering Rate

First, our derived formula is used to calculate the cratering rate spatial
variation at current values of the Earth-Moon system since such variation can
be compared with previous predictions by Le Feuvre & Wieczorek (2011) and Wang
& Zhou (2016). Figure 2 shows the relative spatial variation of impact flux,
normal speed, and cratering rate on the Moon with parameters set at current
values of the Earth-Moon system. The parameters involved in Eq. (21) are set
as
$(a_{m},e,i_{1},i_{2},v_{p})=(60R_{e},0.0549,5.145^{\circ},1.535^{\circ},19km/s)$
($R_{e}$ is the radius of the Earth). In Figure 2, the relative cratering rate
is symmetry about $0^{\circ}N$. This symmetry arises from the symmetry of the
asteroids’s concentration $C_{0}$. The maximum of impact flux occurs at
$(90^{\circ}W,0^{\circ}N)$ and the minimum is at $(90^{\circ}E,\pm
65^{\circ}N)$. The maximum/minimum ratio of impact flux is 1.24. The maximum
of normal speed occurs at $(90^{\circ}W,0^{\circ}N)$ and the minimum is at
$(90^{\circ}E,\pm 47^{\circ}N)$. The maximum normal speed is 13.7 $km/s$ and
the minimum is 12.1 $km/s$. The maximum of cratering rate occurs at
$(90^{\circ}W,0^{\circ}N)$ and the minimum is at $(90^{\circ}E,\pm
53^{\circ}N)$. The maximum/minimum cratering rate ratio is 1.40. The apex/ant-
apex ratio (the cratering rate ratio between $(90^{\circ}W,0^{\circ}N)$ and
$(90^{\circ}E,0^{\circ}N)$) is 1.36 and this ratio is a measure of the
longitudinal variation. The pole/equator ratio is 0.87 and this is a measure
of the latitudinal variation. The impact flux, normal speed, and cratering
rate with $e=0$ (other parameters are set same as Figure 2) are also
calculated. The relative difference of cratering rate ($e=0$) from Figure 2 is
less than 0.2%.

Figure 2: Distribution of impact flux (a), normal speed (b), and cratering
rate (c) on the Moon for the current lunar orbital obliquity, inclination and
Earth-Moon distance. The maximum is set to 1.00.

###  3.2 Influences of Orbital Inclination and Obliquity of the Moon

We next investigate the specific effects of lunar orbital inclination and
obliquity on apex/ant-apex and pole/equator ratios. Figure 3 shows the
apex/ant-apex ratio $r_{1}$ and pole/equator ratio $r_{2}$ with different
lunar inclination and obliquity. $i_{2}$ is not same with the lunar obliquity
to its orbit normal. For Cassini state 2 ($\omega_{1}=\omega_{2}+\pi$), lunar
obliquity relative to the lunar orbit normal is $i_{1}+i_{2}$, while for
Cassini state 1 ($\omega_{1}=\omega_{2}$), lunar obliquity relative to the
lunar orbit normal is $|i_{1}-i_{2}|$ (Ward 1975). Other parameters in Eq.
(21) are set as $(a_{m},e,v_{p})=(60R_{e},0.0549,19km/s)$. $r_{1}$ decreases
with both $i_{1}$ and $i_{2}$ increased, while $r_{2}$ increases with increase
in $i_{2}$ and seems to be independent of $i_{1}$. Based on our calculated
$i_{1}$ and $i_{2}$ distribution, we speculate the correlation of $r_{1}$ or
$r_{2}$ with $\cos{(i_{1}+i_{2})}$ and $\cos{(2i_{2})}$ can be fitted as
linear regressions.

$\displaystyle r_{1}$
$\displaystyle=a_{11}+a_{12}\cos{(i_{1}+i_{2})}+a_{13}\cos{(2i_{2})}$ (33)
$\displaystyle r_{2}$
$\displaystyle=a_{21}+a_{22}\cos{(i_{1}+i_{2})}+a_{23}\cos{(2i_{2})}$ (34)

When $a_{m}=60R_{e}$, the fitting result is

$\displaystyle\left(\begin{matrix}a_{11}&a_{12}&a_{13}\\\
a_{21}&a_{22}&a_{23}\end{matrix}\right)=\left(\begin{matrix}1.12676&0.2469790&0.0089461\\\
0.97137&-0.0029778&-0.0930123\end{matrix}\right)$ (35)

When $i_{1}$ and $i_{2}$ are between $0^{\circ}$ and $45^{\circ}$, the
relative error between fitting result and Figure 3 is less than 2.4% for
$r_{1}$ and 0.15% for $r_{2}$.

Figure 3: The apex/ant-apex (a) and pole/equator (b) ratios with orbital
obliquity and inclination of the Moon. Only the Cassini state 2
($\omega_{1}=\omega_{2}+\pi$) is calculated.

###  3.3 Evolution of the Apex/ant-apex Ratio

The past obliquity has been very high and the lunar inclination is also
different from current value (Ward 1975; Ćuk et al. 2016). We obtain the
evolution of the lunar orbital obliquity and inclination with the Earth-Moon
distance by reproducing the semi-analytical method for the lunar orbital
evolution from Ćuk et al. (2016). This method includes solving the
differential equations of lunar synchronous orbit controlled by Earth and Moon
tidal dissipation, as well as coupling them with the equation to satisfy the
Cassini state Ward (1975). The solutions show that the lunar inclination damps
from the initial high value to its present low value $5.1^{\circ}$ due to
tidal dissipation, and the lunar obliquity first increases and then decreases
to current value $1.5^{\circ}$ with the jump between 29.7$R_{e}$ and 35$R_{e}$
due to the transitions from Cassini state 1 to Cassini state 2, which is
similar to the extended data Figure 1 in Ćuk et al. (2016). We next apply this
evolution in our model to estimate the evolution of apex/ant-apex ratio
(Figure 4). According to Ćuk et al. (2016), the Moon is in non-synchronous
rotation from 29.7$R_{e}$ to about 35$R_{e}$ (gray box in Figure 4). When the
Moon is at Cassini state 1 (the Earth-Moon distance $<29.7R_{e}$), the
apex/ant-apex ratio decreases with $a_{m}$. When the Moon is at Cassini state
2 ($>35R_{e}$), this ratio reaches a maximum between 40$R_{e}$ and 45$R_{e}$.

Figure 4: Evolution of the apex/ant-apex ratio with Earth-Moon distance. The
X-axis represents the Earth-Moon distance and Y-axis represents the apex/ant-
apex ratio. The black solid line is $1.12e^{-0.0529a_{m}/R_{e}}+1.32$ from Le
Feuvre & Wieczorek (2011). The value of $\alpha_{p}\gamma_{p}$ in Le Feuvre &
Wieczorek (2011) ranges from $0.907$ to $1.25$. The other two dashed black
curves represent Eq. (124) of Wang & Zhou (2016) with
$\alpha_{p}\gamma_{p}=0.987$. The red and blue triangles represent results
from Eq. (21) which uses constant obliquity and inclination same as current
values $(i_{1},i_{2})=(5.145^{\circ},1.535^{\circ})$. The red and blue lines
is calculated basing on the evolution of orbital obliquity and inclination of
the Moon from Ćuk et al. (2016) with lunar tidal dissipation number
$Q_{M}=38$.

###  3.4 Cratering Rate Distribution of 3:2 Resonance

Our formulas can also predict cratering rate distributions with various spin-
orbit resonance. When the resonance is 3:2 (applicable to Mercury, Colombo
1965), we have a different transformation matrix in Eq. (6) as

$\displaystyle T^{{}^{\prime}}$
$\displaystyle=R_{z}(\frac{\pi}{2}+\omega_{2})R_{x}(i_{2})R_{z}(\frac{3}{2}M+M_{0})$
(36)

Also because this resonance is 3:2, a full integration interval for Eqs. (16),
(20), and (21) is extended to two periods $(0,4\pi]$. If setting up the
parameters involved in Eq. (21) as those for Mercury
$(\sqrt{GM_{e}/a_{m}},e,i_{1},i_{2},v_{p})=(48.0km/s,0.205,7.0^{\circ},7.0^{\circ},42.2km/s)$.
Further substituting T in Eqs. (7-21) with T’ and using the asteroids
inclination distribution from Le Feuvre & Wieczorek (2008), the maximum and
minimum of cratering rate with 3:2 resonance is at $(\pm
90^{\circ}E,0^{\circ}N)$ and $\pm 90^{\circ}N$, respectively. The
maximum/minimum cratering rate ratio is 3.64. When the orbital eccentric is
degraded to 0.0 (Le Feuvre & Wieczorek 2008; Wang & Zhou 2016), the maximum
and minimum are at $0^{\circ}N$ and $\pm 90^{\circ}N$, respectively. The
maximum/minimum cratering rate ratio is 2.91.

##  4 Discussion

###  4.1 Comparison with Previous Results

In Figure 2(c), this study gives a similar current cratering rate spatial
variation ($a_{m}=60R_{e}$) as Le Feuvre & Wieczorek (2011) in which the
maximum and minimum appear at $(90^{\circ}W,0^{\circ}N)$ and $(90^{\circ}E,\pm
65^{\circ}N)$ respectively. The difference in the location of the minimum
between our result and Le Feuvre & Wieczorek (2011) may be brought by $f_{1}$
and $g_{1}$ used in our calculations. The apex/ant-apex ratio for current Moon
from Le Feuvre & Wieczorek (2011) is 1.37. The apex/ant-apex ratio for current
Moon from Wang & Zhou (2016) with $v_{p}=19km/s,\alpha_{p}\gamma_{p}=0.987$
and this study are 1.32 and 1.36 respectively. As an extension based on Le
Feuvre & Wieczorek (2011) and Wang & Zhou (2016), this study gives a value
between them. The larger relative difference between this study and Wang &
Zhou (2016) is probably caused by either the asteroids inclination or the
orbital obliquity and inclination of the Moon. The pole/equator ratio in
Figure 2(c) is 0.87. This value is higher than 0.80 in Le Feuvre & Wieczorek
(2011). This difference may be brought by $f_{1}$ and $g_{1}$ used in our
calculations. Besides the current cratering rate, this study also gives the
evolution of the apex/ant-apex ratio in Figure 4. The results from Le Feuvre &
Wieczorek (2011) and Wang & Zhou (2016) are also included in Figure 4. The
value of $\alpha_{p}\gamma_{p}$ adopted in Le Feuvre & Wieczorek (2011) is
about $0.907\sim 1.25$, because they used different parameters in crater
scaling law (in non-porous gravity scaling regime $\gamma_{p}=0.564$ while in
porous $\gamma_{p}=0.410$) and a 10th-order polynomial to fit the size
distribution of asteroids ($\alpha_{p}\approx 2.22$). The influences of
$\alpha_{p}\gamma_{p}$ is shown in Figure 5. When $\alpha_{p}\gamma_{p}$ is
between $0.907\sim 1.25$ and $v_{p}$ is between $19\sim 20km/s$, the apex/ant-
apex ratio for current Moon calculated by this model is about $1.33\sim 1.41$.
Although the value of $\alpha_{p}\gamma_{p}$ or $v_{p}$ in this study is
different from Le Feuvre & Wieczorek (2011), if the orbital obliquity and
inclination of the Moon is assumed constant and same with current values
$(i_{1},i_{2})=(5.145^{\circ},1.535^{\circ})$, this study will reproduce the
results predicted by Le Feuvre & Wieczorek (2011) (red and blue triangle in
Figure 4).

If we consider the variation of obliquity and inclination, when the Earth-Moon
distance is more than $\sim 42R_{e}$, this model gives a consistent value with
Le Feuvre & Wieczorek (2011) and a similar trend with Wang & Zhou (2016).
However, when Earth-Moon distance is between $35R_{e}$ and $\sim 42R_{e}$,
this study gives an opposite trend to previous results. According to Ćuk et
al. (2016), the orbital obliquity and inclination of the Moon decrease in this
interval. This evolution trend of apex/ant-apex ratio can be explained by the
influences of orbital obliquity and inclination of the Moon. When Earth-Moon
distance is between $29.7R_{e}$ and $35R_{e}$, the Moon is in non-synchronous
rotation and the apex/ant-apex ratio will be diminished by non-synchronous
rotation. When Earth-Moon distance is less than $29.7R_{e}$, the apex/ant-apex
ratio is calculated under the assumption: the inclination distribution of
asteroids’ velocity is same as current distribution. Although the obliquity
and inclination is very high, the apex/ant-apex ratio is consistent with Le
Feuvre & Wieczorek (2011). We note that the population of asteroids is
dominated by main-belt asteroids during the late heavy bombardment and and
near-Earth objects since $3.8-3.7$ Ga according to Strom et al. (2015). The
population of near-Earth objects have been in steady state for the past $\sim
3$ Ga (Bottke et al. 2002). The evolution of apex/ant-apex ratio for the
Earth-Moon distance $<29.7R_{e}$ may be quite different from that shown in
Figure 4. It is confirmed that the influences of orbital obliquity and
inclination of the Moon are not negligible in analysing lunar cratering
asymmetry.

Figure 5: Apex/ant-apex ratio of cratering rate with different
$\alpha_{p}\gamma_{p}$ and $v_{p}$. Other parameters involved in Eq. (21) is
set as $(a_{m},e,i_{1},i_{2})=(60R_{e},0.0549,5.145^{\circ},1.535^{\circ})$.

###  4.2 Explanation for the Influences of Orbital Obliquity and Inclination
of the Moon

Figure 6:  Sketch for our model in the lunar fixed coordinate system
$OX_{2}Y_{2}Z_{2}$. The origin $O$ is at the center of the Moon. The
$Z_{2}$-axis is parallel to the lunar spin axis $\vec{s}$ and the $X_{2}$-axis
points to the mean sub earth point. $\vec{k}$ is the ecliptic normal.
$\vec{n}$ is the lunar orbit normal. Point A, B, and C are on the intersection
of the lunar surface and the plane $OY_{2}Z_{2}$. $\overline{OB}$ is parallel
to $\vec{v_{m}}$ and perpendicular to $\vec{n}$. Point A is on the lunar
equator plane $C_{2}$ and it is the ant-apex point $(90^{\circ}E,0^{\circ}N)$.
Point C is on the ecliptic plane. Point D is the mean sub-earth point.

When the lunar orbit eccentric $e=0.0$, our model is sketched in Figure 6. In
a lunar rotation period, the relative position between
$\\{\vec{n},\vec{k},\vec{s}\\}$ and the coordinate system $OX_{2}Y_{2}Z_{2}$
is not fixed. The apex or ant-apex point is on the gray cycle $C_{2}$. When
$i_{1}=i_{2}=0$, the gray circle $C_{2}$, red circle, and yellow circle
coincide and $r_{1}$ reaches its maximum. The influence of asteroids
inclination is related to the length of
$\stackrel{{\scriptstyle\frown}}{{AC}}$. The influence of lunar velocity is
related to the length of $\stackrel{{\scriptstyle\frown}}{{AB}}$. The farther
the ant-apex is from $B$ or $C$, the smaller the leading/trailing asymmetry.
In our model, the angular distance between A (ant-apex) and B is in
$[0,max\\{|i_{1}+i_{2}|,|i_{1}-i_{2}|\\}]$ and the angular distance between A
and C is in $[0,i_{2}]$. For Cassini state 2,
$max\\{|i_{1}+i_{2}|,|i_{1}-i_{2}|\\}=|i_{1}+i_{2}|$. This is consistent with
the fitting result: $r_{1}$ is proportional to $\cos{(i_{1}+i_{2})}$ and
$\cos{(2i_{2})}$. For $r_{2}$ which is related to the angular distance between
point D and the red or yellow plane. When $i_{1}=i_{2}=0$, $r_{2}$ reaches its
minimum. The angular distance between point D and the red circle is in
$[0,i_{2}]$ and the angular distance between point D and the yellow circle is
in $[0,max\\{|i_{1}+i_{2}|,|i_{1}-i_{2}|\\}$. The pole/equator asymmetry will
decrease with increasing in those two angular distances. This is also
consistent with the fitting result: $r_{2}$ is proportional to
$\cos{(2i_{2})}$ and $\cos{(i_{1}+i_{2})}$.

###  4.3 Generalization of this Model

The orbital obliquity and inclination of the Moon, Earth-Moon distance, lunar
orbital eccentric, and lunar rotation speed have been included in this model.
In addition to the Moon, this model can be applied to other planets and moons,
especially for the types of spin-orbit resonance. For example, Mercury is
tidally locked with the Sun in a 3:2 resonance. The cratering rate
distribution of 3:2 Resonance is detailed in the Figure 7 Figure 7(a) shows
the distribution of cratering rate for 3:2 resonance. This cratering rate
asymmetry has been reported in Wieczorek et al. (2012). They predict the
cratering asymmetry maximizes at $(0^{\circ}E,0^{\circ}N)$ and
$(180^{\circ}E,0^{\circ}N)$ and minimizes at $(\pm 90^{\circ}E,0^{\circ}N)$.
Both of this study and Wieczorek et al. (2012) predict the distance between
maxima of cratering asymmetry is $180^{\circ}$. The difference between this
study and Wieczorek et al. (2012) may be from the ignorance of the non-
uniformity on the azimuth of asteroids velocity in this study and different
definition of prime meridian between this study and Wieczorek et al. (2012).
In Figure 7(b), the cratering with orbital eccentric $e=0$ shows a different
distribution from $e=0.205$ for 3:2 resonance. The longitudinal variation of
cratering rate is diminished by rotation of planets and moons with eccentric
$e=0$. However, in the cratering rate on the Moon, the difference caused by
eccentric is less than 0.2%. The influences of eccentric is probably related
to the types of spin-orbit resonance and will be investigated in a future
work.

Figure 7: The relative cratering rate for 3:2 resonance. The maximum is set to
1.00. longitudes $0^{\circ}$ and $180^{\circ}$ are subsolar points when $M=0$.
In subfigure (a), the orbital eccentric is 0.205. In subfigure (b), the
orbital eccentric is 0.0.

##  5 Conclusion

In this study, we have presented an extension of Wang & Zhou (2016) and Le
Feuvre & Wieczorek (2011) to calculate the lunar cratering asymmetry with high
obliquity and inclination. Different from previous models, this model is also
able to calculate the cratering asymmetry with different Cassini states, and
synchronous rotating speed. This model gives consistent results with previous
with low obliquity and inclination. When the obliquity, inclination and Earth-
Moon distance are at current values, this model gives an cratering asymmetry
maximizing at $(90^{\circ}W,0^{\circ}N)$ and minimizing at $(90^{\circ}E,\pm
53^{\circ}N)$ using the encountering velocity inclination distribution
calculated in Le Feuvre & Wieczorek (2008). The apex/ant-apex ratio of this
asymmetry is 1.36 and the pole/equator ratio is 0.87. In order to calculate
the cratering rate with high obliquity and inclination, we have assumed the
orbital obliquity and inclination of the Moon don’t affect the asteroids
population encountering with the Moon. Increasing the orbital obliquity and
inclination of the Moon reduces the apex/ant-apex ratio. According to the
evolution of orbital obliquity and inclination of the Moon, this model gives
an increasing trend in apex/ant-apex ratio with the Earth-Moon distance
between $[35R_{e},42R_{e}]$. This ratio was predicted decreased with the
increasing of Earth-Moon distance in previous studies. Besides the cratering
rate, this model also gives the spatial variation of impact flux and impact
normal speed. Our results provide the quantitative information for evaluating
and rectifying the lunar cratering chronology.

###### Acknowledgements.

We thank M. A. Wieczorek and W. Fa for instructive discussion at the early
stage of this study. We also thank J.-L. Zhou for constructive and insightful
suggestions on our research focus. We thank careful reviews by two anonymous
reviewers. Computations were conducted on the High-performance Computing
Platform of Peking University and the Pawsey Supercomputing Centre with
funding from the Australian Government and the Government of Western
Australia. Z.Y. is supported by the B-type Strategic Priority Program of the
Chinese Academy of Sciences, Grant No. XDB41000000 and NSFC 41972321. N.Z. is
grateful for NSFC 41674098, CNSA D020205 and the B-type Strategic Priority
Program of the Chinese Academy of Sciences, Grant No. XDB18010104. This
research has made use of data and/or services provided by the International
Astronomical Union’s Minor Planet Center.

## References

  * Bandermann & Singer (1973) Bandermann, L. W., & Singer, S. F. 1973, Icarus, 19, 108
  * Bottke et al. (2002) Bottke, W. F., Morbidelli, A., Jedicke, R., et al. 2002, Icarus, 156, 399
  * Colombo (1965) Colombo, G. 1965, Nature, 208, 575
  * Fassett et al. (2012) Fassett, C. I., Head, J. W., Kadish, S. J., et al. 2012, Journal of Geophysical Research: Planets, 117
  * Gallant et al. (2009) Gallant, J., Gladman, B., & Ćuk, M. 2009, Icarus, 202, 371
  * Greenberg (1982) Greenberg, R. 1982, The Astronomical Journal, 87, 184
  * GSFC (2008) GSFC. 2008, LRO Project White Paper Version, 4, 1
  * Hartmann (1970) Hartmann, W. K. 1970, Icarus, 13, 299
  * Hartmann & Neukum (2001) Hartmann, W. K., & Neukum, G. 2001, Space Science Reviews, 96, 165
  * Hiesinger et al. (2000) Hiesinger, H., Jaumann, R., Neukam, G., & Head, J. W. 2000, Journal of Geophysical Research, 105, 29239
  * Holsapple & Housen (2007) Holsapple, K. A., & Housen, K. R. 2007, Icarus, 191, 586
  * Horedt & Neukum (1984) Horedt, G. P., & Neukum, G. 1984, Icarus, 60, 710
  * Korycansky & Zahnle (2005) Korycansky, D. G., & Zahnle, K. J. 2005, Planetary and Space Science, 53, 695
  * Le Feuvre & Wieczorek (2008) Le Feuvre, M., & Wieczorek, M. A. 2008, Icarus, 197, 291
  * Le Feuvre & Wieczorek (2011) Le Feuvre, M., & Wieczorek, M. A. 2011, Icarus, 214, 1
  * McGill (1977) McGill, G. E. 1977, Geological Society of America Bulletin, 88, 1102
  * Neukum (1984) Neukum, G. 1984, Meteorite bombardment and dating of planetary surfaces, PhD thesis, National Aeronautics and Space Administration, Washington, DC.
  * Neukum et al. (2001a) Neukum, G., Ivanov, B. A., & Hartmann, W. K. 2001a, Space Science Reviews, 96, 55
  * Neukum et al. (1975) Neukum, G., König, B., & Arkani-Hamed, J. 1975, The moon, 12, 201
  * Neukum et al. (2001b) Neukum, G., Oberst, J., Hoffmann, H., Wagner, R., & Ivanov, B. A. 2001b, Planetary and Space Science, 49, 1507
  * Opik (1951) Opik, E. J. 1951, Proc. R. Irish Acad. Sect. A, 54, 165
  * Strom et al. (2015) Strom, R. G., Malhotra, R., Xiao, Z.-Y., et al. 2015, Research in Astronomy and Astrophysics, 15, 407
  * Wang & Zhou (2016) Wang, N., & Zhou, J. L. 2016, Astronomy & Astrophysics, 594
  * Ward (1975) Ward, W. R. 1975, Science, 189, 377
  * Wetherill (1967) Wetherill, G. W. 1967, Journal of Geophysical Research (1896-1977), 72, 2429
  * Wieczorek et al. (2012) Wieczorek, M. A., Correia, A. C. M., Le Feuvre, M., Laskar, J., & Rambaux, N. 2012, Nature Geoscience, 5, 18
  * Wiesel (1971) Wiesel, W. 1971, Icarus, 15, 373
  * Zahnle et al. (2001) Zahnle, K., Schenk, P., Sobieszczyk, S., Dones, L., & Levison, H. F. 2001, Icarus, 153, 111
  * Ćuk & Burns (2004) Ćuk, M., & Burns, J. A. 2004, The Astronomical Journal, 128, 2518
  * Ćuk et al. (2016) Ćuk, M., Hamilton, D. P., Lock, S. J., & Stewart, S. T. 2016, Nature, 539, 402

##  Appendix A Asteroids encountering with the Moon with high obliquity and
inclination

In section 2.1, we assumed the concentration of asteroids encountering with
the Moon is unaffected by the orbital obliquity and inclination of the Moon.
Although the probability of asteroids encountering with the Moon has been
estimated based on (Opik 1951; Wetherill 1967; Greenberg 1982) and it has been
demonstrated the probability encountering with the Moon is similar with the
Earth when the lunar inclination and obliquity is about 0 in Figure 6 of Le
Feuvre & Wieczorek (2008). But when the inclination and obliquity is high, the
rationality of our assumption is uncertain. In this section we introduce a
different framework to prove this assumption. For any asteroid with semi-major
axis, eccentricity, inclination, longitude of ascending node, argument of
perihelion, true anomaly, mean anomaly and eccentric anomaly are
$(a,e,i,\Omega,\omega,f,M,E)$. Here we use subscript $e$ to represent the
orbit of Earth and $m$ to represent the Moon. In the heliocentric ecliptic
coordinates system, the position of this asteroid is

$\displaystyle\vec{r}$
$\displaystyle=rR_{z}(\Omega)R_{x}(i)R_{z}(\omega+f)[1,0,0]^{T}$ (37)
$\displaystyle r$ $\displaystyle=a(1-e^{2})/(1+e\cos{f})=a(1-e\cos{E})$ (38)
$\displaystyle r_{max}=a(1+e)\ ,\ r_{min}=a(1-e)$

For elliptic trajectory,

$\displaystyle M$ $\displaystyle=E-e\sin{E}$ (39) $\displaystyle\cos{f}$
$\displaystyle=\frac{\cos{E}-e}{1-e\cos{E}}$ (40) $\displaystyle\sin{f}$
$\displaystyle=\frac{\sqrt{1-e^{2}}\sin{E}}{1-e\cos{E}}$ (41)

Using the common assumption: uniform precession of $\Omega$ and $\omega$(Opik
1951; Wetherill 1967; Greenberg 1982), for given $(a,e,i)$ the joint
probability density is

$\displaystyle
P_{\Omega,\omega,M}(\Omega,\omega,M|a,e,i)=\left(\frac{1}{2\pi}\right)^{3}$
(42) $\displaystyle P_{\Omega,\omega,E}(\Omega,\omega,E|a,e,i)$
$\displaystyle=P_{\Omega,\omega,M}(\Omega,\omega,M|a,e,i)|\frac{\partial
M}{\partial E}|=\left(\frac{1}{2\pi}\right)^{3}(1+e\sin{E})$ (43)
$\displaystyle\Omega\in[-\pi,\pi],\ \omega\in[-\pi,\pi],\ M\in[-\pi,\pi],\
E\in[-\pi,\pi]$ (44)

The position of this asteroid can also be expressed as $\vec{r}=[x,y,z]^{T}$,
then we obtain a transformation: $(\Omega,\omega,E)\mapsto(x,y,z)$.

$\displaystyle\left\\{\begin{array}[]{cll}x&=&r(\cos{\Omega}\cos{(\omega+f)}-\sin{\Omega}\sin{(\omega+f)}\cos{i})\\\
y&=&r(\sin{\Omega}\cos{(\omega+f)}+\cos{\Omega}\sin{(\omega+f)}\cos{i})\\\
z&=&r\sin{(\omega+f)}\sin{i}\end{array}\right.$ (48)

Eq. (A.9) has 4 solutions: $(\Omega_{k},\omega_{k},E_{k}),k=1,2,3,4$.

$\displaystyle\left\\{\begin{array}[]{cll}\cos{f_{k}}&=&\cfrac{a(1-e^{2})}{e\sqrt{x^{2}+y^{2}+z^{2}}}-\cfrac{1}{e}\\\
\sin{(\omega_{k}+f_{k})}&=&\cfrac{z}{\sqrt{x^{2}+y^{2}+z^{2}}\sin{i}}\\\
\Omega_{k}&=&atan2\Big{(}x\cos{(\omega_{k}+f_{k})}+y\sin{(\omega_{k}+f_{k})}\cos{i},\\\
&&y\cos{(\omega_{k}+f_{k})}-x\sin{(\omega_{k}+f_{k})}\cos{i}\Big{)}\\\
E_{k}&=&atan2\Big{(}\sqrt{x^{2}+y^{2}+z^{2}}\cos{f_{k}}+ae,\\\
&&\sqrt{x^{2}+y^{2}+z^{2}}\sin{f_{k}}\frac{1}{\sqrt{1-e^{2}}}\Big{)}\end{array}\right.$
(54)

Using Eqs. (A.6-A.10), the joint probability density is

$\displaystyle
P_{x,y,z}(x,y,z|a,e,i)=\sum_{k=1}^{4}\left|\left|\cfrac{\partial(\Omega_{k},\omega_{k},E_{k})}{\partial(x,y,z)}\right|\right|P_{\Omega,\omega,E}(\Omega_{k},\omega_{k},E_{k}|a,e,i)$
$\displaystyle=\frac{1}{2a\pi^{3}}\frac{1}{\sqrt{r^{2}\sin^{2}{i}-z^{2}}}\frac{1}{\sqrt{(r-r_{min})(r_{max}-r)}}\
,\ \ r=\sqrt{x^{2}+y^{2}+z^{2}}$ (56)

In Eq. (A.11),
$\cfrac{\partial(\Omega_{k},\omega_{k},E_{k})}{\partial(x,y,z)}$ is the Jacobi
matrix. Eq. (A.11) is valid when $(x,y,z)\in\mathbf{D}=\\{|z|\leq|\sin{i}|r\
and\ r_{min}\leq r\leq r_{max}\\}$. When $(x,y,z)\notin\mathbf{D}$,
$P_{x,y,z}(x,y,z|a,e,i)=0$. This asteroid encountering with a fix point
$\vec{r_{0}}=[x_{0},y_{0},z_{0}]^{T}$ is defined as
$|\vec{r}-\vec{r_{0}}|\leq\tau$ ($\tau$ is different for the Moon and the
Earth) and $\tau\ll min\\{|\vec{r}|,|\vec{r_{0}}|\\}$. Then we obtain the
probability encountering with a fixed point $P_{1}$ and its error $\delta
P_{1}$.

$\displaystyle P_{1}$
$\displaystyle=\iiint_{|\vec{r}-\vec{r_{0}}|\leq\tau}P_{x,y,z}(x,y,z|a,e,i)dxdydz$
$\displaystyle\approx\frac{4}{3}\pi\tau^{3}P_{x,y,z}(x_{0},y_{0},z_{0}|a,e,i)$
(57) $\displaystyle\delta P_{1}$
$\displaystyle\approx\iiint_{|\vec{r}-\vec{r_{0}}|\leq\tau}|\vec{r}-\vec{r_{0}}||\nabla
P_{x,y,z}(x_{0},y_{0},z_{0}|a,e,i)|dxdydz$ $\displaystyle=\pi\tau^{4}|\nabla
P_{x,y,z}(x_{0},y_{0},z_{0}|a,e,i)|$ (58)

Because $P_{x,y,z}$ is not bounded with
$(r^{2}\sin^{2}{i}-z^{2})(r-r_{max})(r-r_{min})=0$. Eq.(A.12) and Eq. (A.13)
are valid when $min\\{|r\sin{i}\pm
z|,|r-r_{max}|,|r-r_{min}|\\}\geq\varepsilon a\ (\varepsilon>0)$. When
$min\\{|r\sin{i}\pm z|,|r-r_{max}|,|r-r_{min}|\\}<\varepsilon a$, the supremum
of $P_{1}$ can be estimated in spherical coordinates,

$\displaystyle P_{1}\leq\iiint_{E}P_{x,y,z}r^{2}|\sin{\theta}|drd\theta
d\varphi=\frac{1}{2a\pi^{3}}\left(\int_{\varphi}d\varphi\right)$
$\displaystyle\cdot
Re\left(\int_{\theta}\frac{\sin{\theta}d\theta}{\sqrt{\sin^{2}{i}-\cos^{2}{\theta}}}\right)Re\left(\int_{r}\frac{rdr}{\sqrt{(r-r_{min})(r_{max}-r)}}\right)$
(59) $\displaystyle
E=\Big{\\{}|r-r_{0}|\leq\tau,|\theta-\theta_{0}|\leq\arctan{\tau/r_{0}},|\varphi-\varphi_{0}|\leq\arctan{\frac{\tau}{r_{0}\cos{\theta_{0}}}}\Big{\\}}$
(60)

This section is only a qualitative explanation. For simplicity, following
derivations are under the condition: $min\\{|r\sin{i}\pm
z|,|r-r_{max}|,|r-r_{min}|\\}\geq\varepsilon a$. When $\vec{r_{0}}$ is not
fixed, for the Earth
$\vec{r_{0}}=\vec{r_{e}}=r_{e}R_{z}(\Omega_{e})R_{x}(i_{e})R_{z}(\omega_{e}+f_{e})[1,0,0]^{T}$,
this asteroid encounters with the Earth by probability $P_{2}$.

$\displaystyle
P_{2}=\iiint_{(\Omega_{e},\omega_{e},M_{e})}P(\Omega_{e},\omega_{e},M_{e}|a_{e},e_{e},i_{e})P_{1}d\Omega_{e}d\omega_{e}dM_{e}$
$\displaystyle=2\pi\iint_{(\omega_{e},M_{e})}P(\Omega_{e}=0,\omega_{e},M_{e}|a_{e},e_{e},i_{e})P_{1}d\omega_{e}dM_{e}$
(61) $\displaystyle\delta
P_{2}=2\pi\iint_{(\omega_{e},M_{e})}P(\Omega_{e}=0,\omega_{e},M_{e}|a_{e},e_{e},i_{e})\delta
P_{1}d\omega_{e}dM_{e}$ (62)

Eq. (A.16) and Eq. (A.17) use the rotational symmetry about the z axis of
$P_{x,y,z}$. For the Moon
$\vec{r_{0}}=\vec{r_{e}}+\vec{r_{m}}=r_{e}R_{z}(\Omega_{e})R_{x}(i_{e})R_{z}(\omega_{e}+f_{e})[1,0,0]^{T}+r_{m}R_{z}(\Omega_{m})R_{x}(i_{m})R_{z}(\omega_{m}+f_{m})[1,0,0]^{T}$
, this asteroid encounters with the Moon by $P_{3}$.

$\displaystyle
P_{3}=2\pi\iint_{(\omega_{e},M_{e})}P(\Omega_{e}=0,\omega_{e},M_{e}|a_{e},e_{e},i_{e})$
$\displaystyle\iiint_{(\Omega_{m},\omega_{m},M_{m})}P(\Omega_{m},\omega_{m},M_{m}|a_{m},e_{m},i_{m})P_{1}d\omega_{e}dM_{e}d\Omega_{m}d\omega_{m}dM_{m}$
(63) $\displaystyle\delta
P_{3}=2\pi\iint_{(\omega_{e},M_{e})}P(\Omega_{e}=0,\omega_{e},M_{e}|a_{e},e_{e},i_{e})$
$\displaystyle\iiint_{(\Omega_{m},\omega_{m},M_{m})}P(\Omega_{m},\omega_{m},M_{m}|a_{m},e_{m},i_{m})\delta
P_{1}d\omega_{e}dM_{e}d\Omega_{m}d\omega_{m}dM_{m}$ (64)

The difference between $P_{2}$ and $P_{3}$ can be estimated by

$\displaystyle|P_{2}\tau_{e}^{-3}-P_{3}\tau_{m}^{-3}|\leq
2\pi\iint_{(\omega_{e},M_{e})}P(\Omega_{e}=0,\omega_{e},M_{e}|a_{e},e_{e},i_{e})\iiint_{(\Omega_{m},\omega_{m},M_{m})}P(\Omega_{m},\omega_{m},M_{m}|a_{m},e_{m},i_{m})$
$\displaystyle\times\frac{4\pi}{3}|P_{x,y,z}(x_{m}+x_{e},y_{m}+y_{e},z_{m}+z_{e}|a,e,i)-P_{x,y,z}(x_{e},y_{e},z_{e}|a,e,i)|d\omega_{e}dM_{e}d\Omega_{m}d\omega_{m}dM_{m}$
(65)

From Eq. (A.20), $P_{3}$ can be estimated by
$P_{2}\frac{\tau_{m}^{3}}{\tau_{e}^{3}}$. While
$P_{2}\frac{\tau_{m}^{3}}{\tau_{e}^{3}}$ is independent of lunar inclination
and obliquity, therefore we can use the concentration of asteroids
encountering with the Moon with low inclination and obliquity to replace the
concentration with high inclination and obliquity. We note that when
$(a_{e},i_{e},e_{e})=(1,0,0)$, for $87\%$ of the near-earth orbits (the
dataset of near-earth orbts is taken from the International Astronomical
Union’s website), the relative error between
$P_{2}\frac{\tau_{m}^{3}}{\tau_{e}^{3}}$ and $P_{3}$ calculated by Eq. (A.20)
is less than 5%.