33institutetext: Cyber Security Research Group, WMG, University of Warwick
Coventry, United Kingdom
33email<EMAIL_ADDRESS>22institutetext:
<EMAIL_ADDRESS>33institutetext<EMAIL_ADDRESS>
# Leveraging Semantic Relationships to Prioritise Indicators of Compromise in
Additive Manufacturing Systems

Mahender Kumar 11 Gregory Epiphaniou 22 Carsten Maple 33

###### Abstract

Additive manufacturing (AM) offers numerous benefits, such as manufacturing
complex and customised designs quickly and cost-effectively, reducing material
waste, and enabling on-demand production. However, several security challenges
are associated with AM, making it increasingly attractive to attackers ranging
from individual hackers to organised criminal gangs and nation-state actors.
This paper addresses the cyber risk in AM to attackers by proposing a novel
semantic-based threat prioritisation system for identifying, extracting and
ranking indicators of compromise (IOC). The system leverages the heterogeneous
information networks (HINs) that automatically extract high-level IOCs from
multi-source threat text and identifies semantic relations among the IOCs. It
models IOCs with a HIN comprising different meta-paths and meta-graphs to
depict semantic relations among diverse IOCs. We introduce a domain-specific
recogniser that identifies IOCs in three domains: organisation_specific,
regional_source-specific, and regional_target-specific. A threat assessment
uses similarity measures based on meta-paths and meta-graphs to assess
semantic relations among IOCs. It prioritises IOCs by measuring their severity
based on the frequency of attacks, IOC lifetime, and exploited vulnerabilities
in each domain.

###### Keywords:

Indicators of Compromise Cyber-Physical Systems Threat Intelligence, Threat
Prioritisation Heterogeneous Information Networks

##  1 Introduction

Industry 4.0, the fourth industrial revolution, refers to integrating advanced
digital technologies and manufacturing systems to automate and optimize
industrial processes. Additive manufacturing (AM) is a key enabler of Industry
4.0, as it allows for the rapid and flexible production of customized parts
and products [1]. AM is a process that enables the production of complex
devices by applying successive layers of materials. AM offers many advantages,
such as on-demand customisation, enhanced logistics, reduced labour and
production lead times, streamlined production, reduced waste, reduced
inventory, and reduced transportation costs. However, cyber and physical
attacks in AM pose severe concerns and formidable challenges [2], making AM
supply chains susceptible to various attack vectors. As a result, protecting
the security of AM has become increasingly important, and developing robust
security mechanisms that protect against a range of potential attacks has
become a significant challenge for researchers and industry practitioners
alike.

Modern attacks on AM are often sophisticated and can exploit hidden
vulnerabilities that go undetected for long periods [3-7]. A prime example is
Advanced Persistent Threats (APTs), which have been used to target AM
industries for espionage, economic gain, and intellectual property theft. APTs
are commonly described as an extended attack campaign in which one or more
intruders execute a long-term plan to take control of a network or system. In
2020, more than 1000 data breaches were reported in the United States alone,
affecting more than 155.8 million individuals through data exposure [3].
Perhaps the most famous kinetic cyber attack of all time was aimed at Iran’s
nuclear program, considered unprecedented in the industry [4]. The Stuxnet
attack involved a complexly targeted worm that executed zero-day exploits on
operating systems and software for managing programmable logic controllers
(PLCs). The attack resulted in tens of billions of dollars in damage. Another
famous example of a cyberattack is the sewage attack in Maroochy Shire, which
caused a system failure and millions of litres of untreated sewage to leak
into the water supply [5]. Belikovetsky et al. [6] conducted a study on the
vulnerability of additive manufacturing to cyber attacks. He demonstrated the
sabotage attack on a propeller blueprint that can be 3D printed at home. The
findings of their study emphasized the vulnerability of additive constructs to
cyber attacks, which was also confirmed by another recent paper [7]. The
authors of the latter paper identified AM as the second most vulnerable
industry to cyberattacks, second only to the financial sector. With AM’s
growing national and industrial importance, cyberattacks have become more
attractive and induced threat actors are increasingly involved in cybercrime-
as-a-service, commoditising cyberattacks. As a result, APTs are now employing
common attack patterns to compromise targets. Therefore, early identification
of threat exposure and breaches is critical to preventing significant damage
to an organization and providing reliable evidence during prosecution trials.

Cyber Threat Intelligence (CTI) can be a valuable tool for assessing the
threat landscape in AM and developing effective strategies for mitigating
cyber risks. Threat intelligence feeds can help organizations stay informed
about emerging threats and new attack techniques. This can be especially
useful in the fast-paced world of AM, where new technologies and processes are
constantly being developed. CTI involves the extraction of threat intelligence
from threat-related information from multiple sources, utilising several
attributes, including Indicators of Compromise (IOCs), Tactic, Technique, and
Procedure (TTP) and the skills and motive of the threat actor [8]. Some
example of CTI feeds are Structured Threat Information Expression (STIX),
OpenDef, Cybox, and OpenIOC, but a massive amount of information remain
unstructured. IBM X-force [9], Facebook ThreatExchange [10], OpenCTI [11] and
MISP [12] are a few vendors who provide threat intelligence feeds by
extracting threat intelligence from multiple open sources using IOC extraction
methods such as PhishTank, and IOCFinder.

These structured threat information feeds have several disadvantages,
including a limited scope, delayed information, high cost, inflexibility and
false positives, making it very challenging for AM industry to rely on them.
On the other hand, unstructured threat feeds can provide a more comprehensive
and flexible approach to threat intelligence that can be more effective for
many organizations. However, unstructured reports may not be well-organized,
making it hard to identify the relationships between different pieces of
information. Other challenges include errors, inaccuracies, and missing
information. As a result, it requires advanced natural language processing
(NLP) techniques and machine learning algorithms to extract meaningful and
relevant threat information from unstructured reports.

Figure 1: An annotated example of CTI includes IOCs such as attack actor,
vulnerability, time, region, file, attack type, device and platform, and their
relationship.

Consider the following instance of a security-related post: ”In October 2019,
Mr Xing from China exploited the CVE-2019-0114 vulnerability, which affected
multiple Android and Win10 devices in the United States. CVE-2019-0114 is a
remote code execution vulnerability that contains the malicious file abc.bat”.
Figure 1 displays a graphical representation of CTI, including eight IOCs such
as attack actor, vulnerability, time, region, file, attack type, device and
platform, and the relationship between them. Existing methods only consider
IOCs but avoid the relationship between them, and as a result, they cannot
grasp a comprehensive picture of the threat landscape. To overcome the
limitations of existing structured threat feed tools, this paper aims to
automate IOC extraction by exploiting Heterogeneous Information Network (HIN)
that provides insight into the interdependencies between heterogeneous IOCs.

This paper presents a novel semantic-based threat prioritisation framework for
identifying, extracting and ranking IOCs. The summary of the paper is as
follows:

  * •

Recogniser. We propose a recogniser that automatically extracts threat-related
information from multi-source threat text. It also identifies the domains to
which IOCs belong and integrates IOCs with their domain, forming three domain-
specific IOCs such as organisation_domain-specific, regional_source-specific,
and regional_target-specific threat intelligence.

  * •

Threat modelling. We model the range of IOCs with a Heterogeneous Information
Network (HIN), which comprises different meta-paths and meta-graphs that
depicts the semantic relations among diverse IOCs to capture a more
comprehensive landscape of threat events.

  * •

Threat assessment. We present a CTI assessment framework that uses similarity
measures based on meta-paths and meta-graphs and assesses the interdependent
relations among diverse IOCs.

  * •

Prioritisation. We then measure the severity of IOCs by considering the
frequency of attacks, IOC lifetime, and the number of exploited
vulnerabilities in each domain. As a result, they evaluate the ranking
mechanism for each IOC.

The rest of the paper is organized as follows: Section 2 discusses the related
work, and Section 3 provides the conceptual background. The proposed framework
is presented in Section 4. Finally, Section 5 summarizes the paper and
provides directions for future research.

##  2 Related Work

Extracting threat intelligence from the unstructured text of threat-related
information has become an exciting research topic in cyber security. This
section briefly describes key methodologies for identifying cyber threats by
extracting IOCs from multiple sources.

Noor et al. [13] have proposed a model to automate cyber threat attribution by
considering high-level IOCs to determine threat actors. Their technique
extracts high-level IOCs from unstructured CTI reports and then semantically
profiles threat actors with the high-level IOCs taken from MITRE’s ATT&CK.
Zhao et al. [14] present TIMiner, a method to extract and assess domain-
specific CTIs that automatically classify the domains associated with related
CTIs. Gao et al. [15] proposed a cyber threat intelligence method based on the
Hierarchical Information Network (HINCTI) system to identify threat types.
HINCTI provides a threat intelligence-based meta-schema to capture the
semantic relationship between threat nodes leveraging a meta-path and meta-
graph-based similarity method. Zhao et al. [16] proposed a CTI framework,
HINTI, based on HIN that proposed a multi-granular-based IOC recogniser to
enhance the accuracy of the IOC extraction method. HINTI defines different
types of IOCs using meta-paths to identify the relationship between IOCs,
profile the threat events and rank the significance of IOCs to understand the
threat landscape.

Liao et al. [17] proposed a novel automated threat-related information
collection and extraction (iACE) technique from unstructured text that uses a
natural language processing method to extract IOC data from text documents and
then analyse the IOC data using graph mining methods. iACE aims to identify
the grammatical semantic relationships between token threat patterns
associated with IOC in text documents. The method integrates name entity
recognition and relation extraction methods. Gao et al. [18] proposed a
threat-hunting system (THREATRAPTOR) that extracts threat behavioural
information from unstructured CTI reports to facilitate threat hunting.
THREATRAPTOR provides an accurate NLP-based auditing framework to extract
structured threat information from unstructured CTI text and defines a domain-
specific query language to detect malicious system activities. Wang et al.
[19] develop an efficient automated process that recognises and extracts
entities and their relationship from text reports.

