# A Simple and Powerful Framework for
Stable Dynamic Network Embedding

Ed Davis<EMAIL_ADDRESS>
School of Mathematics
University of Bristol
Bristol, United Kingdom Ian Gallagher<EMAIL_ADDRESS>
School of Mathematics
University of Bristol
Bristol, United Kingdom
Daniel John Lawson<EMAIL_ADDRESS>
School of Mathematics
University of Bristol
Bristol, United Kingdom
Patrick Rubin-Delanchy<EMAIL_ADDRESS>
School of Mathematics
University of Bristol
Bristol, United Kingdom

###### Abstract

In this paper, we address the problem of dynamic network embedding, that is,
representing the nodes of a dynamic network as evolving vectors within a low-
dimensional space. While the field of static network embedding is wide and
established, the field of dynamic network embedding is comparatively in its
infancy. We propose that a wide class of established static network embedding
methods can be used to produce interpretable and powerful dynamic network
embeddings when they are applied to the dilated unfolded adjacency matrix. We
provide a theoretical guarantee that, regardless of embedding dimension, these
unfolded methods will produce stable embeddings, meaning that nodes with
identical latent behaviour will be exchangeable, regardless of their position
in time or space.

We additionally define a hypothesis testing framework which can be used to
evaluate the quality of a dynamic network embedding by testing for planted
structure in simulated networks. Using this, we demonstrate that, even in
trivial cases, unstable methods are often either conservative or encode
incorrect structure. In contrast, we demonstrate that our suite of stable
unfolded methods are not only more interpretable but also more powerful in
comparison to their unstable counterparts.

Keywords: dynamic network embedding, graph representation learning, hypothesis
testing, dynamic networks, temporal clustering.

##  1 Introduction

Network embedding is the process of finding a $d$-dimensional vector
representation for each node in a network, upon which downstream tasks such as
regression or clustering can be used to extract information about the network.
This is a thriving field of research with wide applications from natural
language (Mikolov et al., 2013) to neuroscience (Rosenthal et al., 2018).
While there has been a significant body of research behind the embedding of
static networks, we argue that the field of dynamic network embedding is still
in its infancy, despite offering an even wider range of potential
applications. Due to the maturity of the static network embedding literature,
a tantalising goal would be to find a way of adapting established static
network embeddings such that they perform well on dynamic network embedding
tasks.

We consider the problem of embedding discrete-time dynamic networks, i.e.
those that can be represented as a series of symmetric (or undirected)
adjacency matrix “snapshots” over time,
$\mathbf{A}^{(1)},\dots,\mathbf{A}^{(T)}\in\\{0,1\\}^{n\times n}$. A dynamic
network embedding is then the process of finding a representation
$\hat{\mathbf{Y}}_{i}^{(t)}\in\mathbb{R}^{d}$ for each node
$i\in\\{1,\dots,n\\}$ in each snapshot $t\in\\{1,\dots,T\\}$. We will refer to
each
$\hat{\mathbf{Y}}^{(1)},\dots,\hat{\mathbf{Y}}^{(T)}\in\mathbb{R}^{n\times d}$
as embedding time points.

In a recent paper, Gallagher et al. (2021) showed that unfolded adjacency
spectral embedding (UASE) (Jones and Rubin-Delanchy, 2020) provides a
theoretically grounded procedure for obtaining a dynamic network embedding
which is stable in both time and space (these details and others from the
introduction are formally defined below). One implication of the stability
guarantee is that dynamic community structure can be recovered consistently
using spatio-temporal clustering by computing a single spectral embedding on a
column-concatenation of adjacency snapshots,
$\mathbfcal{A}=\left(\mathbf{A}^{(1)},\dots,\mathbf{A}^{(T)}\right)$, referred
to as the unfolded adjacency matrix. This is the only method to our knowledge
that can offer such guarantees. However, the authors leave the question open
of how to achieve such guarantees with methods other than adjacency spectral
embedding. For example, regularised Laplacian spectral embedding can offer
better results in sparse regimes in comparison to unregularised spectral
embedding (Qin and Rohe, 2013) and node2vec (Grover and Leskovec, 2016) is
popular and empirically successful. Generalisation is a complicated problem
however, as Gallagher et al. (2021) proved the stability of UASE using its
central limit theorem, something that cannot be done for a general class of
methods. Up to now, the choice has been to either use UASE to consistently
capture spatial and temporal information while being limited to adjacency
spectral embedding, or to use other embedding methods at the cost of temporal
information. In this paper, we extend these stability guarantees to encompass
a wide class of established network embedding methods, circumventing the need
for a central limit theorem.

The concepts of cross-sectional and longitudinal stability for dynamic
embeddings (which we refer to as spatial and temporal stability), defined by
Gallagher et al. (2021), are hugely important. A stable dynamic embedding is
one where the ordering of input snapshots does not change the probabilistic
outcome of the embedding (up to permutation), and where any two nodes with the
same latent behaviour are exchangeable, even if they are observed at different
times. We build on the work of Gallagher et al. (2021) by extending these
definitions of stability to encompass a wider range of dynamic embedding
methods, including those which are non-deterministic. This extension is
critical for addressing the open question of how to extend the stability
guarantees of UASE to a wider range of methods.

Our main theorem then provides a theoretical guarantee that a wide class of
static network embedding methods will produce a stable dynamic embedding on
the input of the dilated unfolded adjacency matrix. From this, the jump from
static network to dynamic network embedding is straightforward. Further, this
guarantee holds for _any_ embedding dimension $d$, which offers a powerful
advantage over unstable methods whose embedding dimension can greatly affect
the encoding of temporal information. Under our framework, the dimension $d$
is chosen to gain statistical power (see Section 4).

A popular alternative approach to dynamic embedding has been to assume
temporal smoothness, that is, that dynamic networks evolve smoothly and
therefore that subsequent snapshots are similar. Under this assumption,
dynamic embedding methods have been proposed that recurrently align single
network embeddings, for example by Procrustes alignment rotation (Zhou et al.,
2019; Singer et al., 2019) or by initialising the training of node vectors at
time $t$ with the node vectors at $t-1$ (Mahdavi et al., 2018; Zhou et al.,
2018; Goyal et al., 2018; Hou et al., 2020). However, these alignments are
only based on subsequent time points, meaning that the embedding time point
$\hat{\mathbf{Y}}^{(u)}$ is not likely to be related to
$\hat{\mathbf{Y}}^{(t)}$ when $u\ll t$. This means that long-range temporal
structure is lost and stability cannot be guaranteed. Further, the smoothness
assumption may not hold. In contrast, unfolding stands out for its theoretical
guarantees and simplicity.

To formalise a testable notion of stability, we introduce a hypothesis test
that exploits exchangeability over time and/or space and therefore can be used
to quantify the stability of any dynamic embedding. We show that in trivial
settings, dynamic embedding methods without stability often encode incorrect
spatial and temporal relations. Surprisingly, even in settings when snapshots
are independently and identically distributed (the simplest possible dynamic
network), unstable methods produce embeddings which are, at best, conservative
under this hypothesis test. We demonstrate that this does not occur with
stable methods. Due to these findings, we recommend that unstable dynamic
embedding methods should be replaced by stable versions with more
interpretable meaning; something that we show can be done with only a few
lines of code.

As a demonstration of practical power, we consider a dynamic flight network
from 2019 to 2021. We show that dynamic network embeddings based on the
dilated unfolded adjacency matrix can encode the seasonal periodic behaviour
present in the flight network, as well as the two major European waves of the
COVID-19 pandemic as it disrupted air traffic; properties that could not be
encoded using unstable methods.

The structure of the paper is as follows. Section 2 motivates the use of
stable dynamic embedding methods with a trivial example. Section 3 introduces
the problem of dynamic network embedding, introduces definitions of stability
and proves our main theorem. In Section 4, we introduce a hypothesis test and
use it to compare the stability of various dynamic embedding methods on
simulated networks. In Section 5, we compare our stable methods against
existing dynamic network embedding methods on a real-world dynamic flight
network. Finally, Section 6 summarises our findings and recommends the use of
stable dynamic embedding methods.

##  2 Motivation

Let us define a trivial problem. Suppose that we have a network with two
groups of nodes that have similar behaviour, which we call node communities.
Suppose also that as this network evolves over time, one community changes
behaviour, while the other remains static. Then, we expect that the dynamic
embedding of this network should encode both the static and changing behaviour
demonstrated by the respective communities. Despite the simplicity of this
problem, there are almost no dynamic embedding methods that can solve it.

Figure 1 displays the embeddings from spectral- and skip-gram-based methods
with the use of three different temporal mechanisms. The leftmost column of
Figure 1 demonstrates the naïve approach to this problem, where static network
embedding methods have been independently applied to each snapshot. This
independent application places each embedding time point in a different space,
thereby making comparisons over time impossible and the temporal structure
present in the network is lost. This is illustrated by the significant
movement of the (supposedly) stable community.

The middle column of Figure 1 shows dynamic embeddings based on the temporal
smoothness assumption which act to align subsequent time points. While they
mark an improvement over independent embeddings, they still fail to encode the
trivial system; both independent spectral embedding (ISE) with Procrustes
alignment and GloDyNE (Hou et al., 2020) reduce the movement of the stable
community compared to the independent methods, but they still incorrectly
encode some movement. The fact that state-of-the-art methods, like GloDyNE,
cannot solve such a simple problem leads us to claim that the field of dynamic
embedding is very much in its infancy.

We contrast these unstable methods with two stable methods, which have
guarantees to solve this problem. The top right of Figure 1 shows the
aforementioned UASE (Gallagher et al., 2021). Then, as a demonstration of our
main theorem and to contrast with other skip-gram methods, we define a stable
dynamic version of the popular node2vec embedding (Grover and Leskovec, 2016),
which we call unfolded node2vec. This is displayed in the bottom right of
Figure 1. Both stable methods correctly encode the stability of the static
community, and the movement of the other. Note that the difference between
independent node2vec and unfolded node2vec is only correctly structuring the
input matrix and interpretation of the output.

Figure 1: Example dynamic network embeddings from various methods. Points are
coloured based on their community membership and time point. Dashed lines
illustrate the movement of communities. The subplot row indicates the
embedding method type and the subplot column indicates the temporal mechanism
used. The embedding of the stable community at time 1 should be exchangeable
with its embedding at time 2, as achieved by the unfolded methods. For
visualisation, the skip-gram methods are embedded into three dimensions and
then PCA is used to embed down to two.

##  3 Dilated Unfolded Embedding

In this section, we will set up the problem of dynamic network embedding,
explain what makes a dynamic network embedding stable, introduce the dilated
unfolded adjacency matrix, and prove our main theorem.

###  3.1 Dynamic Network Embedding Setup

An $n$ node network is represented using an adjacency matrix,
$\mathbf{A}\in\left\\{0,1\right\\}^{n\times n}$, where $\mathbf{A}_{ij}=1$ if
there exists an edge between nodes $i$ and $j$ and $\mathbf{A}_{ij}=0$
otherwise. In this paper, we consider the discrete-time dynamic network
setting in which we represent an $n$ node dynamic undirected network as a
series of $T$ symmetric adjacency matrix “snapshots” over time,
$\mathbf{A}^{(1)},\dots,\mathbf{A}^{(T)}\in\left\\{0,1\right\\}^{n\times n}$.
Here, we assume that each $\mathbf{A}^{(t)}$ is an inhomogeneous random graph,
$\mathbf{A}^{(t)}_{ij}\overset{\text{ind}}{\sim}\text{Bernoulli}\left(\mathbf{P}_{ij}^{(t)}\right)$,
where $\mathbf{P}^{(t)}\in[0,1]^{n\times n}$ is a fixed matrix giving the edge
probability for each pair of nodes in the network at time $t$. In contrast to
related works (Jones and Rubin-Delanchy, 2020; Gallagher et al., 2021), we do
not assume that $\mathbf{P}^{(t)}$ is low-rank under this framework. As a
consequence of this relaxation, we have no clear notion of a “correct” $d$,
and can select $d$ by preference, e.g. for interpretation (with $d=2$), to
gain statistical power (discussed in section 4), or for computational
constraints. Then, our goal is to find a $d$-dimensional embedding for each
snapshot,
$\hat{\mathbf{Y}}^{(1)},\dots,\hat{\mathbf{Y}}^{(T)}\in\mathbb{R}^{n\times
d}$, with $d\leq n$.

An additional challenge that comes with dynamic network embedding over the
static network case is to ensure that information is well-represented across
different embedding time points. Specifically, we would like to ensure that if
any two nodes $i$ and $j$ have the same edge probability vectors at different
times $t$ and $u$, that is, $\mathbf{P}_{i}^{(t)}=\mathbf{P}_{j}^{(u)}$, then
their embedded positions should be exchangeable. To do this, building on the
work of Gallagher et al. (2021), we introduce the concepts of spatial and
temporal stability. We assume that for a fixed $n$ only a countable number of
embeddings are possible, so that the embedding is a random object with a
discrete distribution. We make this assumption only to help presentation—for
continuous-valued embeddings, the probabilities below can be changed to
densities.

