11institutetext: University of Cincinnati, School of Information Technology,
Cincinnati, Ohio, USA
11email<EMAIL_ADDRESS>
# Deep Learning Algorithm for Threat Detection in Hackers Forum (Deep Web)

Victor Adewopo* 11 Bilal Gonen 11 Nelly Elsayed 11 Murat Ozer 11 Zaghloul Saad
Elsayed 11

###### Abstract

In our current society, the inter-connectivity of devices provides easy access
for netizens to utilize cyberspace technology for illegal activities. The deep
web platform is a consummative ecosystem shielded by boundaries of trust,
information sharing, trade-off, and review systems. Domain knowledge is shared
among experts in hackers forums which contain indicators of compromise that
can be explored for cyberthreat intelligence. Developing tools that can be
deployed for threat detection is integral in securing digital communication in
cyberspace. In this paper, we addressed the use of TOR relay nodes for
anonymizing communications in deep web forums. We propose a novel approach for
detecting cyberthreats using a deep learning algorithm Long Short-Term Memory
(LSTM). The developed model outperformed the experimental results of other
researchers in this problem domain with an accuracy of 94% and precision of
90%. Our model can be easily deployed by organizations in securing digital
communications and detection of vulnerability exposure before cyberattack.

###### Keywords:

Deep Learning Algorithm Deepweb Cyberthreat Cyberattack Threat prediction

##  1 Introduction

In our contemporary society, the existence of interconnected devices in
cyberspace is faced with high vulnerabilities. However, the ubiquitous use of
cyber-technology, the internet of things, and smart technologies steered the
revolution and advancement in the quality of human lives. Cyberspace is one of
the most complex systems ever built by human beings, the resources are widely
used by many, but only a few understand the complexities of cyberspace. The
availability of sophisticated technologies enables criminals to exploit new
ways of committing crimes. Criminogenic activities in cyberspace are a
critical societal problem that requires the swift intervention of
cybersecurity experts [3, 30]. There are often misconceptions with the terms
”Surface web,” ”Deep web,” and ”Dark web,” they are relatively interconnected
but do not mean exactly the same thing. The surface web is web pages
unencrypted and can be accessed using traditional search engines (e.g.,
Google, Bing, Yahoo). The surface web consists of billions of static web
pages, and it occupies only 10% of the internet space.

The dark web is a layer within the deep web and is not accessible using a
standard browser. They are intentionally hidden parts of the web but can be
accessed using a specific URL address. The deep web contains about 90% of the
contents available on the internet [43]. Deep web contents are available on
the web, but cannot be indexed or accessed by using regular search engines.
Several techniques and tools have been designed to understand and crawl the
deep web. A recent report indicated the low harvest rate of the deep web-
about 647,000 distinct web forms was found by sampling 25 million pages using
Google index [43, 12]. Darknet is an encrypted network technology that uses
the internet infrastructure and can only be accessed using special network
configuration and software tools to access its contents [23].

Figure 1: Layers of Internet: Surface Web, Deep Web, And Dark Web Explained

Traditionally, deep web hacker forums are used for criminal activities,
including different types of cybercrime, procurement of illegal drugs, child
pornography, cryptocurrency markets, and hacking [35]. Cryptocurrency (i.e.,
Bitcoins) is adopted as a method of payment in hacker forums because of its’
anonymous characteristics, which make it more difficult to link the users with
any specific bank account details [31]. In the quest for high-quality
security-related information, developing methodologies to crawl the deep web
and excavate wealth of information for cyberthreat intelligence is a paramount
concern [2]. A Cybersecurity firm ”UpGuard” reported 419 million Facebook
users’ data breach in April 2019. Data exposed publicly on the Amazon server
includes; passwords, user IDs, and check-ins [32]. Similarly, in September
2018, over 50 million Facebook users’ breached data was auctioned at a bitcoin
value of 3$ per each user’s data on the Darknet [8]. Nefarious actors have
been able to build a digital ecosystem that provides platforms famously
acclaimed to promote activities such as cybercrime, terrorism, hacking,
procurements of illegal drugs, arms deals, and anonymous markets [18]. The
hidden nature of Darknet websites aided fraudsters in utilizing cyberspace for
criminal activities.

##  2 User’s Motivation In Deep web forums

Why do people want to use the deep web? Over the years, the deep web has been
traditionally used by actors to carry out nefarious activities. In countries
with high censorship of all digital communications, including indiscriminate
censorship of media by oppressive regimes have led to an increase in the use
of hidden services. For example, in March 2019, a darknet website was used to
publish the leaked emails and official documents of the Russian Government to
expose corruption and promote transparency [25]. The darknet is classified as
a platform for private networks and infrastructure for the public systems, in
that the infrastructure is considered as a public entity. When the services
are used as a private platform, users communicate anonymously on the network,
and they consent to store data on their local computers. Most of these stored
data are inaccessible to the local users who harbor the stored data in their
devices. Platforms acquire the existing infrastructure of public networks by
enhancing security features and constraints to privatize communications [33].
Gehl et al. [17] suggest that the platform tends to be of a more private
network while infrastructures lean more toward public entities. The anonymous
nature of the darknet could provide individuals with both well-intended
motives and malicious intentions to carry out various activities on the
platform [17]. Wang et al [41] analyzed the personal characteristics of the
darknet users to identify typical individual user’s profile. They aimed at
developing a tool to track darknet users and reduce cybercrime. A three-step
approach was used in analyzing the characteristics of darknet website users
through Named Entity Recognition and Information Extraction. A. Block
Filtration. B. Attribute candidate generation. C. Attribute candidate
verification.

###  2.1 Use of TOR for Anonymizing Services in Deep Web Forums

TOR is an open-source software that secures anonymous communications and also
allows users to host hidden services on the TOR network [38]. TOR services can
dynamically allow users to access the contents of websites without revealing
their identity, which prevents government censors from monitoring sites
accessed by users. Accordingly, bloggers could host websites and hidden
services by publishing contents on Distributed Hash Tables (DHT) and TOR
services which makes the directory of the database random in which all
participating nodes can access and retrieve information [36]. This makes it
difficult to censor, control, or access the location of the website’s contents
[29, 20]. The downside of such services is that the ubiquitous utilization of
hidden services could pose a significant threat compared to the benefits
mentioned above. Hackers also use TOR services, and the lack of national laws
and international cooperation among relevant agencies are two significant
challenges in combating the investigation of criminal activities committed
using Tor hidden services [15]. TOR software ensures user’s anonymity by
relying on an overlay of networks using three relay nodes: (1) Entry relay
nodes, (2) Middle relay nodes, and (3) Exit relay nodes.

  * •

Entry Relay Nodes
The entry relay nodes serve as the entry point to TOR services, also called
entry guards, which receive the TOR traffic and pass it onto other nodes that
anonymize users as the originator of traffic. This node encrypts communication
between the traffic source and the entry relay nodes. The IP addresses of the
users connected to the TOR network are logged in the list of guard nodes to
keep track and distinguish entry relay nodes [7, 14]

  * •

Middle Relay Nodes
The middle relay nodes transport traffic from an entry relay to the exit relay
to ensure anonymity and bridge gaps between the entry and exit relay. Each
intermediate node can only communicate with its predecessor and successor
relay in Tor services. The nodes are configured in such a way that
communications cannot be linked (send/receive data) without the intermediate
nodes [3].

  * •

Exit Relay Nodes
The exit relay is the final relay that Tor traffic passes through before it
reaches its destination. The IP address of the exit relay is interpreted as
the source of TOR traffic. The exit nodes contain the key of encrypted
information in the Tor network. The original information from the source
traffic can only be recovered through the exit relay.

The exit relay servers are often the sources for complaints and associated
with illegal activities perpetrated by source traffic, which attracts law
enforcement agencies’ attention. Sun et al. [37] discovered two ways of
observing the TOR network and identifying anonymous users, which could be
achieved by either owning multiple TOR relays or manipulating the underlying
network communications. TOR is the most common anonymizing service that runs
over seven thousand relays. TOR uses bridges, which are new TOR routers that
are not listed on the Tor network. Some Internet Service Providers (ISPs)
blocked TOR nodes and access to the TOR network by filtering the IP address of
TOR servers. Bridges are used in circumventing blocked IP addresses of listed
Tor relays [7, 38]

##  3 Pattern of Cyberthreat between 2019 and 2021

A cyberattack can come in the form of malware that is virulent in nature.
Oftentimes, they may look benign on sight but are usually embedded with
malicious infections. Most malware infections originate in a form of
attachment with the ability to hide, replicate themselves on a target system,
and spread to other computers on a networked device, while some malicious
infections such as the trojan virus do not replicate themselves. Malware is
somewhat platform-centric, targeting specific security loopholes to explore
areas with high vulnerability. In a malicious ransomware attack, victims’ file
disks are usually partitioned, changed, or encrypted with a key that can only
be decrypted by a mathematically related key [39]. The existence of
sophisticated technology can be deployed in predicting cyberattacks before the
actual occurrence. The Wall Street Journal reported the likelihood of a
federal judiciary’s systems breached in a SolarWinds hack which prompted the
federal court to store sensitive documents filled in a stand-alone computer
system [40]. In March 2021, Over 30,000 organizations across the U.S were
hacked by cyber espionage focused on siphoning email communications through
the Microsoft exchange flaws. In 2019, IBM security reported that the average
cost of a data breach in the U.S. is about $3.9 million, the average life
cycle of a data breach stands at 314 days, and the impact on an organization
can span over 4 years with 67% of the cost occurring within the first year
[34]. In 2020, over 273% increase in the number of breached records in the
first quarter of 2020 was observed as compared to the first quarter of 2019
[27]. In July 2020, attackers swindle $121,000 worth of bitcoin in 3,000
transactions, by hacking Twitter accounts of 130 high profile US personalities
and resetting 45 user accounts passwords. The phone Spear-phishing method of
attack was targeted towards a few employees to gain access to internal
systems. In April 2020, approximately 500,000 stolen zoom passwords were
available for sale in darknet markets while some account credentials were made
available free. Information compromised includes login credentials, victim
personal meeting URLs, and host keys [27]. Nefarious actors leverage the
virtual ecosystem of cyberspace to perpetrate sophisticated cyberattacks
launched against an unsuspecting individual. In April 2019, there were two
incident reports of Facebook data breach. A database with about 419 million
unique Facebook users’ record, containing ID, phone numbers, and other
personal details linked to users account was exposed on an Amazon cloud
service. Another data breach contains Users’ ID, Facebook friends, likes,
photos, Facebook check-ins, and some other sensitive information that was
leaked in the Amazon S3 bucket cloud server [1]. In June 2019, US medical bill
debt collector, American Medical Collection Agency (AMCA) filed for bankruptcy
after eight months-long data breach, and huge expenses ($2.5 Million loans)
suffered when their system was hacked that led to the exposure of about 20
million Americans’ personal information [28]. In July 2019, Capital One
security breach leaked over 140,000 social security numbers of their clients,
personal information, and around 80,000 bank account numbers were exposed
[28]. The breach is estimated to cost the firm about $100 - $150 million. In
November 2019, personal data belonging to 1.2 billion unique people were
exposed on an Elasticsearch server that contains more than four terabytes of
data [1, 28].

###  3.1 Number of Records Compromised and Types of Data Breach (2005-2019)

In this research, we analyzed the Privacy Rights Clearinghouse (PRC) dataset
containing the report of data breaches for over 14 years. The dynamic nature
of the technology industry creates security loopholes that can be explored for
cyberattacks, some organization suffers a major loss, and often fold up due to
unprecedented events. Most security and data breaches span across all major
industries. The Privacy Rights Clearinghouse (PRC) Chronology of Data Breaches
reported more than ten billion records breached. Figure (2) discloses that
unauthorized access by a non-member of an organization (HACK) led to the
highest records breach with over one billion data compromised in both 2016 and
2017. It is also noteworthy that the pattern of attacks changed from physical
access to information and portable devices to unauthorized access to
information (i.e. Malware, DDOS) and unintended disclosure of information
between 2005 to 2012 and 2013 to 2019. Most cyberattacks are launched through
the use of sophisticated technology; this is also supported by the research
work conducted by Decary et al. [10].

