11institutetext: Timoteo Carletti 22institutetext: naXys, Namur Institute for
Complex Systems, University of Namur, rue de Bruxelles 61, B5000 Namur,
Belgium, 22email<EMAIL_ADDRESS>33institutetext: Duccio Fanelli
44institutetext: Università degli Studi di Firenze, Dipartimento di Fisica e
Astronomia, CSDC and INFN, via G. Sansone 1, 50019 Sesto Fiorentino, Italy
44email<EMAIL_ADDRESS>
# Dynamical systems on hypergraphs

Timoteo Carletti and Duccio Fanelli

###### Abstract

We present a general framework that enables one to model high-order
interaction among entangled dynamical systems, via hypergraphs. Several
relevant processes can be ideally traced back to the proposed scheme. We shall
here solely elaborate on the conditions that seed the spontaneous emergence of
patterns, spatially heterogeneous solutions resulting from the many-body
interaction between fundamental units. In particular we will focus, on two
relevant settings. First, we will assume long-ranged mean field interactions
between populations, and then turn to considering diffusive-like couplings.
Two applications are presented, respectively to a generalised Volterra system
and the Brusselator model.

##  1 Introduction

The study of many body interactions has a long history in science and
technology, and relevant results have been obtained under the assumption of
regularity of the underlying substrates, where the dynamics eventually
develops. When regularity gets lost, general results are scarce and
simplifying assumptions, which implement dedicated approximations, need to be
put forward. It is for instance customary to reduce the many body exchanges
within a pool of simultaneously interacting entities to a vast collection of
pairwise contacts, a working ansatz which drastically reduces the intimate
complexity of the scrutinised dynamics. Governing dynamical systems are hence
cast on top of networks AlbertBarabasi ; BLMCH  with diverse and variegated
topologies: each node contains a replica of the original system, and the
strength of interaction is set by the weight of the associated link.

Despite this crude approximation, relevant results have been obtained which
bear general interest Newmanbook ; Barabasibook ; Latorabook . At the same
time many examples of systems exist for which the above assumption holds true
just as a first order approximation BGL ; LRS2019 . To overcome this intrinsic
limitation, the effect of aggregated structures of nodes, such as cliques,
modules or communities Newmanbook ; fortunatohric2016  has been recently
addressed in the literature. This implies analysing the cooperative
interference within bunches of tightly connected nodes and assessing their
role in shaping the ensuing dynamics, in the framerwok of a generalized
picture which accounts for multiple pairwise exchanges.

There are however several examples where the interactions among individuals,
being them neurons petri2014homological ; LEFGHVD2016 , proteins estradaJTB ,
animals Abrams1983 ; GBMSA  or authors of scientific papers PPV2017 ;
CarlettiEtAl2020 , cannot be reduced to binary interactions. The group action
is indeed the real driver of the dynamics. Starting from this observation,
higher-order models have been developed so as to capture the many body
interactions among individual units. We hereby focus on hypergraphs
berge1973graphs ; estrada2005complex ; GZCN , versatile tools with a broad
potential that is still being fully elucidated. Hypergraphs have been applied
to different fields from social contagion model de2019social ; ATM2020 , to
the modelling of random walks CarlettiEtAl2020 , from the study of
synchronisation Krawiecki2014 ; MKJ2020 ; CarlettiJPC_2020  and diffusion
ATM2020 , to non-linear consensus neuhauser2020multibody , via the emergence
of Turing patterns CarlettiJPC_2020 . It is also worth mentioning an
alternative approach to high-order interactions which exploits the notion of
simplicial complexes DVVM ; BC ; PB . Largely used in the past to tackle
optimisation or algebraic problems, they have been recently invoked to address
problems in epidemic spreading BKS2016 ; IPBL  or synchronisation phenomena
LCB2020 ; GdPGLRCFLB2020 ; PhysRevLett.124.218301 . In this work we will
however adopt the viewpoint of hypergraphs, to represent high-order
interactions.

Hypergraphs constitute indeed a very flexible paradigm. An arbitrary number of
agents are allowed to interact: an hyperedge grouping all the involved agents
encodes for the many body interaction, thus extending conventional network
models beyond the limit of binary contacts. A hypergraph can reproduce, in a
proper limit, a simplicial complex and, in this respect, provides a more
general tool for addressing many body simultaneous interactions.

Based on the above, it can be claimed that many body interactions constitute a
relevant and transversal research field that is still in its embryonic stage,
in particular as concerns studies that relate to hypergraphs. Our contribution
is positioned in this context and aims at systematising the study of dynamical
systems coupled via a hypergraph. For a sake of definitiveness, we will hereby
consider the interactions to be mediated by the hyperedges, that is by the
(hyper)adjacency matrix (see Section 2), or by a diffusive-like process, that
is implemented via a properly engineered Laplace matrix (see Section 3). In
both cases, we will be interested in the emergence of spatially heterogeneous
solutions, i.e. coherent and extended patterns.

##  2 Hypergraphs and high-order interactions.

The aim of this section is to introduce the formalism of (hyper) adjacency
matrix which enables us to account for the high-order interaction among
several identical dynamical systems. We will then present a first study on the
emergence of spatial heterogeneous solutions, i.e. patterns, for systems
interacting via a hypergraph, by assuming that uncoupled individual units do
converge to a (spatially) homogeneous stable solution.

###  2.1 Hypergraphs

An hypergraph $\mathcal{H}(V,E)$ is defined by a set of nodes,
$V=\\{v_{1},\dots,v_{n}\\}$, and a set of $m$ hyperedges
$E=\\{E_{1},\dots,E_{m}\\}$, such that for all $\alpha=1,\dots,m$ :
${E_{\alpha}}\subset V$. If all hyperedges have size $2$ then the hypergraph
reduces to a network. A simplicial complex is recovered if each hyperedge
contains all its subsets.

One can encode the information on how the nodes are shared among hyperedges,
by using the incidence matrix of the hypergraph 111We will adopt the
convention of using roman indexes for nodes and greek ones for edges.,
$e_{i\alpha}$, namely

$e_{i\alpha}=\begin{cases}1&\text{$v_{i}\in E_{\alpha}$}\\\
0&\text{otherwise}\,.\end{cases}$ (1)

Given the latter, one can construct the $n\times n$ hypergraph adjacency
matrix,

$\mathbf{A}=\mathbf{e}\,\mathbf{e}^{\top}\,,\quad
A_{ij}=\sum_{\alpha}e_{i\alpha}e_{j\alpha}\,,$ (2)

thus $A_{ij}$ represents the number of hyperedges containing both nodes $i$
and $j$. Let us observe that often in the literature the adjacency matrix is
defined by imposing a null diagonal. In the following we will adopt a
different notation by defining its diagonal to contain all $1$’s. This in turn
amounts to assume the hypergraph to contain all the trivial hyperedges made of
just a single node. Finally we define the $m\times m$ hyperedges matrix

$\mathbf{C}=\mathbf{e}^{\top}\mathbf{e}\,,\quad
C_{\alpha\beta}=\sum_{i}e_{i\alpha}e_{i\beta}\,,$ (3)

$C_{\alpha\beta}$ counts the number of nodes in $E_{\alpha}\cap E_{\beta}$,
hence $C_{\alpha\alpha}$ is the size of the hyperedge $E_{\alpha}$.

###  2.2 High-order coupling

Let us consider a $d$-dimensional dynamical system described by the ODE :

$\frac{d\mathbf{x}}{dt}(t)=\mathbf{f}(\mathbf{x}(t))\,,$ (4)

where $\mathbf{x}(t)=(x_{1}(t),\dots,x_{d}(t))^{\top}$ denotes the state of
the system at time $t$ and $\mathbf{f}$ is a generic nonlinear function which
describes the rate of variation of $\mathbf{x}$. Assume now to replicate
system (4) into $n$ independent copies, hence yielding a (tensorial) system

$\frac{d\mathbf{x}^{(i)}}{dt}(t)=\mathbf{f}(\mathbf{x}^{(i)}(t))\quad\forall
i=1,\dots,n\,,$ (5)

where $\mathbf{x}^{(i)}(t)=(x^{(i)}_{1}(t),\dots,x^{(i)}_{d}(t))^{\top}$
denotes the state of the $i$-th copy of the generalised system. The whole
system will thus be described by the $n\times d$ vector
$\mathbf{x}=(\mathbf{x}^{(1)},\dots,\mathbf{x}^{(n)})^{\top}$. Finally we
allow each system (5) to simultaneously interact with many others, and
specifically belonging to the same hyperedge.

Let thus $E_{\alpha}$ be an hyperedge containing the $i$-th system. Then the
growth rate associated to this latter will depend on all the systems $j\neq
i$, belonging to the same hyperedge; moreover we assume such interaction to
depend also on the hyperedge size, $\varphi(C_{\alpha\alpha})$, for a generic
function $\varphi$. The system $i$ may belong to several hyperedges
$E_{\alpha}$ and thus all these contributions should be taken into account to
determine its growth rate. In formula

$\frac{d\mathbf{x}^{(i)}}{dt}(t)=\frac{\sum_{\alpha}e_{i\alpha}\sum_{j}e_{j\alpha}\varphi(C_{\alpha\alpha})\mathbf{F}(\mathbf{x}^{(i)}(t),\mathbf{x}^{(j)}(t))}{\sum_{\alpha}e_{i\alpha}\sum_{j}e_{j\alpha}\varphi(C_{\alpha\alpha})}\quad\forall
i=1,\dots,n\,,$ (6)

where we introduced the function $\mathbf{F}$ such that
$\mathbf{F}(\mathbf{x}^{(i)},\mathbf{x}^{(i)})=\mathbf{f}(\mathbf{x}^{(i)})$
and the term at the denominator acts as a normalisation factor. We will show
later on, that different functions $\mathbf{F}$ can be used to return the same
function $\mathbf{f}$.

Let us define the $m\times m$ diagonal matrix $\Phi$ such that
$\Phi_{\alpha\alpha}=\varphi(C_{\alpha\alpha})$ and zero otherwise. Then we
can rewrite Eq. (6) as follows

$\frac{d\mathbf{x}^{(i)}}{dt}(t)=\frac{1}{d_{i}}\sum_{j}D_{ij}\mathbf{F}(\mathbf{x}^{(i)}(t),\mathbf{x}^{(j)}(t))\quad\forall
i=1,\dots,n\,,$ (7)

where we introduced the matrix
$\mathbf{D}=\mathbf{e}\,\Phi\,\mathbf{e}^{\top}$ whose elements read

$D_{ij}=\sum_{\alpha}e_{i\alpha}\Phi_{\alpha\alpha}e_{j\alpha}\quad\forall
i\neq j\text{ and }D_{ii}=\varphi(1)\,.$ (8)

Let us observe that the different definition for the diagonal elements is due
to the inclusion of the trivial hyperedges containing each single node and
thus having size $1$. Finally let use define $d_{i}=\sum_{j}D_{ij}$.

######  Remark 1 (Isolated systems)

In the case $n$ systems are isolated, i.e. all the hyperedges have size $1$,
then $C_{\alpha\alpha}=1$ for all $\alpha$. Observing that a single
$\alpha^{\prime}$ (the one associated to the unique hyperedge containing $i$)
does satisfy $e_{i\alpha^{\prime}}=1$ (all the other ones being zero,
$e_{i\beta}=0$ for all $\beta=\alpha^{\prime}$), we can rewrite equation (6)
by remarking that the sum over $j$ is restricted to $j=i$:

$\frac{d\mathbf{x}^{(i)}}{dt}(t)=\frac{\varphi(1)\mathbf{F}(\mathbf{x}^{(i)}(t),\mathbf{x}^{(i)}(t))}{\varphi(1)}=\mathbf{f}(\mathbf{x}^{(i)}(t))\quad\forall
i=1,\dots,n\,,$

where use has been made of the relation
$\mathbf{F}(\mathbf{x}^{(i)},\mathbf{x}^{(i)})=\mathbf{f}(\mathbf{x}^{(i)})$.
Because our formalism contains the trivial case of isolated systems (5), it
results thus a natural extension of the latter.

######  Remark 2 (Pairwise interacting systems)

In case of systems interacting in pairs, i.e. when all hyperedges have size
$C_{\alpha\alpha}=2$ for all $\alpha$ (but the ones associated to the trivial
hyperedges containing each node), we can show that equation (6) converges back
to the usual setting of a dynamical model anchored on a conventional network
CENCETTI2020109707 , once we assume $\varphi\equiv 1$, namely the same unitary
weight is associated to each link.

First of all, let us observe that
$D_{ii}=(\mathbf{e}\,\Phi\,\mathbf{e}^{\top})_{ii}=\varphi(1){A}_{ii}$ while
for $i\neq j$ we have
$D_{ij}=(\mathbf{e}\,\Phi\,\mathbf{e}^{\top})_{ij}=\varphi(2){A}_{ij}$, where
we used the definition of the adjacency matrix that includes self-loops. Then
Eq. (7) can be rewritten as

$\frac{d\mathbf{x}^{(i)}}{dt}(t)=\frac{\sum_{j}A_{ij}\mathbf{F}(\mathbf{x}^{(i)}(t),\mathbf{x}^{(j)}(t))}{k_{i}}\quad\forall
i=1,\dots,n\,,$

where use has been made of the definition $k_{i}=\sum_{j}A_{ij}$.

###  2.3 Dynamical behaviour