######  Definition 1 (Spatial stability)

A dynamic embedding, $\hat{\mathbf{Y}}^{(1)},\dots,\hat{\mathbf{Y}}^{(T)}$, is
spatially stable if when $\mathbf{P}_{i}^{(t)}=\mathbf{P}_{j}^{(t)}$ the
positions $\hat{\mathbf{Y}}_{i}^{(t)}$ and $\hat{\mathbf{Y}}_{j}^{(t)}$ are
exchangeable,

$\mathbb{P}\left(\hat{\mathbf{Y}}^{(t)}_{i}=v_{1},\hat{\mathbf{Y}}^{(t)}_{j}=v_{2}\right)=\mathbb{P}\left(\hat{\mathbf{Y}}^{(t)}_{i}=v_{2},\hat{\mathbf{Y}}^{(t)}_{j}=v_{1}\right),$

for any $t\in\\{1,\dots,T\\}$ and any $i,j\in\\{1,\dots,n\\}$.

######  Definition 2 (Temporal stability)

A dynamic embedding, $\hat{\mathbf{Y}}^{(1)},\dots,\hat{\mathbf{Y}}^{(T)}$, is
temporally stable if when $\mathbf{P}_{i}^{(t)}=\mathbf{P}_{i}^{(u)}$ the
positions $\hat{\mathbf{Y}}_{i}^{(u)}$ and $\hat{\mathbf{Y}}_{i}^{(t)}$ are
exchangeable,

$\mathbb{P}\left(\hat{\mathbf{Y}}^{(u)}_{i}=v_{1},\hat{\mathbf{Y}}^{(t)}_{i}=v_{2}\right)=\mathbb{P}\left(\hat{\mathbf{Y}}^{(u)}_{i}=v_{2},\hat{\mathbf{Y}}^{(t)}_{i}=v_{1}\right),$

for any $u,t\in\\{1,\dots,T\\}$ and any $i\in\\{1,\dots,n\\}$.

The only existing dynamic network embedding procedure we know to be spatially
and temporally stable is unfolded adjacency spectral embedding (UASE) (Jones
and Rubin-Delanchy, 2020; Gallagher et al., 2021). UASE constructs a time-
invariant _anchor_ embedding,

$\hat{\mathbf{X}}_{\mathbfcal{A}}=\hat{\mathbf{U}}_{\mathbfcal{A}}\hat{\mathbf{\Sigma}}_{\mathbfcal{A}}^{1/2}\in\mathbb{R}^{n\times
d},$

and a time-varying _dynamic_ embedding,

$\hat{\mathbf{Y}}_{\mathbfcal{A}}=\hat{\mathbf{V}}_{\mathbfcal{A}}\hat{\mathbf{\Sigma}}_{\mathbfcal{A}}^{1/2}\in\mathbb{R}^{nT\times
d},$

where
$\hat{\mathbf{U}}_{\mathbfcal{A}}\hat{\mathbf{\Sigma}}_{\mathbfcal{A}}\hat{\mathbf{V}}_{\mathbfcal{A}}^{\top}$
is a rank-$d$ (truncated) singular value decomposition (SVD) of the unfolded
adjacency matrix,
$\mathbfcal{A}=\left(\mathbf{A}^{(1)},\dots,\mathbf{A}^{(T)}\right)\in\\{0,1\\}^{n\times
nT}$ (column concatenation). $\hat{\mathbf{Y}}{A}$ can then be split into $T$
$\mathbb{R}^{n\times d}$ dynamic embedding time points as
$\hat{\mathbf{Y}}{A}=\left(\hat{\mathbf{Y}}{A}^{(1)};\dots;\hat{\mathbf{Y}}{A}^{(T)}\right)$
(row concatenation).

###  3.2 Dilated Unfolded Embedding

We introduce the dilated unfolded adjacency matrix, which is the symmetric
dilation of the unfolded adjacency matrix,

$\mathbf{A}=\begin{bmatrix}\mathbf{0}&\mathbfcal{A}\\\
\mathbfcal{A}^{\top}&\mathbf{0}\end{bmatrix}\in\\{0,1\\}^{(n+nT)\times(n+nT)}.$
(1)

Almost all network embedding algorithms have some random element, e.g., in the
choice of eigenvectors when the corresponding eigenvalues are equal. To cope
with this, we define a network embedding as a function
$\mathcal{F}:\\{\mathbf{a}\in\\{0,1\\}^{m\times
m}:\mathbf{a}=\mathbf{a}^{\top}\\}\times\Omega\rightarrow\mathbb{R}^{m\times
d}$ which takes an adjacency matrix $\mathbf{a}$ and a seed $\omega\in\Omega$
as input and embeds the $i$th row of $\mathbf{a}$ to the $i$th row of
$\mathcal{F}(\mathbf{a},\omega)$. We will use the shorthand
$\mathcal{F}(\mathbf{A})$ to denote the embedding obtained with a random seed
$W$, that is, $\mathcal{F}(\mathbf{A}):=\mathcal{F}(\mathbf{A},W)$. We let

$\begin{bmatrix}\hat{\mathbf{X}}\\\
\hat{\mathbf{Y}}\end{bmatrix}=\mathcal{F}(\mathbf{A}),$ (2)

where $\hat{\mathbf{X}}\in\mathbb{R}^{n\times d}$, and
$\hat{\mathbf{Y}}\in\mathbb{R}^{nT\times d}$, defining _anchor_ and _dynamic_
embeddings respectively, corresponding to $\mathcal{F}$. $\hat{\mathbf{Y}}$ is
split into $T$ $\mathbb{R}^{n\times d}$ embedding time points,
$\hat{\mathbf{Y}}=\left(\hat{\mathbf{Y}}^{(1)};\dots;\hat{\mathbf{Y}}^{(T)}\right)$.

###### Proposition 3

If $\mathcal{F}$ is a rank-$2d$ adjacency spectral embedding then
$\hat{\mathbf{X}}$ and $\hat{\mathbf{Y}}$ contain scaled rank-$d$ anchor and
dynamic UASEs, that is,

$\begin{bmatrix}\hat{\mathbf{X}}\\\
\hat{\mathbf{Y}}\end{bmatrix}=\mathcal{F}(\mathbf{A})=\hat{\mathbf{U}}|\hat{\mathbf{\Lambda}}|^{1/2}=\frac{1}{\sqrt{2}}\begin{bmatrix}\hat{\mathbf{X}}_{\mathbfcal{A}}&\hat{\mathbf{X}}_{\mathbfcal{A}}\\\
\hat{\mathbf{Y}}_{\mathbfcal{A}}&-\hat{\mathbf{Y}}_{\mathbfcal{A}}\end{bmatrix},$

where $\hat{\mathbf{\Lambda}}\in\mathbb{R}^{2d\times 2d}$ is a diagonal matrix
containing the $2d$ largest eigenvalues of $\mathbf{A}$ in terms of magnitude
in appropriate order and $\hat{\mathbf{U}}\in\mathbb{R}^{(n+nT)\times 2d}$ is
a matrix of appropriately chosen orthonormal eigenvectors (as columns).

###  3.3 Provable Stability of Embedding the Dilated Unfolded Adjacency Matrix

We introduce the important but very mild assumption that $\mathcal{F}$ is a
label-invariant network embedding, that is, that $\mathcal{F}$ is invariant to
the ordering of the nodes in its input adjacency matrix.

######  Definition 4 (Label-invariant network embedding)

A label-invariant network embedding satisfies

$\mathbb{P}\left(\mathcal{F}(\mathbf{a},W)=\mathbf{v}\right)=\mathbb{P}\left(\mathcal{F}(\mathbf{\Pi}\mathbf{a}\mathbf{\Pi}^{\top},W)=\mathbf{\Pi}\mathbf{v}\right),$

for any $\mathbf{a}\in\\{0,1\\}^{m\times m}$,
$\mathbf{v}\in\mathbb{R}^{m\times d}$, and permutation matrix
$\mathbf{\Pi}\in\\{0,1\\}^{m\times m}$.

###### Theorem 5

If $\mathcal{F}$ is label-invariant, then the dynamic embedding
$\mathbf{\hat{Y}}$ obtained from $\mathcal{F}(\mathbf{A})$ is spatially and
temporally stable.

Proof  Let $\mathcal{F}$ be a label-invariant network embedding as defined in
Definition 4. By enforcing that the rows and columns of $\mathbf{A}$ are
ordered as stated in Equation 1, and based on the definition of a label-
invariant network embedding,

$\begin{bmatrix}\hat{\mathbf{X}}\\\
\hat{\mathbf{Y}}\end{bmatrix}=\mathcal{F}(\mathbf{A}).$

Suppose that $\mathbf{P}_{i}^{(u)}=\mathbf{P}_{j}^{(t)}$. By the law of total
probability, we can express the joint probability of the embedding of node $i$
at time $u$ and node $j$ at time $t$ as,

$\mathbb{P}\left(\hat{\mathbf{Y}}_{i}^{(u)}=v_{1},\hat{\mathbf{Y}}_{j}^{(t)}=v_{2}\right)=\sum_{\mathbf{a}\in\mathbb{A}}\mathbb{P}\left(\hat{\mathbf{Y}}^{(u)}_{i}=v_{1},\hat{\mathbf{Y}}^{(t)}_{j}=v_{2}\middle|\mathbf{A}=\mathbf{a}\right)\mathbb{P}\left(\mathbf{A}=\mathbf{a}\right),$

where $\mathbb{A}$ is a sequence containing exactly one instance of every
possible $(n+nT)\times(n+nT)$ dilated unfolded adjacency matrix $\mathbf{A}$
in some order. Let $\mathbf{\Pi}$ be a permutation matrix that permutes node
$i$ at time $u$ with node $j$ at time $t$. Then, as we get $\hat{\mathbf{Y}}$
from a label-invariant $\mathcal{F}$,

$\displaystyle\sum_{\mathbf{a}\in\mathbb{A}}$
$\displaystyle\mathbb{P}\left(\hat{\mathbf{Y}}^{(u)}_{i}=v_{1},\hat{\mathbf{Y}}^{(t)}_{j}=v_{2}\middle|\mathbf{A}=\mathbf{a}\right)\mathbb{P}\left(\mathbf{A}=\mathbf{a}\right)$
$\displaystyle=\sum_{\mathbf{a}\in\mathbb{A}}\mathbb{P}\left(\hat{\mathbf{Y}}^{(u)}_{i}=v_{2},\hat{\mathbf{Y}}^{(t)}_{j}=v_{1}\middle|\mathbf{A}=\mathbf{\Pi
a\Pi}^{\top}\right)\mathbb{P}\left(\mathbf{A}=\mathbf{\Pi
a\Pi}^{\top}\right),$
$\displaystyle=\sum_{\mathbf{\Pi}\mathbf{a}\mathbf{\Pi}^{\top}\in\mathbb{A}}\mathbb{P}\left(\hat{\mathbf{Y}}^{(u)}_{i}=v_{2},\hat{\mathbf{Y}}^{(t)}_{j}=v_{1}\middle|\mathbf{A}=\mathbf{a}\right)\mathbb{P}\left(\mathbf{A}=\mathbf{a}\right).$

By then removing this permutation, we sum over the same terms in a different
order,

$\displaystyle\sum_{\mathbf{\Pi}\mathbf{a}\mathbf{\Pi}^{\top}\in\mathbb{A}}$
$\displaystyle\mathbb{P}\left(\hat{\mathbf{Y}}^{(u)}_{i}=v_{2},\hat{\mathbf{Y}}^{(t)}_{j}=v_{1}\middle|\mathbf{A}=\mathbf{a}\right)\mathbb{P}\left(\mathbf{A}=\mathbf{a}\right)$
$\displaystyle=\sum_{\mathbf{a}\in\mathbb{A}}\mathbb{P}\left(\hat{\mathbf{Y}}^{(u)}_{i}=v_{2},\hat{\mathbf{Y}}^{(t)}_{j}=v_{1}\middle|\mathbf{A}=\mathbf{a}\right)\mathbb{P}\left(\mathbf{A}=\mathbf{a}\right),$
$\displaystyle=\mathbb{P}\left(\hat{\mathbf{Y}}^{(u)}_{i}=v_{2},\hat{\mathbf{Y}}^{(t)}_{j}=v_{1}\right).$

This leaves us with

$\mathbb{P}\left(\hat{\mathbf{Y}}^{(u)}_{i}=v_{1},\hat{\mathbf{Y}}^{(t)}_{j}=v_{2}\right)=\mathbb{P}\left(\hat{\mathbf{Y}}^{(u)}_{i}=v_{2},\hat{\mathbf{Y}}^{(t)}_{j}=v_{1}\right),$

which proves that $\hat{\mathbf{Y}}$ is spatially and temporally stable.

##  4 A Hypothesis Test for Stability