Figure 2: Total Number of Compromised Records in Millions (2005-2019)

###  3.2 Occurrences of Cyberattacks Per States in the US

Several states in the US have started paying more attention to secure vital
information and also invested in improving data protection. The New York Penal
Law S. 156.10 punishes the act of knowing use of a computer to gain
unauthorized access to material, otherwise known as computer trespass. It was
reported that 31 states in the US enacted cybersecurity legislation that
significantly tackles issues relating to cybersecurity threats[16], some
states also have specific laws and policies regarding data breach. Most states
have enacted legislation that binds private or governmental entities to notify
individuals of security breaches of information involving personally
identifiable information after a certain period of days. Figure (3) further
illustrates the total number of records breached in different parts of the
United States. The total number of records breached indicated across the bar-
line is measured in billions with clusters in South and Northeast of US which
indicates a higher risk of data breach for organizations located in the
southeast and northeast. There is only a few research that has taken into
consideration the crime distribution rate in the united states using geo-
spatial tools in identifying the pattern of crime ad type of data breach. The
total number of records breached in different states ranges from zero to five
hundred million data. Only a few states have the total number of records
breached between one billion and two billion reports of data breach. The
research of Khey et al. [22] focused on the spatial distribution of data
breaches in the United States and risk profiling of vulnerabilities across
geographical locations.

Figure 3: US Map Showing Numbers of Breached Records Across Various Cities in
US (2005-2019)

##  4 Using the Deep Web (Hackers Forums) as a Source of Cyberthreat
Intelligences

In the era of big data, crawling the deep web is very important to get
insights into information embedded in the deep web; a large amount of data
exists in cyberspace and can be accessed through various interfaces. The
inability of search engines to crawl the deep web limits the information
provided to users. Oftentimes, users have to spend extra time and effort
querying a website manually to access the contents of a website. In their
study, Dixit et al. [24] developed an Ontology-based adaptive crawler for the
hidden web by mapping relationships between webpages using four components
(Ontology Builder, Hidden Web Miner, Result Processor, and Domain-Specific -
Flow Diagram). The research of Dong et al. [11] utilized a lightweight
framework to predict cyberthreats from darknet data. The method applied a
scrappy crawler with proxy VPN to parse the web page for relevant information
which is processed through verification techniques to generate cyberthreat
warnings. Azene et al. [42] utilized descriptive and predictive analytic
machine learning techniques on a data-set from darknet forum to discover
valuable cyberthreat intelligence. IBM Watson analytic tools and WEKA was used
to classify types of exploits targeted by hackers in the forums. Their results
showed that password cracker, key logger, and remote administration tools are
the most frequent tools targeted by influential authors in the forum Research
work in [21, 35] adopted a framework in predicting the future occurrence of a
cyberattack by mapping different data-sets from the dark web, socio-personal
and technical indicators with online cybersecurity report forums. A
combination of human and automated techniques was used to extract live data
from the dark web, insights gathered from online forums and community can be
used to predict malicious acts before they occur. The use of search engines
and spider services on the TOR network was used to collate a list of more than
100 malicious hacker websites.

###  4.1 Extracting information from the Deep web (Hackers Forum)

The deep web is unindexed, making it challenging to analyze the contents or
have access to information on darknet websites. Standard search engines cannot
harvest or crawl the contents of darknet websites. Two methods can be used to
crawl the contents of darknet websites. The first is through prior knowledge
base and queries, and second, a non-prior knowledge base with the use of
greedy queries. The attempts to close darknet websites always result in a
ripple effect. One of the major challenges encountered in crawling the deep
web forums was CAPTCHA authentication, username and password, invitation code,
and filling of query forms for automatic crawlers. The prototype system built-
in [19] tackled the challenge of query generation, empty page filtering, and
URL duplication. In this research, we utilized an open-source deep web dataset
that contains discussion forums from two darknet markets (Silkroad & Wall
street) extracted from the Arizona State University database [13]. The data
set contains over 128,000 posts from different discussion threads. Discussions
were organized in a thread topic, and other users initiate discussions based
on the thread title. The thread title is related to :

  1. 1.

Carding: Carding thread forums involve discussions on the trade of stolen
credit card details, identity theft, and currency counterfeiting. The forums
allow both vendors and customers to trade credit card details, bank account,
and other personal information online.

  2. 2.

Newbie: This is a discussion thread for new users and inexperienced hackers.
The discussion thread focuses mainly on welcoming new members, introducing
basic concepts and rules of engagements. Some other forums limit new users to
pending when they make payments to access advanced features.

  3. 3.

Scamming: Discussion in this threat is centered towards stealing of
information, fraudulent activities, and transactions. One major characteristic
of this thread is dishonest transactions and the listing of fake products.

  4. 4.

Hacking: Hackers forums consist of a wide variety of discussions ranging from
software leak, malware, credit cards, data breach, identity theft, and
financial information.

  5. 5.

Review Thread: The review thread in the forums represents a town hall meeting
where both vendors and customers discuss a variety of topics and share
knowledge base with each other. This thread also serves as a tool/page to
report users who have violated the forum’s policy or scammed other members of
the forum.

###  4.2 Deep Learning Algorithm for Threat Detection

Deep learning is a subset of machine learning methods based on neural networks
with representation learning. Deep Neural Network (DNN) is an artificial
neural network with multiple layers between the input and output layers.
Machine learning algorithms can be supervised, semi-supervised or
unsupervised. The Long Short-Term Memory (LSTM) neural network is a special
variant of Recurrent Neural Networks (RNN), which overcomes stability
bottlenecks encountered in traditional Recurrent Neural Networks [26]. RNN is
a class of Artificial Neural Network (ANN) that allows the use of previous
outputs as the current input with hidden layers. This special type of ANN
makes use of sequential and time-series information, RNN allows output from
previous steps to be used as input to the current step while having hidden
states. The unique feature of RNN is the sharing of parameters across each
layer of the network, RNNs share the same weight parameter within each layer
of the network and are adjusted through the processes of backpropagation and
gradient descent. RNN stores information in memory because their current
output is dependent on the previous computations. However, the major problem
associated with RNNs are vanishing and exploding gradients, LSTMs were
designed to overcome the bottlenecks of RNN by introducing new gates which
allow better control over the gradient flow and enable better preservation of
long-range dependencies [6, 5]. The LSTM architecture is different from other
deep learning architectures, LSTM model contains the memory cell and gates
which are essential to the architecture namely; the input gate, output gate,
and the forget gate. The LSTM regulates the flow of training information
through these gates by selectively adding information.

$\displaystyle i_{t}$ $\displaystyle=\sigma(W_{i}\cdot[h_{t-1},x_{t}]+b_{i})$
(1) $\displaystyle O_{t}$
$\displaystyle=\sigma(W_{o}\cdot[h_{t-1},x_{t}]+b_{o})$ (2) $\displaystyle
f_{t}$ $\displaystyle=\sigma(W_{f}\cdot[h_{t-1},x_{t}]+b_{f})$ (3)
$\displaystyle C_{t}$ $\displaystyle=f_{t}-C_{t-1}+i_{t}-\tilde{C}_{t}$ (4)
$\displaystyle h_{t}$ $\displaystyle=O_{t}-tanh(C_{t})$ (5)
$\displaystyle\tilde{C}_{t}$
$\displaystyle=tanh(W_{C}\cdot[h_{t-1},x_{t}]+b_{C})$ (6)

The gates in LSTM cell enable it to preserve constant error that can be
backpropagated, the input gate is represented by $i$, output gate by $o$, and
forget gate by $f$. The memory cell represented as $C$ accumulates the state
of information, and the cell output is given by $h$, while the cell input at
time $t$ is denoted as $x_{t}$. $W$ are the weights for each of the gates and
$\tilde{C}$ is the updated cell state. These states are propagated ahead
through the network, the forget gate plays a crucial role in reducing over-
fitting by not retaining all information from the previous time steps. See
figure(4)

Figure 4: Simple illustration of LSTM cells

The performance of machine learning methods is heavily dependent on the choice
of data representation (features) to which they are applied. In this study, we
utilized LSTM to train the deep web data because LSTM performs better when
training a large time-series corpus. In this research, we curated a unique
features extraction pipeline to extract time-series data from hackers’ forums.
The extracted features served as the input of the first hidden layer of our
LSTM model architecture.

##  5 Results and Discussion

We used NER for features extraction and tokenizers to convert raw text to
matrices. The average length of each text in the dataset is 150. The maximum
word length is a sequence of 250 words. The input texts were padded to
maintain a uniform sequence length. The model was trained on twelve epochs
with an accuracy of 91% on train and 87% on validation data.

Figure 5: Accuracy of LSTM Model

Figure (5) shows the accuracy of the trained LSTM model on both train and
validation data. We compared the performance of the developed model with other
related work in predicting cyberthreat and identifying vulnerability exposures
using deep web data. See table (1). The research of Dong et al. utilized a
lightweight framework to predict cyberthreats from darknet data. The
researchers monitored eight darknet marketplaces for a period of five months.
Scrappy crawler was used with Polipo and Vidalia proxies to scrape information
related to cybersecurity. Their model was able to predict 145 existing
threats, 35 new threats, and newly developed hacking tools in the darknet
marketplace. The researchers experimented with SVM. SVM has the highest
accuracy of 81% and precision of 90% in detecting threats [11]. The research
of Azene et al. [42] on cyberthreat discovery from the dark web utilized
descriptive and predictive analytic machine learning techniques to discover
valuable cyberthreat intelligence. Their classification model with Naive Bayes
had an accuracy of 65% on web exploit, Random Tree 84% on system exploits, and
Random forest classifier has the highest accuracy of 97% on system exploits.
The research of Ashok et al.[9] developed a novel approach for predicting
cyberthreat leveraging on hackers’ sentiments. The researchers analyzed over
four hundred thousand posts in different hacker forums. Their experimental
results in organizations malicious email had the precision of 76% and recall
of 63%. The deep learning algorithm (LSTM) developed in this research
outperforms every other base model algorithm reported in related research. As
showcased in 1), LSTM had an accuracy of 94%, precision of 90%, and a recall
of 91%, which is the best result of deep learning algorithms in identifying
anomalous cyberthreat texts and predicting vulnerability exposures leveraging
on the discussion in hackers forum.

Table 1: Comparison of other baseline models

Models |  Algorithm |  Accuracy |  Precision |  Recall |  F1/ROC
---|---|---|---|---|---
My Model |  LSTM |  0.94 |  0.90 |  0.91 |  0.91
|  Random Forest (RFC) |  0.80 |  0.95 |  0.95 |  0.75
Dong et al. [11] |  SVM |  0.81 |  0.90 |  - |  -
Azene et al[42] |  Naive Bayes |  0.65 |  - |  - |  0.85
|  Random Tree |  0.84 |  - |  - |  0.71
|  Random Forest |  0.97 |  - |  - |  0.91
Arora et al. [4] |  Random Forest(RFC) |  0.80 |  0.81 |  0.80 |  0.79

##  6 Conclusion

Cyberspace is a global environment where innovation, creativity, and access to
information have almost no limit. Most businesses rely on web-based
applications such as Zoom, WebEx, Microsoft teams, and other communication
platforms. Actors with malicious intent could exploit the vulnerabilities of
these systems to steal personal information and engage in illegal activities.
The occurrence of Cyberattacks has a significant impact on organizational
activities and people’s privacy. Based on the increased rate of occurrence
gathering social and technical indicators useful for threat intelligence from
hackers forum is imperative in developing tools for detecting cyberthreats
before the actual occurrence. The model developed in this research can be
deployed by security analysts and law enforcement agencies in identifying
existing cyberthreats and predicting potential vulnerability exposures using
deep learning algorithms. The results of our experiments indicate that our
methodology is more effective in detecting cyberthreats with 90% of accuracy.
This paper also identifies potential areas for further research in developing
a methodology to scrape information from deep web forums securely using
automatic smart crawlers. As part of our future studies, we will use
geospatial analysis in profiling risk factors of cyberthreat in various
regions in the United States to establish correlations between deep web and
surface web virtual ecosystem in orchestrating criminal activities by
nefarious netizens.

