[subfigure]justification=centering area density Swap Opt opt lru Department of
Computer Science
University of Maryland, College Park, MD, USA<EMAIL_ADDRESS>https://orcid.org/0000-0003-0868-9802 A. Flores -Velazco [500]Theory of
computation Computational geometry 1 2022

# Improved Search of Relevant Points for Nearest-Neighbor Classification

Alejandro Flores-Velazco

###### Abstract

Given a training set $P\subset\mathbb{R}^{d}$, the _nearest-neighbor
classifier_ assigns any query point $q\in\mathbb{R}^{d}$ to the class of its
closest point in $P$. To answer these classification queries, some training
points are more relevant than others. We say a training point is _relevant_ if
its omission from the training set could induce the misclassification of some
query point in $\mathbb{R}^{d}$. These relevant points are commonly known as
_border points_ , as they define the boundaries of the Voronoi diagram of $P$
that separate points of different classes. Being able to compute this set of
points efficiently is crucial to reduce the size of the training set without
affecting the accuracy of the nearest-neighbor classifier.

Improving over a decades-long result by [8], in a recent paper [12] an
_output-sensitive_ algorithm was proposed to find the set of border points of
$P$ in $\mathcal{O}(n^{2}+nk^{2})$ time, where $k$ is the size of such set. In
this paper, we improve this algorithm to have time complexity equal to
$\mathcal{O}(nk^{2})$ by proving that the first steps of their algorithm,
which require $\mathcal{O}(n^{2})$ time, are unnecessary.

###### keywords:

nearest-neighbor classification, nearest-neighbor rule, decision boundaries,
border points, relevant points

###### category:

iven a training set P ⊂ ℝ^d, the nearest-neighbor classifier assigns any query
point q ∈ ℝ^d to the class of its closest point in P. To answer these
classification queries, some training points are more relevant than others. We
say a training point is relevant if its omission from the training set could
induce the misclassification of some query point in ℝ^d. These relevant points
are commonly known as border points, and they define the boundaries of P
between points of different classes. Being able to compute this set of points
efficiently is crucial to reduce the size of the training set without
affecting the accuracy of the nearest-neighbor classifier.

Improving over a decades-long result by Clarkson, in a recent paper by
Eppstein an output-sensitive algorithm was proposed to find the set of border
points of P in O(n^2 + nk^2) time, where k is the size of such set. In this
paper, we further improve this algorithm to have time complexity equal to
O(nk^2) by proving that the first steps of their algorithm, which require
O(n^2) time, are in fact unnecessary.

##  1 Introduction

In the context of non-parametric classification, we are given a training set
$P\subset\mathbb{R}^{d}$ consisting of $n$ _labeled_ points in $d$-dimensional
Euclidean space, where the label of every point in $P$ indicates the _class_
(or _color_) that the point belongs to. The goal of a classifier is to use the
training set $P$ to _predict_ the class for any _unlabeled_ query point
$q\in\mathbb{R}^{d}$, that is, to _classify_ $q$.

The _nearest-neighbor classifier_ (also known as _nearest-neighbor rule_) [13]
stands out as a simple yet powerful method, that works by assigning any query
point $q$ to the class of its closest point in $P$. Despite its simplicity,
the nearest-neighbor classifier is well-known to exhibit good classification
accuracy both experimentally and theoretically [30, 10, 11]. In fact, it is
still frequently used in many applications [5, 22, 18, 25, 27, 29, 26] over
more recent and sophisticated techniques like support-vector machines [9] and
deep neural networks [28].

One of the principal disadvantages of this technique is its high dependency on
the size and dimensionality of the data, especially in light of _big data_
applications. With training sets with billions of points becoming increasingly
common, reducing the nearest-neighbor classifier’s dependency on $n$ and $d$
is one approach to enhance its efficiency. There has been significant progress
towards this goal, mainly focusing on two directions. The first involves the
design of efficient data structures to answer approximate nearest-neighbor
queries [17, 23, 19, 20, 4, 3, 2]. The second direction focuses on reducing
the size of the training set used by the nearest-neighbor classifier, thus
effectively reducing $n$. However, most practical techniques for training set
reduction provide limited guarantees on the effect of this reduction to the
accuracy of the nearest-neighbor classifier [15, 14, 1, 16].

[b].35  [b].35

Figure 1: Original training set $P$ Figure 2: Border/Relevant points of $P$
Figure 3: On the left, a training set $P$ with points of three classes: _red_
, _blue_ and _yellow_. There the black lines highlight the _boundaries_ of $P$
between points of different classes. On the right, a subset of these points
corresponding to the set of border points of $P$. Note that by definition, the
boundaries between points of different classes remain the same for $P$ and for
its set of border points.

Only a handful of works [8, 6, 12] have proposed training set reduction
algorithms that guarantee the same classification of every query point, before
and after the reduction took place. These are called _boundary preserving_
algorithms, and it is the focus of this paper.

The set of _border points_ (or _relevant points_ 111While [12] uses the term
_relevant points_ , the term _border points_ has been the standard in the
literature of this and other related problems [15, 14, 21, 32]. For this
reason, we stick to the term _border points_.) of the training set $P$ are
those that define the boundaries between points of different classes, and
whose omission from the training set would imply the misclassification of some
query points in $\mathbb{R}^{d}$. Formally, two points $p,\hat{p}\in P$ are
border points of $P$ if they belong to different classes, and there exist some
point $q\in\mathbb{R}^{d}$ such that $q$ is equidistant to both $p$ and
$\hat{p}$, and no other point of $P$ is closer to $q$ than these two points
(_i.e._ , the empty ball property of Voronoi Diagrams). See Figure 3 for an
example of a training set $P$ in $\mathbb{R}^{2}$ and its set of border
points. Throughout, we let $k$ denote the total number of border points in the
training set. By definition, if instead of building the nearest-neighbor
classifier with the entire training set $P$ we use the set of border points of
$P$, its dependency is reduced from $n$ to $k$, while still obtaining the same
classification for any query point in $\mathbb{R}^{d}$. This becomes
particularly relevant for applications where $k\ll n$.

In this paper, we improve a recently proposed algorithm by [12] that computes
the set of border points of any training set $P\subset\mathbb{R}^{d}$, where
dimension $d$ is assumed to be constant. While the original algorithm computes
such set in $\mathcal{O}(n^{2}+nk^{2})$ time, where $k$ is the number of
border points of $P$, our new algorithm computes the same set in
$\mathcal{O}(nk^{2})$ time.

###  1.1 Previous Work

Other related problems in the realm of training set reduction are NP-hard [31,
24, 33] to solve exactly (_e.g._ , those of finding minimum cardinality
_consistent_ subsets and _selective_ subsets). However, the problem of
preserving the class boundaries of the nearest-neighbor classifier, or simply,
finding the set of border points of $P$, is tractable.

For training sets $P\subset\mathbb{R}^{2}$ in 2-dimensional Euclidean space,
Bremner et al. [6] proposed an output-sensitive algorithm for finding the set
of border points of $P$ in $\mathcal{O}(n\log{k})$ worst-case time. However,
how to generalize this algorithm for higher dimensions remained unclear.

Until very recently, the best result for the higher dimensional case was that
of Clarkson [8]. He proposed an algorithm to find the set of border points of
$P\subset\mathbb{R}^{d}$, with bounded $d$, that runs in
$\mathcal{O}(\min{(n^{3},kn^{2}\log{n})})$ worst-case time. For almost three
decades, this remained the best result for training sets in $\mathbb{R}^{d}$.
Recently, Eppstein [12] proposed a significantly faster algorithm for the
$d$-dimensional Euclidean case, which runs in $\mathcal{O}(n^{2}+nk^{2})$
worst-case time.

Eppstein’s algorithm is strikingly simple, yet full of interesting ideas (see
Algorithm 1). The algorithm works as follows: it begins by selecting an
initial set of border points of $P$, one point from every class region. From
here, the algorithm uses a series of subroutines which we will group together
and denote as the “ _inversion method_ ”, to find the remaining border points
of $P$. Thus, the algorithm can be naturally split into two steps: the
initialization of $R$ with some border points, and the search process for the
remaining border points of $P$.

Input: Initial training set $P$

Output: The set of border points of $P$

1 Let $M$ be the MST of $P$

2 Initialize $R$ with the end points of every bichromatic edge of $M$

3 foreach  _$p\in R$_ do

4 Let $c$ be $p$’s class and $P_{c}$ be the points of $P$ that belong to class
$c$

5 Let $S_{p}$ be the inverted points of $P\setminus P_{c}$ around $p$

6 Find all extreme points of $S_{p}$ and their corresponding original points
$E_{p}$

7 $R\leftarrow R\cup E_{p}$

8

return _$R$_

Algorithm 1 Recent Eppstein’s algorithm [12] to find the set of border points
of $P$.

The initialization step (lines 1–2 of Algorithm 1) involves finding a subset
of border points such that at least one point for every class region is
selected. Eppstein observes that this can be achieved by computing the Minimum
Spanning Tree (MST) of $P$, identifying the edges of the MST that connect
points of different classes (denoted as _bichromatic_ edges), and selecting
the endpoints of all such edges. This step takes $\mathcal{O}(n^{2})$ time,
but we will prove that it is not necessary.

The search step (lines 3–6 of Algorithm 1) is in charge of finding every
remaining border point of $P$. This step iterates over all selected points,
and for each such point $p$, it performs what we call the inversion method.
This method identifies a subset of border points of $P$, which are added to
$R$. Once the algorithm has done the inversion method on every point of $R$,
it terminates with the guarantee of having selected every border point of $P$.

Given any point $p\in P$, the inversion method on $p$ is described in lines
4–6 of Algorithm 1. Let $c$ be $p$’s class, and $P_{c}$ be the points of $P$
that belong to class $c$, the inversion method on $p$ consists of:  (i)
inverting all points of $P\setminus P_{c}$ around a ball centered at $p$ (call
the set of these inverted points as $S_{p}$ and include $p$ itself in the
set),  (ii) computing the set of extreme points of $S_{p}$, and finally  (iii)
returning the set $E_{p}$ of those points of $P$ that correspond to the
extreme points of $S_{p}$ before inversion.  For a detailed description and
proof of correctness of this method, we refer the reader to Eppstein’s paper
[12]. However, for the purposes of this paper we only need a property
presented in Lemma 3 of [12]: the points in $E_{p}$ reported by the inversion
method are the Delaunay neighbors of $p$ with respect to the set $(P\setminus
P_{c})\cup\\{p\\}$.

Every call of the inversion method takes $\mathcal{O}(nk)$ time by leveraging
well-known output-sensitive algorithms for computing extreme points. Given
that this method is called exclusively on every border point of the training
set, this yields a total of $\mathcal{O}(nk^{2})$ time to complete the search
step of the algorithm. Overall, this implies that Eppstein’s algorithm
computes the entire set of border points of $P$ in $\mathcal{O}(n^{2}+nk^{2})$
worst-case time.