##  3 Conceptual Background

###  3.1 Cyber Threat Intelligence

Modern cybercriminals have developed sophisticated tactics, techniques, and
procedures (TTP) to realise their aim of compromising their targets quickly
and efficiently. Thus, traditional defence mechanisms, such as anti-virus
software, firewalls and intrusion detection methods, struggle to effectively
detect cyber attacks such as advance persistent threats (APTs) and zero-day
attacks. Cyber attacks have successfully compromised systems in a wide range
of sectors. For example, the WannaCry ransomware attack extorted money to
unlock sensitive information and designs across various industries [20].
Security experts have increasingly turned to sharing cyber threat intelligence
(CTI) to combat such emerging cyber threats. CTI is any relevant information
that helps detect, monitor, assess, and respond to cyber threats. CTI
facilitates a comprehensive and significant threat warning and includes
information such as IOCs [21].

Nowadays, a rich source of commercial and free CTI feeds are available, making
it difficult for network defenders to evaluate the quality of information and
select the optimal set of data feeds to pay attention to. Acting on results
from low-quality feeds can give rise to many false alerts while concentrating
on only a few data feeds increases the risk of missing relevant threats.
However, it is challenging to extract IOCs from unstructured form sources.
Several automated methods for extracting IOCs (such as malicious IP addresses,
malware, and file hashes of malicious payloads) are based on the OpenIOC
standard, including PhishTank, IOCFinder, and CleanMX [14]. To facilitate
efficient threat intelligence sharing among organisations, CybOX [22], STIX
[23], and TAXII [24] have emerged as de-facto standards for describing threat
intelligence and are widely consumed by the threat intelligence sharing
platforms, including MISP [12] and AT&T Open Threat eXchange (OTX).

###  3.2 Indicators of Compromise

Cyber Threat Intelligence (CTI) includes IOCs, which organisations can use to
identify possible threats and protect themselves and their customers.
Specifically, IOCs are artefacts observed about an attacker or their
behaviour, such as tactics, techniques and procedures [25]. IOCs can be kept
at a network or host level and help network defenders block malicious traffic
and identity actions or determine if a cyber intrusion has occurred. Security
and forensic analysts prepare reports of in-depth analysis of cyber attacks,
including the IOCs, to be shared with communities, often through public data
sources. Examples of IOC found in reports from data sources include actor
identity behind cyber attacks, the malware used in threat attacks and their
typical behaviour, communication and control server list, and other types of
information. The information used in creating these reports is gathered from
multiple sources, such as host logs, proxy logs and alerts. The reports may be
widely distributed through various channels, including blogs, forums and
social media.

The pyramid of pain (PoP) classifies the common types of IOCs. The PoP
identifies the types of indicators that a system defender might use to detect
an adversary’s activities. The pyramid organises the pain an adversary will
cause when the defender can deny those indicators. At the bottom end, if a
defender identifies hash values of malicious files and then blocks these, it
causes the attacker little pain since making an insignificant change to the
file to produce the same outcome with a different hash is trivial. TTP sit at
the top of the pyramid. When a defender detects and responds at this level,
this disrupts behaviours much more complicated for an adversary to change;
defining new behaviours is a significant challenge for adversaries.

###  3.3 Heterogeneous Information Network

Heterogeneous Information Network (HIN) is a simple way of modelling a problem
as a graph compromising different types of nodes and one or more correlations
between nodes (edges) [26]. The set of node and edge types correspond to the
network scheme. HIN delivers a high conceptualisation of modelling for a
complex collection of data. From the graphical representation of the dataset,
feature vectors can be extracted by defining meta-paths and meta-graphs
corresponding to the graph and implementing a guided random walk over defined
meta-paths and meta-graphs. A meta-path is a path defined within the graph of
network schema, covering a specific sequence of relation types. A meta-graph
[27] can handle the in-depth relationship between nodes by employing a direct
acyclic graph of nodes defined over the HIN from a single source node to a
single target node. The guided random walk generates a sequence of nodes
processed in an embedding model such as word2vec, skip-gram or Continuous Bag-
of-Words (CBOW). Once the nodes are represented numerically, it is possible to
determine a set of nodes and resolve many problems (classification,
clustering, and similarity search).

###  3.4 Overview

We introduce a novel system designed to automatically extract and prioritise
high-level IOCs (Indicators of Compromise) from multiple sources of threat
text. Our system addresses the limitations of existing IOC extraction methods
by considering the semantic relationships among different IOCs. We present a
novel approach to extracting threat-related information and identifying the
domains IOCs belong to. This information is then integrated with their
respective domains to form three domain-specific threat intelligence
categories: the organisational domain, regional-source domain, and regional-
target domain. We also present a threat modelling that utilizes a
Heterogeneous Information Network (HIN) comprising different meta-paths and
meta-graphs. The proposed system captures the interdependent relationships
among diverse IOCs and provides a more comprehensive view of the landscape of
threat events. Our system then utilizes similarity measures based on these
meta-paths and meta-graphs to assess the interdependent relationships among
different IOCs.

To prioritize IOCs, we measure their severity by considering various factors,
including the frequency of attacks, the lifetime of the IOC, and the number of
exploited vulnerabilities in each domain. Our system then evaluates the
ranking mechanism for each IOC, providing a more comprehensive and accurate
view of the threat landscape. Our system significantly contributes to
cybersecurity, providing a more effective and efficient method for
automatically extracting, assessing, and prioritizing high-level IOCs. With
the increasing frequency and complexity of cyber threats, the need for such a
system has become more critical.

##  4 Methodology

The architecture of the proposed method, as shown in Fig 2, comprises of
following phases: Data collection and Preprocessing, Relation Extraction and
Threat Modelling, Domain Recognition and Tag Generation, Domain-specific
threat identification and Tagging, and Severity measure and Threat
Prioritisation. Table 1 summarises the list of notations and abbreviations
used throughout the paper.

Table 1: List of Abbreviations and Notations Notations | Description
---|---
IOC | Indicators of Compromise
AM | Additive manufacturing
HIN | Heterogeneous Information Network
APT | Advanced Persistent Threats
CTI | Cyber Threat Intelligence
TTP | Tactic, technique and procedure
PoP | Pyramid of pain
STIX | Structured threat information exchange
TAXII | Trusted Automated eXchange of Indicator Information

###  4.1 Data Collection and Preprocessing

The system automatically collects threat information identifying IOCs from
multiple resources, including forums, blogs, security news, and bulletins. We
use a breadth-first search to capture the HTML course code and Xpath for data
extraction. We then reduce the dimension of each text report and remove noisy
features by pre-processing. This pre-processing includes the removal of
stopwords, punctuations, and markup characters.

Figure 2: Process flow of the Proposed system

###  4.2 Relation extraction and threat intelligence modelling

Using a heterogeneous information network (HIN) for threat intelligence, we
first build a graph that shows the interdependent (semantic) relationships
between the different IOCs involved in the attack. By denoting nodes (IOCs)
and relationships, we can identify patterns and anomalies that may indicate
the presence of a threat. For example, we can use HINs to identify groups of
attackers that share common attack vectors or targets or to track the
evolution of an attack over time as new entities become involved. To better
understand, we can characterise the nodes and relations as follows.

####  4.2.1 Node Features.

In the context of risk in Additive Manufacturing, it is essential to consider
the domain-specific threat information. For instance, a threat post discussing
the Stuxnet virus and its impact on industrial control systems is more
relevant to manufacturing organisations than those in the finance or
healthcare sectors. This highlights the need for threat intelligence tailored
to an organisation’s domain.

Additionally, geographical location plays a significant role in cyber attacks.
Over 500 geopolitical cyber attacks have been reported worldwide in the past
decade, with 30% originating from China or Russia and 26.3% targeting the USA.
In 2018 alone, 27% of attacks occurred in the USA [28]. Therefore, when
developing threat models for Additive Manufacturing, it is crucial to consider
the regional source and target source of cyber attacks.

To account for these domain-specific and regional factors in our threat
intelligence model for Additive Manufacturing, we define nodes as
organisation-specific, regional_source-specific, and regional_target-specific.
This enables us to capture the complex relationships between entities involved
in cyber attacks, such as attackers, attack vectors, and targets. Moreover, we
consider time-related node features such as attack frequency and IOC
lifecycle, which can provide valuable insight into the TTPs of attackers and
help defenders calculate the level of risk posed by a particular threat.

####  4.2.2 Semantic relation features.

The node features in the HIN represent a specific action, but actions can be
employed multiple times in conjunction with other activities in a campaign.
These complex relationships among nodes can provide more valuable intelligence
for threat identification; therefore, we consider relation-based and node
features. This allows us to analyse highly sophisticated malicious cyber
attacks. To model the interdependent relationship between eight IOCs, we
define the following semantic relationships:

  * •

R1: The relation actor-exploit-vulnerability matrix $A$ represents the link
between the threat actor and vulnerability. For each element,
$A_{i,j}\in\\{0,1\\}$, where $A_{i,j}=1$ means actor $i$ exploits
vulnerability $j$.

  * •

R2: The relation actor-invade-device matrix $B$ represents the link between
the threat actor and device. For each element, $B_{i,j}\in\\{0,1\\}$, where
$B_{i,j}=1$ means actor $i$ invades device $j$.

  * •

R3: The link between two actors is represented by the relation actor-assist-
actor matrix $C$. For each element, $C_{i,j}\in\\{0,1\\}$, where $C_{i,j}=1$
means actor $i$ assists actor $j$.

  * •

R4: The relation attack_type-originate_from-region matrix $D$ represents the
link between the attack type and location. For each element,
$D_{i,j}\in\\{0,1\\}$, where $D_{i,j}=1$ means attack type $i$ originate from
region $j$.

  * •

R5: The relation attack_type-target-region matrix $E$ represents the link
between the attack type and location. For each element, $E_{i,j}\in\\{0,1\\}$,
where $F_{i,j}=1$ means attack type $i$ target to region $j$.

  * •

