# Circularly polarized states and propagating bound states in the continuum
in a periodic array of cylinders

Amgad Abdrabou  Corresponding author<EMAIL_ADDRESS>Department
of Mathematics, City University of Hong Kong, Kowloon, Hong Kong, China  Ya
Yan Lu  Department of Mathematics, City University of Hong Kong, Kowloon, Hong
Kong, China

###### Abstract

Bound states in the continuum (BICs) in a periodic structure sandwiched
between two homogeneous media have interesting properties and useful
applications in photonics. The topological nature of BICs was previously
revealed based on a topological charge related to the far-field polarization
vector of the surrounding resonant states. Recently, it was established that
when a symmetry-protected BIC (with a nonzero topological charge) is destroyed
by a generic symmetry-breaking perturbation, a pair of circularly polarized
resonant states (CPSs) emerge and the net topological charge is conserved. A
periodic structure can also support propagating BICs with a nonzero
wavevector. These BICs are not protected by symmetry in the sense of symmetry
mismatch, but they need symmetry for their robust existence. Based on a highly
accurate computational method for a periodic array of slightly noncircular
cylinders, we show that a propagating BIC is typically destroyed by a
structural perturbation that breaks only the in-plane inversion symmetry, and
when this happens, a pair of CPSs of opposite handedness emerge so that the
net topological charge is conserved. We also study the generation and
annihilation of CPSs when a structural parameter is varied. It is shown that
two CPSs with opposite topological charge and same handedness, connected to
two BICs or in a continuous branch from one BIC, may collapse and become a CPS
with a zero charge. Our study clarifies the important connection between
symmetry and topological charge conservation.

##  I Introduction

Bound states in the continuum (BICs) are trapped or guided modes with their
frequencies in the radiation continua neumann29 ; hsu16 . They exist in a
variety of photonic structures including periodic structures sandwiched
between two homogeneous media bonnet94 ; padd00 ; ochiai01 ; tikh02 ;
shipman03 ; lee12 ; port05 ; mari08 ; hsu13_2 ; bulg14b ; gan16 , waveguides
with local defects evans94 ; bulg08 , waveguides with lateral leakage channels
zou15 ; bezus18 ; nguyen19 ; yu19 ; byk20 , etc. In structures that are
invariant or periodic in one or two spatial directions, a BIC is a special
point with an infinite quality factor ($Q$ factor) in a band of resonant
states hsu13_2 ; yuan17_2 ; yuan18 ; jin19 , and it becomes a high-$Q$
resonance if the structure is perturbed generically kosh18 ; yuan20 . High-$Q$
resonances lead to strong local field enhancement mocella15 ; yoon15 ;
huzhen20  and abrupt features in reflection and transmission spectra yuan17 ,
and are essential for sensing, lasing, switching and nonlinear optics
applications.

For theoretical interest and practical applications, it is important to
understand how a BIC is affected by a perturbation of the structure. If the
BIC is protected by a symmetry bonnet94 ; padd00 ; ochiai01 ; tikh02 ;
shipman03 ; lee12 ; evans94 ; bulg08 , i.e., there is a symmetry mismatch
between the BIC and the compatible radiation modes, it naturally continues its
existence if the perturbation preserves the symmetry bonnet94 ; evans94 ;
shipman07 . Therefore, a symmetry-protected BIC is robust with respect to
symmetry-preserving perturbations. If a perturbation breaks the symmetry, a
symmetry-protected BIC typically, but not always, becomes a resonant state
with a finite $Q$ factor kosh18 ; yuan20 ; yuan20b . In periodic structures
sandwiched between two homogeneous media, there are also BICs with a nonzero
Bloch wavevector, and they propagate in the periodic directions port05 ;
mari08 ; hsu13_2 ; bulg14b ; yuan17 ; maksim ; hu18 . Such a propagating BIC
is not protected by symmetry in the usual sense, but can also be robust with
respect to symmetry-preserving perturbations zhen14 ; yuan17ol ; conrob . More
precisely, in a periodic structure with an up-down mirror symmetry and an in-
plane inversion symmetry, a generic low-frequency propagating BIC (with only
one radiation channel) continues its existence if the structure is perturbed
by a perturbation preserving these two symmetries yuan17ol ; conrob .

The BICs in periodic structures exhibit interesting topological properties.
Zhen el al. zhen14  realized that a BIC in a biperiodic structure is a
polarization singularity in momentum space (the plane of two wavevector
components), and defined a topological charge based on the winding number of
the far-field polarization vector. Since the far field of a resonant state is
typically elliptically polarized, Bulgakov and Makismov refined the definition
using the major polarization vector maksim . The topological charge can be
used to classify the BICs and illustrate their generation, interaction and
annihilation processes when structural parameters are tuned zhen14 ; maksim .
Importantly, the topological charge is a conserved quantity that cannot be
changed by small structural perturbations. However, this does not imply that
the BICs (with a nonzero topological charge) are robust with respect to
arbitrary structural perturbations, because a resonant state can be circularly
polarized and also has a nonzero topological charge. In a recent work, Liu el
al. fudan  showed that symmetry-protected BICs in a photonic crystal slab,
protected by the in-plane inversion symmetry, turn to pairs of circularly
polarized resonant states (CPSs) when the structure is perturbed breaking the
symmetry.

Most propagating BICs are found in periodic structures with both the up-down
mirror symmetry and the in-plane inversion symmetry. It is known that if one
of these two symmetries is broken, a propagating BIC is usually, but not
always, destroyed hu18 ; yuan20b . It has been shown that CPSs can exist in
structures without the up-down mirror symmetry and in-plane inversion symmetry
yin20 . In this paper, we show that CPSs emerge when propagating BICs are
destroyed by generic and arbitrarily small perturbations that break only the
in-plane inversion symmetry. We consider vectorial BICs in a periodic array of
of dielectric cylinders, introduce a small deformation of the cylinder
boundary, and show that pairs of CPSs emerge when propagating BICs with
topological charge $\pm 1$ are destroyed. In addition, we follow the CPSs as
the deformation parameter is further increased, and show that two CPSs with
opposite topological charges and same handedness may merge, and their
dependence on the deformation parameter may be multi-valued with self-
generation and annihilation points.

The rest of this paper is organized as follows. In Sec. II, we recall the
mathematical formulation for vectorial eigenmodes in two-dimensional (2D)
periodic structures and the definition of topological charge. In Sec. III, we
consider a periodic array of circular cylinders and briefly describe our
computational method. In Sec. IV, we introduce boundary deformations to the
cylinders and show the generation and annihilation of CPSs. The paper is
concluded with a brief discussion in Sec. V.

##  II Eigenmodes in periodic structures

We consider a 2D structure that is invariant in $x$, periodic in $y$ with
period $L$, bounded in $z$ by $\left|z\right|<D$ for some $D>0$, and
surrounded by air (for $|z|>D$), where $\\{x,y,z\\}$ is a Cartesian coordinate
system. Let $\varepsilon=\varepsilon(y,z)$ be the relative permittivity of
this structure and its surrounding medium, then $\varepsilon({\bf
r})=\varepsilon(y+L,z)$ for all ${\bf r}=(y,z)$ and $\varepsilon({\bf r})=1$
for $\left|z\right|>D$. We study time-harmonic electromagnetic waves that
depend on time $t$ and variable $x$ as $\exp[i(\alpha x-\omega t)]$, where
$\omega$ is the angular frequency and $\alpha$ is a real wavenumber in the $x$
direction. From the frequency-domain Maxwell’s equations, it is easy to obtain
the following system

$\displaystyle\nabla\cdot\left(\frac{\varepsilon}{\eta}\nabla
E_{x}\right)+\nabla\cdot\left(\frac{\alpha}{k\eta}P\cdot\nabla\widetilde{H}_{x}\right)+\varepsilon
E_{x}=0,\quad$ (1)
$\displaystyle\nabla\cdot\left(\frac{1}{\eta}\nabla\widetilde{H}_{x}\right)-\nabla\cdot\left(\frac{\alpha}{k\eta}P\cdot\nabla
E_{x}\right)+\widetilde{H}_{x}=0,$ (2)

where $E_{x}$ is the $x$ component of the electric field, $\widetilde{H}_{x}$
is the $x$ component of a scaled magnetic field (magnetic field multiplied by
free space impedance), $k=\omega/c$ is the free space wavenumber, $c$ is the
speed of light in vacuum, and

$\eta=k^{2}\varepsilon({\bf
r})-\alpha^{2},\quad\nabla=\begin{bmatrix}\partial_{y}\\\
\partial_{z}\end{bmatrix},\quad P=\begin{bmatrix}0&1\\\ -1&0\end{bmatrix}.$
(3)

The other four field components can be obtained from the following equations:

$\displaystyle\begin{bmatrix}E_{y}\\\
E_{z}\end{bmatrix}=\frac{i}{\eta}\left(\alpha\nabla
E_{x}+k\,P\cdot\nabla\widetilde{H}_{x}\right),$ (4)
$\displaystyle\begin{bmatrix}\widetilde{H}_{y}\\\
\widetilde{H}_{z}\end{bmatrix}=\frac{i}{\eta}\left(\alpha\nabla\widetilde{H}_{x}-k\,\varepsilon
P\cdot\nabla E_{x}\right).$ (5)