Now that we have a way of generating stable dynamic embeddings from a wide
class of methods, this section will demonstrate why stability is so important,
or rather demonstrate exactly why dynamic embedding methods without stability
cannot be interpreted. We will consider a range of trivial simulated networks
with planted spatial or temporal changes. We then define a hypothesis test
that can determine whether an embedding is encoding the correct structure and
quantify how powerful the embedding is (how clearly it encodes changes). We
would expect that any reasonable dynamic embedding method would perform well
on these very simple systems. If a method were to perform poorly, we argue
that when applied to real data, the embedding would not be interpretable.

###  4.1 Generating Dynamic Networks

We use three simple dynamic stochastic block models (DSBMs) (Yang et al.,
2011; Xu and Hero, 2014) to model static, moving and merging node communities
over time. The DSBM is a simple dynamic network model where nodes are assigned
to one of $K$ communities. Then, the edge probability between any two nodes is
determined based on their community memberships. Specifically, we draw our
adjacency snapshots as

$\mathbf{A}_{ij}^{(t)}\sim\text{Bernoulli}\left(\mathbf{B}^{(t)}_{\tau_{i},\tau_{j}}\right),$

where $\mathbf{B}^{(t)}\in[0,1]^{K\times K}$ is a matrix of inter-community
edge probabilities and $\tau\in\\{1,\dots,K\\}^{n}$ is the community
allocation vector. For example, in the Moving system defined in Table 1, the
second community initially has a within-community edge probability of 0.5,
which then changes slightly to 0.53 in the second snapshot. The first
community has the same within- and inter-community edge probabilities in both
snapshots. Therefore, for this system, we have planted a temporal change in
the second community, but not the first.

As real-world networks are generally sparse, we additionally define systems in
the sparse regime. Following recommendations from Leskovec et al. (2009, Sec.
6.2), we use a Chung Lu model,

$\mathbf{A}_{ij}^{(t)}\sim\text{Bernoulli}\left(\frac{w_{i}w_{j}}{\sum_{k}w_{k}}\right),$

where $\mathbf{w}=\left(w_{1},\dots,w_{n}\right)$ is a vector of node weights
that follows a power-law distribution. These systems are summarised in Table
1.

We have designed these systems to be as simple as possible with at most two
communities at any given time. Therefore, we would expect a reasonable dynamic
embedding method to perform well on them.

System |  DSBM Matrix |  Description
---|---|---
Static |  $\mathbf{B}^{(1)}=\begin{bmatrix}0.5&0.5\\\ 0.5&0.4\end{bmatrix}=\mathbf{B}^{(2)}$ |  Both communities are static
Moving |  $\mathbf{B}^{(1)}=\begin{bmatrix}0.5&0.2\\\ 0.2&0.5\end{bmatrix},\mathbf{B}^{(2)}=\begin{bmatrix}0.5&0.2\\\ 0.2&0.53\end{bmatrix}$ |  Community 2 changes, community 1 is static
Merge |  $\mathbf{B}^{(1)}=\begin{bmatrix}0.9&0.2\\\ 0.2&0.1\end{bmatrix},\mathbf{B}^{(2)}=\begin{bmatrix}0.5&0.5\\\ 0.5&0.5\end{bmatrix}$ |  Community 1 and 2 merge over time
System |  Probability Matrix |  Description
Static with Power Dist. |  $\mathbf{P}^{(1)}_{ij}=w_{i}w_{j}/\sum_{k}{w_{k}}=\mathbf{P}^{(2)}_{ij}$ |  Static sparse network
Moving with Power Dist. |  $\mathbf{P}^{(1)}_{ij}=w_{i}w_{j}/\sum_{k}{w_{k}}$, $\mathbf{P}^{(2)}_{ij}=0.97\left(w_{i}w_{j}/\sum_{k}{w_{k}}\right)$ |  Changing sparse network
Table 1: Generated dynamic network systems used to test the spatio-temporal
properties of a dynamic embedding. In the power example the weights $w_{i}$
are drawn from a power distribution as described in Leskovec et al. (2009,
Sec. 6.2).

###  4.2 The Paired Displacement Test

To work out if a dynamic network embedding has encoded the correct structure,
we must compare the embedding distributions either over time or over space,
depending on the type of planted change. To evaluate whether embedding
distributions have changed, we introduce an exact test. Suppose that we
suspect a change in a set of nodes $\mathscr{N}$, with $|\mathscr{N}|\leq n$,
at some time $t_{c}$. Here, we use $|\cdot|$ to denote the cardinality of a
set. We can test the significance of this change by comparing a window of
embedding time points before and after $t_{c}$. Let
$\mathbf{S}^{[1]}=\hat{\mathbf{Y}}^{(t_{c}-r_{1}:t_{c})}_{\mathscr{N}}\in\mathbb{R}^{|\mathscr{N}|r_{1}\times
d}$ be the embeddings of nodes $\mathscr{N}$ from time $t_{c}-r_{1}$ up to but
not including $t_{c}$, and let
$\mathbf{S}^{[2]}=\hat{\mathbf{Y}}^{(t_{c}:t_{c}+r_{2})}_{\mathscr{N}}\in\mathbb{R}^{|\mathscr{N}|r_{2}\times
d}$ be the embeddings of $\mathscr{N}$ from time $t_{c}$ up to but not
including $t_{c}+r_{2}$. Our hypotheses are then

  * •

$\mathbf{H}_{0}$: $\mathbf{S}^{[1]}$ and $\mathbf{S}^{[2]}$ are drawn from the
same distribution,

  * •

$\mathbf{H}_{1}$: $\mathbf{S}^{[1]}$ and $\mathbf{S}^{[2]}$ are drawn from
different distributions.

Here, our null hypothesis $\mathbf{H}_{0}$ implies that there is no change at
$t_{c}$, whereas the alternative $\mathbf{H}_{1}$ implies that there is a
change at $t_{c}$. We emphasise that the appropriate use of null hypothesis
significance testing is when $\mathbf{H}_{0}$ is plausible, which is true for
simulated data as a test for the validity of embedding procedures.

To test this hypothesis, we define the (temporal) paired displacement test,
with test statistic

$t\left(\mathbf{S}^{[1]},\mathbf{S}^{[2]}\right)=\left\|\sum_{i=1}^{|\mathscr{N}|r_{1}}\mathbf{S}^{[1]}_{i}-\sum_{i=1}^{|\mathscr{N}|r_{2}}\mathbf{S}^{[2]}_{i}\right\|_{2},$
(3)

where $\mathbf{S}^{[1]}_{i}$ is the $i$th row of $\mathbf{S}^{[1]}$. This
gives the L2 norm of the mean node displacement across the two embedding sets.
Our observed test statistic is then
$t_{\text{obs}}=t\left(\mathbf{S}^{[1]},\mathbf{S}^{[2]}\right)$.

For each node we then propose to swap its position with one of its
$r_{1}+r_{2}$ past or future representations with uniform probability, giving
us the permuted embedding $\hat{\mathbf{Y}}_{\mathscr{N}}^{\prime}$. We choose
this permutation as under $\mathbf{H}_{0}$ the nodes in $\mathbf{S}^{[1]}$
should be exchangeable with those in $\mathbf{S}^{[2]}$, and so this
permutation will have no effect on the test in this case. Using the permuted
embedding we can form the permuted embedding matrices,
$\mathbf{S}^{\prime[1]}=\hat{\mathbf{Y}}^{\prime(t_{c}-r_{1}:t_{c})}_{\mathscr{N}}$
and
$\mathbf{S}^{\prime[2]}=\hat{\mathbf{Y}}^{\prime(t_{c}:t_{c}+r_{2})}_{\mathscr{N}}$.
The permuted test statistic is then
$t^{\star}=t\left(\mathbf{S}^{\prime[1]},\mathbf{S}^{\prime[2]}\right)$. Our
p-value is given by how extreme our observed test statistic is against
$n_{\text{sim}}$ permuted test statistics,

$\hat{p}=\frac{1}{n_{\text{sim}}+1}\sum_{i=1}^{n_{\text{sim}}+1}\mathbb{I}\left(\mathbf{T}^{\star}_{i}\geq
t_{\text{obs}}\right),$

where $\mathbf{T}^{\star}\in[0,1]^{n_{\text{sim}}+1}$ is a vector of permuted
test statistics with $\mathbf{T}^{\star}_{1}=t_{\text{obs}}$ and $\mathbb{I}$
is the indicator function. This process is summarised in algorithm 1.

For spatial testing, we can use the same procedure, but instead we compare two
different node sets in a single time point, as opposed to comparing the same
node set across time. In this case we set,
$\mathbf{S}^{[1]}=\hat{\mathbf{Y}}^{(t)}_{\mathscr{N}_{1}}$ and
$\mathbf{S}^{[2]}=\hat{\mathbf{Y}}^{(t)}_{\mathscr{N}_{2}}$, where
$\mathscr{N}_{1}$ and $\mathscr{N}_{2}$ are two different node sets.

Algorithm 1 Temporal Paired Displacement Test

Input:

Dynamic network embedding
$\hat{\mathbf{Y}}^{(1)},\dots,\hat{\mathbf{Y}}^{(T)}\in\mathbb{R}^{n\times d}$

Node set $\mathscr{N}$

Proposed changepoint $t_{c}$ and time windows $r_{1}$ and $r_{2}$ before and
after $t_{c}$

Number of permuted test statistics $n_{\text{sim}}$

Initialize $\mathbf{S}^{[1]}=\hat{\mathbf{Y}}^{(t_{c}-r_{1}:t_{c})}$,
$\mathbf{S}^{[2]}=\hat{\mathbf{Y}}^{(t_{c}:t_{c}+r_{2})}$

Compute
$\mathbf{T}^{\star}_{1}=t_{\text{obs}}=t\left(\mathbf{S}^{[1]},\mathbf{S}^{[2]}\right)$,
as defined in equation 3

for $i=1,\dots,n_{\text{sim}}$ permutations do

for $j=1,\dots,n$ nodes do

Set $\hat{\mathbf{Y}}^{\prime}_{j}$ to be the shuffled version of
$\hat{\mathbf{Y}}_{j}=\left\\{\hat{\mathbf{Y}}_{j}^{(t_{c}-r_{1})},\dots,\hat{\mathbf{Y}}_{j}^{(t_{c}+r_{2})}\right\\}$

end for

Set $\mathbf{S}^{\prime[1]}=\hat{\mathbf{Y}}^{\prime(t_{c}-r_{1}:t_{c})}$,
$\mathbf{S}^{\prime[2]}=\hat{\mathbf{Y}}^{\prime(t_{c}:t_{c}+r_{2})}$

Compute
$t^{\star}=t\left(\mathbf{S}^{\prime[1]},\mathbf{S}^{\prime[2]}\right)$

$\mathbf{T}^{\star}_{i+1}=t^{\star}$

end for

Compute
$\hat{p}=\frac{1}{n_{\text{sim}}+1}\sum_{i=1}^{n_{\text{sim}}+1}\mathbb{I}\left(\mathbf{T}^{\star}_{i}\geq
t_{\text{obs}}\right)$

Output: $\hat{p}$

###  4.3 Methods to Compare

From Proposition 3, we have shown that UASE is a special case of our framework
when $\mathcal{F}$ is an adjacency spectral embedding (ASE). However, from
Theorem 5, stable dynamic embeddings are no longer limited to the use of ASE.
Regularised Laplacian spectral embedding (RLSE) is a static network embedding
that has been shown to perform well on sparse networks (Qin and Rohe, 2013).
By setting $\mathcal{F}$ to be RLSE, we can define _unfolded_ regularised
Laplacian spectral embedding (URLSE),

$\begin{bmatrix}\hat{\mathbf{X}}\\\
\hat{\mathbf{Y}}\end{bmatrix}=\mathcal{F}(\mathbf{A})=\text{ASE}\left(\mathbf{D}_{\gamma}^{-1/2}\mathbf{A}\mathbf{D}_{\gamma}^{-1/2}\right),$

where $\text{ASE}(\mathbf{a})$ indicates an ASE of an adjacency matrix
$\mathbf{a}$, and
$\mathbf{D}_{\gamma}=\mathbf{D}+\gamma\mathbf{I}\in\mathbb{R}^{n\times n}$.
Here, $\mathbf{D}$ is a diagonal degree matrix with entries
$\mathbf{D}_{ii}=\sum_{j}\mathbf{A}_{ij}$ and $\gamma\geq 0$ is a
regularisation parameter. Following recommendations from Qin and Rohe (2013),
we will set $\gamma$ to be the average node degree.

We compare URLSE with three other dynamic spectral embedding methods.
Independent spectral embedding (ISE) represents the naïve case, where each
$\hat{\mathbf{Y}}^{(t)}$ is an independent adjacency spectral embedding of
each $\mathbf{A}^{(t)}$. Note that when using ISE we have ensured that the
eigenvector sum is always positive as this makes the method more consistent.
We also include a version of ISE where subsequent embedding time points are
aligned by a Procrustes rotation. The final spectral method we compare is the
omnibus embedding (OMNI) (Levin et al., 2017). This method constructs a single
block matrix containing pairwise-averaged adjacency snapshots,
$\mathbf{M}_{s,t}=(\mathbf{A}^{(s)}+\mathbf{A}^{(t)})/2$, and then computes an
adjacency spectral embedding on this matrix to achieve a dynamic embedding. It
has been shown that OMNI is temporally stable, albeit at the cost of spatial
stability (Gallagher et al., 2021).