##  2 Our Approach

We propose a simple modification to Eppstein’s algorithm, which avoids the
step of computing the MST of the training set $P$, along with the subsequent
selection of bichromatic edges to produce the initial subset of border points.

Instead, we simply start the search process with any arbitrary point of $P$.
The rest of the algorithm remains virtually unchanged (see Algorithm 2 for a
formal description). We show that this new approach is not only correct,
meaning that it only finds border points of $P$, but also complete, as all
border points of $P$ are eventually found by our algorithm. Additionally, by
avoiding the main bottleneck of the original algorithm, our new algorithm
computes the same result in $\mathcal{O}(nk^{2})$ time, eliminating the
$\mathcal{O}(n^{2})$ term.

Input: Initial training set $P$

Output: The set of border points of $P$

1 Let $s$ be any “ _seed_ ” point from $P$

2 $R\leftarrow\phi$

3 foreach  _$p\in R\cup\\{s\\}$_ do

4 Let $c$ be $p$’s class and $P_{c}$ be the points of $P$ that belong to class
$c$

5 Let $S_{p}$ be the inverted points of $P\setminus P_{c}$ around $p$

6 Find all extreme points of $S_{p}$ and their corresponding original points
$E_{p}$

7 $R\leftarrow R\cup E_{p}$

8

return _$R$_

Algorithm 2 New algorithm to find the set of border points of $P$.

Before proceeding, it is useful to explore why Eppstein’s algorithm computes
the MST of the training set $P$. First, note that the original algorithm only
applies the inversion method on border points of $P$. In fact, Eppstein’s
correctness proof relies on it: Lemma 6 in [12] proves that all points in
$E_{p}$ are border points by assuming that point $p$ is also a border point.
From the description of our algorithm, note that we initially apply the
inversion method on a “seed” point $s$, which might not be a border point.
Therefore, we need to generalize Lemma 6 in [12] for the case where $p$ is not
a border point of $P$. Additionally, using the points from all bichromatic
pairs of the MST of $P$ guarantees that Eppstein’s algorithm starts the search
step with at least one point from every boundary of $P$. Eppstein’s
completeness proof shows that the search step can then “move along” any given
boundary and eventually select all its defining points. We show that the
search process is far more powerful, and can even “jump” between nearby
boundaries, thus rendering the MST computation unnecessary.

The following description outlines the necessary steps to prove both the
correctness and completeness of our new algorithm, which are unfolded in the
rest of this section.

  * •

By applying the inversion method to any point of $P$, not necessarily a border
point, all reported points are border points of $P$. This is established in
Lemma 2.1, generalizing the statement of Lemma 6 of [12] for non-border
points.

  * •

For any class boundary of $P$, once the algorithm selects a point from this
boundary, it will eventually select every other point defining the same
boundary. This is originally proved in Lemma 10 [12], however, we provide
simpler proofs in Lemmas 2.3 and 2.5.

  * •

Given two disconnected boundaries separated by a class region, we prove that
if our algorithm selects a defining point from one of the boundaries, it will
eventually select all defining points from both boundaries. This is proved in
Lemma 2.6.

All together, these lemmas are used to prove the main result: the correctness,
completeness, and worst-case time complexity of Algorithm 2, as stated in
Theorems 2.8 and 2.10.

[b].32  [b].32  [b].32

Figure 4: Training set $P$ Figure 5: Points from $(P\setminus
P_{c})\cap\\{p\\}$ Figure 6: $q$ and $\hat{p}$ are border points Figure 7:
Example showing the inversion method from any point $p\in P$. On the left,
training set $P$. The middle figure shows every non red point of $P$, except
for $p$ itself, along with a point $q$ selected from the inversion method on
$p$. On the right, we see evidence that $q$ is a border point of $P$.

######  Lemma 2.1.

Let $p\in P$ be any point of the training set. Then every point selected using
the inversion method on $p$ must be a border point of $P$.

######  Proof 2.2.

Let $E_{p}$ be the points of $P$ corresponding (before inversion) to the
extreme points of $S_{p}$. According to Lemma 3 [12], every point in $E_{p}$
is a neighbor of point $p$ with respect to the Voronoi Diagram of set
$(P\setminus P_{c})\cup\\{p\\}$. This implies that for every point $q\in
E_{p}$ other than p, there exists a ball such that both $p$ and $q$ are on its
surface and no points of $P\setminus P_{c}$ lie inside (see Figures 7 and 7).
We can now leverage similar techniques to the ones described in [6], to find a
“witness” point to the hypothesis that $q$ must be a border point of $P$.

Recall that the empty ball we just described, as illustrated in Figure 7, is
empty from points of $P\setminus P_{c}$. However, there might be points of
$P_{c}$ inside. And moreover, we know that at least one point of $P_{c}$,
point $p$, lies on its surface. Now, let $r$ be the center of this ball, we
grow an empty ball, this time with respect to the entire training set $P$,
such that its center lies on the line $\overline{qr}$ and point $q$ is on its
surface (see Figure 7). This ball will grow until it hits another point
$\hat{p}$ of $P$, which we are guaranteed it will be of the same class as
point $p$, and thus, of different class as point $q$. Finally, we have just
found an empty ball with respect to $P$, which has points $q$ and $\hat{p}$ on
its surface, and were the class of both points differ. Therefore, this implies
that $q$ is a border point of $P$.

[b].32  [b].32

Figure 8:  Figure 9:  Figure 10: By definition, any two adjacent walls $w_{1}$
and $w_{2}$ of the Voronoi Diagram of $P$ hold the empty ball property with
the points that define them. When these walls are part of the class boundaries
of $P$, the points that define them belong to at least two classes.

Before continuing, we need to formally define a few concepts. First, we define
a _wall_ of $P$ as any $(d-1)$-dimensional face of the Voronoi Diagram of $P$.
By known properties of these structures, every wall $w$ is _defined_ by two
distinct points $p,q\in P$ such that any point on $w$ has $p$ and $q$ as its
two equidistant nearest-neighbors in the training set. We say two walls are
_adjacent_ if their intersection is not empty. That is, if there exists a
point in $\mathbb{R}^{d}$ with all the defining points of these two walls as
its equidistant nearest-neighbors in $P$.

Additionally, we define a _class boundary_ (or just _boundary_) of $P$ as the
union of adjacent walls, where each of these walls is defined by two points of
different classes. Similarly, we define a _class region_ of $P$ as the union
of adjacent Voronoi cells whose defining points belong to the same class.
Based on these definitions, note that class boundaries are the ones that
separate different class regions of $P$. Figure 11 illustrates a training set
in $\mathbb{R}^{2}$ with points of three classes, whose Voronoi Diagram
describes five class regions and two class boundaries.

######  Lemma 2.3.

Let $w_{1}$ and $w_{2}$ be two adjacent walls in a class boundary of $P$. If
the algorithm selects one of the points defining one of these walls, it
eventually selects the remaining points defining both walls.

######  Proof 2.4.

Let $\mathcal{W}$ be the set of points defining both walls $w_{1}$ and $w_{2}$
(see Figure 10). By definition, these two walls of the Voronoi Diagram of $P$
are adjacent if there exists an empty ball with all the points of
$\mathcal{W}$ on its surface. Knowing these two walls are part of the class
boundaries of $P$, the set $\mathcal{W}$ must contain at least three points,
and at least two classes.

Let $p_{1}$ be the first point of $\mathcal{W}$ to be selected by the
algorithm. When doing the inversion method on point $p_{1}$, the algorithm
will select all points of $\mathcal{W}$ of different class than $p_{1}$, of
which we know there is at least one. Let $p_{2}$ be one such point. Finally,
when doing the inversion method on point $p_{2}$, the algorithm will select
the remaining points of $\mathcal{W}$ of the same class as $p_{1}$. Therefore,
all points of $\mathcal{W}$ will eventually be selected by the algorithm.

Figure 11: A training set with five class regions (one _blue_ , two _red_ ,
and two _yellow_ regions), along with two disconnected class boundaries that
separate all these regions. On the left, a boundary separating the blue region
and the leftmost yellow and red regions. On the right, a boundary that
separates the same blue region and the remaining red and yellow regions.

######  Lemma 2.5.

Let $\mathcal{A}$ be a class boundary of $P$, and assume that the algorithm
selects one of the defining points of $\mathcal{A}$. Then, the algorithm will
eventually select all defining points of $\mathcal{A}$.

This comes as a direct consequence of Lemma 2.3 and the definition of a class
boundary of the training set $P$. It remains to show what happens with
boundaries that are disconnected.

######  Lemma 2.6.

Let $\mathcal{A}$ and $\mathcal{B}$ be two disconnected boundaries of $P$,
such that there exists a path in space from $\mathcal{A}$ to $\mathcal{B}$
that is completely contained within one color region. Without loss of
generality, say that every point that defines $\mathcal{A}$ has been selected
by the algorithm. Then, every point that defines $\mathcal{B}$ must also be
selected by the algorithm.

######  Proof 2.7.

Given these two disconnected boundaries $\mathcal{A}$ and $\mathcal{B}$, we
assume there exists some path $\mathcal{P}$ in $\mathbb{R}^{d}$ going from a
wall of $\mathcal{A}$ to a wall of $\mathcal{B}$, such that this path passes
exclusively through a single class region (see Figure 15). Without loss of
generality, say this is a _red_ class region. Formally, for every point $r$
along $\mathcal{P}$ we know $r$’s nearest-neighbor in $P$ is red.
Additionally, we assume that every border point defining $\mathcal{A}$ is
selected by the algorithm. Hence, the proof consists of showing that there
exists a sequence of border points $\langle
p_{1},\hat{p}_{1},p_{2},\hat{p}_{2},\dots,p_{m},\hat{p}_{m}\rangle$ such that
(i) $p_{1}$and $\hat{p}_{m}$ are defining points of $\mathcal{A}$ and
$\mathcal{B}$, respectively,  (ii) $\hat{p}_{i}$is retrieved by the inversion
method on $p_{i}$, for every $i\in[1,m]$, and finally  (iii) points $p_{i}$
and $\hat{p}_{i-1}$ are both defining the same boundary, for every
$i\in[2,m]$.  See Figure 15 for a visual description.

[b].32  [b].32  [b].32

Figure 12:  Figure 13:  Figure 14:  Figure 15: On the right, two disconnected
boundaries $\mathcal{A}$ and $\mathcal{B}$ enclosing a red class region. Thus,
there is a path $\mathcal{P}$ completely contained inside such region and
connecting both boundaries. Other boundaries can also be enclosing the same
region and be near path $\mathcal{P}$. On the left, we proof that there exists
a sequence of points that can be retrieved by calls to the inversion method,
such that if points of $\mathcal{A}$ are selected by the algorithm, eventually
points of $\mathcal{B}$ will also be selected.

