# Bayesian Graph Neural Network for Fast identification of critical nodes in
Uncertain Complex Networks

Sai Munikoti, Student Member, IEEE, Laya Das, Balasubramaniam Natarajan,
Senior Member, IEEE This work has been submitted to the IEEE for possible
publication. Copyright may be transferred without notice, after which this
version may no longer be accessible. Preprint submitted to IEEE SMC 2021

###### Abstract

In the quest to improve efficiency, interdependence and complexity are
becoming defining characteristics of modern complex networks representing
engineered and natural systems. Graph theory is a widely used framework for
modeling such complex networks and to evaluate their robustness to
disruptions. Particularly, identification of critical nodes/links in a graph
can facilitate the enhancement of graph (system) robustness and characterize
crucial factors of system performance. Most existing methods of critical node
identification are based on an iterative approach that explores each node/link
of a graph. These methods suffer from high computational complexity and the
resulting analysis is network specific. Additionally, uncertainty associated
with the underlying graphical model further limits the potential value of
these traditional approaches. To overcome these challenges, we propose a
Bayesian graph neural network based node classification framework that is
computationally efficient and systematically incorporates uncertainties.
Instead of utilizing the observed graph for training the model, a MAP estimate
of the graph is computed based on the observed topology and node target
labels. Further, a Monte-Carlo (MC) dropout algorithm is incorporated to
account for the epistemic uncertainty. The fidelity and the gain in
computational complexity offered by the Bayesian framework is illustrated
using simulation results.

###### Index Terms:

Graph uncertainty, GNN, Graph robustness, Bayesian, critical nodes

##  I Introduction

Graph theory is a popular framework to study complex systems and is gaining
popularity with the growth of network size, complexity and data availability.
Hence, various graph specific machine learning algorithms are recently being
proposed in the literature for applications like node classification,
community detection, link prediction, etc. [1, 2]. The proper functioning of a
complex system depends on its constituent sub-systems (nodes) and their
interconnections (links). Usually, there exists a set of critical nodes that
are relatively more crucial to the overall operation of system. The removal of
critical nodes severely disrupts the network functioning. The network
functionality is often studied in the literature with the help of graph
robustness metrics such as effective graph resistance, flow robustness, etc.
[3]. Robustness quantifies the impact of loss of resources (nodes) on the
performance of a system. Due to the inherent topological structure, each node
contributes differently to graph robustness, i.e., its removal/loss affects
the robustness to a different degree. In this regard, the notion of node
criticality score is introduced, which quantifies the decrease in robustness
(such as effective graph resistance) when the corresponding node is removed
from the graph [4, 5]. Critical scores are then employed to classify the
nodes, with the high criticality assigned to the nodes whose removal maximally
decreases the graph robustness.

There are numerous methods proposed in the literature for identifying critical
nodes [6, 4, 7, 8]. However, there are four major drawbacks with these
conventional approaches of identifying critical nodes: (1) Computationally
expensive due to iterative nature of the algorithms; (2) The results obtained
from the iterative approaches are specific to the network considered and
cannot be generalized to other systems/graphs or alternate robustness metrics;
(3) Topological information is not explicitly used in node identification; (4)
Accuracy of the approximation methods decreases as the size of the network
grows. To overcome these challenges, we have proposed a graph neural network
based scalable framework in [5] that efficiently computes criticality scores
and consequently identifies critical nodes/links in large complex networks.
Although [5] can address some of the fundamental challenges of existing
approaches, it still overlooks one of the fundamental problems of graph-based
models, i.e., the uncertainty in the underlying graph (system). Most of the
current works on node criticality consider the observed graph as a perfect
model of the underlying complex system. However, in reality, this is not
always the case. For instance, a node (user) in a social network influences
another node with some probability. In protein-protein interaction networks,
the edges, representing interactions among proteins, are results of noisy and
error-prone measurements.

Uncertainties are present in the graph because of modeling assumptions,
insufficient data and noisy measurements. Uncertainties can lead to the
addition of spurious links/weights in the graph or miss the important links
between the nodes that have a strong relationship. Consequently, any
robustness analysis of the observed graph could lead to inaccurate inferences.
There are a few prior attempts to incorporate the uncertainties while
developing graphical models of the systems. For instance, authors in [9] have
used fuzzy set theory to determine the link weights. However, the approach in
[9] is not very generic. From a learning based perspective, authors in [10,
11] have extended the classic Graph Neural Network (GNN) to a Bayesian
framework. Here, the observed graph is viewed as a realization of a parametric
random graph model, and then the joint posterior estimate of graph parameters
as well as GNN weights is computed.

Inspired by this work, we propose a novel Bayesian inductive learner for graph
robustness (BILGR). The nodes are classified into three classes based on the
criticality scores. BILGR is trained on subset of nodes and then prediction
can be done for any node of the graph. The proposed node classification
approach is scalable, computationally efficient, and systematically accounts
for the uncertainties in link connections. To the best of the authors’
knowledge, this is the first attempt to develop a Bayesian framework of
critical node identification through GNN.

###  I-A Related work

The authors in [4] explore effective graph resistance as a measure of
robustness for complex networks. Robustness is improved by protecting the
node/link whose removal maximally increases the effective graph resistance.
The complexity of an exhaustive search to identify such a link is of order
$O(N^{5})$ for a graph of $N$ nodes. To overcome computational complexity, the
authors in [12], propose a method based on genetic algorithm to enhance
network robustness. Particularly, the authors focus on identifying links whose
removal would severely decrease the effective graph resistance of the graph.
However, the algorithm requires manual fine-tuning of various parameters and
is not scalable. Authors in [5] proposed a learning based framework for fast
identification of critical nodes/links in large graphs. However, the framework
does not handle graph uncertainties. Specifically for GNN, authors in [10,
11], leverage a Bayesian framework to systematically incorporate
uncertainties. The fundamental idea behind the proposed method is similar to
the Bayesian neural network where some prior distribution is assumed for the
unknown weights of the network and then the posterior is computed for a
different set of weight parameters sampled from Markov chain Monte Carlo
(MCMC) algorithm. From various experiments, the Bayesian GNN formulation shows
that, although computationally more demanding, it leads to (1) an ability to
learn more rapidly with less data, (2) a better capacity to represent
uncertainty, and (3) better robustness and resilience to noise or adversarial
attacks.

###  I-B Contributions

This work develops a new framework for estimating node criticality class while
systematically accounting for uncertainties. The key contributions of this
work are as follows:

  * •

A novel critical node identification framework is developed with a GNN-based
inductive learning approach that is applicable for graphs in the presence of
aleatoric and epistemic uncertainties.

  * •

The scalability of the proposed technique is demonstrated by training GNN
models on a small subset of nodes and predicting node classes in a relatively
larger fraction of the graphs that were not used for training.

  * •

Systematic integration of aleatoric and epistemic uncertainty through MAP
estimate of underlying graph and Monte-Carlo (MC) dropout algorithm,
respectively.

  * •

The superior performance is validated on power law cluster graphs with more
than $90$% identification accuracy and credible intervals. The computational
efficiency is demonstrated with execution time three orders smaller compared
to the conventional approaches.

The rest of the paper is organized as follows. Section II describes the
fundamentals of classic GNN followed by Bayesian GNN and graph robustness
metrics. Section III presents the proposed framework with experiments and
results in Section IV. Final conclusions are provided in Section V.

##  II Background

This section summarizes the working principle of classic GNN and Bayesian GNN,
which are then employed to devise our BILGR framework.

###  II-A Fundamentals of GNN