As a demonstration of the flexibility of this framework, we define an
additional unfolded method using a non-spectral and non-deterministic static
network embedding. node2vec (Grover and Leskovec, 2016) is a popular static
network embedding method which uses random walks to generate node
neighbourhoods, then a skip-gram with negative sampling (SGNS) is used to
compute an embedding that preserves these node neighbourhoods. Despite the
lack of statistical rigour behind node2vec, in comparison to that of spectral
embedding, we can still guarantee the stability of unfolded node2vec as static
node2vec satisfies our label-invariance condition. Unfolded node2vec is given
by

$\begin{bmatrix}\hat{\mathbf{X}}\\\
\hat{\mathbf{Y}}\end{bmatrix}=\mathcal{F}(\mathbf{A})=\text{node2vec}(\mathbf{A}).$

Note that we do not imply that node2vec is the best choice of $\mathcal{F}$.
We instead leave this as an example and invite readers to experiment with
different choices of $\mathcal{F}$.

Similar to ISE, we contrast the performance of unfolded node2vec against
independent node2vec, where we compute each $\hat{\mathbf{Y}}^{(t)}$ using
independent static node2vec embeddings. The final dynamic embedding we shall
compare is GloDyNE (Hou et al., 2020). We choose to compare against this
method as, like node2vec, it is based on SGNS, allowing us directly compare
the temporal mechanism used by GloDyNE with unfolding. The method aims to
preserve global topology by only updating a subset of nodes which accumulate
the largest topological changes over subsequent network snapshots. Then, an
incremental SGNS is used to train; to compute an embedding at $t$, its
training is initialised using the pre-trained weights of the SGNS at $t-1$ in
order to keep the time points similar. As discussed in the introduction, this
temporal mechanism is not sufficient for temporal stability as it only
attempts to align subsequent time points. In addition, due to the use of
incremental SGNS, the embedding at $t$ will be distorted to look like the
embedding at $t-1$, which means that its spatial stability is also lost.
Therefore GloDyNE has neither temporal nor spatial stability.

Table 2 provides a summary of each dynamic embedding method that we consider.

Embedding Method |  Embedding Matrix |  Description |  Stability
---|---|---|---
UASE |  $\mathbfcal{A}=\left[\mathbf{A}^{(1)},\dots,\mathbf{A}^{(T)}\right]$ |  Spectral embedding |  Both
URLSE |  $\mathbf{L}=\left(\mathbf{D}_{\text{L}}-\gamma\right)^{-1/2}\mathbfcal{A}\left(\mathbf{D}_{\text{R}}-\gamma\right)^{-1/2}$ |  Spectral embedding |  Both
ISE |  $\mathbf{A}^{(t)}$ |  Spectral embedding |  Spatial
ISE Procrust. |  $\mathbf{A}^{(t)},\mathbf{A}^{(t-1)}$ |  Spectral embedding & Procrustes alignment |  Spatial
OMNI |  $\mathbf{M}_{s,t}=\left(\mathbf{A}^{(s)}+\mathbf{A}^{(t)}\right)/2$ |  Spectral embedding |  Temporal
Unfolded node2vec |  $\mathbf{A}=\begin{bmatrix}\mathbf{0}&\mathbfcal{A}\\\ \mathbfcal{A}^{\top}&\mathbf{0}\end{bmatrix}$ |  SGNS |  Both
Independent node2vec |  $\mathbf{A}^{(t)}$ |  SGNS |  Spatial
GloDyNE |  $\mathbf{A}^{(t)}$, $\mathbf{A}^{(t-1)}$ |  Incremental SGNS & topological updates |  None
Table 2: A description of each dynamic embedding method we consider. To our
knowledge, only unfolded methods have both spatial and temporal stability.

###  4.4 Simulated Testing Setup

We now evaluate how well each dynamic embedding can encode planted spatial and
temporal changes at the full graph, community and node levels. As we can
generate many random graphs with the same underlying distribution using the
systems displayed in Table 1, it is possible to generate a distribution of
p-values for each method. By plotting the cumulative $\hat{p}$ distribution,
we can conclude one of three outcomes. If the cumulative distribution is
uniform (as decided by a Kolmogorov-Smirnov test), then we fail to reject the
null hypothesis; there is no significant change encoded. If it is super-
uniform (above 0.05 at the 5% level), then we reject the null hypothesis; a
change has been encoded. If it is sub-uniform (below 0.05 at the 5% level),
then we fail to reject the null hypothesis, but the test has a high false
negative rate; the embedding is conservative. The height of the curve at the
5% level gives the power of an embedding; a more powerful embedding is better
at encoding changes.

200 p-values were computed for each system. At the full graph level,
$\mathscr{N}$ is the set of all nodes in the network. At the community level,
we set $\mathscr{N}=\mathbb{I}\left(\tau=1\right)$ or
$\mathscr{N}=\mathbb{I}\left(\tau=2\right)$ (the set of all nodes in community
1 or 2) for temporal testing and
$\mathscr{N}_{1}=\mathbb{I}\left(\tau=1\right)$,
$\mathscr{N}_{2}=\mathbb{I}\left(\tau=2\right)$ for spatial testing. Then, the
temporal test at the community and graph level is
$t\left(\hat{\mathbf{Y}}_{\mathscr{N}}^{(1)},\hat{\mathbf{Y}}_{\mathscr{N}}^{(2)}\right)$
and the spatial test is
$t\left(\hat{\mathbf{Y}}_{\mathscr{N}_{1}}^{(2)},\hat{\mathbf{Y}}_{\mathscr{N}_{2}}^{(2)}\right)$.
For node testing, we consider single nodes. As we cannot use a hypothesis test
to compare only two positions, we instead generate 50 time points with the
changepoint $t_{c}=26$. Then, for temporal node testing, we use
$t\left(\hat{\mathbf{Y}}_{i}^{(1,\dots,25)},\hat{\mathbf{Y}}_{i}^{(26,\dots,50)}\right)$
and for spatial node testing we use
$t\left(\hat{\mathbf{Y}}_{i}^{(26,\dots,50)},\hat{\mathbf{Y}}_{j}^{(26,\dots,50)}\right)$
where nodes $i$ and $j$ belong to different communities. In each case, we set
$n_{\text{sim}}=1000$.

For each method, we set $d$ to be the rank of the corresponding noise-free
embedding matrix (see Table 2 for embedding matrices). As GloDyNE requires a
single embedding dimension, we use the same $d$ as used by unfolded methods.
While we believe that these are reasonable decisions, it is still an open
question of exactly what $d$ will be optimum. Whilst we require no concept of
a “correct” $d$, in Appendix C we show some examples of how $d$ can affect the
power of an embedding.

###  4.5 Simulated Testing Results

Tables 3 and 4 display the results of the temporal testing experiments and
tables 5 and 6 display the results of the spatial testing experiments. We
present the results of these experiments as coloured cumulative p-value
distributions. For experiments where we have planted a change, we colour a
distribution green if it encodes the change, grey if the change is not encoded
significantly enough (the distribution is uniform), and red if the method is
conservative. For experiments where there is no planted change, we colour a
distribution green if there is no change encoded (but not conservative), blue
if the method is conservative and red if a change is encoded.

For the moving systems, we can test on either the static community or the
moving community. In these cases, if a method cannot encode the static
community correctly, then we deem the method invalid for testing on the moving
community. This is because if the stability cannot be encoded then the amount
of movement in the moving community becomes immaterial.

We find that all methods without temporal stability fail to achieve a uniform
p-value distribution on the static temporal experiments (Table 3). On the
static graph system (the simplest possible case), these methods are at best
conservative. While a conservative embedding does not represent a failure on
the static temporal systems, this conservative nature significantly reduces
power on the moving temporal experiments; ISE Procrustes even remains
conservative on the moving graph system. In contrast, the methods with
temporal stability (OMNI and unfolded methods) always return a uniform p-value
distribution in the static temporal experiments. As they are never
conservative, the temporally stable methods have more power in comparison to
their unstable counterparts, as demonstrated by the moving systems.

On the power-distributed moving graph, URLSE was the most powerful spectral
embedding. This is a demonstration of how we can take results from the static
network embedding literature, for example that regularised spectral embeddings
perform well on sparse networks (Qin and Rohe, 2013), and apply them to
dynamic embedding through the selection of $\mathcal{F}$.

Experiment |  ISE |  ISE Procrust. |  OMNI |  UASE |  URLSE |  Indep. node2vec |  GloDyNE |  Unfolded node2vec
---|---|---|---|---|---|---|---|---
Static Graph |  |  |  |  |  |  |  |
Static Graph Static Community |  |  |  |  |  |  |  |
Static Graph Static Node |  |  |  |  |  |  |  |
Static Graph with Power Dist. |  |  |  |  |  |  |  |
Moving Graph Static Community |  |  |  |  |  |  |  |
Moving Graph Static Node |  |  |  |  |  |  |  |

Key | Uniform | High False Positive Rate | High False Negative Rate
---|---|---|---

Table 3: Summary of temporal testing experiments on static nodes. The outcome
of each experiment is determined using a plot of the cumulative $\hat{p}$
distribution, which is expected to be uniform for these experiments. We see
that only temporally stable methods (OMNI and unfolded methods) consistently
encode the correct information.

Experiment |  ISE |  ISE Procrust. |  OMNI |  UASE |  URLSE |  Indep. node2vec |  GloDyNE |  Unfolded node2vec
---|---|---|---|---|---|---|---|---
Moving Graph |  |  |  |  |  |  Invalid |  Invalid |
Moving Graph with Power Dist. |  Invalid |  Invalid |  |  |  |  Invalid |  Invalid |
Moving Graph Moving Community |  Invalid |  Invalid |  |  |  |  Invalid |  Invalid |
Moving Graph Moving Node |  Invalid |  |  |  |  |  |  Invalid |

Key | Uniform | High Positive Rate | High False Negative Rate
---|---|---|---

Table 4: Summary of temporal testing experiments with moving nodes. The
outcome of each experiment is determined using a plot of the cumulative
$\hat{p}$ distribution, where a higher curve implies a more powerful
embedding. A cell is left “Invalid” when the method had failed on the static
version of the experiment, listed in Table 3. In the power distributed example
URLSE has the most power, demonstrating its advantage over UASE in sparse
systems.

In the static spatial test, all methods were able to encode a difference
between the two communities. In these small-scale systems, the power of
node2vec was below that of spectral embedding. Despite this, we see that
unfolded node2vec is the method with the most power, being far greater than
that of independent node2vec and GloDyNE. This demonstrates how unfolding can
borrow power from other similar time points, without having to sacrifice
temporal stability. In the merging spatial experiments only the spatially
unstable methods, OMNI and GloDyNE, failed.

Overall, these experiments demonstrate that unstable methods are at best
under-powered or conservative. However, even on trivial systems, they more
often encode incorrect structure. Notably, on the static graph system (the
simplest possible dynamic network), a Procrustes alignment causes the
embedding to be conservative, and GloDyNE encodes incorrect structure. In
contrast, unfolded methods are never conservative. With this consistent
validity, we are free to select an $\mathcal{F}$ that will increase the power
of the embeddings, similar to how we set $\mathcal{F}$ to be RLSE to gain
extra power on the moving sparse system.

Experiment |  ISE |  ISE Procrust. |  OMNI |  UASE |  URLSE |  Indep. node2vec |  GloDyNE |  Unfolded node2vec
---|---|---|---|---|---|---|---|---
Static Graph Community |  |  |  |  |  |  |  |

Key | Uniform | High Positive Rate | High False Negative Rate
---|---|---|---

Table 5: A spatial test on a static system. The outcome of each method is
determined using a plot of the cumulative $\hat{p}$ distribution, with super-
uniform curves implying that the method can correctly separate the two
communities. We see that unfolded node2vec is significantly more powerful than
its skip-gram counterparts. This demonstrates the ability of unfolded methods
to borrow strength across similar snapshots.

Experiment |  ISE |  ISE Procrust. |  OMNI |  UASE |  URLSE |  Indep. node2vec |  GloDyNE |  Unfolded node2vec
---|---|---|---|---|---|---|---|---
Merge Community |  |  |  |  |  |  |  |
Merge Node |  |  |  |  |  |  |  |

Key | Uniform | High False Positive Rate | High False Negative Rate
---|---|---|---

Table 6:  Merged nodes tested for spatial separation. The outcome of each
experiment is determined using a plot of the cumulative $\hat{p}$
distribution, which is expected to be uniform here. OMNI and GloDyNE fail here
as they have temporal mechanisms that compromise spatial stability; the merged
node communities appear separate due to a separation in a previous time point.

##  5 Temporal Analysis of a Flight Network Over the COVID-19 Pandemic

