# A Generalized Neural Diffusion Framework on Graphs

Yibo Li1, Xiao Wang2, Hongrui Liu3, Chuan Shi1 Corresponding author.

###### Abstract

Recent studies reveal the connection between GNNs and the diffusion process,
which motivates many diffusion-based GNNs to be proposed. However, since these
two mechanisms are closely related, one fundamental question naturally arises:
Is there a general diffusion framework that can formally unify these GNNs? The
answer to this question can not only deepen our understanding of the learning
process of GNNs, but also may open a new door to design a broad new class of
GNNs. In this paper, we propose a general diffusion equation framework with
the fidelity term, which formally establishes the relationship between the
diffusion process with more GNNs. Meanwhile, with this framework, we identify
one characteristic of graph diffusion networks, i.e., the current neural
diffusion process only corresponds to the first-order diffusion equation.
However, by an experimental investigation, we show that the labels of high-
order neighbors actually exhibit monophily property, which induces the
similarity based on labels among high-order neighbors without requiring the
similarity among first-order neighbors. This discovery motives to design a new
high-order neighbor-aware diffusion equation, and derive a new type of graph
diffusion network (HiD-Net) based on the framework. With the high-order
diffusion equation, HiD-Net is more robust against attacks and works on both
homophily and heterophily graphs. We not only theoretically analyze the
relation between HiD-Net with high-order random walk, but also provide a
theoretical convergence guarantee. Extensive experimental results well
demonstrate the effectiveness of HiD-Net over state-of-the-art graph diffusion
networks.

## Introduction

Graphs, such as traffic networks, social networks, citation networks, and
molecular networks, are ubiquitous in the real world. Recently, Graph Neural
Networks (GNNs), which are able to effectively learn the node representations
based on the message-passing manner, have shown great popularity in tackling
graph analytics problems. So far, GNNs have significantly promoted the
development of graph analysis towards real-world applications. e.g, node
classifification (Abu-El-Haija et al. 2019; Wu, He, and Xu 2019), link
prediction (Kipf and Welling 2016b; You, Ying, and Leskovec 2019), subgraph
isomorphism counting (Yu et al. 2023), and graph classifification (Gao and Ji
2019; Zhang et al. 2018).

Some recent studies show that GNNs are in fact intimately connected to
diffusion equations (Chamberlain et al. 2021; Wang et al. 2021; Thorpe et al.
2022), which can be considered as information diffusion on graphs. Diffusion
equation interprets GNNs from a continuous perspective (Chamberlain et al.
2021) and provides new insights to understand existing GNN architectures,
which motives some diffusion-based GNNs. For instance, (Wang et al. 2021)
proposes continuous graph diffusion. (Thorpe et al. 2022) utilizes the
diffusion process to handle oversmoothing issue. (Song et al. 2022) considers
a graph as a discretization of a Riemannian manifold and studies the
robustness of the information propagation process on graphs. Diffusion
equation can also build a bridge between traditional GNNs and control theory
(Zang and Wang 2020). Although the diffusion process and graph convolution are
closely related, little effort has been made to answer: Is there a unified
diffusion equation framework to formally unify the current GNN architectures?
A well-informed answer can deepen our understanding of the learning mechanism
of GNNs, and may inspire to design a broad new class of GNNs based on
diffusion equation.

Actually, (Chamberlain et al. 2021) has explained GCN (Kipf and Welling 2016a)
and GAT (Veličković et al. 2017) from diffusion equation. However, with more
proposed GNN architectures, it is highly desired to formally revisit the
relation between diffusion equation and GNNs. In this paper, we discover that
many GNN architectures substantially can be unified with a general diffusion
equation with the fidelity term, such as GCN/SGC [22], APPNP[14], GAT [19],
AMP [16], DAGNN [15]. Basically, the diffusion equation describes that the
change of a node representation depends on the movement of information on
graphs from one node to its neighbors, and the fidelity term constraints that
the change of a node representation depends on the difference with its initial
feature. Furthermore, we show that the unified diffusion framework can also be
derived from an energy function, which explains the whole framework as an
energy minimization process in a global view. Compared with other unified
frameworks (Zhu et al. 2021; Ma et al. 2021), our framework is from the
diffusion perspective, which has many advantages. For example, diffusion-based
methods are able to address the common plights of graph learning models such
as oversmoothing (Chamberlain et al. 2021; Wang et al. 2021). What’s more, the
diffusion equation can be seen as partial differential equations (PDEs)
(Chamberlain et al. 2021), thus introducing many schemes to solve the graph
diffusion equation such as explicit scheme, implicit scheme, and multi-step
scheme, some of which are more stable and converge faster.

Based on the above findings, we can see that the diffusion process employed by
most current GNNs just considers the first-order diffusion equation, which
only diffuses messages among 1-hop neighbors. That is, the first-order
diffusion has the underlying homophily assumption among 1-hop neighbors. While
we empirically discover that the labels of 2-hop neighborhoods actually appear
monophily property (Altenburger and Ugander 2018), i.e., nodes may have
extreme preferences for a particular attribute which are unrelated to their
own attribute and 1-hop neighbors’ attribute, but are more likely to be
similar with the attribute of their 2-hop neighbors. Simply put, monophily can
induce a similarity among 2-hop neighbors without requiring similarity among
1-hop neighbors. So when the 1-hop neighbors are heterophily-dominant or have
noise, the 2-hop neighbors will provide more relevant context. Therefore, a
more practical diffusion process should take both the first-order and second-
order neighbors into account. How can we design a new type of graph diffusion
networks satisfying the above requirement based on our framework?

In this paper, we design a new high-order neighbor-aware diffusion equation in
our proposed diffusion framework, and then derive a High-order Graph Diffusion
Network (HiD-Net). Specifically, our model simultaneously combines the first-
order and second-order diffusion process, then we regularize the diffusion
equation by minimizing the discrepancy between the estimated and the original
graph features. The whole diffusion equation is finally integrated into the
APPNP architecture. With second-order diffusion equation, HiD-Net is more
robust against the attacks and more general on both homophily and heterophily
graphs. We theoretically prove that HiD-Net is essentially related with the
second-order random walk. We also provide the convergence guarantee that HiD-
Net will converge to this random walk’s limit distribution as the number of
layers increases, and meanwhile, the learned representations do not converge
to the same vector over all nodes. The contributions of this paper are
summarized as follows:

  * •

We propose a novel generalized diffusion graph framework, consisting of
diffusion equation and fidelity term. This framework, formally establishing
the relation between diffusion process with a wide variety of GNNs, describes
a broad new class of GNNs based on the discretized diffusion equations on
graphs and provides new insight to the current graph diffusion/neural
networks.

  * •

We discover the monophily property of labels, and based on our diffusion
framework, we propose a high-order graph diffusion network, HiD-Net, which is
more general and robust. We theoretically build the relation between HiD-Net
and second-order random walk, together with the convergence guarantee.

  * •

Our extensive experiments on both the homophily and heterophily graphs clearly
show that HiD-Net outperforms the popular GNNs based on diffusion equation.

## Related Work

Graph convolutional networks. Recently, graph convolutional network (GCN)
models (Bruna et al. 2013; Defferrard, Bresson, and Vandergheynst 2016; Kipf
and Welling 2016a; Veličković et al. 2017; Hamilton, Ying, and Leskovec 2017)
have been widely studied. Based on the spectrum of graph Laplacian, (Bruna et
al. 2013) generalizes CNNs to graph signal. Then (Defferrard, Bresson, and
Vandergheynst 2016) further improves the efficiency by employing the Chebyshev
expansion of the graph Laplacian. (Kipf and Welling 2016a) proposes to only
aggregate the node features from the one-hop neighbors and simplifies the
convolution operation. (Veličković et al. 2017) introduces the attention
mechanisms to learn aggregation weights adaptively. (Hamilton, Ying, and
Leskovec 2017) uses various ways of pooling for aggregation. More works on
GNNs can be found in surveys (Wu et al. 2020b; Zhou et al. 2020).

Diffusion equation on graphs. Graph Heat Equation (GHE) (Chung and Graham
1997), which is a well-known generalization of the diffusion equation on graph
data, models graph dynamics with applications in spectral graph theory. GRAND
(Chamberlain et al. 2021) studies the discretized diffusion PDE on graphs and
applies different numerical schemes for their solution. GRAND++ (Thorpe et al.
2022) mitigates the oversmoothing issue of graph neural networks by adding a
source term. DGC (Wang et al. 2021) decouples the terminal time and
propagation steps of linear GCNs from a perspective of graph diffusion
equation, and analyzes why linear GCNs fail to benefit from deep layers. ADC
(Zhao et al. 2021) strategies to automatically learn the optimal diffusion
time from the data. However, these works focus on specific graph diffusion
network, thus there is not a framework to formally unify the GNNs.

The unified GNN framework. (Zhu et al. 2021) establishes a connection between
different propagation mechanisms with a unified optimization problem, and
finds out that the proposed propagation mechanisms are the optimal solution
for optimizing a feature fitting function over a wide class of graph kernels
with a graph regularization term. (Ma et al. 2021) establishes the connections
between the introduced GNN models and a graph signal denoising problem with
Laplacian regularization. It essentially is still an optimization solving
framework from the perspective of signal denoising. However, our framework is
based on the diffusion equation, where the advantages are two fold: one is
that diffusion-based methods are able to address the oversmoothing problem
(Chamberlain et al. 2021; Wang et al. 2021). The other is that the diffusion
equation can be seen as partial differential equations (PDEs) (Chamberlain et
al. 2021) and thus can introduce many schemes that have many good properties
such as fast convergence rate and high stability.