Recently, new generalizations and operations have been developed for extending
deep learning approaches to graph based data. [13] is one of the first works
that efficiently extended the two dimensional convolutional operations to
graph convolutional operations. In GNNs, the foremost task is to represent a
node by a vector such that it captures both structural as well as the desired
target information. These node embedding vectors can then be used for various
downstream applications like node classification, link prediction, by
appending the node embeddings with the classic feed forward layers. The
standard procedure to learn this embedding vector from initial node features
involves a “ message passing mechanism”. Here, the information (node feature)
is aggregated from the neighbors of a node and then combined with its own
feature to generate a new feature vector. This process is repeated to generate
the final embedding for each node of a graph.

Broadly, there are two methodologies to learn node embeddings. The first
category belongs to transductive based learning algorithms like deepwalk of
[14] and node2vec, where a complete graph is needed to learn the embeddings of
nodes. On the contrary, the second category belongs to an inductive based
learning algorithm, namely the GraphSAGE [14]. This algorithm learns the
mapping (aggregator) function instead of learning the embedding vectors.
Hence, GraphSAGE can induce the embedding of a new node or node unseen during
training, given its features and neighborhood. GraphSAGE learns a
representation for every node based on some combination of its neighboring
nodes, parametrized by $K$. The parameter $K$ controls the number of hops to
be considered as a neighbor. After defining the neighborhood, the aggregate
function is employed to collect neighbor’s embedding for a selected node. The
aggregate function could be as simple as a mean function or as complex as a
multilayer perceptron with trainable weights. For initialization, the node
feature vectors are assigned as embeddings. Thereafter, for each neighborhood
depth until $K$, a neighborhood embedding is generated with the aggregator
function for each node and concatenated with the existing node embedding.
Passing the concatenated vector through a feed forward multilayer layer
perceptron (MLP) layer results in the final node embeddings. GraphSAGE learns
the weights of an aggregator function and MLP by minimizing a relevant loss
function. In a supervised learning setup, the algorithm tries to reduce the
loss associated with node target labels. Once weights are learned, then an
embedding vector for any node can be predicted given the node features and its
neighboring information. Our framework is based on GraphSAGE as it is more
appropriate for large dynamic graphs.

###  II-B Fundamentals of Bayesian GNN

Bayesian GNN is an extension of GNN to account for uncertainties in the
underlying graph. It is possible to take a Bayesian perspective on GNN by
constructing a joint posterior distribution of the graph, the GNN weights and
the node labels. The marginal posterior probability of node labels is given
by,

$\begin{split}p(\bm{Z|Y_{L}},\bm{X},G_{obs})=\int p(\bm{Z|W,X},G_{obs}),\\\
p(\bm{W|Y_{L},X},G_{obs})p(G|G_{obs},\bm{X,Y_{L}})d\bm{W}dG,\end{split}$ (1)

where $\bm{W}$ denotes the random weights of a Bayesian GNN over graph $G$
[11]. As the integral in (1) cannot be computed analytically, MCMC
approximation is formed as,

$\begin{split}p(\bm{Z|Y_{L}},\bm{X},G_{obs})\approx\frac{1}{S}\sum_{s=1}^{S}\frac{1}{N_{G}}\sum_{i=1}^{N_{G}}\int
p(\bm{Z|W_{s,i},X},G_{obs})\end{split}$ (2)

In this approximation, $N_{G}$ graphs $G_{i}$ are sampled from
$p(G|G_{obs},\bm{X,Y_{L}})$ and $S$ weight samples $W_{s,i}$ are drawn from
$p(\bm{W|Y_{L},X,G_{i}})$ by training a Bayesian GNN corresponding to the
graph $G_{i}$. The final prediction result is an average of all the
predictions coming from different samples. Here, one can sample a number of
similar distributions but diverse topologies based on the posterior of the
graph generation model. Thus, in the aggregation step, the GNN model is able
to learn more general node embeddings by incorporating information from
different potential neighbors.

###  II-C Robustness metrics

Robustness of a graph signifies the resistance of a graph against node/link
failure. This is typically measured via the connectivity of the graph post
failure. Thus, higher robustness implies higher connectivity of the residual
graph post node/link failure. There are various network-based surrogate
metrics that approximately quantify the graph robustness against node failure
[3, munikoti2021robustness]. The authors in [3, 15] studied and compared
various graph theoretic metrics for approximating network robustness. The key
take away from [3, 15] is that no single metric can approximate the robustness
for all the different families of graph under all scenarios. Thus, different
metrics are used based on the study and type of graphs. Effective graph
resistance ($R_{g}$) appears to be most accurate and generic metric for graph
robustness [4]. Therefore, we focus on $R_{g}$ in this study. It is important
to note that the proposed approach can be extended to other surrogate
robustness metrics. $R_{g}$ is the sum of the effective resistances over all
pairs of vertices [16]. $R_{g}$ considers both the number of paths between the
nodes and their length (link weight), intuitively measuring the presence and
quality of back-up possibilities. $R_{g}$ can be computed in several ways and
equation (3) shows the spectral form.

$R_{g}=\frac{2}{N-1}\sum_{i=1}^{N-c}\frac{1}{\lambda_{i}},$ (3)

where $\lambda_{i}$ are the eigen values of the laplacian matrix of graph $G$
and $c$ is the number of connected components in the graph, which matches with
the number of zero eigen values. Thus, only non-zero eigen values are utilized
while computing $R_{g}$.

##  III Proposed approach

The proposed BILGR has two key features: (1) It incorporates data (i.e.,
graph) and model (i.e., GNN) uncertainty. (2) It poses identification of
critical nodes as node classification problem where classes denote the
criticality level of the nodes. The framework is described in the following
steps.

### Step 1- Train node embeddings

The embedding vector for nodes is generated by maximally utilizing all the
available information namely the observed graph $G_{obs}$, initial feature
vector $\bm{X}$, and targets $\bm{Y_{L}}$. Here, targets denote the node
criticality score/class. Algorithm 1 summarizes the steps involved in
computing node criticality scores [5]. The criticality scores are used to
classify nodes into three classes as shown in the Table I. High criticality
nodes correspond to class 1. Each node is assigned a feature vector $\bm{X}$,
consisting of node weighted degree as well as its average neighbor degree.

In contrast to [10], here we propose to learn the embeddings in a supervised
way by training a node classification model on $G_{obs}$ with node class
labels. The GNN based classification model consists of a GraphSAGE layer with
depth $K=3$. Once the GraphSAGE layer aggregates features from the neighbors,
a feed forward layer is connected to the output followed by a softmax
activation. The detail architecture is shown in figure 1. Thereafter, the
intermediate output of the trained model, i.e., the final output of graph
aggregation operation ( before the feedforward layer) would serve as node
embeddings.

TABLE I: Class labels of nodes based on criticality scores Criticality Score | class label
---|---
$0-0.3$ | Class 1
$0.3-0.7$ | Class 2
$0.7-1.0$ | Class 3

### Step 2-Compute node distance matrix

The node distance matrix represents the connectivity of nodes in an embedding
space. In other words, it determines the relative position of nodes by
computing the distances between them. The distance between any two arbitrary
nodes is [10],

$\begin{split}Z_{i,j}(\bm{X},G_{obs})=||\boldsymbol{e_{i}-e_{j}}||^{2},\end{split}$
(4)

where, the vectors $\boldsymbol{e_{i}}$ and $\boldsymbol{e_{j}}$ are the
embedding vector of node $i$ and $j$, respectively.

### Step 3-Obtain MAP estimate of graph