In this section, we consider the OpenSky flight network data set where flights
from 2019 to 2021 are tracked, resulting in a dynamic network of 17,388
airports connected by recorded flight paths over 36 network snapshots (Olive
et al., 2022).

The temporal structure of European airports, encoded by a URLSE, is displayed
in Figure 2. For visualisation, we keep only one dimension of the embedding
and then we plot the mean one-dimensional position of Italian and UK airports,
as well as a mean across all European airports. Before the start of the
COVID-19 pandemic (roughly March 2020 for Europe), the network cycles between
a summer/autumn state and a winter/spring state as seasonal destinations
become more and less popular. As well as this periodic structure, it is
possible to see the two main European waves of the COVID-19 pandemic, with the
first wave being shorter and more disruptive in comparison to the second.
After this second wave, the network returns to its normal periodic behaviour.
This return to normality, after the significant disruption of the pandemic, is
possible due to the temporal stability of URLSE. Supplementary Figure S4 shows
the embeddings of other continents, and Supplementary Figure S3 shows the
spatial structure encoded in a time-invariant anchor URLSE.

Figure 2: A one-dimensional URLSE displaying the average temporal structure of
airports in Italy, the UK and across Europe. The embedding encodes the two
main European waves of COVID-19 as well as an annual periodic structure before
and after the pandemic. We also see that the pandemic hit Italy slightly
before the UK and the rest of Europe.

###  5.1 Temporal Clustering

In order to quantify higher-dimensional changes we perform temporal
clustering. We first construct a temporal dissimilarity matrix
$\mathbf{R}\in\mathbb{R}^{T\times T}$ between each embedding time point, where
we use the temporal paired displacement test statistic as a dissimilarity
measure,
$\mathbf{R}_{ij}=t\left(\hat{\mathbf{Y}}_{\mathscr{N}}^{(i)},\hat{\mathbf{Y}}_{\mathscr{N}}^{(j)}\right)$.
We can then recover the temporal structure present in the network by k-means
clustering on $\mathbf{R}$. We emphasise that this is intended as a visual
guide as to which time-points are similar, rather than an assumption that
there are real underlying clusters.

Figure 3 displays the temporal clustering on both URLSE and UASE. For context,
we additionally plot the number of European countries with travel restrictions
in response to the pandemic (Hale et al., 2021), which can be used as a
measure for the magnitude of disruption caused by the pandemic. We see that
both URLSE and UASE cluster months based on the severity of the pandemic,
having clusters that change as the waves of COVID-19 come and go. However, the
embeddings place a larger emphasis on different qualities of the network.
Given five clusters, UASE will prioritise encoding a difference in pre- and
post-pandemic network behaviour. In contrast, URLSE prioritises representing
the winter/summer periodicity of the network, and instead clusters the post-
pandemic network as a return to normality. This demonstrates the flexibility
of our framework; we can select an $\mathcal{F}$ which is best suited to the
features of the network, without worrying about losing temporal stability.
Supplementary Figure S5 shows the effect of using different numbers of
clusters and varying the embedding dimension.

(a) URLSE.

(b) UASE.

Figure 3: The temporal clustering of 50-dimensional unfolded spectral
embeddings, including only European airports, shown with the number of
European countries with restrictions on travel in response to the COVID-19
pandemic. URLSE focuses more on the periodicity of the network, assigning
periodic clusters both before and after the pandemic. UASE focuses more on the
disruption during the pandemic, assigning clusters for the peaks and troughs
of each wave. UASE does not encode a return to normality after the pandemic.

###  5.2 Comparison of Skip-gram Dynamic Embedding Methods

As an illustration of why stability is important in practice, Figure 4
compares our stable skip-gram dynamic embedding, unfolded node2vec, with an
unstable one, GloDyNE. Currently, GloDyNE cannot embed multi-graphs so these
embeddings are computed on binary adjacency matrices for the purpose of
comparison. For visualisation, we reduce the 150-dimensional embeddings to one
using PCA.

GloDyNE is an example of how the temporal smoothness assumption leads to
artificially smoothed embeddings. The embedding fails not only to encode the
repeating periodic structure but also to encode any significant change for the
pandemic. In stark contrast to this, we see that unfolded node2vec is able to
encode both the annual periodic structure of the network as well as the abrupt
changes that occur as the pandemic hits Italy and the UK (with Italy being hit
slightly earlier than the UK). This is a clear demonstration of how much
detail it is possible to encode in a dynamic embedding which is stable in
comparison to a similar, but unstable, method.

Table 7 lists the computation times for each method to embed the network. In
addition to its increased detail, unfolded node2vec was computed almost seven
times faster than GloDyNE with the same skip-gram parameters. We leave in the
appendix Supplementary Figure S2, which shows how embedding computation times
vary depending on the temporal mechanism used for networks of increasing size.

(a) GloDyNE.

(b) Unfolded node2vec.

Figure 4: Two 150-dimensional skip-gram dynamic embedding methods on the binary flight network (i.e. no multi-edges) reduced to one dimension using PCA. Traces are then given by the average embedding position for airports in the UK, Italy, and Europe. Unfolded node2vec encodes the annual periodic nature of the network as well as abrupt changes for Italy and the UK (Italy slightly before the UK) as they are hit by the pandemic. GloDyNE encodes a smooth continuous change over the entire series, encoding no repeated structure or abrupt changes.  Embedding Method |  Embedding Dimension |  Computation Time |  Description | Stability
---|---|---|---|---
UASE |  50 |  9 s |  Spectral embedding | Both
URLSE |  50 |  12 s |  Spectral embedding | Both
Unfolded node2vec |  150 |  54 min 48 s |  SGNS | Both
GloDyNE |  150 |  6 hrs 10 min 50 s |  Incremental SGNS & topological updates | None
Table 7: Computation times to embed the dynamic flight network. The network
has 17,388 nodes and 36 time points. Embeddings were computed using an AMD
Ryzen 5 3600 processor.

##  6 Conclusion

In this paper, we have solved the problem of how to adapt a wide class of
established static network embedding methods to produce stable dynamic network
embeddings. By embedding the dilated unfolded adjacency matrix we gain
stability without the need for a central limit theorem and without requiring
the concept of a “correct” embedding dimension. This can be trivially
implemented for any label-invariant network embedding. Upon the foundation of
stability that is guaranteed by unfolding, we can then select a static
embedding, $\mathcal{F}$, to gain statistical power. For example, we
demonstrated that URLSE was more powerful than UASE on simulated sparse
networks. Due to the flexibility of this framework, we expect that as the
field of static network embedding progresses, the field of dynamic embedding
is able progress at the same rate through the selection of $\mathcal{F}$.

This paper additionally introduced a hypothesis test which can be used to
evaluate how well dynamic embedding methods encode structure. Even in the
simplest possible case, unstable dynamic embedding methods are at best
conservative, but more often encode incorrect structure. This test provides a
previously unseen view of how embeddings encode structure and it highlights
why intuitive alignment-based methods are fundamentally flawed for comparing
embeddings across time. This test also highlights the significant variation
present in unstable embeddings when the selection of $d$ changes, something
that is not a problem with unfolded embedding. We finally demonstrated the
practical advantage that unfolded methods possess over unstable methods by
applying them to a pandemic-disrupted dynamic flight network. Here, only the
unfolded methods could encode the periodicity and abrupt changes present in
the network.

Future research includes experimenting with different selections of
$\mathcal{F}$ and determining what embedding dimension will produce the most
powerful embedding. There is also future research available into what
downstream analysis we can perform on stable dynamic embeddings now that we
have a guarantee of exchangeability.

##  Appendix A Code Availability

Python code for reproducing all embeddings and experiments can be found at
https://github.com/edwarddavis1/universal_dynamic_embedding_with_testing.

##  Appendix B Proof of Proposition 3

Proof  Let $\mathbfcal{A}$ have singular values $\sigma_{i}$ and left and
right singular vectors $u_{i}$, $v_{i}$ respectively for each
$i\in\\{1,\dots,d\\}$. It is then true that its symmetric dilation (the
dilated unfolded adjacency matrix) $\mathbf{A}$ will have eigenvalues
$\pm\sigma_{i}$ with corresponding eigenvectors
$\left\\{\frac{1}{\sqrt{2}}(u_{i},\pm v_{i})\right\\}$ (Horn and Johnson,
2012), meaning that the spectrum of $\mathbfcal{A}$ is completely encoded in
$\mathbf{A}$. Then, for an appropriately chosen rank-$2d$ (truncated)
eigendecomposition of a $\mathbf{A}$, we have

$\displaystyle\mathbf{A}\approx\hat{\mathbf{U}}_{\mathbf{A}}\hat{\mathbf{\Lambda}}_{\mathbf{A}}\hat{\mathbf{U}}_{\mathbfcal{A}}^{\top}=\frac{1}{\sqrt{2}}\begin{bmatrix}\hat{\mathbf{U}}_{\mathbfcal{A}}&\hat{\mathbf{U}}_{\mathbfcal{A}}\\\
\hat{\mathbf{V}}_{\mathbfcal{A}}&-\hat{\mathbf{V}}_{\mathbfcal{A}}\end{bmatrix}\begin{bmatrix}\hat{\mathbf{\Sigma}}_{\mathbfcal{A}}&\mathbf{0}\\\
\mathbf{0}&-\hat{\mathbf{\Sigma}}_{\mathbfcal{A}}\end{bmatrix}\frac{1}{\sqrt{2}}\begin{bmatrix}\hat{\mathbf{U}}_{\mathbfcal{A}}&\hat{\mathbf{U}}_{\mathbfcal{A}}\\\
\hat{\mathbf{V}}_{\mathbfcal{A}}&-\hat{\mathbf{V}}_{\mathbfcal{A}}\end{bmatrix}^{\top}.$

An appropriately chosen spectral embedding of $\mathbf{A}$ into
$\mathbb{R}^{2d}$ then contains scaled rank-$d$ anchor and dynamic UASEs,

$\displaystyle\hat{\mathbf{U}}_{\mathbf{A}}|\hat{\mathbf{\Lambda}}_{\mathbf{A}}|^{1/2}$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{bmatrix}\hat{\mathbf{U}}_{\mathbfcal{A}}&\hat{\mathbf{U}}_{\mathbfcal{A}}\\\
\hat{\mathbf{V}}_{\mathbfcal{A}}&-\hat{\mathbf{V}}_{\mathbfcal{A}}\end{bmatrix}\left|\begin{bmatrix}\hat{\mathbf{\Sigma}}_{\mathbfcal{A}}&\mathbf{0}\\\
\mathbf{0}&-\hat{\mathbf{\Sigma}}_{\mathbfcal{A}}\end{bmatrix}\right|^{1/2}$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{bmatrix}\hat{\mathbf{U}}_{\mathbfcal{A}}\hat{\mathbf{\Sigma}}_{\mathbfcal{A}}^{1/2}&\hat{\mathbf{U}}_{\mathbfcal{A}}\hat{\mathbf{\Sigma}}_{\mathbfcal{A}}^{1/2}\\\
\hat{\mathbf{V}}_{\mathbfcal{A}}\hat{\mathbf{\Sigma}}_{\mathbfcal{A}}^{1/2}&-\hat{\mathbf{V}}_{\mathbfcal{A}}\hat{\mathbf{\Sigma}}_{\mathbfcal{A}}^{1/2}\end{bmatrix}$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{bmatrix}\hat{\mathbf{X}}_{\mathbfcal{A}}&\hat{\mathbf{X}}_{\mathbfcal{A}}\\\
\hat{\mathbf{Y}}_{\mathbfcal{A}}&-\hat{\mathbf{Y}}_{\mathbfcal{A}}\end{bmatrix}.$

##  Appendix C Embedding Dimension and Testing Power

In Section 4, we tested various dynamic embedding methods on simple simulated
networks. In this testing, we demonstrated that unstable methods were
unreliable as they often either encoded incorrect structure or were
conservative. An assumption that we made in this testing was that the
embeddings would be best when the embedding dimension, $d$, matched the rank
of the method’s noise-free embedding matrix. However, when applying these
methods to real-world networks, the rank of the noise-free embedding matrix
will not be known. Further, this choice of $d$ may not optimum for a given
method. Therefore, we repeat two of the experiments from Section 4 using a
range of embedding dimensions. Tables 8 and 9 show the results for the two
temporal smoothing methods, ISE Procrustes and GloDyNE, along with one
unfolded method, URLSE. The rank of these systems is two.

We see that across both experiments, increasing $d$ beyond the rank of the
system reduces power. For the temporal smoothing methods, this reduction
pushes the methods into being conservative, even when a change is present.
While this aids GloDyNE in the static community experiment (in which it fails
at $d=2$), we see that in the moving experiments, both temporal smoothing
methods fail to detect a change when $d$ is sufficiently over-estimated, due
to their conservative nature. This highlights a key strength of our
construction, where the stability of unfolded methods holds at any $d$. We see
that URLSE always maintains stability in the static community system and still
manages to encode a significant change at a greatly over-estimated $d$ in the
moving community system.