By definition, for every point $r$ along path $\mathcal{P}$ we know $r$’s
nearest-neighbor is a red point. Now, let’s delete every red point from
consideration, including the ones defining boundaries $\mathcal{A}$ and
$\mathcal{B}$ (see Figure 15). This immediately implies that $r$’s nearest-
neighbor just became a non-red border point of $P$. The fact that $r$’s new
nearest-neighbor is a border point is easy to proof, using similar arguments
as the ones laid down in Lemma 2.1. Additionally, these border points could be
defining other boundaries apart from $\mathcal{A}$ and $\mathcal{B}$, as seen
in Figure 15.

Let’s start moving along the path $\mathcal{P}$, starting from the end-point
of the path that lies on a wall of boundary $\mathcal{A}$. Then, find all
$r_{i}$ points along the path, where each $r_{i}$ has two equidistant nearest-
neighbors among the remaining non-red points, and both points define two
distinct boundaries of $P$. We say there are $m$ of these points along the
path, and denote $r_{i}$’s two equidistant nearest-neighbors as $q_{i,1}$ and
$q_{i,2}$ for $i\in[1,m]$. Clearly, $q_{i,1}$ and $q_{i-1,2}$ are border
points defining the same boundary, for all $i\in[2,m]$. See Figure 15, where
the three black points along the path are the $r_{i}$ points, and the _yellow_
and _blue_ points on the surface of the balls centered at each $r_{i}$ are the
corresponding $q_{i,1}$ and $q_{i,2}$ points.

For now, let’s fix the analysis on one such $r_{i}$ point, and consider the
ball centered at $r_{i}$ with both $q_{i,1}$ and $q_{i,2}$ on its surface.
There must exist some other point $q_{i,3}$ lying inside of $r_{i}$’s ball,
such that $q_{i,3}$ is one of the deleted red points defining the same
boundary as $q_{i,1}$. It is now easy to see that there exist an empty ball,
with respect to the set $P\setminus P_{\emph{red}}\cup\\{q_{i,3}\\}$, with
both $q_{i,3}$ and $q_{i,2}$ on its boundary. This implies that $q_{i,2}$ is
retrieved by the inversion method on $q_{i,3}$. Therefore, let’s add
$p_{i}\leftarrow q_{i,3}$ and $\hat{p}_{i}\leftarrow q_{i,2}$ to the sequence
of points that we are looking for. Repeat this for every $r_{i}$ with
$i\in[1,m]$ to identify all points in the sequence.

Finally, we have the sequence of border points $\langle
p_{1},\hat{p}_{1},p_{2},\hat{p}_{2},\dots,p_{m},\hat{p}_{m}\rangle$ such that
for any $i\in[1,m]$ assuming that the algorithm selects the points defining
the same boundary as $p_{i}$, it will also select $\hat{p}_{i}$, and
leveraging Lemma 2.5 it will eventually select all other points defining the
same boundary as $\hat{p}_{i}$. Given that $p_{1}$ and $\hat{p}_{m}$ are
defining border points of boundaries $\mathcal{A}$ and $\mathcal{B}$,
respectively, and by the assumption that all points defining $\mathcal{A}$ are
selected by the algorithm, we know that eventually, all points defining
$\mathcal{B}$ will be selected too.

######  Theorem 2.8.

The algorithm selects every border point of $P$ in $\mathcal{O}(nk^{2})$ time.

######  Proof 2.9.

Proving the worst-case time complexity of our algorithm follows directly from
the time complexity of the search step of Eppstein’s algorithm [12]. However,
the correctness and completeness of our algorithm follows from Lemmas 2.1 to
2.6.

First, we know by Lemmas 2.1-2.5 that Algorithm 2 will select the defining
border points of at least one class boundary of $P$. Denote this boundary as
$\mathcal{A}$ and consider any other boundary $\mathcal{B}$ of $P$. Evidently,
we can draw a path $\mathcal{P}$ from $\mathcal{A}$ to $\mathcal{B}$, which
would generally pass through several class regions. Then, let’s split
$\mathcal{P}$ into several subpaths
$\mathcal{P}_{1},\mathcal{P}_{2},\dots,\mathcal{P}_{m}$ such that each subpath
is completely contained within a single class region. From this, we can
directly apply Lemma 2.6 on each of the intermediate boundaries that “cut”
$\mathcal{P}$ into these subpaths. Finally, this implies that our algorithm
will eventually select every defining point of boundary $\mathcal{B}$, and
similarly, it will do the same with all other boundaries of $P$.

######  Theorem 2.10.

Leveraging Chan’s algorithm [7] for finding extreme points, the algorithm
selects every border point of $P$ in randomized expected time
$\mathcal{O}(nk\log{k})$ for $d=3$, and in

$\mathcal{O}\left(k(nk)^{1-\frac{1}{\lfloor
d/2\rfloor+1}}(\log{n})^{\mathcal{O}(1)}\right)$

time for all constant dimensions $d>3$.

Just as with Eppstein’s original algorithm, we can use Chan’s randomized
algorithm [7] for finding extreme points of point sets in $\mathbb{R}^{d}$, in
order to reduce the expected time complexity of our improved algorithm. The
remaining of the proof is the same as for Theorem 2.8.

#### Acknowledgements

Thanks to Prof. David Mount for pointing out Eppstein’s paper [12] and for the
valuable discussions on the results presented in this paper.

## References

  * [1] Fabrizio Angiulli.  Fast nearest neighbor condensation for large data sets classification.  IEEE Transactions on Knowledge and Data Engineering, 19(11):1450–1464, 2007.
  * [2] Sunil Arya, Guilherme D. da Fonseca, and David M. Mount.  Optimal approximate polytope membership.  In Proceedings of the Twenty-Eighth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 270–288. SIAM, 2017.
  * [3] Sunil Arya, Guilherme D. Da Fonseca, and David M. Mount.  Approximate polytope membership queries.  SIAM Journal on Computing, 47(1):1–51, 2018.
  * [4] Sunil Arya, Theocharis Malamatos, and David M. Mount.  Space-time tradeoffs for approximate nearest neighbor searching.  Journal of the ACM (JACM), 57(1):1, 2009.
  * [5] Oren Boiman, Eli Shechtman, and Michal Irani.  In defense of nearest-neighbor based image classification.  In 2008 IEEE Conference on Computer Vision and Pattern Recognition, pages 1–8. IEEE, 2008.
  * [6] David Bremner, Erik Demaine, Jeff Erickson, John Iacono, Stefan Langerman, Pat Morin, and Godfried Toussaint.  Output-sensitive algorithms for computing nearest-neighbour decision boundaries.  In Frank Dehne, Jörg-Rüdiger Sack, and Michiel Smid, editors, Algorithms and Data Structures: 8th International Workshop, WADS 2003, Ottawa, Ontario, Canada, July 30 - August 1, 2003. Proceedings, pages 451–461, Berlin, Heidelberg, 2003. Springer Berlin Heidelberg.  URL: https://doi.org/10.1007/978-3-540-45078-8_39, http://dx.doi.org/10.1007/978-3-540-45078-8_39 doi:10.1007/978-3-540-45078-8_39.
  * [7] Timothy M Chan.  Output-sensitive results on convex hulls, extreme points, and related problems.  Discrete & Computational Geometry, 16(4):369–387, 1996.
  * [8] Kenneth L Clarkson.  More output-sensitive geometric algorithms.  In Proceedings 35th Annual Symposium on Foundations of Computer Science, pages 695–702. IEEE, 1994.
  * [9] Corinna Cortes and Vladimir Vapnik.  Support-vector networks.  In Machine Learning, pages 273–297, 1995.
  * [10] T. Cover and P. Hart.  Nearest neighbor pattern classification.  IEEE Trans. Inf. Theor., 13(1):21–27, January 1967.  URL: http://dx.doi.org/10.1109/TIT.1967.1053964, http://dx.doi.org/10.1109/TIT.1967.1053964 doi:10.1109/TIT.1967.1053964.
  * [11] Luc Devroye.  On the inequality of cover and hart in nearest neighbor discrimination.  Pattern Analysis and Machine Intelligence, IEEE Transactions on, pages 75–78, 1981.
  * [12] David Eppstein.  Finding relevant points for nearest-neighbor classification.  In Symposium on Simplicity in Algorithms (SOSA), pages 68–78. SIAM, 2022.
  * [13] E. Fix and J. L. Hodges.  Discriminatory analysis, nonparametric discrimination: Consistency properties.  US Air Force School of Aviation Medicine, Technical Report 4(3):477+, January 1951.
  * [14] Alejandro Flores-Velazco.  Social distancing is good for points too!  In Proceedings of the 32st Canadian Conference on Computational Geometry, CCCG 2020, August 5-7, 2020, University of Saskatchewan, Saskatoon, Saskatchewan, Canada, 2020.
  * [15] Alejandro Flores-Velazco and David M. Mount.  Guarantees on nearest-neighbor condensation heuristics.  In Proceedings of the 31st Canadian Conference on Computational Geometry, CCCG 2019, August 8-10, 2019, University of Alberta, Edmonton, Alberta, Canada, 2019.
  * [16] Alejandro Flores-Velazco and David M. Mount.  Coresets for the Nearest-Neighbor Rule.  In Fabrizio Grandoni, Grzegorz Herman, and Peter Sanders, editors, 28th Annual European Symposium on Algorithms (ESA 2020), volume 173 of Leibniz International Proceedings in Informatics (LIPIcs), pages 47:1–47:19, Dagstuhl, Germany, 2020. Schloss Dagstuhl–Leibniz-Zentrum für Informatik.  URL: https://drops.dagstuhl.de/opus/volltexte/2020/12913, http://dx.doi.org/10.4230/LIPIcs.ESA.2020.47 doi:10.4230/LIPIcs.ESA.2020.47.
  * [17] Alejandro Flores-Velazco and David M. Mount.  Boundary-sensitive approach for approximate nearest-neighbor classification.  In Petra Mutzel, Rasmus Pagh, and Grzegorz Herman, editors, 29th Annual European Symposium on Algorithms, ESA 2021, September 6-8, 2021, Lisbon, Portugal (Virtual Conference), volume 204 of LIPIcs, pages 44:1–44:15. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2021.  URL: https://doi.org/10.4230/LIPIcs.ESA.2021.44, http://dx.doi.org/10.4230/LIPIcs.ESA.2021.44 doi:10.4230/LIPIcs.ESA.2021.44.
  * [18] Ruiqi Guo, Philip Sun, Erik Lindgren, Quan Geng, David Simcha, Felix Chern, and Sanjiv Kumar.  Accelerating large-scale inference with anisotropic vector quantization.  In International Conference on Machine Learning, 2020.  URL: https://arxiv.org/abs/1908.10396.
  * [19] Ruiqi Guo, Philip Sun, Erik Lindgren, Quan Geng, David Simcha, Felix Chern, and Sanjiv Kumar.  Accelerating large-scale inference with anisotropic vector quantization, 2020.  http://arxiv.org/abs/1908.10396 arXiv:1908.10396.
  * [20] Sariel Har-Peled.  A replacement for Voronoi diagrams of near linear size.  In Proc. 42nd Annu. IEEE Sympos. Found. Comput. Sci., pages 94–103, 2001.
  * [21] Norbert Jankowski and Marek Grochowski.  Comparison of instances selection algorithms I. Algorithms survey.  In Artificial Intelligence and Soft Computing-ICAISC 2004, pages 598–603. Springer, 2004.
  * [22] Jeff Johnson, Matthijs Douze, and Hervé Jégou.  Billion-scale similarity search with GPUs.  arXiv preprint arXiv:1702.08734, 2017.
  * [23] Jeff Johnson, Matthijs Douze, and Hervé Jégou.  Billion-scale similarity search with gpus, 2017.  http://arxiv.org/abs/1702.08734 arXiv:1702.08734.
  * [24] Kamyar Khodamoradi, Ramesh Krishnamurti, and Bodhayan Roy.  Consistent subset problem with two labels.  In Conference on Algorithms and Discrete Applied Mathematics, pages 131–142. Springer, 2018.
  * [25] Marc Khoury and Dylan Hadfield-Menell.  Adversarial training with Voronoi constraints.  CoRR, abs/1905.01019, 2019.  URL: http://arxiv.org/abs/1905.01019, http://arxiv.org/abs/1905.01019 arXiv:1905.01019.
  * [26] Nicolas Papernot and Patrick McDaniel.  Deep $k$-nearest neighbors: Towards confident, interpretable and robust deep learning.  arXiv preprint arXiv:1803.04765, 2018.
  * [27] Neehar Peri, Neal Gupta, W. Ronny Huang, Liam Fowl, Chen Zhu, Soheil Feizi, Tom Goldstein, and John P. Dickerson.  Deep $k$-nn defense against clean-label data poisoning attacks.  In European Conference on Computer Vision, pages 55–70. Springer, 2020.
  * [28] Jürgen Schmidhuber.  Deep learning in neural networks: An overview.  CoRR, abs/1404.7828, 2014.  URL: http://arxiv.org/abs/1404.7828, http://arxiv.org/abs/1404.7828 arXiv:1404.7828.
  * [29] Chawin Sitawarin and David Wagner.  On the robustness of deep $k$-nearest neighbors.  In 2019 IEEE Security and Privacy Workshops (SPW), pages 1–7. IEEE, 2019.
  * [30] Charles J. Stone.  Consistent nonparametric regression.  The annals of statistics, pages 595–620, 1977.
  * [31] Gordon Wilfong.  Nearest neighbor problems.  In Proceedings of the Seventh Annual Symposium on Computational Geometry, SCG ’91, pages 224–233, New York, NY, USA, 1991. ACM.  URL: http://doi.acm.org/10.1145/109648.109673, http://dx.doi.org/10.1145/109648.109673 doi:10.1145/109648.109673.
  * [32] D. Randall Wilson and Tony R. Martinez.  Instance pruning techniques.  In Proceedings of the Fourteenth International Conference on Machine Learning, ICML ’97, pages 403–411, San Francisco, CA, USA, 1997. Morgan Kaufmann Publishers Inc.  URL: http://dl.acm.org/citation.cfm?id=645526.657143.
  * [33] A. V. Zukhba.  NP-completeness of the problem of prototype selection in the nearest neighbor method.  Pattern Recog. Image Anal., 20(4):484–494, December 2010.  URL: http://dx.doi.org/10.1134/S1054661810040097, http://dx.doi.org/10.1134/S1054661810040097 doi:10.1134/S1054661810040097.