If $\alpha=0$, the equations for $E_{x}$ and $\widetilde{H}_{z}$ are
decoupled, we say the corresponding wave is scalar. The case $\alpha\neq 0$ is
referred to as vectorial.

Due to the periodicity in $y$, any eigenmode of the structure is a Bloch mode
with an electric field given by

${\bf E}({\bf r})={\bf F}({\bf r})e^{i\beta y},$ (6)

where $\beta$ is a real Bloch wavenumber satisfying
$\left|\beta\right|\leq\pi/L$, and ${\bf F}$ is a periodic in $y$ with period
$L$. For $\left|z\right|>D$, the field can be expanded in plane waves as

${\bf E}({\bf r})=\sum_{m=-\infty}^{\infty}{\bf
c}_{m}^{\pm}e^{i(\beta_{m}y\pm\gamma_{m}z)},\quad\pm z>D,$ (7)

where $\beta_{0}=\beta$, and

$\beta_{m}=\beta+\frac{2\pi
m}{L},\quad\gamma_{m}=\sqrt{k^{2}-\alpha^{2}-\beta_{m}^{2}}.$ (8)

A guided mode satisfies the condition ${\bf E}\to{\bf 0}$ as $z\to\pm\infty$.
For a positive $k$, guided modes usually exist when
$k<\sqrt{\alpha^{2}+\beta^{2}}$, so that all $\gamma_{m}$ are pure imaginary
and all plane waves in the right hand side of Eq. (7) are evanescent in $z$. A
BIC is also a guided mode, but it satisfies the condition
$k>\sqrt{\alpha^{2}+\beta^{2}}$, thus $\gamma_{0}$ and probably a few other
$\gamma_{m}$ are real positive. We study BICs with a real $k$, a real
$\alpha$, and a real $\beta$ such that

$\sqrt{\alpha^{2}+\beta^{2}}<k<\sqrt{\alpha^{2}+\left(\frac{2\pi}{L}-|\beta|\right)^{2}}.$
(9)

The above condition implies that $\gamma_{0}$ is positive and all other
$\gamma_{m}$ are pure imaginary. Since the plane waves
$\exp[i(\beta_{0}y\pm\gamma_{0}z)]$ can propagate to infinity, the
coefficients ${\bf c}_{0}^{\pm}$ of the BIC must vanish.

A resonant state is also an eigenmode, but it satisfies an outgoing radiation
condition as $z\to\pm\infty$. Because of the radiation loss, the amplitude of
a resonant state decays in time, thus the frequency $\omega$ or $k$ must be
complex with a negative imaginary part. This implies that the real and
imaginary parts of $\gamma_{0}$ are positive and negative, respectively, and
$\exp[i(\beta_{0}y+\gamma_{0}z)]$ is an amplifying outgoing plane wave as
$z\to+\infty$. We are concerned with resonant states with only one radiating
plane wave for either $z>D$ or $z<-D$. Therefore, it is assumed that
$\mbox{Im}(\gamma_{m})>0$ for all $m\neq 0$. Resonant states form bands where
each band corresponds to $k$ being a complex-valued function of $\alpha$ and
$\beta$. The $Q$ factor of a resonant mode is given by
$Q=-0.5\mathrm{Re}(k)/\mathrm{Im}(k)$. A BIC is a special point (with a real
$k$) in a band of resonant states.

We assume the periodic structure has an up-down mirror symmetry, i.e.,
$\varepsilon({\bf r})=\varepsilon(y,-z)$ for all ${\bf r}$, then it is
sufficient to study the field of any eigenmode in the upper half space
($z>0$), because the field components are either even in $z$ or odd in $z$.
For a BIC with a frequency and wavevector satisfying Eq. (9), it is possible
to define a topological charge based on the far field polarization vector of
the surrounding resonant states zhen14 ; maksim . Let $\mathcal{C}$ be a
closed contour in the $\alpha$-$\beta$ plane. Each point on $\mathcal{C}$
corresponds to a resonant state in a band that contains the BIC. The resonant
state contains a far-field outgoing plane wave (for $z\to+\infty$) with a
vector amplitude ${\bf c}_{0}^{+}$. Its projection on the $x$-$y$ plane is

$\widehat{\mathbf{E}}=\begin{bmatrix}c^{+}_{0x}\\\ c^{+}_{0y}\\\
0\end{bmatrix}e^{i(\alpha x+\beta y+\gamma_{0}z)},$ (10)

where $c^{+}_{0x}$ and $c^{+}_{0y}$ are the $x$ and $y$ components of ${\bf
c}_{0}^{+}$. For any fixed $z>D$, the real projected electric field
$\mbox{Re}(\widehat{\bf E}e^{-i\omega t})$ is typically elliptically
polarized, and the major polarization vector (along the major axis of the
polarization ellipse) forms an angle $\theta$ with the $x$ axis. If it is
possible to define $\theta$ continuously as $(\alpha,\beta)$ traverses along
$C$ in counterclockwise direction from a starting point back to the same point
(the ending point), and $\theta|_{\rm end}-\theta|_{\rm start}=2\pi q$ for
some $q$, where $\theta|_{\rm start}$ and $\theta|_{\rm end}$ are the values
of $\theta$ at the starting point and the ending point respectively, then $q$
is the winding number (of the projected major polarization vector) on
$\mathcal{C}$. Alternatively, $q$ can be evaluated by the integral formula

$q=\frac{1}{2\pi}\oint_{\mathcal{C}}d{\bm{\alpha}}\cdot\nabla_{\bm{\alpha}}\theta({\bm{\alpha}}),$
(11)

where ${\bm{\alpha}}=(\alpha,\beta)$ and
$\nabla_{\bm{\alpha}}=(\partial_{\alpha},\partial_{\beta})$. A BIC is a
polarization singularity, since it does not have a far field (${\bf
c}_{0}^{\pm}=0$). The topological charge of a BIC is defined as the winding
number $q$, if $\mathcal{C}$ is sufficiently close to the BIC and encloses the
BIC in the $\alpha$-$\beta$ plane zhen14 ; maksim .

It is important to note that the major polarization vector or the angle
$\theta$ is undefined, if the projected far-field plane wave is circularly
polarized. This implies that a circularly polarized state (CPS), i.e., a
resonant state with a circularly polarized far field, is also a polarization
singularity. If $\mathcal{C}$ contains the wavevector of a CPS, the winding
number on $\mathcal{C}$ is undefined. The definition of topological charge
requires the underlying assumption that no CPSs exist in a small neighborhood
(in the $\alpha$-$\beta$ plane) of the BIC.

Finally, we recall that the polarization state of the plane wave given by Eq.
(10) can be characterized by the Stokes parameters Stokes

$\displaystyle\mathcal{S}_{0}$
$\displaystyle=|c_{0x}^{+}|^{2}+|c^{+}_{0y}|^{2},\quad\mathcal{S}_{1}=|c_{0y}^{+}|^{2}-|c^{+}_{0x}|^{2},$
(12) $\displaystyle\mathcal{S}_{2}$
$\displaystyle=2|c^{+}_{0x}c^{+}_{0y}|\cos(\varphi_{y}-\varphi_{x}),$ (13)
$\displaystyle\mathcal{S}_{3}$
$\displaystyle=2|c^{+}_{0x}c^{+}_{0y}|\sin(\varphi_{y}-\varphi_{x}),$ (14)

where $\varphi_{x}=\arg(c^{+}_{0x})$, $\varphi_{y}=\arg(c^{+}_{0y})$ are the
phases of $c_{0x}^{+}$ and $c_{0y}^{+}$, respectively. Since a BIC has no
radiation, it corresponds to the point
$\left(\mathcal{S}_{0},\mathcal{S}_{1},\mathcal{S}_{2},\mathcal{S}_{3}=0\right)$.
A CPS satisfies the condition
$(\mathcal{S}_{1},\mathcal{S}_{2},\mathcal{S}_{3}/\mathcal{S}_{0})=(0,0,\pm
1)$, where the ellipticity $\mathcal{S}_{3}/\mathcal{S}_{0}=+1$ for the left
CPS and $-1$ for right CPS. A linearly polarized state is any point with
$\mathcal{S}_{3}=0$. The principle value of angle $\theta$ is given by
$\theta=\arg(\mathcal{S}_{1}+i\mathcal{S}_{2})/2$.

##  III Circular cylinders

A periodic array of circular dielectric cylinders surrounded by air, as shown
in Fig. 1,

Figure 1: A periodic array of circular dielectric cylinders surrounded by air.
The array is periodic in $y$ with period $L$. The cylinders are parallel to
the $x$ axis. The radius and dielectric constant of the cylinders are $a$ and
$\varepsilon_{c}$, respectively.

is a simple structure supporting many different BICs shipman03 ; shipman07 ;
bulg14b ; yuan17 ; maksim . In general, a BIC may be a standing wave
($\alpha=\beta=0$), may propagate along the $y$ axis ($\alpha=0$ and
$\beta\neq 0$), along the $x$ axis ($\alpha\neq 0$ and $\beta=0$), or in both
$x$ and $y$ directions ($\alpha\neq 0$ and $\beta\neq 0$). Those BICs with
$\alpha=0$ are scalar ones with either the $H$ or $E$ polarization. The BICs
with $\alpha\neq 0$ are vectorial ones with nonzero $E_{x}$ and
$\widetilde{H}_{x}$. In Table 1,