## Acknowledgment

This material is based upon work supported by the National Science Foundation
under Grant No. (CNS-1801593). Any opinions, findings, and conclusions, or
recommendations expressed in this material are those of the author(s) and do
not necessarily reflect the views of the National Science Foundation.

## References

  * [1] Phone numbers of millions of Facebook users exposed online, https://www.cisomag.com/unprotected-database-exposes-millions-of-facebook-users-contact-numbers/ https://www.cisomag.com/elasticsearch-server-exposed-1-2-billion-people-data/
  * [2] Adewopo, V., Gonen, B., Adewopo, F.: Exploring open source information for cyber threat intelligence. In: 2020 IEEE International Conference on Big Data (Big Data). pp. 2232–2241. IEEE (2020)
  * [3] Adewopo, V., Gonen, B., Varlioglu, S., Ozer, M.: Plunge into the underworld: A survey on emergence of darknet. In: 2019 International Conference on Computational Science and Computational Intelligence (CSCI). pp. 155–159. IEEE (2019)
  * [4] Arora, T., Sharma, M., Khatri, S.K.: Detection of cyber crime on social media using random forest algorithm. In: 2019 2nd International Conference on Power Energy, Environment and Intelligent Control (PEEIC). pp. 47–51. IEEE (2019)
  * [5] Azumah, S.W., Elsayed, N., Adewopo, V., Zaghloul, Z.S., Li, C.: A deep lstm based approach for intrusion detection iot devices network in smart home. In: 2021 IEEE 7th World Forum on Internet of Things (WF-IoT). pp. 836–841 (2021). https://doi.org/10.1109/WF-IoT51360.2021.9596033
  * [6] Bouktif, S., Fiaz, A., Ouni, A., Serhani, M.A.: Optimal deep learning lstm model for electric load forecasting using feature selection and genetic algorithm: Comparison with machine learning approaches. Energies 11(7), 1636 (2018)
  * [7] Chaabane, A., Manils, P., Kaafar, M.A.: Digging into anonymous traffic: A deep analysis of the Tor anonymizing network. In: Proceedings - 2010 4th International Conference on Network and System Security, NSS 2010. pp. 167–174 (2010). https://doi.org/10.1109/NSS.2010.47
  * [8] Cuthbertson, A.: Facebook hack: People’s accounts appear for sale on dark web — The Independent. Tech. rep., https://www.independent.co.uk/life-style/gadgets-and-tech/news/facebook-hack-data-dark-web-login-details-cost-dream-market-a8564671.html
  * [9] Deb, A., Lerman, K., Ferrara, E.: Predicting cyber-events by leveraging hacker sentiment. Information (Switzerland) 9(11), 280 (nov 2018). https://doi.org/10.3390/info9110280, www.mdpi.com/journal/information
  * [10] Décary-Hétu, D., Giommoni, L.: Do police crackdowns disrupt drug cryptomarkets? a longitudinal analysis of the effects of operation onymous. Crime, Law and Social Change 67(1), 55–75 (2017)
  * [11] Dong, F., Yuan, S., Ou, H., Liu, L.: New cyber threat discovery from darknet marketplaces. In: 2018 IEEE Conference on Big Data and Analytics (ICBDA). pp. 62–67 (2018). https://doi.org/10.1109/ICBDAA.2018.8629658
  * [12] Dragut, E.C., Meng, W., Yu, C.T.: Deep web query interface understanding and integration. Synthesis Lectures on Data Management 7(1), 1–168 (2012)
  * [13] Du, P.Y., Zhang, N., Ebrahimi, M., Samtani, S., Lazarine, B., Arnold, N., Dunn, R., Suntwal, S., Angeles, G., Schweitzer, R., et al.: Identifying, collecting, and presenting hacker community data: Forums, irc, carding shops, and dnms. In: 2018 IEEE international conference on intelligence and security informatics (ISI). pp. 70–75. IEEE (2018)
  * [14] Elahi, T., Bauer, K., Alsabah, M., Dingledine, R., Goldberg, I.: Changing of the Guards: A Framework for Understanding and Improving Entry Guard Selection in Tor (2012)
  * [15] Everett, C.: Should the dark net be taken out? Network Security 2015(3), 10–13 (mar 2015). https://doi.org/10.1016/S1353-4858(15)30018-0
  * [16] Fonseca, B., Rosen, J.D.: Cybersecurity in the us: Major trends and challenges. In: The New US Security Agenda, pp. 87–106. Springer (2017)
  * [17] Gehl, R., McKelvey, F.: Bugging out: darknets as parasites of large-scale media objects. Media, Culture and Society 41(2), 219–235 (mar 2019). https://doi.org/10.1177/0163443718818379
  * [18] Guitton, C.: A review of the available content on Tor hidden services: The case against further development. Computers in Human Behavior 29(6), 2805–2815 (2013). https://doi.org/10.1016/j.chb.2013.07.031
  * [19] He, Y., Xin, D., Ganti, V., Rajaraman, S., Shah, N.: Crawling deep web entity pages. In: WSDM 2013 - Proceedings of the 6th ACM International Conference on Web Search and Data Mining. pp. 355–364 (2013). https://doi.org/10.1145/2433396.2433442
  * [20] Irudayaraj, P.: Distributed Hash Tables. pp. 403–403 (2011). https://doi.org/10.1109/icccn.2006.286309, https://blog.keep.network/distributed-hash-tables-49721094403d
  * [21] Kadoguchi, M., Hayashi, S., Hashimoto, M., Otsuka, A.: Exploring the dark web for cyber threat intelligence using machine leaning. In: 2019 IEEE International Conference on Intelligence and Security Informatics, ISI 2019. pp. 200–202. Institute of Electrical and Electronics Engineers Inc. (jul 2019). https://doi.org/10.1109/ISI.2019.8823360
  * [22] Khey, D.N., Sainato, V.A.: Examining the correlates and spatial distribution of organizational data breaches in the united states. Security Journal 26(4), 367–382 (2013)
  * [23] Mansfield-Devine, S.: Darknets. Computer Fraud and Security 2009(12), 4–6 (dec 2009). https://doi.org/10.1016/S1361-3723(09)70150-2
  * [24] Manvi, M., Dixit, A., Bhatia, K.K.: Design of an ontology based adaptive crawler for hidden Web. In: Proceedings - 2013 International Conference on Communication Systems and Network Technologies, CSNT 2013. pp. 659–663 (2013). https://doi.org/10.1109/CSNT.2013.140
  * [25] Mod, D.: Leaked Emails and Messages from Russian Officials Published by the Darknet Site — Dark Web Link (2018), https://www.darkweblink.com/leaked-emails-and-messages-from-russian-officials-published-by-the-darknet-site/
  * [26] Mohan, A., Daniel, D., Chertkov, M., Livescu, D.: Compressed convolutional lstm: An efficient deep learning framework to model high fidelity 3d turbulence. arXiv preprint arXiv:1903.00033 (2019)
  * [27] More, D.: 5 Biggest Data Breaches of 2020 (So Far) - Security Boulevard (2020), https://securityboulevard.com/2020/08/5-biggest-data-breaches-of-2020-so-far/
  * [28] Okazaki, S., Eisend, M., Plangger, K., de Ruyter, K., Grewal, D.: Understanding the strategic consequences of customer privacy concerns: A meta-analytic review. Journal of Retailing (2020)
  * [29] Owen, G., Savage, N.: The Tor Dark Net (2015)
  * [30] Ozer, M., Varlioglu, S., Gonen, B., Adewopo, V., Elsayed, N., Zengin, S.: Cloud incident response: Challenges and opportunities. In: 2020 International Conference on Computational Science and Computational Intelligence (CSCI). pp. 49–54 (2020). https://doi.org/10.1109/CSCI51800.2020.00015
  * [31] Ozer, M., Varlioglu, S., Gonen, B., Bastug, M.: A prevention and a traction system for ransomware attacks. In: 2019 International Conference on Computational Science and Computational Intelligence (CSCI). pp. 150–154. IEEE (2019), https://arxiv.org/pdf/2001.02282.pdf
  * [32] Perez, J.: Unprotected server held 419 million Facebook user records including phone numbers. — Komando.com, https://www.komando.com/happening-now/593814/another-huge-facebook-data-breach-exposes-419-million-user-records
  * [33] Plantin, J.C., Lagoze, C., Edwards, P.N., Sandvig, C.: Infrastructure studies meet platform studies in the age of Google and Facebook. New Media & Society 20(1), 293–310 (jan 2018). https://doi.org/10.1177/1461444816661553, http://journals.sagepub.com/doi/10.1177/1461444816661553
  * [34] Ponemon Institute: 2019 Cost of a Data Breach Report — IBM Security. Tech. rep. (2019), https://databreachcalculator.mybluemix.net/?_ga=2.44171747.616715455.1574429761-1014862204.1574429761 https://www.ibm.com/downloads/cas/ZBZLY7KL
  * [35] Robertson, J., Diab, A., Marin, E., Nunes, E., Paliath, V., Shakarian, J., Shakarian, P.: Darkweb cyber threat intelligence mining (2017). https://doi.org/10.1017/9781316888513
  * [36] Schäfer, M., Fuchs, M., Strohmeier, M., Engel, M., Liechti, M., Lenders, V.: Blackwidow: Monitoring the dark web for cyber security information. In: 2019 11th International Conference on Cyber Conflict (CyCon). vol. 900, pp. 1–21. IEEE (2019)
  * [37] Sun, Y., Edmundson, A., University, P., Vanbever, L., Zürich, E., Li, O., Rexford, J., Chiang, M., Mittal, P., Vanbever ETH Zurich Oscar Li, L.: RAPTOR: Routing Attacks on Privacy in Tor. https://www.usenix.org/conference/usenixsecurity15/technical-sessions/presentation/sun
  * [38] Tsai, H., Harwood, A.: A scalable anonymous server overlay network. In: Proceedings - International Conference on Advanced Information Networking and Applications, AINA. vol. 1, pp. 973–978 (2006). https://doi.org/10.1109/AINA.2006.55
  * [39] Upadhyaya, R., Jain, A.: Cyber ethics and cyber crime: A deep dwelved study into legality, ransomware, underground web and bitcoin wallet. Proceeding - IEEE International Conference on Computing, Communication and Automation, ICCCA 2016 pp. 143–148 (2017). https://doi.org/10.1109/CCAA.2016.7813706
  * [40] Volz, D., McMillan, R.: Federal judiciary’s systems likely breached in solarwinds hack (Jan 2021), https://www.wsj.com/articles/federal-judiciarys-systems-likely-breached-in-solarwinds-hack-11610040175
  * [41] Wang, M., Wang, X., Shi, J., Tan, Q., Gao, Y., Chen, M., Jiang, X.: Who are in the darknet? Measurement and analysis of Darknet person attributes. In: Proceedings - 2018 IEEE 3rd International Conference on Data Science in Cyberspace, DSC 2018. pp. 948–955. Institute of Electrical and Electronics Engineers Inc. (jul 2018). https://doi.org/10.1109/DSC.2018.00151
  * [42] Zenebe, A., Shumba, M., Carillo, A., Cuenca, S.: Cyber threat discovery from dark web. In: EPiC Series in Computing. vol. 64, pp. 174–183. EasyChair (2019). https://doi.org/10.29007/nkfk
  * [43] Zhao, F., Zhou, J., Nie, C., Huang, H., Jin, H.: SmartCrawler: A two-stage crawler for efficiently harvesting deep-web interfaces. IEEE Transactions on Services Computing 9(4), 608–620 (jul 2016). https://doi.org/10.1109/TSC.2015.2414931


estimators $\hat{\eta}_{pos,s}$ and $\hat{\eta}_{neg,s}$ and given as follows:

$\hat{\eta}_{avg,n}=\frac{\hat{\eta}_{pos,s}+\hat{\eta}_{neg,s}}{2}$