## The Generalized Diffusion Graph Framework

Notations. Consider an undirected graph as $G=(V,E)$ with adjacency matrix
$\mathbf{A}\in\mathbb{R}^{n\times n}$, where $V$ contains $n$ nodes
$\\{v_{1},\dots,v_{n}\\}$ and $E$ is the set of edges. The initial node
feature matrix is denoted as $\mathbf{X}^{(0)}\in\mathbb{R}^{n\times q}$,
where $q$ is the dimension of node feature. We denote the neighbors of node
$i$ at exactly $k$ hops/steps away as $N_{k}(i)$. For example,
$N_{1}(i)=\\{j:(i,j)\in E\\}$ are the immediate neighbors of $i$.

Diffusion is a physical process that equilibrates concentration differences
without creating or destroying mass. This physical observation can be easily
cast in the diffusion equation, which is a parabolic partial differential
equation. Fick’s law of diffusion describes the equilibration property
(Weickert 1998):

$J=-f\cdot\nabla u,$ (1)

where $J$ is the diffusion flux, which measures the amount of substance that
flows through a unit area during a unit time interval. $f$ is the diffusivity
coefficient, which can be constant or depend on time and position. $\nabla u$
is the concentration gradient. This equation states that a concentration
gradient $\nabla u$ causes a flux $J$ which aims to compensate for this
gradient. The observation that a change in concentration in any part of the
system is due to the influx and outflux of substance into and out of that part
of the system can be expressed by the continuity equation:

$\dfrac{\partial u}{\partial t}=-\operatorname{div}J,$ (2)

where $t$ denotes the time. Plugging in Fick’s law (1) into the continuity
equation we end up with the diffusion equation:

$\dfrac{\partial u}{\partial t}=\operatorname{div}(f\cdot\nabla u).$ (3)

As $\operatorname{div}$ is the sum of the second derivatives in all
directions, please note that normal first order derivatives and second order
derivatives are on continuous space and can not be generalized directly to
graph which is on discrete space. As (Chamberlain et al. 2021) defined, the
first derivative is the difference between the feature of a node and its
neighbor. And the second derivative can be considered as the difference
between the first derivatives of the node itself and its neighbors.

For better illustration, we provide an example. Consider a chain graph in
Figure 1, where $i$, $i+1$, and $i-1$ are the indexes of the nodes. The
feature of node $i$ is denoted as $x_{i}$.

Figure 1: Chain graph

The first order derivatives on node $i$ is defined as $x_{i+1}-x_{i}$ and
$x_{i}-x_{i-1}$. The diffusion flux from node $j$ to node $i$ at time $t$ on a
graph is:

$J_{ij}^{(t)}=-f\cdot(\nabla x)_{ij}^{(t)}=-f(x_{j}^{(t)}-x_{i}^{(t)}).$ (4)

The second order derivative is the difference of first order derivatives:
$(x_{i+1}-x_{i})-(x_{i}-x_{i-1})=(x_{i+1}-x_{i})+(x_{i-1}-x_{i})$. We notice
that the chain graph only has one dimension, so the divergence of node $i$ on
a chain graph is equal to its second order derivative:
$\operatorname{div}(\nabla x_{i})=(x_{i+1}-x_{i})+(x_{i-1}-x_{i})$.

Thus, on a normal graph, we have the generalized form:
$\operatorname{div}(\nabla x_{i})=\sum_{j\in N_{1}(i)}(x_{j}-x_{i})$.

Here we normalize the diffusion process utilizing the degree of the nodes to
down-weight the high-degree neighbors, and we have $\operatorname{div}(\nabla
x_{i})=\sum_{j\in
N_{1}(i)}\dfrac{\tilde{A}_{ij}}{\sqrt{\tilde{d_{i}}}\sqrt{\tilde{d_{j}}}}(x_{j}-x_{i})$,
where $\tilde{A}_{ij}$ is the element of
$\tilde{\mathbf{A}}=\mathbf{A}+\mathbf{I}$, and
$d_{i}=\sum_{j}\tilde{A}_{ij}$.

So the diffusion equation on node $i$ can be defined as:

$\displaystyle\frac{\partial x_{i}^{(t)}}{\partial t}$
$\displaystyle=-\operatorname{div}J_{ij}^{(t)}=\operatorname{div}[f(\nabla
x)_{ij}^{(t)}]$ (5) $\displaystyle=f\sum_{j\in
N_{1}(i)}\dfrac{\tilde{A}_{ij}}{\sqrt{\tilde{d_{i}}}\sqrt{\tilde{d_{j}}}}(x_{j}^{(t)}-x_{i}^{(t)}).$

The diffusion equation models the change of representation $x_{i}^{(t)}$ with
respect to $t$, which depends on the difference between the nearby nodes,
implying that the greater the difference between a node and its neighbors, the
faster it changes.

However, how fast $x_{i}^{(t)}$ changes should not only depend on the
representation difference between node $i$ and its neighbors, otherwise, it
will cause oversmoothing issue, i.e., as the diffusion process goes by, the
nodes are not distinguishable. Based on this phenomenon, we think that the
representation change of $x_{i}^{(t)}$ should be also related with the node
feature $x_{i}^{(0)}$ itself, i.e., if the difference between $x_{i}^{(t)}$
and $x_{i}^{(0)}$ is small, the change of $x_{i}^{(t)}$ should also be small.
Then we add another fidelity term and obtain our general graph diffusion
framework as follows:

$\dfrac{\partial x_{i}^{(t)}}{\partial
t}=\alpha(x_{i}^{(0)}-x_{i}^{(t)})+\beta\operatorname{div}(f(\nabla
x)_{ij}^{(t)}),$ (6)

where $\alpha$, $\beta$ are coefficients.

Remark 1. (6) can be derived from the energy function:

$E(x):=\int_{\Omega}\left(\alpha\cdot(x_{i}-x_{i}^{(0)})^{2}+\beta\cdot|f(\nabla
x)_{ij}|^{2}\right)d\theta,$ (7)

where $\theta$ represents the position of the nodes, and $\Omega$ represents
the entire graph domain. The corresponding Euler–Lagrange equation, which
gives the necessary condition for an extremum of (7), is given by:

$0=\alpha(x_{i}-x_{i}^{(0)})+\beta\operatorname{div}(f(\nabla x)_{ij}).$ (8)

(8) can also be regarded as the steady-state equation of (6). Based on the
energy function, we can see that (6) constrains space variation and time
variation of the diffusion process, indicating that the representations of the
graph nodes will not change too much between nearby nodes, as well as not
change too much from the initial features.

Remark 2. The framework (6) is closely related to many GNNs, such as GCN/SGC
(Wu et al. 2019), APPNP (Klicpera, Bojchevski, and Günnemann 2018), GAT
(Veličković et al. 2017), AMP (Liu et al. 2021), DAGNN (Liu, Gao, and Ji
2020), as demonstrated by the following propositions. We provide the proofs of
all the subsequent propositions in Appendix.

Proposition 1. With $\alpha=0$, $\beta=1$, $\Delta t=1$ and $f=1$ in (6), the
diffusion process in SGC/GCN is:

$\dfrac{\partial x_{i}^{(t)}}{\partial t}=\operatorname{div}((\nabla
x)_{ij}^{(t)}).$ (9)

Proposition 2. Introducing $\eta$ as coefficient, with $\alpha=1$,
$\beta=1-\dfrac{1}{\eta}$, $\Delta t=1$ and $f=1$ in (6), the diffusion
process in APPNP is:

$\dfrac{\partial x_{i}^{(t)}}{\partial
t}=(x_{i}^{(0)}-x_{i}^{(t)})+(1-\dfrac{1}{\eta})\operatorname{div}((\nabla
x)_{ij}^{(t)}).$ (10)

Proposition 3. With $\alpha=0$, $\beta=1$, $\Delta t=1$ and the learned
similarity coefficient $f_{ij}^{(t)}$ between nodes $i$ and $j$ at time $t$ in
(6), the diffusion process in GAT is:

$\dfrac{\partial x_{i}^{(t)}}{\partial
t}=\operatorname{div}(f_{ij}^{(t)}(\nabla x)_{ij}^{(t)}).$ (11)

Proposition 4. With stepsize $\epsilon$ and coefficient $\lambda$,
$\beta_{i}^{(t)}=\max\left(1-\frac{\epsilon\lambda}{\left\|(1-2\epsilon(1-\lambda))\mathbf{X}_{i}^{(t)}+2\epsilon(1-\lambda)\tilde{\mathbf{A}}\mathbf{X}_{i}^{(t)}-\left(\mathbf{X}^{(0)}\right)_{i}\right\|_{2}},0\right)$.
Let $\alpha=1-\beta_{i}^{(t)}$, $\beta=2\epsilon(1-\lambda)\beta_{i}^{(t)}$,
$\Delta t=1$ and $f=1$ in (6), the diffusion process in AMP is:

$\dfrac{\partial x_{i}^{(t)}}{\partial
t}=(1-\beta_{i}^{(t)})(x_{i}^{(0)}-x_{i}^{(t)})+2\epsilon(1-\lambda)\beta_{i}^{(t)}\operatorname{div}((\nabla
x)_{ij}^{(t)}).$ (12)

Proposition 5. With $\alpha+\beta=1$, $\Delta t=1$ and the learned
coeffiecient $f(t)$ at time $t$ satisfying $\sum_{t=0}^{T}(\beta
f(t))^{t}=\frac{1}{\alpha}$ in (6), the diffusion process in DAGNN is:

$\dfrac{\partial x_{i}^{(t)}}{\partial
t}=\alpha(x_{i}^{(0)}-x_{i}^{(t)})+\beta\operatorname{div}(f(t)((\nabla
x)_{ij}^{(t)})).$ (13)

## The High-order Graph Diffusion Network

### High-order Graph Diffusion Equation

Figure 2: Illustration of the same node pair in different contexts.

In the first-order diffusion process, the diffusion equation only considers
1-hop neighbors. As shown in Figure 2, the nodes in (a), (b), and (c) are the
same, but the structures are different. The first-order diffusion flux from
node $j$ to node $i$ will be the same, even if the local structure of node $i$
and $j$ is very different. Based on the specific local environments, the
diffusion flux should be either different, so as to provide more additional
information and make the learned representations more discriminative. To
better understand the effect of local structures, we conduct an experiment on
six widely used graphs to evaluate the effect of 2-hop neighbors. First, we
have the following definition of $k$-hop neighbor similarity score.

Definition 1. Let $y_{i}$ be the label of node $i$, the $k$-hop neighbor
similarity score $h_{k}=\dfrac{|\sum_{i\in
V}\mathbf{1}_{y_{i}=\mathbf{O}(\\{y_{j},j\in N_{k}(i)\\})}|}{|N_{k}(i)|}$, and
$h_{a+b}=\dfrac{|\sum_{i\in
V}\mathbf{1}_{y_{i}=\mathbf{O}(\\{y_{j},j\in\\{N_{a}(i),N_{b}(i)\\}\\})}|}{|\\{N_{a}(i),N_{b}(i)\\}|}$,
where $\mathbf{O}(\\{y_{j},j\in N_{k}(i)\\})$ represents the element with the
highest frequency in $\\{y_{j},j\in N_{k}(i)\\}$.

The similarity score is based on node labels, and higher similarity score
implies the labels of a node and its $k$-hop neighbors are more consistent.
The scores of the six graphs are shown in Table 1. Interestingly, we find that
the labels of 2-hop neighbors show monophily property (Altenburger and Ugander
2018), i.e., as can be seen from both the homophily graphs (Cora, Citeseer,
Pubmed) and heterophily graphs (Chameleon, Squirrel, Actor), without requiring
the similarity among first-order neighbors, the second-order neighbors are
more likely to have the same labels.

Table 1: The similarity scores of six graphs.

| cora | citeseer | pubmed | chameleon | squirrel | actor
---|---|---|---|---|---|---
$h_{1}$ | 0.8634 | 0.7385 | 0.7920 | 0.2530 | 0.1459 | 0.2287
$h_{2}$ | 0.8696 | 0.8476 | 0.7885 | 0.3131 | 0.1600 | 0.3716
$h_{1+2}$ | 0.8737 | 0.8206 | 0.7880 | 0.3070 | 0.1530 | 0.3363

To take advantage of 2-hop neighbors, we regularize the gradient $\nabla
x_{i}$ utilizing the average gradient of 1-hop neighbors:

$\overline{(\nabla x)_{j}}=\operatorname{avg}(\nabla x_{jk})=\sum_{k\in
N_{1}(j)}\dfrac{\tilde{A}_{jk}}{\sqrt{\tilde{d_{j}}}\sqrt{\tilde{d_{k}}}}(x_{k}-x_{j}).$
(14)

We propose the high-order graph diffusion equation:

$\dfrac{\partial x_{i}^{(t)}}{\partial
t}=\alpha(x_{i}^{(0)}-x_{i}^{(t)})+\beta\operatorname{div}(f(\nabla
x_{ij}^{(t)})+\gamma\overline{(\nabla x)_{j}^{(t)}}),$ (15)

where $\gamma$ is the parameter of the regularization term. The iteration step
of (15) is:

$\begin{split}x_{i}^{t+\Delta t}&=\alpha\Delta tx_{i}^{(0)}+(1-\alpha\Delta
t)x_{i}^{(t)}\\\ &+\beta\Delta t\operatorname{div}(f(\nabla
x_{i}^{(t)}))+\beta\gamma\Delta t\operatorname{div}(\overline{(\nabla
x)_{j}^{(t)}}),\end{split}$ (16)

which is the diffusion-based message passing scheme (DMP) of our model. We can
see that DMP utilizes the 2-hop neighbors’ information, where the advantages
are two-fold: one is that the 2-hop neighbors capture the local environment
around a node, even if there are some abnormal features among 1-hop neighbors,
their negative effect can still be alleviated by considering a larger
neighborhood size, making the learning process more robust. The other is that
the monophily property of 2-hop neighbors provides additional stronger
correlation with labels, thus even if the 1-hop neighbors may be heterophily,
DMP can still make better predictions with information diffused from 2-hop
neighbors.

Comparison with other GNN models. Though existing GNN iteration steps can
capture high-order connectivity through iterative adjacent message passing,
they still have their limitations while having the same time complexity as
DMP. DMP is superior because it can utilize the monophily property, adjust the
balance between first-order and second-order neighbors, and is based on
diffusion equation which has some unique characteristics. More comparisons are
discussed in Appendix.

### Theoretical Analysis

Next, we theoretically analyze some properties of our diffusion-based message
passing scheme.

Definition 2. Consider a surfer walks from node $j$ to $i$ with probability
$P_{ij}$. Let $\mathbb{X}_{t}$ be a random variable representing the node
visited by the surfer at time $t$. The probability $P_{ij}$ can be represented
as a conditional probability
$\mathbb{P}\left[\mathbb{X}_{t}=i\mid\mathbb{X}_{t-\Delta t}=j\right]$. Let