Table 1: A few BICs in periodic arrays of circular cylinders with radius $a$ and dielectric constant $\varepsilon_{c}=15$. | $a/L$ | $\alpha L$ | $\beta L$ | $kL$ | $q$
---|---|---|---|---|---
BIC1 | $0.45042$ | $0$ | $0$ | $3.512894$ | $-1$
BIC2 | $0.45042$ | $0.584168$ | $0$ | $3.523918$ | $+1$
BIC3 | $0.45$ | $0$ | $0.479401$ | $3.471973$ | $+1$
BIC4 | $0.3$ | $1.190401$ | $0.96714$ | $3.196556$ | $+1$

we list four BICs for periodic arrays with a fixed dielectric constant
$\varepsilon_{c}=15$. The topological charges of the BICs are listed in the
last column. More BICs in this periodic array can be found in Ref. maksim .

Normally, the resonant states and the BICs are computed by solving an
eigenvalue problem for the Maxwell’s equations. To take advantage of the
special geometry of the circular cylinders, we use a semi-analytic method
based on cylindrical wave expansions. Since the structure is invariant in $x$
and periodic in $y$, it is sufficient to solve the eigenmodes in one period of
the structure, i.e., a 2D domain $\Omega_{\infty}$ given by $|y|<L/2$ and
$|z|<\infty$. Inside $\Omega_{\infty}$, there is a square $\Omega$ given by
$|y|<L/2$ and $|z|<L/2$. We assume one cylinder is located at the center of
$\Omega$. For given $\alpha$ and $\omega$, the electromagnetic field in
$\Omega$ can be expanded in vectorial cylindrical waves with unknown
coefficients yumao . If $\beta$ is also specified, we can expand the
electromagnetic field in plane waves (also with unknown coefficients) for
$|z|>L/2$. Relating the field at $y=\pm L/2$ by the quasi-periodic condition,
assuming even or odd symmetry in $y$, and imposing continuity conditions at
$z=L/2$, we can obtain an operator $\mathscr{A}$, such that

$\mathscr{A}(k,\alpha,\beta)\,\mathbf{u}|_{z=L/2}={\bf 0},$ (15)

where ${\bf u}$ is a column vector for $E_{x}$ and $\widetilde{H}_{x}$, and
${\bf u}|_{z=L/2}$ denotes ${\bf u}$ at $z=L/2$ for $|y|<L/2$. Notice that
$\mathscr{A}$ is an operator that acts on a vector of two single-variable
functions. Since $\mathscr{A}$ depends on $k$, Eq. (15) is a nonlinear
eigenvalue problem. In practice, $y\in(-L/2,L/2)$ is sampled by $N$ points,
${\bf u}|_{z=L/2}$ becomes a vector of length $2N$, and $\mathscr{A}$ is
approximated by a $(2N)\times(2N)$ matrix. The method can be extended to the
case where the boundary of the cylinders are slightly and smoothly deformed.
Details are given in Appendix.

For computing resonant states, $\alpha$ and $\beta$ are given, we look for a
complex $k$ such that Eq. (15) has a nontrivial solution. One possible
approach is to solve $k$ from

$\lambda_{1}(\mathscr{A})=0,$ (16)

where $\lambda_{1}$ is the eigenvalue of $\mathscr{A}$ with the smallest
magnitude. A BIC is a special resonant state with $\mathrm{Im}(k)=0$. For
propagating BICs, it is more efficient to treat $\alpha$ and/or $\beta$ also
as unknowns. We can solve Eq. (16) for a real $k$, a real $\alpha$, and/or a
real $\beta$.

##  IV Slightly noncircular cylinders

In this section, we consider a periodic array of slightly noncircular
cylinders, study the emergence of CPSs when propagating BICs are destroyed,
and also the annihilation and generation of CPSs. We assume the boundary of
the cylinder centered at the origin is given by

$y=-a\sin(\tau)+\delta\cos(g\tau),\quad z=a\cos(\tau),$ (17)

for $0\leq\tau<2\pi$ and $g=2$ or $4$, where $\delta>0$ is a small deformation
parameter. The dielectric constant of the cylinders is fixed at
$\varepsilon_{c}=15$. The deformed cylinders for $\delta=0.1a$ are shown in
Figs. 2(a) and 2(b)

Figure 2: Cross-sections $\Omega_{1}$ of two deformed cylinders with a
boundary given by Eq. (17) for $\delta=0.1a$ and the cases (a) $g=2$ and (b)
$g=4$. $\Omega_{2}$ is the exterior domain outside the cylinder and inside the
square $\Omega=(-L/2,L/2)\times(-L/2,L/2)$.

for $g=2$ and $g=4$, respectively. The small boundary deformation is a
perturbation that breaks the reflection symmetry in $y$ (the in-plane
inversion symmetry), but preserves the reflection symmetry in $z$ (the up-down
mirror symmetry).

To calculate vectorial resonant states, we extend the method described in the
previous section. Importantly, a general electromagnetic field in $\Omega$
(the square domain containing one cylinder centered at the origin) can still
be expanded in cylindrical waves, although the expansions are more complicated
due to the deformation of the cylinder boundary. As shown in Appendix, the
eigenvalue problem for resonant states is reduced to Eq. (15), where
$\mathscr{A}$ is a $(2N)\times(2N)$ matrix depending on $k$, $\alpha$ and
$\beta$, and $N$ is the number of sampling points for an interval of length
$L$. It is highly desirable to compute a CPS without calculating all nearby
resonant states. To find a left or right CPS, we solve a real $\alpha$, a real
$\beta$ and a complex $k$, from Eq. (16) and

$\mathcal{S}_{1}=0,\quad\mathcal{S}_{3}/\mathcal{S}_{0}=\pm 1,$ (18)

where $\mathcal{S}_{0}$, $\mathcal{S}_{1}$ and $\mathcal{S}_{3}$ are the
Stokes parameters.

First, we show that when the deformation parameter $\delta$ is increased from
zero, all four BICs listed in Table 1 are destroyed and pairs of CPSs emerge.
If the topological charge of the BIC is 1 (or $-1$), a pair of CPSs with
topological charge $1/2$ (or $-1/2$) and different handedness emerge. The net
topological charge is conserved. In Fig. 3,

Figure 3: A pair of CPSs with topological charge $1/2$ emerged from BIC3 for a
periodic array of slightly noncircular cylinders with $g=4$.

we show the emergence of CPSs from BIC3 for a periodic array with $a=0.45L$
and $g=4$. The purple arrows indicate the direction of increasing $\delta$.
The curve shows the wavevector $(\alpha,\beta)$ for a pair of left and right
CPSs, and it is shown in color to indicate the value of $\delta$. These two
CPSs exhibit a symmetry with respect to the $\beta$ axis. For each left CPS
with wavevector $(\alpha,\beta)$, there is also a right CPS with wavevector
$(-\alpha,\beta)$. The complex frequencies of these two CPSs are exactly the
same. In Figs. 4(a) and 4(b),

Figure 4: Pairs of CPSs emerged from BIC4 for periodic arrays of slightly
noncircular cylinders with (a) $g=2$ and (b) $g=4$.

we show CPSs emerged from BIC4 for periodic arrays with $g=2$ and $g=4$,
respectively. The radius of the original cylinders is $a=0.3L$. There is no
apparent symmetry between the two CPSs emerged from the BIC, but there is
still a symmetry with respect to the $\beta$ axis. Since the structure is
invariant in $x$, the resonant states (or BICs) for $(\alpha,\beta)$ and
$(-\alpha,\beta)$ are reflections of each other. The CPSs with wavevectors
$(\alpha,\beta)$ and $(-\alpha,\beta)$ have the opposite handedness.

Next, we show that as $\delta$ is increased, the CPSs emerged from BICs of
opposite topological charge may collapse to a CPS with a zero charge. BIC1 and
BIC2 in Table 1 are found in the same periodic array of circular cylinders
with radius $a=0.45042L$, but their topological charges are $-1$ and $+1$,
respectively. As $\delta$ is increased from zero, both BIC1 and BIC2 are
destroyed. In Figs. 5(a) and 5(b),

Figure 5: Emergence of CPSs from BIC1 and BIC2, and annihilation of CPSs of
opposite charges for periodic arrays of slightly noncircular cylinders with
(a) $g=2$ and (b) $g=4$.

we show the emergence of CPS pairs from BIC1 and BIC2 for periodic arrays of
deformed cylinders with $g=2$ and $g=4$, respectively. The CPSs emerged from
BIC1 and BIC2 carry topological charges $-1/2$ and $1/2$, respectively. As
$\delta$ reaches a critical value $\delta_{\rm AP}$ ($\delta_{\rm AP}\approx
7.586\times 10^{-5}L$ for $g=2$ and $\delta_{\rm AP}\approx 3.736\times
10^{-5}L$ for $g=4$), two CPSs with the same handedness, one from BIC1 and the
other from BIC2, collapse to a CPS with a zero topological charge. Since these
two CPSs cease to exist for $\delta>\delta_{\rm AP}$, we call $\delta_{\rm
AP}$ an annihilation point (AP).

Finally, we show that CPSs connected to a single BIC may encounter self-
generation and self-annihilation points as $\delta$ is increased. For a
periodic array with $a=0.45L$ and $g=2$, a continuous branch of CPSs emerging
from BIC3 is shown in Fig. 6(a).