# Low-Threshold Lasing with Frozen Mode Regime and Stationary Inflection Point
in Three Coupled Waveguide Structure

Kessem Zamir- Abramovich1<EMAIL_ADDRESS>Nathaniel Furman2 Albert
Herrero-Parareda2 Filippo Capolino2 Jacob Scheuer1<EMAIL_ADDRESS>1School of Electrical Engineering Tel- Aviv University, Tel-Aviv, Israel
2Department of Electrical Engineering and Computer Science, University of
California, Irvine, California

(August 28, 2024)

###### Abstract

The frozen mode regime is a unique slow-light scenario in periodic structures,
where the flat-bands (zero group velocity) are associated with the formation
of high-order stationary points (aka exceptional points). The formation of
exceptional points is accompanied by enhancement of various optical properties
such as gain, Q-factor and absorption, which are key properties for the
realization of wide variety of devices such as switches, modulators and
lasers. Here we present and study a new integrated optical periodic structure
consisting of three waveguides coupled via micro-cavities and directional
coupler. We study this design theoretically, demonstrating that a proper
choice of parameters yields a third order stationary inflection point (SIP).
We also show that the structure can be designed to exhibit two almost-
overlapping SIPs at the center of the Brillouin Zone. We study the
transmission and reflection of light propagating through realistic devices
comprising a finite number of unit-cells and investigate their spectral
properties in the vicinity of the stationary points. Finally, we analyze the
lasing frequencies and threshold level of finite structures (as a function of
the number of unit-cells) and show that it outperforms conventional lasers
utilizing regular band edge lasing (such as DFB lasers).

††preprint: APS/123-QED

##  I Introduction

Optical structures operating in the slow-light regime exhibit group velocities
that are substantially lower than the speed of light in vacuum. These optical
structures exhibit numerous interesting properties; hence, they are highly
attractive for many applications. In particular, such structures have been
found to exhibit enhancement of properties such as gain, absorption and
quality factor (Q factor) [1, 2, 3, 4, 5]. The frozen mode regime is a special
case of slow light, which describes a solution of Bloch wave point with zero
group velocity point (stopped light); this point is an outcome of the
coalescence of Bloch waves (both eigenvalues and eigenvectors) at a single
frequency. The points are called stationary points or exceptional points, and
they appear in various types. These types differ in their dispersion
properties (the relation between the $k$ vector and the angular frequency
$\omega$) in the vicinity of these points [6].

In particular, there is a distinct difference between stationary points
corresponding to the coalescence of an even and an odd number of modes. The
latter are called stationary inflection points (SIPs). The lowest order of
this class corresponds to the coalescence of three eigenvalues, and they are
characterized by a cubic dispersion relation:

$\displaystyle\omega-\omega_{SIP}\propto(k-k_{SIP})^{3}$ (1)

Compared to stationary points with a coalescence of an even number of
eigenvalues and eigenvectors, SIPs exhibit several unique properties. More
specifically, SIPs are not formed at the band-edge, but rather within the
Brillouin zone (BZ). Since SIPs correspond to the coalescence of three Bloch
waves propagating in the _same_ direction, they do not form a standing wave.
This property renders SIPs interesting scientifically, as well as attractive
for various applications involving slow and stopped light. For example,
coupling light into and out of an optical structure supporting an SIP is more
efficient, as the excitation of counter-propagating modes can be eliminated
[7]. This is because in the vicinity of the SIP frequency, the $k$ vector
preserves its sign. As a result, counter-propagating waves can be suppressed,
contrary to what occurs for stationary points of even order, such as DBEs.
Moreover, the slow light regime helps match a slowly propagating mode to a
fast mode across the interface of the structure [8]. Furthermore, the SIP
resonance has been shown to be remarkably robust to structural disorder and
perturbation [9].

There are many ways to generate stationary points [10, 11, 12], one of them
being the use of optical periodic structures. In particular, previous studies
have presented several periodic structures supporting SIPs, such as waveguides
with multiple gratings [8, 13], coupled resonators optical waveguides, [14],
three-way periodic microstrip coupled waveguides [15], asymmetric serpentine
optical waveguides [16], and coupled transmission lines [17]. Besides,
previous studies have investigated lasers operating at the frozen mode and
stopped light regimes. Lasers operating at a regular band edge (RBE), such as
DFB laser, have been investigated thoroughly. A laser operating at a
degenerated band edge (DBE) was investigated in [18]. There is limited amount
of research on lasers utilizing odd-order stationary points. In Ref. [19], the
theory of unidirectional lasers operating in the frozen mode regime has been
proposed and studied theoretically. The laser was designed to lase in the
vicinity of an SIP, where non-reciprocity was introduced (by means of a
magnetic layer) to obtain unidirectional lasing. More recently, a new
structure, which also exhibits an SIP, was investigated theoretically as a
potential laser near an SIP [20]. This structure, employing an asymmetric
serpentine optical waveguide was shown to exhibit a lower lasing threshold
level than that of the same structure operating at an RBE instead.

Figure 1: Three periodic waveguides coupled with ring resonator and
directional coupler structure. The boundaries of the unit cell are marked by a
dashed red line.

In this paper, we propose and study a new, integrated, periodic structure that
exhibits SIPs, and can be highly attractive for low-threshold laser
applications. The structure consists of a repeating unit cell comprising three
parallel waveguides coupled through a directional coupler and a ring
resonator, as shown in Fig. 1. The central waveguide is coupled to the top
waveguide with a ring resonator, and to the bottom one with a directional
coupler. By properly setting the parameters of the structure - the resonator
radius ($R$), the length of the unit cell ($2d$), and the power coupling
coefficients of the coupler and the resonators
($\kappa_{1},\;\kappa_{2},\;\kappa_{3}$), it is possible to obtain a
dispersion relation exhibiting an SIP, and to control the properties of the
structure. We show that this structure is versatile, in the sense that it can
support SIPs at various frequencies, depending on the structure parameters. In
contrast to many of the previously proposed structures (such as those
presented in [16]), this structure’s unit cell is robust and relatively simple
to fabricate. It consists of simple and standard integrated optical elements,
whose designed parameters can be readily adjusted and modified to tune the
structures characteristics. We also perform a comprehensive analysis of the
power transmission/reflection characteristics of this structure in two cases:
1) when all six I/O ports are available and 2) when only one input port and
one transmission port are available. This is in contrast to previous studies
that only examined the latter case. In addition, we also analyze finite
segments of the proposed structure when gain material is incorporated. We
calculated the lasing threshold of this geometry, and find that the threshold
gain of our structure decreases proportionally to $N^{-3}$, where $N$ is the
number of unit cells. We learn that SIP-based devices outperform their RBE
counterparts in terms of lasing threshold level. Our results regarding the
properties of the threshold gain reinforce the results obtained in previous
studies [20, 19], especially in the context of SIP-based laser applications.

The rest of the paper is organized as follows: In Sec. II, we calculate the
dispersion relation of the structure by finding a general analytic expression
for the transfer matrix of a unit cell. The structure parameters are then
optimized in order to obtain SIPs. In Sec. III, we numerically calculate the
transmission and reflection for a finite structure. We discuss the spectral
properties of the structure with three input/output (I/O) ports, and also with
a single input and single transmission port. In Sec. IV, we analyze an active
version of this structure, incorporating optical gain material, and in Sec. V
we summarize the results and conclude. The time convention $e^{i\omega t}$ is
used throughout the paper.

##  II Dispersion Relation and SIPs