R6: The relation vulnerability-affect-device matrix $F$ represents the link
between the vulnerability and the device. For each element,
$F_{i,j}\in\\{0,1\\}$, where $F_{i,j}=1$ means vulnerability $i$ affects
device $j$.

  * •

R7: The relation attack_type-associate-vulnerability matrix $G$ represents the
link between the attack type and vulnerability. For each element,
$G_{i,j}\in\\{0,1\\}$, where $G_{i,j}=1$ means attack type $i$ carry
vulnerability $j$.

  * •

R8: The relation vulnerability-held-time matrix $H$ represents the link
between the vulnerability and time. For each element, $H_{i,j}\in\\{0,1\\}$,
where $H_{i,j}=1$ means vulnerability $i$ held in time $j$.

  * •

R9: The relation vulnerability-include -file matrix $B$ represents the link
between the vulnerability and malicious file. For each element,
$I_{i,j}\in\\{0,1\\}$, where $I_{i,j}=1$ means vulnerability $i$ include
malicious file $j$.

R10: The relation vulnerability-evolve-vulnerability matrix $B$ represents the
link between the vulnerabilities. For each element, $J_{i,j}\in\\{0,1\\}$,
where $J_{i,j}=1$ means vulnerability $i$ evolve to vulnerability $j$.

We initiate dependency parsing to leverage the semantic relationships among
the eight IOCs and extract them in a structured format. Using this approach,
we can represent the IOCs as triplets, each consisting of two IOCs and a
relation between them. For instance, if IOC1 is dependent on IOC2, we would
define the relationship as (IOC1-relation-IOC2), where ’relation’ denotes the
nature of the relationship between the two IOCs.

####  4.2.3 Meta-path and Meta-graph.

Figure 3 presents 12 distinct types of meta-paths and meta-graphs denoted by
$\chi_{i}$ that capture interdependent relationships among seven different
IOCs. While the meta-path illustrates the connections between the IOCs, it
falls short in capturing intricate relationships. To address this limitation,
the proposed HIN-based Threat Intelligence (TI) model utilizes a directed
acyclic graph of nodes to handle more complex structures in the HIN
architecture. By learning and analyzing these 12 different meta-paths and
meta-graphs, the model can convey the context of a threat event and offer
threat insights across heterogeneous IOCs. For instance, the $\chi_{1}$ meta-
path is a length-2 meta-path that represents the relatedness of ”threat actors
(A) exploiting the same vulnerability (V).”

Figure 3: Proposed meta-path and meta-graph for threat type identification,
where A denotes threat actor, V denotes vulnerability, D denotes device, R
denotes a region, AT denotes attack type, F denotes file, and T denotes time.

Similarly, $\chi_{8}$ is a meta-path that describes the relationships between
IOCs that “two attack types who leverage the same vulnerability held at the
same time”. Likewise, $\chi_{10}$ is a meta-graph that portrays the
relationship over threat infrastructure with more comprehensive insight that
integrates both external and intrinsic relationships. Meta-graph $\chi_{10}$
depicts the relationship among IOCs: ”two attack types originated from the
same region, and their associated vulnerabilities affect the same device occur
at the same time”.

###  4.3 Domain recognition and tag generation

To extract domain-specific IOCs, it is essential first to identify the domain
of threat information. This initial step helps to ensure that IOCs are
tailored to the specific context of the threat landscape, enabling more
effective threat detection and response. Here, we consider three domains,
organisation_domain-specific, regional_source-specific, and regional_target-
specific. Organisational_domain-specific threat information includes
financial, health, manufacturing, government, and IoT information.
Regional_source-specific and Regional_target-specific threat information
originated and targeted the geographic region, such as China, Russia, India,
Korea, USA, UK, and Europe. We first trained the word2vec model specific to a
threat description embedding that inputs a large corpus (threat description)
and generates a low-dimension vector. Each unique word in the corpus is
allocated an assigned vector in latent space. The convolution function sets a
filter for each word vector to generate a feature called local_feature. Our
model captures the most significant features by leveraging the max-pooling
operation that takes the maximum values of local-feature and runs over a local
feature set. This will generate tags for three domains, i.e., $OD_{t}$,
$SD_{t}$, and $TD_{t}$ denotes the tags corresponding to organisation_domain-
specific, regional_source-specific, and regional_target-specific threat,
respectively.

###  4.4 Domain-Specific Threat Identification and Tagging

After successfully extracting the features of IOCs and their relationships and
identifying the relevant meta-paths and meta-graphs, a meta-path and meta-
graph-based heterogeneous classifier is developed to classify the threat type
of infrastructure nodes in Cyber Threat Intelligence (CTI). The proposed
classification approach integrates node features and explores the semantic
similarities between meta-paths and meta-graphs to represent the nodes
comprehensively. These advanced techniques enable a more comprehensive
depiction of the nodes, enhancing the accuracy of the threat classification.

Given the threat intelligence graph $G=(N,R)$, and meta-path and meta-graph
set $P=\\{\chi_{1},\chi_{2},..\chi_{n}\\}$, the assessment of threat
intelligence includes the following steps:

  * •

Adjacency matrix. The relationships between threat nodes can be explored using
different meta-paths and meta-graphs, which capture the behaviour of threats
in various aspects. To represent these relationships, we propose using an
adjacent weighted matrix, denoted by $Adj_{i}\in R^{N\times N}$, which can be
generated using similarity algorithms such as Euclidean distance, Manhattan
distance, Cosine similarity, word mover distance, and Jaccard similarity. To
assess the similarity between IOCs, we generate a corresponding weighted
adjacency matrix, $Adj_{i}$, based on the meta-path and meta-graph path set,
$P$. The use of weighted adjacency matrices $Adj_{i}$ enables the
identification of the most significant relationships between the different
nodes, which can be used to prioritise threat mitigation efforts.

  * •

Feature matrix. By incorporating attributed information of nodes, we can
construct an attributed feature matrix $F_{i}$ of size $N\times d$, where N
denotes the number of IOCs in $Adj_{i}$, and $d$ is the dimension of the node
feature. This allows us to integrate the attribute information of IOCs and
create a node feature matrix $F_{i}\in R^{N\times d}$. To recognize previously
unnoticed IOCs, we employ the word2vec method to develop a threat intelligence
embedding, which transforms words into a latent vector space. To achieve this,
threat-related texts are pre-processed, accumulated into a word set, and
converted into a latent vector space using word2vec. This approach enables us
to represent threat-related information in a low-dimensional vector space,
facilitating the detection and analysis of potential threats.

  * •

Quantify threat intelligence. After designing an adjacent weighted and
attributed feature matrix, we assess the threat intelligence. Different types
of assessment methods to quantify the proposed HIN-based threat intelligence.
For example, graph convolution network (GCN) and Bidirectional Encoder
Representations from Transformers (BERT). Given the adjacency matrix
$Adj_{i}$, and its corresponding feature matrix $F_{i}$ (low-dimensional
space), we utilise the graph convolution network (GCN) method to quantify the
relationship between IOCs. This will fuse the adjacency matrix $Adj_{i}$, and
feature matrix $F_{i}$ as $Z=(F,Adj)$ and output the predicted labels of IOCs.
Then, the model integrates the domain-specific tags $OD_{t}$, $SD_{t}$ and
$TD_{t}$ to the predicted IOC, representing that the IOCi belongs to the
organisation $OD_{t}$, originates from the country $SD_{t}$ and is targeted to
the country $TD_{t}$ and are considered as the domain-specific IOCs.

###  4.5 Severity Measure and Threat Prioritisation

Utilizing the learned domain-specific IOCs, we can evaluate the severity of
potential threats of various attack vectors within each domain. This motivates
us to develop a quantitative measure to assess threat risks corresponding to
each domain. The proposed severity measure is based on several key
assumptions.

  1. 1.

Firstly, we assume that the frequency of attacks may significantly influence
the severity and scope of the threats manifested.

  2. 2.

Secondly, we postulate that chain exploits, where multiple vulnerabilities use
an attack, can cause considerably more damage.

  3. 3.

Finally, we recognize that the severity of a threat may decrease over time,
particularly during the zero-day risk period.

Consequently, the severity of a threat can be measured by examining the
frequency of attacks, the lifetime of the IOC, and the number of exploited
vulnerabilities in each domain. This approach allows us to develop a more
nuanced and comprehensive understanding of the potential threats facing a
domain, enabling us to take appropriate measures to mitigate risk and enhance
security.

##  5 Conclusion and Future scope

This paper presented a novel semantic-based threat prioritisation system for
AM that intends to expose the comprehensive behaviour of threat events in the
relationships among different IOCs with high accuracy. We proposed an
intelligent IOC acquisition and ranking system based on a Heterogeneous
Information Network (HIN). The proposed system collects threat-related from
multiple sources that automatically extract threat-related information,
measure the severity of IOCs, and quantify them based on the severity. We
considered individual IOCs and one or more relationships among semantically
similar IOCs. We proposed an efficient recogniser to identify domain-specific
IOCs focusing on three domains: organisational_domain-specific,
regional_source-specific, and regional_target-specific threat intelligence.
Further, we evaluated the severity of IOC by exploring the frequency of
attacks, IOC lifetime, and the number of exploited vulnerabilities in each
domain.

The proposed semantic-based threat prioritisation system for AM has potential
future scopes that can be explored, such as:

  * •

Integrating with existing security tools: The proposed system can be combined
with existing security tools to provide real-time threat intelligence and
prioritisation of threats. The integration can help security teams to automate
the detection, investigation, and response to threats and reduce the time to
mitigate them.

  * •

Exploring additional domains: The proposed system focuses on three domains:
organisational_domain-specific, regional_source-specific, and regional_target-
specific. However, other domains, such as industry-specific or technology-
specific, can be explored to provide a more comprehensive view of the threat
landscape.

  * •