Figure 6: (a) A continuous right CPS emerged from BIC3 for a periodic array
with $g=2$, with $\delta$ increased from $0$ to $\delta_{AP}$, decreased to
$\delta_{GP}$, and increased again. The topological charge of the CPS has
changed from $1/2$ to $-1/2$ and back to $1/2$ accordingly. (b) The $Q$
factor, (c) $\alpha$ and (d) $\beta$ of the right CPS as multivalued functions
of $\delta$. The blue and red curves correspond to topological charge $1/2$
and $-1/2$, respectively.

The curve in Fig. 6(a) depicts the wavevector $(\alpha,\beta)$ of a CPS when
$\delta$ is first increased from zero to $\delta_{\rm AP}\approx 0.0077L$,
then decreased to $\delta_{\rm GP}\approx 0.0070L$, finally increased again.
The topological charge of the CPS has changed from $1/2$ to $-1/2$, and back
to $1/2$ accordingly. In the process of an increasing $\delta$, the critical
value $\delta_{\rm AP}$ is annihilation point where two CPSs of opposite
topological charge and same handedness collapse to a CPS with a zero charge,
but these two CPSs are not connected to different BICs. In fact, the CPS with
topological charge $-1/2$ and another CPS with topological charge $1/2$ emerge
from a CPS of zero topological charge at $\delta_{\rm GP}$. We call
$\delta_{\rm AP}$ a self-annihilation point and $\delta_{\rm GP}$ a self-
generation point. However, annihilation and generation are relative terms
depending on how the structure is tuned. For example, $\delta_{\rm GP}$ can
also be regarded as a self-annihilation point, if the structure is tuned by a
decreasing $\delta$. It is clear that the continuous branch of CPSs emerged
from BIC3, as shown in Fig. 6(a), exhibits a multi-valued dependence on
$\delta$. In Figs. 6(b), (c) and (d), we respectively show the $Q$ factor, and
wavevector components $\alpha$ and $\beta$ as multivalued functions of
$\delta$. Despite the generation and annihilation of CPSs as $\delta$ is
varied, the net topological charge is conserved and remains at $1/2$.

##  V Conclusion

Many applications of BICs are realized in periodic structures sandwiched
between two homogeneous media. It is known that the existence and robustness
of BICs, including the propagating BICs in periodic structures, depend
crucially on symmetry hsu13_2 ; zhen14 . More precisely, it has been proved
that some propagating BICs are robust with respect to structural perturbations
that preserve the in-plane inversion symmetry and the up-down reflection
symmetry, even though these BICs do not have a symmetry mismatch with
compatible radiating waves yuan17ol ; conrob . If the perturbation breaks one
of these two symmetries, the propagating BICs are typically destroyed and
become resonant states with a finite $Q$ factor hu18 . The topological charge
of a BIC, defined using the polarization vector of the resonant states
surrounding the BIC in momentum space, is an interesting concept and is useful
for understanding the evolution of BICs as structural parameters are varied
zhen14 ; maksim . The definition of topological charge does not require the
in-plane inversion symmetry. When structural parameters are varied, the
topological charge is always conserved, but this does not imply that the BICs
are robust with respect to arbitrary structural perturbations, because a CPS
is also a polarization singularity in momentum space and it can have a nonzero
topological charge. Therefore, the conservation of topological charge is only
valid when both BICs and CPSs are included. For propagating BICs, the
connection between symmetry-breaking perturbations and the emergence of CPSs
has not been clearly established in existing literature. Using a periodic
array of slightly noncircular cylinders as an example, we show that pairs of
CPSs emerge when propagating BICs are destroyed by arbitrarily small
perturbations that break only the in-plane inversion symmetry. We also study
the generation and annihilation of CPSs when structural parameters are varied.
It is shown that pairs of CPSs of opposite topological charge can collapse at
special CPSs with a zero topological charge, but the net topological charge is
still conserved.

## Acknowledgments

The authors acknowledge support from the Research Grants Council of Hong Kong
Special Administrative Region, China (Grant No. CityU 11305518).

*

##  Appendix A Construction of matrix $\mathscr{A}$

To obtain a matrix $\mathscr{A}$ satisfying Eq. (16) for the case of a
periodic array of slightly noncircular cylinders, we need to find cylindrical
wave solutions that are valid inside and outside a single cylinder, determine
a matrix $\mathscr{C}$ mapping ${\bf u}$ to the normal derivative of ${\bf u}$
on $\partial\Omega$ (the boundary of $\Omega$), and finally construct the
matrix $\mathscr{A}$. As shown in Fig. 2, $\Omega_{1}$ is a subdomain of
$\Omega$ corresponding to the cross section of the cylinder,
$\Omega_{2}=\Omega\backslash\overline{\Omega}_{1}$, $\varepsilon({\bf
r})=\varepsilon_{1}$ or $\varepsilon_{2}$ for ${\bf r}\in\Omega_{1}$ or
$\Omega_{2}$, respectively. The boundary of $\Omega_{1}$ is $\Gamma$.

For given $k$ and $\alpha$ and an integer $p$, we construct a vectorial
cylindrical wave solution that depends on two arbitrary coefficients $c_{p}$
and $r_{p}$. The $x$ components of this solution are assumed to be

$\displaystyle E_{x}^{(p)}$ $\displaystyle=\begin{cases}\displaystyle\sum
a_{pq}\,\frac{J_{q}(\rho_{1}r)}{J_{q}(\rho_{1}a)}\,e^{iq\theta},&{\bf
r}\in\Omega_{1}\\\ \displaystyle
c_{p}\frac{J_{p}(\rho_{2}r)}{J_{p}(\rho_{2}a)}\,e^{ip\theta}+\sum
b_{pq}\frac{Y_{q}(\rho_{2}r)}{Y_{q}(\rho_{2}a)}\,e^{iq\theta},&{\bf
r}\in\Omega_{2}\end{cases}$ $\displaystyle\widetilde{H}_{x}^{(p)}$
$\displaystyle=\begin{cases}\displaystyle\sum
s_{pq}\,\frac{J_{q}(\rho_{1}r)}{J_{q}(\rho_{1}a)}\,e^{iq\theta},&{\bf
r}\in\Omega_{1}\\\ \displaystyle
r_{p}\frac{J_{p}(\rho_{2}r)}{J_{p}(\rho_{2}a)}\,e^{ip\theta}+\sum
t_{pq}\frac{Y_{q}(\rho_{2}r)}{Y_{q}(\rho_{2}a)}\,e^{iq\theta},&{\bf
r}\in\Omega_{2}\end{cases}$

where $r$ and $\theta$ are the polar coordinates of ${\bf r}$,
$\rho_{j}=(k^{2}\varepsilon_{j}-\alpha^{2})^{1/2}$, $a$ is the radius of the
circle to which $\Gamma$ is close, $J_{q}$ and $Y_{q}$ are the first and
second kinds of Bessel functions of order $p$, the sums are for $q$ from
$-\infty$ to $+\infty$, and $a_{pq}$, $b_{pq}$, $s_{pq}$, $t_{pq}$ (for all
$q$) are unknown coefficients.

Let $\nu=(\nu_{y},\nu_{z})$ be the outward unit vector normal to $\Gamma$, and
$\tau=(-\nu_{z},\nu_{y})$ be a unit vector tangential to $\Gamma$, then the
tangential field components of this cylindrical wave solution are

$\displaystyle\widetilde{H}_{\tau}^{(p)}=$
$\displaystyle\frac{i}{\eta}\left[\varepsilon\frac{\partial
E^{(p)}_{x}}{\partial\nu}+\frac{\gamma_{0}}{k}\frac{\widetilde{H}_{x}^{(p)}}{\partial\tau}\right],$
(19) $\displaystyle E_{\tau}^{(p)}=$
$\displaystyle\frac{i}{\eta}\left[\frac{\partial\widetilde{H}_{x}^{(p)}}{\partial\nu}-\frac{\gamma_{0}}{k}\frac{\partial
E_{x}^{(p)}}{\partial\tau}\right].$ (20)

We discretize $\Gamma$ using $M$ points (assuming $M$ is odd) and truncate the
sums by $|q|\leq(M-1)/2$. The continuity conditions of $E_{x}$,
$\widetilde{H}_{x}$, $E_{\tau}$ and $\widetilde{H}_{\tau}$ at those $M$ points
on $\Gamma$ give rise to the following linear system:

$\begin{bmatrix}A_{11}&A_{12}&{\bf 0}&{\bf 0}\cr{\bf 0}&{\bf
0}&A_{23}&A_{24}\cr A_{31}&A_{32}&A_{33}&A_{34}\cr
A_{41}&A_{42}&A_{43}&A_{44}\end{bmatrix}\begin{bmatrix}{\bf a}_{p}\cr{\bf
b}_{p}\cr{\bf s}_{p}\cr{\bf t}_{p}\end{bmatrix}=c_{m}\begin{bmatrix}{\bf
g}_{1}\cr{\bf 0}\cr{\bf g}_{3}\cr{\bf
g}_{4}\end{bmatrix}+r_{m}\begin{bmatrix}{\bf 0}\cr{\bf h}_{2}\cr{\bf
h}_{3}\cr{\bf h}_{4}\end{bmatrix},$