Assume $\mathbf{s}(t)$ to be a solution of the initial system (4), then
$\mathbf{x}^{(i)}(t)=\mathbf{s}(t)$, $i=1,\dots,n$, is trivially also a
homogeneous solution of Eq. (5) but also of Eq. (7). Indeed, for all
$i=1,\dots,n$ one has

$\displaystyle\frac{d\mathbf{x}^{(i)}}{dt}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{d_{i}}{\sum_{j}D_{ij}\mathbf{F}(\mathbf{x}^{(i)}(t),\mathbf{x}^{(j)}(t))}\Big{\rvert}_{\mathbf{x}^{(i)}(t)=\mathbf{s}(t)}=\frac{1}{d_{i}}{\sum_{j}D_{ij}\mathbf{F}(\mathbf{s}(t),\mathbf{s}(t))}$
(9) $\displaystyle=$
$\displaystyle\frac{1}{d_{i}}{\sum_{j}D_{ij}\mathbf{f}(\mathbf{s}(t))}=\mathbf{f}(\mathbf{s}(t))\,,$

where we used the property
$\mathbf{F}(\mathbf{s},\mathbf{s})=\mathbf{f}(\mathbf{s})$ and the definition
of $d_{i}$. By definition of $\mathbf{s}$ the rightmost term equals
$\dot{\mathbf{s}}$ which thus coincides also with the leftmost term.

Consider now a spatially dependent perturbation, i.e. a node depending one,
about the homogeneous solution,
$\mathbf{x}^{(i)}(t)=\mathbf{s}(t)+\mathbf{u}^{(i)}(t)$. Insert this ansatz
into Eq. (7) and determine the evolution of $\mathbf{u}^{(i)}(t)$ by assuming
it to be small (i.e. using a first order expansion), $\forall i=1,\dots,n$:

$\displaystyle\frac{d\mathbf{u}^{(i)}}{dt}(t)+\frac{d\mathbf{s}}{dt}(t)$
$\displaystyle=$
$\displaystyle\frac{1}{d_{i}}{\sum_{j}D_{ij}\mathbf{F}(\mathbf{s}+\mathbf{u}^{(i)},\mathbf{s}+\mathbf{u}^{(j)})}$
$\displaystyle=$
$\displaystyle\mathbf{f}(\mathbf{s})+\frac{1}{d_{i}}{\sum_{j}D_{ij}\left(\sum_{\ell}\partial_{x_{\ell}^{(i)}}\mathbf{F}(\mathbf{s},\mathbf{s}){u}_{\ell}^{(i)}+\sum_{\ell}\partial_{{x}_{\ell}^{(j)}}\mathbf{F}(\mathbf{s},\mathbf{s}){u}_{\ell}^{(j)}\right)}$
$\displaystyle=$
$\displaystyle\mathbf{f}(\mathbf{s})+\sum_{\ell}\partial_{x_{\ell}^{(i)}}\mathbf{F}(\mathbf{s},\mathbf{s}){u}_{\ell}^{(i)}+\frac{1}{d_{i}}{\sum_{j}D_{ij}\sum_{\ell}\partial_{{x}_{\ell}^{(j)}}\mathbf{F}(\mathbf{s},\mathbf{s}){u}_{\ell}^{(j)}}$
$\displaystyle=$
$\displaystyle\mathbf{f}(\mathbf{s})+\mathbf{J}_{1}\mathbf{u}^{(i)}+\frac{1}{d_{i}}{\sum_{j}D_{ij}\mathbf{J}_{2}\mathbf{u}^{(j)}}\,,$

where we defined the Jacobian matrices
$\mathbf{J}_{1}=\partial_{\mathbf{x}_{1}}\mathbf{F}(\mathbf{s},\mathbf{s})$,
i.e. the derivatives are computed with respect to the first group of
variables, and
$\mathbf{J}_{2}=\partial_{\mathbf{x}_{2}}\mathbf{F}(\mathbf{s},\mathbf{s})$,
i.e. the derivatives are performed with respect to the second group of
variables. In both cases the derivatives are evaluated at the reference
solution $\mathbf{s}$.

By using the fact that $\dot{\mathbf{s}}=\mathbf{f}(\mathbf{s})$ and by
slightly rewriting the previous equation, we obtain

$\displaystyle\frac{d\mathbf{u}^{(i)}}{dt}(t)$ $\displaystyle=$
$\displaystyle\mathbf{J}_{1}\mathbf{u}^{(i)}+\frac{1}{d_{i}}{\sum_{j}D_{ij}\mathbf{J}_{2}\mathbf{u}^{(j)}}=\mathbf{J}_{1}\mathbf{u}^{(i)}+\mathbf{J}_{2}\mathbf{u}^{(i)}+\sum_{j}\left(\frac{D_{ij}}{d_{i}}-\delta_{ij}\right)\mathbf{J}_{2}\mathbf{u}^{(j)}$
$\displaystyle=$
$\displaystyle\left(\mathbf{J}_{1}+\mathbf{J}_{2}\right)\mathbf{u}^{(i)}+\sum_{j}\mathcal{L}_{ij}\mathbf{J}_{2}\mathbf{u}^{(j)}\,,$

where we defined the matrix operator

$\mathcal{L}_{ij}=\frac{D_{ij}}{d_{i}}-\delta_{ij}\,.$ (10)

By introducing the $n\times d$ vector
$\mathbf{u}=(\mathbf{u}^{(1)},\dots,\mathbf{u}^{(n)})^{\top}$ we can rewrite
the latter equation in a compact form as:

$\frac{d\mathbf{u}}{dt}(t)=\left[\left(\mathbf{J}_{1}+\mathbf{J}_{2}\right)\otimes\mathbf{I}_{n}+\mathbf{J}_{2}\otimes\mathcal{L}\right]\mathbf{u}\,,$
(11)

where $\mathbf{I}_{n}$ is the $n\times n$ identity matrix and $\otimes$ is the
Kronecker product of matrices.

One can prove that $\mathcal{L}$ is a novel (consensus) high-order Laplace
matrix 222Let us introduce
$\mathcal{L}^{\mathit{sym}}=\mathbf{d}^{-1/2}\mathbf{L}^{H}\mathbf{d}^{-1/2}$,
where $\mathbf{d}$ is the diagonal matrix containing the $d_{i}$’s on the
diagonal and $\mathbf{L}^{H}$ is the high-order (combinatorial) Laplace matrix
defined in CarlettiJPC_2020 . Then
$\mathcal{L}^{\mathit{sym}}=D_{ij}/\sqrt{d_{i}d_{j}}-\delta_{ij}$ from which
it immediately follows that $\mathcal{L}^{\mathit{sym}}$ is symmetric and
nonpositive definite; indeed take any
$\mathbf{x}\in\mathbb{R}^{N}\setminus\\{0\\}$, $N$ standing for the dimension
of the matrices, then
$(\mathbf{x},\mathcal{L}^{\mathit{sym}}\mathbf{x})=(\mathbf{d}^{-1/2}\mathbf{x},\mathbf{L}^{H}\mathbf{d}^{-1/2}\mathbf{x})\leq
0$ where the last inequality follows from the fact that $\mathbf{L}^{H}$ is
nonpositive definite. Finally let us observe that
$\mathcal{L}=\mathbf{d}^{-1}\mathbf{L}^{H}=\mathbf{d}^{-1/2}\mathcal{L}^{\mathit{sym}}\mathbf{d}^{1/2}$,
hence, $\mathcal{L}$ is similar to $\mathcal{L}^{\mathit{sym}}$ and, thus they
display the same non-positive spectrum. Moreover this implies also that
$-2\leq\Lambda^{(\alpha)}\leq 0$. , i.e. it is nonpositive definite, the
largest eigenvalue is $\Lambda^{(1)}=0$ and its is associated to the uniform
eigenvector $\phi^{(1)}\sim(1,\dots,1)^{\top}$.

Recalling the relation
$\mathbf{f}(\mathbf{x})=\mathbf{F}(\mathbf{x},\mathbf{x})$ one can prove that:

$\partial_{\mathbf{x}}\mathbf{f}:=\mathbf{J}=\mathbf{J}_{1}+\mathbf{J}_{2}\,,$

and thus rewrite Eq. (11) as

$\frac{d\mathbf{u}}{dt}(t)=\left[\mathbf{J}\otimes\mathbf{I}_{n}+\mathbf{J}_{2}\otimes\mathcal{L}\right]\mathbf{u}\,.$
(12)

This is a linear system involving matrices with size $nd\times nd$. To
progress with the analytical understanding, we employ the eigenbase of
$\mathcal{L}$, to project the former equation onto each eigendirection

$\frac{d\mathbf{u}^{(\alpha)}}{dt}(t)=\left[\mathbf{J}(\mathbf{s}(t))+\mathbf{J}_{2}(\mathbf{s}(t))\Lambda^{(\alpha)}\right]\mathbf{u}^{(\alpha)}\,,$
(13)

where $\Lambda^{(\alpha)}$ is the eigenvalue relative to the eigenvector
$\phi^{(\alpha)}$. The above equation enables us to infer the stability of the
homogeneous solution, $\mathbf{s}(t)$, by studying the Master Stability
Function, namely the real part of the largest Lyapunov exponent of Eq. (13).
To illustrate the potentiality of the theory we shall turn to considering a
specific application that we will introduce in the following.

###  2.4 Results

In the above analysis we have obtained a one-parameter family (indexed by the
eigenvalues $\Lambda^{(\alpha)}$) of linear but (in general) time dependent
systems (13). For the sake of simplicity we will hypothesise the homogenous
solution to be stationary and stable, $\mathbf{s}(t)=\mathbf{s}_{0}$. In this
way we will hence assume each isolated system to converge to the same
stationary point. This simplifies the study of Eq. (13), by allowing us to
deal with a constant linear system. Let us observe that one could in principle
study the more general setting of a time dependent solution, by using the
Floquet theory in case of a periodic orbit or the full Master Stability
Function in the case of irregular oscillators.

As a concrete application we will consider a Volterra model mckane2005predator
which describes the interaction of prey and predators in an ecological setting
:

$\begin{cases}\dot{x}=-dx+c_{1}xy\\\ \dot{y}=ry-sy^{2}-c_{2}xy\,,\end{cases}$
(14)

here $x$ denotes the concentration of predators, while $y$ stands for the prey
and $\dot{}$ the time derivative. All the parameters are assumed to be
positive; in the following we will make use of the choice $c_{1}=2$,
$c_{2}=13$, $r=1$, $s=1$ and $d=1/2$, but of course our results hold true in
general. The Volterra model (14) admits a nontrivial fixed-point,
$x^{*}=\frac{c_{1}r-sd}{c_{1}c_{2}}$, $y^{*}=\frac{d}{c_{1}}$, which is
positive and stable, provided $c_{1}r-sd>0$. In the case under scrutiny, we
have $x^{*}=3/52\sim 0.0577$ and $y^{*}=1/4$.

Following the above presented scheme, let us now considering $n$ replicas of
the model (14), each associated to a different ecological niche and indexed by
the node index $i$. Assume also that species can sense the remote interaction
with other communities populating neighbouring nodes. For instance, the
competition of prey for food and resources can be easily extended so as to
account for a larger habitat which embraces adjacent patches. At the same
time, predators can benefit from a coordinated action to hunt in team. For a
sake of definitiveness we will study in the following the high-order coupling
(let us stress once again that several “microscopic” high-order models can
give rise to the same network-aggregate model) defined by:

$\begin{cases}\dot{x}_{i}=-dx_{i}+ac_{1}y_{i}\frac{1}{d_{i}}\sum_{j}D_{ij}x_{j}+(1-a)c_{1}x_{i}\frac{1}{d_{i}}\sum_{j}D_{ij}y_{j}\\\
\dot{y}_{i}=ry_{i}-sy_{i}\frac{1}{d_{i}}\sum_{j}D_{ij}y_{j}-c_{2}y_{i}\frac{1}{d_{i}}\sum_{j}D_{ij}x_{j}\,,\end{cases}$
(15)

where the matrix $D_{ij}$ encodes for the high-order interaction among nodes
$i$ and $j$, taking into account the number and size of the hyperedges
containing both nodes (see (8)). The parameters $a\in[0,1]$ describes the
relative strength with which the predators in node $i$ increase because of the
“in-node” predation or because of the interaction among predators in the
hyperedges. The case $a=1$ corresponds to a purely in-node process while if
$a=0$ a coordinated action to hunt in team is assumed to rule the dynamics.
Preys feel the competition for the resources with preys living in nodes
belonging to the same hyperedge (second term on the right hand side of the
second equation of (15)) as well from predators in the same hyperedge
(rightmost terms in the same equation). Birth and death of both species are
local, i.e. due to resources available in-node.

By using the new Laplace matrix (10) we can rewrite the previous model (15)
as:

$\begin{cases}\dot{x}_{i}=-dx_{i}+c_{1}y_{i}x_{i}+ac_{1}y_{i}\sum_{j}\mathcal{L}_{ij}x_{j}+(1-a)c_{1}x_{i}\sum_{j}\mathcal{L}_{ij}y_{j}\\\
\dot{y}_{i}=ry_{i}-sy^{2}_{i}-c_{2}y_{i}x_{i}-sy_{i}\sum_{j}\mathcal{L}_{ij}y_{j}-c_{2}y_{i}\sum_{j}\mathcal{L}_{ij}x_{j}\,,\end{cases}$
(16)