In real world modeling scenarios, the underlying graph may not be known
accurately. There are uncertainties associated with the link connections or
link weights. This type of uncertainty in an underlying graph belongs to the
category of an aleatoric uncertainty. Assuming the observed graph $G_{obs}$ is
a sample from a collection of random graph models, we attempt to a get a
posterior distribution of graph, i.e., $P(G|G_{obs},\bm{X,Y_{L}})$. This
posterior distribution is then plugged into equation ($1$) to generate final
node predictions. The analytical derivation of probability distribution is
intractable. To counter the intractability, one can design a suitable Markov
chain Monte Carlo (MCMC) simulation to sample from the posterior. However,
MCMC is challenging and computationally expensive for the case of graph,
especially when we are dealing with thousands and millions of nodes.
Therefore, a maximum a posteriori estimate (MAP) of the graph is obtained.
Here, the term estimate of a graph refers to estimate of the weighted
adjacency matrix of the graph. If the embedding vectors of the graph nodes are
assumed to be a signal residing on a graph, then one can use the framework of
learning a graph from the smooth signals. The smoothness criterion will of the
estimated graph. The graph $G_{est}$ (equivalently $A_{G_{est}}$) is computed
by solving the following graph learning optimization algorithm [10, 17],

$A_{G_{est}}=\operatorname*{argmin}_{A_{g}}||A_{G}\circ
Z||-\alpha\bm{1}^{T}\log(A_{G}\bm{1})+\beta||A_{G}||_{F}^{2}\\\ $ (5)

where, $A_{G}$ is the set of feasible symmetric matrices; $Z$ is the distance
matrix as computed in Step 2. The Hadamard product of $A_{G}$ and $Z$, i.e.,
$A_{G}\circ Z$, represents the Dirichlet energy of a graph, which quantifies
the smoothness of the signals on a weighted graph. The second term of the
optimization problem in equation (5) ensures that the node degree is never
zero although the individual link weights can be zero. The last term of the
equation constrains the adjacency matrix to be sparse. In a nutshell, the
optimization problem tries to learn a sparse weighted adjacency matrix that is
a better representation of given node embedding vectors.

#### Step 4-Train GNN over the inferred graph $G_{est}$ and employ MCMC to
sample different GNN weights

Once a graph $G_{est}$ is estimated, a new node classification model is
trained on $G_{est}$ with node class labels. This trained model is the final
model that would predict class labels for test nodes. However, here the GNN
model yields point predictions and does not account for epistemic uncertainty.
Therefore, a Monte Carlo dropout algorithm is implemented to account for
epistemic uncertainty [18]. In training time, we run multiple forward passes
with different dropout masks each time. During prediction, we collect
different predictions of test samples with different dropout masks. Then, mean
and variance of these predictions are used to devise the predictive posterior
mean and credible interval, respectively.

Algorithm 2 summarizes all the steps from 1-4. The schematic in figure 1
describes the complete Bayesian framework of GNN for critical node
identification. The details of model training and performance evaluation
experiments is discussed in the next section.

Figure 1: Proposed BILGR framework for critical node identification under
graph/model uncertainty

##  IV Experiments

###  IV-A Baseline approach

To evaluate the performance of our proposed approach, we compare the BILGR
with a conventional method of identifying node criticality. A typical
procedure involves removing a particular node from the graph and computing the
robustness metric of the residual graph. This process is repeated for all the
nodes of the graph. Thereafter, the metrics of all nodes are used to generate
node ranks and identify the high criticality nodes whose removal maximally
decreases the graph robustness. Algorithm $1$ summarizes the existing approach
of identifying high criticality nodes.

Algorithm 1 Conventional algorithm of critical node identification

0: Graph $G$ with $V$ nodes.

0: Node class based on criticality score LOOP Process

1: for $n$ in V do

2: Remove node $n$ from graph G

3: Compute robustness metric of the residual graph ($G-n$)

4: Assign criticality scores to node $n$

5: end for

6: return Assign class to nodes based on criticality scores (refer Table I)

Algorithm 2 Propose BILGR algorithm to identify critical nodes

0: Observed graph $G_{obs}$ with node feature vector $\bm{X}$ and subset node
class $\bm{Y_{L}}$.

0: Node class, i.e. node criticality.

1: Learn node embeddings in a supervised manner using $G_{obs}$, $\bm{X}$ and
$\bm{Y_{L}}$.

2: Compute node distance matrix using eqn. (4).

3: Solve optimization problem in eqn. (5) to obtain adjacency matrix of MAP
estimate of graph

4: Sample multiple networks weights using MC dropout masks by training model
on the estimated graph.

5: Predict multiple class probabilities for each test nodes using different
training dropout masks.

6: return mean and variance of predictions to get posterior mean and credible
interval of node classes.

###  IV-B Model settings and Training:

There are various hyper-parameters in the model that need to be tuned for
obtaining the best possible model. The number of node embedding layers, i.e.,
the depth of the GNN is selected as three. The number of neurons in these
three layers are $64$, $32$, and $16$ respectively. The classification module
consists of two feedforward layers with $10$ and $3$ neurons respectively. The
activation function in all the layers is kept as relu except last layer which
has softmax. The loss function is categorical cross entropy which is optimized
via ADAM optimizer. The training of GNN and Graph generation/manipulation are
carried out in TensorFlow and python NetworkX, respectively.

###  IV-C Results

The node identification problem is framed as node classification problem where
three classes denote the three levels of node criticality, i.e., high, medium
and low. However, in any graph, critical nodes is only $2\%$ to $5$ $\%$ of
the total nodes, thereby inducing a class imbalance problem. The imbalance is
addressed by employing weighted version of GNN also known as cost sensitive
learning. In the loss function, the weights are assigned to a class
proportional to the inverse of the class distribution present in the training
dataset. The weights are $100$, $50$ and $1$ for classes $1$,$2$ and $3$,
respectively. The performance of the models is reported with accuracy,
precision and recall. Here, recall is an important metric as we want to
correctly predict all class-1 nodes, i.e., high criticality nodes which are
usually small in number.

The model is trained on subset of $3000$ nodes in a power law cluster graph of
$5000$ nodes and $5997$ edges. Power law cluster is a family of graphs that
follow a power law degree distribution and also exhibit “clustering” by
requiring that, in some fraction of cases ($p$), a new node connects to a
random selection of the neighbors of the node to which the new node was
connected last time. This graph family is selected for illustration as many
real-world networks have shown to topologically resemble it. However, our
framework is generic for any family of the graph as well as any robustness
metric. Table II tabulates the scores of the trained model when evaluated on a
unseen $2000$ test nodes. As we are more interested in the class $1$
performance since they represent high criticality nodes, the scores are
tabulated for class-1 as well. It can be seen that the trained model is $90$ %
accurate in detecting all the classes of nodes. More importantly, the model is
correctly identifying $95$ % of the top critical nodes (i.e., $41$ out of
$43$) with high recall score of $0.95$.

TABLE II: Performance of the GNN model Model | Accuracy | Recall
---|---|---
Overall | 0.902 $\pm$ 0.009 | 0.900
Class-1 | 0.953 $\pm$ 0.01 | 0.950

Apart from high identification accuracy, the proposed BILGR framework offers
significant advantage in execution time compared to a conventional iterative
approach. Specifically, BILGR identifies $43$ high criticality nodes in $17$
seconds compared to $64600$ seconds with a conventional approach in a $5000$
nodes graph. Thus, the proposed method is multiple orders faster than the
conventional approach, and this gap will increase as the network size grows.
The training and experiments are conducted on a system with an Intel i7
processor running at 3.4GHz with 6GB Nvidia RTX 2070 GPU.