This estimator integrates to 1 but for test data points it is prone to lying
outside the interval $[0,1]$.

The estimators which does not lie inside the interval $[0,1]$ can be truncated
to lie inside the interval. However, there performance might be affected. A
compact form of the above estimators is available in Table 15.

$\hat{\eta}$ |  $\hat{\eta}(\mathbf{x})\in[0,1]$ for $\mathbf{x}$ in test set
---|---
$\hat{\eta}_{pos}(\mathbf{x})=\hat{w}_{pos}(\mathbf{x})\hat{\pi}$ | May not
$\hat{\eta}_{neg}(\mathbf{x})=1-\hat{w}_{neg}(\mathbf{x})(1-\hat{\pi})$ | May not
$\hat{\eta}_{norm}=\frac{\hat{\eta}_{pos}}{\hat{\eta}_{pos}+\hat{\eta}_{neg}}$ | Yes
$\hat{\eta}_{pos,s}=\frac{\hat{\eta}_{pos}}{\hat{\eta}_{max,p}}$ | May not
$\hat{\eta}_{neg,s}=\frac{\hat{\eta}_{neg}}{\hat{\eta}_{max,n}}$ | May not
$\hat{\eta}_{norm,s}=\frac{\hat{\eta}_{pos,s}}{\hat{\eta}_{pos,s}+\hat{\eta}_{neg,s}}$ | Yes
$\hat{\eta}_{Cnorm,s}=\frac{\hat{\eta}_{neg,s}}{\hat{\eta}_{neg,s}+\hat{\eta}_{pos,s}}$ | Yes
$\hat{\eta}_{avg,n}=\frac{\hat{\eta}_{pos,s}+\hat{\eta}_{neg,s}}{2}$ | May not
Table 15: Details about various KLIEP based $\eta$ estimates and whether they
are valid probability or not. We only check for $\hat{\eta}\in[0,1]$ as all
the estimates integrate to 1.

The comparison between these estimators w.r.t. various measures is done in
Section 0.F.7. In our implementation, we have used Gaussian kernels as the
basis functions with centres as the sample points from the numerator density
(positive/negative class). The kernel width $\sigma$ is tuned by KLIEP’s
inbuilt CV procedure where the criterion is maximum objective value of the KL
based optimization problem. The possible candidates for kernel width $\sigma$
to be used in KLIEP’s CV procedure were generated through Centile method
described in Pseudocode Centile method.

###  0.F.7 An experimental comparative study of various in-class probability
estimation methods

In this section, we compare the performance of the $\eta$ estimation methods
described in Section 0.F.1, 0.F.2 and 0.F.6. We implemented Algorithm
$\boldsymbol{(f,\hat{\eta})}$ with squared loss, logistic loss and modified
squared loss and LSPC, KLIEP and $k$-NN on 3 synthetic datasets so that we can
comment on the quality of $\eta$ estimates. In case of $k$-NN, the value of
number of neighbours $k$ to be used is selected by cross validating w.r.t. the
measure to be evaluated on. The dataset generation scheme is adapted from [4].
We consider 9 performance measures out of which first 5 (MSE, RMSE, MAD, MD,
KL) are purely for $\eta$ estimation, next 2 (Acc, BS) are for the scenario
when these estimates are used in label prediction and the last 2 measure the
algorithms estimation capability at the boundary i.e., maximum and minimum.
For an estimate $\hat{\eta}$, when the train data is
$\\{(\mathbf{x}_{i},y_{i})\\}_{i=1}^{m_{tr}}$ and test data set is
$\\{(\mathbf{x}_{i},y_{i})\\}_{i=1}^{m_{te}}$, the details about all these
measures are given below:

  1. 1.

Mean squared error (MSE):
$MSE=\frac{1}{m_{te}}\sum\limits_{i=1}^{m_{te}}(\hat{\eta}(\mathbf{x}_{i})-\eta(\mathbf{x}_{i}))^{2}$

  2. 2.

Root mean squared error (RMSE): $RMSE=\sqrt[]{MSE}$

  3. 3.

Mean absolute deviation (MAD):
$MAD=\frac{1}{m_{te}}\sum\limits_{i=1}^{m_{te}}|\hat{\eta}(\mathbf{x}_{i})-\eta(\mathbf{x}_{i})|$

  4. 4.

Mean deviation (MD):
$MD=\frac{1}{m_{te}}\sum\limits_{i=1}^{m_{te}}(\hat{\eta}(\mathbf{x}_{i})-\eta(\mathbf{x}_{i}))$

  5. 5.

Averaged Kullback-Leibler divergence (KL): Since $\eta$ is a density, we can
compute the KL between its estimate and true value as follows:

$\displaystyle KL$ $\displaystyle=$
$\displaystyle\frac{1}{m}\sum\limits_{i=1}^{m_{te}}KL[\eta(\mathbf{x}_{i})\|\hat{\eta}(\mathbf{x}_{i})]$
$\displaystyle=$
$\displaystyle\frac{1}{m}\sum\limits_{i=1}^{m_{te}}\eta(\mathbf{x}_{i})\log\frac{\eta(\mathbf{x}_{i})}{\hat{\eta}(\mathbf{x}_{i})}+(1-\eta(\mathbf{x}_{i}))\log\frac{1-\eta(\mathbf{x}_{i})}{1-\hat{\eta}(\mathbf{x}_{i})}$

  6. 6.

Accuracy (Acc):
$Acc=\frac{1}{m_{te}}\sum\limits_{i=1}^{m_{te}}\mathbf{1}_{[(2\hat{\eta}(\mathbf{x}_{i})-1)y_{i}\leq
0]}$

  7. 7.

Brier Score (BS): Brier score is a proper score function which measures the
accuracy of probabilistic predictions and is defined below:

$BS=\frac{1}{m_{te}}\sum\limits_{i=1}^{m_{te}}(\hat{\eta}(\mathbf{x}_{i})-y_{i})^{2}$

  8. 8.

DiffMax: In some applications, one requires how good the estimates are at the
boundary. In such cases, we have to compare the maximum value of true $\eta$
to the maximum value of $\hat{\eta}$ on the training set. This measure plays
an important role when one is interested in estimating noise rates as some
noise estimating schemes use corrupted ${\eta}$ values [9]. Therefore,

$DiffMax:=\max_{i=1,\ldots,m_{tr}}\hat{\eta}(\mathbf{x}_{i})-\max_{i=1,\ldots,m_{tr}}\eta(\mathbf{x}_{i})$

  9. 9.

DiffMin: Similar to DiffMax, this measure is also a deciding factor for
quality of noise rates.

$DiffMin:=\min_{i=1,\ldots,m_{tr}}\hat{\eta}(\mathbf{x}_{i})-\min_{i=1,\ldots,m_{tr}}\eta(\mathbf{x}_{i})$

The synthetic dataset is partitioned into training and test set with $80$-$20$
split. The first 7 performance measures are computed on the test set and the
last two on the training set. To make the comparisons more robust, this
process is repeated for 10 trials. Finally, the reported values are
performance measures averaged over the 10 trials along with the standards
deviation.

####  0.F.7.1 Synthetic dataset 1 ($Syn_{2D}$)

This is a 2 dimensional dataset. We first generate 1000 binary class labels
$Y$ from Bernoulli distribution with parameter $p=0.5$. Then a 2-dimensional
feature vector $\mathbf{X}$ for each label is drawn from two different
Gaussian distributions: $\mathbf{X}|Y=1\sim N([1.2,0],\Sigma)$ $\&$
$\mathbf{X}|Y=-1\sim N([-1.2,0],\Sigma)$ where $\Sigma=[[1,0.4];[0.4,1]]$.
Table 16, 17 and 18 present the values of various measures for this dataset.

| ${\eta}$ estimated from
---|---
| Performance
---
measures
| Algorithm $\boldsymbol{(f,\hat{\eta})}$
---
with $l_{sq}$
| Algorithm $\boldsymbol{(f,\hat{\eta})}$
---
with $l_{log}$
| Algorithm $\boldsymbol{(f,\hat{\eta})}$
---
with $l_{msq}$
LSPC | k-NN
MSE |  0.02$\pm$0.002 |  0.0008$\pm$0.0007 |  0.001$\pm$0.0003 |  0.003$\pm$0.002 |  0.013 $\pm$ 0.003
RMSE |  0.144$\pm$0.007 |  0.025$\pm$0.012 |  0.034$\pm$0.004 |  0.054$\pm$0.02 |  0.114 $\pm$ 0.014
MAD |  0.123$\pm$0.005 |  0.018$\pm$0.010 |  0.025$\pm$0.002 |  0.037$\pm$0.015 |  0.066 $\pm$ 0.008
MD |  -0.0032$\pm$0.011 |  -0.0076$\pm$0.001 |  -0.007$\pm$0.002 |  -0.007$\pm$0.004 |  -0.007 $\pm$ 0.005
KL | NA |  0.004$\pm$0.004 | NA | NA | NA
Acc |  0.911$\pm$0.016 |  0.91$\pm$0.017 |  0.91$\pm$0.018 |  0.91$\pm$0.01 |  0.906 $\pm$ 0.016
BS |  0.0804$\pm$0.004 |  0.0651$\pm$0.010 |  0.0657$\pm$0.007 |  0.06809$\pm$0.0068 |  0.077 $\pm$ 0.010
DiffMax |  0.531$\pm$0.04 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0
DiffMin |  -0.497$\pm$0.042 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0
Table 16: Performance of estimates obtained by Algorithm $\boldsymbol{(f,\hat{\eta})}$ with squared loss, logistic loss and modified squared loss and LSPC and $k$-NN on $Syn_{2D}$ dataset. The cells with NA means that some estimate was 0 and the KL couldn’t be defined. Cells with $0.0$ entry means that the value is less than $1\times 10^{-5}$. | KLIEP based estimator used
---|---
| Performance
---
measures
$\mathbf{\hat{\eta}_{pos}}$ | $\mathbf{\hat{\eta}_{neg}}$ | $\mathbf{\hat{\eta}_{norm}}$
MSE |  0.012$\pm$ 0.003 |  0.014$\pm$0.007 |  0.001$\pm$ 0.00
RMSE |  0.108$\pm$ 0.015 |  0.118$\pm$0.027 |  0.034$\pm$0.009
MAD |  0.062$\pm$ 0.006 |  0.069$\pm$0.016 |  0.021$\pm$0.005
MD |  -0.006$\pm$0.008 |  –0.006$\pm$0.013 |  -0.003$\pm$0.005
KL | NA | NA |  0.005$\pm$ 0.002
Acc |  0.89 $\pm$ 0.02 |  0.90 $\pm$ 0.01 |  0.91 $\pm$ 0.02
BS |  0.075 $\pm$ 0.011 |  0.074$\pm$ 0.006 |  0.065$\pm$0.006
DiffMax |  0.16$\pm$0.02 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0
DiffMin |  0.0 $\pm$0.0 |  -0.16 $\pm$0.06 |  0.0 $\pm$0.0
Table 17: Performance of $\hat{\eta}_{pos},\hat{\eta}_{neg},\hat{\eta}_{norm}$, KLIEP based estimates on $Syn_{2D}$ dataset. The cells with NA means that some estimate was $0$ and the KL couldn’t be defined. Cells with $0.0$ entry means that the value is less than $1\times 10^{-5}$. | KLIEP based estimator used
---|---
| Performance
---
measures
$\mathbf{\hat{\eta}_{pos,s}}$ | $\mathbf{\hat{\eta}_{neg,s}}$ | $\mathbf{\hat{\eta}_{norm,s}}$ | $\mathbf{\hat{\eta}_{Cnorm,s}}$ | $\mathbf{\hat{\eta}_{avg,s}}$
MSE |  0.009$\pm$0.004 |  0.019$\pm$0.01 |  0.001$\pm$0.00 |  0.178$\pm$0.004 |  0.009$\pm$0.002
RMSE |  0.096$\pm$0.019 |  0.133$\pm$0.035 |  0.036$\pm$0.009 |  0.133$\pm$0.013 |  0.097$\pm$0.01
MAD |  0.095$\pm$0.008 |  0.072$\pm$0.027 |  0.021$\pm$0.005 |  0.105$\pm$0.014 |  0.072$\pm$0.01
MD |  -0.08$\pm$0.008 |  0.06$\pm$0.027 |  -0.003$\pm$0.007 |  -0.019$\pm$0.018 |  -0.012$\pm$0.014
KL |  0.073 $\pm$0.002 | NA |  0.005 $\pm$ 0.002 | NA |  0.047$\pm$0.008
Acc |  0.87 $\pm$ 0.02 |  0.89 $\pm$ 0.02 |  0.90 $\pm$ 0.01 |  0.90 $\pm$0.01 |  0.90 $\pm$ 0.02
BS |  0.086$\pm$0.010 |  0.084$\pm$0.009 |  0.066$\pm$0.006 |  0.083$\pm$0.006 |  0.0748$\pm$0.005
DiffMax |  0.0 $\pm$0.0 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0 |  -0.003$\pm$0.00
DiffMin |  0.0 $\pm$0.0 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0 |  -0.005$\pm$0.004
Table 18: Performance of
$\hat{\eta}_{pos,s},\hat{\eta}_{neg,s},\hat{\eta}_{norm,s},\hat{\eta}_{Cnorm,s},\hat{\eta}_{avg,s}$,
KLIEP based estimates on $Syn_{2D}$ dataset. The cells with NA means that some
estimate was $0$ and the KL couldn’t be defined. Cells with $0.0$ entry means
that the value is less than $1\times 10^{-5}$.