Improving the ranking system: The proposed system ranks the IOCs based on
severity. However, the ranking system can be improved to consider the evolving
threat landscape and real-time threat intelligence data to enhance the
accuracy of the prioritisation system.

## References

  * [1] Haleem, Abid, and Mohd Javaid. ”Additive manufacturing applications in industry 4.0: a review.” Journal of Industrial Integration and Management 4.04 (2019): 1930001.
  * [2] M. Dietz and G. Pernul, “Unleashing the digital twin’s potential for ICS security,” IEEE Secur. Priv., vol. 18, no. 4, pp. 20–27, 2020.
  * [3] R. Sobers, “Data Breach Response Times: Trends and Tips — Varonis,” Varonis, 2020. https://www.varonis.com/blog/data-breach-response-times/ (accessed Jun. 17, 2020).
  * [4] R. Langner, “Stuxnet: Dissecting a cyberwarfare weapon,” IEEE Secur. Priv., vol. 9, no. 3, pp. 49–51, 2011.
  * [5] J. Slay and M. Miller, “Lessons learned from the maroochy water breach,” in International conference on critical infrastructure protection, 2007, pp. 73–82.
  * [6] S Belikovetsky, M. Yampolskiy, J. Toh, J. Gatlin, and Y. Elovici. ”dr0wned–cyber-physical attack with additive manufacturing.” In 11th USENIX Workshop on Offensive Technologies (WOOT). 2017.
  * [7] Kumar, Mahender, and Satish Chand. ”A provable secure and lightweight smart healthcare cyber-physical system with public verifiability.” IEEE Systems Journal 16.4 (2021): 5501-5508.
  * [8] T. D. Wagner, K. Mahbub, E. Palomar, and A. E. Abdallah, “Cyber threat intelligence sharing: Survey and research directions,” Comput. Secur., vol. 87, p. 101589, 2019.
  * [9] I. B. M. X-Force, “IBM X-force threat intelligence report 2016,” 2016.
  * [10] Github, “Facebook ThreatExchange Report.” https://github.com/facebook/ThreatExchange (accessed Jul. 22, 2021).
  * [11] “OpenCTI.” https://www.opencti.io/en/ (accessed Jul. 22, 2021).
  * [12] C. Wagner, A. Dulaunoy, G. Wagener, and A. Iklody, “Misp: The design and implementation of a collaborative threat intelligence sharing platform,” in Proceedings of the 2016 ACM on Workshop on Information Sharing and Collaborative Security, 2016, pp. 49–56.
  * [13] U. Noor, Z. Anwar, T. Amjad, and K.-K. R. Choo, “A machine learning-based FinTech cyber threat attribution framework using high-level indicators of compromise,” Futur. Gener. Comput. Syst., vol. 96, pp. 227–242, 2019.
  * [14] J. Zhao, Q. Yan, J. Li, M. Shao, Z. He, and B. Li, ”TIMiner: Automatically extracting and analysing categorised cyber threat intelligence from social data,” Comput. Secur., vol. 95, p. 101867, 2020.
  * [15] Y. Gao, L. I. Xiaoyong, P. Hao, B. Fang, and P. Yu, “Hincti: A cyber threat intelligence modeling and identification system based on heterogeneous information network,” IEEE Trans. Knowl. Data Eng., 2020.
  * [16] J. Zhao, Q. Yan, X. Liu, B. Li, and G. Zuo, “Cyber threat intelligence modeling based on heterogeneous graph convolutional network,” in 23rd International Symposium on Research in Attacks, Intrusions and Defenses (RAID 2020), 2020, pp. 241–256.
  * [17] X. Liao, K. Yuan, X. Wang, Z. Li, L. Xing, and R. Beyah, “Acing the ioc game: Toward automatic discovery and analysis of open-source cyber threat intelligence,” in Proceedings of the 2016 ACM SIGSAC Conference on Computer and Communications Security, 2016, pp. 755–766.
  * [18] P. Gao et al., “Enabling efficient cyber threat hunting with cyber threat intelligence,” in 2021 IEEE 37th International Conference on Data Engineering (ICDE), 2021, pp. 193–204.
  * [19] X. Wang et al., “A Method for Extracting Unstructured Threat Intelligence Based on Dictionary Template and Reinforcement Learning,” in 2021 IEEE 24th International Conference on Computer Supported Cooperative Work in Design (CSCWD), 2021, pp. 262–267.
  * [20] J. M. Ehrenfeld, “Wannacry, cybersecurity and health information technology: A time to act,” J. Med. Syst., vol. 41, no. 7, p. 104, 2017.
  * [21] B. Akram and D. Ogi, “The Making of Indicator of Compromise using Malware Reverse Engineering Techniques,” in 2020 International Conference on ICT for Smart Society (ICISS), 2020, pp. 1–6.
  * [22] S. Barnum, R. Martin, B. Worrell, and I. Kirillov, “The cybox language specification,” MITRE Corp., 2012.
  * [23] S. Barnum, ”Standardising cyber threat intelligence information with the structured threat information expression (stix),” Mitre Corp., vol. 11, pp. 1–22, 2012.
  * [24] J. Connolly, M. Davidson, and C. Schmidt, “The trusted automated exchange of indicator information (taxii),” MITRE Corp., pp. 1–20, 2014.
  * [25] A. Iklody, G. Wagener, A. Dulaunoy, S. Mokaddem, and C. Wagner, “Decaying indicators of compromise,” arXiv Prepr. arXiv1803.11052, 2018.
  * [26] T. Fu, W.-C. Lee, and Z. Lei, “Hin2vec: Explore meta-paths in heterogeneous information networks for representation learning,” in Proceedings of the 2017 ACM on Conference on Information and Knowledge Management, 2017, pp. 1797–1806.
  * [27] H. Zhao, Q. Yao, J. Li, Y. Song, and D. L. Lee, “Meta-graph based recommendation fusion over heterogeneous information networks,” in Proceedings of the 23rd ACM SIGKDD international conference on knowledge discovery and data mining, 2017, pp. 635–644.
  * [28] J Robinson, “Cyberwarfare statistics: A decade of geopolitical attacks,” Privacy Affairs, 2019, https://www.privacyaffairs.com/ geopolitical-attacks/ Accessed 5 September 2021.


# Calculating the polarization in bi-partite lattice models: application to an
extended Su-Schrieffer-Heeger model

Balázs Hetényi1,2, Yetkin Pulcu1, and Serkan Doǧan3 1 Department of Physics,
Bilkent University, TR-06800 Bilkent, Ankara, Turkey
2MTA-BME Quantum Dynamics and Correlations Research Group, Department of
Physics, Budapest University of Technology and Economics, H-1111 Budapest,
Hungary
3 Department of Mathematics, Bilkent University, TR-06800 Bilkent, Ankara,
Turkey

###### Abstract

We address the question of different representation of Bloch states for
lattices with a basis, with a focus on topological systems. The
representations differ in the relative phase of the Wannier functions
corresponding to the diffferent basis members. We show that the phase can be
chosen in such a way that the Wannier functions for the different sites in the
basis both become eigenstates of the position operator in a particular band. A
key step in showing this is the extension of the Brillouin zone. When the
distance between sites within a unit cell is a rational number, $p/q$, the
Brillouin extends by a factor of $q$. For irrational numbers, the Brillouin
zone extends to infinity. In the case of rational distance, $p/q$, the Berry
phase “lives” on a cyclic curve in the parameter space of the Hamiltonian, on
the Brillouin zone extended by a factor of $q$. For irrational distances the
most stable way to calculate the polarization is to approximate the distance
as a rational sequence, and use the formulas derived here for rational
numbers. The use of different bases are related to unitary transformations of
the Hamiltonian, as such, the phase diagrams of topological systems are not
altered, but each phase can acquire different topological characteristics when
the basis is changed. In the example we use, an extended Su-Schrieffer-Heeger
model, the use of the diagonal basis leads to toroidal knots in the
Hamiltonian space, whose winding numbers give the polarization.

††preprint: APS/123-QED

##  I Introduction

Recently, several papers investigated Bena09 ; Watanabe18  the different
representation of the wave function and physical observables when the
underlying model is a lattice with a basis. This fundamental question is
important per se, moreover, such models have been the subject of intense
research in the last decades, exemplified by graphene, the kagome, or Lieb
lattices, and others. The canonical models for topological insulators
Bernevig13 ; Kane13 ; Asboth15  are all lattices with bases: the Haldane,
Haldane88  Kane-Mele, Kane05a ; Kane05b  Su-Schrieffer-Heeger Su79  models are
some of the most fundamental ones.

Bena and Montambaux Bena09  have studied the two commonly used bases of the
graphene lattice. The difference between the two is that in one (basis $I$),
the different position of the sites within a unit cell are not explicit (they
are taken to be the same, the position of the unit cell itself), whereas in
the other (basis $II$), the positions of the different sites within a unit
cell are explicit. What is meant by this is that in $I$, the Wannier functions
corresponding to the different basis members have the same phase, while in
$II$ their relative phase depends on the distance between them. Bena and
Montambaux Bena09  show that in the case of the graphene lattice, although the
Hamiltonian and operators take different forms in the two bases, important
physical quantities, such as the density or the density of states, as well as
the low-energy theory, are unchanged. However, care must be taken in using the
correct construction of observables, once a representation of the wave
functions is chosen.

Watanabe and Oshikawa Watanabe18  have taken up this question more recently,
and they show that, while the precise form of the current operator does depend
on the choice of basis, the quantized charge pumped during an adiabatic cycle
Thouless82  is unchanged. Note that the current operator is intimately related
to the phase of the wave function, in particular the relative phase wave
functions centered on different sites may exhibit. Also note that the charge
pumped is a topological invariant, which remains the same in different bases.

In this paper we reconsider the question of basis dependence in lattice models
in the context of modern polarization theory with a focus on topological
systems. We compare the two bases used in Ref. Bena09 , and show that one
corresponds to a linear combination of Wannier functions which are
eigenfunctions of the position operator. When the distance between the basis
sites is a rational number $p/q$, the Brillouin zone has to be extended by a
factor of $q$ to derive the above result. If the distance is irrational, then
the Brillouin zone has to be extended to the whole real line.