The Bayesian formulation of the proposed framework provides confidence in the
predictions which is useful when test nodes differs significantly from the
training nodes in terms of topological properties such as neighborhood sub-
graph, assortativity, etc. The confidence around the mean predictions of
accuracy is shown in the Table II. Further, to demonstrate the applicability
of confidence interval, the prediction of two class-1 test samples are shown
in the Figure 2. It can be observed that the model is more confident with test
sample $1$ (left hand side) with a tighter confidence interval compared to
test sample $2$. This type of extra inference is only possible in Bayesian
framework, which further increases the usability of the model.

Due to the systematic incorporation of graph and model uncertainty, the model
is also robust to noisy or adversarial graphs. To demonstrate this, noisy
versions of graph are generated by randomly adding links at $5$% of the nodes.
Table III lists the scores of the model when evaluated on two noisy graphs
(i.e., Gnoise1 and Gnoise2) with $6096$ and $6196$ links. The number of new
randomly added links in Gnoise1 and Gnoise2 are $100$ and $200$, respectively.
It can be observed that model prediction remains consistent to a very high
degree even though spurious links are added to the input graph. This
demonstrates the robustness of the BILGR framework to noisy or adversarial
test samples (nodes). Furthermore, with MC dropout algorithm we took $100$
predictions of each test node and reported the posterior mean and confidence
interval of the predictions in tables II and III. Compared to point
prediction, the accuracy with ensemble prediction has increased by approx $1$%
for each of the test scenarios.

Figure 2: Probability distributions of class probabilities for test samples of class-1  TABLE III: Robustness of Trained model to graph noise Test graph | Measure | Accuracy | Recall
---|---|---|---
Gnoise1 6096 links | Ovverall | 0.896 $\pm$ 0.09 | 0.866
Class-1 | 0.931 $\pm$ 0.01 | 0.930
Gnoise2 6196 links | Ovverall | 0.876$\pm$0.08 | 0.886
Class-1 | 0.890$\pm$0.02 | 0.89

As mentioned earlier, all the results are demonstrated for power law cluster
graph. A different model is needed for each family of synthetic graphs, i.e.,
power law, power law cluster, Stochastic block model, among others. This is
because, different families of graphs vary in their overall structure and link
connections, i.e., degree distributions, assortativity, average clustering
coefficient, etc. However, model trained on one family can be employed for
prediction of other graphs (even larger one) of the same family. Indeed,
models trained on synthetic graphs can be used for identification of critical
nodes in real networks [5]. This is possible due to the generalizability of
our BILGR framework ( inductive learning) and topological resemblance of real
networks with synthetic graphs.

##  V Conclusions

This work proposes a Graph neural network based BILGR framework to quickly
identify critical nodes under uncertain conditions. The criticality is related
to node robustness metrics such as effective graph resistance, weighted
robustness, among others. BILGR addresses data uncertainty by computing a MAP
estimate of a graph. Then, a GNN architecture based on node classification is
designed with node labels belonging to one of the three classes. The model is
trained on the estimated graph and tested on unseen nodes. To incorporate
model uncertainty, MC dropout is included in the GNN architecture. The model
has sufficient accuracy to identify critical nodes as well as it provides a
confidence interval around its predictions. The significant computational
advantage over classical approach of critical node identification is
demonstrated through algorithm running times. As a part of future work, the
BILGR will be extended to critical Link identification and other graph
families.

## Acknowledgement

This material is based upon work supported by National Science Foundation
under award number $1855216$.

## References

  * [1] Z. Lu, W. Lv, Y. Cao, Z. Xie, H. Peng, and B. Du, “Lstm variants meet graph neural networks for road speed prediction,” _Neurocomputing_ , vol. 400, pp. 34–45, 2020.
  * [2] P. Bongini, M. Bianchini, and F. Scarselli, “Molecular generative graph neural networks for drug discovery,” _Neurocomputing_ , 2021.
  * [3] M. J. Alenazi and J. P. Sterbenz, “Comprehensive comparison and accuracy of graph metrics in predicting network resilience,” in _2015 11th International Conference on the Design of Reliable Communication Networks (DRCN)_. IEEE, 2015, pp. 157–164.
  * [4] X. Wang, E. Pournaras, R. E. Kooij, and P. Van Mieghem, “Improving robustness of complex networks via the effective graph resistance,” _The European Physical Journal B_ , vol. 87, no. 9, p. 221, 2014.
  * [5] S. Munikoti, L. Das, and B. Natarajan, “Scalable graph neural network-based framework for identifying critical nodes and links in complex networks,” _arXiv preprint arXiv:2012.11831_ , 2020.
  * [6] V. Boginski and C. W. Commander, “Identifying critical nodes in protein-protein interaction networks,” in _Clustering challenges in biological networks_. World Scientific, 2009, pp. 153–167.
  * [7] X. Wang, Y. Koç, R. E. Kooij, and P. Van Mieghem, “A network approach for power grid robustness against cascading failures,” in _2015 7th international workshop on reliable networks design and modeling (RNDM)_. IEEE, 2015, pp. 208–214.
  * [8] P. Van Mieghem, D. Stevanović, F. Kuipers, C. Li, R. Van De Bovenkamp, D. Liu, and H. Wang, “Decreasing the spectral radius of a graph by link removals,” _Physical Review E_ , vol. 84, no. 1, p. 016101, 2011.
  * [9] C. Lam and K. Tai, “Modeling infrastructure interdependencies by integrating network and fuzzy set theory,” _International Journal of Critical Infrastructure Protection_ , vol. 22, pp. 51–61, 2018.
  * [10] Y. Zhang, S. Pal, M. Coates, and D. Ustebay, “Bayesian graph convolutional neural networks for semi-supervised classification,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 33, 2019, pp. 5829–5836.
  * [11] S. Pal, F. Regol, and M. Coates, “Bayesian graph convolutional neural networks using non-parametric graph learning,” _arXiv preprint arXiv:1910.12132_ , 2019.
  * [12] C. Pizzuti and A. Socievole, “A genetic algorithm for enhancing the robustness of complex networks through link protection,” in _International Conference on Complex Networks and their Applications_. Springer, 2018, pp. 807–819.
  * [13] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” _arXiv preprint arXiv:1609.02907_ , 2016.
  * [14] W. Hamilton, Z. Ying, and J. Leskovec, “Inductive representation learning on large graphs,” in _Advances in neural information processing systems_ , 2017, pp. 1024–1034.
  * [15] M. J. Alenazi and J. P. Sterbenz, “Evaluation and comparison of several graph robustness metrics to improve network resilience,” in _2015 7th International Workshop on Reliable Networks Design and Modeling (RNDM)_. IEEE, 2015, pp. 7–13.
  * [16] W. Ellens, F. Spieksma, P. Van Mieghem, A. Jamakovic, and R. Kooij, “Effective graph resistance,” _Linear algebra and its applications_ , vol. 435, no. 10, pp. 2491–2506, 2011.
  * [17] V. Kalofolias and N. Perraudin, “Large scale graph learning from smooth signals,” _arXiv preprint arXiv:1710.05654_ , 2017.
  * [18] Y. Gal and Z. Ghahramani, “Dropout as a bayesian approximation: Representing model uncertainty in deep learning,” in _international conference on machine learning_ , 2016, pp. 1050–1059.


# Topological metal and high-order Dirac point in cubic Rashba model

Haijiao Ji  Harbin Institute of Technology, Shenzhen, 518055, P. R. China
Ning Zhang  Harbin Institute of Technology, Shenzhen, 518055, P. R. China
Noah F. Q. Yuan<EMAIL_ADDRESS>Harbin Institute of Technology,
Shenzhen, 518055, P. R. China

###### Abstract