Method |  $d=2$ |  $d=3$ |  $d=4$ |  $d=5$
---|---|---|---|---
ISE Procrust. |  |  |  |

Method |  $d=2$ |  $d=4$ |  $d=6$ |  $d=8$
---|---|---|---|---
GloDyNE |  |  |  |

Method |  $d=2$ |  $d=4$ |  $d=8$ |  $d=20$
---|---|---|---|---
URLSE |  |  |  |

Key | Uniform | High False Positive Rate | High False Negative Rate
---|---|---|---

Table 8: A repeat of the moving graph static community experiment for a range
of $d$. At higher dimensions, the temporal smoothing methods, ISE Procrustes
and GloDyNE, become more conservative. In contrast, the stability guarantee
for URLSE holds for any $d$.

Method |  $d=2$ |  $d=3$ |  $d=4$ |  $d=5$
---|---|---|---|---
ISE Procrust. |  |  |  |

Method |  $d=2$ |  $d=4$ |  $d=6$ |  $d=8$
---|---|---|---|---
GloDyNE |  |  |  |

Method |  $d=2$ |  $d=4$ |  $d=8$ |  $d=20$
---|---|---|---|---
URLSE |  |  |  |

Key | High Positive Rate | High False Negative Rate
---|---|---

Table 9: Testing at other embedding dimensions on the moving graph moving
community system. Similar to Table 8, the alignment-based embedding methods
become more conservative at higher dimensions to the point where changes
cannot be detected when $d$ is sufficiently over-estimated. In contrast,
unfolded methods have more power, allowing changes to be detected even at a
significantly over-estimated $d$.

While unfolded methods will be stable at any $d$, it is clear that there
exists an optimum $d$ which produces an embedding with the most power.
However, as our construction does not assume that $\mathbf{P}^{(t)}$ is low-
rank, the notion of what the optimum $d$ is becomes less clear. Consider a
DSBM system with form,

$\mathbf{B}^{(1)}=\begin{bmatrix}p_{1}&0.2&\dots&0.2\\\
0.2&p_{2}&\ddots&\vdots\\\ \vdots&\ddots&\ddots&0.2\\\
0.2&\dots&0.2&p_{K}\end{bmatrix},\mathbf{B}^{(2)}=\begin{bmatrix}p_{1}&0.2&\dots&0.2\\\
0.2&p_{2}+0.03&\ddots&\vdots\\\ \vdots&\ddots&\ddots&0.2\\\
0.2&\dots&0.2&p_{K}\end{bmatrix},$

where $p_{1},\dotsm p_{K}$ are the within-community connection probabilities.
This is an example of a $K$-dimensional random network which has a planted
temporal change in the second community. Using DSBMs of this form allows us to
investigate how the optimum $d$ varies with systems of increasing rank. Figure
S1 displays the testing power of two temporally stable dynamic embedding
methods, UASE and OMNI, on two $8$-dimensional systems with

  * •

system (a): $p_{1},\dots,p_{K}=0.5,\dots,0.5$,

  * •

system (b): $p_{1},\dots,p_{K}=0.3,\dots,0.9$.

While the choice of these two systems is quite arbitrary, the relationship
between the embedding power and $d$ for the two systems is quite different,
particularly in the region of $d<K$. In addition to this, neither system had
the rank of the respective noise-free embedding matrices as the optimum $d$.
In fact, the optimum $d$ in both systems was below $K$. Therefore, it is
generally unclear which $d$ will be best for a given network. We leave the
precise explanation of this behaviour as an open question for future research.

(a) Using a DSBM where the within-community edge probabilities are identical
for each stable cluster in both time points.

(b) Using a DSBM where within-community edge probabilities are different for
each stable cluster in both time points.

Figure S1: Testing for change by a moving community in two 8-community DSBMs.
In the region below the rank of the embedding matrices, the behaviour of this
curve is wildly different for the two systems. In both cases the optimum $d$
was below the rank of the noise-free embedding matrices. While at $d<K$, UASE
and OMNI have similar testing power, the testing power from OMNI drops off
more gradually than UASE.

##  Appendix D Supplimentary Figures

Figure S2: Computation times for three different dynamic skip-gram embedding
methods on the power-distributed moving system for an increasing number of
nodes. While GloDyNE is initially the fastest for low $n$, its computation
time increases significantly for larger $n$. As $n$ increases, unfolded
node2vec is the fastest here. Figure S3: A five-dimensional time-invariant
anchor embedding of the flight network using URLSE. The embedding features
rays from the origin, where continents lie on rays of different angles. The
further from the origin a node is, the higher its degree in the network.
Figure S4: The temporal structure of the flight network for the different
continents encoded by a one-dimensional URLSE. Each continent displays a
significant change around the start of the COVID-19 pandemic. From this
initial change, the continents display various behaviour, for example, Europe
and Oceania had two clear waves of disruption, with the second Oceania wave
much later than the second European wave.

(a) URLSE d=50.

(b) UASE d=50.

(c) URLSE d=100.

(d) UASE d=100.

Figure S5: Temporal clustering on flight network URLSEs and UASEs at 50 and
100 dimensions using different numbers of clusters, $K$. Each month is
assigned to a cluster, which is indicated by a coloured bar. URLSE focuses
more on the periodicity of the network; in each case it encodes periodic
behaviour before the pandemic and often returns to normality after the
pandemic. In contrast, UASE focuses more on the disruption of the pandemic; it
consistently allocates single clusters to the peak and trough of each wave of
the pandemic.

## References

  * Gallagher et al. (2021) Ian Gallagher, Andrew Jones, and Patrick Rubin-Delanchy.  Spectral embedding for dynamic networks with stability guarantees.  _Advances in Neural Information Processing Systems_ , 34:10158–10170, 2021.
  * Goyal et al. (2018) Palash Goyal, Nitin Kamra, Xinran He, and Yan Liu.  Dyngem: Deep embedding method for dynamic graphs.  _arXiv preprint arXiv:1805.11273_ , 2018.
  * Grover and Leskovec (2016) Aditya Grover and Jure Leskovec.  node2vec: Scalable feature learning for networks.  In _Proceedings of the 22nd ACM SIGKDD international conference on Knowledge discovery and data mining_ , pages 855–864, 2016.
  * Hale et al. (2021) Thomas Hale, Noam Angrist, Rafael Goldszmidt, Beatriz Kira, Anna Petherick, Toby Phillips, Samuel Webster, Emily Cameron-Blake, Laura Hallas, Saptarshi Majumdar, et al.  A global panel database of pandemic policies (oxford covid-19 government response tracker).  _Nature human behaviour_ , 5(4):529–538, 2021\.
  * Horn and Johnson (2012) Roger A. Horn and Charles R. Johnson.  _Positive Definite and Semidefinite Matrices_ , page 425–516.  Cambridge University Press, 2 edition, 2012.  doi: 10.1017/CBO9781139020411.010.
  * Hou et al. (2020) Chengbin Hou, Han Zhang, Shan He, and Ke Tang.  Glodyne: Global topology preserving dynamic network embedding.  _IEEE Transactions on Knowledge and Data Engineering_ , 2020.
  * Jones and Rubin-Delanchy (2020) Andrew Jones and Patrick Rubin-Delanchy.  The multilayer random dot product graph.  _arXiv preprint arXiv:2007.10455_ , 2020.
  * Leskovec et al. (2009) Jure Leskovec, Kevin J Lang, Anirban Dasgupta, and Michael W Mahoney.  Community structure in large networks: Natural cluster sizes and the absence of large well-defined clusters.  _Internet Mathematics_ , 6(1):29–123, 2009.
  * Levin et al. (2017) Keith Levin, Avanti Athreya, Minh Tang, Vince Lyzinski, Youngser Park, and Carey E Priebe.  A central limit theorem for an omnibus embedding of multiple random graphs and implications for multiscale network inference.  _arXiv preprint arXiv:1705.09355_ , 2017.
  * Mahdavi et al. (2018) Sedigheh Mahdavi, Shima Khoshraftar, and Aijun An.  dynnode2vec: Scalable dynamic network embedding.  In _2018 IEEE international conference on big data (Big Data)_ , pages 3762–3765. IEEE, 2018.
  * Mikolov et al. (2013) Tomas Mikolov, Kai Chen, Greg Corrado, and Jeffrey Dean.  Efficient estimation of word representations in vector space.  _arXiv preprint arXiv:1301.3781_ , 2013.
  * Olive et al. (2022) Xavier Olive, Martin Strohmeier, and Jannis Lübbe.  Crowdsourced air traffic data from The OpenSky Network 2020 [CC-BY], January 2022.  URL https://doi.org/10.5281/zenodo.5815448.
  * Qin and Rohe (2013) Tai Qin and Karl Rohe.  Regularized spectral clustering under the degree-corrected stochastic blockmodel.  _Advances in neural information processing systems_ , 26, 2013.
  * Rosenthal et al. (2018) Gideon Rosenthal, František Váša, Alessandra Griffa, Patric Hagmann, Enrico Amico, Joaquín Goñi, Galia Avidan, and Olaf Sporns.  Mapping higher-order relations between brain structure and function with embedded vector representations of connectomes.  _Nature communications_ , 9(1):2178, 2018.
  * Singer et al. (2019) Uriel Singer, Ido Guy, and Kira Radinsky.  Node embedding over temporal graphs.  _arXiv preprint arXiv:1903.08889_ , 2019.
  * Xu and Hero (2014) Kevin S Xu and Alfred O Hero.  Dynamic stochastic blockmodels for time-evolving social networks.  _IEEE Journal of Selected Topics in Signal Processing_ , 8(4):552–562, 2014.
  * Yang et al. (2011) Tianbao Yang, Yun Chi, Shenghuo Zhu, Yihong Gong, and Rong Jin.  Detecting communities and their evolutions in dynamic social networks—a bayesian approach.  _Machine learning_ , 82:157–189, 2011.
  * Zhou et al. (2018) Lekui Zhou, Yang Yang, Xiang Ren, Fei Wu, and Yueting Zhuang.  Dynamic network embedding by modeling triadic closure process.  In _Proceedings of the AAAI conference on artificial intelligence_ , volume 32, 2018.
  * Zhou et al. (2019) Yujing Zhou, Weile Liu, Yang Pei, Lei Wang, Daren Zha, and Tianshu Fu.  Dynamic network embedding by semantic evolution.  In _2019 International Joint Conference on works (IJCNN)_ , pages 1–8. IEEE, 2019.


(130.0130) Integrated optics, infrared; (130.6010) Semiconductors; (050.5524)
Subwavelength structures; (050.2770) Gratings.

# Suspended silicon waveguides for long-wave infrared wavelengths

J. Soler Penadés  Optoelectronics Research Centre, University of Southampton,
Southampton SO17 1BJ, United Kingdom  Corresponding author:
<EMAIL_ADDRESS>A. Sánchez-Postigo  Universidad de Málaga, Dept. de
Ingeniería de Comunicaciones, ETSI Telecomunicación, Campus de Teatinos, 29071
Málaga, Spain  M. Nedeljkovic  Optoelectronics Research Centre, University of
Southampton, Southampton SO17 1BJ, United Kingdom  A. Ortega-Moñux
Universidad de Málaga, Dept. de Ingeniería de Comunicaciones, ETSI
Telecomunicación, Campus de Teatinos, 29071 Málaga, Spain  J. G. Wangüemert-
Pérez  Universidad de Málaga, Dept. de Ingeniería de Comunicaciones, ETSI
Telecomunicación, Campus de Teatinos, 29071 Málaga, Spain  Y. Xu
Optoelectronics Research Centre, University of Southampton, Southampton SO17
1BJ, United Kingdom  R. Halir  Universidad de Málaga, Dept. de Ingeniería de
Comunicaciones, ETSI Telecomunicación, Campus de Teatinos, 29071 Málaga, Spain
Z. Qu  Optoelectronics Research Centre, University of Southampton, Southampton
SO17 1BJ, United Kingdom  A. Z. Khokhar  Optoelectronics Research Centre,
University of Southampton, Southampton SO17 1BJ, United Kingdom  A. Osman
Optoelectronics Research Centre, University of Southampton, Southampton SO17
1BJ, United Kingdom  W. Cao  Optoelectronics Research Centre, University of
Southampton, Southampton SO17 1BJ, United Kingdom  C.G. Littlejohns
Optoelectronics Research Centre, University of Southampton, Southampton SO17
1BJ, United Kingdom  NOVITAS, School of Electrical and Electronic Engineering,
Nanyang Technological University, Singapore 639798, Republic of Singapore  P.
Cheben  National Research Council Canada, Building M-50, Ottawa, K1A 0R6
Canada  I. Molina-Fernández  Universidad de Málaga, Dept. de Ingeniería de
Comunicaciones, ETSI Telecomunicación, Campus de Teatinos, 29071 Málaga, Spain
G. Z. Mashanovich  Optoelectronics Research Centre, University of Southampton,
Southampton SO17 1BJ, United Kingdom

###### Abstract

