# Cream Skimming the Underground: Identifying Relevant Information Points from
Online Forums

Felipe Moreno-Vera, Mateus Nogueira, Daniel S. Menasché, Miguel Bicudo, Ashton
Woiwood, Enrico Lovat, Anton Kocheturov, and Leandro Pfleger de Aguiar  Felipe
Moreno-Vera1, Mateus Nogueira1, Cainã Figueiredo1, Daniel S. Menasché1, Miguel
Bicudo1,
Ashton Woiwood2, Enrico Lovat3, Anton Kocheturov3, and Leandro Pfleger de
Aguiar4  1Federal University of Rio de Janeiro (UFRJ), 3Siemens Corporation,
2ESO, 4Amazon.com

###### Abstract

This paper proposes a machine learning-based approach for detecting the
exploitation of vulnerabilities in the wild by monitoring underground hacking
forums. The increasing volume of posts discussing exploitation in the wild
calls for an automatic approach to process threads and posts that will
eventually trigger alarms depending on their content. To illustrate the
proposed system, we use the CrimeBB dataset, which contains data scraped from
multiple underground forums, and develop a supervised machine learning model
that can filter threads citing CVEs and label them as Proof-of-Concept,
Weaponization, or Exploitation. Leveraging random forests, we indicate that
accuracy, precision and recall above 0.99 are attainable for the
classification task. Additionally, we provide insights into the difference in
nature between weaponization and exploitation, e.g., interpreting the output
of a decision tree, and analyze the profits and other aspects related to the
hacking communities. Overall, our work sheds insight into the exploitation of
vulnerabilities in the wild and can be used to provide additional ground truth
to models such as EPSS and Expected Exploitability.

###### Index Terms:

Cybersecurity, online forums, data mining.

##  I Introduction

The exploitation of vulnerabilities in the wild poses significant threats to
the Internet ecosystem, being a concern to end users, companies, and, more
generally, to the stability of the Internet itself. In essence, exploitation
refers to the use of a weaponized exploit to attack a target. In this stage,
the attacker uses the weapon to take advantage of a vulnerability and gain
unauthorized access to a system or steal sensitive information. Therefore,
early detection of weaponization and tentative exploitation is key for
defending against attacks.††This paper appears at IEEE International
Conference on Cyber Security and Resilience (IEEE CSR), 2023. The first two
authors contributed equally to the work. Corresponding authors:
{felipe.moreno.vera<EMAIL_ADDRESS><EMAIL_ADDRESS>
While public databases on weapons, such as ExploitDB, are continuously updated
with information about how to exploit vulnerabilities, underground hacking
forums still contain privileged and more up-to-date information about the
availability and development of exploits and, more importantly, about the
tentative use of those exploits in the wild [4, 5, 14]. In particular, certain
forums contain information about the prices of exploits, and instructions on
how to make attacks Fully UnDetectable (FUD). In this context, figuring out
what users are discussing in those forums is instrumental to detecting and
neutralizing the exploitation of vulnerabilities in the wild. Monitoring the
discussion in these forums also allows for tracking exploit prices, their
usage, demand, and main targets.

We use the CrimeBB dataset, made available by Cambridge Cybercrime Centre,
which contains data scraped from multiple underground forums [14]. We focus on
activity related to hacking, noting that the increasing volume of posts
discussing exploitation in the wild calls for an automatic approach to process
threads that will trigger alarms depending on their content. To that aim, we
developed a supervised machine learning model, which filters threads citing a
Common Vulnerabilities and Exposures (CVE) identifier and labels them as
Proof-of-Concept (PoC), Weaponization, or Exploitation. Then, we indicate
rules that can be automatically derived from data, providing insights into the
difference between weaponization and exploitation.

Prior art.  Weaponization and exploitation are two of the key stages involved
in the development of a cyberattack. Most of the literature has focused on
weaponization [8, 2], i.e., the process of building exploits for
vulnerabilities. Much less attention has been given to exploitation in the
wild, i.e., the actual use of a weaponized exploit to attack a target, or to
gain unauthorized access into a system or steal sensitive information [4]. In
part, this occurs because the study of exploitation in the wild involves
sensitive data and stringent non-disclosure agreements.

Exploit Prediction Scoring System (EPSS) [9] and Expected Exploitability [16]
are two examples of systems that aim at determining exploitability in the
wild. Whereas EPSS uses private sources to derive its parameters, Expected
Exploitability uses public artifacts. However, to the best of our knowledge,
there is no prior work using CrimeBB for the purpose of understanding
exploitation in the wild. Our work serves to close this gap and can be used to
provide additional ground truth to previous models such as EPSS and Expected
Exploitability.

Contributions. In summary, our key contributions are twofold. First, we
provide an analysis of exploitation of vulnerabilities in the wild, using the
CrimeBB dataset. We conduct a longitudinal analysis of profits and other
aspects related to the hacking communities, e.g., indicating the prices
associated with exploits and the distribution of delays between discussions
about vulnerabilities on those forums and the release of information at the
National Vulnerability Database (NVD). Second, we present a classifier for
assessing eminent threats based on underground forums.

Paper structure.  The remainder of this paper is structured as follows. In
Section II we discuss related work and Section III presents our dataset, with
some general statistics. Section IV reports our empirical findings, in Section
V we discuss our thread classifier, and Section VI concludes.

##  II Related Work and Background

In what follows, we discuss related work and background pertaining to the main
themes of our work.

###  II-A NLP and threat intelligence (TI)

The use of Natural Language Processing (NLP) for the analysis of hacker forums
has been considered in [15, 13, 14]. In this work, we complement such body of
literature by focusing on discussions about software vulnerabilities within
CrimeBB forums, which have been previously considered for the analysis of
eWhoring [13] and other cybercrimes [14].

Threat Miner [6] is a system to identify threats based on hacker forums. The
authors of Threat Miner classify notifications or reports as “good” if they
represent a cyber threat that can be linked to a known CVE. In this study, we
focus specifically on analyzing threads within CrimeBB forums that can be
linked to known CVEs. By leveraging CVEs, we relate data from CrimeBB forums
against other sources such as the Common Vulnerability Scoring Systems (CVSS)
and EPSS to gain new insights into the lifecycle of vulnerabilities.

###  II-B Blackhat forums

Blackhat forums comprise unstructured posts. All posts include their content,
author, and subject. Leveraging a public dataset collected by the CrimeBB
project, we analyze posts from multiple blackhat forums. In the posts, we find
references to vulnerabilities, IP addresses, and products that are being
exploited in the wild. Blackhat forums provide a way for researchers as well
as badly-intentioned users to trade knowledge about hacking. These forums
supply information ranging from beginner hacking skills to functional hacking
tools that anyone can easily get access, sometimes for free. The so-called
_script-kiddies_ , i.e., home users with limited computing skills, for
instance, can leverage those tools to initiate cyberattacks. One of our goals
is to distinguish between research activity that poses a potential threat
against exploitation in the wild, wherein criminals chat about threats.

Figure 1: Proposed framework composed of two main steps: (1) threads and posts
pre-processing including feature extraction and labeling; (2) three-class
classifier to classify threads based on their content.

##  III Dataset

Cambridge Cybercrime Centre makes available sixteen underground forums through
CrimeBB. In these CrimeBB underground forums, we have 54,460,134 posts under
5,270,587 discussion threads. Those posts were filtered to extract information
about software vulnerabilities. In this work, we will focus on the largest
forum, Hackforums, which has many boards for exchanging hacking knowledge,
encompassing topics ranging from hacking to games. In what follows, we provide
further details on the dataset used in this paper (Table I). 111All the
material to reproduce the results presented in this paper is available at
https://tinyurl.com/crimebbpaper

###  III-A Producing the dataset

To produce the dataset, we consider the following steps listed in Figure 1.
First, we filter all posts citing at least one CVE (details in Section
III-C1). Each of those posts is contained in a thread. Then, we group all the
posts in each of these threads, along with the thread title. Finally, for each
thread, we proceed with feature extraction. The features correspond to the
presence of words in threads, e.g., through Bag-of-Words (BoW), Term
Frequency-Inverse Document Frequency (TF-IDF), or doc2vec (Section III-C2).
Then, our classifier takes all features extracted from each thread as input
and classifies the thread into one of the considered target classes. By
classifying threads as opposed to individual posts, the proposed approach
provides a greater amount of contextual information to the classifier, which
enhances classification accuracy.

TABLE I: Number of posts (threads) citing CVEs in the top 10 Hackforums
boards, ranked by number of tagged posts

Board | Number of posts (threads) citing CVEs
---|---
| Posts tagged as |
| PoC | Weapon | Exploit | All posts
Pentesting and Forensics | 271 | (55) | 210 | (57) | 11 | (3) | 557 | (166)
Premium Tools and Programs | 198 | (1) | 28 | (3) | 142 | (4) | 433 | (20)
Website and Forum Hacking | 93 | (34) | 139 | (43) | 16 | (12) | 333 | (132)
Hacking Tools and Programs | 10 | (7) | 57 | (28) | 174 | (7) | 260 | (59)
Premium Sellers Section | – | – | 81 | (28) | 89 | (26) | 210 | (66)
Beginner Hacking | 86 | (43) | 58 | (47) | 6 | (6) | 219 | (143)
Botnets, IRC, and Zombies | 24 | (4) | 85 | (34) | 22 | (5) | 160 | (62)
Hacking Tutorials | 58 | (21) | 8 | (4) | 3 | (3) | 74 | (33)
Secondary Sellers Market | 8 | (4) | 33 | (21) | – | – | 91 | (40)
News and Happenings | 9 | (9) | 11 | (5) | 1 | (1) | 75 | (54)
Total, all boards | 757 | (244) | 710 | (397) | 464 | (102) | 3,037 | (1,162)

###  III-B Target classes and manual labeling

Our search for vulnerabilities involved using a case-insensitive regular
expression cve-[0-9]{4}-[0-9]{4,} (slightly more specific than cve(-id)?(?i)
used in [2]) to search for posts referring to vulnerabilities by their CVE
identifiers. Across all CrimeBB forums, we found 4,116 posts citing 1,498
unique CVEs, under 1,700 discussion threads. This aligns with previous
research on marketed exploits, which considered similar quantities of posts
[1, 10, 3]. We discard the second most relevant forum, Antichat, due to the
posts being primarily in Russian. We highlight that Antichat provides
citations to around 90 additional unique CVEs. For Hackforums, we found 3,037
posts explicitly referring to 1,068 unique CVEs, under 1,162 discussion
threads between December 2007 and October 2019 (to be contrasted against 194
discussion threads and around 3,000 posts considered in [2]). From these 1,162
threads, a total of 1,067 were manually labeled by experts. The experts used
the following code book to manually label the threads:222Accounting for slang
and abbreviations that are typical in those communities is left as a subject
for future work.

Figure 2: CDF of hacking tools prices: prices at CrimeBB are relatively low
compared against the Russian market – some prices correspond to subscriptions,
and others to repackaging and FUD. Price statistics: CrimeBB (Min: 1, Median:
100, Max: 4400), Russian market (Min: 100, Median: 2000, Max: 8000)

  * •

PoC:  (1) contain keywords such as PoC, tutorial, guide (given the appropriate
context of producing tools in a lab or controlled environment); (2) provide a
tutorial description about how to build a PoC or (3) discuss vulnerabilities
without signs of using exploits in the wild.

  * •

Weaponization:  (1) contain keywords such as vulnerability and exploit (given
the appropriate context of weaponization); (2) discuss the availability of
fully functional or highly mature exploits, providing references or source
code.

  * •

Exploitation:  (1) mention a well-known hacker group; (2) contain references
to cryptocurrencies and keywords such as bitcoin, exploitation, and attack
(given the appropriate context of attacks in the wild); (3) discuss approaches
to make exploits fully undetectable; or (4) involve markets of exploits.

In addition to the above categories, the experts also labeled a few threads as
Scam, when it was identified the selling of an exploit that was a posteriori
recognized as non-functional. A total of 244, 397, 102 and 10 threads were
labeled as PoC, Weaponization, Exploitation and Scam, respectively (see last
line of Table I). Note that the remainder 314 threads did not fit into any of
the above categories. Scams and the latter threads were not considered in this
study.

Figure 3: CDF of the difference in days between CrimeBB citation and NVD
publish date. Negative values correspond to citations to CVEs that occurred
before NVD published the corresponding vulnerability. Age statistics: CrimeBB
(Min: -396, Median: 132, Max: 7181), Russian market (Min: 1, Median: 95,5,
Max: 2610)

###  III-C Blackhat forums and markets statistics

####  III-C1 General statistics

In Table I, we show the boards that contain most of the posts citing CVEs. In
the top two boards, we find users selling and buying exploits, which indicates
that discussions about vulnerabilities are generally about exploits already
available on the market. Furthermore, Table I also shows the distribution of
posts across different classes over the different boards at Hackforums. Note
that in the board of pentesting, for instance, we find significant activity
related to weaponization and exploitation. In contrast, few posts explicitly
cite CVE identifiers in the board of hacking tutorials.

####  III-C2 Features

The blackhat forums are unstructured. The intrinsic features contained in the
posts are the words that can be encoded using different strategies such as
BoW, TF-IDF, and doc2vec. In this work, we compare the three encoding
techniques, noting that BoW is more interpretable, whereas TF-IDF and doc2vec
yield higher accuracy (see Section IV). 333In future work, we consider
leveraging additional features, such as CVSS and EPSS scores of
vulnerabilities and prices of hacking tools.

####  III-C3 NVD data

We use data from NVD to determine properties of the considered
vulnerabilities, such as severity level (CVSS). In particular, NVD provides a
brief description of each vulnerability, together with its publish date,
products affected, and external resources.

##  IV Empirical Findings

In this section, we report empirical findings from CrimeBB forums, including
exploit prices, delays, and risks.

###  IV-A Prices

We discuss the prices of artifacts cited by users from CrimeBB forums. Figure
2 shows the CDF of prices in dollars. For comparative purposes, we also plot
the CDF of prices of exploits reported at the Russian market studied in [2].
Whereas in CrimeBB forums the minimum, median and maximum values were 1 USD,
100 USD, and 4,400 USD, in the Russian market the corresponding values were
100, 2,000, and 8,000 USD, respectively. Besides, we note that more than 80%
of the references to hacking tools correspond to prices less than 1,000 USD.
The larger prices observed in the Russian market when compared to CrimeBB
forums can be explained by the fact that the Russian market requires explicit
admission by market administrators. Indeed, admission to the market is
conditioned on the user being active in related communities [2]. For this
reason, in the Russian market users tend to discuss more mature, hence more
expensive, artifacts.

In the CrimeBB forums, in contrast, we observed that users tend to propose the
repackaging of already existing exploits, e.g., under new FUD versions [17].
Alternatively, some of the prices refer to subscriptions to websites that tend
to be naturally cheaper than exploits. Despite the differences between prices,
we also observe some similarities. In both platforms, the maximum prices did
not surpass USD 8,000, the majority of prices are below USD 2,000, and roughly
20% of the prices are close to USD 100. Together, those numbers indicate that
the activity in those forums can be monetarily rewarding, with rewards aligned
with most bug bounty programs that offer up to USD 3,000 for a critical bug.
Nonetheless, those numbers are still far from the million-dollar bug bounties
that were recently reported in the
literature.444https://portswigger.net/daily-swig/million-dollar-bug-bounties-
the-rise-of-record-breaking-payouts

###  IV-B Delays

Knowing how long it takes for information about vulnerabilities to appear on
online forums is key, e.g., to assess risks associated with vulnerabilities
and for patch management purposes [7, 12]. In this section, we explore the
delay between the publication of vulnerabilities at NVD and posts appearing at
the forums. For CrimeBB forums, we compute the post age as the difference
between the day of the post and the day on which the corresponding
vulnerability was published at NVD, PostAge = PostPubDate - CVEPubDate.
Similarly, the exploit age reported by [2] is the difference between the day
on which an exploit was published at the Russian market and the day on which
the corresponding vulnerability was published at NVD, ExplAge = ExplPubDate -
CVEPubDate. Figure 3 shows the CDF of PostAge and ExplAge, for CrimeBB forums
and the Russian market, respectively.

Note that more than 50% of exploits discussed in CrimeBB are about
vulnerabilities that were disclosed over the previous 69 days before being
cited at CrimeBB. Considering that 50% of Industrial Control Systems (ICS) are
not patched 60 days after vulnerability disclosure[18], the use of blackhat
forums is imperative to estimate risks associated with vulnerabilities. Among
the similarities between CrimeBB forums and the Russian black market, we
observe that roughly 60% of the activity occurs very close to CVE publish
date. Discussion tends to phase out for virtually all vulnerabilities 6 years
after they are released.

We observe that in the Russian market, we have only positive age values,
whereas under CrimeBB we have a small fraction of negative values. This is
explained by the different nature of the two forums, as discussed in the
previous section: whereas CrimeBB also counts with messages querying about
vulnerabilities and discussing strategies to produce proof-of-concept weapons,
the Russian market contains mostly discussion of mature exploits to be sold at
higher values, and typically being released only after the CVE has already
been published at NVD. Fully functional or high-maturity exploits are rarely
produced before the vulnerability publish date, i.e., ExplPubDate is larger
than CVEPubDate. Discussions about vulnerabilities, however, can initiate
before they are released, as CVE identifiers are announced to the public
before they are published at NVD.

###  IV-C Risks

---
(a) CVSS
(b) EPSS
Figure 4: Distribution of CVSS and EPSS scores across different classes. Note
that 91% of posts refer to CVEs whose CVSS score is higher than the mean CVSS
across all NVD CVEs (not shown in the figure).

Across all vulnerabilities, Figure 4 shows the distribution of CVSS and EPSS
scores conditioned on the thread class (PoC, Weaponization and Exploitation)
and across all threads. Gray boxplots consider one sample per CVE citation
(counting repeated citations to a CVE within posts or threads as distinct
citations), and red boxplots consider a single sample per distinct CVE
(without counting repetitions). Figure 4(a) displays the distribution of CVSS
scores for CVSS versions 2.0 and 3.1 (the latter is available for a subset of
CVEs). The overall median CVSS value exceeds 7.0. Accounting for CVSS 2.0,
exploitation usually corresponds to higher CVSS values when compared against
PoC and weaponization. Additionally, still accounting for CVSS 2.0, CVEs with
higher CVSS values are the most frequently cited. This is evidenced by the
fact that the medians of the red boxplots are all above 9.0, indicating that
the risk scores are magnified by repeated citations. A similar trend is
observed in Figure 4(b) for EPSS. EPSS was released in 2021, and we used its
latest version, as of February 28, 2023, to represent the probability of
exploitation in the wild in the next 30 days. Despite the gap between the post
citations and the release of EPSS scores, we observe that EPSS scores are able
to capture high risks for the vulnerabilities in our dataset.

|
---|---
(a) All | (b) PoC
|
(c) Weaponization | (d) Exploitation
Figure 5: Word clouds showing the most frequent keywords appearing across
posts: (a) all posts; (b) PoC; (c) weaponization and (d) exploitation.

###  IV-D Natural Language Processing (NLP) and word clouds

Figure 5 shows the word clouds obtained from CrimeBB posts grouped by PoC,
exploitation, and weaponization classes and across all posts. The presence of
the word “multisploit” in the PoC category indicates that Trillium Security
MultiSploit Tool is a popular tool used by hackers. The presence of the word
“thanks”, in the exploitation cloud, refers to the behavior of users that
confirm that a given exploitation worked.

Note that the term “clean” is also prevalent across posts about exploitation
in the wild. Indeed, it appears in the context of FUD exploits. Users
typically indicate that they were able to run the exploit and that it was not
detected by any antivirus. As an example, in one post, we found 35 occurrences
of the word “clean” when referring to a “silent” exploit to CVE-2011-3544. The
exploit authors share its MD5 hash and report a detection ratio of 0/35. To
evidence the lack of detection, the authors present a sample output: AVG Free:
Clean; ArcaVir: Clean; Avast 5: Clean. The list continues with 35 antivirus
programs.

##  V Classifiers: Distinguishing Potential Threat Against Eminent Threat

###  V-A Feature extraction and Training of classifier

To encode textual information, we consider three alternatives: BoW, TF-IDF,
and doc2vec. For BoW and TF-IDF, we consider the following four parameters: 1)
the top 30,000 most frequently occurring words, 2) that appear at least 5
times, and 3) in at least 90% of the posts in the corpus are considered for
analysis. In addition, we consider 4) $n$-grams, ranging from one word to
three words. For doc2vec, we encode posts into 5000-dimensional vectors. We
use standard NLP pre-processing techniques, e.g., filtering English language
stop words and punctuation from the posts. Then, we proceed with the use of
decision trees and random forests for classification.