We investigate the properties of the two-dimensional model with Rashba-type
spin-orbit coupling cubic in electron momentum. In the normal phase, edge
states emerge on open boundaries. In the superconducting phase, edge states
could evolve into gapped fermionic edge states. Applications to realistic
materials of interface superconductors are also discussed.

Introduction.— The Rashba effect is one of the well-known manifestations of
spin-orbit couplings (SOCs) in solids where inversion symmetry is broken by an
electric field normal to the heterointerface [1, 2]. Another well-known type
of SOCs due to the inversion symmetry breaking is the so-called Dresselhaus
SOC [3]. In both cases, the SOC Hamiltonian is linear in electron momentum
$\bm{k}$, and the electron spins exhibit one full winding in moving around the
closed Fermi contour. Usually such $\bm{k}$-linear SOCs are leading order
effects in inversion symmetry breaking systems [4, 5, 6].

In certain systems the leading SOCs will not be $\bm{k}$-linear [7, 8, 9, 10,
11]. For example in transition metal dichalcogenides (TMDs) or the Kane-Mele
model, due to the special point group symmetry, the leading order SOC near
$\Gamma$ point is $\bm{k}$-cubic, while the electron spin is pinned to the
out-of-plane direction.

Other examples of nonlinear SOCs can be found in interfaces such as (111)
surface of GaAs quantum well [12, 13], (001) surface of oxide SrTiO3 (STO)
[14, 15, 16, 17, 18], and Si-terminated (001) surface of the rare-earth
antiferromagnet RRh2Si2 [19, 20, 21] (R denotes the rare-earth element such as
Yb, Tb, etc.), where the $\bm{k}$-cubic Rashba SOC can well describe the
dominant band structure, as verified by experiments [8, 19]. In GaAs quantum
well, fine tuning is needed to cancel $\bm{k}$-linear Rashba and Dresselhaus
effects [22]. When the dominant electrons are $d$\- (STO) or $f$-orbitals
(TbRh2Si2), and the leading order SOCs could become cubic [14, 21].

Superconductivity has been found experimentally in some of the above systems
with $\bm{k}$-cubic SOCs, including bulk and two-dimensional (2D) TMDs [23,
24, 25, 26, 27], LaAlO3/SrTiO3 (LAO/STO) interfaces [28, 29, 30], and RRh2Si2
[31]. Correspondingly, unconventional superconductivity such as Ising
superconductivity [32, 33, 34] and topological nodal superconductivity [35,
36, 37] has been found in 2D TMDs. It is therefore interesting to investigate
the properties of superconductors with $\bm{k}$-cubic SOCs such as LAO/STO
interfaces and RRh2Si2.

Figure 1: (a) Spin texture of the Fermi contours of the cubic model Eq. (1).
electron spins have three full windings in moving around the Fermi contours.
(b) The distribution of spin field $(S_{x},S_{y})$ in the Brillouin zone,
where the spin winding number $W=-3$ at $\Gamma=(0,0)$ point (colored blue),
while $W=+1$ at $X=(\pi,0)$, $Y=(0,\pi)$ and $M=(\pi,\pi)$ points (colored
red). The spin winding number $W$ is defined in Eq. (3) and the spin field
$(S_{x},S_{y})$ is defined underneath. (c) Phase diagram of the cubic Rashba
model Eq. (1), where edge states can be found explicitly in the pink region,
for example the black (d) and green (e) dots, while edge states are mixed with
bulk states in the white region such as the red dot [Fig. 4(a)]. (d) and (e)
are corresponding energy spectra with open boundary conditions along
$y$-direction with 200 lattice sites. (f) is the wavefunction of the edge
states at momentum indicated by black dagger in (e) with 50 lattice sites.

In this work, we focus on the abstract two-dimensional model equipped with the
cubic Rashba SOC, and study the superconducting phase. In this cubic Rashba
model, the SOC Hamiltonian is cubic in electron momentum $\bm{k}$, and the
electron spins exhibit three full windings in moving around the closed Fermi
contour, in contrast to the one full winding of the linear Rashba case. In the
linear Rashba model, the so-called topological metal has been proposed to
describe its normal phase [41, 38, 39, 40], and the corresponding
superconducting phase will also be affected to host finite-momentum
superconductivity. In the cubic Rashba model, the notion of topological metal
also applies to its normal phase, with edge states emerging on open edges. In
the superconducting phase, these edge states could evolve into gapped
fermionic edge states. Such exotic behaviors may have potential impacts on the
experimental measurements of superconductors such as LAO/STO interfaces and
RRh2Si2.

Figure 2: Zero contours $S_{x}=0$ (blue) and $S_{y}=0$ (red) of the spin field
$(S_{x},S_{y})$, and the intersecting points of blue and red lines are Dirac
points with corresponding spin winding numbers. (a) One cubic Dirac point at
zero field. (b) Three linear Dirac points under in-plane Zeeman field
$B_{x}=0.1\alpha$ along $x$-axis. (c) Five linear Dirac points under linear
Rashba $\lambda=0.1\alpha$. (d) Nine linear Dirac points under symmetry-
breaking perturbations.

Topological metal and edge states.— It is theoretically proposed and
experimentally verified that the following cubic Rashba model could describe
the low-energy band structure of interface superconductor LAO/STO [8],

${H}_{\rm
cR}(\mathbf{k})=\eta_{0}k^{2}+\eta_{1}k^{4}+\frac{1}{2}i\alpha(k_{+}^{3}\sigma_{-}-k_{-}^{3}\sigma_{+})$
(1)

with electron momentum $\mathbf{k}=(k_{x},k_{y})$, $k=|\mathbf{k}|$,
$k_{\pm}=k_{x}\pm ik_{y}$, $\sigma_{\pm}=\sigma_{x}\pm i\sigma_{y}$, and Pauli
matrices $\bm{\sigma}$ denoting spin. Here $\eta_{0,1}$ are coefficients for
kinetic terms, and $\alpha$ is the cubic Rashba SOC coefficient.

The unitary symmetries of the cubic Rashba model form the point group $D_{2d}$
with three generators: Twofold in-plane rotation $C_{2z}$, twofold rotation
$C_{2x^{\prime}}$ along the diagonal direction
$\bm{x}^{\prime}\equiv(\hat{\bm{x}}+\hat{\bm{y}})/\sqrt{2}$, and in-plane
mirror reflection $M_{x}:(k_{x},k_{y})\to(-k_{x},k_{y})$,
$(\sigma_{x},\sigma_{y},\sigma_{z})\to(\sigma_{x},-\sigma_{y},-\sigma_{z})$.
Though the point group is anisotropic, the energy bands of the cubic Rashba
model are isotropic

$E_{\pm}(\mathbf{k})=\eta_{0}k^{2}+\eta_{1}k^{4}\pm\alpha k^{3},$ (2)

and the two Fermi contours are hence circular and concentric. The inner $(+)$
and outer $(-)$ Fermi contours correspond to the upper $(+)$ and lower $(-)$
bands.

Due to the cubic Rashba SOC, electron spins have three full windings in moving
around the Fermi contours as shown in Fig. 1(a), which can be characterized by
the Berry phase $\pi$ of both Fermi contours as enforced by the parity
symmetries including time-reversal symmetry $\mathcal{T}$, mirror symmetries
$M_{x}$, $M_{y}$ and combined symmetry $C_{2z}\mathcal{T}$. The metallic phase
with quantized nonzero Berry phase $\pi$ of the Fermi contours is known as the
topological metal [41, 38, 39, 40], which is consistent with quantum
oscillation data of LAO/STO interfaces [17].