####  0.F.7.2 Synthetic dataset 2 ($Syn_{3D}$)

This is a 3 dimensional dataset. We first generate 1000 binary class labels
$Y$ from Bernoulli distribution with parameter $p=0.5$. Then a 3-dimensional
feature vector $\mathbf{X}$ for each label is drawn from two different
Gaussian distributions: $\mathbf{X}|Y=1\sim N([1.3,0,0],\Sigma)$ $\&$
$\mathbf{X}|Y=-1\sim N([-1.3,0,0],\Sigma)$ where
$\Sigma=[[1,0.1,0.1];[0.1,1,0.1];[0.1,0.1,1]]$. Table 19, 20 and 21 present
the values of various measures for this dataset.

| ${\eta}$ estimated from
---|---
| Performance
---
measures
| Algorithm $\boldsymbol{(f,\hat{\eta})}$
---
with $l_{sq}$
| Algorithm $\boldsymbol{(f,\hat{\eta})}$
---
with $l_{log}$
| Algorithm $\boldsymbol{(f,\hat{\eta})}$
---
with $l_{msq}$
LSPC | k-NN
MSE |  0.02$\pm$0.001 |  0.0008$\pm$0.0007 |  0.001$\pm$0.0003 |  0.003$\pm$0.002 |  0.0108 $\pm$ 0.0022
RMSE |  0.145$\pm$0.007 |  0.026$\pm$0.007 |  0.035$\pm$0.003 |  0.064$\pm$0.003 |  0.103 $\pm$ 0.0102
MAD |  0.126$\pm$0.003 |  0.018$\pm$0.007 |  0.025$\pm$0.002 |  0.045$\pm$0.015 |  0.0612 $\pm$ 0.005
MD |  -0.016$\pm$0.008 |  -0.0078$\pm$0.002 |  -0.006$\pm$0.0036 |  -0.009$\pm$0.005 |  -0.008 $\pm$ 0.007
KL | NA |  0.004$\pm$0.003 | NA | NA | NA
Acc |  0.92$\pm$0.014 |  0.923$\pm$0.015 |  0.924$\pm$0.014 |  0.923$\pm$0.016 |  0.929 $\pm$ 0.0197
BS |  0.0820$\pm$0.006 |  0.0638$\pm$0.008 |  0.0663$\pm$0.008 |  0.0692$\pm$0.008 |  0.0667 $\pm$ 0.0127
DiffMax |  0.628$\pm$0.01 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0
DiffMin |  -0.458$\pm$0.01 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0
Table 19: Performance of estimates obtained by Algorithm $\boldsymbol{(f,\hat{\eta})}$ with squared loss, logistic loss and modified squared loss and LSPC and $k$-NN on $Syn_{3D}$ dataset. The cells with NA means that some estimate was 0 and the KL couldn’t be defined. Cells with $0.0$ entry means that the value is less than $1\times 10^{-5}$. | KLIEP based estimator used
---|---
| Performance
---
measures
$\mathbf{\hat{\eta}_{pos}}$ | $\mathbf{\hat{\eta}_{neg}}$ | $\mathbf{\hat{\eta}_{norm}}$
MSE |  0.022$\pm$ 0.003 |  0.025$\pm$0.005 |  0.002$\pm$ 0.00
RMSE |  0.146$\pm$ 0.010 |  0.158$\pm$0.015 |  0.050$\pm$0.005
MAD |  0.097$\pm$ 0.007 |  0.103$\pm$0.007 |  0.040$\pm$0.005
MD |  -0.0004$\pm$0.009 |  –0.009$\pm$0.019 |  -0.006$\pm$0.004
KL | NA | NA |  0.014$\pm$ 0.003
Acc |  0.90 $\pm$ 0.02 |  0.89 $\pm$ 0.0 |  0.92 $\pm$ 0.016
BS |  0.0793 $\pm$ 0.0077 |  0.084$\pm$ 0.010 |  0.067$\pm$0.007
DiffMax |  0.34$\pm$0.02 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0
DiffMin |  0.0002 $\pm$0.0 |  -0.344 $\pm$0.02 |  0.006 $\pm$0.002
Table 20: Performance of $\hat{\eta}_{pos},\hat{\eta}_{neg},\hat{\eta}_{norm}$, KLIEP based estimates on $Syn_{3D}$ dataset. The cells with NA means that some estimate was $0$ and the KL couldn’t be defined. Cells with $0.0$ entry means that the value is less than $1\times 10^{-5}$. | KLIEP based estimator used
---|---
| Performance
---
measures
$\mathbf{\hat{\eta}_{pos,s}}$ | $\mathbf{\hat{\eta}_{neg,s}}$ | $\mathbf{\hat{\eta}_{norm,s}}$ | $\mathbf{\hat{\eta}_{Cnorm,s}}$ | $\mathbf{\hat{\eta}_{avg,s}}$
MSE |  0.014$\pm$0.004 |  0.044$\pm$0.001 |  0.03$\pm$0.005 |  0.041$\pm$0.003 |  0.024$\pm$0.002
RMSE |  0.118$\pm$0.014 |  0.209$\pm$0.021 |  0.051$\pm$0.005 |  0.203$\pm$0.01 |  0.154$\pm$0.009
MAD |  0.164$\pm$0.008 |  0.141$\pm$0.014 |  0.042$\pm$0.005 |  0.183$\pm$0.010 |  0.136$\pm$0.009
MD |  -0.151$\pm$0.004 |  0.122$\pm$0.0151 |  -0.008$\pm$0.015 |  -0.018$\pm$0.009 |  -0.014$\pm$0.008
KL |  0.150 $\pm$0.011 | NA |  0.015 $\pm$ 0.003 | NA |  0.102$\pm$0.009
Acc |  0.8535 $\pm$ 0.01 |  0.83 $\pm$ 0.02 |  0.92 $\pm$ 0.01 |  0.92 $\pm$0.01 |  0.92 $\pm$ 0.02
BS |  0.1147$\pm$0.006 |  0.1144$\pm$0.011 |  0.0677$\pm$0.0.007 |  0.108$\pm$0.005 |  0.0905$\pm$0.005
DiffMax |  0.0 $\pm$0.0 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0 |  -0.006$\pm$0.02
DiffMin |  0.00 $\pm$0.0 |  0.0 $\pm$0.0 |  0.0006 $\pm$0.0 |  0.0 $\pm$0.0 |  -0.008$\pm$0.002
Table 21: Performance of
$\hat{\eta}_{pos,s},\hat{\eta}_{neg,s},\hat{\eta}_{norm,s},\hat{\eta}_{Cnorm,s},\hat{\eta}_{avg,s}$,
KLIEP based estimates on $Syn_{3D}$ dataset. The cells with NA means that some
estimate was $0$ and the KL couldn’t be defined. Cells with $0.0$ entry means
that the value is less than $1\times 10^{-5}$.

####  0.F.7.3 Synthetic dataset 3 ($Syn_{10D}$)

This is a 10 dimensional dataset. We first generate 1000 binary class labels
$Y$ from Bernoulli distribution with parameter $p=0.5$. Then, a 10-dimensional
feature vector $\mathbf{X}$ for each label is drawn from two different
Gaussian distributions: $\mathbf{X}|Y=1\sim N(\boldsymbol{\mu}_{+},\Sigma)$
$\&$ $\mathbf{X}|Y=-1\sim N(\boldsymbol{\mu}_{-},\Sigma)$ where $\Sigma$ is
such that all the diagonal elements (variances) are $1$ and the non-diagonal
elements (covariances) are $10$. Here, the means are
$\boldsymbol{\mu}_{+}=[2,1.2,2,2.1,2,2,2,0.2,2,3.6]$ and
$\boldsymbol{\mu}_{+}=[2,-1.2,2,-2.1,2,2,2,-0.2,2,-3.6].$ Table 22, 23 and 24
present the values of various measures for this dataset.