We started with an imbalanced dataset. Recall from the last line of Table I
that our dataset comprises 244 threads labeled as PoC, 397 as weaponization,
and 102 as exploitation. To address this imbalance, we used the Random Over
Sampling heuristic [11] to produce a balanced dataset. We split the dataset
into 75% for training+validation and 25% for testing. We then conducted a grid
search to find the optimal hyperparameters, such as tree depth, the number of
features to consider at each tree split, minimum samples required to split an
internal node, and maximum node degree. To maintain the proportion of
categories during training, we employed 5-fold stratified cross-validation on
the training+validation set.

TABLE II: Decision Tree (DT) and Random Forest (RF) performance.

| Text encoding | Target classes | Accuracy | Precision | Recall | F1
---|---|---|---|---|---|---
DT | BoW | PoC, Weaponization, Exploitation | 0.71 | 0.71 | 0.72 | 0.70
DT | TF-IDF | PoC, Weaponization, Exploitation | 0.73 | 0.73 | 0.74 | 0.72
DT | doc2vec | PoC, Weaponization, Exploitation | 0.74 | 0.74 | 0.74 | 0.73
DT | BoW | Exploitation vs Non-exploitation | 0.85 | 0.86 | 0.85 | 0.85
DT | TF-IDF | Exploitation vs Non-exploitation | 0.91 | 0.91 | 0.91 | 0.91
DT | doc2vec | Exploitation vs Non-exploitation | 0.92 | 0.93 | 0.92 | 0.92
DT | BoW | PoC vs Non-PoC | 0.75 | 0.75 | 0.75 | 0.75
DT | TF-IDF | PoC vs Non-PoC | 0.77 | 0.78 | 0.77 | 0.77
DT | doc2vec | PoC vs Non-PoC | 0.70 | 0.71 | 0.70 | 0.70
DT | BoW | Weaponization vs Non-weapon. | 0.68 | 0.68 | 0.68 | 0.68
DT | TF-IDF | Weaponization vs Non-weapon. | 0.63 | 0.64 | 0.63 | 0.62
DT | doc2vec | Weaponization vs Non-weapon. | 0.59 | 0.59 | 0.59 | 0.59
RF | BoW | PoC, Weaponization, Exploitation | 0.85 | 0.84 | 0.85 | 0.84
RF | TF-IDF | PoC, Weaponization, Exploitation | 0.86 | 0.87 | 0.86 | 0.86
RF | doc2vec | PoC, Weaponization, Exploitation | 0.86 | 0.90 | 0.86 | 0.86
RF | BoW | Exploitation vs Non-exploitation | 0.98 | 0.98 | 0.98 | 0.98
RF | TF-IDF | Exploitation vs Non-exploitation | 0.98 | 0.98 | 0.98 | 0.98
RF | doc2vec | Exploitation vs Non-exploitation | 0.99 | 0.99 | 0.99 | 0.99
RF | BoW | PoC vs Non-PoC | 0.84 | 0.84 | 0.84 | 0.84
RF | TF-IDF | PoC vs Non-PoC | 0.87 | 0.87 | 0.87 | 0.87
RF | doc2vec | PoC vs Non-PoC | 0.88 | 0.90 | 0.88 | 0.87
RF | BoW | Weaponization vs Non-weapon. | 0.67 | 0.67 | 0.67 | 0.67
RF | TF-IDF | Weaponization vs Non-weapon. | 0.69 | 0.70 | 0.69 | 0.69
RF | doc2vec | Weaponization vs Non-weapon. | 0.58 | 0.58 | 0.58 | 0.58