In the topological metal, there are odd numbers of Dirac points enclosed by
the Fermi contours. However, the detailed properties of the Dirac point are
not captured by the Berry phase, hence the properties of the topological metal
is yet to be revealed.

For this purpose, we introduce the spin winding number to describe the
topological property of the Dirac point

$\displaystyle W=\frac{1}{2\pi}{\rm Im}\oint_{\rm
FS}\frac{dS_{x}+idS_{y}}{S_{x}+iS_{y}},$ (3)

where FS denotes the Fermi surface (or Fermi contour in 2D) and
$\bm{S}=(S_{x},S_{y},S_{z})=\langle\bm{\sigma}\rangle$ is the spin expectation
value of the lower band $E_{-}$. The spin winding number of the cubic Dirac
point in Eq. (1) can be worked out as $W=3$. For linear Dirac point
$H=\sum_{ij}k_{i}A_{ij}\sigma_{j}$, the spin winding number is $W={\rm
sgn}(\det A)$.

The cubic Rashba model Eq. (1) is defined near the $\Gamma$ point. When this
model is extended to the entire Brillouin zone (BZ), the Poincaré-Hopf index
theorem dictates that the total spin winding number is zero

$\sum_{i}W_{i}=0.$ (4)

The BZ compatible with point group $D_{2d}$ would be that of a square lattice.
Without loss of generality, we assume the lattice constant of the square
lattice to be 1. When the time-reversal-invariant points $X=(\pi,0)$,
$Y=(0,\pi)$ and $M=(\pi,\pi)$ are linear Dirac points $|W_{X,Y,M}|=1$ while
$\Gamma$ is cubic $W_{\Gamma}=3$, the only choice of spin winding numbers is
$W_{X,Y,M}=-1$ due to the Poincaré-Hopf index theorem, which is shown in the
distribution of spin field $(S_{x},S_{y})$ in the BZ in Fig. 1(b). Details of
the tight-binding cubic Rashba model of Eq. (1) can be found in the Appendix.

The nontrivial spin winding number of the cubic Dirac point indicates the
nontrivial bulk topology and hence topological edge states in the tight-
binding model. In the phase diagram of Fig. 1(c), when the parameters
$\eta_{0,1}/\alpha$ fall in the topological region (colored pink), the two
energy bands are separated by local gaps as shown in Fig. 1(d) and (e). Under
the open boundary condition of $x$ or $y$-direction, edge states emerge, whose
dispersions are within local gaps [red lines in Fig. 1(d) and (e)], and whose
wavefunctions are localized on two edges with opposite spin polarization along
$z$-axis as shown in Fig. 1(f). Notice that the edge states are spin
degenerate.

To figure out the origin of such edge states, we can consider the chiral limit
$\eta_{0}=\eta_{1}=0$ as shown in Fig. 1(d). In this chiral limit, the cubic
model Eq. (1) has the chiral symmetry $\\{H_{\rm cR},\sigma_{z}\\}=0$ and
hence hosts zero energy modes on its open boundaries. Classified in the BDI
class, nodal points of the bulk energy spectrum are characterized by the spin
winding number $W$, and under open boundary condition, flat bands will emerge
between projected nodal points with opposite spin winding numbers. In
particular, when $y$-direction is open, the nodal points are projected to
$k_{x}=0$ with net spin winding number $W(k_{x}=0)=+3-1=+2$, and
$k_{x}=\pm\pi$ with net spin winding number $W(k_{x}=\pm\pi)=-1-1=-2$. Thus a
doubly degenerate flat band emerges connecting $k_{x}=0$ and $k_{x}=\pi$, and
another doubly degenerate flat band emerges connecting $k_{x}=0$ and
$k_{x}=-\pi\equiv\pi({\rm mod}\ 2\pi)$. Furthermore, zero-energy edge states
are eigenstates of chiral symmetry $\sigma_{z}$, and are hence spin polarized
along $z$-axis. As shown in Fig. 1(f), we hence find spin-degenerate edge
states.

Away from the chiral limit, one may treat $\eta_{0,1}$ as perturbations, under
which the nodal points evolve into Dirac points connecting upper and lower
bands, and the edge states become dispersive. As long as the local gaps are
not closed, the corresponding edge states are protected by the bulk topology
and separated from bulk states. As a result, we obtain phase diagram of Fig.
1(c).

Topological metal under perturbations.— We consider the following
perturbations in the metallic phase

${H}_{\rm
per}(\mathbf{k})=\frac{1}{2}i\lambda(k_{+}\sigma_{+}-k_{-}\sigma_{-})+\bm{B}\cdot\bm{\sigma}$
(5)

where $\lambda$ is the linear Rashba coefficient and $\bm{B}$ is the Zeeman
field. These perturbations are allowed under point group $D_{2d}$, under which
the number and positions of Dirac points of the band structure will change.

Under an in-plane Zeeman field, the single cubic Dirac point in Fig. 2(a) is
split into three linear Dirac points as shown in Fig. 2(b), forming a
threefold structure. Under the linear Rashba effect, the cubic Dirac point is
split into five linear Dirac points as shown in Fig. 2(c), forming a fourfold
structure. During the Dirac point splitting process, the total spin winding
number in the neighborhood of the $\Gamma$ point is conserved due to topology
[42], as shown in Fig. 2(b) and (c). In fact, the two energy bands can be
expressed in terms of the polar coordinates of the electron momentum
$\mathbf{k}=k(\cos\theta,\sin\theta)$ as

$\displaystyle E_{\pm}(\mathbf{k})=\eta_{0}k^{2}+\eta_{1}k^{4}\pm
D(k,\theta),$ (6)

where the energy bands are threefold symmetric under an in-plane Zeeman field
$\bm{B}=B(\cos\phi,\sin\phi)$, while fourfold symmetric under linear Rashba
effect $\lambda$:

$D(k,\theta)=\begin{cases}\sqrt{B^{2}+\alpha^{2}k^{6}-2B\alpha
k^{3}\sin(3\theta-\phi)}&\lambda=0,\\\
k\sqrt{\lambda^{2}+\alpha^{2}k^{4}-2\lambda\alpha k^{2}\cos
4\theta}&B=0.\end{cases}$ (7)

With both in-plane Zeeman field and linear Rashba effect, the in-plane
rotation symmetry of the energy bands is then lost.

We introduce the multiplicity (also known as the Milnor number) $\mu$ of the
Dirac point, which is defined as the maximal number of split Dirac points
under perturbations [42]. Here perturbations are defined as the Hamiltonian
terms with lower orders in momentum than that of the Dirac point. The cubic
Dirac point can be split into at most five linear Dirac points under the
symmetry-allowed perturbations in Eq. (5), hence the symmetry-constraint
multiplicity of the cubic Dirac point is $\mu_{\rm sym}=5$. Without symmetry
constraints, the cubic Dirac point can be split into at most $\mu=W^{2}=9$
linear Dirac points, as shown in Fig. 2(d) where the time-reversal symmetry
$\mathcal{T}$ is broken and the point group is reduced to
$C_{1v}=\\{1,M_{y}\\}$.

We have discussed the effects of perturbations on the bulk spectrum especially
the Dirac points. Due to the bulk-edge correspondence, the edge states will
also be affected. The general mechanism indicates that edge states will emerge
between projected nodal points with opposite spin winding numbers. For
simplicity we consider the chiral limit first. As shown in Fig. 3, edge states
are found under open boundary conditions, which are flat bands in the chiral
limit. Notice that under in-plane Zeeman field (a, b) and linear Rashba effect
(c), the net spin winding numbers of the nodal points near projected $\Gamma$
point can become 0, $\pm 1$ besides $\pm 2$. Notice that under an out-of-plane
field $B_{z}$ (d), the spin-degenerate edge band will get split into two spin
polarized bands, whose wavefunctions are similar to those shown in Fig. 1(f).