where $A_{jk}$ for $1\leq j,k\leq 4$ are $M\times M$ matrices, ${\bf g}_{j}$
and ${\bf h}_{j}$ for $1\leq j\leq 4$ are column vectors of length $M$, ${\bf
a}_{p}$, ${\bf b}_{p}$, ${\bf s}_{p}$ and ${\bf t}_{p}$ are also column
vectors of length $M$, ${\bf a}_{p}$ is the vector for all $a_{pq}$, etc.
Solving the above linear system, we obtain column vectors ${\bf f}_{jk}^{(p)}$
for $1\leq j,k\leq 2$, such that

${\bf b}_{p}=c_{p}{\bf f}_{11}^{(p)}+r_{p}{\bf f}_{12}^{(p)},\quad{\bf
t}_{p}=c_{p}{\bf f}_{21}^{(p)}+r_{p}{\bf f}_{22}^{(p)}.$ (21)

Therefore, the cylindrical wave solution can be written as ${\bf
u}_{p}=c_{p}{\bf u}_{p}^{(1)}+r_{p}{\bf u}_{p}^{(2)}$, where ${\bf
u}_{p}^{(1)}$ and ${\bf u}_{p}^{(2)}$ are completely determined.

The general field in $\Omega$ can be expanded in the above cylindrical waves
as

${\bf u}=\sum_{p=-\infty}^{\infty}{\bf
u}_{p}=\sum_{p=-\infty}^{\infty}(c_{p}{\bf u}_{p}^{(1)}+r_{p}{\bf
u}_{p}^{(2)}).$ (22)

If $\partial\Omega$ is sampled by $4N$ points, we can truncate the index $p$
in the sum above to $-2N\leq p\leq 2N-1$, and evaluate ${\bf u}$ and the
normal derivative of ${\bf u}$ at the $4N$ points on $\partial\Omega$. This
leads to $(8N)\times(8N)$ matrices $\mathscr{D}_{1}$ and $\mathscr{D}_{2}$
such that

${\bf u}|_{\partial\Omega}=\mathscr{D}_{1}\begin{bmatrix}{\bf c}\\\ {\bf
r}\end{bmatrix},\quad\left.\frac{\partial{\bf
u}}{\partial\nu}\right|_{\partial\Omega}=\mathscr{D}_{2}\begin{bmatrix}{\bf
c}\\\ {\bf r}\end{bmatrix},$ (23)

where ${\bf u}|_{\partial\Omega}$ and $\partial_{\nu}{\bf
u}|_{\partial\Omega}$ are column vectors of length $8N$, ${\bf c}$ and ${\bf
r}$ are column vectors of length $4N$ with entries $c_{p}$ and $r_{p}$,
respectively. For simplicity, $\partial_{\nu}$ is simply taken to be
$\partial_{z}$ or $\partial_{y}$ on the horizontal and vertical sides of
$\Omega$. Therefore, we have matrix
$\mathscr{C}=\mathscr{D}_{2}\mathscr{D}_{1}^{-1}$, such that

$\left.\frac{\partial{\bf
u}}{\partial\nu}\right|_{\partial\Omega}=\mathscr{C}{\bf
u}|_{\partial\Omega}.$ (24)

If $\beta$ is given, ${\bf u}$ and $\partial_{y}{\bf u}$ satisfies the
following quasi-periodic condition:

$\displaystyle{\bf u}|_{y=L/2}$ $\displaystyle=e^{i\beta L}{\bf u}|_{y=-L/2},$
(25) $\displaystyle\left.\frac{\partial{\bf u}}{\partial y}\right|_{y=L/2}$
$\displaystyle=e^{i\beta L}\left.\frac{\partial{\bf u}}{\partial
y}\right|_{y=-L/2}.$ (26)

Since the structure has the up-down mirror symmetry, we have either ${\bf
u}(y,z)={\bf u}(y,-z)$ or ${\bf u}(y,z)=-{\bf u}(y,-z)$. Substituting these
conditions into Eq. (24), we obtain a matrix $\mathscr{B}_{0}$ such that

$\left.\frac{\partial{\bf u}}{\partial z}\right|_{z=L/2}=\mathscr{B}_{0}{\bf
u}|_{z=L/2}.$ (27)

For $z>L/2$, the field can be expanded in plane waves as in Eq. (7). Using the
plane wave expansion for ${\bf u}$, evaluating ${\bf u}$ and $\partial_{z}{\bf
u}$ at $z=L/2$, and eliminating the unknown coefficients, we obtain a matrix
$\mathscr{B}_{1}$, such that

$\left.\frac{\partial{\bf u}}{\partial z}\right|_{z=L/2}=\mathscr{B}_{1}{\bf
u}|_{z=L/2}.$ (28)

This leads to Eq. (15) for $\mathscr{A}=\mathscr{B}_{1}-\mathscr{B}_{0}$.

## References

  * (1) J. von Neumann and E. Wigner, “Über merkwürdige diskrete Eigenwerte,” Phys. Z. 30, 465-467 (1929).
  * (2) C. W. Hsu, B. Zhen, A. D. Stone, J. D. Joannopoulos, and M. Soljačić, “Bound states in the continuum,” Nat. Rev. Mater. 1, 16048 (2016).
  * (3) A.-S. Bonnet-Bendhia and F. Starling, “Guided waves by electromagnetic gratings and nonuniqueness examples for the diffraction problem,” Math. Methods Appl. Sci. 17, 305-338 (1994).
  * (4) P. Paddon and J. F. Young, “Two-dimensional vector-coupled-mode theory for textured planar waveguides,” Phys. Rev. B 61, 2090-2101 (2000).
  * (5) T. Ochiai and K. Sakoda, “Dispersion relation and optical transmittance of a hexagonal photonic crystal slab,” Phys. Rev. B 63, 125107 (2001).
  * (6) S. G. Tikhodeev, A. L. Yablonskii, E. A Muljarov, N. A. Gippius, and T. Ishihara, “Quasi-guided modes and optical properties of photonic crystal slabs,” Phys. Rev. B 66, 045102 (2002).
  * (7) S. P. Shipman and S. Venakides, “Resonance and bound states in photonic crystal slabs,” SIAM J. Appl. Math. 64, 322-342 (2003).
  * (8) J. Lee, B. Zhen, S. L. Chua, W. Qiu, J. D. Joannopoulos, M. Soljačić, and O. Shapira, “Observation and differentiation of unique high-Q optical resonances near zero wave vector in macroscopic photonic crystal slabs,” Phys. Rev. Lett. 109, 067401 (2012).
  * (9) R. Porter and D. Evans, “Embedded Rayleigh-Bloch surface waves along periodic rectangular arrays,” Wave Motion 43, 29-50 (2005).
  * (10) D. C. Marinica, A. G. Borisov, and S. V. Shabanov, “Bound states in the continuum in photonics,” Phys. Rev. Lett. 100, 183902 (2008).
  * (11) C. W. Hsu, B. Zhen, J. Lee, S.-L. Chua, S. G. Johnson, J. D. Joannopoulos, and M. Soljačić, “Observation of trapped light within the radiation continuum,” Nature 499, 188–191 (2013).
  * (12) E. N. Bulgakov and A. F. Sadreev, “Bloch bound states in the radiation continuum in a periodic array of dielectric rods,” Phys. Rev. A 90, 053801 (2014).
  * (13) R. Gansch, S. Kalchmair, P. Genevet, T. Zederbauer, H. Detz, A. M. Andrews, W. Schrenk, F. Capasso, M. Lončar, and G. Strasser, “Measurement of bound states in the continuum by a detector embedded in a photonic crystal,” Light: Science & Applications 5, e16147 (2016).
  * (14) D. V. Evans, M. Levitin and D. Vassiliev, “Existence theorems for trapped modes,” J. Fluid Mech. 261, 21-31 (1994).
  * (15) E. N. Bulgakov and A. F. Sadreev, “Bound states in the continuum in photonic waveguides inspired by defects,” Phys. Rev. B 78, 075105 (2008).
  * (16) C.-L. Zou, J.-M. Cui, F.-W. Sun, X. Xiong, X.-B. Zou, Z.- F. Han, and G.-C. Guo, “Guiding light through optical bound states in the continuum for ultrahigh-$Q$ microresonantors,” Laser Photonics Rev. 9, 114-119 (2015).
  * (17) E. A. Bezus, D. A. Bykov, and L. L. Doskolovich, “Bound states in the continuum and high-$Q$ resonances supported by a dielectric ridge on a slab waveguide,” Photonics Research 6, 1084-1093 (2018).
  * (18) T. G. Nguyen, G. Ren, S. Schoenhardt, M. Knoerzer, A. Boes, and A. Mitchell, “Ridge resonance in silicon photonics harnessing bound states in the continuum,” Laser Photonics Rev. 13, 1900035 (2019).
  * (19) Z. Yu, X. Xi, J. Ma, H. K. Tsang, C.-L. Zou, and X. Sun, “Photonic integrated circuits with bound states in the continuum,” Optica 6, 1342-1348 (2019).
  * (20) D. A. Bykov, E. A. Bezus, and L. L. Doskolovich, “Bound states in the continuum and strong phase resonances in integrated Gires-Tournois interferometer,” Nanophotonics 9(1), 83-92 (2020).
  * (21) L. Yuan and Y. Y. Lu, “Strong resonances on periodic arrays of cylinders and optical bistability with weak incident waves,” Phys. Rev. A 95, 023834 (2017).
  * (22) L. Yuan and Y. Y. Lu, “Bound states in the continuum on periodic structures surrounded by strong resonances,” Phys. Rev. A 97, 043828 (2018).
  * (23) J. Jin, X. Yin, L. Ni, M. Soljacic, B. Zhen, and C. Peng, “Topologically enabled unltrahigh-$Q$ guided resonances robust to out-of-plane scattering,” Nature 574, 501-504 (2019).
  * (24) K. Koshelev, S. Lepeshov, M. Liu, A. Bogdanov, and Y. Kivshar, “Asymmetric metasurfaces with high-$Q$ resonances governed by bound states in the contonuum,” Phys. Rev. Lett. 121, 193903 (2018)
  * (25) L. Yuan and Y. Y. Lu, “Perturbation theories for symmetry-protected bound states in the continuum on two-dimensional periodic structures,” Phys. Rev. A 101, 043827 (2020).
  * (26) V. Mocella and S. Romano, “Giant field enhancement in photonic lattices,” Phys. Rev. B 92, 155117 (2015).
  * (27) J. W. Yoon, S. H. Song, and R. Magnusson, “Critical field enhancement of asymptotic optical bound states in the continuum,” Sci. Rep. 5, 18301 (2015).
  * (28) Z. Hu, L. Yuan, and Y. Y. Lu, “Resonant field enhancement near bound states in the continuum on periodic structures,” Phys. Rev. A 101, 043825 (2020).
  * (29) L. Yuan and Y. Y. Lu, “Propagating Bloch modes above the lightline on a periodic array of cylinders,” J. Phys. B: Atomic, Mol. and Opt. Phys. 50, 05LT01 (2017).
  * (30) S. Shipman and D. Volkov, “Guided modes in periodic slabs: existence and nonexistence,” SIAM J. Appl. Math. 67, 687–713 (2007).
  * (31) L. Yuan and Y. Y. Lu, “Parametric dependence of bound states in the continuum on periodic structures,” Phys. Rev. A 102, 033513 (2020).
  * (32) E. N. Bulgakov and D. N. Maksimov, “Bound states in the continuum and polarization singularities in periodic arrays of dielectric rods,” Phys. Rev. A 96, 063833 (2017).
  * (33) Z. Hu and Y. Y. Lu, “Resonances and bound states in the continuum on periodic arrays of slightly noncircular cylinders,” J. Phys. B: At. Mol. Opt. Phys. 51, 035402 (2018).
  * (34) B. Zhen, C. W. Hsu, L. Lu, A. D. Stone, and M. Soljačič, “Topological nature of optical bound states in the continuum,” Phys. Rev. Lett. 113, 257401 (2014).
  * (35) L. Yuan and Y. Y. Lu, “Bound states in the continuum on periodic structures: perturbation theory and robustness,” Opt. Lett. 42(21), 4490-4493 (2017).
  * (36) L. Yuan and Y. Y. Lu, “Conditional robustness of propagating bound states in the continuum on biperiodic structures,” arXiv preprint arXiv:2001.00832 (2020).
  * (37) W. Liu, B. Wang, Y. Zhang, J. Wang, M. Zhao, F. Guan, X. Liu, L. Shi, and J. Zi, “Circularly Polarized States Spawning from Bound States in the Continuum,” Phys. Rev. Lett. 123, 116104 (2019).
  * (38) T. Yoda and M. Notomi, “Generation and Annihilation of Topologically Protected Bound States in the Continuum and Circularly Polarized States by Symmetry Breaking,” Phys. Rev. Lett. 125, 053902 (2020).
  * (39) X. Yin, J. Jin, M. Soljačić, C. Peng, and B. Zhen, “Observation of topologically enabled unidirectional guided resonances,” Nature (London) 580, 467-471 (2020).
  * (40) W. H. McMaster, “Polarization and the Stokes parameters,” Am. J. Phys. 22(6), 351-362 (1954).
  * (41) Y. Wu and Y. Y. Lu, “Dirichlet-to-Neumann map method for analyzing periodic arrays of cylinders with oblique incident waves,” J. Opt. Soc. Am. B 26, 1442-1449 (2009).