###  V-B Evaluation and interpretation of classifiers

In this section, we report the performance of the considered classifiers. To
that aim, we account for four metrics: accuracy, precision, recall, and F1.
Table II shows the obtained results, considering the best set of
hyperparameters for each configuration, as described above. In particular, our
configurations vary as a function of the text encoding strategy and target
classes. We observed increasing accuracy when switching from the simpler and
more interpretable encoding (BoW) to the most complex but less interpretable
one (doc2vec). Indeed, doc2vec outperforms BoW and TF-IDF except for the cases
“PoC versus Non-PoC” and “Weaponization versus Non- Weaponization”.
Nonetheless, BoW is instrumental to produce the interpretable tree presented
in Figure 6.

With respect to the target classes, we consider PoC vs Weaponization vs
Exploitation and three additional one-against-all classifiers, in which each
binary classifier separates members of a class from members of other classes.
The best results were obtained when filtering exploitation in the wild from
the rest of the threads, which is arguably the first step towards identifying
relevant information at underground forums, as exploitation poses the most
eminent risk. With respect to the classifier model, decision trees are simpler
than random forests, producing less accurate predictions but being amenable to
interpretation, as illustrated below.

Figure 6: Decision tree to classify PoC, weaponization, and exploitation.

Figure 6 illustrates the decision tree used to classify between PoC,
weaponization, and exploitation (first line in Table II). Each internal node
in the tree contains a rule that splits the dataset, and each leaf indicates
the most prominent class at that split and its frequency. Despite the fact
that not all splitting rules that appear in Figure 6 are interpretable, we can
already extract interesting insights from it. In the root of the tree, we find
the rule with the highest splitting power, according to the Gini index
criterion. Indeed, the root together with the leaf immediately below it
indicate that if the thread contains the keyword “poc”, with a 74.2% chance,
it is actually a proof-of-concept. The following rule indicates that posts
wherein users are concerned about privacy, i.e., containing the keyword “pm”,
which stands for “private message” in the black forum jargon, correspond to
exploitation in the wild. Finally, we also observe that JavaScript is a common
language used to produce exploits, e.g., that inject code through unverified
input fields.

##  VI Conclusion

“Data is power as long as you know how to wield it.” In this work, we
leveraged CrimeBB and machine learning methods to learn textual content and
distinguish between: (1) potential threat (proof of concept), (2) eminent
threat (weaponization), and (3) criminals chatting about a threat
(exploitation in the wild). Among our empirical findings obtained by relating
CVSS and EPSS against CrimeBB, we found that the most cited CVEs are typically
related to higher risks and that it is feasible to automatically filter
exploitation threads, with an accuracy above 99%.

We believe that this work opens up interesting avenues for future research,
including the use of transformers such as ChatGPT to distill data from online
forums, and the analysis of additional labels and elements, such as the
maturity level of discussions and its correlation against EPSS scores. We also
aim to expand the number of posts considered in our study, accounting for
vulnerabilities cited in forums by their names, for “named vulnerabilities”,
as opposed to CVE-ids, e.g., Bleichenbacher as opposed to CVE-2018-12404.

## Acknowledgment

This project was sponsored by CAPES, CNPq, and FAPERJ (315110/2020-1,
E-26/211.144/2019 and E-26/201.376/2021).

## References

  * [1] Ablon, L., Libicki, M.C., Golay, A.A.: Markets for cybercrime tools and stolen data: Hackers’ bazaar. Rand Corporation (2014)
  * [2] Allodi, L.: Economic factors of vulnerability trade and exploitation. In: Proceedings of the 2017 ACM SIGSAC conference on computer and communications security. pp. 1483–1499 (2017)
  * [3] Allodi, L., Massacci, F.: Comparing vulnerability severity and exploits using case-control studies. ACM TISSEC 17(1), 1–20 (2014)
  * [4] Basheer, R., Alkhatib, B.: Threats from the dark: a review over dark web investigation research for cyber threat intelligence. Journal of Computer Networks and Communications 2021, 1–21 (2021)
  * [5] Campobasso, M., Allodi, L.: Threat/crawl: a trainable, highly-reusable, and extensible automated method and tool to crawl criminal underground forums. In: APWG eCrime 2022 (2022), arXiv:2212.03641
  * [6] Deguara, N., et al.: Threat miner: A text analysis engine for threat identification using dark web data. In: Big Data. pp. 3043–3052 (2022)
  * [7] Figueiredo, C., Lopes, J.G., Azevedo, R., Vieira, D., Miranda, L., Zaverucha, G., de Aguiar, L.P., Menasché, D.S.: A statistical relational learning approach towards products, software vulnerabilities and exploits. IEEE Transactions on Network and Service Management (2023)
  * [8] Hanks, C., Maiden, M., Ranade, P., Finin, T., Joshi, A., et al.: Recognizing and extracting cybersecurity entities from text. In: ICML (2022)
  * [9] Jacobs, J., Romanosky, S., et al.: Exploit prediction scoring system (EPSS). Digital Threats: Research and Practice 2(3), 1–17 (2021)
  * [10] Kotov, V., Massacci, F.: Anatomy of exploit kits: Preliminary analysis of exploit kits as software artefacts. In: Engineering Secure Software and Systems. pp. 181–196. Springer (2013)
  * [11] Lemaître, G., Nogueira, F., Aridas, C.K.: Imbalanced-learn: A Python Toolbox to Tackle the Curse of Imbalanced Datasets in Machine Learning. Journal of Machine Learning Research 18(17), 1–5 (2017)
  * [12] Miranda, L., Vieira, D., de Aguiar, L.P., Menasché, D.S., Bicudo, M.A., Nogueira, M.S., Martins, M., Ventura, L., Senos, L., Lovat, E.: On the flow of software security advisories. IEEE Transactions on Network and Service Management 18(2), 1305–1320 (2021)
  * [13] Pastrana, S., Hutchings, A., et al.: Measuring ewhoring. In: Proceedings of the Internet Measurement Conference. pp. 463–477 (2019)
  * [14] Pastrana, S., Thomas, D.R., et al.: CrimeBB: Enabling cybercrime research on underground forums at scale. In: Proceedings of the 2018 World Wide Web Conference. pp. 1845–1854 (2018)
  * [15] Rahman, M.R., et al.: What are the attackers doing now? automating cyberthreat intelligence extraction from text on pace with the changing threat landscape: A survey. ACM Computing Surveys (2021)
  * [16] Suciu, O., Nelson, C., Lyu, Z., Bao, T., Dumitraș, T.: Expected exploitability: Predicting the development of functional vulnerability exploits. In: 31st USENIX Security Symposium. pp. 377–394 (2022)
  * [17] Valeros, V., Garcia, S.: Growth and commoditization of remote access trojans. In: 2020 IEEE EuroS&PW. pp. 454–462. IEEE (2020)
  * [18] Wang, B., Li, X., de Aguiar, L.P., et al.: Characterizing and modeling patching practices of industrial control systems. Proceedings of the ACM on Measurement and Analysis of Computing Systems 1(1), 1–23 (2017)


# On an analytical method for the satellite problem revisited

Niccolò Pilloni and Claudio Saccon 111Dipartimento di Matematica, Largo Bruno
Pontecorvo 5, I56127 Pisa, ITALY, e-mail<EMAIL_ADDRESS>
##  1 Introduction

In this note we discuss an analitical method, which goes back to Dirk Brouwer
(1959 see [1]), to compute the trajectories of an artificial satellite
rotating around the Earth when the oblate shape is taken into account. Despite
being somewhat old these kind of techniques are still used nowadays in
software packages like SGP4 (see [5]) used for the tracking of space debris.
Actually the models on which these softwares are based are the result of later
generalizations (see [2, 3, 4]) which include the effect of the atmosferic
drag. We don’t cover the general theory here – we just remark that the problem
with drag is treated as a pertubation of the dragless case.

As well known the Hamiltonian associated with the “non spherical Earth” can be
written as a series of “harmonics” of three types: zonal (only depending on
the latitude), sectorial (only depending on the longitude), and tesseral
(depending on both) each having a suitable coefficient. The most important
term (apart from the zero-th one which corresponds to the perfectly spherical
case) is the coefficient $J_{2}$ of the first zonal harmonic. The main idea in
[1] is looking for a canonical tranformation to a new set of variables such
that the new Hamiltonian only depends on the momenta. If this is possible then
the new Hamilton equations become trivial: the momenta are constant while the
coordinates are linear (this goes with a price: the transformation is usually
complicated and made up of a very long list of terms). To accomplish this a
Brouwer looks for a generating function of the form of a power series in the
coefficient $J_{2}$. Using some calculus involving composition of power series
he derives some conditions on the coefficients of such a series so that the
associated canonical tranformation has the desired property. In the second
part of [1] also takes into account the coefficient $J_{4}$ of the second
zonal harmonic and in the last part some sectorial harmonic are considered. It
has to be said that Brouwer credits Hugo von Zeipel in [6] for the above
mentioned idea of the power series expansion. Unfortunately we were not able
to find the original article so our main source is [1]: nevetheless we use the
term “von Zeipel method” to indicate the technique we want to discuss.

The porpose of the present work is to give a general and coincise presentation
of the von Zeipel method by using a more up-to-date vectorial notation. We use
a multi dimensional power series expansion, corresponding to a generic set of
parameters (not just $J_{2}$). In this way we can derive a more transparent
general formula, which allows to replicate the computations of [1] “all at
once”. Since such computation involve a lot of derivatives, we have used
Wolfram Mathematica to perform the actual work. In this way we could also
compute a second order term for the generating function that [1] omits. The
corresponding Mathematica notebook that we wrote for this is publicly
avaliable at the link:
https://www.wolframcloud.com/obj/claudiosaccon/VonZeipel2.nb