In this paper we report suspended silicon waveguides operating at a wavelength
of 7.67 $\mu$m with a propagation loss of $3.1\pm 0.3$ dB/cm. To our knowledge
this is the first demonstration of low loss silicon waveguides at such a long
wavelength, with loss comparable to other platforms that use more exotic
materials. The suspended Si waveguide core is supported by a sub-wavelength
grating that provides lateral optical confinement, while also allowing access
to the buried oxide layer so that it can be wet etched using hydrofluoric
acid. We also demonstrate low loss waveguide bends and s-bends.

††journal: ol

Mid-infrared (MIR) photonics research has experienced rapid development in the
last few years, which has been sparked by increased interest in the
commercial, military and security applications of the field, and by much
improved availability of characterization and testing equipment that is able
to reach longer MIR wavelengths. Some of these applications, such as
absorption spectroscopy of gases and chemicals in the fingerprint region,
require the use of wavelengths above 7 $\mu$m [1] (e.g. Sulfur dioxide has
strong absorption near 7.6 $\mu$m [2]), and there is therefore a need to
develop photonic integrated circuits that can work throughout the MIR.

Silicon-on-insulator (SOI) is the dominant waveguide platform in silicon
photonics, but in the MIR alternative platforms are needed because SiO2 has
extremely high absorption at wavelengths above 4 $\mu$m (i.e. higher than 2.5
$\times 10^{3}$ dB/cm at 7.6 $\mu$m, rising to $\sim$2.5 $\times 10^{4}$ dB/cm
at 8.0 $\mu$m) [3]. Several alternative MIR group-IV material waveguide
platforms with wider transparency ranges have been reported in the last few
years, such as silicon-on-sapphire [4, 5, 6, 7], silicon-on-nitride [8],
germanium-on-nitride [9], graded Si/Ge [10, 11], and germanium-on-silicon
(GOS) [12, 13, 14]. So far group-IV material waveguiding has been demonstrated
up to the 7.5-8.5 $\mu$m range, where graded Si/Ge [15, 11] and GOS [16] have
exhibited minimum propagation losses below 3 dB/cm.

An alternative approach to changing the material is to use SOI wafers, but to
limit the optical mode overlap with SiO2, such as by suspending the Si
waveguide core or by substantially increasing the Si core size [13]. Compared
to other platforms SOI wafers are easily available and the fabrication
techniques are very well known, and as such at MIR wavelengths where Si is
still transparent it may be preferred. In [17, 18, 19] two dry-etch steps were
used to first define a Si (or Ge) rib waveguide, and then to create an array
of holes away from the waveguide mode that enable the buried oxide (BOX) to be
locally removed by wet etching in hydrofluoric acid (HF). In [20, 21], we
presented an alternative approach, consisting of a strip waveguide with
subwavelength grating (SWG) holes that provide lateral optical confinement as
well as access to the BOX. The concept of subwavelength grating (SWG)
engineering of Si-wire waveguides was first proposed in [22, 23] and, more
recently, extensively used in a variety of photonic devices [24]. With this
approach we achieved a minimum loss of 0.82 dB/cm at 3.8 $\mu$m. This design
allows simpler fabrication consisting of a single dry etch step and it has the
further advantage of reducing the necessary suspended region width and
improving the mechanical stability, allowing the definition of wider devices
like tapers and MMIs [20]. The same approach has recently been implemented to
create suspended Si slot waveguides [25] at 2.25 $\mu$m.

We note that a further advantage of this waveguide platform compared to GOS
and graded Si/Ge is the symmetrical high refractive index contrast between the
core and upper and lower claddings. For example, GOS exhibits high index
contrast between the Ge core (n = 4.0) and air upper cladding (n = 1.0), but
low contrast with the Si lower cladding (n = 3.4), which results in the need
for single mode waveguides to be rib waveguides with thickness on the order of
several micrometers, and also results in very little mode overlap with the air
upper cladding. Suspended Si waveguides would have a thinner waveguide core
and larger evanescent field, and therefore much increased sensitivity for
sensing applications where high mode overlap with an analyte placed on the
waveguide surface is desirable.

To our knowledge, suspended Si waveguides have not yet been demonstrated at
wavelengths beyond 3.8 $\mu$m. In this paper we report suspended silicon
waveguides with 3.1 dB/cm propagation losses, 90o bends and s-bends at
$\lambda$ = 7.67 $\mu$m. At this wavelength bulk Si has a material absorption
of $\sim$2.1 dB/cm [26], increasing to 4.6 dB/cm at 9 $\mu$m and 9.3 dB/cm at
11 $\mu$m, thus the demonstrated waveguides in this letter approach the
transparency limit of Si.

A schematic representation of the proposed suspended waveguide is shown in
Fig. 1. This guiding structure, based on a SOI platform, comprises a solid
silicon core of thickness $t_{Si}$ and width $W_{core}$ with a SWG structure
of width $W_{clad}$ acting as the waveguide cladding [27]. This lattice with
holes and silicon strips of lengths $L_{hole}$ and $L_{Si}$, respectively,
behaves as a homogeneous medium with an equivalent refractive index $n_{SWG}$,
offering the index contrast required for waveguiding. In addition, the
cladding is used to anchor the waveguide core to the unsuspended lateral
silicon areas and to allow the entrance of the HF acid needed to remove the
$t_{BOX}$-thick buried oxide layer.

Figure 1: Schematic representation (front and top views) of the proposed
suspended Si waveguide with subwavelength grating cladding.

In our previous design working at $\lambda$ = 3.8 $\mu$m [20], we used a SOI
wafer with $t_{Si}$ = 500 nm and $t_{BOX}$ = 3 $\mu$m. The waveguide core
width was $W_{core}$ = 1.3 $\mu$m. The SWG cladding parameters ($W_{clad}$ =
2.5 $\mu$m, $L_{Si}$ = 100 nm and $L_{hole}$ = 450 nm) were optimized to
provide mechanical stability to the entire structure and, simultaneously,
facilitate the HF flow through the lattice holes. Scaling all dimensions by a
factor of 2 could be acceptable to move from $\lambda$ = 3.8 $\mu$m to
$\lambda$ = 7.67 $\mu$m, but this would imply $t_{BOX}$ = 6 $\mu$m. However,
as in the long-wave range we use a SOI wafer with the same SiO2 layer as
before, $t_{BOX}$ thickness remains fixed to 3 $\mu$m. If all dimensions
except $t_{BOX}$ are doubled, the power leakage to the Si substrate could lead
to excessive propagation losses. Hence, the waveguide must be completely
redesigned.

We follow a design methodology based on two separate steps. Firstly, the
silicon thickness $t_{Si}$ is chosen to have vertical single-mode operation
with negligible leakage to the substrate. Secondly, both the core width
$W_{core}$ and the SWG cladding ($W_{clad}$, $L_{Si}$ and $L_{hole}$) are
adjusted to mitigate power leakage to the lateral Si without under-etching.

To study the vertical and lateral leakage losses of the TE-polarized modes
supported by the waveguides, we used Rsoft™ and Photon Design® simulation
suites in combination with our 2D in-house electromagnetic tool FEXEN [28],
which is especially appropriate for designing periodic structures such as
grating couplers and the SWG cladding of the waveguides.

Choosing the silicon thickness is critical to enable waveguiding with
minimized propagation losses. A high value of $t_{Si}$ is beneficial to
vertically confine the fundamental mode with low leakage to the substrate. On
the other hand, a reduced silicon thickness is desirable to cut off the
second-order vertical guided mode, or, at least, to ensure it leaks towards
the substrate in the shortest possible propagation distance. In our initial
simulations we kept the same SWG equivalent refractive index $n_{SWG}=1.72$ as
in our original design [20]. For the fundamental mode, we used a waveguide
width $W_{core}$ = 2.6 $\mu$m, which coincides with the scaled value we
employed at $\lambda$ = 3.8 $\mu$m. This value of 2.6 $\mu$m will be later
optimized in the last step of the design process. Regarding the second-order
vertical mode, it must be noticed that this mode should be avoided not only in
interconnecting waveguides, but also in much wider optical devices like
multimode interferometer (MMI) couplers. Since the vertical confinement for a
wide waveguide is always higher than for a narrow waveguide, we have decided
to calculate the effective index of the second-order vertical mode for a slab
(i.e. $W_{core}\rightarrow\infty$) as a worst case. This provides us an upper
bound, since any finite-width waveguide will exhibit a weaker guidance and it
will therefore suffer higher substrate leakage losses. Figure 2(a) shows the
variation of the real part of the effective indices of the fundamental and
second-order vertical modes as a function of the waveguide core thickness
$t_{Si}$. For $t_{Si}$ > $\sim$1 $\mu$m, the fundamental quasi-TE mode is well
confined. As is observed, the second-order vertical mode is only weakly guided
for Si thicknesses between $t_{Si}$ = 1.2 $\mu$m and $t_{Si}$ = 1.5 $\mu$m, as
its low effective index, near 1 (refractive index of air), denotes the field
is substantially expanded outside the waveguide silicon core. So, the second-
order mode suffers from high vertical leakage in the depicted range. Figure
2(a) (inset) shows, for $t_{Si}$ = 1.4 $\mu$m, the transverse distributions of
both the fundamental mode ($W_{core}$= 2.6 $\mu$m) and second-order vertical
mode ($W_{core}\rightarrow\infty$).

Figure 2: (a) Real part of the effective indices ($n_{eff}$) vs. the Si
thickness ($t_{Si}$) for the fundamental mode ($W_{core}$ = 2.6 $\mu$m) and
the second-order vertical mode ($W_{core}\rightarrow\infty$) supported by the
suspended waveguides. Simulation data: $t_{BOX}$ = 3 $\mu$m, $n_{SWG}$ = 1.72.
The real part of the effective index of the final designed waveguide with
$W_{core}$ = 2.9 $\mu$m is also included. Insets: transverse field
distributions for the first two modes ($W_{core}$ = 2.6 $\mu$m,
$W_{core}\rightarrow\infty$) when $t_{Si}$ = 1.4 $\mu$m. (b) Power loss of the
fundamental vertical mode of the suspended waveguide with $W_{core}$ = 2.6
$\mu$m and the finally designed waveguide as a function of the waveguide
silicon thickness $t_{Si}$.

To design an exact $t_{Si}$ value, a further examination of the propagation
losses $\alpha$ is required. They have two contributions: the intrinsic
material loss, which is $\alpha_{Si}\approx$ 2.1 dB/cm at $\lambda$ = 7.67
$\mu$m [26], and the vertical power leakage to the substrate. In our
simulations, we included silicon losses in the waveguide core and substrate.
The total loss of the fundamental mode of a waveguide with $W_{core}$ = 2.6
$\mu$m is shown in Fig. 2(b) as a function of the core thickness $t_{Si}$. For
Si thicknesses smaller than $t_{Si}$ $\sim$1 $\mu$m, substrate leakage losses
are dominant and $\alpha$ exhibits an exponential increase. On the contrary,
for larger thicknesses, propagation losses asymptotically tend to the material
absorption $\alpha_{Si}$, as the mode is well confined for higher $t_{Si}$
values (> 1.3 $\mu$m). Furthermore, to increase the robustness of the design,
the sensitivity of $\alpha$ with respect to $t_{Si}$ at the selected nominal
Si thickness should be close to zero; otherwise, propagation losses could be
rapidly increased if undesired reduction in the silicon thickness is produced
by a HF overexposure during the fabrication process. Taking these
considerations into account, we selected $t_{Si}$ = 1.4 $\mu$m. For this
value, the fundamental mode leakage to substrate is negligible, more than ten
times smaller than the intrinsic silicon losses. Although for the chosen
$t_{Si}$ value the vertical second-order mode is indeed guided (see Fig.
2(a)), we can consider single-mode operation, as its propagation losses are
greater than 50 dB/cm. Higher $t_{Si}$ values would barely benefit the
fundamental mode loss, but undesirably decrease the second-order mode leakage.

The waveguide core width and the SWG cladding geometry are chosen to avoid
lateral leakage and to enhance the mechanical stability. Note that a scaled
version of the original design in [20], i.e. $W_{core}$ = 2.6 $\mu$m and
$W_{clad}$ = 5 $\mu$m, would compromise waveguide stability because of the
thicker Si thickness of 1.4 $\mu$m. Therefore, we made the waveguide core
wider ($W_{core}$ = 2.9 $\mu$m), increasing the lateral mode confinement and
allowing a reduction in the cladding width up to $W_{clad}$ = 3 $\mu$m, with
the intention of reducing the torque without degrading the lateral leakage. It
is also important to mention that, even though all dimensions are
substantially enlarged to migrate to $\lambda$ = 7.67 $\mu$m, the cladding
width hardly increases with respect to the design at $\lambda$ = 3.8 $\mu$m,
resulting in a remarkably more stable waveguide. At $\lambda$ = 7.67 $\mu$m,
the Bragg period is moved away up to $\Lambda_{Bragg}\approx$ 1300 nm for a
wide waveguide, thereby extending the subwavelength regime in relation to the
design in [20]. Due to the higher $t_{Si}$ value, larger holes are needed to
assure an effective dry etching of the cladding and let HF penetrate through
them. Therefore, we selected $L_{hole}$ = 900 nm. Finally, the fundamental
mode lateral confinement and the waveguide mechanical stability were slightly
improved by choosing $L_{Si}$ = 250 nm.