where one can easily recognise the in-node Volterra model (14) and the
corrections stemming from high-order contributions.

As previously shown, in the general setting (see (9)) the homogenous solution
$(x^{*},y^{*})$ is also a solution of the coupled system (15), that is
$x_{i}=x^{*}$ and $y_{i}=y^{*}$ solves the latter. In the following we will
prove that such solution can be destabilised due to the high-order coupling so
driving the system towards a new heterogenous, spatially dependent, solution.
To prove this claim, we will linearise system (14) about the homogeneous
equilibrium by setting $u_{i}=x_{i}-x^{*}$ and $v_{i}=y_{i}-y^{*}$ and then
make use of the eigenbase of the Laplace matrix $\mathcal{L}$,
$(\Lambda^{(\alpha)},\phi^{(\alpha)})$, to project the linear system onto each
eigenmode, that is $u_{i}=\sum_{\alpha}u^{\alpha}\phi_{i}^{(\alpha)}$ and
$v_{i}=\sum_{\alpha}v^{\alpha}\phi_{i}^{(\alpha)}$:

$\displaystyle\frac{d}{dt}\binom{u^{\alpha}}{v^{\alpha}}$ $\displaystyle=$
$\displaystyle\left[\left(\begin{matrix}0&c_{1}x^{*}\\\
-c_{2}y^{*}&-sy^{*}\end{matrix}\right)+\Lambda^{(\alpha)}\left(\begin{matrix}ac_{1}y^{*}&(1-a)c_{1}x^{*}\\\
-c_{2}y^{*}&-sy^{*}\end{matrix}\right)\right]\binom{u^{\alpha}}{v^{\alpha}}$
(17) $\displaystyle=$
$\displaystyle\left(\mathbf{J}+\Lambda^{(\alpha)}\mathbf{J}_{2}\right)\binom{u^{\alpha}}{v^{\alpha}}=:\mathbf{J}^{(\alpha)}\binom{u^{\alpha}}{v^{\alpha}}\,.$

The homogenous solution will prove unstable if (at least) one eigenmode
$\bar{\alpha}$ exists for which the largest real part of the eigenvalues of
$\mathbf{J}^{(\bar{\alpha})}$ is positive. The real part of the largest
eigenvalue $\lambda$ as function of $\Lambda^{(\alpha)}$ is called the
dispersion relation. One can easily realise that $\lambda$ is the solution
with the largest real part of the second order equation

$\lambda^{2}-\mathrm{tr}\mathbf{J}^{(\alpha)}\lambda+\det\mathbf{J}^{(\alpha)}=0\,.$

Hence the required condition for the instability is

$\mathrm{tr}\mathbf{J}^{(\alpha)}>0\text{ or
}\mathrm{tr}\mathbf{J}^{(\alpha)}<0\text{ and }\det\mathbf{J}^{(\alpha)}<0\,.$
(18)

A straightforward computation returns

$\mathrm{tr}\mathbf{J}^{(\alpha)}=-sy^{*}+\Lambda^{(\alpha)}\left(-s+ac_{1}\right)\text{
and
}\det\mathbf{J}^{(\alpha)}=c_{1}y^{*}\left(1+\Lambda^{(\alpha)}\right)\left[\Lambda^{(\alpha)}\left(c_{2}x^{*}(1-a)-asy^{*}\right)+c_{2}x^{*}\right]\,.$

Let us recall that the homogenous equilibrium is stable for the decoupled
system corresponding to setting $\Lambda^{(1)}=0$. Indeed
$\mathrm{tr}\mathbf{J}^{(1)}=-sy^{*}<0$ and
$\det\mathbf{J}^{(1)}=c_{1}c_{2}x^{*}y^{*}>0$. We have thus to determine the
existence of (at least one) $\bar{\alpha}\geq 2$ for which the conditions for
instability (18), allowing us to prove the positivity of
$\lambda\left(\Lambda^{(\bar{\alpha})}\right)$. In Fig. 1 we report a case
where the high-order coupling is able to destabilise the homogenous solution
(panel b), thus returning a patchy solution (panels c and d) for the involved
species. Finally let us observe that interestingly some niches ($6$ over $20$)
become empty, that is deprived of any species.

t]

Figure 1: Patterns in the Volterra model with high-order interactions (I). In
panel a) we represent the hypergraph used to model the high-order interactions
among species living in different niches. The hypergraph is composed of $n=20$
nodes and it has been generated using a random attachment process and it is
composed by $20$ trivial hyperedges of size $1$, $11$ hyperedges of size $2$,
$10$ hyperedges of size $3$ and $1$ hyperedge of size $4$. In panel b) we
report the dispersion relation for the Volterra model (15), the red symbols
refer to $\lambda\left(\Lambda^{(\alpha)}\right)$, $\alpha\in\\{1,\dots,n\\}$,
while the blue line denotes the dispersion relation for the Volterra model
reformulated on a continuous support. In panel c) we show the time evolution
of the predator density in each node as a function of time, $x_{i}(t)$; let us
observe that in (almost) each node the density of predators is much larger
than the corresponding homogenous equilibrium $x^{*}\sim 0.0577$ (blue). Panel
d) report the time evolution of the prey density in each node as a function of
time, $y_{i}(t)$; let us observe that in (almost) each node the density of
preys is much lower than the corresponding homogenous equilibrium $y^{*}=1$
(green). The model parameters have been set to $c_{1}=2$, $c_{2}=13$, $r=1$,
$s=1$, $d=1/2$ and $a=1/2$. We fix $\varphi(c)=c^{\sigma}$ with $\sigma=1.5$.

Another even more interesting case is reported in Fig. 2. In this case the
uncoupled homogeneous equilibrium yields $\tilde{x}=0$ and $\tilde{y}=r/s$.
When extending the study to account for multi body interactions, predators do
survive in each niche while the preys go through extinction in a few location
($6$ nodes over $20$). Generally the density of preys is lower than the
equilibrium value found in the isolated case.

t]

Figure 2: Patterns in the Volterra model with high-order interactions (II).
Using the same hypergraph shown in Fig. 1 we study the emergence of patterns
close to the homogeneous equilibrium $\tilde{x}=0$ and $\tilde{y}=r/s=1$. We
report in panel a) the dispersion relation for the Volterra model (15), the
red symbols refer to $\lambda\left(\Lambda^{(\alpha)}\right)$,
$\alpha\in\\{1,\dots,n\\}$, while the blue line denotes the dispersion
relation for the Volterra model computed on a continuous support. In panel b)
we show the time evolution of the predator density in each node as a function
of time, $x_{i}(t)$; let us observe that in each node the density of predators
is positive in striking contrast with it happens for the uncoupled system.
Panel c) reports the time evolution of the prey density in each node as a
function of time, $y_{i}(t)$; let us observe that in each node the density of
preys is much lower than the homogenous equilibrium $y^{*}=1$ (green) and in
$8$ niches the preys have gone through extinction. The model parameters have
been set to $c_{1}=2$, $c_{2}=13$, $r=1$, $s=1$, $d=1/2$ and $a=1/2$. We fix
$\varphi(c)=c^{\sigma}$ with $\sigma=1.5$.

##  3 Hypergraph and high-order diffusive-like coupling

In the previous section we have introduced and studied the problem of the
emergence of a spatially heterogenous solution in a system of several
identical dynamical units coupled together via the (hyper) adjacency matrix of
the hypergraph. In particular the microscopic units defining the system are
constrained to stay anchored to the node where they interact with those
sharing the same location and those belonging to nodes of the incident
hyperedges. In this section we will present a modified framework based on the
assumption that the basic units can travel across the hypergraph jumping from
node to node via the available hyperedges.

Starting from the definition of hyper adjacency matrix, Eq. (2), the notion of
(combinatorial) Laplace matrix for networks can be straightforwardly
generalised to the case of hypergraphs JM2019 ; MKJ2020 , by defining
$k_{i}\delta_{ij}-A_{ij}$, where $k_{i}=\sum_{j}A_{ij}$. Let us however
observe that the latter does not account in full for the higher-order
structures encoded in the hypergraph. Notably, the sizes of the incident
hyperedges are neglected.

To overcome this limitation, authors of CarlettiEtAl2020  studied a random
walk process defined on a generic hypergraph using a new (random walk) Laplace
matrix. It is worth mentioning that the transition rates of the associated
process, linearly correlates with the size of the involved hyperedges. Stated
differently, exchanges are favoured among nodes belonging to the same
hyperedge (weighted according to its associated size). Note that a similar
construction has been proposed in Evans_2010  to extract a $n$-clique graph
from a network. The main difference in the present case is that hyperedges can
have an heterogeneous size distribution and thus provide a more flexible
framework for tackling a wide range of problems.

For the sake of completeness, let us briefly recall the construction of the
random walk process on a hypergraph and invite the interested reader to
consult CarlettiEtAl2020  for further details. The agents are located on the
nodes and hop between them. In a general setting, the walkers may weight
hyperedges depending on their size, introducing a bias in their moves that we
shall encode into a function $\varphi$ of the hyperedge size. This yields the
weighted adjacency matrix $\mathbf{D}=\mathbf{e}\,\Phi\,\mathbf{e}^{\top}$,
already defined in Eq. (8) and hereby recalled:

$D_{ij}=\sum_{\alpha}e_{i\alpha}\Phi_{\alpha\alpha}e_{j\alpha}\quad\forall
i\neq j\text{ and }D_{ii}=\varphi(1)\,,$

where $\Phi$ is the diagonal matrix whose elements read
$\varphi(C_{\alpha\alpha})$. The transition probabilities of the examined
process are then obtained by normalising the columns of the weighted adjacency
matrix $T_{ij}=\frac{D_{ij}}{d_{i}}$ for all $i$, where again
$d_{i}=\sum_{j}D_{ij}$.

Let us briefly observe that assuming $\varphi(c)=c^{\sigma}$ allows to cover
several existing models of random walks on hypergraphs. For $\sigma=1$, we get
the random walk defined in CarlettiEtAl2020 , while for $\sigma=-1$ we obtain
the one introduced by Zhou zhou2007learning . Finally, the case $\sigma=0$
returns a random walk on the so called clique reduced multigraph. The latter
is a multigraph where each pair of nodes is connected by a number of edges
equal to the number of hyperedges containing that pair in the hypergraph.

From the above introduced transition probabilities one can define the random
walk Laplacian generalising that of standard networks,
$L_{ij}=\delta_{ij}-T_{ij}$, and eventually derive the (combinatorial) Laplace
matrix,

$\mathbf{L}^{H}=\mathbf{D}-\mathbf{d}\,,$ (19)

this latter will be employed in the following to model diffusion on higher-
order structures. In the above equation, matrix $\mathbf{d}$ displays, on the
diagonal, the values $d_{i}=\sum_{j}D_{ij}$ and zeros otherwise. It is clear
from its very definition that $\mathbf{D}$ takes into account both the number
and the size of the hyperedges incident with the nodes. It can also be noted
that $\mathbf{D}$ can be considered as a weighted adjacency matrix whose
weights have been self-consistently defined so as to account for the higher-
order structures encoded in the hypergraph.

Consider again the $d$-dimensional system Eq. (4) described by local, i.e.
aspatial, equations:

$\frac{d\mathbf{x}}{dt}=\mathbf{f}(\mathbf{x})\quad\mathbf{x}\in\mathbb{R}^{d}\,,$
(20)

and assume further $n$ identical copies of the above system coupled through a
hypergraph. In this way each copy of the system attached to a node of a
hypergraph belonging to one (or more) hyperedge. Units sharing the same
hyperedge are tightly coupled, due to existing many body interactions. In
formulas:

$\frac{d\mathbf{x}_{i}}{dt}=\mathbf{f}(\mathbf{x}_{i})+\varepsilon\sum_{\alpha:i\in
E_{\alpha}}\sum_{j\in
E_{\alpha}}\varphi(C_{\alpha\,\alpha})\left(\mathbf{G}(\mathbf{x}_{j})-\mathbf{G}(\mathbf{x}_{i})\right)\,,$

where $\mathbf{x}_{i}$ denotes the state of the $i$-th unit, i.e. anchored to
the $i$-th node, $\varepsilon$ the strength of the coupling, $\varphi$ is the
function encoding the bias due to the hyperedge size and $\mathbf{G}$ a
generic nonlinear coupling function. From the definition of $e_{i\alpha}$ one
can rewrite the previous formula as

$\displaystyle\frac{d\mathbf{x}_{i}}{dt}$ $\displaystyle=$
$\displaystyle\mathbf{f}(\mathbf{x}_{i})+\varepsilon\sum_{\alpha,j}e_{i\alpha}e_{j\alpha}\varphi(C_{\alpha\,\alpha})\left(\mathbf{G}(\mathbf{x}_{j})-\mathbf{G}(\mathbf{x}_{i})\right)$
(21) $\displaystyle=$
$\displaystyle\mathbf{f}(\mathbf{x}_{i})+\varepsilon\sum_{j}D_{ij}\left(\mathbf{G}(\mathbf{x}_{j})-\mathbf{G}(\mathbf{x}_{i})\right)=\mathbf{f}(\mathbf{x}_{i})+\varepsilon\sum_{j}\left(D_{ij}-d_{i}\delta_{ij}\right)\mathbf{G}(\mathbf{x}_{j})$
$\displaystyle=$
$\displaystyle\mathbf{f}(\mathbf{x}_{i})+\varepsilon\sum_{j}L^{H}_{ij}\mathbf{G}(\mathbf{x}_{j})\,,$