We present the derivation of the unit-cell transfer matrix, and obtain the
dispersion relation of the eigenmodes in the infinitely-long waveguiding
structure. To calculate the Bloch wavenumber dispersion relation of the
structure shown in Fig. 1, we use the transfer matrix method. For simplicity,
we divide the unit cell into two sections: two coupled waveguides with a third
separate waveguide seen in Fig. 2(a), and an Add-Drop multiplexer of length
$2d$ seen in Fig. 2(b). Here we assume the coupling sections are infinitely
small (e.g. point coupling, $\epsilon\rightarrow 0$).

Following the unit cell separation, obtaining the transfer matrices for each
section is straightforward. The transfer matrix for the three waveguides with
a single directional coupler $\kappa_{1}$ is in Eq. 2.

$\displaystyle\begin{bmatrix}E_{1}^{+}(\epsilon)\\\ E_{1}^{-}(\epsilon)\\\
E_{2}^{+}(\epsilon)\\\ E_{2}^{-}(\epsilon)\\\ E_{3}^{+}(\epsilon)\\\
E_{3}^{-}(\epsilon)\\\ \end{bmatrix}=\begin{bmatrix}t_{1}&0&r_{1}&0&0&0\\\
0&t_{1}&0&-r_{1}&0&0\\\ 0&0&t_{1}&0&r_{1}&0\\\ 0&0&0&t_{1}&0&-r_{1}\\\
0&0&0&0&1&0\\\ 0&0&0&0&0&1\\\ \end{bmatrix}\begin{bmatrix}E_{1}^{+}(0)\\\
E_{1}^{-}(0)\\\ E_{2}^{+}(0)\\\ E_{2}^{-}(0)\\\ E_{3}^{+}(0)\\\
E_{3}^{-}(0)\\\ \end{bmatrix}$ (2)

where $t_{1}\overset{\Delta}{=}\sqrt{1-\kappa_{1}}$,
$r_{1}\overset{\Delta}{=}-i\sqrt{\kappa_{1}}$. Here, $\kappa_{1}$ is the
intensity coupling between the two waveguides. $E_{n}$ is the electric field
phasor in the $n^{\mathrm{th}}$ waveguide.

Similarly, the transfer matrix for the section of the add-drop multiplexer is
written in Eq. 3, where $T_{1,2},\;D_{1,2}$ are the field transmission
functions of the Through and Drop port of two add-drop multiplexers (ADMs).
$T_{1,2},\;D_{1,2}$ are written in Eq. 4. Index one ($T_{1},\;D_{1}$)
corresponds to an ADM with coupling coefficients $\kappa_{3}$ at the through
port and $\kappa_{2}$ at the drop port. Index two ($T_{2},\;D_{2}$)
corresponds to the opposite case of ADM ($\kappa_{2}$ at the through port and
$\kappa_{3}$ at the drop port).

$\displaystyle\begin{bmatrix}E_{1}^{+}(2d)\\\ E_{1}^{-}(2d)\\\
E_{2}^{+}(2d)\\\ E_{2}^{-}(2d)\\\ E_{3}^{+}(2d)\\\ E_{3}^{-}(2d)\\\
\end{bmatrix}=\begin{bmatrix}e^{-2i\phi}&0&0&0&0&0\\\ 0&e^{2i\phi}&0&0&0&0\\\
0&0&(T_{2}-\frac{D_{1}D_{2}}{T_{1}})e^{-2i\phi}&0&0&\frac{D_{1}}{T_{1}}\\\
0&0&0&\frac{1}{T_{2}}e^{2i\phi}&-\frac{D_{1}}{T_{2}}&0\\\
0&0&0&\frac{D_{2}}{T_{2}}&(T_{1}-\frac{D_{1}D_{2}}{T_{2}})e^{-2i\phi}&0\\\
0&0&-\frac{D_{2}}{T_{1}}&0&0&\frac{1}{T_{1}}e^{-2i\phi}\\\
\end{bmatrix}\begin{bmatrix}E_{1}^{+}(0)\\\ E_{1}^{-}(0)\\\ E_{2}^{+}(0)\\\
E_{2}^{-}(0)\\\ E_{3}^{+}(0)\\\ E_{3}^{-}(0)\\\ \end{bmatrix}$ (3)
$\displaystyle T_{1}$
$\displaystyle=\frac{\sqrt{1-\kappa_{3}}-\sqrt{(1-\kappa_{2})(1-\alpha)}e^{-i\varphi}}{1-\sqrt{(1-\kappa_{3})(1-\kappa_{2})(1-\alpha)}e^{-i\varphi}}$
(4) $\displaystyle D_{1}$
$\displaystyle=\frac{-\sqrt{\kappa_{3}\kappa_{2}}(1-\alpha)^{\frac{1}{4}}e^{-i\frac{\varphi}{2}}}{1-\sqrt{(1-\kappa_{3})(1-\kappa_{2})(1-\alpha)}e^{-i\varphi}}$
$\displaystyle T_{2}$
$\displaystyle=\frac{\sqrt{1-\kappa_{2}}-\sqrt{(1-\kappa_{3})(1-\alpha)}e^{-i\varphi}}{1-\sqrt{(1-\kappa_{2})(1-\kappa_{3})(1-\alpha)}e^{-i\varphi}}$
$\displaystyle D_{2}$ $\displaystyle=D_{1}$

$\phi=k_{0}n_{w}\,d$, $\varphi=k_{0}n_{r}\,2\pi R$, are phase accumulations in
a strait waveguide with length $d$, and a ring resonator with radius $R$,
respectively. $n_{w}$ is the effective refractive index of the straight
waveguide, and $n_{r}$ is the effective refractive index of the curved ring
waveguide. $\alpha$ represents intensity loss per revolution inside the
resonator ring. For the dispersion relation analysis we assume $\alpha=0$. The
complete transfer matrix of the unit cell, M, is obtained by multiplying the
matrices of the two sections. By invoking the Bloch theorem, the dispersion
relation of the periodic structure is obtained from Eq. 5.

$\displaystyle|\textbf{M}-\textbf{I}e^{-ik2d}|=0$ (5)

where $||$ denotes the determinant operation, and I is the 6x6 identity
matrix.

(a)

(b)

Figure 2: The two sections of a unit cell: (a) the directional coupler, and
(b) the Add- drop multiplexer of length $2d$.

###  II.1 Separated Stationary Inflection Points

Eq. 5 provides the foundation for analyzing the eigenmodes in the structure
depicted in Fig. 1, leading to the dispersion relation of an infinitely-long
structure for a specific set of parameters
($d$,$R$,$\kappa_{1}$,$\kappa_{2}$,$\kappa_{3}$). Obtaining an SIP requires a
specific choice of parameters. Fig. 3(a) depicts the magnitude of the
determinant of the matrix $\textbf{M}-\textbf{I}\,e^{-ik\,2d}$ for the set of
parameters: ‘set 1’ listed in Table 1.

Table 1: Two sets of parameters for studied structure. Parameter | $\kappa_{1}$ | $\kappa_{2}$ | $\kappa_{3}$ | $\alpha$ | $2d[\upmu\textrm{m}]$ | $R[\upmu\textrm{m}]$ | $n_{w}$ | $n_{r}$
---|---|---|---|---|---|---|---|---
Set 1 | 0.30 | 0.16 | 0.56 | 0 | 44 | 7 | 2.21 | 2.21
Set 2 | 0.20 | 0.42 | 0.73 | 0 | 20 | 4 | 3.42 | 3.42

The blue-green lines in the figure, indicate magnitude close to zero, i.e.
solutions of Eq. 5, thus indicating the $k-\omega$ dispersion relation. The
red circles denote a potential SIP at $\lambda\approx 1.54\,\upmu\textrm{m}$,
related to frequency $\omega_{s}=\frac{2\pi\,c}{\lambda}$. Note that at this
point we obtain $\frac{\partial\omega}{\partial k}=0$, and a third-degree
polynomial in its vicinity, which is indicative of a frozen mode regime.
Nevertheless, this does not guaranties an SIP, as the coalescence of both
eigenvalues and eigenvectors still needed to be verified.

(a)

(b)

Figure 3: ‘set 1’ parameters dispersion relation (a) Logarithmic scale (using
base of ten) of the determinant $|\textbf{M}-\textbf{I}\,e^{-ik\,2d}|$, for
varying angular frequency and Bloch wavenumber. The SIP point is denoted with
a red circle. (b) normalized $k$’s related to the eigenvalues of the transfer
matrix M at a range of frequencies near the SIP. The magenta lines indicates
the propagating modes, marked as $\\{k_{p}\\}$. The light blue lines indicates
the decaying modes, marked as $\\{k_{d}\\}$. Black-dotted lines and solid
lines represent the imaginary and real parts of the normalized wave numbers,
respectively.

To verify that this is indeed a stationary point of order three, we consider
the eigenvalues of the transfer matrix M. Any input field in this structure
can be spanned by a basis of six propagating and evanescent waves which are
the eigenvectors of M. At a given frequency, the eigenvalues of M are
designated as $\\{\gamma\\}_{i=1}^{6}$, where the relation between the
eigenvalues and the corresponding wavenumbers $k$ is given by Eq. 6.

$\displaystyle\gamma=e^{-ik\,2d}$ (6)

A stationary point of third order (i.e., an SIP) is characterized by the
coalescence of three eigenvalues and three eigenvectors. Fig. 3(b) plots the
dispersion relation of the structure (normalized complex wavenumber
$\frac{k\,2d}{\pi}$ vs. angular frequency). At each frequency, there are six
Bloch modes with six eigenvalues. The magenta lines in Fig. 3(b) indicate
propagating modes. The solid magenta lines correspond to the real parts of the
eigenvalues, while the doted ones indicate the imaginary parts. These
eigenvalues are purely real, and the real parts of the eigenvalues match
exactly the green lines near the SIP in Fig. 3(a). The light blue lines in
Fig. 3(b) indicate the non-propagating (decaying) modes. Similar to the
propagating modes, the solid and dotted lines indicate the real and imaginary
parts of the wavenumbers, respectively. Note that there are four decaying
modes: two that are complex conjugates. In other words, two modes have
eigenvalues with a positive real part, and the other two with a negative real
part. From Fig. 3(b) we identify the SIP frequency as the frequency at which
three of the wavenumbers coalesce into a single real value. It corresponds to
the frequency marked by a red circle in Fig. 3(a), $\omega=\omega_{s}$. Due to
reciprocity, two SIPs are obtained (at the same frequency), corresponding to
two triply- degenerate counter-propagating modes.

We emphasize that an SIP is a condition where both eigenvalues and
eigenvectors coalesce. However, Fig. 3(b) only shows the coalescence of the
eigenvalues. The coalescence of the eigenvectors is verified through the
coalescence parameter that was defined in [15, 16] for the SIP (see Appendix A
for details). Appendix A also describes the structure design and optimization
process used for obtaining the parameters in Table 1. Consequently, it can be
inferred that Fig. 3(b) verifies that the marked points in Fig. 3(a) are
indeed SIPs, meaning that $\omega_{s}=\omega_{SIP}$.