We then turn to the question of topological phase transitions in the two
bases. As an example system we consider an extended Su-Schrieffer-Heeger Su79
(SSH) model. When the $k$-space Hamiltonian is written considering distance
dependence explicitly the three-dimensional curve $d_{x},d_{y},d_{z}$ (see Eq.
7) is not a closed curve. When the distance between basis sites is a rational
number ($p/q$) extending the Brillouin zone by a factor of $q$ closes the
curve, and an expression for the polarization Berry phase Berry84 ; Zak89 ;
King-Smith93 ; Resta98 ; Resta00  can be derived. In the irrational case a
closed $d_{x},d_{y},d_{z}$-curve results only if the Brillouin zone is
extended to infinity. For this reason, calculation of the irrational case is
not straightforward. We show that a stable calculation results if the
irrational distance is approximated as a rational sequence.

Implementing the dependence of the Hamiltonian on the distance between basis
sites proceeds through a unitary transformation. As such, the Hamiltonian can
acquire auxiliary topological features. In the model we study, the
$d_{x},d_{y},d_{z}$-curve traces out a toroidal knot, and the ratio of the
winding numbers of the knot give the polarization. However, the phase diagram
itself is not altered, since these topological characteristics are merely
different representations of the phases of the original extended SSH model.

In the next section we discuss the two different bases used in lattices with a
basis. In section III the extended SSH model is described. In section IV the
change of basis is discussed in the context of the extended SSH model. In
section V our results are presented, and in section VI we conclude our work.

##  II Polarization and choice of representation

It is well-known Blount62 ; Kivelson82  that Wannier functions in systems
without a basis are eigenstates of the position operator within a given band.
Here we show that in a lattice with a basis, one of the representations of the
wave-function studied by Bena and Montambaux Bena09  leads to position
eigenstates.

We consider a one-dimensional system whose potential obeys $V(x)=V(x+1)$ (in
other words we take the length of a unit cell as unity). The unit cell
consists of a two-site basis. We chose the distance between the two internal
sites to be $p/q$, a rational number, and we extend the Brillouin zone by a
factor of $q$. This step is justified below.

In this case, the two commonly used bases Bena09  can be written as

$\displaystyle u_{k}^{I}(x)=$
$\displaystyle\sum\limits_{m=-\infty}^{\infty}e^{ik(m-x)}\left[w_{a}(x-m)+w_{b}\left(x-m-\frac{p}{q}\right)\right]$
$\displaystyle u_{k}^{II}(x)=$
$\displaystyle\sum\limits_{m=-\infty}^{\infty}\left[e^{ik(m-x)}w_{a}(x-m)+e^{ik(m+p/q-x)}\times\right.$
(1) $\displaystyle\left.\times w_{b}\left(x-m-\frac{p}{q}\right)\right].$

The functions $u_{k}^{I}(x)$ and $u_{k}^{II}(x)$ are periodic Bloch functions.
$u_{k}^{I}(x)$ is the more standard approach, and it appears Ashcroft76  in
textbooks. $w_{a}(x),w_{b}(x)$ are two different localized Wannier functions.
In the extended SSH model we study below, they differ only in having their
centers displaced with respect to each other.

Following Zak Zak89  we derive the geometric phase which arises upon
integrating the connection across the Brillouin zone. We define a geometric
phase of the Zak type, in the extended Brillouin zone as:

$\gamma=\frac{i}{2\pi q}\int_{-q\pi}^{q\pi}dk\langle
u_{k}|\partial_{k}|u_{k}\rangle.$ (2)

Using $u_{k}^{I}$ the Zak phase becomes

$\gamma_{I}=\int_{-\infty}^{\infty}dxx\left|w_{a}(x)+w_{b}\left(x-\frac{p}{q}\right)\right|^{2},$
(3)

whereas with $u_{k}^{II}$

$\gamma_{II}=\int_{-\infty}^{\infty}dxx\left(|w_{a}(x)|^{2}+|w_{b}(x)|^{2}\right).$
(4)

In deriving $\gamma_{II}$, we made use of the fact that

$\frac{1}{2\pi q}\int_{-\pi q}^{\pi
q}dk\exp\left(ik(n-n^{\prime})/q\right)=\delta_{nn^{\prime}}.$ (5)

Comparing $\gamma_{I}$ and $\gamma_{II}$ we see that the second Bloch basis
function, $u_{k}^{II}$ is a diagonal basis for the position operator in a
periodic system, since the result is a simple sum of Wyckoff positions for the
two different Wannier functions. It is due to Eq. (5) that the cross terms in
$\gamma_{II}$ disappear, but only of the Brillouin zone is extended by a
factor of $q$.

It is instructive to consider the effect of inversion symmetry on either
$\gamma$. In the original work of Zak Zak89 , where the sites are equivalent,
two types of inversion symmetries were considered, site-inversion and bond-
inversion. The former gives a Zak phase of zero, while the latter gives $\pi$,
which corresponds to a polarization of $qa/2$ ($a$ is the length of the unit
cell). In our case, there are also two types of inversion. If the sites of the
model are a distance $p/q$ apart, one can consider inverting around the
midpoint of the $p/q$ bond, or the $1-p/q$ bond. In the former case it holds
that $w_{c}(-x+p/q)=\pm w_{c}(x)$ for both $c=a,b$, which leads to
$\gamma=\frac{p}{2q}$ (for both $I$ and $II$). In deriving this one can follow
the stame steps as ZakZak89 . Inverting around the other bond center leads to
$w_{c}(-x-1+p/q)=\pm w_{c}(x)$, from which it follows that
$\gamma=\frac{p-q}{2q}$ (again for both $I$ and $II$).

In the case of irrational distance, the Brillouin zone has to be extended to
the range $(-\infty,\infty)$. When this is done, Eq. (5) eliminates the cross
terms and leads to $\gamma_{II}$.

##  III Model

The full model we consider is an extension of the well-known Su-Schrieffer-
Heeger model Su79 . It is a bi-partite lattice model ($A$ and $B$
sublattices), whose Hamiltonian has the form,

$\hat{H}=-\sum_{i=1}^{L}[Jc_{i}^{\dagger}d_{i}+J^{\prime}d_{i}^{\dagger}c_{i+1}+iKc_{i}^{\dagger}c_{i+1}-iKd_{i}^{\dagger}d_{i+1}+\mbox{H.c.}],$
(6)

where $L$ denotes the number of unit cells, $c_{i}$($c_{i}^{\dagger}$) denote
the annihiliation(creation) operators on the $A$ site of the $i$th unit cell,
and $d_{i}$($d_{i}^{\dagger}$) denote the same for the $B$ site of the $i$th
unit cell. We take the length of a unit cell to be unity. An extended SSH
model of a different kind, with real second nearest neighbor hoppings, was
invesigated by Li et al. Li14 .

When $K=0$ the SSH model is recovered. The third and fourth terms in Eq. (6)
correspond to second nearest neighbor hoppings (hoppings within one
sublattice). A Peierls phase of $\pi/2$ and $-\pi/2$ is applied along these
hoppings, the former for sublattice $A$, the latter for sublattice $B$. These
second nearest neighbor hoppings play a similar role to those in the Haldane
model Haldane88 . There and in Eq. (6) the Peierls phases of second nearest
neighbor hoppings point in opposite directions on the two sublattices. The
consequence of this term is that for gapped phases and finite $K$ there is
persistent current in the system, but the currents on the different
sublattices cancel, and insulation results. If $K$ changes sign, the direction
of the persistent currents on the different sublattices also change sign. In
this sense, the model can be viewed as an example of Kohn’s tenet Kohn64 .
This tenet states that insulation in a quantum system is not a function of
localization of individual charge carriers, but ultimately results from many-
body localization, or the localization of the center of mass of all the charge
carriers. Individually, charge carriers can be quite mobile, as long as the
center of mass is localized.

Figure 1: Curves traced out by the parameters of the Hamiltonian
$d_{x},d_{y},d_{z}$ as $k$ traverses Brillouin zone for the usual SSH and the
extended SSH model. Upper panel: SSH model and extended SSH model in the usual
basis. Lower panel: extended SSH model in the distance dependent basis
comparing the first Brilllouin zone to the extended Brillouin zone.

According to the symmetry classification Altland97 ; Schnyder08  of
topological insulators, the SSH model, which exhibits time reversal, particle-
hole, and chiral symmetries, falls into the BDI class, which in one dimension
is a $\mathbb{Z}$-topological system. After the addition of the second nearest
neighbor hopping term in Eq. (7) only the particle-hole symmetry remains, and
the model falls in the D symmetry class, which is a $\mathbb{Z}_{2}$ system.
The gap closure points of the SSH remain (they occur at $k=0,\pi$), and we
show below that topological edge states still arise.

Fourier transforming Eq. (6) results in

$\hat{H}=\sum_{k}H_{k}=\sum_{k}[d_{x}(k)\sigma_{x}+d_{y}(k)\sigma_{y}+d_{z}(k)\sigma_{z}],$
(7)

where $\sigma_{x},\sigma_{y},\sigma_{z}$ denote the Pauli matrices, and

$\displaystyle d_{x}(k)$ $\displaystyle=$ $\displaystyle-J-J^{\prime}\cos(k),$
(8) $\displaystyle d_{y}(k)$ $\displaystyle=$
$\displaystyle-J^{\prime}\sin(k),$ $\displaystyle d_{z}(k)$ $\displaystyle=$
$\displaystyle 2K\sin(k).$

Written this way, we see that gap closure occurs at $k=\pi$ if $J=J^{\prime}$.
The upper panel Fig. 1 shows the curves traced out by $d_{x},d_{y},d_{z}$ as
$k$ traverses the Brillouin zone for the usual SSH model
($J=1.25,J^{\prime}=1.75$) with $K=0$ and an example of the extended one, with
$K=2$. The former forms a circle in the $xy$-plane, the latter is a circle
with a tilted axis (axis along the $yz$-plane). Both circles cross the
$x$-axis itself in two places at the same points. This means, that by varying
the ratio of $J/J^{\prime}$ it is possible to close the gap in both models,
and encounter topological phase transitions.