where we have used the above definitions for $d_{i}$ and $L^{H}_{ij}$. Let us
stress once again that the whole high-order structure is encoded in a $n\times
n$ matrix. Hence there is no need for tensors and this simplifies the
resulting analysis.

By exploiting the fact that $\sum_{j}L^{H}_{ij}=0$ for all $i=1,\dots,n$, it
is immediate to conclude that the aspatial reference solution $\mathbf{s}(t)$,
i.e. the time dependent function solving Eq. (20), is also a solution of Eq.
(21). A natural question hence arises: what can we say of the stability of the
homogeneous solution for the system in its diffusive-like coupled variant?

To answer to this question one introduces again the deviations from the
reference orbit, i.e. $\mathbf{u}_{i}=\mathbf{x}_{i}-\mathbf{s}$. Assuming
this latter to be small, one can derive a self-consistent set of linear
differential equations for tracking the evolution of the perturbation in time.
To this end, we make use of the expression in the above Eq. (21) and perform a
Taylor expansion to the linear order of approximation, to eventually get:

$\frac{d\mathbf{u}_{i}}{dt}=\mathbf{J}(\mathbf{s}(t))\mathbf{u}_{i}+\varepsilon\sum_{j}{L}^{H}_{ij}\mathbf{J}_{\mathbf{G}}(\mathbf{s}(t))\mathbf{u}_{j}\,,$
(22)

where $\mathbf{J}(\mathbf{s}(t))$ (resp.
$\mathbf{J}_{\mathbf{G}}(\mathbf{s}(t))$) denotes the Jacobian matrix of the
function $\mathbf{f}$ (resp. $\mathbf{G}$) evaluated on the trajectory
$\mathbf{s}(t)$.

We can improve on our analytical understanding of the problem by employing
again the eigenbase of the Laplace matrix $\mathbf{L}^{H}$. Being the latter
symmetric there exists a basis of orthonormal eigenvectors,
$\phi_{H}^{(\alpha)}$, associated to the eigenvalues $\Lambda_{H}^{(\alpha)}$.
We can then project $\mathbf{u}_{i}$ on this basis and obtain, for all
$\alpha$:

$\frac{d\mathbf{y}_{\alpha}}{dt}=\left[\mathbf{J}(\mathbf{s}(t))+\varepsilon{\Lambda^{(\alpha)}_{H}}\mathbf{J}_{\mathbf{G}}(\mathbf{s}(t))\right]\mathbf{y}_{\alpha}\,,$
(23)

where $\mathbf{y}_{\alpha}$ is the projection of $\mathbf{u}_{i}$ on the
$\alpha$-th eigendirection.

The (in)stability of the homogenous solution $\mathbf{s}(t)$ can be checked by
looking at the eigenvalue of the linear system (23), and more specifically the
eigenvalue with the largest real part. In a general framework, where i.e.
$\mathbf{s}(t)$ depends on time, we are dealing with a time dependent
eigenvalue problem that can be tackled by using the Master Stability Function
Pecora ; HCLP . For simplicity we will hereby solely consider the case of a
stationary reference orbit, i.e. $\mathbf{s}(t)=\mathbf{s}_{0}$. In this way
Eq. (23) can be directly solved by using spectral methods. We invite the
interested reader to refer to CarlettiJPC_2020  where the general case of a
periodic or even a chaotic $\mathbf{s}(t)$ has been analysed.

###  3.1 Turing patterns on hypergraphs

The problem introduced in the previous section opens up the perspective to
address the notion of a Turing instability on hypergraphs. Indeed, according
to the Turing instability mechanism, a stable homogeneous equilibrium becomes
unstable upon injection of a heterogeneous, i.e. spatially dependent,
perturbation once diffusion and reaction terms are simultaneously at play. The
Turing phenomenon is exemplified with reference to $2$ dimensional systems. In
the following we will consequently assume $d=2$ and rewrite
$\mathbf{x}_{i}=(u_{i},v_{i})$ as well as
$\mathbf{f}(\mathbf{x}_{i})=\left(f(u_{i},v_{i}),g(u_{i},v_{i})\right)$, where
the index $i=1,\dots,n$ refers to the specific node to which the dynamical
variables are bound. Hence Eq. (21) becomes

$\begin{cases}\dot{u_{i}}&=f(u_{i},v_{i})+D_{u}\sum_{j}L^{H}_{ij}u_{j}\\\
\dot{v_{i}}&=g(u_{i},v_{i})+D_{v}\sum_{j}L^{H}_{ij}v_{j}\end{cases}\,,$ (24)

where $D_{u}$ and $D_{v}$ replace the diffusion coefficients of species $u$
and $v$ in the case of network and can thus be called generalised diffusion
coefficients. At first sight, the above model seems to solely account for
binary interactions. However, higher-order interactions are also present, as
encoded in the matrix $\mathbf{L}^{H}$. Finally, let us observe that if the
hypergraph is a network, then $\mathbf{L}^{H}$ reduces to the standard Laplace
matrix and thus Eqs. (24) converges to the usual reaction-diffusion system
defined on a network.

The condition for the emergence of a Turing instability can be assessed by
performing a linear stability analysis about the homogeneous equilibrium
NM2010 ; Asllani2013 ; Asllani2014NC ; Asllani2014 , as previously shown.
Assuming $\mathbf{G}$ to be the identity function and the reference orbit to
coincide with a stable stationary equilibrium $\mathbf{s}_{0}=(u_{0},v_{0})$,
Eq. (22) simplifies into:

$\begin{cases}\dot{\delta u_{i}}&=\partial_{u}f(u_{0},v_{0})\delta
u_{i}+\partial_{v}f(u_{0},v_{0})\delta v_{i}+D_{u}\sum_{j}L^{H}_{ij}\delta
u_{j}\\\ \dot{\delta v_{i}}&=\partial_{u}g(u_{0},v_{0})\delta
u_{i}+\partial_{v}g(u_{0},v_{0})\delta v_{i}+D_{v}\sum_{j}L^{H}_{ij}\delta
v_{j}\,,\end{cases}$

where $\delta u_{i}=u_{i}-u_{0}$ and $\delta v_{i}=v_{i}-v_{0}$. By exploiting
again the eigenbasis of the Laplace matrix we can write $\delta
u_{i}(t)=\sum_{\alpha}\hat{u}^{\alpha}(t)\phi^{\alpha}_{i}$ and $\delta
v_{i}(t)=\sum_{\alpha}\hat{v}^{\alpha}(t)\phi^{\alpha}_{i}$. Finally the
ansatz, $\hat{u}^{\alpha}(t)\sim e^{\lambda_{\alpha}t}$ and
$\hat{v}^{\alpha}(t)\sim e^{\lambda_{\alpha}t}$, allows us to compute the
dispersion relation, i.e. the linear growth rate
$\lambda_{\alpha}=\lambda(\Lambda_{H}^{\alpha})$ of the eigenmode $\alpha$, as
a function of the Laplacian eigenvalue $\Lambda_{H}^{\alpha}$.

As it can be straightforwardly proved, the linear growth rate is the largest
real part of the roots of the second order equation

$\lambda_{\alpha}^{2}-\lambda_{\alpha}\left[\mathrm{tr}\mathbf{J}_{0}+\Lambda_{H}^{\alpha}(D_{u}+D_{v})\right]+\det\mathbf{J}_{0}+\Lambda_{H}^{\alpha}(D_{u}\partial_{v}g+D_{v}\partial_{u}f)+D_{u}D_{v}(\Lambda_{H}^{\alpha})^{2}=0\,,$
(25)

where $\mathbf{J}_{0}=\left(\begin{smallmatrix}\partial_{u}f&\partial_{v}f\\\
\partial_{u}g&\partial_{v}g\end{smallmatrix}\right)$ is the Jacobian matrix of
the reaction part evaluated at the equilibrium $(u_{i},v_{i})=(u_{0},v_{0})$.
In Eq. (25), $\mathrm{tr}(\cdot)$ and $\det(\cdot)$ stand respectively for the
trace and the determinant. The existence of at least one eigenvalue
$\Lambda_{H}^{\tilde{\alpha}}$ for which the dispersion relation takes
positive values, implies that the system goes unstable via a typical path
first identified by Alan Turing in his seminal work. At variance, if the
dispersion relation is negative the system cannot undergo a Turing
instability: any tiny perturbation fades away and the system settles back to
the homogeneous equilibrium.

To proceed further with a concrete example we selected the Brusselator
reaction system PrigogineNicolis1967 ; PrigogineLefever1968 . This is a
nonlinear model defined by $f(u,v)=1-(b+1)u+cu^{2}v$ and $g(u,v)=bu-cu^{2}v$,
where $b$ and $c$ act as tunable parameters. In Fig. 3 we report the results
for a choice of the model parameters giving rise to Turing patterns ($b=4$,
$c=6$, $D_{u}=0.02$ and $D_{v}=0.17$) and the same hypergraph previously used
in Figs. 1 and 2. The dispersion relation (panel a) is clearly positive for a
selection of $\Lambda_{H}^{(\alpha)}$ (red points). The homogeneous solution
becomes hence unstable and the ensuing patterns are displayed in panels b) and
c).

t]

Figure 3: Turing patterns in the Brusselator model with high-order diffusive-
like couplling. Using the same hypergraph shown in Fig. 1 we study the Turing
patterns emerging from the homogeneous equilibrium $({u}_{0},{v}_{0})$. We
report in panel a) the dispersion relation for the Brusselator model defined
by the reaction terms $f(u,v)=1-(b+1)u+cu^{2}v$ and $g(u,v)=bu-cu^{2}v$; the
red symbols refer to $\lambda\left(\Lambda_{H}^{(\alpha)}\right)$,
$\alpha\in\\{1,\dots,n\\}$, while the blue line denotes the dispersion
relation for the Brusselator model defined on a continuous support. In panel
b) we show the time evolution of the $u$ variable in each node as a function
of time, $u_{i}(t)$. Panel c) reports the time evolution of the $v$ variable
in each node as a function of time, $v_{i}(t)$. The model parameters have been
set to $b=4$, $c=6$, $D_{u}=0.02$ and $D_{v}=0.17$. Hence $u_{0}=1$ and
$v_{0}=b/c=2/3$. We fix $\varphi(c)=c^{\sigma}$ with $\sigma=1.5$.

##  4 Conclusions

Complex systems are composed of a large number of simple units, mutually
interacting via nonlinear exchanges. Many-body interactions sit hence at the
root of a large plethora of spontaneously emerging phenomena, as exhibited by
complex systems. The former are often reduced to a vast collection of pairwise
interactions, involving agents interacting in pairs. This enables one to model
the inspected problem as a dynamical system flowing on a conventional binary
network, a powerful approximation that allows for progresses to be made. In
many cases of interest, this reductionist choice constitutes a rough first
order approximation to the examined dynamics and more precise models are to be
invoked which encompass for the high-order interactions being at play.

In this work, we presented a general framework which allows one to account for
multi-body interacting systems coupled via a hypergraph. This materialises in
a natural extension of the conventional network paradigm. More specifically,
we considered the problem of the emergence of heterogeneous stable solutions
in interconnected systems, under the assumption that, once isolated, all units
converge to the same, and thus globally homogenous, solution. The high-order
interaction is the driver of the resulting patchy states, which emerge as
follow a symmetry breaking instability caused by the injection of a tiny non
homogeneous perturbation. This can be though as a generalisation of the Turing
instability on hypergraphs. In particular, we considered the interaction
mediated by the number of interacting neighbouring units, namely the size of
the hyperedge, and a diffusive-like process, again biased by the number of
neighbours. In both cases we provided sufficient conditions for the emergence
of spatial patterns.

Our findings have been corroborated by numerical simulations applied to two
reference models. A Volterra model that describes the interaction among
predators and prey in ecological niches, and the Brusselator model, a
prototype model of nonlinear dynamics, that describes the interaction among
reacting and diffusing chemicals.

The proposed framework goes beyond the examples hereby presented and, because
of its generality, it could prove useful in tackling those problems were
simultaneous many-body interactions within a complex environment are to be
properly accounted for.