# Using machine learning algorithms to determine the post-COVID state of a
person by his rhythmogram

Sergey Stasenko
Institute of Applied Physics of RAS,Russia
Lobachevsky University, Russia
<EMAIL_ADDRESS>&Andrey Kovalchuk
Institute of Applied Physics of RAS,Russia
<EMAIL_ADDRESS>
&Eremin Evgeny
Lobachevsky University, Russia
<EMAIL_ADDRESS>&Natalya Zarechnova
Privolzhsky District Medical Center, Russia
&Maria Tsirkova
Privolzhsky District Medical Center, Russia
&Sergey Permyakov
Lobachevsky University, Russia
&Sergey Parin
Lobachevsky University, Russia
<EMAIL_ADDRESS>
&Sofia Polevaya
Lobachevsky University, Russia
<EMAIL_ADDRESS>

###### Abstract

In this study we applyed machine-learning algorithms to determine the post-
COVID state of a person. During the study, a marker of the post-COVID state of
a person was found in the electrocardiogram data. We have shown that this
marker in the patient’s ECG signal can be used to diagnose a post-COVID state.

_K_ eywords Machine-learning algorithms $\cdot$ Electrocardiogram $\cdot$
Post-COVID state $\cdot$ Covid-19 $\cdot$ Data analysis

##  1 Introduction

The use of cardiac rhythmography to assess the functional state began in the
middle of the twentieth century [1]. However, for a number of decades, the
possibilities of the method were limited by stationary registration
conditions. As a rule, measurements were carried out before or after a
specific load, in a horizontal position, at rest [2]. Rare exceptions were
studies performed directly in the process of loading samples: in space and
sports medicine, using simulators [3, 4]. In the 21st century, with the
expansion of the arsenal of telemetric systems, it became possible to register
cardiorhythmograms in the contexts of various types of natural activity.

Methods and technologies have been developed that allow non-invasively,
remotely, without interference in human behavior, to assess its functional
state during sports competitions, public speaking, driving, simultaneous
translation, fire fighting and other extreme loads [5, 6, 7, 8, 9, 10].

The COVID-19 pandemic has actualized a complex of pathological conditions
associated with risks of cognitive impairment: severe acute respiratory
syndrome (SARS), chronic stress, multisystem inflammatory syndrome,
disseminated intravascular coagulation syndrome [11]. A significant proportion
of patients with COVID-19 have extremely low blood oxygen saturations, but,
remarkably, there are disproportionately few symptoms of cerebral, or "happy"
hypoxia. Oxygen starvation of the brain can provoke neurological disorders,
especially in areas of the brain that are very sensitive to hypoxia. Cognitive
dysfunction is promoted by ischemic or hypoxic lesions of the hippocampus,
basal ganglia, cerebellum, and impaired functional connections typical of
SARS. Specific changes in cognitive processes may be associated with specific
neurotropic manifestations of the activity of coronoviruses [12]. Properties
of SARS-CoV-2 as brain protein aggregation catalyst and accelerator contribute
to serious damage to the structure and function of the central nervous system
(CNS), including infections of immune macrophages, microglia or astrocytes,
severe encephalitis, toxic encephalopathy and severe acute demyelinating
lesions. The fMRI data indicate that the foci of destruction extend to the
complex of subcortical structures, capture the thalamus, basal ganglia, and
neocortical zones included in the limbic system. The neuroanatomical scheme of
brain lesions is in good agreement with the neuroarchitecture of dopaminergic
pathways. It is no coincidence that one of the manifestations of the action of
SARS-CoV-2 is a decrease in the activity of dopamine in the nigrostriatal
complex [13]. The mechanisms of retrograde or anterograde neuronal transport
ensure the migration of viruses to motor and sensory terminals, which can
provoke significant distortions in sensorimotor reactions of any level of
complexity. Thus, SARS-CoV-2 affects the main components of the neural
platform that supports key cognitive processes:

  1. 1.

Damage to exteroceptive and interoceptive sensory channels provokes a
violation of perceptual processes of mapping objective signals into signs
subjective information images and back afferent signals;

  2. 2.

Destruction of the thalamus provokes a violation of the process of
concentration of information resources on the most significant objects and
events, that is, selective attention;

  3. 3.

Hypoxic lesions of the hippocampus provoke disturbances in the dynamic memory
system and provoke distortions in the processes of preservation and
reproduction information images;

  4. 4.

Degradation of the dopaminergic system distorts evaluative functions, disrupts
locomotor processes, and provokes the reduction of emotions, including vital
ones, such as “pain”.

From patient reports, a repertoire of cognitive problems associated with
COVID-19 is known: difficulty concentrating, reduction in motor activity,
impaired coordination of movements, loss of smell and taste, decrease in
visceral sensitivity, a sharp reduction in evaluative functions and
motivations [14]. The researchers note that about 18% of patients who
underwent coronavirus infection noted a decrease in concentration and memory,
and the effect persisted for more than 3 years [15]. At the same time,
disorders cognitive processes are not associated with either the age of
patients or the severity of the disease. This type of disorder can manifest
itself both in the elderly patients who, in addition to age and the presence
of chronic diseases, had a severe infection, and in younger patients who had
COVID-19 in a fairly mild form [16, 17, 18, 19]. There are studies showing the
relationship between various inflammatory markers and cognitive impairment
[16]. It should be noted that in recent years, studies have increasingly
mentioned another important symptom that may indicate the manifestation of
COVID-19 – neurological disorders manifested in the form of delirium [17, 20].
There are references to a significant deterioration in the state of cognitive
functions of patients who have come out of a state of delirium [21].