###  II.2 Almost-Overlapping Stationary Inflection Points

The proposed structure can exhibit SIP dispersion properties in a range of
frequencies and locations in the BZ, controlled by the design parameters.
Modifying the parameters ($d$,$R$,$\kappa_{1}$,$\kappa_{2}$,$\kappa_{3}$),
depicted in Fig. 1, enables us to generate an SIP in different locations in
the dispersion relation. This section considers the set of parameters
designated as ‘set 2’ in Table 1. As in Sec. II.1, Fig. 4(a) plots the
colormap of the determinant of $\textbf{M}-\textbf{I}\,e^{-ik\,2d}$. Fig. 4(b)
depicts the Bloch wavenumbers of the transfer matrix as a function of
frequency. At first glance, it seems that a flat band with zero group velocity
is formed at the center of the Brillouin Zone (BZ) with wavelength
$\lambda\approx 1.55\,\upmu\textrm{m}$, represented with a red circle. We
define the frequency related to this wavelength as before, $\omega_{s}$. In a
scenario of a flat band, the transfer matrix of the structure should exhibit a
degeneracy of six eigenvalues related to $k=0$. However, the eigenvalues
depicted in Fig. 4(b) indicate that this so-called ”flat band” consists, in
fact, of two very close SIPs located at both sides of the BZ center. This
should not be surprising as it is well known that an intersection of different
spectral branches at $k=0$ is only possible for a 1D periodic structure that
exhibits glide plane symmetry [21, 22, 23, 24]. However, as the structure
depicted in Fig. 1 does not exhibit glide symmetry, the two SIPs cannot fully
intersect at the center of the BZ and therefore cannot support degenerate
SIPs.

Thus, we understand that the studied structure can support the formation of
SIPs over large domains in the BZ, and at different wavelengths. The
frequencies and wavenumbers of such SIPs can be controlled by modifying and
optimizing the set of parameters that describe the structure. In Sec. II.1 we
showed an SIP that is located at $\frac{k\,2d}{\pi}\approx 0.2$, which makes
this SIP separated from its backward propagating counterpart. In this section
we described an SIP that is formed almost at the center of the BZ, at
$\frac{k\,2d}{\pi}\approx 0.05$. This shows that this structure is incredibly
versatile, since it can generate SIPs at various locations in the BZ.

(a)

(b)

Figure 4: ‘set 2’ parameters dispersion relation (a) Logarithmic scale (using
base of ten) of the determinant $|\textbf{M}-\textbf{I}\cdot e^{-ik\,2d}|$, ,
for varying angular frequency and Bloch wavenumber. The SIP is marked with a
red circle (b) normalized $k$’s related to the eigenvalues of the transfer
matrix $M$ at a range of frequencies near the SIP. The magenta lines indicate
propagating modes, marked as $\\{k_{p}\\}$. The light blue lines indicate
decaying modes, marked as $\\{k_{d}\\}$. Black-dotted lines and solid lines
correspond to the imaginary and real parts of the normalized wave numbers,
respectively.

##  III Finite Length Structure

###  III.1 Single Input - Multiple Output Transmission and Reflection

The dispersion relation and frozen mode regimes discussed in Sec. II are
obtained only for infinitely-long structures. However, any realistic structure
is finite in length, thus leading to resonances with optical properties that
may differ a bit from those of the infinitely-long structure. Resonances in
such a realistic finite-length structure may not support a “perfect” frozen
light mode. However, as more unit cells are added, the closer the SIP
resonance is to the SIP frequency, and the smaller the group velocity is at
the resonance. Thus, the longer the structure, the closer the characteristics
of the propagating waves in the structure become to those expected at an SIP.
Therefore, we expect to obtain enhancement of properties, such as the Q
factor, by using finite-length structures with a sufficient number of unit
cells.

In order to investigate the transmission and reflection properties of a
finite-length structure, we analyze a waveguide with length $l=2d\,N$ where
$N$ is the number of unit cells. Any resonance response of such structure in
the vicinity of an SIP is expected to be enhanced due to the small group
velocity. In addition, the larger the number of unit cells, the closer the
resonance frequency to the SIP. Recall that the structure consists of three
parallel waveguides, (see Fig. 1) with three I/O ports. Consequently, the
structure can be excited in many ways. We analyze the power transmission and
reflection properties, for the case where an input field excites only the
upper waveguide (port 3) from the left. We use the transfer matrix approach to
obtain the intensity of the fields exiting from the three left-hand side
waveguides (i.e. reflection), and from the three right-hand side waveguides
(transmission) of the structure. Eq. 7 defines the transmission and reflection
power, which are denoted as $T_{m}^{2}$ and $R_{m}^{2}$ respectively. The plus
and minus signs indicate forward and backward propagating fields, and the
subscript $m$ denotes at which port the field is calculated. The number in the
brackets is the point on the propagation axes where the field is calculated.
Fig. 2 shows the port numbers.

$\begin{matrix}T_{m}^{2}=\left|\frac{E_{m}^{+}(l)}{E_{3}^{+}(0)}\right|^{2}&R_{m}^{2}=\left|\frac{E_{m}^{-}(0)}{E_{3}^{+}(0)}\right|^{2}&m=1,2,3\end{matrix}$
(7)

(a)

(b)

Figure 5: (a) Transmission and reflection from six ports of a finite-length
waveguide structure with 60 unit cells. (b) Transmission from input port 3 to
output port 3, near the SIP frequency. The parameters in ‘set 1’ are used.

Fig. 5(a) shows the transmission and reflection spectra of a finite-length
structure with $N=60$, excited through the upper waveguide input (port 3). In
this figure, the SIP frequency obtained for ‘set 1’ is designated as
$\omega_{SIP}$. It can be seen that the transmitted signal through port 3 and
those reflected through ports 1 and 2 exchange power as a function of
frequency, while the signals at the other ports are zero. At frequencies
around the SIP, the spectral profiles in the transmission and reflection ports
become more oscillatory and exhibit narrower peaks. This is clearly seen in
Fig. 5(b), which is a zoom-in of $T^{2}_{3}$ (marked in the black rectangle in
Fig. 5(a)). At $\omega\approx 0.999\omega_{SIP}$, the fast oscillations in the
spectral responses stop and a resonance in the reflection from port
$R^{2}_{1}$ is formed, accompanied by zero transmission through ports
$T^{2}_{3}$ and a complementary decrease in $R^{2}_{2}$. At $\omega_{SIP}$,
all the power is reflected through the lower waveguide ($R^{2}_{1}$). The
observed increase in the density of the spectral peaks in the vicinity of the
SIP indicates that the group velocity is indeed smaller close to the SIP
frequency. It also implies that by introducing optical gain into the
structure, it should be possible to obtain large optical amplification and
even lasing.

###  III.2 Single Input - Single Output Transmission Port

Figure 6: Finite-length model, closed with reflective mirrors. This
arrangement has one input port and one transmission port.

In order to obtain an SIP, it is necessary for the structure to support at
least three modes in each direction, thus necessitating the three-path
waveguide arrangement. However, for any practical application such as high-
gain amplifier or laser, it is advantageous to have a single input and two
outputs, reflection (from the input port) and transmission from a second port.
This can be achieved by closing the I/O ports at the end of the (finite)
structure as shown in Fig. 6: Input ports 1 and 2 and output ports 2 and 3 are
closed by introducing reflectors which can be implemented as distributed Bragg
reflectors (DBRs). A DBR is a natural choice for implementing an integrated
waveguide reflector. In this scenario, we chose the reflectors to exhibit
perfect reflectivity. Different DBR designs facilitate high reflectivity as
well as control over the reflected spectral band [25, 26]. The closed
structure possesses two ports which can be considered as I/O for optical
amplification purposes or two output ports for lasing (similar to a Fabry-
Perot laser). Consequently, the reflection and transmission ports of the
structures are at port 3 (left) and at port 1 (right), respectively. The
choice between different waveguides for the I/O ports stems from the fact that
an SIP formation necessitates coupling between forward and backward
propagating waves. The position of the DBRs ensures that the input signal
propagates through all the waveguides in both directions. Nevertheless,
additional configurations for “closing” the structure (i.e. obtaining single
input and single output) are possible. The study of such configurations is
beyond the scope of this paper. In this finite-length model, we assume no
losses, therefore due to energy conservation the sum of the transmission and
the reflection powers equals one.

Fig. 7 plots the transmission spectrum $T^{2}$, for two structure lengths of
30 and 70 unit cells. There are several things to note in the transmission
spectrum of the device. The transmission exhibits oscillations similar to a
Fabry-Perot (FP) cavity (compared to the spectral properties of the ‘open’
structure shown in Fig. 5(a)). This result is not surprising, since the
closing of some of the ports in the structure introduces feedback that couples
between forward and backward propagating waves (as in an FP cavity). However,
in the vicinity of the SIP frequency, these oscillations corresponding to the
FP resonances become dense and sharper (higher Q factor). This is attributed
to the slow-light effect in the vicinity of the SIP which effectively reduces
the local free spectral range (FSR) of the FP cavity. The resonances of the
structure consisting of $N=70$ unit cells are denser and sharper than those of
the $N=30$ device. This is due to the longer FP cavity, which, similar to the
slow-light effect, reduces the FSR. This effect is particularly strong near
the SIP frequency (shown in the inset of Fig. 7), where the resonances of the
longer structure are much sharper than those of the shorter one due to the
combination of the longer device and slow-light effect.

Figure 7: Transmission of closed model with mirrors, for 30 and 70 unit cells.

##  IV Analyzing Gain Properties

The sharp resonances associated with high quality factor and the SIP
frequency, indicate that such a structure could be useful for the realization
of low-threshold lasers [27, 28, 20]. To identify the lasing threshold of the
structure we introduce an imaginary part to the refractive index of the
resonators and waveguides, and consider the dependence of the threshold level
on the number of unit cells. The lasing condition is a singular point, which
is characterized as a pole in the transmission spectral response. In other
words, at the lasing threshold, the output power calculated using the transfer
matrix, approaches infinity regardless of the input intensity level, as shown
in the book [29], chapter 5. In practice, gain saturation effects, which are
often not considered in a linear transfer matrix analysis, limit the actual
output power. Nevertheless, such effects do not affect the lasing threshold.

As the structure supports many resonances that can potentially lase, we
consider the resonance which is closest to the SIP frequency because it is the
one that retains the SIP properties the most. As the group velocity in the
vicinity of this frequency is lowest, the intensity buildup at this frequency
is expected to be enhanced substantially and exhibit the lowest threshold
level. We continue to focus on the SIP obtained with ‘set 1’ parameters. Fig.
8(a) depicts the resonance (i.e., the lasing) angular frequency (denoted by
magenta dots) as function of the number of unit cells. The resonance frequency
approaches $\omega_{s}$, which corresponds to $\omega_{SIP}$ frequency for the
magenta points, as the number of the unit-cells increases.