## References

  * [1] P. A. Abrams.  Arguments in favor of higher order interactions.  The American Naturalist, 121:887, 1983.
  * [2] Réka Albert and Albert-László Barabási.  Statistical mechanics of complex networks.  Reviews of modern physics, 74(1):47, 2002.
  * [3] M. Asllani, D.M. Busiello, T. Carletti, D. Fanelli, and G. Planchon.  Turing patterns in multiplex networks.  Physical Review E, 90:042814, 2014.
  * [4] Malbor Asllani, Tommaso Biancalani, Duccio Fanelli, and Alan J. McKane.  The linear noise approximation for reaction-diffusion systems on networks.  The European Physical Journal B, 86(11):476, 2013.
  * [5] Malbor Asllani, Joseph D. Challenger, Francesco Saverio Pavone, Leonardo Sacconi, and Duccio Fanelli.  The theory of pattern formation on directed networks.  Nature Communications, 5(1):4517, 2014.
  * [6] Albert-László Barabási et al.  Network science.  Cambridge university press, 2016.
  * [7] Austin R Benson, David F Gleich, and Jure Leskovec.  Higher-order organization of complex networks.  Science, 353(6295):163–166, 2016.
  * [8] Claude Berge.  Graphs and hypergraphs.  North-Holland Pub. Co. American Elsevier Pub. Co, 1973.
  * [9] Stefano Boccaletti, Vito Latora, Yamir Moreno, Martin Chavez, and D-U Hwang.  Complex networks: Structure and dynamics.  Physics Reports, 424(4-5):175–308, 2006.
  * [10] Á. Bodó, G.Y. Katona, and P.L. Simon.  Sis epidemic propagation on hypergraphs.  Bull. Math. Biol., 78(4):713, 2016.
  * [11] T Carletti, F. Battiston, G. Cencetti, and D. Fanelli.  Random walks on hypergraphs.  Physical Review E, 101:022308, 2020.
  * [12] Timoteo Carletti, Duccio Fanelli, and Sara Nicoletti.  Dynamical systems on hypergraphs.  Journal of Physics: Complexity, 1(3):035006, aug 2020.
  * [13] Giulia Cencetti, Federico Battiston, Timoteo Carletti, and Duccio Fanelli.  Generalized patterns from local and non local reactions.  Chaos, Solitons and Fractals, 134:109707, 2020.
  * [14] Owen T Courtney and Ginestra Bianconi.  Generalized network structures: The configuration model and the canonical ensemble of simplicial complexes.  Physical Review E, 93(6):062311, 2016.
  * [15] Guilherme Ferraz de Arruda, Giovanni Petri, and Yamir Moreno.  Social contagion models on hypergraphs.  Phys. Rev. Research, 2:023032, 2020.
  * [16] Karel Devriendt and Piet Van Mieghem.  The simplex geometry of graphs.  Journal of Complex Networks, 7(4):469–490, 2019.
  * [17] Ernesto Estrada and Juan A Rodríguez-Velázquez.  Complex networks as hypergraphs.  arXiv preprint physics/0505137, 2005.
  * [18] Ernesto Estrada and G.J. Ross.  Centralities in simplicial complexes. applications to protein interaction networks.  J. Their. Biol., 438:46, 2018.
  * [19] T S Evans.  Clique graphs and overlapping communities.  Journal of Statistical Mechanics: Theory and Experiment, 2010(12):P12037, dec 2010.
  * [20] Guilherme Ferraz de Arruda, Michele Tizzani, and Yamir Moreno.  Phase transitions and stability of dynamical processes on hypergraphs.  arXiv preprint arXiv:2005.10891, 2020.
  * [21] Santo Fortunato and Darko Hric.  Community detection in networks: A user guide.  Physics Reports, 2016.
  * [22] L. V. Gambuzza, F. Di Patti, Gallo L., S. Lepri, M. Romance, R. Criado, M. Frasca, V. Latora, and S. Boccaletti.  The master stability function for synchronization in simplicial complexes.  arXiv preprint arXiv:2004.03913v1, 2020.
  * [23] Gourab Ghoshal, Vinko Zlatić, Guido Caldarelli, and Mark EJ Newman.  Random hypergraphs and their applications.  Physical Review E, 79(6):066118, 2009.
  * [24] Jacopo Grilli, György Barabás, Matthew J Michalska-Smith, and Stefano Allesina.  Higher-order interactions stabilize dynamics in competitive network models.  Nature, 548(7666):210, 2017.
  * [25] Liang Huang, Qingfei Chen, Ying-Cheng Lai, and Louis M Pecora.  Generic behavior of master-stability functions in coupled nonlinear dynamical systems.  Physical Review E, 80:036204, 2009.
  * [26] Iacopo Iacopini, Giovanni Petri, Alain Barrat, and Vito Latora.  Simplicial models of social contagion.  Nature communications, 10(1):2485, 2019.
  * [27] J. Jost and R. Mulas.  Hypergraph laplace operators for chemical reaction networks.  Advances in Mathematics, 351:870, 2019.
  * [28] A. Krawiecki.  Chaotic synchronization on complex hypergraphs.  Chaos, Solitons and Fractals, 65:44, 2014.
  * [29] R. Lambiotte, M. Rosvall, and I. Scholtes.  From networks to optimal higher-order models of complex systems.  Nat. Phys., 15:313, 2019.
  * [30] Vito Latora, Vincenzo Nicosia, and Giovanni Russo.  Complex networks: principles, methods and applications.  Cambridge University Press, 2017.
  * [31] L.-D Lord, P. Expert, H.M. Fernandes, G. Petri, T.J. Van Hartevelt, F. Vaccarino, G. Deco, F. Turkheimer, and M.L. Kringelbach.  Insights into brain architectures from the homological scaffolds of functional connectivity networks.  Front. Syst. Neurosci., 10:85, 2016.
  * [32] Maxime Lucas, Giulia Cencetti, and Federico Battiston.  A multi-order laplacian framework for the stability of higher-order synchronization.  arXiv preprint arXiv: 2003.09734v1, 2020.
  * [33] Alan J McKane and Timothy J Newman.  Predator-prey cycles from resonant amplification of demographic stochasticity.  Physical review letters, 94(21):218102, 2005.
  * [34] Ana P. Millán, Joaquín J. Torres, and Ginestra Bianconi.  Explosive higher-order kuramoto dynamics on simplicial complexes.  Phys. Rev. Lett., 124:218301, 2020.
  * [35] Raffaella Mulas, Christian Kuehn, and Jürgen Jost.  Coupled dynamics on hypergraphs: Master stability of steady states and synchronization.  Physical Review E, 101:062313, 2020.
  * [36] Hiroya Nakao and Alexander S. Mikhailov.  Turing patterns in network-organized activator-inhibitor systems.  Nature Physics, 6:544, 2010.
  * [37] Leonie Neuhäuser, Andrew Mellor, and Renaud Lambiotte.  Multibody interactions and nonlinear consensus dynamics on networked systems.  Physical Review E, 101(3):032310, 2020.
  * [38] Mark EJ Newman.  Networks: An Introduction.  Oxford University Press, 2010.
  * [39] A. Patania, G. Petri, and F. Vaccarino.  The shape of collaborations.  EPJ Data Sci., 6:18, 2017.
  * [40] Louis M Pecora and Thomas L Carroll.  Master stability functions for synchronized coupled systems.  Physical Review Letters, 80(10):2109, 1998.
  * [41] Giovanni Petri and Alain Barrat.  Simplicial activity driven model.  Physical Review Letters, 121(22):228301, 2018.
  * [42] Giovanni Petri, Paul Expert, Federico Turkheimer, Robin Carhart-Harris, David Nutt, Peter J Hellyer, and Francesco Vaccarino.  Homological scaffolds of brain functional networks.  Journal of The Royal Society Interface, 11(101):20140873, 2014.
  * [43] I. Prigogine and R. Lefever.  Symmetry breaking instabilities in dissipative systems. ii.  J. Chem. Phys., 48:1695, 1968.
  * [44] I. Prigogine and G. Nicolis.  Symmetry breaking instabilities in dissipative systems.  J. Chem. Phys., 46:3542, 1967.
  * [45] Dengyong Zhou, Jiayuan Huang, and Bernhard Schölkopf.  Learning with hypergraphs: Clustering, classification, and embedding.  In Advances in neural information processing systems, pages 1601–1608, 2007.


# Backscattering of topologically protected helical edge states by line
defects

Mohadese Karimi  Mohsen Amini<EMAIL_ADDRESS>Morteza Soltani
<EMAIL_ADDRESS>Mozhgan Sadeghizadeh  Department of Physics, Faculty
of Physics, University of Isfahan, Isfahan 81746-73441, Iran

###### Abstract

The quantization of conductance in the presence of non-magnetic point defects
is a consequence of topological protection and the spin-momentum locking of
helical edge states in two-dimensional topological insulators. This protection
ensures the absence of backscattering of helical edge modes in the quantum
Hall phase of the system. However, our study focuses on exploring a novel
approach to disrupt this protection. We propose that a linear arrangement of
on-site impurities can effectively lift the topological protection of edge
states in the Kane-Mele model. To investigate this phenomenon, we consider an
armchair ribbon containing a line defect spanning its width. Utilizing the
tight-binding model and non-equilibrium Green’s function method, we calculate
the transmission coefficient of the system. Our results reveal a suppression
of conductance at energies near the lower edge of the bulk gap for positive
on-site potentials. To further comprehend this behavior, we perform analytical
calculations and discuss the formation of an impurity channel. This channel
arises due to the overlap of in-gap bound states, linking the bottom edge of
the ribbon to its top edge, consequently facilitating backscattering. Our
explanation is supported by the analysis of the local density of states at
sites near the position of impurities.

###### pacs:

##  I Introduction

The emergence of two-dimensional topological insulators (2D TIs) has
revolutionized the field of condensed matter physics by introducing a
fascinating phenomenon known as bulk-boundary correspondence topol1 . This
principle dictates the existence of metallic gap states that are confined to
the one-dimensional edges of the material. These edge states exhibit a
distinct property of spin polarization, where the spin of an electron becomes
rigidly coupled to its momentum. The intriguing consequence of this phenomenon
is that electrons traveling along the edge exhibit remarkable immunity to
backscattering caused by non-magnetic defects, leading to conductance
quantization. This phenomenon is known as topological edge state protection
and lies at the heart of the celebrated quantum spin Hall effect topol1 ;
topol2 ; rev2020 .

The study of topological protection in 2D TIs raises a crucial question
concerning the mechanisms that can lift this protection. This issue holds a
dual significance. Firstly, it pertains to the puzzling deviations from
conductance quantization observed in experiments on 2D TIs PZ1 ; PZ2 ; PZ3 ,
necessitating a deeper understanding of the underlying mechanisms that induce
the lifting of protection US1 ; US2 ; US3 . Secondly, in many practical
scenarios, the ability to control and manipulate edge states becomes paramount
Yazdani ; APP1 ; APP2 ; APP3 ; APP4 ; APP5 ; APP6 . Depending on the specific
device and application, there arises a need to modulate or suppress these edge
states, underscoring the importance of comprehending the mechanisms behind the
lifting of protection.

So far, numerous proposals have been put forth to explain the mechanisms
behind the lifting of topological protection and the subsequent occurrence of
backscattering in 2D TIs. Possible sources of backscattering include
mechanisms that violate time-reversal symmetry, such as the influence of an
external magnetic field MF , the presence of charge puddles PRL2013 , the
interaction with embedded nuclear spins SP1 ; SP2 , coupling to phonons
PRL2012 , and the impact of electromagnetic noise PRL2018 . However, a more
direct and controllable approach for lifting topological protection is
achieved through tunneling between opposite edges of a 2D TI APP4 ; APP5 ;
APP6 ; T1 ; T2 ; T3 ; T4 . The tunneling between opposite edges of a 2D TI
enables the coupling of electrons moving in one direction with their
counterparts of the same spin orientation on the opposing edge. As a result, a
small gap is opened at the Fermi level Amini1 , and more significantly, a
channel for electron backscattering is created without the need to break time-
reversal symmetry. A potential scenario for achieving this condition involves
the creation of spatially extended defects within a 2D TI, which can be
engineered through methods such as nano-patterning or the introduction of
specific line defects in the atomic lattice Nan1 .

Recently, a promising achievement in the observation of edge coupling due to
the presence of an extended linear defect in a 2D TI is reported
experimentally NP . In this experimental study, bismuthene, which serves as a
prototypical 2D TI is utilized. Interestingly, it is observed that narrow
constrictions spontaneously manifest themselves within the material in the
form of line defects. By employing scanning tunneling microscopy/spectroscopy
(STM/STS) measurements and analyzing the local density of states (LDOS), it
has been discovered that the presence of the line defect leads to a spatial
overlap of the edge states localized along both edges of the line defect. This
spatial overlap induces a phenomenon known as hybridization, which involves
the mixing of the edge states from each edgeNP . As a consequence of this
hybridization, inter-edge scattering occurs, giving rise to a gap in the
energy spectrum. The presence of this gap creates a channel through which back
scattering can take place.

In this study, we introduce a novel mechanism for the disruption of
conductance quantization in 2D TIs. Our proposed mechanism involves the
application of line defects, which result in the breakdown of conductance
quantization without the need for direct coupling of edge states. Importantly,
this breakdown occurs without the opening of a gap in the energy spectrum. In
this study, we focus on a 2DTI ribbon that contains a line of on-site
impurities (line defect) along its longitudinal direction. Our objective is to
investigate the robustness of the transport properties of the system using the
non-equilibrium Green’s functions formalismDatta . By employing this
framework, we aim to gain insights into how the presence of line defects
affects the transport characteristics of the 2DTI ribbon. Utilizing a
combination of numerical simulations and analytical calculations, we have
discovered that the strength of the on-site potential plays a crucial role in
the emergence of quasi one-dimensional eigenstates along the line defect in
the 2DTI ribbon. This intriguing phenomenon transforms the line defect into a
quantum wire that facilitates the transportation of current across the width
of the ribbon. Moreover, the presence of the line defect introduces a
mechanism for backscattering, which affects the transport properties of the
system. Our findings shed light on the intricate interplay between on-site
potential, line defects, and the transport behavior in 2DTI ribbons, offering
new insights into the manipulation and control of quantum states in these
systems.

The rest of the paper is organized as follows. In Section II, we provide a
comprehensive description of the model system under consideration and outline
the methodology employed to study the transmission through this system. Moving
to Section III, we present the results obtained from our quantum transport
study of the system in the presence of the line defect. We thoroughly discuss
the mechanism responsible for the breakdown of topological protection through
analytical calculations, elucidating the formation of the impurity channel
resulting from bound states in the bulk gap region. Finally, in Section IV, we
conclude our findings.