We can also develop a continuum approximation, the usual way Jackiw76 ,

$H=it\partial_{x}\sigma_{x}+iK\partial_{x}\sigma_{y}+2\delta t\sigma_{z},$ (9)

where we used the parametrization:

$\displaystyle t$ $\displaystyle=$ $\displaystyle\frac{J+J^{\prime}}{2}$ (10)
$\displaystyle\delta t$ $\displaystyle=$ $\displaystyle\frac{J^{\prime}-J}{2}$

We can solve for zero energy eigenstates of the form

$\begin{pmatrix}\Psi_{A}(x)\\\ \Psi_{B}(x)\end{pmatrix}.$ (11)

The zero energy solution is

$\Psi_{A/B}^{\prime\prime}(x)=\frac{4\delta t^{2}}{t^{2}+K^{2}}\Psi_{A/B}(x).$
(12)

The solutions are exponentials. If the system is extended, these solutions
diverge, hence they are unnormalizable. If the system has open boundary
conditions, exponential solutions can be normalized, since they can start at
one of the boundaries and decay towards the bulk. Therefore zero-energy states
are possile. The sign of $\delta t$ determines the direction of decay, in
other words, the side of the finite system on which the edge state is located.
The signs of $t$ and $K$ are irrelevant.

It is interesting to compare the above results to those of Jackiw and Rebbi
Jackiw76 . There a first order differential equation is solved for a function
multiplying a $\sigma_{z}$ eigenstate, while in our case each component of the
Pauli spinor satisfies the same second order differential equation. Our
calculations below on the lattice model (see Fig. 3) bear out the predictions
of our continuum theory.

##  IV Implementation of distance dependence

We take the shortest distance between the two sublattices to be $p/q$, where
$p$ and $q$ are co-prime integers. The Hamiltonian in $k$-space can be shown
to be:

$\displaystyle d_{x}(k)$ $\displaystyle=$
$\displaystyle-J\cos\left[k\left(\frac{p}{q}\right)\right]-J^{\prime}\cos\left[k\left(1-\frac{p}{q}\right)\right]$
(13) $\displaystyle d_{y}(k)$ $\displaystyle=$ $\displaystyle
J\sin\left[k\left(\frac{p}{q}\right)\right]-J^{\prime}\sin\left[k\left(1-\frac{p}{q}\right)\right]$
$\displaystyle d_{z}(k)$ $\displaystyle=$ $\displaystyle 2K\sin(k).$

This Hamiltonian can also be derived from the original one (Eq. (7)) by

$H_{k}\rightarrow\exp\left(ik\frac{\sigma_{z}}{2}\frac{p}{q}\right)H_{k}\exp\left(-ik\frac{\sigma_{z}}{2}\frac{p}{q}\right),$
(14)

a rotation by angle $kp/q$ around the $z$ axis. We also note that the vector
function $[d_{x},d_{y},d_{z}]$ has a period of $2\pi q$, rather than the usual
$2\pi$.

Due to the fact that the distance dependence coincides with a rotation, a
qualitative difference between the case of rational and irrational distances
arises. As is known Asboth15 ; Kane13 , the curve traced out by the SSH model
on the $d_{x}-d_{y}$ plane is an ellipse whose center, in general, lies on the
$d_{x}$ axis. Implementing the distance dependence will generate a different
curve. A closed curve can be obtained if the distance between sublattice sites
is taken to be a rational number ($p/q$ with $p$ and $q$ coprimes), and the
Brillouin zone is extended to $2\pi q$. Observables related to the Berry (Zak)
phase (polarization, topological invariants) depend on the curve being closed,
since these physical quantities are integrals over a connection over a closed
curve (usually the Brillouin zone).

Figure 2: Curves traced out by the parameters of the Hamiltonian
$d_{x},d_{y},d_{z}$ as $k$ traverses the extended Brillouin zone for the
$\frac{p}{q}=\frac{1}{2}$ case. The Hamiltonian parameters are $K=1$ in all
cases, and $J=1,J^{\prime}=2$ (upper panel), $J=1.5,J^{\prime}=1.5$ (middle
panel), $J=2,J^{\prime}=1$ (lower panel). The curves of the upper and lower
panels appear equivalent, but the orientation of the curves changes as the gap
closure point (the figure-8 of the middle panel) is crossed.

Let us rewrite the $d$-curve of Eq. (13) with a scaled wave vector,
$k=q\kappa$, as

$\displaystyle d_{x}(\kappa)$ $\displaystyle=$ $\displaystyle-\cos(\kappa
p)\left[J+J^{\prime}\cos(\kappa q)\right]$ $\displaystyle-\sin(\kappa
p)J\sin(\kappa q),$ $\displaystyle d_{y}(\kappa)$ $\displaystyle=$
$\displaystyle+\sin(\kappa p)\left[J+J^{\prime}\cos(\kappa q)\right]$
$\displaystyle-\cos(\kappa p)J^{\prime}\sin(\kappa q),$ $\displaystyle
d_{z}(\kappa)$ $\displaystyle=$ $\displaystyle+2K\sin(\kappa q).$

This curve exists on the surface of a torus, albeit, not a torus of the usual
parametrization, and the curve itself is a toroidal knot (see examples in
Figs. 2, 4, 5).

The fundamental group Nakahara03  of the torus is
$\pi_{1}(T^{2})=\pi_{1}(S^{1})\oplus\pi_{1}(S^{1})=\mathbb{Z}\oplus\mathbb{Z}$,
meaning that homotopically equivalent classes of loops on the torus can be
characterized by two integers. In Eq. (IV) the two integers are $p,q$, and in
this parametrization they are toroidal winding numbers. A torus knot is
obtained by a trajectory which loops around the hole of the torus an integer
number of times ($w_{1}=p$), while it makes an integer number of rotations
around the “body” of the torus itself ($w_{2}=q$). As argued in the previous
section, the ratio of winding numbers is proportial to the polarization. A
curve parametrized as Eq. (IV) traces out a toroidal knot for a given integer
pair $p,q$, even if $p$ and $q$ are not co-primes. If they are co-primes then
there is only one curve, if they are not, then the same curve is traced more
than once, as many times as the common divisor of $p$ and $q$. In the original
parametrization, Eq. (13), the rescaling leads to the same curves for a given
$p,q$.

Alternatively, one can define $r=q-p$ in Eq. (13), and arrive at a different
curve, similar in form to Eq. (IV). In this case the winding numbers will be
$w_{1}=r$ and $w_{2}=q$. Thus, we can write the two different fractionally
quantized polarization results of the previous section as:

$\mathfrak{P}=-\frac{e}{2}\frac{w_{1}}{w_{2}}.$ (16)

The lower panel of Fig. 1 shows two curves for this extended model, one using
a normal Brillouin zone ($-\pi,\pi$), the other an extended one ($-q\pi,q\pi$)
(model parameters are $J=1.25,J^{\prime}=1.75,K=2$ and $p/q=1/3$ in both
cases). The former curve is open, the latter is a closed curve, a toroidal
knot.

In our calculations below the Brillouin zone is stretched by a factor of $q$,
and we take every $q$th point in $k$ space in this enlarged Brillouin zone.
The starting point for the polarization expression is the scalar product

$Z_{q}=\langle\Psi|\exp\left(i\frac{2\pi\hat{X}}{L}q\right)|\Psi\rangle.$ (17)

Following the steps of Resta Resta98 , in a band insulator we have

$Z_{q}=\prod_{s=0}^{L-1}\mbox{det}S(k_{qs},k_{q(s+1)}).$ (18)

where

$S_{m,m^{\prime}}(k_{qs},k_{q(s+1)})=\int_{0}^{L}dx\hskip
2.84544ptu^{*}_{k_{qs},m}(x)u_{k_{q(s+1)},m^{\prime}}(x),$ (19)

where $u_{k_{qs},m}(x)$ denote the periodic Bloch functions of band $m$, and
$k_{qs}=2\pi qs/L$. In our numerics we calculate the polarization via

$\mathfrak{P}=-\frac{e}{2\pi}\mbox{Im}\ln\prod_{s=0}^{L-1}[\mbox{det}S(k_{qs},k_{q(s+1)})]^{\frac{1}{q}}.$
(20)

This definition is important from the point of view of implementation, because
the $\frac{1}{q}$ power keeps the terms in the Berry phase between $0$ and
$2\pi$, and provides the correct fractional polarization.

Again the question remains, whether irrational distances can be handled. The
difficulty in this case is that the Brillouin zone has to be extended to
infinity. Two methods suggest themselves. One, use the irrational distance
itself in place of $p/q$, and compare larger and larger Brillouin zones. In
this case, the curve in the $d_{x},d_{y},d_{z}$ curve will be open. Another
way would be to approximate the irrational distance as the limit of a sequence
of rational numbers, and investigate the limit of Eq. (20) along that
sequence. Below a comparison is presented for the case when the distance is
the inverse of the golden ratio.

Figure 3: Calculation for a system with open boundary conditions with $L=100$
(meaning also that there are $100$ states). The parameters of the Hamiltonian
are $J=1.25,J^{\prime}=1.75,K=2.0$. The minimum distance between sublattices
is $p/q=1/3$. The figure show: (a) the zero energy state number $50$, (b) the
zero energy state number $51$, (c) the phase of the wave function as a
function of lattice site for state number $50$, and (d) the phase of the wave
function as a function of lattice site for state number $51$.

##  V Results