As a comparison, we also calculate the threshold level of a structure
exhibiting a regular band edge (RBE) of a CROW structure [30] comprising ring
resonators with radius $R=7\mu\textrm{m}$ (which is the same as the one in
‘set 1’ Table 1 here considered) coupled to each other with the coupling
coefficient $\kappa_{r}=\kappa_{2}$ from set 1. The unit cell length of the
CROW is $l=4R$, which is smaller than the unit cell length of our structure.
Any realistic structure exhibits losses that should be included in the
calculation of the lasing threshold. Thus, we introduce a loss parameter of
$\alpha=0.01$ to the ring resonators, equivalent to a 1 percent power loss per
revolution.

Although the comparison is not perfect, it does allow us to get a better
understanding of the lasing properties of a device operating in the vicinity
of an SIP. The process of calculating the threshold gain for a finite
structure is detailed in Appendix B.

(a)

(b)

Figure 8: Dependence of the threshold gain in the waveguides, and the lasing
frequency versus the number of unit cells of the structure. (a) The angular
resonance frequency $\omega_{res}$, which is the closest to the stationary
point $\omega_{s}$, as a function of the number of unit cells. $\omega_{s}$
refers to the stationary point frequency, either $\omega_{SIP}$ (for the
magenta dots) or $\omega_{RBE}$ (for the blue dots). Magenta points: change in
resonance frequency normalized by the SIP frequency in our structure. Blue
points: change in resonance normalized by the RBE frequency in the CROW. (b)
The threshold gain for RBE in CROW model, and SIP in the proposed model versus
the number of unit cells. The y-axis is in logarithmic scale. Magenta points:
SIP threshold gain, Blue points: RBE threshold gain. The two insets show
polynomial fitting of the SIP (magenta) and the RBE (blue). Note that the y-
axis in those figures is in linear scale.

Fig. 8(a) also shows the lasing frequency of a CROW structure as a function of
the number of unit-cells, $N$ (blue dots). The lasing frequencies are
normalized to either $\omega_{SIP}$ or to $\omega_{RBE}$, according to the
relevant structure. Note that in this case $\omega_{RBE}\approx
1224(\textrm{rad/ps})$. It can be seen, that the number of unit cells $N$ has
more impact on the lasing frequency of the SIP supporting structure than on
that of the RBE. However, although the lasing frequency of the RBE supporting
structure seems to be independent of $N$, this is not the case. The inset of
the figure shows a zoom-in on the lasing frequency of the RBE structure,
indicating that it also approaches $\omega_{RBE}$ as $N$ increases. In both
cases, the dependence of the resonances on $N$ can be described as
$\omega_{th}\propto\alpha+\beta N^{-1}+\gamma N^{-3}$. Fig. 8(b) depicts a
comparison between the threshold gain (in the waveguides) near the SIP (in
pink) of the structure depicted in Fig. 6, and that corresponding to the RBE
in the CROW. The parameters of ’set 1’ in Table 1 have been used in the
comparison. The threshold gain is presented as a function of the number of
unit cells in each structure. The scale of the $y$ axis is in logarithmic
scale. There are a couple of important points to be noted: first, the lasing
threshold gain of the SIP supporting structure decreases rapidly with the
number of unit cells, reaching a minimal level of $\approx 4.5$ dB/cm at
$N=64$, and second is that the lasing threshold saturates for longer
structures. We attribute the threshold gain saturation to the fact that the
loss in structure is also proportional to the number of unit cells. Thus, for
a sufficiently long structure, this loss mechanism dominates all other
mechanisms in the structure (e.g. output coupling loss), and the lasing
threshold is obtained when the gain becomes equal to the loss at each unit
cell. In contrast, the dependence of the threshold gain on $N$ for the RBE
case seems to be approximately constant. However, zooming in on the RBE lasing
threshold $N$ dependence (shown as the left inset in Fig. 8(b)) indicates that
this is not the case. The threshold gain near the RBE of the CROW is not
constant but rather decreases slowly with $N$, at a slope that is
substantially slower than that of the SIP supporting structure.

Second, by fitting the two calculated curves (see insets of Fig. 8(b)) we find
that the threshold gain dependence on $N$ in both cases can be described by a
$3^{rd}$ order polynomial in $N^{-1}$: $a+b\,N^{-1}+c\,N^{-3}$. This
dependence agrees with previous results [20, 14, 6]. The main difference
between the curves is the value of the fitting coefficients $b$ and $c$
corresponding to the dependence on $N^{-1}$ and $N^{-3}$ respectively. In the
CROW structure we find that $|b_{RBE}|\approx 10^{-1}\textrm{dB/cm}$ and
$c_{RBE}\approx 2\cdot 10^{4}\textrm{dB/cm}$, while in the SIP supporting
structure we find that $|b_{SIP}|\approx 2\cdot 10^{-4}\textrm{dB/cm}$ and
$c_{SIP}\approx 2\cdot 10^{7}\textrm{dB/cm}$. There is a substantial
difference (of orders of magnitude) between these coefficients, which has a
dramatic impact on the properties of the two stationary points (the SIP and
RBE). The decrease of the lasing threshold as a function of $N$ is
significantly faster for the structure supporting an SIP. This is seen by the
larger fitting coefficient of the $N^{-3}$ term and the smaller coefficient of
the $N^{-1}$ term for that structure. Consequently, the dominant dependence of
the lasing threshold of the SIP structure and the CROW is $N^{-3}$ and
$N^{-1}$, respectively. It should also be noted that although the threshold
value of the SIP structure is larger than that of the CROW for $N<54$, once
the structure length exceeds 54 unit cells, the SIP lasing threshold becomes
smaller. Although a direct comparison between the structures is difficult
because they are very different from each other, the overall trend indicates
the advantages of operating near a SIP. Clearly, a laser operating at an SIP
can potentially exhibit a lower lasing threshold than that of conventional
lasers.

##  V Conclusion

In this paper we introduced a new periodic structure composed of three
parallel waveguides that are coupled to each other by ring resonators and
directional couplers. We calculated the dispersion relation of the structure
by means of the transfer matrix method, and showed that a proper choice of the
structure parameters (coupling coefficients, ring radius and the length of the
unit cell), leads to the formation of SIPs in the dispersion relation. We also
showed that it is possible to control the frequency and wavenumber of the
SIPs. Moreover, we studied the properties of finite-length structures (i.e.
with a finite number of unit cells). We calculated the spectral transmission
and reflection at each port, and studied their properties in the vicinity of
the SIP frequencies. We then studied the transmission and reflection
properties of a closed structure, exhibiting only two ports, by introducing
reflectors at the other four ports. This is only one of the possibilities of
closing such a structure, and other closed-structure configurations might lead
to different spectral transmission and reflection. Specifically, we focused on
two structures with $N=30$ and $N=70$ unit cells. We found that longer
structures yield sharper resonances exhibiting higher Q-factors. We attribute
this to the fact that the longer the structure, the closer the resonance gets
to the zero group-velocity point in the dispersion relation. Finally, we
analyzed the ability of this structure to serve as a laser by introducing
optical gain into the structure. Specifically, we focused on the dependence of
the lasing threshold gain on the number of unit cells, $N$. We compared the
lasing threshold gain at the SIP frequency and that of a structure supporting
a second-order stationary point (i.e., RBE). We found that although the
dependence of the threshold gain on $N$ is similar in both cases and can be
described by a $3^{rd}$ order polynomial in $N^{-1}$, the threshold level of
SIP-supporting structure decreases faster with increasing $N$. Thus, we
conclude that this structure is highly attractive for low-threshold laser
devices.

##  Appendix A Appendix A: Optimization Process

This appendix describes the design and optimization process of the SIP-
supporting structure. For this design, there are five independent parameters
that need to be determined (and optimized): the power coupling coefficients
($\kappa_{1},\kappa_{2},\kappa_{3}$), the ring radius ($R$) and a length
parameter, which is related to the length of the unit cell ($d$). There is an
additional degree of freedom, which is the position of the directional coupler
(in the longitudinal direction) between the middle and the lowest waveguide.
Here, we choose to position it at the center of the unit cell, equally
distanced from the two micro-rings. Nevertheless, this is just a specific
choice. Changing the position of the directional coupler brakes a mirror
symmetry, in the z direction inside the unit cell, which might allow us to get
another SIP, because it is a non-symmetric stationary point. Note that this
structure does not exhibit any symmetry in the vertical axes, which is a
requirement for obtaining an SIP.

Since there are many parameters needed to be determined, a numerical
optimization approach is required. To determine the design parameters we
minimize the coalescence parameter in Eq. 8, introduced in [31]. This
parameter quantifies the degree of coalescence between a number of complex
vectors, determined by the order of the stationary point -$SP_{order}$. For
example, in the case of an SIP, the degree of coalescence of three
eigenvectors of the transfer matrix $\\{\phi_{k}\\}$, having the same
eigenvalue, is minimized. When the angles between all of the eigenvectors
approach zero (i.e. all the eigenvectors are parallel), the coalescence
parameter approaches zero as well, i.e., $D_{h}\rightarrow 0$. At the start of
the optimization process, the length of the unit cell and the radius of the
micro-ring are chosen randomly. The coupling coefficients are then varied
manually by means of trial and error until a potential SIP is formed in the
dispersion relation at a certain wavelength $\lambda_{p}$. Then, the coupling
coefficients are optimized numerically in an attempt to minimize $D_{h}$ for
the wavelength $\lambda_{p}$. Here we employ the Nelder-Mead simplex algorithm
for the optimization [32]. Our objective function is:
$f(\\{\phi_{k}\\})=D_{h}\left(\\{\phi_{k}\\}\right)+10^{2}\cdot(degeneracy\left(\gamma\right)-degeneracy_{desired})$.
The term $D_{h}$ corresponds to the coalescence parameter of a group of
eigenvectors having the same eigenvalue $\gamma$. The _‘degeneracy’_ parameter
corresponds to the degeneracy (i.e. algebraic multiplicity) of the eigenvalue
$\gamma$, related to the group of eigenvectors $\\{\phi_{k}\\}$. The
‘$degeneracy_{desired}$’ parameter is the desired order of the degeneracy of
the eigenvalue $\gamma$ (e.g. $degeneracy_{desired}=3$ for an SIP). Therefore,
the use of this objective function favors the formation of
$degeneracy_{desired}$ identical eigenvectors and eigenvalues.

$\displaystyle D_{h}=\frac{1}{{SP_{order}\choose
2}}\;\;\sum_{\mathclap{\begin{subarray}{c}\;\;\;\;\;\;m=1,n=1\\\
m>n\end{subarray}}}^{SP_{order}}\sin(\theta_{n,m})\phantom{\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
}\cos(\theta_{n,m})=\frac{\langle\phi_{n}|\phi_{m}\rangle}{\|\phi_{n}\|\|\phi_{m}\|}$
(8)