##  II Model and method

Figure 1:  Schematic representation of a two-terminal device with an
infinitely long Kane-Mele ribbon with armchair boundaries connected to
macroscopic leads on the left and right sides. The scattering region contains
a line defect of on-site impurities arranged in a zigzag chain across the
width of the ribbon (yellow circles surrounded by dashed orange rectangle).
The sign of $\nu_{ij}$ for the spin-orbit coupling term is shown in the
scattering region. The width of the ribbon is $W=6$ armchair chains.

###  II.1 Model

In this section, we present a setup designed to investigate the transport
properties of a 2D TI ribbon with an embedded line defect. Our device consists
of two conducting leads that enclose a scattering region, as depicted in Fig.
1. The scattering region includes a line of sites with defects, which are
visually represented by distinct colors in Fig. 1. As it is shown we consider
a two-dimensional honeycomb lattice with armchair edges along the
$x$-direction. In order to emulate the transport through an infinitely long
ribbon in the $x$-direction, the ribbon under consideration in this model has
a finite width $W$ in the $y$-direction. To achieve this, two semi-infinite
leads are connected to the central region, as illustrated in Fig. 1. The
Hamiltonian governing the scattering region as well as the conducting leads in
the tight-binding approximation is expressed as KM1 ; KM2 :

$\displaystyle\mathcal{H}_{\text{KM}}$ $\displaystyle=$ $\displaystyle
t\sum_{\langle
i,j\rangle,\alpha}c_{i\alpha}^{\dagger}c_{j\alpha}+i\lambda_{\text{SO}}\sum_{\langle\langle
i,j\rangle\rangle,\alpha\beta}v_{ij}c_{i\alpha}^{\dagger}s^{z}_{\alpha\beta}c_{j\beta}$
$\displaystyle+$
$\displaystyle\lambda_{\nu}\sum_{i\alpha}\xi_{i}c_{i\alpha}^{\dagger}c_{i\alpha}+\text{h.c.}$

This Hamiltonian consists of several terms. The first term corresponds to the
nearest-neighbor (NN) hopping with amplitude $t$, where
$c_{i\alpha}^{\dagger}(c_{i\alpha})$ denotes the creation (annihilation)
operator for an electron with spin $\alpha$ at site $i$. The summation
$\langle i,j\rangle$ runs over all the NN sites. The second term represents
the intrinsic spin-orbit coupling (SOC) with coupling strength
$\lambda_{\text{SO}}$ between next-nearest-neighbor (NNN) sites. The summation
$\langle\langle i,j\rangle\rangle$ in the index indicates that it runs over
all pairs of NNN sites and the Pauli matrices $\bm{s}=(s^{x},s^{y},s^{z})$ are
associated with the physical spins. The factor
$v_{ij}=\frac{\bm{d}_{i}\times\bm{d}_{j}}{|\bm{d}_{i}\times\bm{d}_{j}|}=\pm 1$
in the second term depends on the hopping path between NNN sites $i$ and $j$.
It is determined by the cross product of the vectors $\bm{d}_{i}$ and
$\bm{d}_{j}$ connecting the NNN sites, and takes values of $\pm 1$ as shown in
Fig. 1. The third term in the Hamiltonian corresponds to the staggered
sublattice potential with strength $\lambda_{\nu}$. It introduces a potential
difference between the two sublattices of the honeycomb lattice. The
sublattice index $\xi_{i}$ takes values of $\pm 1$ depending on the sublattice
of site $i$. The term $h\cdot c$ represents the Hermitian conjugate of the
previous terms.

To introduce a line of impurities, we incorporate a zigzag chain of sites
spanning across the width of the ribbon in the scattering region. These sites
possess an additional on-site potential $U$, creating a line defect indicated
by a red rectangle in Fig. 1. The corresponding Hamiltonian for this line
defect can be written as follows:

$\mathcal{H}_{\text{LD}}=U\sum_{i\in\text{LD},\alpha}c_{i\alpha}^{\dagger}c_{i\alpha}.$
(2)

###  II.2 Method

We employ the non-equilibrium Green’s functions formalism to investigate the
electronic transport properties of our system Datta . This approach allows us
to analyze the flow of electrons and calculate quantities such as the
transmission coefficient and conductance. The transmission coefficient
($\tau$) at zero temperature within this formalism can be computed using the
Landauer-Buttiker formula, given by:

$\tau(E)=\text{Tr}[\Gamma_{L}(E)G_{r}(E)\Gamma_{R}(E)G_{a}(E)].$ (3)

In this context, the retarded Green’s function in the site representation,
denoted as $G_{r}(E)$, is given by
$G_{r}(E)=[E-\mathcal{H}_{C}-\Sigma_{R}(E)-\Sigma_{L}(E)]^{-1}$, where
$\mathcal{H}_{C}=\mathcal{H}_{KM}+\mathcal{H}_{LD}$ represents the Hamiltonian
of the scattering region, which incorporates the influence of the line defect.
Similarly, the advanced Green’s function, denoted as $G_{a}(E)$, is defined as
the Hermitian conjugate of $G_{r}(E)$, i.e., $G_{a}(E)=[G_{r}(E)]^{\dagger}$.
The self-energy terms $\Sigma_{R(L)}$ correspond to the embedding self-energy,
which relies on the retarded contact Green’s functions and the coupling
between the leads $R(L)$ and the central (scattering) region. The right and
left line-width functions, denoted as $\Gamma_{R(L)}(E)$, can be expressed as
$\Gamma_{R(L)}(E)=i[\Sigma_{R(L)}-(\Sigma_{R(L)})^{\dagger}]$.

##  III Results and discussion

Figure 2:  (a) Electronic band structure of the pristine Kane-Mele model on a
ribbon with armchair terminations and a width of $W=30$. The spin-orbit
coupling strength is $\lambda_{\text{SO}}=0.2t$. (b) The corresponding zero-
temperature transmission coefficient $\tau(E)$ as a function of energy for
both absence $(U=0)$ and presence $(U=1.0)$ of the line defect in the
scattering region. The transmission coefficient curve demonstrates the
quantization plateaus in the absence of the line defect, indicating
topological protection, and the suppressed transmission regions caused by the
line defect, leading to dips in the transmission coefficient curve. Two
specific energy points with large and small suppression are highlighted at
$E=-0.88t$ and $E=-0.68t$, respectively.

In this section, we analyze the electronic transport properties of the
structure mentioned above. In our calculations, we take the NN hopping
amplitude $t$ as the energy unit. Additionally, we set the width of the ribbon
to be $W=30$, a value chosen to be sufficiently large to avoid coupling of
edge states Amini1  and fix $\lambda_{\text{SO}}=0.2t$. To numerically
calculate the transmission coefficients $\tau(E)$, we will utilize the PYQULA
library Pyqula .

###  III.1 Breakdown of the conductance quantization

To begin, let us examine the band structure of the pristine system ($U=0$),
which is illustrated in Fig. 2(a). As expected, the band structure displays a
bulk band gap of approximately $\Delta=6\sqrt{3}\lambda_{\text{SO}}$ KM1 ; KM2
, along with the presence of edge bands near the gap region. The corresponding
transmission coefficient $\tau(E)$ for the edge states is depicted by the blue
lines in Fig. 2(b). Notably, in the vicinity of the gap region, the
transmission coefficient of the edge states is $\tau(E)=1$, indicating their
protected nature.

However, when we introduce a line defect with $U\neq 0$, the situation becomes
significantly different. Our calculations reveal that for $U=1.0t$, the
transmission coefficient of the edge states is suppressed for energies near
the bottom of the bulk bands ($E\sim-0.9t$). This leads to the backscattering
of the corresponding edge states and consequently results in the lifting of
their protection. In the subsequent analysis, we will delve into the
underlying reasons for this intriguing phenomenon.

###  III.2 Analysis of the impurity channel formation

In this subsection, we will demonstrate that the lifting of the topological
protection of the edge states in the presence of the line defect can be
attributed to the formation of an additional channel that connects the top
edge of the ribbon to the bottom edge. To address this phenomenon
analytically, we will examine the in-gap impurity states, which play a crucial
role in the creation of this conducting channel. To achieve this, we will
begin by analyzing the in-gap impurity states of a Su-Schrieffer-Heeger (SSH)
SSH  model. Subsequently, we will extend and generalize these findings to our
Kane-Mele ribbon. By studying the SSH model first, we can gain insights into
the formation of in-gap impurity states, which will provide a foundation for
our analysis of the Kane-Mele ribbon’s behavior in the presence of a line
defect.

####  III.2.1 In-gap impurity state in the SSH model

Figure 3:  The geometry of an SSH chain with two types of sub-lattices,
labeled as $A$ and $B$. The chain’s unit cell is enclosed by a dashed
rectangle, and the intracell and intercell hopping amplitudes are represented
by $t_{1}$ and $t_{2}$, respectively.

In this section, we investigate an SSH chain with an on-site impurity, as
depicted in Fig. 3. The SSH chain consists of two nonequivalent sites, labeled
as $A$ and $B$, within each periodic cell. The Hamiltonian of the system is
given by the sum of two terms, the SSH term and the impurity term, expressed
as:

$\mathcal{H}=\mathcal{H}_{\text{SSH}}+\mathcal{H}_{\text{impurity}},$ (4)

where

$\mathcal{H}_{\text{SSH}}=\sum_{i}t_{1}a_{i}^{\dagger}b_{i}+t_{2}b_{i}^{\dagger}a_{i+1}+\text{h.c.},$
(5)

and

$\mathcal{H}_{\text{impurity}}=Ub_{0}^{\dagger}b_{0}.$ (6)

Here, $a_{i}^{\dagger}(a_{i+1})$ and $b_{i}^{\dagger}(b_{i})$ are the creation
(annihilation) operators at the corresponding sites in the unit cell, with
different coupling amplitudes denoted as $t_{1}$ and $t_{2}$ for the intra-
cell and inter-cell couplings respectively. The parameter $U$ represents the
strength of the on-site impurity potential which is considered to be located
at site $i=0$.

Let us begin by considering the case of the clean system, where $U=0$. We are
primarily interested in the bulk properties of the system which implies
considering a sufficiently long chain Book . To analyze the bulk properties of
the system in momentum space, we apply the Fourier transformation to the SSH
term $\mathcal{H}_{\text{SSH}}$ of Eq. (5), leading to the corresponding
momentum-space Hamiltonian:

$\mathcal{H}_{\text{SSH}}(k)=\sum_{k}\begin{bmatrix}a_{k}^{\dagger}&b_{k}^{\dagger}\end{bmatrix}\begin{bmatrix}h(k)\end{bmatrix}\begin{bmatrix}a_{k}\\\
b_{k}\end{bmatrix},$ (7)

where

$\displaystyle h(k)$ $\displaystyle=$
$\displaystyle\begin{bmatrix}0&t_{1}+t_{2}e^{ik}\\\
t_{1}+t_{2}e^{-ik}&0\end{bmatrix}$ (8) $\displaystyle=$
$\displaystyle(t_{1}+t_{2}\cos(k))\sigma_{x}+(t_{2}\sin(k))\sigma_{y}.$

In the above equation, $\sigma_{x}$ and $\sigma_{y}$ represent the Pauli
matrices. Now, the energy spectrum of $h(k)$ can be obtained by calculating
its eigenvalues, which are given by:

$E^{\pm}(k)=\pm\sqrt{t_{1}^{2}+t_{2}^{2}+2t_{1}t_{2}\cos k}.$ (9)

As a result, the energy gap between the two energy bands is given by:

$E_{g}=2|t_{2}-t_{1}|.$ (10)

In this model, when the intercell hopping parameter $(t_{1})$ is smaller than
the intracell hopping parameter $(t_{2})$, the chain exhibits a dimerized
pattern, leading to a topologically nontrivial phase characterized by a
winding number $(W)$ of 1 Book . Conversely, when $t_{1}$ is larger than
$t_{2}$, the chain becomes uniform, resulting in a trivial phase with $W=0$.
We will now focus on the topologically nontrivial phase where $t_{2}>t_{1}$
and attempt to find the end modes, also known as edge states, which satisfy
the Schrödinger equation
$\mathcal{H}_{\text{SSH}}|\psi_{\text{edge}}\rangle=0$. Due to the chiral
symmetry of the system, the corresponding edge state can be expressed in such
a way that it has support on only one sub-lattice Book , for example, sub-
lattice $A$. Therefore, we have $\langle 0|b_{i}|\psi_{\text{edge}}\rangle=0$
and the solution of the zero-energy eigenstate can be expressed as:

$|\psi_{\text{edge}}\rangle=\sum_{i}\psi^{A}_{i}a_{i}^{\dagger}|0\rangle=\sum_{i}\left(-\frac{t_{1}}{t_{2}}\right)^{i}\psi^{A}_{1}a_{i}^{\dagger}|0\rangle,$
(11)

where the coefficient $\psi^{A}_{1}$ can be determined using the normalization
condition as $(\psi^{A}_{1})^{2}=1-\left(\frac{t_{1}}{t_{2}}\right)^{2}$.