Studying the impact of COVID-19 on the human body is one of the most pressing
problems in medicine at the moment. The complexity of diagnosing COVID-19 and
the impossibility operational monitoring of the condition outside the clinic
contribute to the deterioration of the epidemiological situation and increase
the load on the system health care in general. Due to the development mobile
devices - rhythm recorders heart and telecommunication technologies, it
becomes possible to carry out remote screening of dangerous conditions at
different phases of the disease. Screening capabilities are provided by
knowledge of COVID-specific regimes of regulatory systems.

Within the framework of cooperation between UNN named after N.I. Lobachevsky
(Nizhny Novgorod) and the Volga Regional Medical Center (POMC) monitoring of
the functional state of patients diagnosed with COVID-19 was carried out. The
key method was the technology of event-related heart rate telemetry [22] in
the process of various clinical events. On the 86% of the rhythmogram records
of patients in the red zone of the covid hospital, the same type of low-
amplitude fluctuations of RR-intervals of the r-shaped form, called
cardiospike, were revealed 6. The term "Cardiospike" was first introduced by
the authors in the work [23]. In the database of rhythmograms for 2015–2018
(Internet service UNN "Kognit", cogni-nn.ru) r-shaped spikes are present only
in 3% of records. This suggests that cardio spikes can be considered as
markers of COVID-specific modes of operation of the regulator systems.

##  2 Methods

The diagram of the system for recording cardiac rhythmogram and analisys is
shown in Figure 1.

Figure 1: The diagram of the system for recording cardiac rhythmogram and
analisys

The technology of event-related heart rate telemetry provides mobility (does
not limit the subject in mobility and allows you to reliably record a signal
on significant distance of its source from the receiving unit), no breaks in
recording, autonomy of measurement, resistance to exogenous interference, the
possibility of obtaining and processing data in on-line mode.

To achieve these indicators, we have chosen a sensor platform - ZephyrTM HxMTM
Smart - Zephyr BIO PACH BH3-M1 (HxM, Zephyr Technology), which includes a
microprocessor, a radio signal receiving and transmitting unit and miniature
ECG, acceleration and distance sensors.

The transmission of data packets to a Smart-phone or a personal computer is
carried out via the Bluetooth SPP -2.4 GHz channel with an interval of 1 s
(Figure1). Each packet contains a unique sensor platform identifier, the last
15 R-R intervals, the time is fixed relative to the start records. Temporary
accumulation and pre-processing of data takes place on the Smart Phone mobile
communication device with the Android operating system, after which the
processed data is transmitted via GSM channels to the network Internet to a
specialized server of the system. Export of experimental data for further
processing is carried out in TXT and CSV formats.

###  2.1 Data collection