We simulated 90°-bends. Negligible losses are achieved for bend radius of 35
$\mu$m. S-bends were also designed, with propagation losses smaller than 0.04
dB/bend for an offset of 5 $\mu$m and a length of 75 $\mu$m.

To couple light from the optical fiber to the chip, we designed, for the first
time, a focusing suspended surface grating coupler at a nominal wavelength of
$\lambda$ = 7.67 $\mu$m. It is composed of a silicon membrane of thickness
$t_{Si}$ where gaps are etched in the area where light should be coupled. The
gaps again have a double purpose: they allow the removal of the silicon
dioxide and constitute a periodic radiating waveguide formed by the succession
of silicon and air wide strips. We achieve a simulated coupling efficiency of
58% at the central wavelength and a 1-dB bandwidth of 230 nm. Very low back-
reflections around 0.1% are also achieved by including a thinner air gap at
the beginning of the grating. The radiation angle is $\sim$19.3°. More details
about the theoretical design process of this grating coupler and a full
experimental characterization will be given in a future work.

Fabrication was carried out on a 6-inch SOI wafer. We implemented the same
fabrication process as described in [20]: e-beam lithography was used for
patterning with ZEP-520A as positive resist, and etching was accomplished with
an ICP tool using SF6/C4F8 chemistry. The wafer was then submerged into a 1:7
HF bath for 2.5 hours, to ensure removal of the BOX beneath the widest
devices. The resist was only removed after the HF etching in order to protect
the Si top surface, as it was observed that this long wet etch produced a
significant roughening of the Si surfaces. In our previous work [20] it was
observed that the dry and wet etches both resulted in some lateral etching,
reducing the width of the silicon supports and the waveguide core thickness.
To compensate for this a positive biasing of 220 nm (i.e. increase of the
dimensions) was applied to the designed layout of all the unetched features.

Figure 3(a) shows a scanning electron microscope (SEM) image of a fabricated
waveguide, and Fig. 3(b) shows a fabricated grating coupler. It can be seen
that both devices are fully suspended and that there is no significant bending
of the free standing Si. However, it can be seen that in some places roughness
has been created along the edges of the etched structures in places where the
HF was able to reach between the Si and photoresist layers during wet etching.
Atomic force microscope measurements show that in these places 40 nm deep
grooves have been etched into the Si. The final thickness of the Si waveguide,
as measured using the SEM, was approximately 1420 nm.

Figure 3: a) SEM image of a cleaved waveguide facet after fabrication. The BOX
can be seen to have been completely removed below the waveguide, but the HF
etching has created roughness along the etched edges. b) SEM image of the
focusing grating coupler and taper with the BOX removed. It has been cleaved
perpendicular to the waveguide, part way through the grating. The grating was
the widest structure on the chip at 40 $\mu$m width.

The experimental setup employed a single mode continuous wave distributed
feedback quantum cascade laser (Thorlabs QD7500CM1) with 106 mW maximum output
power at a wavelength of 7.67 $\mu$m. A black diamond-2 lens with a focal
length of 1.9 mm was used to collimate the light from the laser, a chopper
wheel modulated the light and another black diamond-2 lens with a focal length
of 6.0 mm was used to couple it into a single mode As2Se3 fiber (Coractive
IRT-SE-28/170) with a measured output of $\sim$ 3.4 mW. Light was coupled into
the waveguides with grating couplers (Fig. 3(b)) and the output was collected
via another single mode fiber and coupled to a liquid nitrogen cooled HgCdTe
detector (Infrared Associates Inc. MCT-13-1.00). The signal from the detector
was passed through a pre-amplifier before being sent to a lock-in amplifier to
improve the signal to noise ratio. In order to characterize the dynamic range
of the setup we placed a gold mirror in the position of a chip, finding that
the maximum transmission was 60 dB above the noise floor.

The waveguide propagation loss was measured with the effective "cut-back"
method. Figures 4(a) and 4(b) show the transmissions of waveguides of
different lengths in samples before and after BOX removal, respectively. The
figures show that removing the SiO2 reduces the waveguide loss from $62.3\pm
9.6$ dB/cm to $3.1\pm 0.3$ dB/cm. According to simulations 2.2 dB/cm of the
suspended waveguide loss can be attributed to material loss ($\sim$2.1 dB/cm)
and substrate leakage ($\sim$0.1 dB/cm), therefore approximately 0.9 dB/cm of
the loss is likely due to scattering, predominantly from the roughness created
by HF etching. This waveguide loss is comparable to other group-IV material
waveguide demonstrations in this wavelength range [11, 16, 15], and we expect
that it can be reduced further through optimization of the HF etching process.

The 90° bend loss was 0.08$\pm$0.02 dB/bend and the loss of the s-bends was
0.06$\pm$0.02 dB/bend.

Figure 4: Waveguide cut-back loss measurement results at 7.67 $\mu$m
wavelength for (a) Si waveguides with the BOX still present, and (b) suspended
Si waveguides. For each set of points the transmissions have been normalized
to the transmission through the shortest waveguide (0.02 cm) in each set. The
red lines are linear fits to the measured points (cross symbols).

We have designed, fabricated and characterized suspended silicon waveguides,
90o bends and s-bends at a wavelength of 7.67 $\mu$m. The propagation loss for
the waveguides was $3.1\pm 0.3$ dB/cm while the bend loss was 0.08$\pm$0.02
dB/bend and 0.06$\pm$0.02 dB/bend for the 90o bends and s-bends respectively.
The material loss of Si contributes by $\sim$2.1 dB/cm to the propagation loss
figure, and it is expected that the remainder of the loss predominantly comes
from scattering at the roughness that prolonged HF etching creates at the
Si/air interfaces. Our future work will focus on refining the fabrication
process in order to reduce roughness and to approach propagation loss values
closer to the material loss. These results demonstrate that suspended silicon
waveguides designed using the sub-wavelength grating concept and fabricated
from SOI wafers can be used throughout the mid-infrared transparency window of
silicon, and unlock potential sensing applications in this range.

Acknowledgements. The data used to create Figs. 2 and 4 can be accessed online
at http://doi.org/10.5258/SOTON/D0322

Funding. Engineering and Physical Sciences Research Council (EPSRC) (MIGRATION
EP/L01162X/1, National Hub in High Value Photonic Manufacturing EP/N00762X/1,
CORNERSTONE EP/L021129/1); Royal Academy of Engineering (RAEng) (M.
Nedeljkovic fellowship RF201617/16/33); Ministerio de Economía y
Competitividad (MINECO) (TEC2016-80718-R); Ministerio de Educación, Cultura y
Deporte (MECD) (FPU14/06121); Fondo Europeo de Desarrollo Regional – FEDER;
Universidad de Málaga.

## References

  * [1] R. A. Soref, S. J. Emelett, and W. R. Buchwald, Journal of Optics A: Pure and Applied Optics 8, 840 (2006).
  * [2] “NIST standard reference data program,” http://webbook.nist.gov/cgi/cbook.cgi?ID=C7446095&Type=IR-SPEC&Index=0#.
  * [3] R. Kitamura, L. Pilon, and M. Jonasz, Appl. Opt. 46, 8118 (2007).
  * [4] T. Baehr-Jones, A. Spott, R. Ilic, A. Spott, B. Penkov, W. Asher, and M. Hochberg, Opt. Express 18, 12127 (2010).
  * [5] Y. Liu, A. Spott, T. Baehr-Jones, R. Ilic, and M. Hochberg, “Silicon waveguides and ring resonators at 5.5 $\mu$m,” in “7th IEEE International Conference on Group IV Photonics,” (2010), pp. 377–379.
  * [6] F. Li, S. D. Jackson, C. Grillet, E. Magi, D. Hudson, S. J. Madden, Y. Moghe, C. O’Brien, A. Read, S. G. Duvall, P. Atanackovic, B. J. Eggleton, and D. J. Moss, Opt. Express 19, 15212 (2011).
  * [7] N. Singh, D. D. Hudson, and B. J. Eggleton, Opt. Express 23, 17345 (2015).
  * [8] P. T. Lin, V. Singh, L. Kimerling, and A. M. Agarwal, Applied Physics Letters 102, 251121 (2013).
  * [9] W. Li, P. Anantha, S. Bao, K. H. Lee, X. Guo, T. Hu, L. Zhang, H. Wang, R. Soref, and C. S. Tan, Applied Physics Letters 109, 241101 (2016).
  * [10] J. M. Ramirez, V. Vakarin, C. Gilles, J. Frigerio, A. Ballabio, P. Chaisakul, X. L. Roux, C. Alonso-Ramos, G. Maisons, L. Vivien, M. Carras, G. Isella, and D. Marris-Morini, Opt. Lett. 42, 105 (2017).
  * [11] M. Brun, P. Labeye, G. Grand, J.-M. Hartmann, F. Boulila, M. Carras, and S. Nicoletti, Opt. Express 22, 508 (2014).
  * [12] M. Nedeljkovic, J. S. Penadés, C. J. Mitchell, A. Z. Khokhar, S. Stanković, T. D. Bucio, C. G. Littlejohns, F. Y. Gardes, and G. Z. Mashanovich, IEEE Photonics Technology Letters 27, 1040 (2015).
  * [13] A. Malik, M. Muneeb, S. Pathak, Y. Shimura, J. V. Campenhout, R. Loo, and G. Roelkens, IEEE Photonics Technology Letters 25, 1805 (2013).
  * [14] Y.-C. Chang, V. Paeder, L. Hvozdara, J.-M. Hartmann, and H. P. Herzig, Opt. Lett. 37, 2883 (2012).
  * [15] V. Vakarin, J. M. Ramírez, J. Frigerio, A. Ballabio, X. L. Roux, Q. Liu, D. Bouville, L. Vivien, G. Isella, and D. Marris-Morini, Opt. Lett. 42, 3482 (2017).
  * [16] M. Nedeljkovic, J. S. Penades, V. Mittal, G. S. Murugan, A. Z. Khokhar, C. Littlejohns, L. G. Carpenter, C. B. E. Gawith, J. S. Wilkinson, and G. Z. Mashanovich, Opt. Express 25, 27431 (2017).
  * [17] Y. Xia, C. Qiu, X. Zhang, W. Gao, J. Shu, and Q. Xu, Opt. Lett. 38, 1122 (2013).
  * [18] Z. Cheng, X. Chen, C. Y. Wong, K. Xu, and H. K. Tsang, IEEE Photonics Journal 4, 1510 (2012).
  * [19] T.-H. Xiao, Z. Zhao, W. Zhou, M. Takenaka, H. K. Tsang, Z. Cheng, and K. Goda, Opt. Lett. 42, 2882 (2017).
  * [20] J. S. Penadés, A. Ortega-Moñux, M. Nedeljkovic, J. G. Wangüemert-Pérez, R. Halir, A. Z. Khokhar, C. Alonso-Ramos, Z. Qu, I. Molina-Fernández, P. Cheben, and G. Z. Mashanovich, Opt. Express 24, 22908 (2016).
  * [21] J. S. Penadés, C. Alonso-Ramos, A. Z. Khokhar, M. Nedeljkovic, L. A. Boodhoo, A. Ortega-Moñux, I. Molina-Fernández, P. Cheben, and G. Z. Mashanovich, Opt. Lett. 39, 5661 (2014).
  * [22] P. Cheben, D.-X. Xu, S. Janz, and A. Densmore, Opt. Express 14, 4695 (2006).
  * [23] P. Cheben, P. J. Bock, J. H. Schmid, J. Lapointe, S. Janz, D.-X. Xu, A. Densmore, A. Delâge, B. Lamontagne, and T. J. Hall, Opt. Lett. 35, 2526 (2010).
  * [24] R. Halir, P. J. Bock, P. Cheben, A. Ortega-Moñux, C. Alonso-Ramos, J. H. Schmid, J. Lapointe, D.-X. Xu, J. G. Wangüemert-Pérez, I. Molina-Fernández, and S. Janz, Laser & Photonics Reviews 9, 25 (2015).
  * [25] W. Zhou, Z. Cheng, X. Wu, B. Zhu, X. Sun, and H. K. Tsang, Opt. Lett. 42, 1245 (2017).
  * [26] D. Chandler-Horowitz and P. M. Amirtharaj, Journal of Applied Physics 97, 123526 (2005).
  * [27] R. Halir, P. Cheben, S. Janz, D.-X. Xu, I. Molina-Fernández, and J. G. Wangüemert-Pérez, Opt. Lett. 34, 1408 (2009).
  * [28] L. Zavargo-Peche, A. O.-M. nux, J. G. Wangüemert-Pérez, and I. Molina-Fernández, Prog. Electromagn. Res. 123, 447 (2012).

sample