Figure 4:  Graphical representation of the solutions to the equation set (16).
In (a), we observe the variation of the bound state energy $E_{s}$ (measured
in units of $t$) as a function of the impurity potential $U/t$ for a single
impurity located on the SSH chain with $t_{2}/t1=0.5$. Notably, for extremely
large values of the on-site potential, the bound state energy tends to zero,
which aligns with the energy of the edge states in the system. (b) illustrates
the corresponding coefficient $\alpha$, as defined in the equation set
(LABEL:E13).

Now, let us consider the case where $U\neq 0$. We aim to demonstrate that the
presence of an on-site impurity potential induces a bound state with energy
$E_{s}$ in the gap region, $-E_{g}/2<E_{s}<E_{g}/2$. To achieve this, we
assume that the impurity is located on a $B$ site as depicted in Fig. 3. We
begin by making an initial assumption about the form of the wave function
associated with such an in-gap bound state as

$|\psi_{\text{s}}\rangle=|\psi_{A}^{r}\rangle+|\psi_{A}^{l}\rangle+|\psi_{B}\rangle.$
(12)

Here $|\psi_{A}^{r}\rangle$ and $|\psi_{A}^{l}\rangle$ represent the
contributions involving solely the A sub-lattice sites on the right and left
sides, respectively, while $|\psi_{B}\rangle$ denotes the contribution
involving solely the B sub-lattice sites and are defined as:

$\displaystyle|\psi_{A}^{r}\rangle$ $\displaystyle=$
$\displaystyle\sum_{i>0}\psi_{i}^{A_{r}}a^{\dagger}_{i}|0\rangle=\sum_{i>0}\psi^{A_{r}}_{1}(-\alpha)^{i}a^{\dagger}_{i}|0\rangle,$
$\displaystyle|\psi_{A}^{l}\rangle$ $\displaystyle=$ $\displaystyle\sum_{i\leq
0}\psi_{i}^{A_{l}}a^{\dagger}_{i}|0\rangle=\sum_{i\leq
0}\psi^{A_{l}}_{0}(-\alpha)^{i}a^{\dagger}_{i}|0\rangle,$
$\displaystyle|\psi_{B}\rangle$ $\displaystyle=$
$\displaystyle\sum_{i}\psi_{i}^{B}b^{\dagger}_{i}|0\rangle=\sum_{i}\psi^{B}_{0}|-\alpha|^{i}b^{\dagger}_{i}|0\rangle,$

where $\psi^{{A_{r}}({A_{l}})}_{j}$ and $\psi^{B}_{j}$ represent the
amplitudes of the wave function on the $j$th site of sub-lattice A on the
right (left) side of the impurity and on sub-lattice B, respectively. Indeed,
the reason for considering two distinct parts, $|\psi_{A}^{r}\rangle$ and
$|\psi_{A}^{l}\rangle$, in the case of the A sub-lattice components is due to
the lack of symmetry in the wave function on the left and right sides of the
impurity. The impurity introduces an asymmetry in the system, leading to
different contributions from the A sub-lattice sites on each side. Therefore,
to accurately describe the behavior of the wave function, we need to consider
separate components for the A sub-lattice sites on the right and left sides of
the impurity. Subsequently, we insert this assumed solution into the equation

$\mathcal{H}|\psi_{\text{s}}\rangle=(\mathcal{H}_{\text{SSH}}+\mathcal{H}_{\text{impurity}})|\psi_{\text{s}}\rangle=E_{s}|\psi_{\text{s}}\rangle,$
(14)

enabling us to solve for the coefficients that render the equation valid.
These considerations lead to the following set of equations for the amplitudes
of the wave function:

$\displaystyle
t_{1}\psi_{0}^{A_{l}}+t_{2}\psi_{1}^{A_{r}}=(E_{s}-U)\psi_{0}^{B},$
$\displaystyle(t_{2}-\alpha t_{1})\psi_{0}^{B}=E_{s}\psi_{1}^{A_{r}},$
$\displaystyle(t_{1}-\alpha t_{2})\psi_{1}^{A_{r}}=-\alpha E_{s}\psi_{0}^{B},$
$\displaystyle(t_{1}-\alpha t_{2})\psi_{0}^{B}=E_{s}\psi_{0}^{A_{l}},$
$\displaystyle(t_{2}-\alpha t_{1})\psi_{0}^{A_{l}}=-\alpha E_{s}\psi_{0}^{B}.$
(15)

Besides the three unknown amplitudes of the wave function
$\psi_{0}^{A_{l}},\psi_{1}^{A_{r}}$ and $\psi_{0}^{B}$ , we have two
additional unknown quantities, namely $E_{s}$ and $\alpha$. To determine $E$
and $\alpha$, we need to solve the set of five equations. However, our main
interest lies in finding $E_{s}$ and $\alpha$ directly. Therefore, we can
eliminate the other three unknowns from the equations through some algebraic
manipulations, resulting in the following expressions:

$\displaystyle t_{1}T_{1}+t_{2}T_{2}=E_{s}(E_{s}-U)$ $\displaystyle
T_{1}T_{2}=-\alpha E_{s}^{2}.$ (16)

Here, we have defined $T_{1}=t_{1}-\alpha t_{2}$ and $T_{2}=t_{2}-\alpha
t_{1}$. By solving these two equations, we can obtain the desired values of
$E_{s}$ and $\alpha$. This set of equations has two sets of solutions: one
with $E_{s}$ inside the gap region and the other with $E_{s}$ outside. We have
depicted the solutions for $E_{s}$ inside the gap region and its corresponding
$\alpha$ graphically in Fig. 4. As shown, regardless of the strength of the
impurity potential $U$, a bound state emerges within the energy gap. For
extremely large values of $U$ $(U\rightarrow\infty)$, the energy of the bound
state approaches zero $(E_{s}/t_{1}\rightarrow 0)$, and the value of $\alpha$
converges to $\frac{t_{2}}{t_{1}}$. As we decrease the impurity strength $U$,
the energy of the bound state gradually decreases until it eventually merges
with the bulk states at $U=0$.

####  III.2.2 Generalization of the analysis for the Kane-Mele ribbon embedded
with a line defect

Figure 5:  The figure illustrates the schematic representation of the two-leg
ladder model, obtained through a dimensional reduction procedure outlined in
Ref. Amini2 , for the Kane-Mele ribbon. The unit cell of the model consists of
four sites (labeled as $a,b,c,$ and $d$) denoted by a dashed rectangle,
forming a plaquette. Following the analogy with the SSH chain, each plaquette
is divided into two blocks, $A$ and $B$, shown by dotted lines, with each
block containing two sites. The hopping parameters $t^{\prime\prime}(k_{y})$
and $-t^{\prime\prime}(k_{y})$ are indicated by solid and dashed vectors,
respectively, and are defined as
$t^{\prime\prime}(k_{y})=2\lambda_{\text{SO}}\sin\big{(}\dfrac{k_{y}}{2}\big{)}$.

We can extend the previous calculations to the case of a Kane-Mele ribbon with
a line defect along the $y$-direction, as illustrated in Fig. 1. To simplify
the analysis, we assume that the ribbon is wide enough along the $y$-direction
so that the momentum $k_{y}$ can be treated as a good quantum number. This
assumption allows us to employ the method and notation introduced in Ref.
Amini2 , where the Kane-Mele Hamiltonian of Eq. (II.1) is mapped onto a two-
leg ladder system with a generalized SSH Hamiltonian in the $k_{y}$-space,
$\mathcal{H}_{KM}(k_{y})=\mathcal{H}_{0}(k_{y})+\mathcal{H}_{1}(k_{y})$. In
this ladder model, the unit cell consists of a plaquette composed of four
sites, denoted as $a$, $b$, $c$, and $d$. To establish an analogy with the SSH
chain, we label the sites $a$ and $b$ as block $A$, and the sites $c$ and $d$
as block $B$, as depicted in Fig. 5. In this notation, the transformed
Hamiltonian for the ladder system consists of two terms,
$\mathcal{H}_{0}(k_{y})$ and $\mathcal{H}_{1}(k_{y})$.
$\mathcal{H}_{0}(k_{y})$ represents a zero-energy flat band, while
$\mathcal{H}_{1}(k_{y})$ is responsible for the dispersion of the edge band.
To proceed, we first aim to obtain the impurity bound state for the ladder
system in the absence of $\mathcal{H}_{1}(k_{y})$. A schematic representation
of this ladder system is depicted in Fig. 5, and we can write the Hamiltonian
$\mathcal{H}_{0}(k_{y})$ as follows Amini1 :

$\mathcal{H}_{0}(k_{y})=\sum_{n}[a_{n}^{\dagger}b_{n}^{\dagger}]{\widehat{t}_{1}}\begin{bmatrix}c_{n}\\\
d_{n}\end{bmatrix}+\sum_{n}[c_{n}^{\dagger}d_{n}^{\dagger}]{\widehat{t}_{2}}\begin{bmatrix}a_{n+1}\\\
b_{n+1}\end{bmatrix}+h\cdot c,$ (17)

where

${\widehat{t}}_{1}=\begin{bmatrix}t&-2\lambda_{\text{SO}}\sin\big{(}\dfrac{k_{y}}{2}\big{)}\\\
2\lambda_{\text{SO}}\sin\big{(}\dfrac{k_{y}}{2}\big{)}&0\end{bmatrix},$ (18)
${\widehat{t}}_{2}=\begin{bmatrix}0&2\lambda_{\text{SO}}\sin\big{(}\dfrac{k_{y}}{2}\big{)}\\\
-2\lambda_{\text{SO}}\sin\big{(}\dfrac{k_{y}}{2}\big{)}&t\end{bmatrix},$ (19)

and the summation on $n$ runs over all the plaquettes in the ladder model.

Figure 6:  Electronic band structure of a Kane-Mele ribbon with zigzag
boundaries, consisting of $20$ zigzag chains. The ribbon includes a zigzag
chain of on-site impurities, positioned away from the edges, leading to the
formation of impurity bands. The impurity bands in the gap region are
highlighted in black, while the edge bands have not been shown to avoid
confusion.

In the presence of the line defect, the Hamiltonian given by Eq. (2) alters
the on-site potential of the corresponding block of the ladder system (based
on the location of the defected zigzag chain) to the on-site energy $U$. For
instance, in Fig. 5, the line defect is shown on block B of the $n$-th
plaquette. Thus, there exists an analogy between the SSH ladder and SSH chain,
and we can derive the corresponding set of equations given in Eq. (15) for the
ladder system by simply replacing the wave function amplitudes in the SSH
chain model as two-component vectors. That is,
$|\psi^{A_{r}}\rangle=\begin{bmatrix}\psi^{a_{r}}\\\
\psi^{b_{r}}\end{bmatrix}$,
$|\psi^{A_{l}}\rangle=\begin{bmatrix}\psi^{a_{l}}\\\
\psi^{b_{l}}\end{bmatrix}$, and $|\psi^{B}\rangle=\begin{bmatrix}\psi^{c}\\\
\psi^{d}\end{bmatrix}$. This immediately results in the following set of
equations:

$\displaystyle{\widehat{t}_{1}}^{\dagger}\begin{bmatrix}\psi^{a_{l}}_{0}\\\
\psi^{b_{l}}_{0}\end{bmatrix}+{\widehat{t}_{2}}\begin{bmatrix}\psi^{a_{r}}_{1}\\\
\psi^{b_{r}}_{1}\end{bmatrix}=(E_{s}-U)\begin{bmatrix}\psi^{c}_{0}\\\
\psi^{d}_{0}\end{bmatrix},$
$\displaystyle({\widehat{t}}_{2}^{\dagger}-\alpha{\widehat{t}}_{1})\begin{bmatrix}\psi^{c}_{0}\\\
\psi^{d}_{0}\end{bmatrix}=E_{s}\begin{bmatrix}\psi^{a_{r}}_{1}\\\
\psi^{b_{r}}_{1}\end{bmatrix},$
$\displaystyle({\widehat{t}}_{1}^{\dagger}-\alpha{\widehat{t}}_{2})\begin{bmatrix}\psi^{a_{r}}_{1}\\\
\psi^{b_{r}}_{1}\end{bmatrix}=-\alpha E_{s}\begin{bmatrix}\psi^{c}_{0}\\\
\psi^{d}_{0}\end{bmatrix},$
$\displaystyle({\widehat{t}}_{1}-\alpha{\widehat{t}}_{2}^{\dagger})\begin{bmatrix}\psi^{c}_{0}\\\
\psi^{d}_{0}\end{bmatrix}=E_{s}\begin{bmatrix}\psi^{a_{l}}_{0}\\\
\psi^{b_{l}}_{0}\end{bmatrix},$
$\displaystyle({\widehat{t}}_{2}-\alpha{\widehat{t}}_{1}^{\dagger})\begin{bmatrix}\psi^{a_{l}}_{0}\\\
\psi^{b_{l}}_{0}\end{bmatrix}=-\alpha E_{s}\begin{bmatrix}\psi^{c}_{0}\\\
\psi^{d}_{0}\end{bmatrix}.$ (20)

By the same token, this will allow us to obtain the unknown quantities
$E_{s}(k_{y})$ and $\alpha(k_{y})$ through the solution of the following
equations:

$\displaystyle{\widehat{t}}_{1}{\widehat{T}}_{1}+{\widehat{t}}_{2}{\widehat{T}}_{2}=E_{s}(E_{s}-U){\widehat{I}}$
$\displaystyle{\widehat{T}}_{1}{\widehat{T}}_{2}=-\alpha
E_{s}^{2}{\widehat{I}}.$ (21)