The database on which the neural network algorithm was trained to recognize
the post-COVID state according to ECG data contains ECG data of both subjects
who had Covid-19 (their condition is documented in the discharge records of
the Clinical Hospital No. - any visible deviations in the state of health
(diseases that can cause extrasystole cannot occur without visible signs or
complaints from patients - such diseases have clearly defined chronic or acute
conditions that will be noticeable to the subject). At the moment, the
database contains ECG records of 18 patients who have had documented Covid-19
and more than 1000 records of healthy subjects.

The database does not contain patients suffering from extrasystole (the
diseases that cause it do not go unnoticed, so they would be in the epicrisis
if they were present), in connection with this, anomalies called "cardiospike"
were found in the ECG data of patients who underwent Covid-19 , which became a
marker for training a neural network algorithm to recognize similar anomalies
in other ECG data.

All ECG data are presented in the COGNITOM Web platform (cogni-nn.ru). The
uploaded data from the system, on which the neural network algorithm was
trained to recognize the post-COVID state according to ECG data, is attached
as a zip archive. There is a CSV file in the archive, in which all covid
records are collected, in which there are 5 columns:

  * •

end-to-end identifier;

  * •

RR interval value;

  * •

markup, 0 - no maximum spike, 1 - maximum spike;

  * •

measurement time in milliseconds of a certain covid record.

###  2.2 Machine-learning algorithms

To solve the problem of segmenting the flow of RR intervals, a dilated
convolution network was chosen. Convolutional networks have a higher
performance and learning speed than recurrent networks (RNNs) that are
commonly used in time series detection and segmentation problems.
Convolutional networks, on the other hand, require more convolutions to get a
sufficient receptive field, i.e. capture long enough time sequences. To
eliminate this shortcoming, stretched convolutions are used, an example for 3
layers of which is shown in Figure 2.

Figure 2: Dilated convolution with kernel=3 and its receptive field.

The stretch for each subsequent layer can be defined as $d=k^{n-1}$, where $n$
is the number of the convolutional layer. As can be seen from the scheme is
shown in Figure 2, the receptive field is 27 and is defined as
$r=(k-1)\sum_{i=1}^{L}k^{i-1}+1$, where $L$ is the number of layers. Of
course, such a multilayer scheme can be equivalent to one convolutional layer
with a kernel of 27, but the calculation speed of the multilayer version is
higher.

Figure 3: The Redisual Block.

As shown in [24], the speed of the convolutional layer can also be increased
by dividing the convolution into horizontal channel-by-channel convolutions
and vertical inter-channel mixing. This scheme is borrowed for the base block
of the RR interval segmenter. The Redisual Block (see Figure 3) is a variant
of the MobileNetv3 base block adapted for time series. It retains the basic
principle of expanding the feature space from size (T, C) to (T, H) (Expansion
Layer), per-channel convolution with gaps, non-linearity based on GELU [25]
projection of the input (T, H) based on approach Squeeze-and-Excitation [26],
as well as a layer of inter-channel mixing to the original size (T, C)
(Compression Layer), to the output of which the input (Residual) is mixed. To
speed up convergence, parameterized dimension reset (Ts, S) (Skip) was used,
which allows more deep layers to be involved in learning at the earliest stage
[27]. It is worth noting that the temporal dimension of the Skip output is
truncated on the left and right $Ts=T-2P$, where $P$ is the size of the
padding. The value of P is determined by the truncation of the receptive field
of the extreme elements, caused by the addition of convolutions with boundary
values to preserve the dimension, as well as the maximum duration of the
desired object. Figure 4 shows the schematic of the RR interval detector. It
consists of F filters, each with L base blocks connected in series. The Skip
branches of each of the layers are summed up and fed into the final segmenting
unit. The diagram of this block is shown in Figure 5. The block converts the
dimension input (Ts, S) into the output (Ts, M), where M is the number of
classes (in the work M = 1).

Figure 4: Spike detector circuit. Figure 5: Head Block.

##  3 Results

The heart in conjunction with the regulation system rhythm can be considered a
phase sensor of the functional state of the organism [28]. The optimal signal
for displaying the operating modes of the physiological system is a
rhythmogram or a diagram of sequentially recorded time intervals between
myocardial ventricular systoles (QRS complexes). The wide diagnostic
capabilities of rhythmograms make it possible to objectify the functional
state in various fields, including sports medicine and psychophysiology. In
this regard, for patients with COVID-19 during treatment and rehabilitation in
the clinic, in addition to the mandatory set of examinations and procedures,
rhythmograms are registered with subsequent processing on the Cognite web
platform. An experimental database of records of 74 patients with a total size
of 110330 samples was collected. A typical rhythmogram of a patient with acute
respiratory disease, which COVID-19 is characterized by rigidity of RR
intervals [29, 30]. However, only for patients with an established diagnosis
of COVID-19 in the signal from the ZephyrSmart sensor platform, there are
special low-amplitude anomalies - cardio spikes (Figure 7).

Figure 6: An example of an ensemble of normalized spikes of various subjects
diagnosed with COVID-19 Figure 7: Section of the rhythmogram of a patient with
COVID-19 in 400 counts, containing 9 episodes of spike anomalies. Anomalies
have a characteristic repeating pattern and different amplitudes in the range
of ±100 ms from the average value

The spike pattern distinguishes two successive (by RR number) jumps from the
average value: a longer RR is followed by a shorter one with further slight
relaxation (see Figure 6). The data obtained indicate that cardiospike are of
an endogenous nature and can be considered as markers of COVID-specific modes
of heart rhythm regulation.

To train the model, we used labeled records of RR intervals, which are
sequences of delays between adjacent QRS impulses on a cardiogram. The records
are divided into segments of length $T$ with an overlap of $2P$, the response
of the $Ts$ model corresponds to the central part of the segment, which is
$T=P+Ts+P$. On the records of RR intervals, the maximum values of the desired
spike were previously marked, which are the goal of the study and
unambiguously describe the position and size spike. The loss function was
chosen based on the well-known FocalLoss [31]. The AdamW function [32] was
used as an optimizer.

Sample length T=32, P=4 was used for experiments. The choice of parameters (L,
F, C, H, S) was based on the analysis of cross-validation experiments when
dividing the base into 10 subsets and testing on one part, when the rest are
used for training. The model with the number of base channels L=32, hidden
channels H=40, side channels S=72, the number of layers in the block L=4 and
blocks F=2 became the optimal parameters in terms of speed and quality. At the
same time, the average quality on the test sets was at least 89% according to
F-score [33].

##  4 Conclusion and Discussion

Due to high-precision cardiointervalography based on heart rate telemetry,
specific rhythm anomalies – cardiospikes – were detected in the rhythmograms
of COVID-19 patients. The substantiation of the possibility to consider
cardiospike as an effective marker of violations of regulatory processes in
different phases of COVID-19 is given. A method for recording a post-state
based on the detection of cardio spikes in electrocardiogram records is
proposed.

##  5 Acknowledgements

The work in terms of data preprocessing was supported by project
0729-2021-013, which is carried out within the framework of the State task for
the performance of research work by laboratories that have passed the
competitive selection within the framework of the national project "Science
and Universities", in respect of which the decision of the Budget Commission
of the Ministry of Education and Science of Russia (dated September 14, 2021
No. BC-P/23) on the provision of subsidies from the federal budget for the
financial support of the state task for the implementation of research work
and in term of data analisys was supported by the Ministry of Education and
Science of Russian Federation under project No. 075-15-2021-634.

## References

  * [1] Parin V.V., Baevsky P.M. Introduction to medical cybernetics. Medicine. Moscow pp. 220 (1966)
  * [2] Kaznacheev V.P., Baevsky R.M. , Berseneva A.P., Domakhina G.M. , Polyakov Ya.V. On some features of the adaptation of the organism in connection with the profession and age. Labor hygiene and occupational diseases. 2 pp. 21-26 (1978)
  * [3] Grigoriev A.I., Bayevsky R.M. The concept of health and the problem of the norm in space medicine. Slovo. Moscow pp. 95 (2001)
  * [4] Shlyk N.I., Sapozhnikova E.N., Kirillova E.N., Semenov V.G. Typological features of the functional state of regulatory systems in schoolchildren and young athletes (according to the analysis of heart rate variability). Human Physiology. 35 pp. 85-93 (2009)
  * [5] Nekrasov M.M., Karatushina D.I., Parin S.B., Polevaya S.A. Application of information technologies for assessment of professional risks for high-altitude assemblers during periodic medical examination. Medical almanac. 3 pp. 26-31 (2011)
  * [6] Parin S.B., Vetyugov V.V., Bakhchina A.V., Polevaya S.A. The role of the endogenous opioid system in management heart rate variability in the context of cognitive loads different levels. The role of the endogenous opioid system in management heart rate variability in the context of cognitive loads different levels. 6 pp. 116-126 (2014)
  * [7] Polevaya S.A., Parin S.B., Runova E.V., Nekrasova M.M., Fedotova I.V., Bakhchina A.V., Kovalchuk A.V., Shishalov I.S. Telemetric and information technologies for monitoring of the functional state of athletes. Modern technologies in medicine. 4 pp. 94-98 (2012)
  * [8] Runova E.V., Grigorieva V.N., Bakhchina A.V., Parin S.B., Shishalov I.S., Kozhevnikov I.S., Nekrasov M.M., Karatushina D.I., Grigorieva D.I., Polevaya S.A. Vegetative correlates of arbitrary mappings emotional stress. International Journal of Psychophysiology. 4 pp. 69-77 (2013)
  * [9] Chernigovskaya T., Parin S., Parina I., Konina A., Urikh D., Yachmonina Yu., Chernova M., Polevaya S.A. Simultaneous interpreting and stress: pilot experiment. International Journal of Psychophysiology. 108 pp. 165 (2016)
  * [10] Chernigovskaya T.V., Parina I.S., Alekseeva S.V., Konina A.A., Urich D.K., Mansurova Yu.O., Parin S.B. Simultaneous interpreting and stress: pilot experiment. Modern technologies v medicine. 11 pp. 132-140 (2019)
  * [11] Miskowiak, K., Johnsen, S., Sattler, S., Nielsen, S., Kunalan, K., Rungby, J., Lapperre, T. & Porsberg, C. Cognitive impairments four months after COVID-19 hospital discharge: Pattern, severity and association with illness variables. European Neuropsychopharmacology. (2021)
  * [12] Beaud, V., Crottaz-Herbette, S., Dunet, V., Vaucher, J., Bernard-Valnet, R., Du Pasquier, R., Bart, P. & Clarke, S. Pattern of cognitive deficits in severe COVID-19. Journal Of Neurology, Neurosurgery And Psychiatry. (2021)
  * [13] Mukaetova-Ladinska, E., Kronenberg, G. & Raha-Chowdhury, R. COVID-19 and neurocognitive disorders. Current Opinion In Psychiatry. (2021)
  * [14] Devita, M., Bordignon, A., Sergi, G. & Coin, A. The psychological and cognitive impact of Covid-19 on individuals with neurocognitive impairments: research topics and remote intervention proposals. Aging Clinical And Experimental Research. (2021)
  * [15] Rogers, J., Chesney, E., Oliver, D., Pollak, T., McGuire, P., Fusar-Poli, P., Zandi, M., Lewis, G. & David, A. Psychiatric and neuropsychiatric presentations associated with severe coronavirus infections: a systematic review and meta-analysis with comparison to the COVID-19 pandemic. The Lancet Psychiatry. 7, 611-627 (2020,7), https://linkinghub.elsevier.com/retrieve/pii/S2215036620302030
  * [16] Zhou, H., Lu, S., Chen, J., Wei, N., Wang, D., Lyu, H., Shi, C. & Hu, S. The landscape of cognitive function in recovered COVID-19 patients. Journal Of Psychiatric Research. 129 pp. 98-102 (2020,10), https://linkinghub.elsevier.com/retrieve/pii/S0022395620308542
  * [17] Helms, J., Kremer, S., Merdji, H., Clere-Jehl, R., Schenck, M., Kummerlen, C., Collange, O., Boulay, C., Fafi-Kremer, S., Ohana, M., Anheim, M. & Meziani, F. Neurologic Features in Severe SARS-CoV-2 Infection. New England Journal Of Medicine. 382, 2268-2270 (2020,6), http://www.nejm.org/doi/10.1056/NEJMc2008597
  * [18] Chaumont, H., San-Galli, A., Martino, F., Couratier, C., Joguet, G., Carles, M., Roze, E. & Lannuzel, A. Mixed central and peripheral nervous system disorders in severe SARS-CoV-2 infection. Journal Of Neurology. (2020)
  * [19] Woo, M., Malsy, J., Pöttgen, J., Seddiq Zai, S., Ufer, F., Hadjilaou, A., Schmiedel, S., Addo, M., Gerloff, C., Heesen, C., Schulze Zur Wiesch, J. & Friese, M. Frequent neurocognitive deficits after recovery from mild COVID-19. Brain Communications. (2020)
  * [20] O’Hanlon, S. & Inouye, S. Delirium: A missing piece in the COVID-19 pandemic puzzle. Age And Ageing. (2020)
  * [21] Tsoi L.V. COVID-19-associated delirium: mechanisms of development and features currents. Medical Bulletin of the South of Russia. 12, 33-37 (2021), DOI 10.21886/2219-8075- 2021-12-1-33-37
  * [22] Polevaya, S., Eremin, E., Bulanov, N., Bakhchina, A., Kovalchuk, A. & Parin, S. Event-related telemetry of heart rate for personalized remote monitoring of cognitive functions and stress under conditions of everyday activity. Modern Technologies In Medicine. 11, 109-114 (2019)
  * [23] Permyakov S.A., Polevaya S.A., Tsirkova M.M., Eremin E.V., Parin S.B. Endogenous cardiac rhythm abnormalities in patients with COVID-19. Proceedings of the VII All-Russian Conference "Nonlinear Dynamics in Cognitive Research -2021", Nizhny Novgorod. https://nd-cogsci.ipfran.ru/2021/img/NDKI-2021.pdf
  * [24] Howard, A., Sandler, M., Chu, G., Chen, L., Chen, B., Tan, M., Wang, W., Zhu, Y., Pang, R., Vasudevan, V. & Others Searching for mobilenetv3. Proceedings Of The IEEE/CVF International Conference On Computer Vision. pp. 1314-1324 (2019)
  * [25] Hendrycks, D. & Gimpel, K. Gaussian error linear units (gelus). ArXiv Preprint ArXiv:1606.08415. (2016)
  * [26] Hu, J., Shen, L. & Sun, G. Squeeze-and-excitation networks. Proceedings Of The IEEE Conference On Computer Vision And Pattern Recognition. pp. 7132-7141 (2018)
  * [27] He, K., Zhang, X., Ren, S. & Sun, J. Deep residual learning for image recognition. Proceedings Of The IEEE Conference On Computer Vision And Pattern Recognition. pp. 770-778 (2016)
  * [28] Kuznetsov A.A. Processing methods ki and analysis of electrocardiographic information during prenosological research: textbook. Vladimir: Publishing House of VlGU. 237 (2012)
  * [29] Shirshov Yu.A., Govorin A.I. Vegetative disorders in patients with influenza A (H1N1). Siberian medical journal.5, 41–44 (2011)
  * [30] Zufarov A.A. Indicators of heart rate variability in acute respiratory syndrome in children. Young scientist.3, pp. 98–103 (2020), URL: https://moluch.ru/archive/293/66153/
  * [31] Lin, T., Goyal, P., Girshick, R., He, K. & Dollár, P. Focal loss for dense object detection. Proceedings Of The IEEE International Conference On Computer Vision. pp. 2980-2988 (2017)
  * [32] Loshchilov, I. & Hutter, F. Decoupled weight decay regularization. ArXiv Preprint ArXiv:1711.05101. (2017)
  * [33] Rijsbergen (CJ) Information retrieval. (Butterworth,1979)