| ${\eta}$ estimated from
---|---
| Performance
---
measures
| Algorithm $\boldsymbol{(f,\hat{\eta})}$
---
with $l_{sq}$
| Algorithm $\boldsymbol{(f,\hat{\eta})}$
---
with $l_{log}$
| Algorithm $\boldsymbol{(f,\hat{\eta})}$
---
with $l_{msq}$
LSPC | k-NN
MSE |  0.026$\pm$0.002 |  0.002$\pm$0.0008 |  0.003$\pm$0.0008 |  0.011$\pm$0.002 |  0.0232 $\pm$ 0.002
RMSE |  0.161$\pm$0.006 |  0.047$\pm$0.009 |  0.005$\pm$0.007 |  0.106$\pm$0.010 |  0.152 $\pm$ 0.007
MAD |  0.138$\pm$0.005 |  0.026$\pm$0.005 |  0.031$\pm$0.004 |  0.077$\pm$0.011 |  0.0911 $\pm$ 0.00756
MD |  0.004$\pm$0.006 |  0.002$\pm$0.006 |  0.003$\pm$0.006 |  0.003$\pm$0.003 |  0.01434 $\pm$ 0.005
KL | NA |  0.010$\pm$0.003 | NA | NA | NA
Acc |  0.912$\pm$0.018 |  0.911$\pm$0.016 |  0.908$\pm$0.016 |  0.913$\pm$0.017 |  0.9195 $\pm$ 0.0181
BS |  0.0795$\pm$0.008 |  0.0645$\pm$0.012 |  0.0636$\pm$0.012 |  0.0733$\pm$0.0102 |  0.0834 $\pm$ 0.011
DiffMax |  0.463$\pm$0.03 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0
DiffMin |  -0.551$\pm$0.024 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0 |  0.0 $\pm$0.0
Table 22: Performance of estimates obtained by Algorithm $\boldsymbol{(f,\hat{\eta})}$ with squared loss, logistic loss and modified squared loss and LSPC and $k$-NN on $Syn_{10D}$ dataset. The cells with NA means that some estimate was 0 and the KL couldn’t be defined. Cells with $0.0$ entry means that the value is less than $1\times 10^{-5}$. | KLIEP based estimator used
---|---
| Performance
---
measures
$\mathbf{\hat{\eta}_{pos}}$ | $\mathbf{\hat{\eta}_{neg}}$ | $\mathbf{\hat{\eta}_{norm}}$
MSE |  0.06$\pm$ 0.003 |  0.062$\pm$0.006 |  0.035$\pm$ 0.005
RMSE |  0.250$\pm$ 0.006 |  0.249$\pm$0.012 |  0.189$\pm$0.014
MAD |  0.208$\pm$ 0.007 |  0.207$\pm$0.013 |  0.177$\pm$0.014
MD |  0.0031$\pm$0.01 |  0.006$\pm$0.012 |  0.008$\pm$0.005
KL | NA | NA |  0.014$\pm$ 0.017
Acc |  0.848 $\pm$ 0.03 |  0.846 $\pm$ 0.017 |
BS |  0.125 $\pm$ 0.010 |  0.123$\pm$ 0.009 |  0.0995$\pm$0.006
DiffMax |  0.32$\pm$0.098 |  -0.02 $\pm$0.006 |  -0.04 $\pm$0.001
DiffMin |  0.015 $\pm$0.005 |  -0.422 $\pm$0.12 |  0.031 $\pm$0.009
Table 23: Performance of $\hat{\eta}_{pos},\hat{\eta}_{neg},\hat{\eta}_{norm}$, KLIEP based estimates on $Syn_{10D}$ dataset. The cells with NA means that some estimate was $0$ and the KL couldn’t be defined. Cells with $0.0$ entry means that the value is less than $1\times 10^{-5}$. | KLIEP based estimator used
---|---
| Performance
---
measures
$\mathbf{\hat{\eta}_{pos,s}}$ | $\mathbf{\hat{\eta}_{neg,s}}$ | $\mathbf{\hat{\eta}_{norm,s}}$ | $\mathbf{\hat{\eta}_{Cnorm,s}}$ | $\mathbf{\hat{\eta}_{avg,s}}$
MSE |  0.035$\pm$0.005 |  0.061$\pm$0.003 |  0.104$\pm$0.014 |  0.036$\pm$0.005 |  0.084$\pm$0.006
RMSE |  0.247$\pm$0.006 |  0.321$\pm$0.022 |  0.191$\pm$0.013 |  0.290$\pm$0.01 |  0.260$\pm$0.006
MAD |  0.245$\pm$0.009 |  0.263$\pm$0.009 |  0.178$\pm$0.014 |  0.275$\pm$0.011 |  0.247$\pm$0.007
MD |  -0.181$\pm$0.03 |  0.157$\pm$0.035 |  0.019$\pm$0.007 |  0.018$\pm$0.012 |  -0.019$\pm$0.010
KL |  0.262 $\pm$0.022 | NA |  0.143 $\pm$ 0.016 | NA |  0.228$\pm$0.009
Acc |  0.750 $\pm$ 0.044 |  0.736 $\pm$ 0.044 |  0.913 $\pm$ 0.016 |  0.913 $\pm$0.016 |  0.913 $\pm$ 0.016
BS |  0.158$\pm$0.0159 |  0.165$\pm$0.0158 |  0.099$\pm$0.006 |  0.149$\pm$0.008 |  0.1318$\pm$0.006
DiffMax |  0.0 $\pm$0.0 |  -0.01 $\pm$0.006 |  -0.040 $\pm$0.012 |  0.0 $\pm$0.0 |  -0.08$\pm$0.018
DiffMin |  0.011 $\pm$0.005 |  0.0 $\pm$0.0 |  0.033 $\pm$0.009 |  0.0 $\pm$0.0 |  0.081$\pm$0.02
Table 24: Performance of
$\hat{\eta}_{pos,s},\hat{\eta}_{neg,s},\hat{\eta}_{norm,s},\hat{\eta}_{Cnorm,s},\hat{\eta}_{avg,s}$,
KLIEP based estimates on $Syn_{10D}$ dataset. The cells with NA means that
some estimate was $0$ and the KL couldn’t be defined. Cells with $0.0$ entry
means that the value is less than $1\times 10^{-5}$.

Observations We observe that as far as label prediction is concerned Algorithm
$\boldsymbol{(f,\hat{\eta})}$ with squared loss, logistic loss and modified
squared loss, LSPC, $k$-NN and KLIEP’s $\hat{\eta}_{norm}$ estimators perform
equally well as their accuracy is the highest in comparison to others. This
trend has been observed across all 3 synthetic datasets.

With respect to DiffMax and DiffMin measures, Algorithm
$\boldsymbol{(f,\hat{\eta})}$ with logistic loss and modified squared loss,
LSPC, $k$-NN and KLIEP’s $\hat{\eta}_{norm}$ estimators are good.

##  Appendix 0.G Risk minimization under weighted 0-1 loss in the extended
cost space

In view of the negative result obtained with the $\alpha$-weighted $0$-$1$
loss function, we explore a larger cost space. In a similar kind of
representation as that of $C_{1},C_{-1}$, one can extend the parameter space
from $(C_{1},C_{-1})$ to $(C_{1},C_{-1},c_{1},c_{-1})$ given as follows.

  * $\bullet$

$C_{1}$ (or $C_{-1}$): loss when a positive (or negative) class point is
misclassified

  * $\bullet$

$c_{1}$ (or $c_{-1}$): loss when a positive (or negative) class point is
correctly classified

Let the extended cost sensitive $0$-$1$ loss be defined as

$l_{c,4}(f(\mathbf{x}),y)=\mathbf{1}_{\\{f(\mathbf{x})>0,y=1\\}}c_{1}+\mathbf{1}_{\\{f(\mathbf{x})>0,y=-1\\}}C_{-1}+\mathbf{1}_{\\{f(\mathbf{x})\leq
0,y=1\\}}C_{1}+\mathbf{1}_{\\{f(\mathbf{x})\leq 0,y=-1\\}}c_{-1}$ (55)

Then, the corresponding risk on $D$ and $\tilde{D}$ be denoted by
$R_{D,l_{c,4}}(f)$ and $R_{\tilde{D},l_{c,4}}(f)$ respectively.

One can reduce the extended $(C_{1},c_{1},C_{-1},c_{1})$ based $4$ cost space
again to $(C_{1}-c_{1},C_{-1}-c_{-1})$ based $2$ cost space as the threshold
of the optimal classifiers does not change by reducing the parameter space.
This has also been shown in [8]. However, we still compute the minimizers of
the risks $R_{D,l_{c,4}}(f)$ and $R_{\tilde{D},l_{c,4}}(f)$ and show that they
need not be equal.

Using the same technique as in $(C_{1},C_{-1})$ cost case, one can show that
the minimizers of $R_{D,l_{c,4}}$ and $R_{\tilde{D},l_{c,4}}$ are as follows:

$f^{*}_{4c}(\mathbf{x})=sign\left(\eta(\mathbf{x})-\frac{C_{-1}-c_{-1}}{C_{1}-c_{1}+C_{-1}-c_{-1}}\right)$
(56)
$\tilde{f}^{*}_{4c}(\mathbf{x})=sign\left(\eta(\mathbf{x})-\frac{C_{-1}-c_{-1}-\rho(C_{1}-c_{1}+C_{-1}-c_{-1})}{(1-2\rho)(C_{1}-c_{1}+C_{-1}-c_{-1})}\right)$
(57)

The threshold of $\eta$ in equation (56) and (57) cannot be same unless
$C_{1}-c_{1}=s=C_{-1}-c_{-1}$ or $\rho=0$. Therefore, even in the extended
cost space, under uniform noise and differential cost,
$f^{*}_{4c}\neq\tilde{f}^{*}_{4c}$, i.e., learning a classifier from $l_{c,4}$
which is both noise robust and cost sensitive is not always possible. This is
independent of the extended sufficient condition (E-SC) that we now introduce;
as we will see, this condition doesn’t help us for joint uniform noise
robustness and cost sensitivity.

A natural extension of sufficient condition (Symmetry condition SC [6]) in
cost sensitive case when used with weighted $0$-$1$ loss can be written as
follows. For an arbitrary classifier $f(\mathbf{x})$ and $K>0$,

$\displaystyle\text{{E-SC}}\hskip 71.13188pt\text{when
}f(\mathbf{x})>0~{}\text{ SC becomes }~{}C_{-1}+c_{1}=K$ (58)
$\displaystyle\text{when }f(\mathbf{x})\leq 0~{}\text{ SC becomes
}~{}C_{1}+c_{-1}=K$

Clearly, $C_{-1}+c_{1}=K=C_{1}+c_{-1}$ which can be equivalently written as
$C_{-1}-c_{-1}=s=C_{1}-c_{1}$. This implies that the effective cost for each
class is same, i.e. $s$ and hence no differential costs for the classes.
Therefore, E-SC cannot be sufficient condition for SLN robustness if there is
differential costing of $(C_{1},c_{1},C_{-1},c_{-1}).$


# Optically Induced Avoided Crossing in Graphene

Sören Buchenau<EMAIL_ADDRESS>Benjamin Grimm-Lebsanft
Florian Biebl  Tomke Glier  Lea Westphal  Institute of Nanostructure and Solid
State Physics, University of Hamburg, 22761, Germany  Janika Reichstetter
Dirk Manske  Max Planck Institute for Solid State Research, 70569 Stuttgart,
Germany  Michael Fechner  Andrea Cavalleri  Max Planck Institute for the
Structure and Dynamics of Matter, 22761 Hamburg, Germany  Sonja Herres-Pawlis
Institute of Inorganic Chemistry, RWTH Aachen University, 52074 Aachen,
Germany  Michael Rübhausen<EMAIL_ADDRESS>Institute of
Nanostructure and Solid State Physics, University of Hamburg, 22761, Germany

###### Abstract

Degenerate states in condensed matter are frequently the cause of unwanted
fluctuations, which prevent the formation of ordered phases and reduce their
functionalities. Removing these degeneracies has been a common theme in
materials design, pursued for example by strain engineering at interfaces.
Here, we explore a non-equilibrium approach to lift degeneracies in solids. We
show that coherent driving of the crystal lattice in bi- and multilayer
graphene, boosts the coupling between two doubly-degenerate modes of E1u and
E2g symmetry, which are virtually uncoupled at equilibrium. New vibronic
states result from anharmonic driving of the E1u mode to large amplitdues,
boosting its coupling to the E2g mode. The vibrational structure of the driven
state is probed with time-resolved Raman scattering, which reveals laser-field
dependent mode splitting and enhanced lifetimes. We expect this phenomenon to
be generally observable in many materials systems, affecting the non-
equilibrium emergent phases in matter.

Phononics, Raman, Pump-Probe, Graphene

##  I Introduction

In condensed matter systems and at interfaces, static symmetry breaking has
been used to create new states of matter and to introduce new functionalities
in equilibrium [1, 2, 3]. However, using coherent and dynamic optical control
offers an innovative approach to expand this concept beyond the limitations of
static control [4, 5] and on short or even ultrashort time scales. There are
many examples where coherent dynamical control of matter by photons has
triggered phase transitions, tuned interactions, and generated new forms of
matter [6]. Even though some of these states exist only for picoseconds, it is
possible to investigate them using pulsed lasers and pulsed currents in
transport measurement [7, 8]. The structural degrees of freedom, represented
by phonon modes, are essential for the thermodynamic properties in crystals.
Thus, materials with degenerate phonons can exhibit strong anharmonic
couplings leading to new vibronic states resulting in new forms of matter.

Graphene systems have been previously excited mostly by high-photon-energy
excitations, in resonance with electronic transitions over a broader energy
range of 1-2 eV [9, 10]. Here, however, we drive structural degrees of freedom
in resonance with a low-energy degenerate IR phonon of bilayer graphene. The
E1u IR phonon and E2g Raman active phonon are uncoupled in equilibrium. We
resonantly drive the infrared vibration at 196 meV ($\sim$ 48 THZ) using a
mid-IR pump. The vibronic state of the driven system is then monitored via the
corresponding transient Raman response. We explore the physics of lifting
degeneracies by optical control in bi- and multilayer graphene and demonstrate
the emergence of an optically induced avoided crossing that can be tuned by
the fluence of the pump laser.