In order to make this note more self–contained, we added an appendix where we
recall the definitions and the main properties of canonical transformations
and generating functions.

We remind that, if $f:{}^{N}\to{}^{M}$, the _Jacobian matrix_ of $f$ is
defined by:

$J_{f}(x)=\frac{\partial f}{\partial x}(x):=\begin{pmatrix}\dfrac{\partial
f_{1}}{\partial x_{1}}(x)&\cdots&\dfrac{\partial f_{1}}{\partial x_{N}}(x)\\\
\vdots&\ddots&\vdots\\\ \dfrac{\partial f_{M}}{\partial
x_{1}}(x)&\cdots&\dfrac{\partial f_{M}}{\partial x_{N}}(x)\end{pmatrix}.$

If $M=1$, i.e. if $f$ is a scalar function, we sometimes consider the
_gradient_ of $f$:

$\nabla f(x):=J_{f}(x)^{t}=\begin{pmatrix}\dfrac{\partial f}{\partial
x_{1}}(x)\\\ \vdots\\\ \dfrac{\partial f}{\partial x_{N}}(x)\end{pmatrix}.$

When more than two indices are needed (for instance when we have to deal with
second order derivatives of vectors or matrices) we use Einstein’s convention
of repeated indices. In this note elements of N are regarded as “vectors”,
i.e. as $N$ by $1$ (column) matrices.

##  2 The satellite problem

The satellite problem can be written as:

$\ddot{\mathbf{r}}=-\frac{\mu\mathbf{r}}{|\mathbf{r}|^{3}}-\nabla
U(\mathbf{r})+\mathbf{F}(\mathbf{r})$ (1)

The term $\mu$ stands for $GM$ where $M$ is the mass of the Earth and $G$ is
the Gravitational Constant. In (1) $\nabla U$ is an additional conservative
term taking into account the non spherical shape of the Earth. In the case of
the “oblated Earth”

$U(\mathbf{r})=\sum_{n=2}^{\infty}\frac{\mu}{|\mathbf{r}|}J_{n}\left(\frac{R}{|\mathbf{r}|}\right)^{n}P_{n}\left(\sin(\beta)\right)$

In this formula $\beta$ denotes the latitude, $R$ is the radius oth the Earth,
and $P_{n}$ are the Legrendre polinomials. In this description the potential
is still rotationally symmetric (only “zonal harmonics are considered). The
even terms are also symmetric with respect to the equator, while the odd ones
introduce an asymmetry between the north and the south emisphere. In [1] terms
up to $J_{5}$ are considered (with different names). The methods presented in
this note rely on the fact that the coefficient $J_{n}$ and are small.

The term $\mathbf{F}$ allows to take into account external, nonconservative
forces acting on the satellite. The most important effect of this type is the
resistance of the athmosphere, usually called _drag_.

In this note we drop the term $\mathbf{F}$. In [2, 3, 4] the problem with drag
is considered (as a perturbation of the case without drag), but the technique
we are going to present here only works in the case $\mathbf{F}=0$. In the
Appendix, however, we show how the full problem gets transformed under a
canonical tranformation, so the reader can understand how this technique can
be relevant in the problem with drag.

##  3 The Von Zeipel method

Now we take a more general point of view. Let
$\mathcal{H}:{}^{N}\times{}^{N}\times{}^{M}\to\real$,
$\mathcal{H}=\mathcal{H}(p,q,\mathcal{J})=\mathcal{H}(v,\mathcal{J})$ be a
Hamiltonian depending on M parameters $J_{1},\dots,J_{M}$ which form the
vector $\mathcal{J}$ ($J_{i}$ will be supposed to be small). We remind that
for an $M\times M$ matrix the expression $A\mathcal{J}^{2}$ denotes
$\mathcal{J}^{t}A\mathcal{J}$. We suppose that

$\mathcal{H}(v,\mathcal{J})=\mathcal{H}_{0}(v)+\mathcal{H}_{1}(v)\cdot\mathcal{J}+\mathcal{H}_{2}(v)\mathcal{J}^{2}+0(|k^{2}|).$