where ${\widehat{T}}_{1}={\widehat{t}}_{2}^{\dagger}-\alpha{\widehat{t}}_{1}$
and ${\widehat{T}}_{2}={\widehat{t}}_{1}^{\dagger}-\alpha{\widehat{t}}_{2}$.
Using the definition of ${\widehat{T}}_{1}$ and ${\widehat{T}}_{2}$, one can
simplify the second equation in the equation set (21) as
${\widehat{T}}_{1}{\widehat{T}}_{2}=(1-\alpha)^{2}{t^{\prime\prime}}^{2}(k_{y})-\alpha
t^{2}(k_{y})=-\alpha E_{s}^{2}(k_{y})$. Indeed, the resulting energy of these
bound states associated with the impurity potential $U$ explicitly depends on
the wave number $k_{y}$. This observation indicates the formation of an
impurity band Amini3 , which leads to the creation of an impurity channel
responsible for the coupling between the top and bottom edges of the ribbon.
The above-mentioned set of equations has four solutions, and two of them lie
in the gap region, which we will refer to later on. The final challenge lies
in understanding the effect of $\mathcal{H}_{1}(k_{y})$. To tackle this, we
need to employ a perturbative approach following the method described in Ref.
Amini2 . This entails using the unperturbed wave functions to obtain the
expectation values of $\mathcal{H}_{1}(k_{y})$ in between. Due to the
complexity of these calculations, analytical solutions are not readily
available, necessitating the use of numerical computations.

To obtain more accurate results and minimize perturbation errors, we can take
a different approach. We consider a zigzag ribbon of the Kane-Mele model with
a line defect, represented by a zigzag chain of on-site impurities with an
energy $U$. By numerically computing the band structure associated with this
ribbon, we can gain deeper insights. For example, we can analyze the behavior
of a ribbon with $N=20$ zigzag chains, where we set $\lambda_{\text{SO}}=0.2$
and $U=1.0$. The resulting band structure is presented in Fig. 6. As
previously mentioned, the presence of the line defect gives rise to two
additional bands within the energy gap. In Fig. 6, we have highlighted these
in-gap impurity bands using black colors for clarity. The edge bands are not
shown in the plot to avoid confusion and better focus on the behavior of the
impurity bands.

###  III.3 Impurity channel analysis

Figure 7:  The figure presents spatial profiles of the local density of states
on sites adjacent to the line defect at energies specified in Fig. 2. The
scattering region is a Kane-Mele ribbon (depicted in Fig. 1) with a width of
$W=26$ and spin-orbit coupling $\lambda_{\text{SO}}=0.2t$ that hosts a line
defect with on-site potential strength $U=1.0t$. (a) shows the LDOS at energy
$E/t=-0.88$ with complete suppression of the transmission coefficient. (b)
shows the LDOS at energy $E/t=-0.68$ with partial suppression of the
transmission coefficient. (c) displays the LDOS at $E/t=0.0$ with an unchanged
transmission coefficient.

Another approach to comprehend the formation of the impurity channel is by
investigating the LDOS at energies where topological protection is broken. By
analyzing the LDOS, which is related to the matrix elements of the retarded
Green’s function as $\rho_{i}(E)=-\frac{1}{\pi}\text{Im}[\langle
i|G_{r}(E)|i\rangle]$, we can gain insights into how the presence of the
impurity potential affects the electronic states in real space and leads to
the emergence of the impurity channel. Fig. 7 displays the spatial
distribution of the LDOS for the scattering region depicted in Fig. 1 near the
line defects at energies specified in Fig. 2. We are interested in the
energies near the lower gap edge where the presence of line defects has a
visible impact. As before, we consider the parameters
$\lambda_{\text{SO}}=0.2t$ and $U=1.0t$ for our calculations. These values
remain consistent with our previous analysis and investigations of the
system’s behavior in the presence of the line defect. Specifically, Fig. 7(a)
illustrates the LDOS map at a specific energy value of $E/t=-0.88$, showcasing
a complete suppression of transmission, indicating the formation of the
impurity channel. Additionally, Fig. 7(b) depicts the LDOS at $E/t=-0.68$,
where the transmission is not entirely suppressed. However, due to the
presence of impunity channel connecting the top and bottom edges of the ribbon
partial backscattering existence. For comparison, Fig. 7(c) shows the LDOS at
$E/t=0.0$, which exclusively displays only the presence of edge states
(without formation of any impurity channel), as anticipated. These results
provide valuable insights into the impact of the line defect on the edge
states’ transmission and the formation of the impurity channel in the system.

##  IV CONCLUDING REMARKS

In conclusion, we have investigated the impact of a line defect in the form of
on-site impurities, arranged in a zigzag chain connecting the top and bottom
edges of a Kane-Mele armchair ribbon. Our study focused on the scattering of
edge states by this line defect and its effect on the topological protection
of these states. We found that the presence of the line defect can
significantly break the topological protection of the edge states, leading to
backscattering of the helical electrons. As a consequence, the conductance of
the system becomes lower than that of the pristine ribbon. Our analysis
involved a combination of numerical and analytical calculations to elucidate
the underlying mechanism behind this backscattering phenomenon. The key result
of our study is the observation that the line defect induces an impurity
channel that spans the width of the ribbon. This impurity channel facilitates
the backscattering of helical electrons, causing a suppression of the
quantized transmission within specific energy windows. In summary, our work
sheds light on the role of line defects in disrupting the topological
protection of edge states in the Kane-Mele model. These findings have
implications for the design and understanding of topological systems and their
robustness against imperfections and defects.

###### Acknowledgements.

We extend our sincere gratitude to the office of graduate studies at the
University of Isfahan(UI) for their generous support and for providing us with
research facilities, which were instrumental in carrying out this study.
Additionally, MA would like to acknowledge the support received from the Abdus
Salam (ICTP) associateship program.

## References

  * (1) M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, Rev. Mod. Phys. 82, 3045 (2010).
  * (2) X.-L. Qi and S.-C. Zhang, Topological insulators and superconductors, Rev. Mod. Phys. 83, 1057 (2011).
  * (3) D. Culcer, A. C. Keser, Y. Li, and G. Tkachov, Transport in two-dimensional topological materials: recent developments in experiment and theory, 2D Mater. 7, 022007 (2020).
  * (4) S. Wu, V. Fatemi, Q. D. Gibson, K. Watanabe, T. Taniguchi, R. J. Cava, and P. Jarillo-Herrero, Observation of the quantum spin Hall effect up to 100 kelvin in a monolayer crystal, Science 359, 76 (2018).
  * (5) L. Lunczer, P. Leubner, M. Endres, V. L. Muller, C. Brune, H. Buhmann, and L. W. Molenkamp, Approaching Quantization in Macroscopic Quantum Spin Hall Devices through Gate Train- ing, Phys. Rev. Lett. 123, 047701 (2019).
  * (6) C.-H. Hsu, P. Stano, J. Klinovaja, and D. Loss, Helical liquids in semiconductors, Semicond. Sci. Technol. 36, 123003 (2021).
  * (7) L. Vannucci, T.Olsen, and K. S. Thygesen, Conductance of quantum spin Hall edge states from first principles: The critical role of magnetic impurities and inter-edge scattering, Phys. Rev. B 101 (2020).
  * (8) L.R. F. Lima and C. Lewenkopf, Breakdown of topological protection due to nonmagnetic edge disorder in two-dimensional materials in the quantum spin Hall phase ,Phys. Rev. B 106, 245408 (2022).
  * (9) P. Novelli, F. Taddei, A. K. Geim, and M. Polini, Failure of Conductance Quantization in Two-Dimensional Topological Insulators due to Nonmagnetic Impurities, Phys. Rev. Lett. 122, 016601 (2019).
  * (10) B. Jäck, Y. Xie, B. A. Bernevig, and A. Yazdani, Observation of backscattering induced by magnetism in a topological edge state. Proc. Natl Acad. Sci. USA 117, 16214 (2020).
  * (11) J. L. Collins, et al. Electric-field-tuned topological phase transition in ultra-thin Na3Bi, Nature 564, 390 (2018).
  * (12) X. Qian, J. Liu, L. Fu, and J. Li, Quantum spin Hall effect in two-dimensional transition metal dichalcogenides, Science 346, 1344 (2014).
  * (13) W. Vandenberghe, and M. Fischetti, Imperfect two-dimensional topological insulator field-effect transistors, Nat. Commun. 8, 14184 (2017).
  * (14) V. Krueckl, and K. Richter, Switching spin and charge between edge states in topological insulator constrictions, Phys. Rev. Lett. 107, 086803 (2011).
  * (15) R. P. Maciel, A. L. Araújo, C. H. Lewenkopf, and G. J. Ferreira, Fabry-Pérot resonant vortices and magnetoconductance in topological insulator constrictions with magnetic barriers. Phys. Rev. B 103, 205124 (2021).
  * (16) H. Ishida, and A. Liebsch, Engineering edge-state currents at the interface between narrow ribbons of two-dimensional topological insulators, Phys. Rev. Res. 2, 023242 (2020).
  * (17) L. N. Zheng, L. Qi, L. Y. Cheng, H. F. Wang, and S. Zhang, Defect-induced controllable quantum state transfer via a topologically protected channel in a flux qubit chain, Phys. Rev. A 102, 012606 (2020).
  * (18) M. Konig, S. Wiedmann, C. Brune, A. Roth, H. Buhmann, L.W. Molenkamp, X.-L. Qi, and S.-C. Zhang, Science 318, 5851 (2007).
  * (19) J. I. Vayrynen, M. Goldstein, and L. I. Glazman, Helical Edge Resistance Introduced by Charge Puddles, Phys. Rev. Lett. 110, 216402 (2013).
  * (20) C.-H. Hsu, P. Stano, J. Klinovaja, and D. Loss, Effects of nuclear spins on the transport properties of the edge of two-dimensional topological insulators, Phys. Rev. B 97, 125432 (2018).
  * (21) A. A. Bagrov, F. Guinea, and M. I. Katsnelson, Suppressing backscattering of helical edge modes with a spin bath, Phys. Rev. B 100, 195426 (2019).
  * (22) J. C. Budich, F. Dolcini, P. Recher, and B. Trauzettel, Phonon-Induced Backscattering in Helical Edge States, Phys. Rev. Lett. 108, 086602 (2012).
  * (23) J. I. Vayrynen, D. I. Pikulin, and J. Alicea, Noise-Induced Backscattering in a Quantum Spin Hall Edge, Phys. Rev. Lett. 121, 106601 (2018).
  * (24) B. Zhou, H.-Z. Lu, R.-L. Chu, S. Q. Shen, and Q. Niu, Finite size effects on helical edge states in a quantum spin-Hall system. Phys. Rev. Lett. 101, 246807 (2008).
  * (25) F. Romeo, R. Citro, D. Ferraro, and M. Sassetti, Electrical switching and interferometry of massive Dirac particles in topological insulator constrictions. Phys. Rev. B 86, 165418 (2012).
  * (26) Takagaki, Y. Backscattering from width variations in quasi-one-dimensional strips of topological insulators. J. Phys. Condens. Matter 24, 435301 (2012).
  * (27) Takagaki, Y. Strong reflection and periodic resonant transmission of helical edge states in topological-insulator stub-like resonators. J. Appl. Phys. 118, 054304 (2015).
  * (28) M. Sadeghizadeh, M. Soltani, and M. Amini, Rigorous analysis of the topologically protected edge states in the quantum spin Hall phase of the armchair ribbon geometry, arXiv:2306.03690 (2023).
  * (29) E. N. Lima, T. Schmidt, and R. W. Nunes, Topologically protected metallic states induced by a one-dimensional extended defect in the bulk of a 2D topological insulator. Nano Lett. 16, 4025 (2016).
  * (30) R. Stuhler, A. Kowalewski, F. Reis, D. Jungblut, F. Dominguez, B. Scharf, G. Li, J. Schafer, E. M. Hankiewicz, and R. Claessen Nature Communications 13, 3480 (2022).
  * (31) S. Datta, Electronic Transport in Mesoscopic Systems (Cam- bridge University Press, Cambridge, 1995).
  * (32) C. L. Kane and E. J. Mele, $Z_{2}$ Topological Order and the Quantum Spin Hall, Effect, Phys. Rev. Lett. 95, 146802 (2005).
  * (33) C. L. Kane and E. J. Mele, Quantum Spin Hall Effect in Graphene, Phys. Rev. Lett. 95, 226801 (2005).
  * (34) W. P. Su, J. R. Schrieffer, and A. J. Heeger, Phys. Rev. Lett. 42, 1698 (1979).
  * (35) pyqula Library https://github.com/joselado/pyqula.
  * (36) J. K. Asboth, L. Oroszlany, & A. Palyi, A Short Course on Topological Insulators: Band Structure and Edge States in One and Two Dimensions. Springer International Publishing (2016).
  * (37) F. Rahmati, M. Amini, M. Soltani, & M. Sadeghizadeh, Explicit derivation of the chiral and generic helical edge states for the Kane-Mele model: Closed expressions for the wave function, dispersion relation, and spin rotation, Phys. Rev. B 107 (20), 205408 (2023).
  * (38) M. Rezaei, H. Karbaschi, M. Amini, M. Soltani, & G. Rashedi, Thermoelectric properties of armchair phosphorene nanoribbons in the presence of vacancy-induced impurity band, Nanotechnology 32 (37), 375704 (2021).