Figure 3: Energy spectra near $\Gamma$ point in the chiral limit under
external Zeeman fields (a,b,d) and linear Rashba effect (c), where edge states
are in red color. Parameters are the same as Fig. 2, and $B_{z}=0.1\alpha$ in
(d).

The introduction of kinetic terms $\eta_{0,1}$ will gradually close the local
gaps, which disperse and eventually diminish edge states. In the following, we
would like to show that even when kinetic terms close the local gaps and
diminish edge states in the normal phase, fermionic edge states will get
revealed under appropriate pairing potentials in the superconducting phase.

Topological metal under pairings.— We consider the cubic Rashba model with
pairing potential, whose Bogouliubov-de Gennes (BdG) Hamiltonian reads

${H}_{\rm BdG}(\mathbf{k})=\left(\begin{array}[]{cc}{H}_{\rm
cR}(\mathbf{k})-\varepsilon_{F}&H_{\Delta}(\mathbf{k})\\\
H_{\Delta}^{\dagger}(\mathbf{k})&-{H}^{*}_{\rm
cR}(-\mathbf{k})+\varepsilon_{F}\end{array}\right),$ (8)

where $\varepsilon_{F}$ is the Fermi energy and the generic pairing is

$\displaystyle H_{\Delta}=(\psi+\bm{d}\cdot\bm{\sigma})i\sigma_{y},$ (9)

with the spin-singlet $\psi(\mathbf{k})$ and spin-triplet $\bm{d}(\mathbf{k})$
pairings. Due to the Pauli exclusion principle, the spin-singlet is even in
momentum $\psi(-\mathbf{k})=\psi(\mathbf{k})$ and spin-triplet is odd
$\bm{d}(-\mathbf{k})=-\bm{d}(\mathbf{k})$. According to the symmetry group
$D_{2d}$, the possible spin-singlet and spin-triplet pairings are calculated
as shown in Table. 1.

Under pairings, edge states will be revealed, as long as

$\displaystyle\Delta_{+}<\Delta_{-}$ (10)

where $\Delta_{\pm}$ denote the pairing gaps of the inner $(+)$ and outer
$(-)$ Fermi contours respectively. For example, one can consider the $s$-wave
pairing with constant $(\Delta_{0})$ and extended $(\Delta_{1})$ parts
$\psi=\Delta_{0}+\Delta_{1}k^{2}$. As shown in Fig. 4(a), when $\Delta_{1}>0$
and the condition Eq. (10) is satisfied, the edge states (colored red) will
emerge joining the inner $(+)$ and outer $(-)$ Fermi contours with pairing
gaps $\Delta_{+}$ and $\Delta_{-}$ respectively, whose wavefunction is also
localized on two open edges with opposite spin polarization along $z$-axis as
shown in Fig. 4(b), where electron spin $\bm{\sigma}$ is equivalent to hole
spin $-\bm{\sigma}^{*}$. Importantly, the condition Eq. (10) applies when bulk
and edge states are separated [Fig. 1(c) pink region] or mixed [Fig. 1(c)
white region].

Figure 4: (a) Superconducting spectrum of the BdG Hamiltonian Eq. (8) with $\psi=\Delta_{0}+\Delta_{1}k^{2}$ and $\bm{d}=\bm{0}$. The red lines denote the gapped fermionic edge states. Parameters are $(\eta_{0},\eta_{1})=(1,\frac{1}{2})\alpha$ [red dot in Fig. 1(c)] and $(\Delta_{0},\Delta_{1})=(0,\frac{1}{4})\alpha$. (b) Edge state wavefunctions at the black dagger in (a), with electron/hole (e/h) and spin ($\uparrow/\downarrow$) components. Inset: Logarithm of edge state wavefunctions, showing the exponential localization of edge states. Table 1: Possible pairings in the cubic Rashba model. $D_{2d}$ | $\psi$ | $\bm{d}$
---|---|---
$A_{1}$ | $1$ | $k_{y}\hat{\bm{x}}+k_{x}\hat{\bm{y}}$
$k^{2}$ | $k_{x}^{2}k_{y}\hat{\bm{x}}+k_{y}^{2}k_{x}\hat{\bm{y}}$
$k^{4},{\rm Re}k_{+}^{4}$ | $k_{y}^{3}\hat{\bm{x}}+k_{x}^{3}\hat{\bm{y}}$
$A_{2}$ | ${\rm Im}k_{+}^{4}$ | $k_{x}\hat{\bm{x}}-k_{y}\hat{\bm{y}}$
$k_{x}k_{y}^{2}\hat{\bm{x}}-k_{y}k_{x}^{2}\hat{\bm{y}}$
$k_{x}^{3}\hat{\bm{x}}-k_{y}^{3}\hat{\bm{y}}$
$B_{1}$ | ${\rm Im}k_{+}^{2}$ | $\bm{k}$
$k_{x}k_{y}^{2}\hat{\bm{x}}+k_{y}k_{x}^{2}\hat{\bm{y}}$
$k_{x}^{3}\hat{\bm{x}}+k_{y}^{3}\hat{\bm{y}}$
$B_{2}$ | ${\rm Re}k_{+}^{2}$ | $\hat{\bm{z}}\times\bm{k}$
$k_{x}^{2}k_{y}\hat{\bm{x}}-k_{y}^{2}k_{x}\hat{\bm{y}}$
$k_{y}^{3}\hat{\bm{x}}-k_{x}^{3}\hat{\bm{y}}$
$E$ | 0 |  $\left(\begin{array}[]{c}k_{x}\hat{\bm{z}}\\\ k_{y}\hat{\bm{z}}\end{array}\right)$, $\left(\begin{array}[]{c}k_{x}k_{y}^{2}\hat{\bm{z}}\\\ k_{x}^{2}k_{y}\hat{\bm{z}}\end{array}\right)$, $\left(\begin{array}[]{c}k_{x}^{3}\hat{\bm{z}}\\\ k_{y}^{3}\hat{\bm{z}}\end{array}\right)$

Discussions.— As experimental candidates, we mainly consider LAO/STO interface
[8]. The conduction band of bulk SrTiO3 originates from Ti 3$d$ orbitals,
whose bottom is composed of $t_{2g}$ states ($d_{xy},d_{yz}$ and $d_{xz}$-like
orbitals) due to crystal field splitting [43]. Including the atomic SOC, the
lowest energy states for bulk SrTiO3 consist of fourfold degenerate bands,
which correspond to total angular momentum $J=3/2$ labeled by $z$-component
$J_{z}=\pm 1/2$ and $\pm 3/2$. Due to the interface potential induced by
LaAlO3, it turns out that $d_{xy}$-like orbitals ($J_{z}=\pm 3/2$) dominate
the bottom conduction bands of LAO/STO interface. Under point group $D_{2d}$
on the basis of $J_{z}=\pm 3/2$ states, the cubic Rashba model Eq. (1) could
be derived, and the cubic Rashba coefficient $\alpha$ in Eq. (1) can be tuned
by the electric gate voltage [8].

The edge states and bulk states can mix at the same energy, which prevents
differentiating bulk and edge states in the topological metallic phase. One
possible way is to apply an out-of-plane magnetic field to create Landau
levels, so that bulk states become gapped. Another way is to consider the
superconducting phase under magnetic field. In that case, the edge states and
bulk states give rise to different spatial profiles of supercurrent density.
To be more precise, one may prepare a Josephson junction in terms of LAO/STO,
and measure the Fraunhofer pattern under an out-of-plane field. By appropriate
transforms of the Fraunhofer pattern [44], the supercurrent density
distribution can be obtained, where the edge localized supercurrent
corresponds to edge states [45].