where $\mathcal{H}_{0}:{}^{2N}\to\real$ (scalar),
$\mathcal{H}_{1}:{}^{2N}\to{}^{M}$ (vector valued), and
$\mathcal{H}_{2}:{}^{2N}\to{}^{M^{2}}$ (matrix valued). Of course (by Taylor’s
expansion:

$\mathcal{H}_{0}(v)=\mathcal{H}(v,0),\quad\mathcal{H}_{1}(v)=\nabla_{k}\mathcal{H}(v,0),\quad\mathcal{H}_{2}(v)=2\nabla_{k}^{2}\mathcal{H}(v,0)$

($\nabla_{k}$ and $\nabla^{2}_{k}$ denote the gradient and the Hessian matrix
with respect to the $k$ variable).

We assume that:

$\displaystyle\mathcal{H}_{0}\mbox{ does not depend on $q$};$ (2)
$\displaystyle\mathcal{H}_{1}\mbox{ and }\mathcal{H}_{2}\mbox{ are $2\pi$
periodic in each }q_{i}\quad i=1,\dots,N.$ (3)

The periodicity assumption means that the variables $q_{1},\dots,q_{N}$ are
_angles_. Given any function $f=f(p,q)$ we can define the _secular part_ of
$f$ as:

$f_{sec}(p):=\frac{1}{(2\pi)^{N}}\int_{0}^{2\pi}\cdots\int_{0}^{2\pi}f(p,q)\,dq_{1}\dots
dq_{N}$

(i.e the average in $q$ over the cube $[0,2\pi]\times\cdots\times[0,2\pi]$)
and the _periodic part_ of $f$:

$f_{per}(p,q):=f(p,q)-f_{sec}(p).$

It is clear that $f\mapsto f_{sec}$ and $f\mapsto f_{per}$ are linear, and
that $\left[f_{sec}\right]_{sec}=f_{sec}$,
$\left[f_{per}\right]_{per}=f_{per}$,
$\left[f_{sec}\right]_{per}=0=\left[f_{per}\right]_{sec}$. Moreover
$f=f_{sec}$ if and only if $f$ does not depend on $q$.

The purpose of this section is finding a canonical transformation such that
the new Hamiltonian $\mathcal{K}$ only depends on the (new) momenta $P$. To
this aim we look for a generating function of the form:

$S(P,q,\mathcal{J})=S_{0}(P,q)+S_{1}(P,q)\cdot\mathcal{J}+S_{2}(P,q)\mathcal{J}^{2}+o(|\mathcal{J}|^{2})$

(again with $S_{0}$ is scalar, $S_{1}$ an $M$ vector, and $S_{2}$ an $M\times
M$ matrix. For the moment we assume that such an $S$ is given and derive some
relationships between the old Hamiltonian $\mathcal{H}$ and the new one
$\mathcal{K}$. We then use such relationships to choose $S_{0}$, $S_{1}$, and
$S_{2}$ in such a way that $\mathcal{K}$ exhibis the desired property.

If $V(v)=\begin{pmatrix}P(p,q)\\\ Q(p,q)\end{pmatrix}$ be the new variables
induced by $S$, and $\mathcal{K}=\mathcal{K}(V)=\mathcal{K}(P,Q)$ be the new
Hamiltonian. We can write

$\mathcal{K}(V,k)=\mathcal{K}_{0}(V)+\mathcal{K}_{1}(V)\cdot\mathcal{J}+\mathcal{K}_{2}(V)\mathcal{J}^{2}+0(|\mathcal{J}^{2}|).$

where $\mathcal{K}_{0}$, $\mathcal{K}_{1}$, and $\mathcal{K}_{2}$ have the
same properties of $\mathcal{H}_{0}$, $\mathcal{H}_{1}$, and
$\mathcal{H}_{2}$. We know that:

$\mathcal{K}\left(P(p,q,\mathcal{J}),\frac{\partial
S(P(p,q),q,\mathcal{J})}{\partial P},k\right)=\mathcal{H}\left(\frac{\partial
S(P(p,q),q,\mathcal{J})}{\partial q},q,\mathcal{J}\right)\quad\forall
p,q,\mathcal{J}.$

In the following we are going to keep things coincise by skipping some of the
“inner” variables. By Taylor’s expansion:

$\mathcal{H}\left(\frac{\partial S}{\partial
q},q,\mathcal{J}\right)=\mathcal{H}_{0}\left(\frac{\partial S}{\partial
q},q\right)+\mathcal{H}_{1}\left(\frac{\partial S}{\partial
q},q\right)^{t}\mathcal{J}+\mathcal{H}_{2}\left(\frac{\partial S}{\partial
q},q\right)\mathcal{J}^{2}+0(|\mathcal{J}^{2}|)=\\\
\mathcal{H}_{0}+\frac{\partial\mathcal{H}_{0}}{\partial p_{i}}\frac{\partial
S_{1,k}}{\partial q_{i}}J_{k}+\frac{\partial\mathcal{H}_{0}}{\partial
p_{i}}\frac{\partial S_{2,k,h}}{\partial
q_{i}}J_{k}J_{h}+\frac{1}{2}\frac{\partial^{2}\mathcal{H}_{0}}{\partial
p_{i}\partial p_{j}}\frac{\partial S_{1,k}}{\partial q_{i}}\frac{\partial
S_{1,h}}{\partial q_{j}}J_{k}J_{h}+\\\
\mathcal{H}_{1,k}J_{k}+\frac{\partial\mathcal{H}_{1,h}}{\partial
p_{i}}\left(\frac{\partial S_{0}}{\partial q},q\right)\frac{\partial
S_{1,k}}{\partial
q_{i}}J_{k}J_{h}+\mathcal{H}_{2,k,h}J_{k}J_{h}+o(|\mathcal{J}|^{2}).$

($\mathcal{H}_{i}$ are computed in $\left(\frac{\partial S_{0}}{\partial
q},q\right)$). In the same way:

$\mathcal{K}\left(P,\frac{\partial S}{\partial
P},k\right)=\mathcal{K}_{0}+\frac{\partial\mathcal{K}_{0}}{\partial
Q_{i}}\frac{\partial S_{1,k}}{\partial
P_{i}}J_{k}+\frac{\partial\mathcal{K}_{0}}{\partial Q_{i}}\frac{\partial
S_{2,k,h}}{\partial P_{i}}J_{k}J_{h}+\\\
\frac{1}{2}\frac{\partial^{2}\mathcal{K}_{0}}{\partial Q_{i}\partial
Q_{j}}\frac{\partial S_{1,k}}{\partial P_{i}}\frac{\partial S_{1,h}}{\partial
P_{j}}J_{k}J_{h}+\mathcal{K}_{1,k}J_{k}+\frac{\partial\mathcal{K}_{1}}{\partial
Q_{i}}\frac{\partial S_{1,k}}{\partial
P_{i}}J_{k}J_{h}+\mathcal{K}_{2,k,h}J_{k}J_{h}+o(|\mathcal{J}|^{2})$

($\mathcal{K}_{i}$ are computed in $\left(P,\frac{\partial S_{0}}{\partial
P}\right)$). We take:

$S_{0}(P,q)=P\cdot q$

(if no other terms were present, then the associated canonical transformation
would be the identity) and we try to find $S_{1}$ and $S_{2}$ such that
$\mathcal{K}(V,k)=\mathcal{K}(P,k)$, that is $\mathcal{K}$ does not depend on
$Q$. By equating the terms having corresponding powers of $\mathcal{J}$, we
get:

$\displaystyle\mathcal{K}_{0}(P,Q)=\mathcal{K}_{0}(P)=\mathcal{H}_{0}(P)$ (4)
$\displaystyle\frac{\partial\mathcal{K}_{0}}{\partial Q_{i}}\frac{\partial
S_{1,k}}{\partial
P_{i}}+\mathcal{K}_{1,k}=\frac{\partial\mathcal{H}_{0}}{\partial
p_{i}}\frac{\partial S_{1,k}}{\partial q_{i}}+\mathcal{H}_{1,k}$ (5)
$\displaystyle\frac{\partial\mathcal{K}_{0}}{\partial Q_{i}}\frac{\partial
S_{2,k,h}}{\partial
P_{i}}+\frac{1}{2}\frac{\partial^{2}\mathcal{K}_{0}}{\partial Q_{i}\partial
Q_{j}}\frac{\partial S_{1,k}}{\partial P_{i}}\frac{\partial S_{1,h}}{\partial
P_{j}}+\frac{\partial S_{1,k}}{\partial
P_{i}}\frac{\partial\mathcal{K}_{1,h}}{\partial Q_{i}}+\mathcal{K}_{2,k,h}$
$\displaystyle=$ (6) $\displaystyle\frac{\partial\mathcal{H}_{0}}{\partial
p_{i}}\frac{\partial S_{2,k,h}}{\partial
q_{i}}+\frac{1}{2}\frac{\partial^{2}\mathcal{H}_{0}}{\partial p_{i}\partial
p_{j}}\frac{\partial S_{1,k}}{\partial q_{i}}\frac{\partial S_{1,h}}{\partial
q_{j}}+\frac{\partial S_{1,k}}{\partial
q_{i}}\frac{\partial\mathcal{H}_{1,k}}{\partial p_{i}}+\mathcal{H}_{2,k,h}$

(with $k$ and $h$ ranging between $1$ and $M$). Using (4) we have found
$\mathcal{K}_{0}$ (which only depends on $P$). Then all derivatives of
$\mathcal{K}_{0}$ in (5) and (6) are zero. If $\mathcal{K}_{1}$ and
$\mathcal{K}_{2}$ are to be undependent of $Q$, then both (5) and (6) split
into a pair of conditions:

$\displaystyle\mathcal{K}_{1,k}(P)=\left[\frac{\partial S_{1,k}}{\partial
q_{i}}\frac{\partial\mathcal{H}_{0}}{\partial
p_{i}}+\mathcal{H}_{1,k}\right]_{sec}$ (7) $\displaystyle
0=\left[\frac{\partial S_{1,k}}{\partial
q_{i}}\frac{\partial\mathcal{H}_{0}}{\partial
p_{i}}+\mathcal{H}_{1,k}\right]_{per}$ (8)
$\displaystyle\mathcal{K}_{2,k,h}(P)=\left[\frac{\partial\mathcal{H}_{0}}{\partial
p_{i}}\frac{\partial S_{2,k,h}}{\partial
q_{i}}+\frac{1}{2}\frac{\partial^{2}\mathcal{H}_{0}}{\partial p_{i}\partial
p_{j}}\frac{\partial S_{1,k}}{\partial q_{i}}\frac{\partial S_{1,h}}{\partial
q_{j}}+\frac{\partial S_{1,k}}{\partial
q_{i}}\frac{\partial\mathcal{H}_{1,h}}{\partial
p_{i}}+\mathcal{H}_{2,k,h}\right]_{sec}$ (9) $\displaystyle
0=\left[\frac{\partial\mathcal{H}_{0}}{\partial p_{i}}\frac{\partial
S_{2,k,h}}{\partial
q_{i}}+\frac{1}{2}\frac{\partial^{2}\mathcal{H}_{0}}{\partial p_{i}\partial
p_{j}}\frac{\partial S_{1,k}}{\partial q_{i}}\frac{\partial S_{1,h}}{\partial
q_{j}}+\frac{\partial S_{1,k}}{\partial
q_{i}}\frac{\partial\mathcal{H}_{1,h}}{\partial
p_{i}}+\mathcal{H}_{2,k,h}\right]_{per}$ (10)

A possible way to verify (7) and (8) is by imposing:

$\displaystyle\mathcal{K}_{1,k}(P):=\left[\mathcal{H}_{1,k}(P)\right]_{sec}$
(11) $\displaystyle\frac{\partial\mathcal{H}_{0}}{\partial
p_{i}}(P)\frac{\partial S_{1,k}}{\partial
q_{i}}(P,q)=-\left[\mathcal{H}_{1,k}(P,q)\right]_{per}$ (12)

In order to solve equation (12): let:

$\displaystyle w(P):=\frac{\partial\mathcal{H}_{0}}{\partial
p}(P),\quad,\quad\hat{w}(P):=\left\lVert\dfrac{\partial\mathcal{H}_{0}}{\partial
p}(P)\right\rVert^{-1}\dfrac{\partial\mathcal{H}_{0}}{\partial
p}(P)=\frac{w(P)}{\|w(P)\|}$

(notice that the “coefficient” $w$ neither depends on $q$ nor on $k$) and
define $S_{1}$ by:

$S_{1,k}(P,q):=-\frac{1}{\|w(P)\|}\int_{0}^{q\cdot\hat{w}}\left[\mathcal{H}_{1,k}\left(P,q-(1-t)(\hat{w}(P)\cdot
q)\hat{w}(P)\right)\right]_{per}\,dt=\\\
\frac{\left[\mathcal{H}_{1,k}\right]_{sec}(P)}{\left\lVert
w(P)\right\rVert}\hat{w}\cdot q-\frac{1}{\left\lVert
w(P)\right\rVert}\int_{0}^{q\cdot\hat{w}}\mathcal{H}_{1,k}\left(P,q-(1-t)(\hat{w}(P)\cdot
q)\hat{w}(P)\right)\,dt.$

It is then simple to check that $S_{1}$ solves Equation (12). So we can find
$\mathcal{K}_{1}$ and $S_{1}$. We can use the same idea to find
$\mathcal{K}_{2}$ and $S_{2}$:

$\displaystyle\mathcal{K}_{2,k,h}(P)=\left[\frac{1}{2}\frac{\partial^{2}\mathcal{H}_{0}}{\partial
p_{i}\partial p_{j}}\frac{\partial S_{1,k}}{\partial q_{i}}\frac{\partial
S_{1,h}}{\partial q_{j}}+\frac{\partial S_{1,k}}{\partial
q_{i}}\frac{\partial\mathcal{H}_{1,h}}{\partial
p_{i}}+\mathcal{H}_{2,k,h}\right]_{sec}$ (13)
$\displaystyle\frac{\partial\mathcal{H}_{0}}{\partial p_{i}}\frac{\partial
S_{2,k,h}}{\partial
q_{i}}=-\left[\frac{1}{2}\frac{\partial^{2}\mathcal{H}_{0}}{\partial
p_{i}\partial p_{j}}\frac{\partial S_{1,k}}{\partial q_{i}}\frac{\partial
S_{1,h}}{\partial q_{j}}+\frac{\partial S_{1,k}}{\partial
q_{i}}\frac{\partial\mathcal{H}_{1,h}}{\partial
p_{i}}+\mathcal{H}_{2,k,h}\right]_{per}$ (14)

Again (13) provides the expression of $\mathcal{K}_{2}$ (undependent of $Q$),
while (14) can be used to find $S_{2,k,h}$ (although in [1] $S_{2}$ is
neglected ????).

##  4 The solution of the satellite problem without drag

Let $r:=|\mathbf{r}|$, $a$ denote the osculating semi-major axis, $e$ denote
the eccentricity, and $\nu$ denote the true anomaly.

We consider the _Delaunay variables_ :

$\displaystyle(p_{1}=)\quad L$ $\displaystyle:=(\mu a)^{1/2},$
$\displaystyle(q_{1}=)\quad l$ $\displaystyle:=\mbox{ mean anomaly},$
$\displaystyle(p_{2}=)\quad G$ $\displaystyle:=L(1-e^{2})^{1/2},$
$\displaystyle(q_{2}=)\quad g$ $\displaystyle:=\mbox{ argument of the
pericenter},$ $\displaystyle(p_{3}=)\quad H$ $\displaystyle:=G\cos(I),$
$\displaystyle(q_{3}=)\quad h$ $\displaystyle:=\mbox{ longitude of ascending
node}.$

We first consider the problem with $J_{n}=0$ if $n\neq 2$ (only $J_{2}$ is
present). In terms of $(L,G,H,l,g,h)$ the Hamiltonian takes the form:

$\mathcal{H}(L,G,H,l,g,h):=\mathcal{H}_{0}(L)+\mathcal{H}_{1}(L,G,H,l,g)J_{2}$

where:

$\displaystyle\mathcal{H}_{0}(L)=$ $\displaystyle\frac{\mu^{2}}{2L^{2}}$ (15)
$\displaystyle\mathcal{H}_{1}(L,G,H,l,g)=$
$\displaystyle\frac{\mu^{4}R^{2}}{4L^{6}G^{2}}\frac{a^{3}}{\mathbf{r}^{3}}\left[\left(-G^{2}+3H^{2}\right)+\right.$
(16) $\displaystyle\qquad+\left.3\left(G^{2}-H^{2}\right)\cos(2g+2\nu)\right]$

Notice that (2) and (3) hold. To maintain the notations of [1] we use doubly
primes to denote the new variables: $P_{1}=L^{\prime\prime}$,
$P_{2}=G^{\prime\prime}$, $P_{2}=H^{\prime\prime}$, $Q_{1}=l^{\prime\prime}$,
$Q_{2}=g^{\prime\prime}$, $Q_{3}=h^{\prime\prime}$ (notice that we are
skipping the single prime variables since we are doing the two steps at once).
So
$\mathcal{K}=\mathcal{K}(L^{\prime\prime},G^{\prime\prime},H^{\prime\prime})$
and $S=S(L^{\prime\prime},G^{\prime\prime},H^{\prime\prime},l,g,h)$.

Of course the term $\dfrac{a^{3}}{\mathbf{r}^{3}}$ and the variable $\nu$ need
to be expressed in terms of the Delaunay variables. For this we can use the
Fourier series expansions:

$\displaystyle\frac{a^{3}}{\mathbf{r}^{3}}$
$\displaystyle=\frac{L^{3}}{G^{3}}+\sum_{n=1}^{\infty}2P_{n}(e)\cos(nl)$ (17)
$\displaystyle\frac{a^{3}}{\mathbf{r}^{3}}\cos(2g+2\nu)$
$\displaystyle=\sum_{n=-\infty}^{\infty}Q_{n}(e)\cos(2g+nl)$ (18)

where $P_{n}$ and $Q_{n}$ are suitable _Hansen coefficients_

Notice that $e=\left(1-\dfrac{G^{2}}{L^{2}}\right)^{1/2}$ so
$P_{n}=P_{n}(L,G)$, $Q_{n}=Q_{n}(L,G)$.

From (17) and (18) we infer that:

$\displaystyle\mathcal{H}_{1,sec}(L,G,H)=$
$\displaystyle\frac{\mu^{4}R^{2}}{4L^{3}G^{5}}\left(-G^{2}+3H^{2}\right)$
$\displaystyle\mathcal{H}_{1,per}(L,G,H)=$
$\displaystyle\frac{\mu^{4}R^{2}}{4L^{6}G^{2}}\left[\left(-G^{2}+3H^{2}\right)\left(\frac{a^{3}}{\mathbf{r}^{3}}-\frac{L^{3}}{G^{3}}\right)+\right.$
$\displaystyle\qquad+\left.\left(3G^{2}-3H^{2}\right)\frac{a^{3}}{\mathbf{r}^{3}}\cos(2g+2\nu)\right]$

Using the above formula we get from conditions (11)

$\mathcal{K}_{1}(L^{\prime\prime},G^{\prime\prime},H^{\prime\prime})=\frac{\mu^{4}R^{2}}{4{L^{\prime\prime}}^{3}{G^{\prime\prime}}^{5}}\left(-{G^{\prime\prime}}^{2}+3{H^{\prime\prime}}^{2}\right)$
(19)

while (12) turns into:

$\frac{\partial S_{1}}{\partial
l}(L^{\prime\prime},G^{\prime\prime},H^{\prime\prime},l,g)=\frac{\mu^{2}R^{2}}{4{L^{\prime\prime}}^{3}{G^{\prime\prime}}^{2}}\times\\\
\times\left[\left(-{G^{\prime\prime}}^{2}+3{H^{\prime\prime}}^{2}\right)\left(\frac{a^{3}}{\mathbf{r}^{3}}-\frac{{L^{\prime\prime}}^{3}}{{G^{\prime\prime}}^{3}}\right)+3\left({G^{\prime\prime}}^{2}-{H^{\prime\prime}}^{2}\right)\frac{a^{3}}{\mathbf{r}^{3}}\cos(2g+2\nu)\right].$
(20)

If we use the Fuorier expansions mentioned above, we get:

$S_{1}(L^{\prime\prime},G^{\prime\prime},H^{\prime\prime},l,g)=\frac{\mu^{2}R^{2}}{4{L^{\prime\prime}}^{3}{G^{\prime\prime}}^{2}}\left[\left(3{H^{\prime\prime}}^{2}-{G^{\prime\prime}}^{2}\right)\sum_{n=1}^{\infty}\frac{2P_{n}}{n}(L^{\prime\prime},G^{\prime\prime})\sin(nl)\right.+\\\
\left.3\left({G^{\prime\prime}}^{2}-{H^{\prime\prime}}^{2}\right)\sum_{n=-\infty}^{\infty}\frac{Q_{n}}{n}(L^{\prime\prime},G^{\prime\prime})\sin(2g+nl)\right]$

(notice that $S_{1}$ does not depend on $h$). However, using the formulas:

$dl=\frac{L}{G}\frac{r^{2}}{a^{2}}d\nu\quad,\quad\frac{a}{r}=\frac{1+e\cos(\nu)}{1-e^{2}}=\frac{L^{2}}{G^{2}}\left(1+\sqrt{1-\frac{G^{2}}{L^{2}}}\cos(\nu)\right)$
(21)

we can also find a closed formula for $S_{1}$. From (20) we get:

$S_{1}(L^{\prime\prime},G^{\prime\prime},h^{\prime\prime},l,g)=\frac{\mu^{2}R^{2}}{4{G^{\prime\prime}}^{5}}({G^{\prime\prime}}-3{H^{\prime\prime}}^{2})l+\\\
\frac{\mu^{2}R^{2}}{4{L^{\prime\prime}}^{3}{G^{\prime\prime}}^{2}}\int\left[\left(-{G^{\prime\prime}}^{2}+3{H^{\prime\prime}}^{2}\right)\left(\frac{L}{G}\right)\frac{1+e\cos(\nu)}{1-e^{2}}+\right.\\\
\left.3\left({G^{\prime\prime}}^{2}-{H^{\prime\prime}}^{2}\right)\left(\frac{L}{G}\right)\frac{1+e\cos(\nu)}{1-e^{2}}\cos(2g+2\nu)\right]d\nu$
(22)

which yields:

$S_{1}(L^{\prime\prime},G^{\prime\prime},H^{\prime\prime},l,g)=\frac{\mu^{2}R^{2}}{4{G^{\prime\prime}}^{5}}\left[({G^{\prime\prime}}^{2}-3{H^{\prime\prime}}^{2})(l-\nu-e\sin(\nu))\right.+\\\
\left.({G^{\prime\prime}}^{2}-{H^{\prime\prime}}^{2})\left(\frac{3}{2}\sin(2g+2\nu)+\frac{3e}{2}\sin(2g+\nu)+\frac{1}{2}\sin(2g+3\nu)\right)\right]$
(23)

In the same way from (13) we get:

$\mathcal{K}_{2}(L^{\prime\prime},G^{\prime\prime},H^{\prime\prime})=\left[\underbrace{\frac{1}{2}\frac{\partial^{2}\mathcal{H}_{0}}{\partial
L^{2}}\left(\frac{\partial S_{1}}{\partial
l}\right)^{2}+\frac{\partial\mathcal{H}_{1}}{\partial L}\frac{\partial
S_{1}}{\partial l}+\frac{\partial\mathcal{H}_{1}}{\partial G}\frac{\partial
S_{1}}{\partial
g}}_{=:\bar{\mathcal{H}}(L^{\prime\prime},G^{\prime\prime},H^{\prime\prime},l,g)}\right]_{sec}$
(24)

(every term on the R.H.S. is computed in
$(L^{\prime\prime},G^{\prime\prime},H^{\prime\prime},l,h,g)$)

and from (14):

$\frac{\partial S_{2}}{\partial
l}(L^{\prime\prime},G^{\prime\prime},H^{\prime\prime},l,g)=-\frac{2\mu^{2}R^{2}}{{L^{\prime\prime}}^{3}}\left[\frac{1}{2}\frac{\partial^{2}\mathcal{H}_{0}}{\partial
L^{2}}\left(\frac{\partial S_{1}}{\partial
l}\right)^{2}+\frac{\partial\mathcal{H}_{1}}{\partial L}\frac{S_{1}}{\partial
l}+\frac{\partial\mathcal{H}_{1}}{\partial G}\frac{S_{1}}{\partial
g}\right]_{per}$ (25)

If we compute the term $\bar{\mathcal{H}}$ on the right hand side of (13), we
get:

$\bar{\mathcal{H}}=\bar{\mathcal{H}}_{0}+\bar{\mathcal{H}}_{3}\left(\frac{a}{\mathbf{r}}\right)^{3}+\bar{\mathcal{H}}_{6}\left(\frac{a}{\mathbf{r}}\right)^{3}$

with:

$\displaystyle\bar{\mathcal{H}}_{0}(L,G,H)=$
$\displaystyle\frac{3\mu^{6}R^{4}\left(G^{2}-3H^{2}\right)^{2}}{32G^{10}L^{4}}$
$\displaystyle\bar{\mathcal{H}}_{3}\nu,g)=$
$\displaystyle-\frac{3\mu^{6}R^{4}}{16G^{8}L^{7}}\left(-12eG^{2}H^{2}L\sin^{2}(g+\nu)\cos(2g+\nu)\right.$
$\displaystyle\left.-4eG^{2}H^{2}L\sin^{2}(g+\nu)\cos(2g+3\nu)+\right.$
$\displaystyle\left.12eH^{4}L\sin^{2}(g+\nu)\cos(2g+\nu)+\right.$
$\displaystyle\left.4eH^{4}L\sin^{2}(g+\nu)\cos(2g+3\nu)+\right.$
$\displaystyle\left.3G^{5}\cos(2(g+\nu))-12G^{3}H^{2}\cos(2(g+\nu))\right.$
$\displaystyle\left.-12G^{2}H^{2}L\sin^{2}(g+\nu)\cos(2(g+\nu))+\right.$
$\displaystyle\left.9GH^{4}\cos(2(g+\nu))+12H^{4}L\sin^{2}(g+\nu)\cos(2(g+\nu))\right.$
$\displaystyle\left.-G^{5}+6G^{3}H^{2}-9GH^{4}\right)$
$\displaystyle\bar{\mathcal{H}}_{6}(L,G,H,\nu,g)=$
$\displaystyle-\frac{9\mu^{6}R^{4}}{32G^{4}L^{10}}\left(9G^{4}\cos^{2}(2(g+\nu))-6G^{4}\cos(2(g+\nu))\right.$
$\displaystyle\left.-18G^{2}H^{2}\cos^{2}(2(g+\nu))+24G^{2}H^{2}\cos(2(g+\nu))+\right.$
$\displaystyle\left.9H^{4}\cos^{2}(2(g+\nu))-18H^{4}\cos(2(g+\nu))+\right.$
$\displaystyle\left.G^{4}-6G^{2}H^{2}+9H^{4}\right)$

Using (21) we can take the average of $\bar{\mathcal{H}}$ with respect to $l$
and $g$ (by passing to the integral in $\nu$). We hence get:

$\mathcal{K}_{2}(L^{\prime\prime},G^{\prime\prime},H^{\prime\prime})=-\frac{3\mu^{6}R^{4}}{512{G^{\prime\prime}}^{13}{L^{\prime\prime}}^{5}}\left(99{G^{\prime\prime}}^{8}-48{G^{\prime\prime}}^{7}{L^{\prime\prime}}\right.\\\
\left.-2{G^{\prime\prime}}^{6}\left(167{H^{\prime\prime}}^{2}+495{L^{\prime\prime}}^{2}\right)+288{G^{\prime\prime}}^{5}{H^{\prime\prime}}^{2}{L^{\prime\prime}}+\right.\\\
\left.{G^{\prime\prime}}^{4}\left(307{H^{\prime\prime}}^{4}+2860{H^{\prime\prime}}^{2}{L^{\prime\prime}}^{2}+1155{L^{\prime\prime}}^{4}\right)\right.\\\
\left.-432{G^{\prime\prime}}^{3}{H^{\prime\prime}}^{4}{L^{\prime\prime}}-70{G^{\prime\prime}}^{2}\left(37{H^{\prime\prime}}^{4}{L^{\prime\prime}}^{2}+45{H^{\prime\prime}}^{2}{L^{\prime\prime}}^{4}\right)+2835{H^{\prime\prime}}^{4}{L^{\prime\prime}}^{4}\right)$

All the computations above have been checked in the Mathematica notebook
avaliable at https://www.wolframcloud.com/obj/claudiosaccon/VonZeipel2.nb .
The notations inside the notebook should be self–explanatory:

$\displaystyle\texttt{H0}=\mathcal{H}_{0},\quad\texttt{H1}=\mathcal{H}_{1},\quad\texttt{H1sec}=\mathcal{H}_{1,sec},\quad\texttt{K1}=\mathcal{K}_{1},\quad\texttt{K2}=\mathcal{K}_{2},$
$\displaystyle\texttt{S1}=S_{1},\quad\texttt{DS1l}=\frac{\partial
S_{1}}{\partial l},\quad\texttt{DS1g}=\frac{\partial S_{1}}{\partial
g},\quad\texttt{S2}=S_{2},$

while the terms RH* are sum up to form $\mathcal{H}_{2}$. The rows of the
notebook should be evaluated in sequence to provide the required expressions
of $S=S_{0}+S_{1}+S_{2}$ and the corresponding derivatives. As said before
also $S_{2}$ is computed according to our formulas. Since $S_{2}$ is made up
by a very long list of terms, we do not insert it here. Notice that [1]
considers $S_{2}$ as an upper order quantity and therefore neglectes it
(making no attempts to compute it).

##  Appendix A The Hamiltonian approach. Canonical transformations and
generating functions

Let $\mathcal{H}:{}^{N}\times{}^{N}\to\real$ be a smooth function (the
_Hamiltonian_), we we say that $p,q:\real\to{}^{N}$ are solutions to the
Hamiltonian equations, if:

$\left\\{\begin{aligned} &\dot{p}=-\frac{\partial}{\partial
q}\mathcal{H}(p,q)^{t},\\\ &\dot{q}=\frac{\partial}{\partial
p}\mathcal{H}(p,q)^{t}\end{aligned}\right.$ (26)

As well known (26) describe the dynamics of a system only effected by
conservative forces, according the laws of Classical Mechanics. The variables
$q$ are called “coordinates” while $p$ are the “momenta”. of the system.

In case we need to consider the presence of dissipative forces we add two two
terms to (26):

$\left\\{\begin{aligned} &\dot{p}=-\frac{\partial}{\partial
q}\mathcal{H}(p,q)^{t}+f_{1}(p,q),\\\ &\dot{q}=\frac{\partial}{\partial
p}\mathcal{H}(p,q)^{t}+f_{2}(p,q).\end{aligned}\right.$ (27)

where $f_{1},f_{2}:{}^{N}\times{}^{N}\to{}^{N}$. It will be sometimes
convenient to indicate $v=\begin{pmatrix}p\\\ q\end{pmatrix}$ and rewrite (27)
as:

$\dot{v}=J\nabla\mathcal{H}(v)+f,$ (28)

where:

$\nabla\mathcal{H}(p,q)=\begin{pmatrix}\dfrac{\partial}{\partial
p}\mathcal{H}(p,q)\\\ \ \\\ \dfrac{\partial}{\partial
q}\mathcal{H}(p,q)\end{pmatrix}\ ,\ J:=\begin{pmatrix}0&-I_{N}\\\ \\\
I_{N}&0\end{pmatrix}\ ,\ f(v)=f(p,q):=\begin{pmatrix}f_{1}(p,q)\\\ \ \\\
f_{2}(p,q)\end{pmatrix}$

($I_{N}$ denotes the $N$-dimensional identity matrix). Notice that
$J^{2}=-I_{N}$.

######  A.1 Definition.

Let $\Phi:{}^{2N}\to{}^{2N}$ be a diffeomorphism i.e. a bijection such that
$\Phi$ and $\Phi^{-1}$ are smooth. We say that $\Phi$ is a _canonical
transformation_ , if for any solution $v$ of $\dot{v}=J\nabla\mathcal{H}(v)$,
the transformed function $V:=\Phi(v)$ is a solution of
$\dot{V}=J\nabla\mathcal{K}(V),$ where
$\mathcal{K}(V):=\mathcal{H}(\Phi^{-1}(V))$. This is usually expressed by
saying that “$\Phi$ preserves the form of Hamilton’s equations”.

######  A.2 Theorem.

Let $\Phi$ be a diffeomorphism. If $\dfrac{\partial\Phi}{\partial v}$ is
symplectic, that is if:

$\frac{\partial\Phi}{\partial v}J\frac{\partial\Phi}{\partial v}^{t}=J$ (29)

holds, then $\Phi$ is a canonical transformation.

###### Proof.

From the definition of $\mathcal{K}$ we have:
$\mathcal{K}(\Phi(v))=\mathcal{H}(v)$. Then:

$\frac{\partial\mathcal{H}(v)}{\partial
v}=\frac{\partial\mathcal{K}(V)}{\partial V}\frac{\partial\Phi(v)}{\partial
v}.$

Now, if $\dot{v}=J\nabla\mathcal{H}(v)$ and $V=\Phi(v)$, we have:

$\dot{V}=\frac{\partial\Phi(v)}{\partial
v}\dot{v}=\frac{\partial\Phi(v)}{\partial
v}J\frac{\partial\mathcal{H}(v)}{\partial
v}^{t}=\frac{\partial\Phi(v)}{\partial v}J\frac{\partial\Phi(v)}{\partial
v}^{t}\frac{\partial\mathcal{K}(V)}{\partial V}^{t}.$ (30)

We can conclude that, if (29) holds, then $\dot{V}=J\nabla\mathcal{K}(V).$ ∎

######  A.3 Remark.

As shown in 1 of (B.2) the equality (29) is equivalent to:

$\frac{\partial P}{\partial p}\frac{\partial P}{\partial q}^{t}=\frac{\partial
P}{\partial q}\frac{\partial P}{\partial p}^{t}\qquad\frac{\partial
Q}{\partial p}\frac{\partial Q}{\partial q}^{t}=\frac{\partial Q}{\partial
q}\frac{\partial Q}{\partial p}^{t}\qquad\frac{\partial P}{\partial
p}\frac{\partial Q}{\partial q}^{t}-\frac{\partial P}{\partial
q}\frac{\partial Q}{\partial p}^{t}=I_{N}$ (31)

or equivalently to:

$\frac{\partial P}{\partial p}^{t}\frac{\partial Q}{\partial p}=\frac{\partial
Q}{\partial p}^{t}\frac{\partial P}{\partial p}\qquad\frac{\partial
P}{\partial q}^{t}\frac{\partial Q}{\partial q}=\frac{\partial Q}{\partial
q}^{t}\frac{\partial P}{\partial q}\qquad\frac{\partial P}{\partial
p}^{t}\frac{\partial Q}{\partial q}-\frac{\partial Q}{\partial
p}^{t}\frac{\partial P}{\partial q}=I_{N}$ (32)

######  A.4 Remark.

If $\dfrac{\partial\Phi}{\partial v}$ is symplectic, then the change of
variables $V=\Phi(v)$ transforms the nonconservative equation (27) into:

$\dot{V}=J\nabla\mathcal{K}(V)+F(V)$

where $\mathcal{K}=\mathcal{H}\circ\Phi^{-1}$ and:

$F(V)=\frac{\partial\Phi}{\partial
v}(\Phi^{-1}(V))f(\Phi^{-1}(V))=\left(\frac{\partial\Phi^{-1}}{\partial
V}(V)\right)^{-1}f(\Phi^{-1}(V)).$

But since $\dfrac{\partial\Phi}{\partial v}$ is symplectic, so is
$\dfrac{\partial\Phi^{-1}}{\partial V}$ (see section B and we have:

$\frac{\partial\Phi^{-1}(V)}{\partial V}=\begin{pmatrix}\dfrac{\partial
p(P,Q)}{\partial P}&\dfrac{\partial p(P,Q)}{\partial Q}\\\ \dfrac{\partial
q(P,Q)}{\partial P}&\dfrac{\partial q(P,Q)}{\partial
Q}\end{pmatrix}^{-1}=\begin{pmatrix}\dfrac{\partial q(P,Q)}{\partial
Q}^{t}&-\dfrac{\partial p(P,Q)}{\partial Q}^{t}\\\ -\dfrac{\partial
q(P,Q)}{\partial P}^{t}&\dfrac{\partial p(P,Q)}{\partial P}^{t}.\end{pmatrix}$

If $\Phi^{-1}(P,Q)=\begin{pmatrix}p(P,Q)\\\ q(P,Q)\end{pmatrix}$, we have
finally:

$\displaystyle F_{1}(P,Q)=$ $\displaystyle\dfrac{\partial q(P,Q)}{\partial
Q}^{t}f_{1}(\Phi^{-1}(P,Q))-\dfrac{\partial p(P,Q)}{\partial
Q}^{t}f_{2}(\Phi^{-1}(P,Q))$ (33) $\displaystyle F_{2}(P,Q)=$
$\displaystyle-\dfrac{\partial q(P,Q)}{\partial
P}^{t}f_{1}(\Phi^{-1}(P,Q))+\dfrac{\partial p(P,Q)}{\partial
P}^{t}f_{2}(\Phi^{-1}(P,Q))$ (34)

A possible way of obtaining a canonical transformation is via a generating
function.

######  A.5 Definition (generating functions).

Let $S:{}^{N}\times{}^{N}\to\real$ be a smooth function. We see $S$ as
$S(P,q)$ (a combination of “new variables” $P$ and “old ones” $q$). It is
simple to check that the equations:

$p=\frac{\partial S}{\partial q}(P,q)\quad,\quad Q=\frac{\partial S}{\partial
P}(P,q)$ (35)

implicitely define a canonical transformation $(P,Q)=\Phi(p,q)$, in the sense
that $\Phi:{}^{N}\times{}^{N}\to{}^{N}\times{}^{N}$ is canonical and:

$(P,Q)=\Phi(p,q)\Leftrightarrow\mbox{ \eqref{eqn:definition-canonical} holds
}\qquad\qquad\forall\ p,q,P,Q.$

Such an $S$ will be called a _generating function_ for the canonical
transformation $\Phi$.

######  A.6 Theorem.

Let $\Phi(p,q)=\begin{pmatrix}P(p,q)\\\ Q(p,q)\end{pmatrix}$ be a
diffeomorphism and $\dfrac{\partial\Phi}{\partial v}$ be symplectic.

Let $(p_{0},q_{0})$ be a point such that
$\det\left(\dfrac{\partial\Phi_{1}}{\partial p}(p_{0},q_{0})\right)\neq 0$ and
let $(P_{0},Q_{0})=\Phi(p_{0},q_{0})$. Then there exist a neighboorhood $W$ of
$(P_{0},q_{0})$, a neighboorhood $W_{1}$ of $(p_{0},Q_{0})$, and a smooth
function $S:W\to\real$ such that $S$ is a generating function for $\Phi$ in
$W$.

###### Proof.

Let $G:{}^{N}\times{}^{N}\times{}^{N}\times{}^{N}$ be defined by:

$G(P,Q,p,q):=\begin{pmatrix}P\\\ Q\end{pmatrix}-\Phi(p,q)=\begin{pmatrix}P\\\
Q\end{pmatrix}-\begin{pmatrix}\Phi_{1}(p,q)\\\ \Phi_{2}(p,q)\end{pmatrix}.$

We have:

$\frac{\partial
G}{\partial(P,Q,p,q)}=\begin{pmatrix}I_{N}&0&-\dfrac{\partial\Phi_{1}}{\partial
p}&-\dfrac{\partial\Phi_{1}}{\partial q}\\\ 0&I_{N}\
&-\dfrac{\partial\Phi_{2}}{\partial p}&-\dfrac{\partial\Phi_{2}}{\partial
q}\end{pmatrix}.$

From the assumption we derive that:

$\det\left(\dfrac{\partial
G}{\partial(Q,p)}\right)=\det\left(\begin{pmatrix}0&-\dfrac{\partial\Phi_{1}}{\partial
p}\\\ I_{N}&-\dfrac{\partial\Phi_{2}}{\partial
p}\end{pmatrix}\right)=\det\left(\frac{\partial\Phi_{1}}{\partial
p}\right)\neq 0$

By the implicit function Theorem there exist a neighboorhood $W$ of
$(P_{0},q_{0})$, a neighboorhood $W_{1}$ of $(Q_{0},p_{0}$ and a map
$\Psi:W\to{}^{N}\times{}^{N}$ such that:

$(P,Q)=\Phi(p,q)\quad\Leftrightarrow\quad(Q,p)=\Psi(P,q).$

Moreover:

$\frac{\partial\Psi}{\partial(P,q)}=-\begin{pmatrix}0&-\dfrac{\partial\Phi_{1}}{\partial
p}\\\ I_{N}&-\dfrac{\partial\Phi_{2}}{\partial
p}\end{pmatrix}^{-1}\begin{pmatrix}I_{N}&-\dfrac{\partial\Phi_{1}}{\partial
q}\\\ 0&-\dfrac{\partial\Phi_{2}}{\partial q}\end{pmatrix}=\\\
\begin{pmatrix}\dfrac{\partial\Phi_{2}}{\partial
p}\left(\dfrac{\partial\Phi_{1}}{\partial p}\right)^{-1}&-I_{N}\\\
\left(\dfrac{\partial\Phi_{1}}{\partial
p}\right)^{-1}&0\end{pmatrix}\begin{pmatrix}I_{N}&-\dfrac{\partial\Phi_{1}}{\partial
q}\\\ 0&-\dfrac{\partial\Phi_{2}}{\partial q}\end{pmatrix}=\\\
\begin{pmatrix}\dfrac{\partial\Phi_{2}}{\partial
p}\left(\dfrac{\partial\Phi_{1}}{\partial
p}\right)^{-1}&-\dfrac{\partial\Phi_{2}}{\partial
p}\left(\dfrac{\partial\Phi_{1}}{\partial
p}\right)^{-1}\dfrac{\partial\Phi_{1}}{\partial
q}+\dfrac{\partial\Phi_{2}}{\partial q}\\\
\left(\dfrac{\partial\Phi_{1}}{\partial
p}\right)^{-1}&\left(\dfrac{\partial\Phi_{1}}{\partial
p}\right)^{-1}\dfrac{\partial\Phi_{1}}{\partial q}\end{pmatrix}.$

We claim that the above matrix is symmetric. The upper left block is symmetric
since:

$\dfrac{\partial\Phi_{2}}{\partial p}\left(\dfrac{\partial\Phi_{1}}{\partial
p}\right)^{-1}=\left(\dfrac{\partial\Phi_{1}}{\partial
p}\right)^{-t}\left(\dfrac{\partial\Phi_{2}}{\partial
p}\right)^{t}\Leftrightarrow\left(\dfrac{\partial\Phi_{1}}{\partial
p}\right)^{t}\dfrac{\partial\Phi_{2}}{\partial
q}=\left(\dfrac{\partial\Phi_{2}}{\partial
p}\right)^{t}\dfrac{\partial\Phi_{1}}{\partial p}$

and the latter is exactly the first equality in (32). In the same way the
lower right block is symmetric:

$\left(\dfrac{\partial\Phi_{1}}{\partial
p}\right)^{-1}\dfrac{\partial\Phi_{1}}{\partial
q}=\left(\dfrac{\partial\Phi_{1}}{\partial
q}\right)^{t}\left(\dfrac{\partial\Phi_{1}}{\partial
p}\right)^{-t}\Leftrightarrow\dfrac{\partial\Phi_{1}}{\partial
q}\left(\dfrac{\partial\Phi_{1}}{\partial
p}\right)^{t}=\dfrac{\partial\Phi_{1}}{\partial
p}\left(\dfrac{\partial\Phi_{1}}{\partial q}\right)^{t}$

which is the first equality in (31). Finally:

$\left(\dfrac{\partial\Phi_{1}}{\partial
p}\right)^{-t}=-\dfrac{\partial\Phi_{2}}{\partial
p}\left(\dfrac{\partial\Phi_{1}}{\partial
p}\right)^{-1}\dfrac{\partial\Phi_{1}}{\partial
q}+\dfrac{\partial\Phi_{2}}{\partial q}\Leftrightarrow\\\
I_{N}=-\dfrac{\partial\Phi_{2}}{\partial
p}\left(\dfrac{\partial\Phi_{1}}{\partial
p}\right)^{-1}\underbrace{\dfrac{\partial\Phi_{1}}{\partial
q}\left(\dfrac{\partial\Phi_{1}}{\partial
p}\right)^{t}}_{=\frac{\partial\Phi_{1}}{\partial
p}\left(\frac{\partial\Phi_{1}}{\partial
q}\right)^{t}}+\dfrac{\partial\Phi_{2}}{\partial
q}\left(\dfrac{\partial\Phi_{1}}{\partial p}\right)^{t}\Leftrightarrow\\\
I_{N}=-\dfrac{\partial\Phi_{2}}{\partial
p}\left(\dfrac{\partial\Phi_{1}}{\partial
q}\right)^{t}+\dfrac{\partial\Phi_{2}}{\partial
q}\left(\dfrac{\partial\Phi_{1}}{\partial p}\right)^{t}\Leftrightarrow\\\
I_{N}=(I_{N}^{t}=)-\dfrac{\partial\Phi_{1}}{\partial
q}\left(\dfrac{\partial\Phi_{2}}{\partial
p}\right)^{t}+\dfrac{\partial\Phi_{1}}{\partial
p}\left(\dfrac{\partial\Phi_{2}}{\partial q}\right)^{t}$

and the last equality coincides with the third one in (31).

Since $\dfrac{\partial\Psi}{\partial(P,q)}$ is symmetric, there exists a
“potential” $S:\to\real$ such that:

$\Psi(P,q)=\nabla S(P,q)\Leftrightarrow\Psi_{1}(P,q)=\frac{\partial
S}{\partial P}(P,q)\ ,\ \Psi_{2}(P,q)=\frac{\partial S}{\partial q}(P,q)$

for all $(P,q)$ in $W$. This concludes the proof. ∎

##  Appendix B Simplectic matrices

######  B.1 Definition.

A $2N\times 2N$ matrix $M$ is said to be _symplectic_ if $MJM^{t}=M$.

######  B.2 Proposition (Properties of symplectic matrices).

Let $M,M_{1}$ be $2N\times 2N$ matrices and suppose that
$\displaystyle{M=\begin{pmatrix}A&B\\\ C&D\end{pmatrix}}$, where $A,B,C,D$ are
four $N\times N$ matrices.

  1. 1.

$M$ is symplectic if and only if:

$AB^{t}=BA^{t}\quad,\quad CD^{t}=DC^{t}\quad,\quad AD^{t}-BC^{t}=I_{N}$

if and only if:

$A^{t}C=C^{t}A\quad,\quad B^{t}D=D^{t}B\quad,\quad A^{t}D-C^{t}B=I_{N}$

  2. 2.

If $M$ and $M_{1}$ are symplectic, then $MM_{1}$ is symplectic.

  3. 3.

If $M$ is symplectic, than $M$ is invertible and

$M^{-1}=-JM^{t}J=\begin{pmatrix}D^{t}&-B^{t}\\\ -C^{t}&A^{t}\end{pmatrix}.$

  4. 4.

If $M$ is symplectic, than $M^{t}$ is symplectic.

###### Proof.

  1. 1.

We have:

$MJM^{t}=\begin{pmatrix}A&B\\\ C&D\end{pmatrix}\begin{pmatrix}0&-I_{N}\\\
I_{N}&0\end{pmatrix}\begin{pmatrix}A^{t}&C^{t}\\\ B^{t}&D^{t}\end{pmatrix}=\\\
\begin{pmatrix}B&-A\\\ D&-C\end{pmatrix}\begin{pmatrix}A^{t}&C^{t}\\\
B^{t}&D^{t}\end{pmatrix}=\begin{pmatrix}BA^{t}-AB^{t}&BC^{t}-AD^{t}\\\
DA^{t}-CB^{t}-&DC^{t}-CD^{t}.\end{pmatrix}$

So the condition $MJM^{t}=J$ turns out to be equivalent to the first set of
equalities above (the possible fourth equality is easily derived from the
third one, by taking the transpose). The second set of equalities follows from
the fact that $M^{t}$ is symplectic, as shown in 4.

  2. 2.

We have:

$MM_{1}J(MM_{1})^{t}=M\underbrace{M_{1}JM_{1}^{t}}_{=J}M^{t}=MJM^{t}=J.$

  3. 3.

We have:

$M(-JM^{t}J)=-MJM^{t}J=-JJ=I_{N}.$

This implies that $-JM^{t}J$ is a left inverse. A standard argument shows that
$M$ is invertible and $M^{-1}=-JM^{t}J$.

  4. 4.

From the previous statement we have $-JM^{t}JM=I_{N}$. By multiplying by $J$,
we get that $M^{t}$ is symplectic.

∎

## References

  * [1] Dirk Brouwer.  Solution of the problem of artificial satellite theory without drag.  Astron. J., 64:378–96, November 1959.
  * [2] Dirk Brouwer and Gen-ichiro Hori.  Theoretical evaluation of atmospheric drag effects in the motion of an artificial satellite.  Astron. J., 66:193–225, 1961.
  * [3] M. Lane.  The development of an artificial satellite theory using a power-law atmospheric density representation.  pages 65–35, 01 1965.
  * [4] M. Lane and K. H. Crawford.  An improved analytical drag theory for the artificial satellite problem.  pages 69–925, 1969.
  * [5] Hoots Felix R. and Roehrich Ronald L.  Models for propagation of norad element sets.  https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&ved=2ahUKEwibku653Zz8AhUtQ_EDHUWcBSwQFnoECA0QAQ&url=https%3A%2F%2Fcelestrak.org%2FNORAD%2Fdocumentation%2Fspacetrk.pdf&usg=AOvVaw2liZ-7R20gznLO1S-74IZ2, 1980\.
  * [6] H. von Zeipel.  La loi des luminosites dans l’amas globulaire M 3. , v. 11, p. 22. .  Ark. Mat. Astr. Phys., 11:22, 1916.


* Zhang et al. (2019c) Cheng Zhang, Chao Fan, Wenlin Yao, Xia Hu, and Ali Mostafavi. 2019c.  Social media for intelligent public information and warning in disasters: An interdisciplinary review.  _International Journal of Information Management_ 49 (2019), 190–207.  https://doi.org/10.1016/j.ijinfomgt.2019.04.004
  * Zhang et al. (2019a) Yixuan Zhang, Kartik Chanana, and Cody Dunne. 2019a.  IDMVis: Temporal Event Sequence Visualization for Type 1 Diabetes Treatment Decision Support.  _IEEE Transactions on Visualization and Computer Graphics_ 25, 1 (2019), 512–522.  https://doi.org/10.1109/TVCG.2018.2865076
  * Zhang et al. (2019b) Yixuan Zhang, Sara Di Bartolomeo, Fangfang Sheng, Holly Jimison, and Cody Dunne. 2019b.  Evaluating Alignment Approaches in Superimposed Time-Series and Temporal Event-Sequence Visualizations. In _2019 IEEE Visualization Conference (VIS)_. IEEE, 501 Hoes Ln #3, Piscataway, NJ 08855, 1–5.  https://doi.org/10.1109/VISUAL.2019.8933584
  * Zhang et al. (2020) Yixuan Zhang, Nurul Suhaimi, Rana Azghandi, Mary Amulya Joseph, Miso Kim, Jacqueline Griffin, and Andrea G. Parker. 2020.  Understanding the Use of Crisis Informatics Technology among Older Adults. In _Proceedings of the 2020 CHI Conference on Human Factors in Computing Systems_ (Honolulu, HI, USA) _(CHI ’20)_. Association for Computing Machinery, New York, NY, USA, 1–13.  https://doi.org/10.1145/3313831.3376862