$P_{ij}=\left\\{\begin{array}[]{lr}1-(\alpha+\beta)\Delta t,&i=j\\\
(\beta-\beta\gamma)\Delta t\hat{\tilde{A}}_{ij},&j\in N_{1}(i)\\\
\beta\gamma\Delta tB_{ij},&j\in N_{2}(i)\\\ \alpha\Delta t,&restart,\\\
\end{array}\right.$ (17)

where
$\hat{\tilde{\mathbf{A}}}=\tilde{\mathbf{D}}^{-\frac{1}{2}}\tilde{\mathbf{A}}\tilde{\mathbf{D}}^{-\frac{1}{2}}$,
$B_{ij}$ is the element of $\mathbf{B}=\hat{\tilde{\mathbf{A}}}^{2}$, and
restart means that the node $i$ will teleport back to the initial root node
$i$. Based on the definition, we have the following propositions.

Proposition 6. Given the probability
$\mathbb{H}^{(t)}_{ij}=\mathbb{P}\left[\mathbb{X}_{t}=i\mid\mathbb{X}_{0}=j\right]$,
DMP (16) is equivalent to the second-order random walk with the transition
probability $P_{ij}$ in (17):

$x_{i}^{(t)}=\sum_{j\in V}\mathbb{H}^{(t)}_{ij}x_{j}^{(0)}.$ (18)

Proposition 7. With $f=1$, $\alpha,\beta,\gamma,\Delta t\in\left(0,1\right]$,
DMP (16) converges, i.e., when $t\rightarrow\infty$,
$\mathbf{X}^{(\infty)}=\alpha((\alpha+\beta)\mathbf{I}-\beta(1-\gamma)\hat{\tilde{\mathbf{A}}}-\beta\gamma\hat{\tilde{\mathbf{A}}}^{2})^{-1}\mathbf{X}^{(0)}$.

Proposition 8. When $t\rightarrow\infty$, the representations of any two nodes
on the graph will not be the same as long as the two nodes have different
initial features, i.e., $\forall i,j\in V$, if $x_{i}^{(0)}\neq x_{j}^{(0)}$,
then $x_{i}^{(t)}\neq x_{j}^{(t)}$ as $t\rightarrow\infty$.

The proofs of the above propositions are in Appendix.

### Our Proposed HiD-Net

To incorporate the high-order graph diffusion DMP (16) into deep neural
networks, we introduce High-Order Graph Diffusion Network (HiD-Net). In this
work, we follow the decoupled way as proposed in APPNP (Klicpera, Bojchevski,
and Günnemann 2018):

$\mathbf{Y}^{\prime}=\mathbf{DMP}\left(l_{\omega}\left(\mathbf{X}^{(0)}\right),t,\Delta
t,\alpha,\beta,\gamma\right),$ (19)

where $l_{\omega}$ is a representation learning model such as an MLP network,
$\omega$ is the learnable parameters in the model. The training objective is
to minimize the cross entropy loss defined by the final prediction
$\mathbf{Y}^{\prime}$ and labels for training data. Because of DMP, HiD-Net is
more robust and works well on both homophily and heterophily graphs in
comparison with other graph diffusion networks.

Time complexity. The time complexity of HiD-Net can be optimized as
$\mathcal{O}(n^{2}\zeta)$, which is the same as the propagation step of GCN,
where $n$ is the number of the nodes, $\zeta$ is the dimension of the feature
vector. We provide the proof in Appendix.

## Experiments

### Node Classification

Datasets.

For comprehensive comparison, we use seven real-world datasets to evaluate the
performance of node classification. They are three citation graphs, i.e.,
Cora, Citeseer, Pubmed (Kipf and Welling 2016a), two Wikipedia networks, i.e.,
Chameleon and Squirrel (Pei et al. 2020), one Actor co-occurrence network
Actor (Pei et al. 2020), one Open Graph Benchmark(OGB) graph ogbn-arxiv(Hu et
al. 2020). Among the seven datasets, Cora, Citeseer, Pubmed and ogbn-arxiv are
homophily graphs, Chameleon, Squirrel, and Actor are heterophily graphs.
Details of datasets are in Appendix.

Table 2:  Quantitative results (%$\pm\sigma$) on node classification. (bold: best) Datasets | Metric | GCN | GAT | APPNP | GRAND | GRAND++ | DGC | ADC | HiD-Net
---|---|---|---|---|---|---|---|---|---
Cora | F1-macro | 81.5$\pm 0.6$ | 79.7$\pm 0.4$ | 82.2 $\pm 0.5$ | 79.4$\pm 2.4$ | 81.3$\pm 3.3$ | 82.1$\pm 0.1$ | 80.0$\pm 1.0$ | 82.8$\mathbf{\pm 0.6}$
F1-micro | 82.5 $\pm 0.6$ | 80.1$\pm 0.8$ | 83.2$\pm 0.2$ | 80.1$\pm 2.7$ | 82.95 $\pm 1.4$ | 83.1$\pm 0.1$ | 81.0$\pm 0.7$ | 84.0$\mathbf{\pm 0.6}$
AUC | 97.3 $\pm 0.1$ | 96.4 $\pm 0.5$ | 97.5$\pm 0.1$ | 96.0$\pm 0.3$ | 97.3$\pm 0.5$ | 97.2$\pm 0.0$ | 97.1$\pm 0.1$ | 97.6$\mathbf{\pm 0.0}$
Citeseer | F1-macro | 66.4 $\pm 0.4$ | 68.5 $\pm 0.3$ | 67.7$\pm 0.6$ | 64.9$\pm 1.5$ | 66.4$\pm 2.6$ | 68.3$\pm 0.4$ | 47.0$\pm 1.4$ | 69.5$\mathbf{\pm 0.6}$
F1-micro | 69.9 $\pm 0.5$ | 72.2 $\pm 0.3$ | 71.0 $\pm 0.4$ | 68.6$\pm 1.7$ | 70.9$\pm 2.3$ | 72.5$\pm 0.4$ | 53.7$\pm 1.5$ | 73.2$\mathbf{\pm 0.2}$
AUC | 89.9 $\pm 0.4$ | 90.2 $\pm 0.1$ | 90.3 $\pm 0.0$ | 89.5$\pm 0.8$ | 91.2$\pm 2.8$ | 91.0 $\pm 0.0$ | 87.1$\pm 1.1$ | 91.5$\mathbf{\pm 0.1}$
Pubmed | F1-macro | 78.4 $\pm 0.2$ | 76.7 $\pm 0.5$ | 79.3$\pm 0.2$ | 77.5$\pm 3.2$ | 78.9$\pm 2.5$ | 78.4$\pm 0.1$ | 73.7$\pm 2.3$ | 80.1$\mathbf{\pm 0.1}$
F1-micro | 79.1 $\pm 0.4$ | 77.3 $\pm 0.4$ | 79.9$\pm 0.3$ | 78.0$\pm 3.2$ | 79.8$\pm 1.6$ | 79.2$\pm 0.1$ | 74.3$\pm 2.3$ | 81.1$\mathbf{\pm 0.1}$
AUC | 91.2 $\pm 0.2$ | 90.3$\pm 0.5$ | 92.2$\mathbf{\pm 0.1}$ | 90.7$\pm 1.6$ | 91.5$\pm 2.2$ | 92.0$\pm 0.0$ | 89.1$\pm 1.9$ | 92.2$\mathbf{\pm 0.1}$
Chameleon | F1-macro | 38.5 $\pm 2.1$ | 45.0 $\pm 1.0$ | 57.5$\pm 1.0$ | 35.7$\pm 1.8$ | 46.3$\pm 2.4$ | 58.0$\pm 0.1$ | 32.6$\pm 0.6$ | 61.0$\mathbf{\pm 0.3}$
F1-micro | 41.8 $\pm 1.2$ | 44.4 $\pm 1.9$ | 57.1$\pm 1.4$ | 37.7$\pm 1.5$ | 45.7$\pm 3.4$ | 58.2$\pm 0.1$ | 33.2$\pm 0.5$ | 60.8$\mathbf{\pm 0.7}$
AUC | 69.8$\pm 0.5$ | 75.5 $\pm 1.0$ | 85.0$\pm 0.6$ | 69.0 $\pm 1.3$ | 74.8 $\pm 2.8$ | 82.4$\pm 0.0$ | 63.7$\pm 0.8$ | 85.2$\mathbf{\pm 0.3}$
Squirrel | F1-macro | 25.2 $\pm 1.2$ | 26.5 $\pm 1.3$ | 41.1$\pm 1.1$ | 24.7$\pm 2.0$ | 30.5$\pm 3.7$ | 42.1$\pm 0.4$ | 24.7$\pm 1.2$ | 47.5$\mathbf{\pm 0.9}$
F1-micro | 25.8 $\pm 0.8$ | 27.3. $\pm 0.7$ | 43.2$\pm 1.0$ | 28.6 $\pm 1.0$ | 34.6 $\pm 2.5$ | 43.1$\pm 0.3$ | 25.4$\pm 1.0$ | 48.4$\mathbf{\pm 0.8}$
AUC | 57.5 $\pm 0.5$ | 58.2$\pm 1.1$ | 78.9$\pm 0.3$ | 60.2$\pm 1.0$ | 65.6$\pm 1.4$ | 74.3$\pm 0.0$ | 55.0$\pm 2.4$ | 79.4$\mathbf{\pm 0.3}$
Actor | F1-macro | 21.5$\pm 0.4$ | 19.7 $\pm 0.8$ | 30.3$\pm 4.7$ | 28.0$\pm 1.1$ | 30.4$\pm 1.1$ | 31.6$\mathbf{\pm 0.0}$ | 20.0$\pm 0.6$ | 25.7$\pm 0.44$
F1-micro | 29.2 $\pm 0.6$ | 27.1 $\pm 0.5$ | 33.2$\pm 0.6$ | 32.5$\pm 1.0$ | 33.7$\pm 2.3$ | 34.1$\pm 0.0$ | 25.5$\pm 0.5$ | 34.7$\mathbf{\pm 0.4}$
AUC | 58.0 $\pm 0.5$ | 55.8$\pm 0.4$ | 64.8$\pm 0.1$ | 56.2$\pm 2.0$ | 60.8.2$\pm 0.6$ | 64.7$\pm 0.0$ | 53.8$\pm 0.1$ | 68.1$\mathbf{\pm 0.2}$
ogbn-arxiv | Accuracy | 71.5$\pm 0.3$ | 71.6$\pm 0.5$ | 71.2 $\pm 0.3$ | 71.7 $\pm 0.1$ | 71.9 $\pm 0.6$ | 70.9$\pm 0.2$ | 70.0$\pm 0.1$ | 72.2$\mathbf{\pm 0.1}$

Baselines. The proposed HiD-Net is compared with several representative GNNs,
including three traditional GNNs: GCN (Kipf and Welling 2016a), GAT
(Veličković et al. 2017), APPNP (Klicpera, Bojchevski, and Günnemann 2018),
and four graph diffusion networks: GRAND (Chamberlain et al. 2021), GRAND++
(Thorpe et al. 2022), ADC (Zhao et al. 2021), DGC (Wang et al. 2021). They are
implemented based on their open repositories, where the code can be found in
Appendix.

Experimental setting. We perform a hyperparameter search for HiD-Net on all
datasets and the details of hyperparameter can be seen in Appendix. For other
baseline models: GCN, GAT, APPNP, GRAND, GRAND++, DGC, and ADC, we follow the
parameters suggested by (Kipf and Welling 2016a; Veličković et al. 2017;
Klicpera, Bojchevski, and Günnemann 2018; Chamberlain et al. 2021; Thorpe et
al. 2022; Wang et al. 2021; Zhao et al. 2021) on Cora, Citeseer, and Pubmed,
and carefully fine-tune them to get optimal performance on Chameleon,
Squirrel, and Actor. For all methods, we randomly run 5 times and report the
mean and variance. More detailed experimental settings can be seen in
Appendix.

Results. Table 2 summarizes the test results. Please note that OGB prepares
standardized evaluators for testing results and it only provides accuracy
metric for ogbn-arxiv. As can be seen, HiD-Net outperforms other baselines on
seven datasets. Moreover, in comparison with the graph diffusion networks, our
HiD-Net is generally better than them with a large margin on heterophily
graphs, which indicates that our designed graph diffusion process is more
practical for different types of graphs.

(a) Cora

(b) Citeseer

(c) Squirrel

Figure 3: Results of different models under random edge addition.

(a) Cora

(b) Citeseer

(c) Squirrel

Figure 4: Results of different models under random edge deletion.

(a) Cora

(b) Citeseer

(c) Squirrel

Figure 5: Results of different models under random feature perturbation.

(a) Cora

(b) Chameleon

Figure 6: Analysis of parameter $\gamma$.

### Robustness Analysis

Utilizing the information from 2-hop neighbors, our model is more robust in
abnormal situations. We comprehensively evaluate the robustness of our model
on three datasets (Cora, Citeseer and Squirrel) in terms of attacks on edges
and features, respectively.

Attacks on edges. To attack edges, we adopt random edge deletions or additions
following (Chen, Wu, and Zaki 2020; Franceschi et al. 2019). For edge
deletions and additions, we randomly remove or add 5%, 10%, 15%, 20%, 25%,
30%, 35%, 40% of the original edges, which retains the connectivity of the
attacked graph. Then we perform node classification task. All the experiments
are conducted 5 times and we report the average accuracy. The results are
plotted in Figure 3 and Figure 4. From the figures, we can see that as the
addition or deletion rate rises, the performances on three datasets of all the
models degenerate, and HiD-Net consistently outperforms other baselines.

Attacks on features. To attack features, we inject random perturbation into
the node features as in (Wu et al. 2020a). Firstly, we sample a noise matrix
$\mathbf{M}\in\mathbb{R}^{n\times q}$, where each entry in $\mathbf{M}$ is
sampled from the normal distribution $N(0,1)$. Then, we calculate reference
amplitude $r$, which is the mean of the maximal value of each node’s feature.
We add Gaussian noise $\mu\cdot r\cdot\mathbf{M}$ to the original feature
matrix, and get the attacked feature matrix, where $\mu$ is the noise ratio.
The results are reported in Figure 5. Again, HiD-Net consistently outperforms
all other baselines under different perturbation rates by a margin for three
datasets.

### Non-over-smoothing with Increasing Steps

(a) Chameleon

(b) Actor

Figure 7: Non-over-smoothing with increasing steps.

To demonstrate that our model solves the oversmoothing problem compared with
other graph diffusion networks, we test different graph diffusion models with
increasing propagation step $k$ from 2 to 20. Baselines include DGC, ADC and
GRAND. The results are plotted in Figure 7. We can see that with the increase
of $k$, HiD-Net consistently performs better than other baselines.

### Parameter Study

In this section, we investigate the sensitivity of parameters on all datasets.

(a) $\alpha$

(b) $\beta$

Figure 8: Analysis of parameter $\alpha$ and $\beta$.

Analysis of $\alpha$. We test the effect of $\alpha$ in (16), and vary it from
0 to 1. From Figure 8(a) we can see that with the increase of $\alpha$, the
performances of citation graphs rise first and then start to drop slowly, the
performances of Chameleon and Squirrel have not changed too much, and the
performance of Actor first rises and then remains unchanged. As citation
graphs are more homophily, we need to focus less on the node itself, implying
a small $\alpha$, while on heterophily graphs, we need to focus more on the
node itself.

Analysis of $\beta$. In order to check the impact of the diffusion term, we
study the performance of HiD-Net with $\beta$ varying from 0 to 1. The results
are shown in Figure 8(b). We can see that as the value of $\beta$ increases,
the accuracies generally increase, while the accuracy on Actor remains
relatively stable, implying that the features diffused from 1-hop and 2-hop
neighbors are very informative.

Analysis of $\gamma$. Finally we test the effect of $\gamma$ in (16) and vary
it from 0 to 0.6. With the increase of $\gamma$, the accuracies on different
datasets do not change much, so we just separately plot each dataset for a
clearer illustration here. As can be seen in Figure 6, with the increase of
$\gamma$, the performance on Cora and Chameleon rises first and then drops,
and different graphs have different best choices of $\gamma$. The results on
other datasets are shown in Appendix.

## Conclusion

In this paper, we propose a generalized diffusion graph framework, which
establishes the relation between diffusion equation with different GNNs. Our
framework reveals that current graph diffusion networks mainly consider the
first-order diffusion equation, then based on our finding of the monophily
property of labels, we derive a novel high-order diffusion graph network (HiD-
Net). HiD-Net is more robust and general on both homophily and heterophily
graphs. Extensive experimental results verify the effectiveness of HiD-Net.
One potential issue is that our model utilizes a constant diffusivity
coefficient, and a future direction is to explore a learnable diffusivity
coefficient depending on time and space. Our work formally points out the
relation between diffusion equation with a wide variety of GNNs. Considering
that previous GNNs are designed mainly based on spatial or spectral
strategies, this new framework may open a new path to understanding and
deriving novel GNNs. We believe that more insights from the research community
on the diffusion process will hold great potential for the GNN community in
the future.

## Acknowledgments

This work is supported in part by the National Natural Science Foundation of
China (No. U20B2045, 62192784, U22B2038, 62002029, 62172052, 62322203,
62172052).

## References

  * Abu-El-Haija et al. (2019) Abu-El-Haija, S.; Perozzi, B.; Kapoor, A.; Alipourfard, N.; Lerman, K.; Harutyunyan, H.; Ver Steeg, G.; and Galstyan, A. 2019.  Mixhop: Higher-order graph convolutional architectures via sparsified neighborhood mixing.  In _international conference on machine learning_ , 21–29. PMLR.
  * Altenburger and Ugander (2018) Altenburger, K. M.; and Ugander, J. 2018.  Monophily in social networks introduces similarity among friends-of-friends.  _Nature human behaviour_ , 2(4): 284–290.
  * Bruna et al. (2013) Bruna, J.; Zaremba, W.; Szlam, A.; and LeCun, Y. 2013.  Spectral networks and locally connected networks on graphs.  _arXiv preprint arXiv:1312.6203_.
  * Chamberlain et al. (2021) Chamberlain, B.; Rowbottom, J.; Gorinova, M. I.; Bronstein, M.; Webb, S.; and Rossi, E. 2021.  Grand: Graph neural diffusion.  In _International Conference on Machine Learning_ , 1407–1418. PMLR.
  * Chen, Wu, and Zaki (2020) Chen, Y.; Wu, L.; and Zaki, M. 2020.  Iterative deep graph learning for graph neural networks: Better and robust node embeddings.  _Advances in Neural Information Processing Systems_ , 33: 19314–19326.
  * Chung and Graham (1997) Chung, F. R.; and Graham, F. C. 1997.  _Spectral graph theory_.  92\. American Mathematical Soc.
  * Defferrard, Bresson, and Vandergheynst (2016) Defferrard, M.; Bresson, X.; and Vandergheynst, P. 2016.  Convolutional neural networks on graphs with fast localized spectral filtering.  _Advances in neural information processing systems_ , 29.
  * Franceschi et al. (2019) Franceschi, L.; Niepert, M.; Pontil, M.; and He, X. 2019.  Learning discrete structures for graph neural networks.  In _International conference on machine learning_ , 1972–1982. PMLR.
  * Gao and Ji (2019) Gao, H.; and Ji, S. 2019.  Graph u-nets.  In _international conference on machine learning_ , 2083–2092. PMLR.
  * Hamilton, Ying, and Leskovec (2017) Hamilton, W.; Ying, Z.; and Leskovec, J. 2017.  Inductive representation learning on large graphs.  _Advances in neural information processing systems_ , 30.
  * Hu et al. (2020) Hu, W.; Fey, M.; Zitnik, M.; Dong, Y.; Ren, H.; Liu, B.; Catasta, M.; and Leskovec, J. 2020.  Open graph benchmark: Datasets for machine learning on graphs.  _Advances in neural information processing systems_ , 33: 22118–22133.
  * Kipf and Welling (2016a) Kipf, T. N.; and Welling, M. 2016a.  Semi-supervised classification with graph convolutional networks.  _arXiv preprint arXiv:1609.02907_.
  * Kipf and Welling (2016b) Kipf, T. N.; and Welling, M. 2016b.  Variational graph auto-encoders.  _arXiv preprint arXiv:1611.07308_.
  * Klicpera, Bojchevski, and Günnemann (2018) Klicpera, J.; Bojchevski, A.; and Günnemann, S. 2018.  Predict then propagate: Graph neural networks meet personalized pagerank.  _arXiv preprint arXiv:1810.05997_.
  * Liu, Gao, and Ji (2020) Liu, M.; Gao, H.; and Ji, S. 2020.  Towards deeper graph neural networks.  In _Proceedings of the 26th ACM SIGKDD international conference on knowledge discovery & data mining_, 338–348.
  * Liu et al. (2021) Liu, X.; Ding, J.; Jin, W.; Xu, H.; Ma, Y.; Liu, Z.; and Tang, J. 2021.  Graph Neural Networks with Adaptive Residual.  _Advances in Neural Information Processing Systems_ , 34.
  * Ma et al. (2021) Ma, Y.; Liu, X.; Zhao, T.; Liu, Y.; Tang, J.; and Shah, N. 2021.  A unified view on graph neural networks as graph signal denoising.  In _Proceedings of the 30th ACM International Conference on Information & Knowledge Management_, 1202–1211.
  * Pei et al. (2020) Pei, H.; Wei, B.; Chang, K. C.-C.; Lei, Y.; and Yang, B. 2020.  Geom-gcn: Geometric graph convolutional networks.  _arXiv preprint arXiv:2002.05287_.
  * Song et al. (2022) Song, Y.; Kang, Q.; Wang, S.; Kai, Z.; and Tay, W. P. 2022.  On the robustness of graph neural diffusion to topology perturbations.  _arXiv preprint arXiv:2209.07754_.
  * Thorpe et al. (2022) Thorpe, M.; Nguyen, T. M.; Xia, H.; Strohmer, T.; Bertozzi, A.; Osher, S.; and Wang, B. 2022.  GRAND++: Graph neural diffusion with a source term.  In _International Conference on Learning Representations_.
  * Veličković et al. (2017) Veličković, P.; Cucurull, G.; Casanova, A.; Romero, A.; Lio, P.; and Bengio, Y. 2017.  Graph attention networks.  _arXiv preprint arXiv:1710.10903_.
  * Wang et al. (2021) Wang, Y.; Wang, Y.; Yang, J.; and Lin, Z. 2021.  Dissecting the diffusion process in linear graph convolutional networks.  _Advances in Neural Information Processing Systems_ , 34.
  * Weickert (1998) Weickert, J. 1998.  _Anisotropic diffusion in image processing_ , volume 1.  Teubner Stuttgart.
  * Wu et al. (2019) Wu, F.; Souza, A.; Zhang, T.; Fifty, C.; Yu, T.; and Weinberger, K. 2019.  Simplifying graph convolutional networks.  In _International conference on machine learning_ , 6861–6871. PMLR.
  * Wu, He, and Xu (2019) Wu, J.; He, J.; and Xu, J. 2019.  Net: Degree-specific graph neural networks for node and graph classification.  In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, 406–415.
  * Wu et al. (2020a) Wu, T.; Ren, H.; Li, P.; and Leskovec, J. 2020a.  Graph information bottleneck.  _Advances in Neural Information Processing Systems_ , 33: 20437–20448.
  * Wu et al. (2020b) Wu, Z.; Pan, S.; Chen, F.; Long, G.; Zhang, C.; and Philip, S. Y. 2020b.  A comprehensive survey on graph neural networks.  _IEEE transactions on neural networks and learning systems_ , 32(1): 4–24.
  * You, Ying, and Leskovec (2019) You, J.; Ying, R.; and Leskovec, J. 2019.  Position-aware graph neural networks.  In _International Conference on Machine Learning_ , 7134–7143. PMLR.
  * Yu et al. (2023) Yu, X.; Liu, Z.; Fang, Y.; and Zhang, X. 2023.  Learning to Count Isomorphisms with Graph Neural Networks.  In Williams, B.; Chen, Y.; and Neville, J., eds., _Thirty-Seventh AAAI Conference on Artificial Intelligence, AAAI 2023, Thirty-Fifth Conference on Innovative Applications of Artificial Intelligence, IAAI 2023, Thirteenth Symposium on Educational Advances in Artificial Intelligence, EAAI 2023, Washington, DC, USA, February 7-14, 2023_ , 4845–4853. AAAI Press.
  * Zang and Wang (2020) Zang, C.; and Wang, F. 2020.  Differential Deep Learning on Graphs and its Applications.
  * Zhang et al. (2018) Zhang, M.; Cui, Z.; Neumann, M.; and Chen, Y. 2018.  An end-to-end deep learning architecture for graph classification.  In _Thirty-second AAAI conference on artificial intelligence_.
  * Zhao et al. (2021) Zhao, J.; Dong, Y.; Ding, M.; Kharlamov, E.; and Tang, J. 2021.  Adaptive Diffusion in Graph Neural Networks.  _Advances in Neural Information Processing Systems_ , 34.
  * Zhou et al. (2020) Zhou, J.; Cui, G.; Hu, S.; Zhang, Z.; Yang, C.; Liu, Z.; Wang, L.; Li, C.; and Sun, M. 2020.  Graph neural networks: A review of methods and applications.  _AI Open_ , 1: 57–81.
  * Zhu et al. (2021) Zhu, M.; Wang, X.; Shi, C.; Ji, H.; and Cui, P. 2021.  Interpreting and unifying graph neural networks with an optimization framework.  In _Proceedings of the Web Conference 2021_ , 1215–1226.


# Soliton blockade in bidirectional microresonators

Zhiwei Fan  Department of Physics, University of Bath, Bath BA2 7AY, UK
Dmitry V. Skryabin  Department of Physics, University of Bath, Bath BA2 7AY,
UK  Corresponding author<EMAIL_ADDRESS>
###### Abstract

We report a method to control - disrupt and restore, a regime of the
unidirectional soliton generation in a bidirectionally pumped ring
microresonator. This control, i.e., the soliton blockade, is achieved by
tuning pump frequency of the counterrotating field. The blockade effect is
correlated with the emergence of a dark-bright nonlinear resonance of the cw-
states.

††journal: ol

Refractive index control with light is an important tool of achieving the
functionality of photonic circuits and optical devices, see, e.g., [2]. One of
its applications is inducing shifts of the resonance frequencies and thereby
controlling reflectivity and transmission of resonators. Ultimate efficiency
of this approach is achieved when one photon interacting with one atom is
enough to alter the system transmission [4, 3, 5]. A famous example of this is
photon blockade [4, 3]. In this work, we demonstrate that tuning frequency of
one of the counterrotating fields in a high-Q ring microresonator can disrupt
and restore soliton transmission in the other. We term this phenomenon -
soliton blockade.

Frequency conversion and soliton generation in high-Q microresonators impact
many areas of modern photonics from precision spectroscopy to optical
information processing [6]. Bi-directionally pumped resonator have recently
became an active sub-area in this field [7, 8, 9]. Rich symmetry breaking
phenomenology [10, 11, 12, 13] and applications in velocity measurements [14,
15, 16] have been some of the underpinning driving forces. A key property of
microresonators is their high finesse which can provide three to six orders of
magnitude boost to the circulating powers relative to the input one. This
creates a variety of opportunities for efficient nonlinear control of optical
signals, and, in particular, the soliton blockade effect is utilising high
sensitivity of the nonlinear response of one of the intracavity fields towards
changes in the driving frequency of the counterrotating one.

One of the bottlenecks in studies of the bidirectional microresonators has
been the absence of a sufficiently rigorous mathematical model to describe
their operation in the multimode regime. A variety of qualitative models [8,
9] have been used to treat problems of the nonlinear crosscoupling between the
modes, backscattering, and also of the opposing group velocities of the two
waves in different ways, see [17] for a comparative overview. These issues
have been reconciled by demonstrating an equivalence of the coupled-mode model
derived from the ab-initio Maxwell system to a certain form of the envelope
equations [17, 18, 19]. Below, we are applying this methodology to study the
soliton blockade and accompanying effects.

Following Ref. [17], we define total electric field, ${\cal E}$, in a
microresonator as a superposition of the clockwise, ’+’, and counterclockwise,
’$-$’, rotating modes. If $\psi^{\pm}_{\mu}(t)$ are the complex mode
amplitudes, $\mu=-N/2+1,\cdots,0,\cdots N/2$ is the relative mode number, $t$
is time, and $\theta$ is angle varying along the resonator circumference, then
${\cal
E}=\big{[}e^{iM\theta-i\omega_{+}t}\sum_{\mu}\psi_{\mu}^{+}e^{i\mu(\theta-
D_{1}t)}+c.c.\big{]}+\big{[}e^{iM\theta+i\omega_{+}t}\sum_{\mu}\psi_{\mu}^{-*}e^{i\mu(\theta+D_{1}t)}+c.c.\big{]}$.
$D_{1}/2\pi$ is the repetition rate, i.e., free spectral range. $M$ is the
absolute mode number corresponding to the resonance frequency
$\omega_{\mu=0}=\omega_{0}$. The first square bracket in the ${\cal
E}$-equation corresponds to the clockwise rotating field, and the second to
the counterclockwise one. $\omega_{+}$ is the laser frequency pumping the ’+’
wave, which is detuned from $\omega_{0}$ by $\delta=\omega_{0}-\omega_{+}$,
and by $\varepsilon=\omega_{+}-\omega_{-}$ from the pump frequency of the
’$-$’ wave, $\omega_{-}$.

If $\psi^{\pm}(t,\theta)=\sum_{\mu}\psi^{\pm}_{\mu}(t)e^{i\mu\theta}$, then
coupled mode equations for $\psi^{\pm}_{\mu}$ are equivalent to the following
equations for the $\psi^{\pm}$ envelopes [17],

$\displaystyle i\partial_{t}\psi^{+}$
$\displaystyle=-\tfrac{1}{2}D_{2}\partial_{\theta}^{2}\psi^{+}$ (1a)
$\displaystyle+\delta\psi^{+}-\gamma|\psi^{+}|^{2}\psi^{+}-2\textsl{g}^{-}\psi^{+}-i\tfrac{1}{2}\kappa(\psi^{+}-{\cal
H}),$ $\displaystyle i\partial_{t}\psi^{-}$
$\displaystyle=-\tfrac{1}{2}D_{2}\partial_{\theta}^{2}\psi^{-}$ (1b)
$\displaystyle+\delta\psi^{-}-\gamma|\psi^{-}|^{2}\psi^{-}-2\textsl{g}^{+}\psi^{-}-i\tfrac{1}{2}\kappa(\psi^{-}-{\cal
H}e^{i\varepsilon t}),$ $\displaystyle\textsl{g}^{\pm}$
$\displaystyle=\gamma\int_{0}^{2\pi}|\psi^{\pm}|^{2}\frac{d\theta}{2\pi}=\gamma\sum_{\mu}|\psi_{\mu}^{\pm}|^{2}.$
(1c)

In Eqs. (1), the offset between the pump frequencies, $\varepsilon\neq 0$, is
the only parameter that breaks the symmetry between the two directions and it
is used below as the main control parameter.

Figure 1: (a) Nonlinear resonances of the single mode states, i.e. cw-cw
states, $\textsl{g}^{\pm}=\gamma|\widehat{\psi}^{\pm}|^{2}$, vs $\delta$ for a
set of fixed $\varepsilon$. Same colours show matching branches on the
$\textsl{g}^{+}$ and $\textsl{g}^{-}$ plots. (b) $\textsl{g}^{\pm}$ vs
$\varepsilon$ for $\delta$ fixed. Black lines show the dark-bright resonance
around $\varepsilon=0$. Colours and line styles in (b) are independent from
(a). Here, they show how branches of the $\textsl{g}^{+}$ (left) and
$\textsl{g}^{-}$ (right) cw-cw states match at their turning points. (c)
Soliton blockade effect in the ’$+$’ field, when $\varepsilon$ is
adiabatically scanned from being large negative, across zero, and towards
large positive values. Colour density shows $|\psi^{\pm}|^{2}(\theta)$ along
the resonator circumference, $\theta\in[-\pi,\pi)$. Long and sparse straight
lines correspond to the soliton-cw states. (d) Mode spectrum corresponding to
the real space waveforms in (c). ${\cal H}^{2}=16$W, $\kappa/2\pi=1.5$MHz,
$D_{1}/2\pi=15$GHz, $\gamma/2\pi=0.4$MHz/W in (a)-(d), $\delta/2\pi=4.5$MHz in
(b)-(d), and $D_{2}/2\pi=10$kHz in (c), (d). Figure 2:  (a,b) Soliton
components of the soliton-cw (a) and of the cw-soliton (b) states as computed
from the time-independent models, Eqs. (3), (4), and corresponding to the
solitons found in the dynamical simulations, see Fig. 1(c). (c,d) is the
soliton-soliton state. Panel (a) shows the soliton blockade. The blue color in
(b) marks the $\varepsilon$ interval where the ’$-$’ soliton is sustained by
the chaotic state in the ’$+$’ wave, that provides relatively large XPM
induced resonance shift, $-2g^{+}$, compensating for $\varepsilon>0$.

$-2\textsl{g}^{\pm}$ are the negative frequency shifts of the cavity resonance
$\omega_{0}$ induced by the nonlinear cross-coupling, i.e., by the XPM (cross-
phase modulation) effects. These shifts are different for ’$+$’ and ’$-$’
waves. In the high repetition rate devices, the phase dependent four-wave
mixing terms inside the XPM part of the Kerr response,
$|\psi^{\mp}|^{2}\psi^{\pm}$, oscillate with the fast rate $\sim 2D_{1}$ [17,
18, 19]. If finesse ${\cal F}=D_{1}/\kappa\gg 1$, then the fast oscillations
can be disregarded and this is why the XPM nonlinearity in Eq. (1c), i.e.,
$2\textsl{g}^{\mp}\psi^{\pm}$, contains only powers of the interacting modes
and not their phase differences. This approximation was recently formally
justified and used for Fabry-Perot resonators [19], and soon after applied for
the ring ones [17, 18].

$D_{2}$ is the 2nd order dispersion coefficient, so that the resonator
spectrum is $\omega_{\mu}=\omega_{0}+D_{1}\mu+D_{2}\mu^{2}/2!$. $\gamma$ is
the nonlinear parameter [17]. $\kappa$ is the resonance linewidth. The
intracavity pump amplitude, ${\cal H}$, is linked to the laser power, $\cal
W$, as ${\cal H}^{2}$ = $\frac{\eta}{\pi}{{\cal F}}\cal W$, where $\eta<1$ is
the coupling coefficient.

To study a resonator operating in the single-mode regime in both directions we
set $\psi^{+}=\widehat{\psi}^{+}_{cc}$,
$\psi^{-}=\widehat{\psi}^{-}_{cc}e^{i\varepsilon t}$,
$\textsl{g}^{\pm}=\gamma|\widehat{\psi}^{\pm}_{cc}|^{2}$,
$\partial_{\theta}\widehat{\psi}^{\pm}_{cc}=\partial_{t}\widehat{\psi}^{\pm}_{cc}=0$.
These are cw-cw states describing counterrotating fields with their carrier
frequencies locked to the respective pump frequencies. In this regime Eqs. (1)
can be reduced to a pair of real equations,

$\displaystyle\gamma{\cal H}^{2}$
$\displaystyle=\textsl{g}^{+}+4\textsl{g}^{+}\big{(}\delta-\textsl{g}^{+}-2\textsl{g}^{-}\big{)}^{2}/\kappa^{2},$
(2a) $\displaystyle\gamma{\cal H}^{2}$
$\displaystyle=\textsl{g}^{-}+4\textsl{g}^{-}\big{(}\delta+\varepsilon-\textsl{g}^{-}-2\textsl{g}^{+}\big{)}^{2}/\kappa^{2}.$
(2b)

Large detunings between the two pump frequencies, $|\varepsilon|/\kappa\gg 1$,
effectively decouple nonlinear resonances of the two fields. The latter
implies that the $\textsl{g}^{+}$ vs $\delta$ plot shows the nonlinearity
tilted resonance originating from $\delta=0$, and $\textsl{g}^{-}$ has a
similar resonance at $\delta=-\varepsilon$. Taking several representative
values of $\varepsilon/2\pi$ and other parameters as in Fig. 1(a) shows that
varying $\varepsilon$ from the relatively large negative to large positive
values drags the $\textsl{g}^{-}$ resonance across the effectively ’immobile’
$\textsl{g}^{+}$ one, see Fig. 1(a). Two resonances overlap and interact
strongly for $|\varepsilon|/\kappa\lesssim 1$.

For $\varepsilon=0$, Eq. (1a) and Eq. (1b) become symmetric relative to
$\psi^{+}\to\psi^{-}$. Hence, there exists a symmetric cw-cw state,
$\widehat{\psi}^{+}_{cc}=\widehat{\psi}^{-}_{cc}$, see black lines in the
$\varepsilon=0$ panels of Fig. 1(a). As $\delta$ varies, this state undergoes
a symmetry breaking bifurcation into a pair of asymmetric solutions, see blue
and green lines in the same panels. Linear stability of
$\widehat{\psi}^{\pm}_{cc}$ relative to noise in the $\mu=0$ mode itself and
in the $\mu\neq 0$ modes has been evaluated using a substitution
$\psi^{\pm}=\big{(}\widehat{\psi}^{\pm}_{cc}+\widetilde{\psi}^{\pm}_{\mu}e^{i\mu\theta+\lambda_{\mu}t}+\widetilde{\psi}^{\pm*}_{-\mu}e^{-i\mu\theta+\lambda_{\mu}^{*}t}\big{)}e^{i\varepsilon^{\pm}t}$,
$\varepsilon^{+}=0$, $\varepsilon^{-}=\varepsilon$, and linearisation of Eqs.
(1), $|\widehat{\psi}^{\pm}_{\mu}|/|\widetilde{\psi}^{\pm}_{cc}|\ll 1$.

What is critical for the blockade effect of the ’$+$’ solitons induced by the
’$-$’ wave, is how $\textsl{g}^{\pm}$ in the cw-cw regime vary with
$\varepsilon$ for fixed $\delta$, i.e. when the frequency of the ’$+$’ pump is
fixed and frequency of the ’$-$’ pump is tuned. The respective
$\textsl{g}^{\pm}$ vs $\varepsilon$ plots change qualitatively depending on
how $\delta$ is chosen. It is critical for what follows to take $\delta$
within the bistability interval of the $\textsl{g}^{+}$ resonance, e.g.,
$\delta/2\pi=4.5$MHz as is in Figs. 1(b)-(d). This choice of $\delta$ ensures
that the ’$+$’ wave can be excited into a strongly nonlinear state
independently if $\varepsilon/\kappa$ is the large positive or large negative.
While the ’$-$’ wave can respond to the pump in a bistable manner only for
$|\varepsilon|/\kappa\sim 1$ or less. From this one can anticipate a resonance
feature in how the ’$-$’ wave impacts transmission of the ’$+$’ wave at
$\varepsilon$ around zero.

Three cw-cw states coexisting for $\varepsilon/2\pi\lesssim-7$MHz are well
approximated by the $\textsl{g}^{+}$ bistability loop with
$\textsl{g}^{-}\approx 0$. Since dispersion is anomalous, $D_{2}>0$, the low
amplitude state is stable and the high amplitude one is unstable. $\mu\neq 0$
instabilities of the latter generate solitons in the ’$+$’ wave, while the
’$-$’ wave remains in the cw-state, see Fig. 1(c). The low amplitude state
serves as the soliton background. As $\varepsilon$ is moved towards zero, the
$\textsl{g}^{-}$ power starts building up, and the first effect it brings is
that the soliton background becomes unstable relative to the $\mu=0$ noise,
see green points in Fig. 1(b). This instability initiates background
oscillations which do not make much impact on the soliton cores. Frequency of
the oscillations is the order of $10$ to $10^{2}$kHz, and they imply that the
counterrotating waves are becoming unlocked from the pump frequencies and tend
to lock between themselves. Soon, the background becomes unstable relative to
the broadband $\mu\neq 0$ perturbations, see blue points in Fig. 1(b). This
brings a more detrimental impact on solitons, that is becoming destroyed by
the chaotic multi-mode oscillations, see Figs. 1(c).

Around $\varepsilon/2\pi\simeq-1$MHz, the low amplitude cw (full blue line)
merges with the middle branch (dotted blue line) and both disappear, see Fig.
1(b). This implies that the ’$+$’ solitons become strictly forbidden. Instead,
one can see the formation of the pronounced quasi-stationary roll pattern in
the ’$-$’ component, see Fig. 1(c). At the same time, $\textsl{g}^{+}$
corresponding to the high power branch of the cw-cw state starts going down
sharply, while the respective $\textsl{g}^{-}$ is rising, which creates the
dark-bright resonance, see black lines in Fig. 1(b). Large values of
$\textsl{g}^{-}$ make the net effective detuning of the $\textsl{g}^{+}$
field, i.e., $\delta-2\textsl{g}^{-}$, negative and hence the ’$+$’ only
solitons impossible. Contrary $\textsl{g}^{+}$ is now relatively small, and as
$\varepsilon$ is increased little further the solitons emerge inside the ’$-$’
wave, see $3\lesssim\varepsilon/2\pi\lesssim 7$MHz in Fig. 1(c). For
$\varepsilon/2\pi\gtrsim 7$MHz they disappear and the $\textsl{g}^{-}$
resonance again has no impact on the ’$+$’ one, and therefore the soliton
states in the ’$+$’ field are restored.

The above scenario, see the left panel of Fig. 1(c), constitutes the blockade
effect for the ’$+$’ solitons created by tuning the resonance of the ’$-$’
field across the ’$+$’ one, see Figs. 1(a),(c). An approximate analytical
expression for the interval of $\varepsilon=\omega_{+}-\omega_{-}$ where the
blockade is taking place can be worked out using following considerations. For
$|\varepsilon|/\kappa\gg 1$, $\textsl{g}^{-}\approx 0$ and $\gamma{\cal
H}^{2}\approx\textsl{g}^{+}+4\textsl{g}^{+}\left(\delta-\textsl{g}^{+}\right)^{2}/\kappa^{2}$.
Values of $\textsl{g}^{+}$ and $\delta$ at the tip of the nonlinear resonance
are connected one to the other by
$\textsl{g}^{+}=\tfrac{2}{3}\delta+\tfrac{1}{3}\sqrt{\smash[b]{\delta^{2}-\tfrac{3}{4}\kappa^{2}}}\approx\delta-\kappa^{2}/8\delta$
[19]. Hence $\gamma{\cal H}^{2}\approx\delta-\kappa^{2}/16\delta$, and
$\delta$ at the tip is well approximated by
$\delta_{tip}=\tfrac{1}{2}\gamma{\cal
H}^{2}+\tfrac{1}{2}\sqrt{\smash[b]{\gamma^{2}{\cal
H}^{4}+\tfrac{1}{4}\kappa^{2}}}$. Thus, $|\varepsilon|\lesssim\delta_{tip}$
provides quasi-resonant interaction of the counterrotating waves and can be
used as an estimate for the width of the blockade interval. Using parameters
as in Fig. 1 gives $\delta_{tip}/2\pi\approx 6.4$MHz, and
$-6.4\lesssim\varepsilon/2\pi\lesssim 6.4$MHz is in the remarkably close
agreement with numerical data in Fig. 1(c).

Steady-state soliton solutions of Eqs. (1) can be divided into three broad
categories of the soliton-cw, cw-soliton, and soliton-soliton states. Here,
the order of words indicates if the soliton component belongs to the ’$+$’ or
to the ’$-$’ wave. The soliton-cw states, which blockade effect we report
here, have been observed directly in the numerical scan data shown Fig. 1(c).
The two remaining states are complimentary to the blockade effect. The state
close to the cw-soliton one appears in the scan with its ’+’ component
distorted by the modulations, cf., [20, 21].

The soliton-cw state is sought in the form
$\psi^{+}=\widehat{\psi}^{+}_{sc}(\theta),\psi^{-}=\widehat{\psi}^{-}_{sc}\leavevmode\nobreak\
e^{i\varepsilon t}$, and satisfies

$\displaystyle\tfrac{1}{2}D_{2}\partial_{\theta}^{2}\widehat{\psi}^{+}_{sc}$
$\displaystyle=\big{(}\delta-\gamma|\widehat{\psi}^{+}_{sc}|^{2}-2\textsl{g}^{-}\big{)}\widehat{\psi}^{+}_{sc}-\tfrac{i\kappa}{2}\big{(}\widehat{\psi}^{+}_{sc}-{\cal
H}\big{)},$ (3a) $\displaystyle\gamma{\cal H}^{2}$
$\displaystyle=\textsl{g}^{-}+4\textsl{g}^{-}\big{(}\delta+\varepsilon-\textsl{g}^{-}-2\textsl{g}^{+}\big{)}^{2}/\kappa^{2}.$
(3b)

Similarly, the cw-soliton states are
$\psi^{+}=\widehat{\psi}^{+}_{cs},\psi^{-}=\widehat{\psi}^{-}_{cs}(\theta)\leavevmode\nobreak\
e^{i\varepsilon t}$, and satisfy

$\displaystyle\gamma{\cal H}^{2}$
$\displaystyle=\textsl{g}^{+}+4\textsl{g}^{+}\big{(}\delta-\textsl{g}^{+}-2\textsl{g}^{-}\big{)}^{2}/\kappa^{2},$
(4a)
$\displaystyle\tfrac{1}{2}D_{2}\partial_{\theta}^{2}\widehat{\psi}^{-}_{cs}$
$\displaystyle=\big{(}\delta+\varepsilon-\gamma|\widehat{\psi}^{-}_{cs}|^{2}-2\textsl{g}^{+}\big{)}\widehat{\psi}^{-}_{cs}-\tfrac{i\kappa}{2}\big{(}\widehat{\psi}^{-}_{cs}-{\cal
H}\big{)}.$ (4b)

The soliton part of the steady-state soliton-cw solution plotted vs
$\varepsilon$, see Fig. 2(a), shows the blockade interval. The latter is
narrower, than the one in the dynamical scan, because the scan picks the
dynamic non-stationary regime for $\varepsilon<0$, and drags it adiabatically
to $\varepsilon>0$, see Fig. 1(c).

The cw-soliton state in Fig. 2(b) and the soliton-soliton state in Figs.
2(c),(d) also exist. The latter are found by jointly solving Eqs. (3a), (4b),
and can be anticipated to emerge from the $\varepsilon=0$ point, when this
pair of equations is reduced to the single-component case. A range of
existence of the ’-’ solitons in the dynamical scan is wider than is found
from the time-independent Eqs. (4), see the blue interval in Fig. 2(b). This
is because the scan, see Fig. 1(c), picks the modulated solution in the ’$+$’
component. The modulated wave pattern provides larger nonlinear shift
$2g^{+}$, for the ’$-$’ wave than the exact cw component of the cw-soliton
state. Larger negative $-2g^{+}$ compensate for the growing $\varepsilon>0$
and keep the ’$-$’ soliton afloat for longer, thereby, extending the blockade
range. All three soliton families in Fig. 2 have been found to have accessible
regimes of stable dynamical evolution, also confirmed by the linear stability
analysis. However, the oscillatory instabilities have also been seen for
$|\varepsilon|$ close to zero, which could be a subject of a separate
investigation. Interplays of the soliton-blockade with the power imbalance,
backscattering, higher-order dispersions and Raman scattering also represent
interesting problems. E.g., the Raman effect could be used as a natural
frequency shifter to the soliton channel [8, 20].

In summary, we have reported how tuning frequency of one of the two
counterrotating fields causes the dark-bright nonlinear resonance in the
single-mode microresonator operation, which translates to the soliton blockade
effect in the multi-mode regime. The soliton blockade constitutes a disruption
of the soliton transmission in the signal field by changing the frequency of
the control field.

Funding: EU Horizon 2020 Framework Programme (812818, MICROCOMB).

Disclosures: The authors declare no conflicts of interest.

## References

  * [1]
  * [2] V. R. Almeida, C. A. Barrios, R. R. Panepucci, and M. Lipson, Nature 431, 1081–1084 (2004).
  * [3] D. E. Chang, V. Vuletić, and M.D. Lukin, Nat. Phot. 8, 685–694 (2014).
  * [4] K. M. Birnbaum, A. Boca, R. Miller, A. D. Boozer, T. E. Northup, and H.J. Kimble, Nature 436, 87–90 (2005).
  * [5] T. Aoki, B. Dayan, E. Wilcut, W. P. Bowen, A. S. Parkins, T. J. Kippenberg, K.J. Vahala, and H.J. Kimble, Nature 443, 671–674 (2006).
  * [6] T.J. Kippenberg, A.L. Gaeta, M. Lipson, and M. Gorodetsky, Science 361, eaan8083 (2018).
  * [7] C. Joshi, A. Klenner, Y. Okawachi, M. Yu, K. Luke, X. Ji, M. Lipson, and A.L. Gaeta, Opt. Lett. 43, 547 (2018).
  * [8] Q.F. Yang, X. Yi, K.Y. Yang, and K. Vahala, Nat. Photon. 11, 560-564 (2017).
  * [9] W. Weng, R. Bouchand, E. Lucas, and T.J. Kippenberg, Phys. Rev. Lett. 123, 253902 (2019).
  * [10] M.T.M. Woodley, J.M. Silver, L. Hill, F. Copie, L.D. Bino, S. Zhang, G. Oppo, and P. Del’Haye, Phys. Rev. A 98, 053863 (2018).
  * [11] F. Copie, M.T.M. Woodley, L. Del Bino, J. M. Silver, S. Zhang, and P. Del’Haye, Phys. Rev. Lett. 122, 013905 (2019).
  * [12] L. Hill, G. L. Oppo, M.T.M. Woodley, P. Del’Haye, Phys. Rev. A 101, 013823 (2020).
  * [13] Q. Cao, H. Wang, C. Dong, H. Jing, R. Liu, X. Chen, L. Ge, Q. Gong, and Y. Xiao, Phys. Rev. Lett. 118, 033901 (2017).
  * [14] J. Li, M.-G. Suh, and K. Vahala, Optica 4, 346-348 (2017).
  * [15] Y.H. Lai, Y.K. Lu, M.G. Suh, Z.Q. Yuan, and K. Vahala, Nature 576, 65 (2019).
  * [16] W. Liang, V.S. Ilchenko, A.A. Savchenkov, E. Dale, D. Eliyahu, A.B. Matsko, and L. Maleki, Optica 4, 114-117 (2017).
  * [17] D. V. Skryabin, OSA Continuum 3, 1364-1375 (2020).
  * [18] N.M. Kondratiev and V.E. Lobanov, Phys. Rev. A 101, 013816 (2020).
  * [19] D.C. Cole, A. Gatti, S.B. Papp, F. Prati, and L. Lugiato, Phys. Rev. A 98, 013831 (2018).
  * [20] M. Anderson, Y. Wang, F. Leo, S. Coen, M. Erkintalo, and S.G. Murdoch, Phys. Rev. X 7, 031031 (2017).
  * [21] A.U. Nielsen, B. Garbin, S. Coen, S.G. Murdoch, and M. Erkintalo, Phys. Rev. Lett. 123, 013902 (2019).