##  I Acknowledgements

We thank Dong-Hui Xu for bringing us the cubic Rashba model. We thank K. T.
Law, B. Q. Lv and Yuanfeng Xu for helpful discussions. This work is supported
by the National Natural Science Foundation of China (Grant. No. 12174021).

## Appendix: Tight-binding cubic Rashba model

The continuum cubic Rashba model reads

${H}_{\rm
cR}(\mathbf{k})=\eta_{0}k^{2}+\eta_{1}k^{4}+\frac{1}{2}i\alpha(k_{+}^{3}\sigma_{-}-k_{-}^{3}\sigma_{+}),$
(11)

where

$k^{2}=k_{x}^{2}+k_{y}^{2},\quad k^{4}=(k_{x}^{2}+k_{y}^{2})^{2},$ (12)
${k}_{+}^{3}=(k_{x}+ik_{y})^{3}=k_{x}^{3}+3ik_{x}^{2}k_{y}-3k_{x}k_{y}^{2}-ik_{y}^{3},$
(13)
${k}_{-}^{3}=(k_{x}-ik_{y})^{3}=k_{x}^{3}-3ik_{x}^{2}k_{y}-3k_{x}k_{y}^{2}+ik_{y}^{3}.$
(14)

Via the following substitution

$\begin{split}&k_{x(y)}\to\sin k_{x(y)},\\\ &k_{x(y)}^{2}\to 2(1-\cos
k_{x(y)}),\\\ &k_{x(y)}^{3}\to 2\sin k_{x(y)}(1-\cos k_{x(y)}).\\\
\end{split}$ (15)

The tight-binding Hamiltonian becomes

$\begin{split}H&=-2\eta_{0}(\cos{k_{x}}+\cos k_{y}-2)+4\eta_{1}(\cos
k_{x}+\cos k_{y}-2)^{2}\\\ &+\alpha(3\cos 2k_{y}-4\sin k_{x}-\sin
2k_{x})\sigma_{y}\\\ &-2\alpha\sin k_{y}(-3\cos k_{x}+\cos
k_{y}+2)\sigma_{x}.\end{split}$ (16)

## References

  * [1] E. I. Rashba and V. I. Sheka, Fiz. Tverd. Tela: Collected Papers 2, 162 (1959).
  * [2] Y. A. Bychkov and E. I. Rashba, J. Phys. C 17, 6039 (1984).
  * [3] G. Dresselhaus, Phys. Rev. 100, 580 (1955).
  * [4] J. Nitta, T. Akazaki, H. Takayanagi and T. Enoki, Phys. Rev. Lett. 78, 1335 (1997).
  * [5] A.Manchon, H. C. Koo, J. Nitta, S.M. Frolov and R. A. Duine, Nat. Mater. 14, 871 (2015).
  * [6] G. Bihlmayer, O. Rader, and R. Winkler, New J. Phys. 17, 050202 (2015).
  * [7] R. Winkler, Phys. Rev. B, 62, 4245 (2000).
  * [8] H. Nakamura, T. Koga, and T. Kimura, Phys. Rev. Lett. 108, 206601 (2012).
  * [9] H. Liu, E. Marcellina, A. R. Hamilton, and D. Culcer, Phys. Rev. Lett. 121, 087701 (2018).
  * [10] H. J. Zhao, et al. Phys. Rev. Lett. 125, 216405 (2020).
  * [11] D. Y. Usachov, et al. Phys. Rev. B 101, 245140 (2020).
  * [12] T.Koga, J. Nitta, T. Akazaki, and H. Takayanagi, Phys. Rev. Lett. 89, 046801 (2002).
  * [13] G. M. Minkov, A. A. Sherstobitov, A. V. Germanenko, O. E. Rut, V. A. Larionova, and B. N. Zvonkov, Phys. Rev. B 71, 165312 (2005).
  * [14] A. D. Caviglia, M. Gabay, S. Gariglio, N. Reyren, C. Cancellieri, and J. M. Triscone, Phys. Rev. Lett. 104, 126803 (2010).
  * [15] P. D. C. King, et al., Nat. Commun. 5, 3414 (2014).
  * [16] H. Liang, et al., Phys. Rev. B 92, 075309 (2015).
  * [17] M. J. Veit, et al., Nat. Commun. 9, 1458 (2018).
  * [18] W. Lin, et al., Nat. Commun. 10, 3052 (2019).
  * [19] D. Y. Usachov, et al. Phys. Rev. Lett. 124, 237202 (2020).
  * [20] I. A. Nechaev, and E. E. Krasovskii, Phys. Rev. B 98, 245415 (2018).
  * [21] S. Schulz, et al., Phys. Rev. B 103, 035123 (2021).
  * [22] S. D. Ganichev, and L. E. Golub, Phys. Status Solidi, 251, 1801 (2014).
  * [23] X. Qian, J. Liu, L. Fu, and J. Li, Science, 346, 1344 (2014).
  * [24] E. Sajadi, et al., Science 362, 922 (2018)
  * [25] W. Wang, S. Kim, M. Liu, F. A. Cevallos, R. J. Cava, and N. P. Ong, Science, 368, 534 (2020) .
  * [26] S. Kezilebieke, et al., Nature, 588, 424 (2020).
  * [27] E. Zhang, et al., Nat. Phys. 19, 106 (2023).
  * [28] L. Li, C. Richter, J. Mannhart, and R. C. Ashoori, Nat. Phys. 7, 762 (2011).
  * [29] J. A. Bert, et al., Nature physics, 7, 767 (2011).
  * [30] Y. L. Han, et al., Appl. Phys. Lett. 105, 19 (2014).
  * [31] E. Schuberth, et al., Science, 351, 485 (2016).
  * [32] J. M. Lu, O. Zheliuk, I. Leermakers, N. F. Yuan, U. Zeitler, K. T. Law, and J. Ye, Science, 350, 1353 (2015).
  * [33] B. T. Zhou, N. F. Yuan, H. L. Jiang, and K. T. Law, Phys. Rev. B, 93, 180501 (2016).
  * [34] S. C. De la Barrera, et al., Nat. Commun. 9, 1427 (2018).
  * [35] W. Y. He, B. T. Zhou, J. J. He, N. F. Yuan, T. Zhang, and K. T. Law, Commun. Phys. 1, 40 (2018).
  * [36] A. K. Nayak, et al., Nat. Phys. 17, 1413 (2021).
  * [37] A. Hamill, et al., Nat. Phys. 17, 949 (2021).
  * [38] A. Alexandradinata, C. Wang, W. Duan, L. Glazman, Phys. Rev. X 8, 011027 (2018).
  * [39] C. Wang, W. Duan, L. Glazman, A. Alexandradinata, Phys. Rev. B 100, 014442 (2019).
  * [40] S. Sun, Z. Song, H. Weng, X. Dai, Phys. Rev. B 101, 125118 (2020).
  * [41] N. F. Yuan, and L. Fu, Proc. Nat. Acad. Sci. 118, e2019063118 (2021).
  * [42] N. F. Yuan, and L. Fu, Phys. Rev. B 101, 125120 (2020).
  * [43] Z. S. Popovic, S. Satpathy, and R. M. Martin, Phys. Rev. Lett. 101, 256801 (2008).
  * [44] R. C. Dynes and T. A. Fulton, Phys. Rev. B 3, 3015 (1971).
  * [45] W. Wang, et al. Science 368, 534 (2020).