Figure 9: Convergence of optimization process. The $x,y,z$ axes denote the
intensity coupling coefficients. The size of the dots and their color indicate
the value of $D_{h}$ in dB. The bigger and darker the dots are the smaller
$D_{h}$ is.

Fig. 9 depicts the convergence of the optimization process that yielded the
coupling coefficients of ‘set 1’. In this case, $\lambda_{p}$ is the
wavelength of the SIP $\lambda\approx 1.54\upmu$m, and the ring radius and the
unit cell length are defined in Table 1. The dots in the figure indicate the
position in the 3D coupling coefficients space
$\\{\kappa_{1},\kappa_{2},\kappa_{3}\\}$ of the structure during the
optimization process. The size and the color of the dots, indicate the value
of the coalescence parameter. The small yellow dot located at
$\\{\kappa_{1},\kappa_{2},\kappa_{3}\\}=\\{0.32,0.15,0.55\\}$ indicates
$D_{h}$ in the initial point (obtained manually as described above) of the
optimization process. At this point, $D_{h}>0.3$, indicating that the
eigenvectors have not coalesced yet. The dots show the progress of the
optimization process, indicating the trajectory in the coupling coefficients
space and the corresponding coalescence parameter. At the end of the process
we get $D_{h}\approx-30\,\textrm{dB}$ (dB stands for
$10\,\textrm{log}_{10}(D_{h})$).

##  Appendix B Appendix B: Calculating the threshold gain

In this appendix, we present the process for calculating the threshold gain
level and the resonance frequency for a finite-length structure with single
input and transmission ports, as depicted in Fig. 6. A uniformly distributed
gain is added to the structure, by introducing an imaginary part ($n_{i}$) to
the refractive index of the waveguides and the ring resonators. First, we find
the amplitude transmission of the structure as a function of the frequency and
gain. We then plot the absolute value of the amplitude transmission
$\textrm{Log}_{10}|T^{2}|$ as a function of these parameters, as depicted in
Fig. 10 for a structure comprising $56$ unit cells (in logarithmic scale).

Figure 10: Log scale of the transmission power for a finite-length structure
with $56$ unit cells, and single Input port and transmission port. The
transmission is plotted as a function of the angular frequency and the
imaginary part of the effective refractive index of a uniform (single mode)
waveguide. Three lasing conditions are visible as the narrow frequency sharp
peaks.

This plot allows for the identification of the resonance (i.e., lasing)
frequencies and and their corresponding lasing thresholds. A threshold
condition corresponds to a pole in the transmission function, which is
manifested as a peak in the transmission function. Next, we identify the
resonance frequency closest to the SIP frequency (the SIP frequency was
obtained by using the optimization process described in Appendix A). While the
resonance frequency does not coincide with the SIP frequency, it approaches it
for a very large number of unit cells, i.e.,
$\omega_{res}\xrightarrow{N\rightarrow\infty}\omega_{SIP}$ [5]. The gain level
at which the transmission is maximal is the threshold gain of the structure at
the relative resonance (i.e., lasing) frequency. The lasing threshold is
identified as the gain at which the transmission is maximal at the resonance
frequency of interest.

## References

  * Kiyota  _et al._ [2006] K. Kiyota, T. Kise, N. Yokouchi, T. Ide, and T. Baba, Various low group velocity effects in photonic crystal line defect waveguides and their demonstration by laser oscillation, Applied physics letters 88, 201904 (2006).
  * Ek  _et al._ [2014] S. Ek, P. Lunnemann, Y. Chen, E. Semenova, K. Yvind, and J. Mork, Slow-light-enhanced gain in active photonic crystal waveguides, Nature Communications 5, 5039 (2014).
  * Mortensen and Xiao [2007] N. A. Mortensen and S. Xiao, Slow-light enhancement of beer-lambert-bouguer absorption, Applied Physics Letters 90, 141108 (2007).
  * Song  _et al._ [2005] B. S. Song, S. Noda, T. Asano, and Y. Akahane, Ultra-high-q photonic double-heterostructure nanocavity, nature materials 4, 207 (2005).
  * Noh  _et al._ [2010] H. Noh, J.-K. Yang, I. Vitebskiy, A. Figotin, and H. Cao, Giant resonances near the split band edges of two-dimensional photonic crystals, Phys. Rev. A 82, 013801 (2010).
  * Figotin and Vitebskiy [2011] A. Figotin and I. Vitebskiy, Slow wave phenomena in photonic crystals, Laser & Photonics Reviews 5, 201 (2011).
  * Figotin and Vitebskiy [2006] A. Figotin and I. Vitebskiy, Slow light in photonic crystals, Waves in Random and Complex Media 16, 293 (2006).
  * Gutman  _et al._ [2012] N. Gutman, C. M. de Sterke, A. A. Sukhorukov, and L. C. Botten, Slow and frozen light in optical waveguides with multiple gratings: Degenerate band edges and stationary inflection points, Physical Review A 85, 033804 (2012).
  * Li  _et al._ [2017] H. Li, I. Vitebskiy, and T. Kottos, Frozen mode regime in finite periodic structures, Physical Review B 96, 180301 (2017).
  * Baba and Mori [2007] T. Baba and D. Mori, Slow light engineering in photonic crystals, Journal of Physics D: Applied Physics 40, 2659 (2007).
  * Lukin and Imamogulu [2001] M. D. Lukin and A. Imamogulu, Controlling photons using electromagnetically induced transparency, NATURE 413, 273 (2001).
  * Kash  _et al._ [1999] M. M. Kash, V. A. Sautenkov, A. S. Zibrov, L. Hollberg, G. R. Welch, M. D. Lukin, Y. Rostovtsev, E. S. Fry, and M. O. Scully, Ultraslow group velocity and enhanced nonlinear optical effects in a coherently driven hot atomic gas, Physical Review Letters 82, 5229 (1999).
  * Furman  _et al._ [2023] N. Furman, T. Mealy, M. S. Islam, I. Vitebskiy, R. Gibson, R. Bedford, O. Boyraz, and F. Capolino, Frozen mode regime in an optical waveguide with a distributed bragg reflector, J. Opt. Soc. Am. B 40, 966 (2023).
  * Nada  _et al._ [2017] M. Y. Nada, M. A. Othman, and F. Capolino, Theory of coupled resonator optical waveguides exhibiting high-order exceptional points of degeneracy, Physical Review B 96, 184304 (2017).
  * Nada  _et al._ [2020] M. Y. Nada, T. Mealy, and F. Capolino, Frozen mode in three-way periodic microstrip coupled waveguide, IEEE Microwave and Wireless Components Letters 31, 229 (2020).
  * Herrero-Parareda  _et al._ [2022] A. Herrero-Parareda, I. Vitebskiy, J. Scheuer, and F. Capolino, Frozen mode in an asymmetric serpentine optical waveguide, Advanced Photonics Research 3, 2100377 (2022).
  * Othman  _et al._ [2016a] M. A. K. Othman, F. Yazdi, and F. Capolino, Exceptional points of degeneracy in coupled-mode periodic structures, IEEE Xplore (2016a).
  * Veysi  _et al._ [2018] M. Veysi, M. A. Othman, A. Figotin, and F. Capolino, Degenerate band edge laser, Physical Review B 97, 195107 (2018).
  * Ramezani  _et al._ [2014] H. Ramezani, S. Kalish, I. Vitebskiy, and T. Kottos, Unidirectional lasing emerging from frozen light in nonreciprocal cavities, Physical review letters 112, 043904 (2014).
  * Herrero-Parareda  _et al._ [2023] A. Herrero-Parareda, N. Furman, T. Mealy, R. Gibson, R. Bedford, I. Vitebskiy, and F. Capolino, Lasing at a stationary inflection point, Optical Materials Express 13, 1290 (2023).
  * Hessel  _et al._ [1973] A. Hessel, M. H. Chen, R. C. Li, and A. A. Oliner, Propagation in periodically loaded waveguides with higher symmetries, Proceedings of the IEEE 61, 183 (1973).
  * QUEVEDO-TERUEL  _et al._ [2021] O. QUEVEDO-TERUEL, Q. CHEN, F. MESA, N. J. G. FONSECA, and G. VALERIO, On the benefits of glide symmetries for microwave devices, IEEE Journal of Microwaves 1, 457 (2021).
  * Patil  _et al._ [2022] C. M. Patil, G. Arregui, M. Mechlenborg, X. Zhou, H. Alaeian, P. D. García, and S. Stobbe, Observation of slow light in glide-symmetric photonic-crystal waveguides, Opt. Express 30, 12565 (2022).
  * Tamayo-Dominguez  _et al._ [2019] A. Tamayo-Dominguez, J.-M. Fernandez-Gonzalez, and O. Quevedo-Teruel, One-plane glide-symmetric holey structures for stop-band and refraction index reconfiguration, Symmetry 11 (2019).
  * Thomas and Pang [1996] I. Thomas, S. and S. W. Pang, Dry etching of horizontal distributed bragg reflector mirrors for waveguide lasers, Journal of Vacuum Science and Technology B: Microelectronics and Nanometer Structures Processing, Measurement, and Phenomena 14, 4119 (1996).
  * Pruessner  _et al._ [2008] M. W. Pruessner, T. H. Stievater, and W. S. Rabinovich, In-plane microelectromechanical resonator with integrated fabry–pérot cavity, Applied Physics Letters 92 (2008).
  * Othman  _et al._ [2016b] M. A. K. Othman, F. Yazdi, A. Figotin, and F. Capolino, Giant gain enhancement in photonic crystals with a degenerate band edge, PHYSICAL REVIEW B 93, 024301 (2016b).
  * Nada  _et al._ [2018] M. Y. Nada, M. A. K. Othman, O. Boyraz, and F. Capolino, Giant resonance and anomalous quality factor scaling in degenerate band edge coupled resonator optical waveguides, JOURNAL OF LIGHTWAVE TECHNOLOGY 36, 3031 (2018).
  * Numai [2004] T. Numai,  _Fundamentals of semiconductor lasers_ (Springer, 2004).
  * Poon  _et al._ [2004] J. K. S. Poon, J. Scheuer, S. Mookherjea, G. T. Paloczi, Y. Huang, and A. Yariv, Matrix analysis of microring coupled-resonator optical waveguides, Opt. Express 12, 90 (2004).
  * Abdelshafy  _et al._ [2019] A. F. Abdelshafy, M. A. Othman, D. Oshmarin, A. T. Almutawa, and F. Capolino, Exceptional points of degeneracy in periodic coupled waveguides and the interplay of gain and radiation loss: Theoretical and experimental demonstration, IEEE Transactions on Antennas and Propagation 67, 6909 (2019).
  * Lagarias  _et al._ [1998] J. C. Lagarias, J. A. Reeds, M. H. Wright, and P. E. Wright, Convergence properties of the nelder–mead simplex method in low dimensions, SIAM Journal on optimization 9, 112 (1998).