Monolayer graphene exhibits the prominent G-phonon [11, 12], which is only
Raman-active. On the other hand, from two layers and up, the material exhibits
fully degenerate IR- and Raman-active phonons with orthogonal E1u and E2g
symmetry [12, 13, 14]. These modes are ideal candidates for strong anharmonic
coupling due to their strongly related eigenvectors and degenerate energies.
Figure 1 (a) shows the IR- and Raman-active oscillators on each side. Without
coupling in equlibrium, the phonons are at the same eigenfrequency. Figure 1
(b), driving the IR mode leads to a lifting of the degeneracy via anharmonic
coupling of the two phonon modes [15], which is indicated in the figure
symbolically as an additional spring coupling the IR and Raman active phonons.
This optically induced anharmonic coupling out of equilibrium lifts the
degeneracy resulting in two new vibrational states $\beta^{+}$ and
$\beta^{-}$, representing two coupled harmonic oscillators. The vibration
$\beta^{-}$ is antisymmetric, i.e., the masses move in opposite directions,
whereas the vibration $\beta^{+}$ is symmetric, i.e. the masses move in the
same direction. From those, the IR- and Raman-active modes can be obtained by
$\beta^{+}+\beta^{-}$ and $\beta^{+}-\beta^{-}$, respectively. Hence the
transient Raman signal vanishes in the degenerate state, i.e.
$\beta^{+}=\beta^{-}$ and is, therefore, a precise tool to study the
anharmonic coupling out of equilibrium.

Figure 1: Optically induced avoided crossing shown with classical springs and
masses. (a) The IR oscillator XIR and the Raman active oscillator XR are
degenerate at $\omega_{0}$, but exhibit different symmetries, i.e., odd and
even. Therefore, without any coupling, an oscillation of either mode does not
impact the other one. (b) When the IR mode is strongly driven out of the
harmonic approximation, the anharmonic coupling $\alpha$ couples the
oscillators and the vibronic systems form the two novel vibrational states
$\beta^{-}$ and $\beta^{+}$. The vibration $\beta^{-}$ is antisymmetric, i.e.
the masses move in opposite directions, whereas the vibration $\beta^{+}$ is
symmetric, i.e. the masses move in the same direction. The introduced energy
difference between the two states, $\Delta$, is given by the fluence of the
pump laser. The Raman-active superposition $\beta^{+}-\beta^{-}$ can be
experimentally measured and exhibits a subharmonic frequency, tunable by the
field-dependent splitting $\Delta$.

If one considers a relative shift of the respective eigenvalues in the limit
of small Rabi frequencies, the energy-shifted values $\beta^{-}$ and
$\beta^{+}$ can be expressed by
$\omega_{\pm}\approx\omega_{0}\left(1\pm\Delta/2\omega_{0}^{2}\right)$, with
$\omega_{0}$ as the degenerate frequency of the phonons in the equilibrium
state and $\Delta$ the representing lifting of the degeneracy due to the
anharmonic coupling out of equilibrium. Thus, the orthogonal harmonic motions
of the IR- and Raman-active vibrations are
$\sim\text{e}^{i\omega_{0}t}\sin{\left(\Delta\cdot t/2\omega_{0}\right)}$ and
$\sim\text{e}^{i\omega_{0}t}\cos{\left(\Delta\cdot t/2\omega_{0}\right)}$,
respectively. Accordingly, the IR and Raman coordinates vibrate with the
original degenerate eigenfrequency $\omega_{0}$, modulated by a subharmonic
frequency and are phase-shifted by $\pi/2$. In static quantum systems, the
avoided crossing [16, 17] phenomenon leads to lower eigenvalues, which
increases the stability of bonded states in molecular systems. In optically
driven non-equilibrium systems, the magnitude of degeneracy lifting $\Delta$
depends on the strength of the electric field as it tunes the anharmonic
coupling. As a result, the populations of the lower and higher energy states
undergo time-dependent oscillations. The transient spontaneous Raman
scattering of a Raman-active phonon, driven by an IR-active phonon, is the
ideal method to observe these states. However, this observation is only
possible in case of a symmetry breaking yielding a finite response due to an
optically driven avoided crossing, i.e., when $\beta^{+}\neq\beta^{-}$.

##  II Experiment

The general outline of the experiment is shown in Figure 2 (a). A tunable mid-
IR pulse drives the E1u IR phonon mode in resonance at 196 meV. At different
time delays, transient spontaneous Raman scattering probes the Raman-active
E2g mode. In bi- and multilayer systems, the presence of the polar IR-active
phonon results in a strong transient response in the excited state, as seen in
Figure 2 (b). The transients are obtained by subtracting the unpumped signal
from the pumped signal, isolating the response of the IR-driven Raman mode.
Note that monolayer graphene lacks the polar vibrational mode [11] and,
therefore, does not exhibit a transient response (see SI Appendix, Figure S2).
In multilayer graphene, the spontaneous Raman signal in the pumped state is
enhanced on the Stokes and Anti-Stokes side.

Figure 2: Outline of the transient spontaneous Anti-Stokes and Stokes Raman
scattering experiment. (a) A mid-IR pump pulse excites the IR vibration in the
sample, which is probed by a Raman pulse. The delay, pump energy, and fluence
are controlled. The insets show the movements of the carbon atoms during the
vibrations. (b) Resonant mid-IR pumping causes a strong enhancement of the
signal in bi- and multilayer graphene. (c) Selective driving of the E1u
vibration by tuning the pump energy through the sharp resonance. Note that
small visible features (see asterisk) might be due to sum-frequency
scattering. Horizontal grey dashed lines indicate the zero line for each
spectrum and vertical ones the energy of the unpumped IR phonon.