The curves traced out by the parameters of the Hamiltonian
($d_{x},d_{y},d_{z}$) are shown in Fig. 2 for three cases: $J=1,J^{\prime}=2$;
$J=1.5,J^{\prime}=1.5$; $J=2,J^{\prime}=1$ ($K=1$ in all three cases) for
$p=1,q=2$. The topological transition occurs at $J=1.5,J^{\prime}=1.5$, where
the curve traces out a figure$-8$, with a crossing point at the origin. The
curves on the different sides of this phase diagram (the uppermost and
lowermost panels) can be transformed into each other via a rotation by
$\pi/2$. Since the orientation of the curves (the sense of rotation around the
$d_{z}$-axis) does change when crossing the phase boundary the winding number
goes from one to minus one. This corresponds to a polarization reversal from
$\mathfrak{P}=e/4$ to $\mathfrak{P}=-e/4$, as expected.

Figure 4: Curves traced out by the parameters of the Hamiltonian
$d_{x},d_{y},d_{z}$ as $k$ traverses the extended Brillouin zone. The
upper(lower) set of plots show the case
$\frac{p}{q}=\frac{1}{3}$($\frac{p}{q}=\frac{2}{3}$), where the parameters of
the Hamiltonian are, left: $J^{\prime}=1.75,J=1.25$,
center:$J=1.5,J^{\prime}=1.5$, and right:$J^{\prime}=1.25,J=1.75$ ($K=1$ in
all cases).

Fig. 3 presents a study of how edge states arise in a system with open
boundary conditions. The parameters of the Hamiltonian are
$J=1.25,J^{\prime}=1.75,K=2.0$ with $p=1$ and $q=3$. The upper two panels ((a)
and (b)) show the two zero energy wave functions which are localized at the
edges. The system is fairly small $L=100$, and the paramater $K$ is fairly
large, giving rise to edge state wavefunctions with a sizable finite component
even halfway through the lattice, however, as the size is increased and/or if
$K$ is decreased the bulk compoment of the state decreases. These results are
consistent with our previous derivation (Eq. (9)). For a system with
$J^{\prime}<J$ edge-states are not found. The bottom panels ((c) and (d)) show
the phase of the zero energy wave functions. The phases of both edge states
oscillate when going towards the center from the edges. The left and right
sides mirror each other, and a sudden change occurs at and within a few sites
near the center of the system.

In Figs. 4 and 5 the $d_{x},d_{y},d_{z}$ curves are shown for various values
of $p$ and $q$. The Hamiltonian parameters are as follows: left
$J^{\prime}=1.75,J=1.25$, center $J^{\prime}=1.5,J=1.5$, and right
$J^{\prime}=1.25,J=1.75$. $K=1.0$ for all cases. The plots are knot diagrams
Adams01 , two dimensional knot representations, arrived at by projecting the
curve onto the $d_{x}-d_{y}$ plane, but indicating which strand is below the
other at crossings. Other than the plots in the center in both figures,
neither of which are knots, the left and right plots show torus knots. For
example, in the top left figure of Fig. 5, the curve winds around $w_{1}=4$
times while it makes $w_{2}=5$ revolutions around the body of the torus.

Figure 5: Curves traced out by the parameters of the Hamiltonian
$d_{x},d_{y},d_{z}$ as $k$ traverses the extended Brillouin zone. The
upper(lower) set of plots show the case
$\frac{p}{q}=\frac{1}{5}$($\frac{p}{q}=\frac{2}{5}$). The parameters of the
Hamiltonian are, left: $J^{\prime}=1.75,J=1.25$,
center:$J=1.5,J^{\prime}=1.5$, and right:$J^{\prime}=1.25,J=1.75$ ($K=1$ in
all cases).

The center plots in both Figs. 4 and 5 show curves for the topological phase
transition. These curves do not form knots, they cross at the origin, which is
the gap closure point. It is obvious that these center plots “connect” the two
topologically distinct phases on either side of the phase transition. At the
phase transition the torus surface on which the curves “live” is a horn torus.
On either side of the horn torus, the curves reconnect, in a manner
reminiscent of “partner switching” of edge states in the Kane-Mele model
Kane05a ; Kane05b .

Figure 6: The trajectory of $\phi_{k}$ a parameter which defines the wave-
function (Eq. (21)) across the extended Brillouin zone for a trivial (a) and a
topological (b) case. As $k$ sweeps across the extended Brillouin zone
$\phi_{k}$ increases by $2\pi$ in the trivial, and $-4\pi$ in the topological
system.

The upper plots in Fig. 4 are for $p=1$, the lower ones are $p=2$. The upper
left plot shows a trefoil knot with $w_{1}=2$, $w_{2}=-3$. The upper right
plot shows an un-knot with $w_{1}=1$, $w_{2}=3$. Even though there are
crossings on this curve, they can be eliminated by type I Reidemeister moves
Adams01 , one of three modifications on two-dimensional representations of
knots which leave knot invariants unaltered. The knot invariants (in this case
the two winding numbers) are not altered by these modifications (if all three
twists are eliminated).

The topological phase transition connects the two states. In the insulating
phases Eq. (16) is recovered. In the lower set the $J^{\prime}>J$ curve
exhibits the unknot ($w_{1}=1,w_{2}=-3$), and the $J>J^{\prime}$ the trefoil
knot ($w_{1}=2,w_{2}=3$), again consistent with Eq. (16).

In Fig. 5 the $d_{x},d_{y},d_{z}$ curves are shown for various $q=5$ cases.
The upper plots are $p=1$, the lower ones are $p=2$. The phase transition in
the upper set is a topological one connecting a $w_{1}=4$ phase
($J^{\prime}>J$) with a $w_{1}=1$ phase ($J>J^{\prime}$), while the bottom set
shows a transition between a $w_{1}=3$ phase ($J^{\prime}>J$) and a $w_{1}=2$
phase ($J>J^{\prime}$), both accompanied by sign changes in $w_{2}$ whose
absolute value is five. The cases $p=3$ and $p=4$ (not shown) give winding
numbers of $w_{1}=2$, $w_{1}=3$ and $w_{1}=1$, $w_{1}=4$, respectively, again
with changes in sign in $w_{2}$ across the quantum phase transition. In all
cases, our polarization formula, Eq. (16), is recovered.

To summarize the information of Figs. 4 and 5, it is useful to connect the
results to what is known from the “usual” representation Su79 , in which the
Wannier functions corresponding to different basis sites are taken to have the
same phase (basis $I$). There, there are two distinct phases, one with winding
number zero (trivial phase), one with winding number one (topological phase),
separated by gap closure. Performing the unitary transformation in Eq. (14)
leads to a different representation of the same topological phases, and it is
these different representations which are shown in the figures. In a given
system with fixed $p/q$, the insulating phases are knots on ring or spindle
tori. A ring torus can be deformed into a spindle torus via a horn torus. The
horn torus occurs if $J=J^{\prime}$, that is when the gap closes, and where
the phase transition separating different topological phases occurs (center
panels in Figs. 4 and 5). Note that in the $p/q=1/3$ case, the
topological(trivial) phase corresponds to the spindle(ring) torus, whereas the
reverse is true for $p/q=2/3$.

When $K=0$ (the usual SSH model), the above curves remain on the $d_{x},d_{y}$
plane, hence one cannot speak of a torus. In the “usual” representation
Asboth15 ; Kane13  of this model the $J^{\prime}>J$ is is the topological
state with winding number of one, while the $J^{\prime}<J$ is the trivial site
with winding number zero, depending on whether the $d_{x},d_{y}$ curve
encircles the origin. Implementing distance dependence, turns each of those
phases into phases with distinct quantized polarizations. Gradually changing
the sign of $K$ changes the sign of both winding numbers for a given toroidal
knot (which goes into its chiral counterpart), leading to no change in the
polarization (since the polarization is the ratio of winding numbers). This is
consistent with the result derived in section III: the sign of $K$ does not
determine the presence of edge states (nor whether the state is topological or
trival), only the sign of the $J^{\prime}-J$ does.

Figure 7: Semilogarithmic plot of the absolute value of the difference between
the exact value of the polarization $P_{exact}$ and the polarization
calculated via two different methods. The solid like shows the use of an
irrational distance (inverse of the golden ratio), the diamonds show a
calculation based on a rational approximation to the inverse golden ratio.

Another way of seeing how the fractional polarization arises is to look at the
phase variable of the wave function. The wave functions that are the solution
of the Hamiltonian in Eq. (7) with the $d$-components given by Eq. (13) can be
written in the form

$\begin{pmatrix}\sin\left(\frac{\theta_{k}}{2}\right)\\\
\cos\left(\frac{\theta_{k}}{2}\right)\exp(i\phi_{k})\end{pmatrix},$ (21)

where $\theta_{k}$ and $\phi_{k}$ depend on $d_{x},d_{y},d_{z}$ (Eq. (13)).
The phase $\phi_{k}$ is shown in Fig. 6, for a system with $p=1$ and $q=3$. In
the trivial phase, the phase $\phi_{k}$ changes by $2\pi$ as $k$ covers the
extended Brillouin zone ($2\pi q=6\pi$). Accordingly, the polarization is
$\mathfrak{P}=-\frac{e}{2}\frac{1}{3}$. Meanwhile, in the topological phase,
the variable $\phi_{k}$ changes by $-4\pi$ along the Brillouin zone, and the
polarization is $\mathfrak{P}=\frac{e}{2}\frac{2}{3}$.

In Fig. 7 we present calculations for an irrational distance, namely the
inverse of the golden ratio. In one calculation, we used the irrational number
itself to derive a Bloch Hamiltonian (in place of $p/q$ in Eq. (13)), and
Brillouin zones of size $-\pi q,\pi q$ were used. In another calculation
ratios of the Fibonacci sequence were used to approximate the golden ratio,
and Eq. (20) was used. The results in Fig. 7 show that, while the polarization
based on using an irrational distance is reasonably close to the exact result,
using a rational approximation is a much more stable procedure, which
uniformly approaches the exact result.

If an adiabatic charge pumping experiment Atala13 ; Nakajima16  undergoing a
full cycle is considered the amount of charge pumped would be one unit of
charge per cell. This coincides with the results of Watanabe and Oshikawa
Watanabe18  who showed that the charge pumped is independent of how the
Hamiltonian is written. Fig. 8a shows the Fourier transform of $Z_{q}$,
defined as

$P_{x}=\sum_{q}\exp(i2\pi xq/L)Z_{q};\hskip 5.69046ptx=1,...,L;$ (22)

which corresponds to the underlying polarization distribution Souza00 ;
Resta99 ; Yahyavi17 ; Hetenyi19 . The case $p/q=1/3$ is considered. Shown are
two different parameter sets for each ordered state, as well as the phase
transition, where $J=J^{\prime}$. The distributions of the ordered states show
maxima corresponding to the values of the polarization calculated above. The
tuning of the parameters of the Hamiltonian within one phase change the shape
of the distribution (variance and other cumulants), but not the position of
the maximum. At the transition the distribution becomes nearly flat, the
system does not exhibit a well-defined polarization due to gap closure. For a
finite system the curve representing $J=J^{\prime}$ is not entirely flat,
because the $k-$space sampling does not precisely hit the points where the
gaps occur. However, as the thermodynamic limit is taken, the curve
progressively flattens out.

Figure 8: (a) Fourier transform of the polarization amplitude, or polarization
distribution shown for different cases of the system with $p/q=1/3$; two sets
in each ordered state and at the gap closure point. The curves are normalized
to one over the lattice of $L=100$. Legend is valid only for the upper panel.
(b) Reconstructed polarization distributions along an adiabatic charge pumping
process. For further explanation see the text.

In the cases treated heretofore the polarization can be easily determined by
the ratio of the winding numbers, hence the Berry phase formula is not
necessary. If other terms are added to the Hamiltonian, for example, an on-
site potential, then the Berry phase based formula in Eq. (20) is needed. In
the lower panel of Fig. 8b an adiabatic charge pumping process Nakajima16  is
shown. The Hamiltonian is reparametrized according to Eq. (10) and an
alternating on-site potential of strength $\Delta$ is also added. In the
calculations of Fig. 8b $t=K=1$. The plots shown are along total polarization
distribution curves along the circle in the $\delta t-\Delta$ plane. The
equation of the circle is $0.2\cos(\alpha),0.2\sin(\alpha)$, and the curves
labeled $A,B,C,D,E,F,G,H$ correspond to values of
$\alpha=0,\pi/4,\pi/2,3\pi/4,\pi,5\pi/4,3\pi/2,7\pi/4$ respectively. As
$\alpha$ increases, the maximum of the polarization distribution function
shifts towards the right within the unit cell. For the curves $G$ and $H$ the
maxima of the curves is to the left of the $A$-curve. This can be interpreted
as a charge which entered from a neighboring unit cell (to the left of the one
shown), as is expected to occur for adiabatic charge pumping. Note that at
points $A$ and $E$ the maxima correspond to the ones shown in Fig. 8a, and the
maximum of the charge distribution appears to travel smoothly in between. The
part of the adiabatic charge pump process between $\alpha=0$ and $\alpha=\pi$
occurs between systems which exhibit fractionally quantized polarization. If
the $I$ basis was used, the maximu would be located at the edge of and halfway
through the unit cell.

##  VI Conclusion

Lattices with a basis can be solved by various representations of the wave
function. The representations differ in the relative phases of the Wannier
functions associated with different sites within the unit cell. The different
representations are related by unitary transformations. The transformation can
be used to diagonalize the position operator, if t he Brillouin zone is
extended. When the nearest distance between basis sites is a rational number
$p/q$, the Brillouin zone has to be extended by a factor of $q$. When this
distance is an irrational number, the Brillouin zone extends to the entire
real line. This latter state of affairs raises the question of how to
calculate the polarization, to which we suggested the use of a rational
approximation to the distance between intra-unit cell sites.

Applying such a transformation to the Hamiltonian can give rise to auxiliary
topological features. In our example calculations, based on an extension of
the SSH model, the $d_{x},d_{y},d_{z}$-curve forms a toroidal knot, and the
ratio of the two winding numbers of the torus is proportional to the
polarization. Since the modification here is a mere unitary transformation of
the Hamiltonian, or a change in the relative phase of the Wannier functions
corresponding to the different sites within a unit cell, no new phases are
found, the phase diagram is unaltered. If the extended SSH model is studied in
the “usual” basis, where the phases of the Wannier functions are the same,
then there is no torus, the $d_{x},d_{y},d_{z}$ curve is an ellipse, and its
winding number around the $d_{z}$ axis defines the polarization (the distinct
phases have winding numbers zero or one). When the Hamiltonian is transformed
into the basis which diagonalizes the position, and the Brillouin zone is
extended, the $d_{x},d_{y},d_{z}$ curve becomes a knot on a torus. Within each
phase the torus on which the “curve” lives is of a different type, a ring or a
spindle torus. The quantum phase transition occurs when the torus is a horn
torus.

These conclusions fit into the recent findings of Bena and Montambaux Bena09
and those of Watanabe and Oshikawa Watanabe18 . The unitary transformation
which relates different representation choices can change quantities which
depend on the phase of the wave fuction. In Ref. Watanabe18  this is shown for
the current, in our case it is shown for the polarization (Berry-Zak phase),
and for the $d_{x},d_{y},d_{z}$ curves of the system. However, as argued in
both Ref. Bena09  and Watanabe18 , phase independent quantities Watanabe18
will not change. The phase independent quantity we investigated was the charge
pumped in an adiabatic process, and indeed it does not depend on the
representation used. The path of the polarization maximum during the process,
however, does depend on it.

As for experimental tests of our work, the model suggested here can be
realized in a setting of cold atoms trapped in optical lattices Atala13 ;
Nakajima16 . We anticipate that a partial change in the polarization, of the
type between maxima in Fig. 8 can be measured in experiments in which the
polarization is changed, as long as the process does not constitute a full
adiabatic cycle, for example, a polarization switch. One possible way to
realize the distance dependence explicitly is the following. For instance a
$p=1$, $q=3$ system could be realized by constructing first a tri-partite
periodic lattice, and approach our bi-partite system as a limit by decreasing
the coupling of one of the sites in each unit cell to zero.

## Acknowledgments

BH would like to thank J. K. Asbóth and B. Dóra for helpful discussions. BH
was supported by the National Research, Development and Innovation Fund of
Hungary within the Quantum Technology National Excellence Program (Project Nr.
2017-1.2.1-NKP-2017-00001) and by the BME-Nanotechnology FIKP grant (BME FIKP-
NAT).

## References

  * (1) C. Bena and G. Montambaux, New J. Phys., 11 095003 (2009).
  * (2) H. Watanabe and M. Oshikawa, Phys. Rev. X 8 021065 (2018).
  * (3) B. A. Bernevig and T. L. Hughes, Topological Insulators and Superconductors, Princeton University Press (2013).
  * (4) C. L. Kane, in Topological Insulators, vol. 6 Eds. M. Franz and L. Molenkamp, Chapter 1, Elsevier (2013).
  * (5) J. K. Asbóth, L. Oroszlány, and A. Pályi, A Short Course on Topological Insulators: Band Structure Topology and Edge states in One and Two Dimensions, Lecture Notes in Physics, vol. 919, Springer (2018).
  * (6) F. D. M. Haldane, Phys. Rev. Lett. 61 2015 (1988).
  * (7) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95 226801 (2005).
  * (8) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95 146802 (2005).
  * (9) W. P. Su, J. R. Schrieffer and A. J. Heeger, Phys. Rev. Lett. 42 1698 (1979).
  * (10) D. J. Thouless, M. Kohmoto, M. P. Nightingale and M. den Nijs Phys. Rev. Lett. 49 405 (1982).
  * (11) M. V. Berry, Proc. Roy. Soc. London A392 45 (1984).
  * (12) J. Zak, Phys. Rev. 62 2747 (1989).
  * (13) R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 47 R1651 (1993).
  * (14) R. Resta, Phys. Rev. Lett. 80 1800 (1998).
  * (15) R. Resta, J. Phys.: Cond. Mat. 12 R107 (2000).
  * (16) E. I. Blount, Solid State Physics 13 305 (1962).
  * (17) S. Kivelson, Phys. Rev. B 26 4269 (1982).
  * (18) N. W. Ashcroft and M. D. Mermin, Solid State Physics (Harcourt, 1976).
  * (19) L. Li, Z. Xu, and S. Chen, Phys. Rev. B 89 085111 (2014).
  * (20) W. Kohn, Phys. Rev. 133 A171 (1964).
  * (21) A. Altland and M. R. Zirnbauer, Phys. Rev. B 55 1142 (1997).
  * (22) A. P. Schnyder, S. Ryu, A. Furusaki, A. W. W. Ludwig, Phys. Rev. B 78 195125 (2008).
  * (23) R. Jackiw and C. Rebbi, Phys. Rev. D 13 3398 (1976).
  * (24) M. Nakahara, Geometry, Topology, and Physics, Graduate Student Series in Physics, CRC Press, 2nd Edition (2003).
  * (25) B. Hetényi and B. Dóra, Phys. Rev. B  99 085126 (2019).
  * (26) C. C. Adams, The Knot Book, American Mathematical Society, Providence, Rhode Island (2001).
  * (27) I. Souza, T. Wilkens, and R. M. Martin, Phys. Rev. B 62 1666 (2000).
  * (28) R. Resta and S. Sorella, Phys. Rev. Lett. 82 370 (1999).
  * (29) M. Yahyavi and B. Hetényi, Phys. Rev. A 95 062104 (2017).
  * (30) M. Atala, M. Aidelsburger, J. T. Barreiro, D. Abanin, T. Kitagawa, E. Demler, I. Bloch, Nat. Phys. 9 795 (2013).
  * (31) S. Nakajima, T. Tomita, S. Taie, T. Ichinose, H. Ozawa, L. Wang, M. Troyer, Y. Takahashi, Nat. Phys. 12 296 (2016).