In the following, we show the transient spontaneous Raman data on the Anti-
Strokes side with the expected stronger transient as the pump excites the
coupled vibronic system out of equilibrium. The transient Anti-Stokes Raman
response as a function of the energy of the mid-IR pump is shown in Figure 2
(c). The first apparent observation is the ultra-sharp resonance when tuning
the mid-IR laser through the energy of the polar phonon. We observe a peaking
at 196 meV with a width of about 6 meV. This observation, together with the
absence of a transient Raman response in monolayer graphene (see SI Appendix,
Figure S2) and the absence of any transients above and below the energy of the
polar IR phonon (see SI Appendix, Figure S4 (a) and Figure 2 (c), highlights
the direct coupling of the transient Raman E2g phonon to the IR-active E1u
phonon. Note that we do not observe any electronic background, as might be
expected when pumping at optical energies (see SI Appendix, Figure S4 (b)).
Additionally, we observe changes in the shape of the transient Anti-Stokes
Raman response as the energy of the mid-IR pump is changed and slightly
detuned from the energy of the E1u phonon. Thus, the data show a purely
phononic mechanism for driving graphene systems.

The incident mid-IR pump photon drives the E1u phonon, which alters the
response of the Raman signal of the E2g phonon. The degeneracy of the two
phonons leads to strong anharmonic coupling [15, 18]. Without the pump
interaction, the Raman response is given by the Raman susceptibility, yielding
a typical Lorentzian lineshape. The anharmonic coupling $\alpha$ between a
Raman and IR-active phonon requires a quadratic coupling to the Raman-active
phonon in the equation of motion as a consequence of a cubic contribution to
the potential, i.e. the anharmonic coupling. This leads to a set of equations
of motions describing the electric field-driven IR-active phonon
$X_{\text{IR}}(t)$ and Raman-active phonon $X_{\text{R}}(t)$

$\displaystyle\widetilde{D}(\omega_{0},\gamma_{0})X_{\text{IR}}(t)$
$\displaystyle=\alpha X_{\text{R}}(t)X_{\text{IR}}(t)+F_{\text{L}}(t)$ (1a)
$\displaystyle\widetilde{D}(\omega_{0},\gamma_{0})X_{\text{R}}(t)$
$\displaystyle=\alpha X^{2}_{\text{IR}}(t).$ (1b)

Here,
$\widetilde{D}\left(\omega_{0},\gamma_{0}\right)=\partial_{t}^{2}+\gamma_{0}\partial_{t}+\omega_{0}^{2}$
is the operator describing a damped harmonic oscillator and
$F_{\mathrm{L}}\left(t\right)$ represents the Gaussian laser beam with
frequency $\omega_{\mathrm{L}}$ and width $\sigma$. Assuming that IR- and
Raman-active phonons are degenerate, the above set of coupled nonlinear
differential equations can be solved analytically. For this procedure, one
adds and subtracts both equations from each other and introduces common
coordinates, linking the coordinates of the IR- and Raman-active phonons by
$\beta^{+}=X_{\mathrm{IR}}+X_{\mathrm{R}}$ and
$\beta^{-}=X_{\mathrm{IR}}-X_{\mathrm{R}}$. The inverse Fourier transformation
of $\beta^{\pm}$ then yields the expressions for these effective phonon
propagators after applying the convolution theorem for the quadratic parts of
the equation

$\displaystyle\beta^{+}(\omega)$
$\displaystyle=\frac{\alpha\left[(\beta^{-}*\beta^{+})(\omega)\right]+2F_{\text{L}}(\omega)}{2(\omega_{-}^{2}-\omega^{2}+i\gamma_{-}\omega)}$
(2a) $\displaystyle\beta^{-}(\omega)$
$\displaystyle=\frac{-\alpha\left[(\beta^{+}*\beta^{-})(\omega)\right]+2F_{\text{L}}(\omega)}{2(\omega_{+}^{2}-\omega^{2}+i\gamma_{+}\omega)},$
(2b)

with the shifted eigenfrequencies $\omega_{\pm}^{2}=\omega_{0}^{2}\pm\Delta$
and modified damping constants $\gamma_{\pm}=\gamma_{0}\pm\delta/\omega$.
These changes scale with the anharmonic coupling constant $\alpha$ and the
square of the electric field $E_{0}^{2}$. The modifications $\Delta$ and
$\delta$ represent the frequency-dependent contributions from the real and
imaginary parts of the convolution integrals. The contribution from the non-
mixed convolution integrals
$\left(\beta^{+}\ast\beta^{+}\right)\left(\omega\right)$ and
$\left(\beta^{-}\ast\beta^{-}\right)\left(\omega\right)$ vanishes and they are
omitted in eq. (2a) and (2b). The mixed terms
$\left(\beta^{-}\ast\beta^{+}\right)\left(\omega\right)$ have to be evaluated
self-consistently. In the following, we assume $\Delta$ and $\delta$ to be
frequency independent parameters to model the Raman response. However, both
parameters depend on time and the fluence of the pump $F_{\mathrm{L}}$. The
leading order contribution results in the modified propagators with the
driving pump laser in the numerators, shown in eq. (2a) and (2b). The Green’s
function of the Raman-active vibration can now be calculated by
$X_{\mathrm{R}}\left(\omega\right)=\frac{1}{2}\left(\beta^{+}-\beta^{-}\right)$
and the Raman scattering intensity is given by
$-\mathrm{Im}\left[X_{\mathrm{R}}\left(\omega\right)\right]$

$\displaystyle I($
$\displaystyle\omega)=\frac{F_{\text{L}}(\omega)\Delta\omega(\gamma_{+}(\omega_{-}^{2}\\!-\\!\omega^{2})\\!+\\!\gamma_{-}(\omega_{+}^{2}\\!-\\!\omega^{2}))}{\mathcal{D}}+$
$\displaystyle\frac{F_{\text{L}}(\omega)\delta(\gamma_{-}\gamma_{+}\omega^{2}\\!-\\!(\omega_{-}^{2}\\!-\\!\omega^{2})(\omega_{+}^{2}\\!-\\!\omega^{2}))}{\mathcal{D}}$
(3a)
$\displaystyle\begin{split}\mathcal{D}=((\omega^{2}_{+}-\omega^{2})(\omega^{2}_{-}-\omega^{2})-\gamma_{-}\gamma_{+}\omega^{2})^{2}\\\\[6.45831pt]
+\omega^{2}(\gamma_{+}(\omega_{-}^{2}-\omega^{2})+\gamma_{-}(\omega_{+}^{2}-\omega^{2}))^{2}\end{split}$
(3b)
$\displaystyle\begin{split}F_{\text{L}}(\omega)\\!=\\!\frac{f_{\text{L}}}{\sqrt{2\pi}}(e^{-\frac{1}{2}\sigma^{2}(\omega-\omega_{\text{L}})^{2}}\\!+\\!e^{-\frac{1}{2}\sigma^{2}(\omega+\omega_{\text{L}})^{2}}).\end{split}$
(3c)

Figure 3: Experimental and calculated transients for mid-IR resonant driving.
(a) Exemplary transient spectra around the resonance energy of the E1u mode
for multilayer graphene. All non-trivial transients are shown in Figure 2 (c).
Note that small visible features (see asterisk) might be due to sum-frequency
scattering. (b) Simulation results obtained by eq. (3) for the values of
$\Delta=4$ meV2 and $\delta=0.12$ meV2. (c) Intensity scaling of the transient
as a function of pump energy. (d) Frequency tuning of the transient as a
function of pump energy. (e) Change of the width as a function of pump energy.
In (c) to (e) black squares represent experimental results and blue circles
the calculation. Horizontal grey dashed lines indicate the frequency and width
of the unpumped Raman phonon and vertical ones the energy of the IR phonon.
Solid grey and blue lines are guides to the eye for each parameter.

As expected, the transient Raman response vanishes for
$\omega_{+}=\omega_{-}=\omega_{0}$, i.e. $\Delta\rightarrow 0$ and
$\gamma_{+}=\gamma_{-}=\gamma_{0}$, i.e., $\delta\rightarrow 0$, representing
the case of two non-interacting oscillators at identical frequencies and
$\alpha=0$. The response from the above expression exhibits a classical anti-
resonance contribution in the first term and a sharpened Lorentzian
contribution from the second term. It is, therefore, the unusual case of a
Fano-like response resulting from the interaction of two discrete states,
which is quite different from the typical case of a Fano lineshape where a
discrete excitation interacts with a continuum of states as in electron-phonon
coupled systems [19, 20, 21, 22]. It is also remarkable that the strength of
the change in frequency and width, i.e., the changes in real and imaginary
parts of the phonon susceptibility, directly influences the lineshape and the
strength of the response. Both are proportional to the square of the electric
field, which drives the IR phonon. Stronger changes in the frequency lead to a
more asymmetric Fano-like lineshape, and larger changes in the width lead to a
more symmetric Lorentzian response. We can now model the observed resonance
interaction between IR- and Raman-active phonons, as shown in Figure 3. The
frequencies and widths $\omega_{0}$, $\gamma_{0}$, $\omega_{L}$, and $\sigma$
are known. The only two model parameters are the modifications to the
frequency $\Delta$ and the width $\delta$ as outlined above. The lineshape
obtained is only dependent on the ratio of the parameters $\Delta$ and
$\delta$, so these values are subject to an unknown scaling parameter that is
given by the enhancement of the response. The simulation results are shown in
Figure 3 (a) and (b) for the experimental and calculated transients. The data
are reproduced very well, considering that the frequency-dependent quantities
$\Delta$ and $\delta$, are modeled by constants. The signal enhancement, the
shifting of mode frequencies, and the sharpening of the transient response are
recovered by the simulation where the ratio of $\Delta/\delta$ was kept
constant and set to $\Delta=4$ meV2 and $\delta=0.12$ meV2. The two anti-
resonances in the transient in Figure 3 (b) are not observed in the
experiment. We attribute these differences to higher-order coupling terms that
we have neglected in our calculation and to the fact that $\Delta$ and
$\delta$ should exhibit a frequency dependence. The shape of the resonance
tuning of the intensity and the absolute frequency shifts are well recovered,
as can be seen in Figure 3 (c) and (d). Even the effect of the sharpening as
such is present in our calculation in Figure 3 (e). However, the details of
the broadening when tuning out of the resonance are not captured. The FWHM of
the Raman mode measured in the probe spectra is 3 meV, the width of the
observed transient in resonance is 2 meV, and the width of the calculated
transient is 1.4 meV. The linewidth variation with the pump energy is well
reproduced but takes place on a smaller scale than the experiment.

Additionally, assuming a time-dependent dephasing between the IR- and Raman-
active modes, the ratio between $\Delta$ and $\delta$ in eq. (3) will increase
with time. The calculation recovers this delay dependence between pump and
probe pulses, as shown in the SI Appendix Figure S3, demonstrating a transient
state with a lifetime of several picoseconds. This lifetime is remarkable
since, in the static limit, the widths of the phonons correspond to a lifetime
of about 200 fs. The observation of the transient Raman response on the time
scale of several picoseconds implies an increased robustness of the non-
equilibrium coupled driven state compared to the equilibrium. This is further
supported by the fact that the transient Raman response is substantially
sharper than its static spontaneous Raman response.

Figure 4: Fluence-dependence of the lifting of the degeneracy. (a) Fluence-
dependent transient signals. (b) Calculated transient signals obtained by
equation (3a). (c) DFT calculation of the splitting of the frequencies as a
function of the fluence of the driving IR field. (d) Amplitude of the
$\beta^{+}$ and $\beta^{-}$ oscillations at 0.60 mJcm-2. The data are derived
from the simulation and by comparison to the DFT calculated fluence-dependent
splitting of the frequencies. (e) Total energy from the mixed vibrational
state for 0.24 and 0.60 mJcm-2. The magnitude tsubh scales with the lifting of
the degeneracy $\Delta$. (f) Schematic of the field-dependent splitting of the
oscillations.

Increasing the fluence at constant mid-IR pump frequency and delay leads to a
linearly increasing transient signal, as shown in Figure 4 (a), which is well
recovered by the simulation in Figure 4 (b). The symmetry breaking of the
degenerate energy levels is directly proportional to the fluence of the mid-IR
pulse, and so is the transient Raman signal. Figure 4 (c) shows a DFT
calculation of the magnitude of the energy splitting of the common coordinates
$\beta^{+}$ and $\beta^{-}$ as a function of the field strength of the driving
mid-IR field. The DFT calculation allows now to determine real values of
$\Delta$ and $\delta$, and to calculate the effective vibrational patterns as
shown in Figure 4 (d), clearly resembling a novel non-equilibrium system of
two coupled quantum states. Figure 4 (e) shows the total energy of the 0.24
mJcm-2 and 0.60 mJcm-2 fluence data. It shows the period of the subharmonic
Rabi oscillation tsubh, which depends on the fluence. Accordingly, it is
inversely proportional to the fluence-dependent splitting of the avoided
crossing $\Delta$, as shown in Figure 4 (f). The splitting separates the
antisymmetric and symmetric vibrations $\beta^{+}$ and $\beta^{-}$ as shown in
Figure 1 (b).

In summary, we have shown that degenerate phonons in bi- and multilayer
graphene can be used to generate novel vibronic states due to an induced
anharmonic coupling, tuned by driving a polar phonon mode. Coherent control
over the vibronic states can be obtained by tuning the electric field of the
IR laser driving the phonon. The lower and higher energy states correspond to
novel symmetric and antisymmetric vibronic states, resembling the physics of
the avoided crossing. These states exchange population density at a
subharmonic Rabi frequency, defined by the applied electric field, effectively
leading to a breaking of continuous time-translation symmetry. The theoretical
as well as experimental framework can be extended towards other 2D materials
with degenerate or nearly degenerate phonon modes evolving the field of
nonlinear phononics.

###### Acknowledgements.

We acknowledge funding by the Deutsche Forschungsgemeinschaft via RU 773/8-1
and the Bundesministerium für Bildung und Forschung via 05K19GU5 and 05K22GU2.

## References

  * Von Boehm and Bak [1979] J. Von Boehm and P. Bak, Devil’s stairs and the commensurate-commensurate transitions in cesb, Physical Review Letters 42, 122 (1979).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Colloquium: topological insulators, Reviews of modern physics 82, 3045 (2010).
  * Dong  _et al._ [2015] S. Dong, J.-M. Liu, S.-W. Cheong, and Z. Ren, Multiferroic materials and magnetoelectric physics: symmetry, entanglement, excitation, and topology, Advances in Physics 64, 519 (2015).
  * Juraschek  _et al._ [2017] D. M. Juraschek, M. Fechner, A. V. Balatsky, and N. A. Spaldin, Dynamical multiferroicity, Physical Review Materials 1, 014401 (2017).
  * McIver  _et al._ [2020] J. W. McIver, B. Schulte, F.-U. Stein, T. Matsuyama, G. Jotzu, G. Meier, and A. Cavalleri, Light-induced anomalous hall effect in graphene, Nature physics 16, 38 (2020).
  * Bloch  _et al._ [2022] J. Bloch, A. Cavalleri, V. Galitski, M. Hafezi, and A. Rubio, Strongly correlated electron–photon systems, Nature 606, 41 (2022).
  * Wessels  _et al._ [2016] P. Wessels, A. Vogel, J.-N. Tödt, M. Wieland, G. Meier, and M. Drescher, Direct observation of isolated damon-eshbach and backward volume spin-wave packets in ferromagnetic microstripes, Scientific reports 6, 22117 (2016).
  * Disa  _et al._ [2020] A. S. Disa, M. Fechner, T. F. Nova, B. Liu, M. Först, D. Prabhakaran, P. G. Radaelli, and A. Cavalleri, Polarizing an antiferromagnet by optical engineering of the crystal field, Nature Physics 16, 937 (2020).
  * Yan  _et al._ [2009] H. Yan, D. Song, K. F. Mak, I. Chatzakis, J. Maultzsch, and T. F. Heinz, Time-resolved raman spectroscopy of optical phonons in graphite: Phonon anharmonic coupling and anomalous stiffening, Physical Review B 80, 121403 (2009).
  * Kang  _et al._ [2010] K. Kang, D. Abdula, D. G. Cahill, and M. Shim, Lifetimes of optical phonons in graphene and graphite by time-resolved incoherent anti-stokes raman scattering, Physical Review B 81, 165405 (2010).
  * Malard  _et al._ [2009] L. Malard, M. A. Pimenta, G. Dresselhaus, and M. Dresselhaus, Raman spectroscopy in graphene, Physics reports 473, 51 (2009).
  * Nemanich  _et al._ [1977] R. Nemanich, G. Lucovsky, and S. Solin, Infrared active optical vibrations of graphite, Solid State Communications 23, 117 (1977).
  * Ferrari and Basko [2013] A. C. Ferrari and D. M. Basko, Raman spectroscopy as a versatile tool for studying the properties of graphene, Nature nanotechnology 8, 235 (2013).
  * Davydov [1964] A. S. Davydov, The theory of molecular excitons, Soviet Physics Uspekhi 7, 145 (1964).
  * Först  _et al._ [2011] M. Först, C. Manzoni, S. Kaiser, Y. Tomioka, Y.-n. Tokura, R. Merlin, and A. Cavalleri, Nonlinear phononics as an ultrafast route to lattice control, Nature Physics 7, 854 (2011).
  * Landau [1932] L. D. Landau, A theory of energy transfer ii, Phys. Z. Sowjetunion 2, 19 (1932).
  * Zener [1932] C. Zener, Non-adiabatic crossing of energy levels, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 137, 696 (1932).
  * Nicoletti and Cavalleri [2016] D. Nicoletti and A. Cavalleri, Nonlinear light–matter interaction at terahertz frequencies, Advances in Optics and Photonics 8, 401 (2016).
  * Fano [1961] U. Fano, Effects of configuration interaction on intensities and phase shifts, Physical Review 124, 1866 (1961).
  * Hasdeo  _et al._ [2014] E. H. Hasdeo, A. R. Nugraha, M. S. Dresselhaus, and R. Saito, Breit-wigner-fano line shapes in raman spectra of graphene, Physical Review B 90, 245140 (2014).
  * Zhang  _et al._ [2015] W.-L. Zhang, H. Li, D. Xia, H. Liu, Y.-G. Shi, J. Luo, J. Hu, P. Richard, and H. Ding, Observation of a raman-active phonon with fano line shape in the quasi-one-dimensional superconductor k 2 cr 3 as 3, Physical Review B 92, 060502 (2015).
  * Bock  _et al._ [1999] A. Bock, S. Ostertun, R. D. Sharma, M. Rübhausen, K.-O. Subke, and C. Rieck, Anomalous self-energy effects of the b 1 g phonon in y 1- x (pr, ca) x ba 2 cu 3 o 7 films, Physical Review B 60, 3532 (1999).