# A direct demonstration of bulk-boundary correspondence in higher-order
topological superconductors with chiral symmetry

Xiaoyu Zhu  School of Physics, MOE Key Laboratory for Non-equilibrium
Synthesis and Modulation of Condensed Matter, Xi’an Jiaotong University,
Xi’an, Shaanxi 710049, China

###### Abstract

A higher-order topological superconductor can experience topological phase
transitions driven by variations in a bulk parameter without closing the bulk
gap. This presents a challenge in establishing a direct bulk-boundary
correspondence, as conventional bulk invariants change only upon the closure
of the bulk gap. Our study of two-dimensional higher-order phases in DIII and
BDI symmetry classes, both characterized by chiral symmetry, demonstrates that
Fermi level crossings facilitate a direct connection between the bulk
Hamiltonian and Majorana zero modes at corners. These crossings, emerging as
boundary conditions vary, can be identified from the bulk Hamiltonian. For
both classes, we introduce a pair of topological invariants derived from these
Fermi level crossings to characterize the higher-order topology. Phases in
which at least one invariant assumes a nonzero value are anticipated to host
Majorana corner modes. Moreover, these invariants may change with the closure
of either bulk or edge gaps, thereby providing a clear and direct
demonstration of bulk-boundary correspondence in higher-order phases. Our
findings offer a promising framework for systematically exploring higher-order
topology through boundary condition modulation.

##  I Introduction

Topology in superconducting systems is typically characterized by invariants
derived from bulk states under the assumption of periodic boundary condition
(PBC) [1, 2]. In most scenarios, a nontrivial topological invariant guarantees
the emergence of localized Majorana modes [3, 4, 5, 6, 7, 8, 9, 10] at the
boundaries (edges or surfaces), which are introduced by “cutting” the periodic
system. This phenomenon, known as bulk-boundary correspondence, is a defining
feature of topological states of matter. The recognition of higher-order
topology [11, 12, 13, 14, 15, 16, 17] extends this correspondence, leading to
the expectation that Majorana modes in superconductors may also appear at
intersections (corners or hinges) of adjacent boundaries. This implies that a
nontrivial higher-order topological superconductor [18, 19, 20, 21, 22, 23,
24, 25, 26] manifests Majorana modes when cut at least twice along different
directions. While this higher-order extension broadens the potential material
base within the topological family, it also presents challenges in fully
understanding the connection between bulk properties and boundary signatures.

A pertinent question is whether one may identify topologically protected
Majorana modes directly from the bulk Hamiltonian in higher-order topological
superconductors. The answer to this question becomes nuanced when considering
the influence of crystalline symmetries. It’s well known that crystalline
symmetry, such as mirror symmetry, can protect higher-order topology unless
the bulk gap closes [27]. States exhibiting this characteristic are termed
intrinsic higher-order topological states [28], resembling topological
crystalline states, wherein characterizing nontrivial topology from the bulk
states of a periodic system, using symmetry indicators for instance [29, 30,
31, 32, 33, 34, 35, 36, 37, 38, 39, 40], is plausible. However, higher-order
states also distinguish themselves from topological crystalline states in that
their boundary signatures can persist even under perturbations breaking
related crystalline symmetries [41, 42, 43, 44, 45]. In this circumstance,
variations in a bulk parameter may drive the system across topological phase
transitions without closing the bulk gap, making it a challenging task to
establish a direct bulk-boundary correspondence.

In our previous work [46] on two-dimensional (2D) superconductors in the D
symmetry class, we discovered that Fermi level crossings — zero-energy band
crossings in the energy spectrum — serve as robust indicators for higher-order
topology. These crossings manifest as the 2D system transitions from toroidal
to cylindrical geometry. Topological phases transitions, occurring where
either the bulk or edge gap closes, are accompanied with the emergence or
disappearance of Fermi level crossings. Importantly, these crossings could be
readily identified from bulk Hamiltonian, thereby acting as a bridge
connecting the bulk Hamiltonian and Majorana corner modes. In this study, we
focus on the DIII and BDI symmetry classes, and examine the higher-order bulk-
boundary correspondence within the same framework. In contrast to the D class,
the topological charges associated with these crossings — well defined due to
the presence of chiral symmetry in the two classes — also play a pivotal role
in determining the higher-order topological invariants.

The rest of this article is organized as follows: In Sec. II, we set up the
general framework and introduce the boundary-modulated Hamiltonian that is
central to our study. In Sec. III, we demonstrate how a pair of higher-order
topological invariants, which are derived from Fermi level crossings,
characterize the higher-order topology in the DIII class, and then take a
specific model as an example for illustration. In Sec. IV we investigate the
higher-order bulk-boundary correspondence in the BDI class following the same
procedure as in the DIII class. Finally, in Sec. V we provide additional
discussions and a brief summary to our work.

Figure 1: Illustration of the equivalence between higher-order topology in 2D
systems and first-order topology in quasi-1D systems. (a) Majorana zero modes,
indicated by the orange ovals, reside at two opposite corners in the 2D
system. In this configuration, the quasi-1D systems along the $x$ and $y$
direction both exhibit nontrivial first-order topology. (b) Majorana zero
modes appear at two adjacent corners. In this scenario, only the quasi-1D
system oriented along the $y$ direction demonstrates nontrivial first-order
topology. Lattice sites enclosed by blue (red) rectangles correspond to a
single unit cell in the quasi-1D system extending along the $x$ $(y)$
direction.

##  II Boundary-modulated Hamiltonian

We consider a generic 2D superconducting system with boundary conditions being
modulated in the $x$ or $y$ direction [31, 47]. The Hamiltonian is composed of
two parts, given by

$\mathcal{\mathcal{H}}_{a}(\lambda_{a})=\mathcal{H}-(1-\lambda_{a})\mathcal{B}_{a},$
(1)

where $a=x(y)$ and $\bar{a}=y(x)$. The first part, $\mathcal{H}$, represents
the bulk Hamiltonian, describing a system with periodic boundary condition in
both directions. The second part, $\mathcal{B}_{a}$, introduces boundary
terminations and is referred to as _boundary Hamiltonian_. $\mathcal{B}_{a}$
includes all inter-cell terms crossing the terminations (edges) that extend
along the $\bar{a}$ direction. It is important to note that this definition of
the boundary Hamiltonian does not involve intra-cell terms on boundary sites.
In this paper, we only consider boundary terminations that are commensurate
with unit cells. As such, the boundary Hamiltonian can be readily inferred
from the bulk Hamiltonian. The real parameter $\lambda_{a}$ controls the
boundary condition in the $a$ direction and modulates only the amplitude of
the boundary Hamiltonian. In particular, $\lambda_{a}=1(0)$ corresponds to
periodic (open) boundary condition, with $\mathcal{H}_{a}(\lambda_{a}=1(0))$
describing a toroidal (cylindrical) system. Varying $\lambda_{a}$ continuously
in the range $[0,1]$ is akin to gradually cutting the periodic system along
the $\bar{a}$ direction.

For the boundary-modulated Hamiltonian $\mathcal{H}_{a}$, the momentum along
the $\bar{a}$ direction remains a good quantum number, allowing us to express
its second-quantized form as

$\mathcal{H}_{a}(\lambda_{a})=\sum_{k_{\bar{a}}}\Psi_{a}^{\dagger}(k_{\bar{a}})H_{a}(k_{\bar{a}},\lambda_{a})\Psi_{a}(k_{\bar{a}}),$
(2)

where $\Psi_{a}(k_{\bar{a}})$ is the Nambu spinor, and
$H_{a}(k_{\bar{a}},\lambda_{a})$ is the Bloch Bogoliubov-de Gennes (BdG)
Hamiltonian, defined in the space parameterized by
$(k_{\bar{a}},\lambda_{a})$, with $k_{\bar{a}}\in[-\pi,\pi)$ and
$\lambda_{a}\in[0,1]$. We may alternatively view $\mathcal{H}_{a}$ as a series
of quasi-one dimensional (1D) Hamiltonian with varying $\lambda_{a}$, which
all extend along the $\bar{a}$ direction, as depicted in Fig. 1. The unit cell
for this quasi-1D system has a dimension proportional to the number of lattice
sites in the $a$ direction. This boundary-modulated Hamiltonian is central to
our study of higher-order topology in both symmetry classes. Unless otherwise
stated, we will omit the subscript of $k_{\bar{a}}$ and $\lambda_{a}$, with
the understanding that the Bloch Hamiltonian $H_{a}(k,\lambda)$ is defined in
the $(k_{\bar{a}},\lambda_{a})$ parameter space.

To determine the higher-order topology, we begin with the cylindrical
Hamiltonian $\mathcal{H}_{a}(\lambda=0)$, which features two edges stretching
along the $\bar{a}$ direction. When the 2D system supports Majorana corner
modes, distributed as in Fig. 1, the quasi-1D cylindrical system exactly
resembles a nontrivial first-order topological superconductor in 1D, with
stable Majorana zero modes at its ends. In this case, we may directly apply
the topological invariants used for diagnosing first-order topology to the
cylindrical system. When assuming nontrivial values, these invariants indicate
Majorana zero modes emerging at the corners, formed by cutting the cylindrical
system along the direction perpendicular to its edges. As the invariants are
derived from cylindrical Hamiltonian, they fail to reveal the connection
between Majorana corner modes and the 2D bulk (toroidal) Hamiltonian. To
explore the bulk-boundary correspondence in higher-order phases, we shall
investigate how the boundary-modulated Hamiltonian responds to the
continuously varying $\lambda$.

Figure 2: Fermi level crossings in nontrivial higher-order phases of three
symmetry classes. Higher-order topology in each class is characterized by a
pair of topological invariants. (a) In the D class, only the crossings (marked
by orange crosses) on high-symmetry lines contribute to the $\mathbb{Z}_{2}$
type topological invariants $(\mathcal{M}_{x},\mathcal{M}_{y})$. (b) In the
DIII class, Fermi level crossings occur in pairs at $(\pm k_{F},\lambda_{F})$,
with each pair possessing opposite topological charges, indicated by crosses
in different colors. The system with one or both of the two $\mathbb{Z}_{2}$
invariants $(\nu_{x},\nu_{y})$ taking nonzero values is expected to support
Majorana Kramers pairs at corners. (c) In the BDI class, Fermi level crossings
off high-symmetry lines appear in pairs at $(\pm k_{F},\lambda_{F})$. Unlike
in the DIII class, the two crossings in each pair may not have opposite
topological charges. The higher-order topological invariants
$(\omega_{x},\omega_{y})$ in this class are of $\mathbb{Z}$ type. In panels
(b) and (c), the colors of the dots indicate the chiralities ($\pm 1$) of
Majorana zero modes.

In our previous study of D-class superconductors, we established that Fermi
level crossings, which emerge during variations of $\lambda$, could be
utilized to identify nontrivial higher-order phases. For a D-class system, it
suffices to consider Fermi level crossings on the high-symmetry lines, namely
$K=0$ and $\pi$, in the $(k,\lambda)$ parameter space. The higher-order
topology is characterized by a pair of $\mathbb{Z}_{2}$ invariants, denoted as
$(\mathcal{M}_{x},\mathcal{M}_{y})$, which relates to the total count of Fermi
level crossings. When either of these invariants assumes the value $-1$, the
system is anticipated to harbor Majorana corner states, as depicted in Fig.
2(a). In this article, we focus on the DIII and BDI symmetry classes, and aim
to demonstrate that, Fermi level crossings emerging at any point in the
parameter space may significantly influence the higher-order topology in these
classes. The topological invariants are determined by the cumulative
topological charges of all crossings.

##  III DIII class

In the tenfold classification scheme [2], a DIII-class system is characterized
by three intrinsic symmetries, namely time-reversal ($T$), particle-hole ($P$)
and chiral symmetry ($S$). They act on the Hamiltonian $H_{a}(k,\lambda)$
according to

$\displaystyle TH_{a}(k,\lambda)T^{-1}=H_{a}(-k,\lambda),\ \ T^{2}=-1$
$\displaystyle PH_{a}(k,\lambda)P^{-1}=-H_{a}(-k,\lambda),\ \ P^{2}=1$
$\displaystyle SH_{a}(k,\lambda)S^{-1}=-H_{a}(k,\lambda),\ \ S^{2}=1$ (3)

where $T$ and $P$ are anti-unitary operators and $S$ is a unitary operator.
From Eq.(3) we observe that time-reversal and particle-hole transformations
act as reflections in the $(k,\lambda)$ parameter space, with the high-
symmetry lines at $K=0,\pi$ serving as the mirror lines. As a result, each
crossing at $(k_{F},\lambda_{F})$ has its partner at $(-k_{F},\lambda_{F})$.
Additionally, chiral symmetry enforces that zero energy levels at any point in
the parameter space are degenerate, with the degeneracy being $2n$, where $n$
represents the number of Fermi level crossings at that point. On the high-
symmetry lines, if Fermi level crossings exist, they are anticipated to come
in pairs due to Kramers degeneracy, thereby ensuring that the total count of
Fermi level crossings on these lines is always even. In D-class systems, the
parity of this number is closely related to the higher-order topological
invariant $\mathcal{M}_{a}$. The emergence of Majorana corner states is
expected when the parity is odd, i.e., $\mathcal{M}_{x}$ or $\mathcal{M}_{y}$
equals $-1$, as illustrated in Fig. 2(a). Apparently, these invariants always
take trivial values in DIII-class systems, and hence could not be used to
characterize the higher-order topology therein. In our subsequent discussions,
we introduce a pair of $\mathbb{Z}_{2}$ invariants, which also relate to Fermi
level crossings, to diagnose the nontrivial higher-order topology of the DIII
class.

###  III.1 $\mathbb{Z}_{2}$ invariants

Following the same procedure as in the D class, we relate the higher-order
topology in the 2D system to the first-order topology of the corresponding
cylindrical Hamiltonian $\mathcal{H}_{a}(0)$, which describes a quasi-1D
system with periodic boundary condition in the $\bar{a}$ direction. Due to
Kramers degeneracy, Majorana zero modes are always paired in real space,
forming Majorana Kramers pairs. When such pairs manifest at two of the four
corners in a square sample, as shown in Figs. 1 and 2(b), we can expect the
quasi-1D system, extending along the $x$ or $y$ direction, to exhibit
nontrivial first-order topology. In DIII class, the first-order topology is
characterized by topological invariants of $\mathbb{Z}_{2}$ type.

To further our analysis, we note that the Bloch BdG Hamiltonian
$H_{a}(k,\lambda)$ can be transformed into an off-diagonal form in the
eigenbasis of the chiral symmetry operator. This is represented by

$U_{S}^{\dagger}H_{a}(k,\lambda)U_{S}=\begin{pmatrix}0&D_{a}(k,\lambda)\\\
D_{a}^{\dagger}(k,\lambda)&0\end{pmatrix},$ (4)

where $U_{S}$ is the unitary matrix that diagonalizes chiral symmetry operator
$S$. For the DIII class, with an appropriate selection of $U_{S}$, the $D$
matrix in Eq.(4) can always be written in a form that satisfies
$D_{a}(k,\lambda)=-D_{a}^{T}(-k,\lambda)$ [48], implying that the $D$ matrix
on high-symmetry lines, $D_{a}(K,\lambda)$, is antisymmetric. Employing
singular value decomposition, we write
$D_{a}(k,\lambda)=U_{a}(k,\lambda)\Lambda_{a}(k,\lambda)V_{a}^{\dagger}(k,\lambda)$,
where $U_{a}$ and $V_{a}$ are both unitary matrices, and $\Lambda_{a}$ is a
diagonal matrix with nonnegative entries. We then define a unitary $q$ matrix,
$q_{a}(k,\lambda)=U_{a}(k,\lambda)V_{a}^{\dagger}(k,\lambda)$. The first-order
topology of the quasi-1D Hamiltonian $\mathcal{H}_{a}(\lambda)$ is
characterized by topological invariant $\tilde{\nu}_{a}(\lambda)$ in terms of
the $q$ matrix [49], which reads

$\displaystyle(-1)^{\tilde{\nu}_{a}(\lambda)}=$
$\displaystyle\frac{\text{Pf}[q_{a}(0,\lambda)]}{\text{Pf}[q_{a}(\pi,\lambda)]}\frac{\sqrt{\det[q_{a}(\pi,\lambda)]}}{\sqrt{\det[q_{a}(0,\lambda)]}}=\frac{\text{Pf}[q_{a}(0,\lambda)]}{\text{Pf}[q_{a}(\pi,\lambda)]}$
$\displaystyle\times$
$\displaystyle\exp\left\\{\frac{1}{2}\int_{0}^{\pi}dk\partial_{k}\ln\det[q_{a}(k,\lambda)]\right\\},$
(5)

where “Pf” represents Pfaffian. In Eq.(5), $\sqrt{\det[q_{a}(k,\lambda)]}$ is
required to be in the same branch for $k\in[0,\pi]$. Notably,
$\tilde{\nu}_{a}(0)$ and $\tilde{\nu}_{a}(1)$ represent the first-order
topological invariants of cylindrical and toroidal systems, respectively. We
denote $\nu_{a}=\tilde{\nu}_{a}(0)$, and the pair of invariant,
$(\nu_{x},\nu_{y})$, can precisely characterize nontrivial higher-order
phases. Specifically, when either $\nu_{x}$ or $\nu_{y}$ equals 1, the 2D
system is expected to host two Majorana Kramers pairs distributed at two
corners, as shown in Fig. 2(b). Alternatively, utilizing the two equalities,
$\text{Pf}(D_{a})=\sqrt{|\det(D_{a})|}\text{ Pf}(q_{a})$ and
$\det(D_{a})=|\det(D_{a})|\det(q_{a})$, the topological invariants can also be
expressed using the $D$ matrix, which reads

$\displaystyle(-1)^{\tilde{\nu}_{a}(\lambda)}=$
$\displaystyle\frac{\text{Pf}[D_{a}(0,\lambda)]}{\text{Pf}[D_{a}(\pi,\lambda)]}\frac{\sqrt{\det[D_{a}(\pi,\lambda)]}}{\sqrt{\det[D_{a}(0,\lambda)]}}=\frac{\text{Pf}[D_{a}(0,\lambda)]}{\text{Pf}[D_{a}(\pi,\lambda)]}$
$\displaystyle\times$
$\displaystyle\exp\left\\{\frac{1}{2}\int_{0}^{\pi}dk\partial_{k}\ln\det[D_{a}(k,\lambda)]\right\\}.$
(6)

Equations (5) and (6), both formulated using cylindrical Hamiltonian
$H_{a}(0)$, do not elucidate the link between corner modes and bulk
Hamiltonian. To achieve this, our subsequent analysis aims to connect the pair
of invariants to the Fermi level crossings that occur as the parameter
$\lambda$ varies. These crossings can be derived from bulk Hamiltonian of the
2D system, thus allowing us to establish a direct bulk-boundary correspondence
in higher-order phases.

Figure 3: Fermi level crossings on and off high-symmetry lines in DIII class.
(a) Variations of energy spectra with $\lambda$ for fixed $k$. On high-
symmetry lines (right panel), Kramers degeneracy requires all energy levels to
be degenerate and hence Fermi level crossings always come in pairs, whereas
there is no such restriction for the crossings off high-symmetry lines (left
panel). The red (solid) and blue (dashed) lines are deliberately shifted for
demonstration of the degeneracy. (b) The first case where the crossings only
occur off high-symmetry lines. (c) The second case with crossings on high-
symmetry lines. The loop $\mathcal{C}_{2}$ now excludes these crossings. The
red and blue crosses indicate Fermi level crossings with opposite topological
charges. A $\mathbb{Z}_{2}$ topological charge can be associated to each pair
of crossings, by considering a loop (black dashed lines in (b) and (c)) that
encloses the pair of crossings.

###  III.2 Fermi level crossings

Before delving into the relationship between Fermi level crossings and higher-
order topological invariants, we first take a look at the physical meaning of
$\tilde{\nu}_{a}(1)$, which is associated with the torus Hamiltonian
$\mathcal{H}_{a}(1)$. When interpreting the torus system as a quasi-1D system
along the $\bar{a}$ direction, a nontrivial $\tilde{\nu}_{a}(1)$ implies the
presence of Majorana Kramers pairs at both ends. This translates to the
existence of gapless Majorana modes propagating along the edges in the $a$
direction of the 2D system, a feature of first-order topological
superconductors in DIII class. However, our interest lies in topological
superconductors that exhibit nontrivial higher-order topology while
maintaining trivial first-order topology, characterized by gapped edge
spectra. Therefore, in our context, $\tilde{\nu}_{a}(1)$ is always assumed to
be zero. This suggests that, in the nontrivial higher-order phase, the
$\mathbb{Z}_{2}$ invariant $\tilde{\nu}_{a}(\lambda)$ is expected to change at
some point as $\lambda$ is continuously varied from 1 to 0. This transition
point is exactly where Fermi level crossings occur. In this regard, Fermi
level crossings reveal the topological distinctions between the toroidal and
the cylindrical system, which is crucial to our understanding of the bulk-
boundary correspondence in higher-order topological phases.

Combining the expression of $\tilde{\nu}_{a}(\lambda)$ for $\lambda=0$ and
$\lambda=1$ from Eq.(6), we arrive at an equivalent formula for the higher-
order topological invariant $\nu_{a}$, namely $\tilde{\nu}_{a}(0)$, which is
given by

$\displaystyle(-1)^{\nu_{a}}=$
$\displaystyle\frac{\text{Pf}[D_{a}(0,0)]}{\text{Pf}[D_{a}(0,1)]}\frac{\text{Pf}[D_{a}(\pi,1)]}{\text{Pf}[D_{a}(\pi,0)]}$
$\displaystyle\times$
$\displaystyle\exp\left\\{\frac{1}{2}\int_{0}^{\pi}dk\partial_{k}\ln\frac{\det[D_{a}(k,0)]}{\det[D_{a}(k,1)]}\right\\}.$
(7)

Our next step is to show that the ratio of Pfaffian terms in Eq.(7) can be
reformulated as an integral. It is important to note that while the $D$ matrix
on high-symmetry lines ($K=0$ or $\pi$), denoted as $D_{a}(K,\lambda)$, is
well defined for all $\lambda\in[0,1]$, its determinant may become zero.
Following this observation, our subsequent analysis will distinguish between
two cases.

In the first case, we assume $\det[D_{a}(K,\lambda)]\neq 0$ for all values of
$\lambda$, indicating that there are no Fermi level crossings on the high-
symmetry lines. We then select $N_{\lambda}$ equally spaced points along these
lines, with $\lambda_{n}=n/N_{\lambda}$ for integers $n$ in the range
$[0,N_{\lambda}]$. This allows us to express the ratio of Pfaffians in Eq.(7)
as the product of ratios for adjacent points, which reads

$\displaystyle\frac{\text{Pf}[D_{a}(K,1)]}{\text{Pf}[D_{a}(K,0)]}$
$\displaystyle=\prod_{n=0}^{N_{\lambda}-1}\frac{\text{Pf}[D_{a}(K,\lambda_{n+1})]}{\text{Pf}[D_{a}(K,\lambda_{n})]}$
(8)

Under the condition that $\det[D_{a}(K,\lambda)]\neq 0$, the phase of
$\text{Pf}[D_{a}(K,\lambda)]$ can be made to vary continuously with $\lambda$,
and so is the phase of $\sqrt{\det[D_{a}(K,\lambda)]}$. Due to this continuity
we have

$\displaystyle\frac{\text{Pf}[D_{a}(K,\lambda_{n+1})]}{\text{Pf}[D_{a}(K,\lambda_{n})]}=\frac{\sqrt{\det[D_{a}(K,\lambda_{n+1})]}}{\sqrt{\det[D_{a}(K,\lambda_{n})]}}.$
(9)

With this equality, we may further express Eq.(8) in an integral form in the
limit of infinitely large $N_{\lambda}$, given by

$\displaystyle\frac{\text{Pf}[D_{a}(K,1)]}{\text{Pf}[D_{a}(K,0)]}=\exp\left\\{\frac{1}{2}\int_{0}^{1}d\lambda\partial_{\lambda}\ln\det[D_{a}(K,\lambda)]\right\\}.$
(10)

From Eq.(10) it follows that the ratio of Pfaffian in the left side of Eq.(10)
is related to line integrals of the determinant of $D$ matrix along high-
symmetry lines, similar to the integral in Eq.(6), which is calculated along
lines with a constant $\lambda$ in the parameter space. Substituting Eq.(10)
into Eq.(7), we obtain

$\displaystyle(-1)^{\nu_{a}}=$
$\displaystyle\exp\left\\{\frac{1}{2}\oint_{\mathcal{C}_{1}}dl\nabla_{l}\ln\det[D_{a}(k,\lambda)]\right\\}$
$\displaystyle=$
$\displaystyle\exp\left[i\pi\sum_{(k_{F},\lambda_{F})}n_{a}(k_{F},\lambda_{F})\right],$
(11)

where the integration along the loop $\mathcal{C}_{1}$, illustrated in Fig.
3(a), represents the winding number of the $D$ matrix. In the second equality
of Eq.(11), we identify the winding number as the sum of topological charges
for all Fermi level crossings $(k_{F},\lambda_{F})$ enclosed by
$\mathcal{C}_{1}$, with $n_{a}(k_{F},\lambda_{F})$ representing the
corresponding charge. In practice, $n_{a}$ is determined by calculating the
winding number along a loop that encloses the crossing, such as the circle $l$
in Fig. 3(b). Time-reversal symmetry guarantees that each crossing at
$(k_{F},\lambda_{F})$ has its partner at $(-k_{F},\lambda_{F})$. Due to the
relation $D_{a}(k_{F},\lambda_{F})=-D_{a}(-k_{F},\lambda_{F})^{T}$, the two
partners have opposite charges, i.e.,
$n_{a}(k_{F},\lambda_{F})=-n_{a}(-k_{F},\lambda_{F})$. Thus, the total
topological charges of Fermi level crossings on either half of the parameter
space could give the same $\nu_{a}$.

Now let us consider the second case, where Fermi level crossings also occur on
high-symmetry lines. Considering that Pfaffian of the $D$ matrix becomes zero
at the crossing $(K,\lambda_{F})$, we need to exclude them while selecting the
sample points in Eq.(8). Also, Eq.(9) may fail if there exists crossings
between the two adjacent points, which we denote as
$(K,\lambda_{F}\pm\delta\lambda)$. The phase difference of the Pfaffian
between these two points is not necessarily infinitely small while
$\delta\lambda$ approaches zero. In this case, the loop $\mathcal{C}_{1}$ in
Eq.(11) is replaced by $\mathcal{C}_{2}$ shown in Fig. 3(b), which exactly
circumvents the crossings. All the sample points still lie on the high-
symmetry lines. On the vertical line segments of $\mathcal{C}_{2}$ ,
$\sqrt{\det[D_{a}(k,\lambda)]}$ can be made to vary continuously, thus
ensuring
$\sqrt{\det[D_{a}(K,\lambda_{F}+\delta\lambda)]}/\sqrt{\det[D_{a}(K,\lambda_{F}-\delta\lambda)]}$
to be infinitely small. This leads to the relation

$\displaystyle\frac{\text{Pf}[D_{a}(K,\lambda_{F}+\delta\lambda)]}{\text{Pf}[D_{a}(K,\lambda_{F}-\delta\lambda)]}$
$\displaystyle=(-1)^{\eta_{a}(K,\lambda_{F})}$ (12)
$\displaystyle\times\frac{\sqrt{\det[D_{a}(K,\lambda_{F}+\delta\lambda)]}}{\sqrt{\det[D_{a}(K,\lambda_{F}-\delta\lambda)]}},$

where $\eta_{a}(K,\lambda_{F})$ is an integer that depends on the crossing at
$(K,\lambda_{F})$. Comparing with Eq.(9), there is an additional phase factor
in Eq.(12). An odd $\eta_{a}$ suggests a $\pi$-phase shift in
$\text{Pf}[D_{a}(K,\lambda)]$ across the Fermi level crossings.

In order to determine the precise value of $\eta_{a}$, we decompose the
$2N\times 2N$ antisymmetric matrix $D_{a}(K,\lambda)$ into the form
$D_{a}=Q_{a}^{T}M_{a}Q_{a}$, with $Q_{a}$ being a unitary matrix, and $M_{a}$
taking the block diagonal form as follows,

$M_{a}=\text{diag}\left\\{\begin{pmatrix}0&m_{a,1}\\\
-m_{a,1}&0\end{pmatrix},...,\begin{pmatrix}0&m_{a,N}\\\
-m_{a,N}&0\end{pmatrix}\right\\}.$ (13)

For a non-singular $D$ matrix, $\pm m_{a,j}$, which represent eigenvalues of
Hamiltonian $\mathcal{H}_{a}(K,\lambda)$, all take finite values. Assuming
$|m_{a,1}|$ to be the smallest nonnegative energy level near the crossing
$(K,\lambda_{F})$, we then have $m_{a,1}=0$ exactly at the crossing. In the
case of linear band crossing, we have $m_{a,1}\propto(\lambda-\lambda_{F})$
near the crossing. As $\lambda$ crosses $\lambda_{F}$, a level crossing occurs
between $m_{a,1}$ and $-m_{a,1}$, as Fig. 3(a) shows in the energy spectrum,
and hence $m_{a,1}$ is expected to change its sign. We should note that, there
are two degenerate crossings in the energy spectrum due to Kramers degeneracy,
but $\pm m_{a,j}$ in Eq.(13) only account for half of the energy levels. While
$\lambda$ crosses $\lambda_{F}$, $Q_{a}$ — specifically, the phase of its
determinant — changes continuously. Alternatively, one may require $m_{a,1}$
to remain nonnegative, by swapping the positions of $\pm m_{a,1}$ in Eq.(13)
after crossing $\lambda_{F}$, albeit at the expense of altering the sign of
$\det[Q_{a}]$. Here, to clarify how Fermi level crossings affect the Pfaffian
of $D$ matrix, we require $\det[Q_{a}]$ to be a continuous function of
$\lambda$ near the crossing. By expressing the Pfaffian of $D$ matrix as

$\text{Pf}[D_{a}(K,\lambda)]=\det[Q_{a}(K,\lambda)]\prod_{j=1}^{N}m_{a,j}(K,\lambda),$
(14)

we readily find that the sign change of $m_{a,1}$ leads to an abrupt
$\pi$-phase shift in the Pfafian of $D$ matrix. At $\lambda_{F}$, it’s
possible for more than one pair of $\pm m_{j,a}$ to cross zero energy. The
number of such pairs is precisely half that of zero singular values of
$D_{a}(K,\lambda_{F})$ if we only consider linear band crossings. According to
Eq.(12), each pair contributes to a $\pi$-phase shift, and hence this number
is also equal to $\eta_{a}(K,\lambda_{F})$. Since only the parity of
$\eta_{a}$ matters, we may define the parity of $\eta_{a}$ to be the
$\mathbb{Z}_{2}$ charge of the crossing pair on the high-symmetry lines. In
this case, Eq.(8) still holds, as the sample points $\lambda_{n}$ do not
include the crossing points. For the two sample points adjacent to the
crossing, Eq.(9) needs to be replaced by Eq.(12). It then follows from Eq.(7)
that, the higher-order topological invariant in the second case is given by

$\displaystyle(-1)^{\nu_{a}}$
$\displaystyle=\prod_{(K,\lambda_{F})}(-1)^{\eta_{a}(K,\lambda_{F})}$
$\displaystyle\times\exp\left[i\pi\sum_{(k_{F}\neq
K,\lambda_{F})}n_{a}(k_{F},\lambda_{F})\right].$ (15)

Comparing with Eq.(11), Eq.(15) also includes contributions from Fermi level
crossings on high-symmetry lines. The summation in the second line of Eq.(15)
involves crossings in either half of the parameter space.

In both cases aforementioned, the pair of higher-order topological invariants
$(\nu_{x},\nu_{y})$ are determined by the topological charges of Fermi level
crossings in the parameter space. For the crossings off the high-symmetry
lines, their topological charges, $n_{a}$, are of $\mathbb{Z}$ type, in
contrast to the $\mathbb{Z}_{2}$ charge, namely $\eta_{a}$, for crossings on
the high-symmetry lines. Although $n_{a}$ may take any integer values, the
topological invariants only depend on its parity, as Eq.(11) and (15)
suggests. Besides, the time-reversal partner of each crossing actually doesn’t
provide additional information. With these observations, we may treat each
pair of Fermi level crossings at $(\pm k_{F},\lambda_{F})$ as a whole, and
associate a $\mathbb{Z}_{2}$ charge to them in a unified fashion, regardless
of whether they appear on or off the high-symmetry lines. The topological
charge, denoted as $\tilde{n}_{a}(k_{F},\lambda_{F})$, is defined according to
[48]

$(-1)^{\tilde{n}_{a}(k_{F},\lambda_{F})}=\frac{\text{Pf}[D_{a}(K,\lambda_{1})]}{\text{Pf}[D_{a}(K,\lambda_{2})]}\frac{\sqrt{\det[D_{a}(K,\lambda_{2})]}}{\sqrt{\det[D_{a}(K,\lambda_{1})]}},$
(16)

where the two points, $(K,\lambda_{1})$ and $(K,\lambda_{2})$, are the
intersections of the high-symmetry line and the dashed square loop, as
depicted in Figs. 3(b) and (c). The loop encloses only the pair of Fermi level
crossings considered. In Eq.(16), $\sqrt{\det[D_{a}(k,\lambda)]}$ is required
to vary continuously along the left or right half of the loops, instead on the
high-symmetry lines, where Fermi level crossings may appear. Equation (16)
reduces to Eq.(12) when considering an infinitely small loop and setting
$\lambda_{1}=\lambda+\delta\lambda$, $\lambda_{2}=\lambda-\delta\lambda$.
Winding numbers of the $D$ matrix along these loops are always zero, but the
$\mathbb{Z}_{2}$ charge defined in Eq.(16) can be nontrivial. To be specific,
$\tilde{n}_{a}(k_{F},\lambda_{F})=1$ if $\eta_{a}(K,\lambda_{F})$ (for
crossings on the high-symmetry lines) or $n_{a}(k_{F},\lambda_{F})$ (for
crossings off high-symmetry lines) is odd; otherwise
$\tilde{n}_{a}(k_{F},\lambda_{F})=0$. Consequently, the higher-order
topological invariant $\nu_{a}$ simply reflects the parity of the total
$\mathbb{Z}_{2}$ charges for all pairs of Fermi level crossings, i.e.

$\nu_{a}=\sum_{(k_{F}\geq
0,\lambda_{F})}\tilde{n}_{a}(k_{F},\lambda_{F})\text{ mod }2.$ (17)

Figure 4: The flow of Fermi level crossings in the parameter space under the
variations of $\delta m$ and $m_{1}$. (a) In the $(k_{y},\lambda_{x})$ space,
a pair of crossings, marked as black circles, reside on the high-symmetry line
$K=0$ when $m_{1}=0$. Finite $m_{1}$ separates the pair and drive them off the
line in opposite directions, as indicated by the arrows. With the increase of
$\delta m$, the crossings move towards $\lambda=0$ and vanish. (b)-(d)
$\lambda$ profiles generated by solving Eq.(24) for $k\in[-\pi,\pi)$. For
sufficiently large $\delta m$, the characteristic equation has no real root,
signaling the disappearance of Fermi level crossings. (e)-(h) A similar
phenomenon is observed in the $(k_{x},\lambda_{y})$ space. The other
parameters are set as $\mu=3t$, $\Delta=t$, $m=0.4t$, $\theta=\pi/4$, which,
unless otherwise specified, remain the same in subsequent plots.

###  III.3 Bulk-boundary correspondence

Having related the Fermi level crossings to higher-order topology, we now
demonstrate how these crossings and their topological charges are effectively
identified from the bulk Hamiltonian, so as to make a direct bulk-boundary
correspondence in higher-order phases. This approach is applicable to systems
in both DIII and BDI classes.

To begin with, we note that at the crossing point $(k_{F},\lambda_{F})$,
$\lambda_{F}$ is a root of characteristic equation $\det[H_{a}(k,\lambda)]=0$
[50]. Its multiplicity is even due to the presence of chiral symmetry. In
terms of $D$ matrix, the characteristic equation reduces to

$\det[D_{a}(k,\lambda)]=0,$ (18)

except that the root’s multiplicity is halved. In particular, for linear
crossings on the high-symmetry lines, we have
$\det[D_{a}(K,\lambda)]\propto(\lambda-\lambda_{F})^{p_{a}}$ near the
crossing, from which it follows that the root’s multiplicity $p_{a}$ is equal
to the number of zero singular values of $D_{a}(K,\lambda_{F})$ ($\pm m_{a,1}$
in Eq.(13)), which exactly equals to $2\eta_{a}(K,\lambda_{F})$, i.e.,
$p_{a}/2=\eta_{a}$. Hence, the root’s multiplicity determines the
$\mathbb{Z}_{2}$ charge of Fermi level crossings on the high-symmetry lines.
It’s possible that the crossings are of quadratic or cubic type, and there may
exist a number of crossings with different types at the same point
$\lambda_{F}$. We will discuss this generic case in Sec. V, where it is
demonstrated that the $\mathbb{Z}_{2}$ charge is still given by the parity of
$p_{a}/2$.

Following Eq.(11), the topological charge for each individual crossing off the
high-symmetry lines is equal to the winding number of the $D$ matrix along a
loop that encloses the crossing. No matter where the crossings occur, we need
to compute the determinant of the $D$ matrix. To demonstrate how it is
obtained from the bulk Hamiltonian, we first write down $H_{a}(k,\lambda)$ as

$H_{a}(k,\lambda)=H_{a}(k,1)-(1-\lambda)B_{a}(k),$ (19)

in the Bloch basis $\Psi_{a}(k)$, with $B_{a}(k)$ being the Bloch Hamiltonian
of boundary Hamiltonian $\mathcal{B}_{a}$ in Eq.(1). We denote the components
of $\Psi_{a}(k)$ as $[c_{a}(k)]_{i\alpha}$, where the Latin letter $i$ labels
the unit cell along the $a$ direction and the Greek letter $\alpha$ labels the
degree of freedom within each cell. The entries of the Bloch Hamiltonian
matrix are thus represented by $[H_{a}(k,\lambda)]_{i\alpha,j\beta}$. Since
chiral symmetry operates within the inner space, the corresponding $D$ matrix
has a similar relation as Eq.(19), which reads

$D_{a}(k,\lambda)=D_{a}(k,1)A_{a}(k,\lambda),$ (20)

with

$A_{a}(k_{,}\lambda)=\mathbb{I}-(1-\lambda)D_{a}^{-1}(k,1)D^{b}_{a}(k).$ (21)

Here, $\mathbb{I}$ represents the identity matrix, and $D^{b}_{a}(k)$ is the
$D$ matrix of $B_{a}(k)$.

For higher-order phases with gapped bulk spectrum, $\det[D_{a}(k,1)]$ always
takes a finite value, and hence Eq.(18) is equivalent to

$\det[A_{a}(k,\lambda)]=0.$ (22)

The boundary Hamiltonian $B_{a}(k)$ only involves boundary degrees of freedom,
which we denote by $i_{b}\alpha$. Consequently, the only possible nonzero
entries of $B_{a}(k)$ are $[B_{a}(k_{\bar{a}})]_{i_{b}\alpha,j_{b}\beta}$. As
such, we demonstrate in Appendix A that $\det[A_{a}(k,\lambda)]$ is equal to
the determinant of a $r_{a}\times r_{a}$ matrix, denoted by
$\tilde{A}_{a}(k,\lambda)$, with $r_{a}$ representing the rank of
$D_{a}^{b}(k)$. The specific form of $\tilde{A}_{a}(k,\lambda)$ is given by

$\displaystyle\tilde{A}_{a}(k_{\bar{a}},\lambda_{a})$
$\displaystyle=\mathbb{I}_{r_{a}\times r_{a}}-\frac{1}{2\pi}(1-\lambda_{a})$
$\displaystyle\times\int dk_{a}\tilde{V}_{a}^{b\dagger}[F(k_{a})\otimes
D^{-1}(\bm{k})]\bar{D}^{b}_{a}(k_{\bar{a}})\tilde{V}_{a}^{b},$ (23)

where we recover the subscripts of $k$ and $\lambda$ to avoid confusion. In
Eq.(23), $F(k_{a})$ is related to Fourier transformation with
$[F(k_{a})]_{i_{b},j_{b}}=e^{ik_{a}(i_{b}-j_{b})}$, $D(\bm{k})$ is the $D$
matrix of bulk Hamiltonian in 2D momentum space, and
$\bar{D}^{b}_{a}(k_{\bar{a}})$ is the block that corresponds to boundary cells
in matrix $D^{b}_{a}(k_{\bar{a}})$. $\tilde{V}_{a}^{b}$ is a semi-unitary
matrix related to compact singular value decomposition of
$\bar{D}^{b}_{a}(k_{\bar{a}})$.

As a result, Eq.(22) reduces to

$\displaystyle\det[\tilde{A}_{a}(k,\lambda)]=0.$ (24)

The loop integral in Eq.(15) may also be reexpressed in terms of
$\tilde{A}_{a}(k,\lambda)$, by noting that with trivial first-order topology
we always have

$\displaystyle\oint_{\mathcal{C}}dl\nabla_{l}\ln\det[D_{a}(k,1)]=0$ (25)

for any loop $\mathcal{C}$ in the parameter space. The topological charge of
each Fermi level crossing is then determined according to

$\displaystyle n_{a}(k_{F},\lambda_{F})=\frac{1}{2\pi
i}\oint_{\mathcal{C}}dl\nabla_{l}\ln\det[\tilde{A}_{a}(k,\lambda)],$ (26)

with $\mathcal{C}$ being any loop that encloses the crossing
$(k_{F},\lambda_{F})$ only. Therefore, all calculations involving the
determinant of $D_{a}(k,\lambda)$ can now be replaced by that of
$\tilde{A}_{a}(k,\lambda)$, which, according to Eq.(23), is obtained from the
$D$ matrix of bulk Hamiltonian, $D(\bm{k})$. After finding all the Fermi level
crossings and their topological charges, we immediately obtain the topological
invariants according to Eq.(15).

So we have established that Fermi level crossings serve as a bridge that
connects bulk Hamiltonian and Majorana corner modes. Although Fermi level
crossings are determined from the bulk Hamiltonian, their numbers may change
when bulk gap remains open. This happens when topological phases transitions
are driven by the closure of edge gaps, as we shall demonstrate in the
following example.

###  III.4 An example

We consider a spinful toy model in 2D with two orbital degrees in each unit
cell. The Bloch BdG Hamiltonian takes the form

$\displaystyle H_{\text{DIII}}(\bm{k})=\epsilon_{\bm{k}}\tau_{z}+2\Delta(\sin
k_{x}\tau_{y}-\sin k_{y}\tau_{x}s_{z}\sigma_{z})$ (27)
$\displaystyle+m(\sin\theta\tau_{z}\sigma_{x}-\cos\theta
s_{z}\sigma_{y})-\delta
m\tau_{x}s_{z}\sigma_{y}-m_{1}\tau_{x}s_{x}\sigma_{y},$

where the kinetic energy $\epsilon_{\bm{k}}=\mu-2t(\cos k_{x}+\cos k_{y})$,
$\Delta$ represents the $p$ wave pairing amplitude, and Pauli matrices $\tau$,
$s$ and $\sigma$ act in particle-hole, spin and orbital space, respectively.
The three symmetry operators defined in Eq.(3) are represented as
$T=-is_{y}\mathcal{K}$, $P=\tau_{x}\mathcal{K}$ and $S=\tau_{x}s_{y}$, where
$\mathcal{K}$ denotes the complex conjugation operator.

The first two terms in Eq.(27) describe two copies ($\sigma_{z}=\pm 1$) of
topological superconductors with $p\pm ip$ pairing. Each copy features two
counter-propagating helical Majorana modes on edges in the nontrivial phase.
The three onsite terms $m$, $\delta m$ and $m_{1}$ could gap out the helical
modes. In particular, the $m$ term introduces anisotropic mass gaps to the
four edges, controlled by $\theta$. This is crucial to the emergence of
Majorana corner modes, which appear when two adjacent edges exhibit mass gaps
of opposite signs. In the absence of $m_{1}$, $s_{z}$ is a good quantum
number, and the model simply reduces to two independent copies ($s_{z}=\pm 1$)
of the D-class model we introduced in Ref.[46]. In the D-class model, the
nontrivial higher-order topology is captured by Fermi level crossings on high-
symmetry lines. Crossings off these lines are not protected due to the absence
of chiral symmetry. In contrast, for the DIII class, Fermi level crossings
both on and off these lines contribute to the higher-order topology. Time-
reversal symmetry-preserving perturbations, such as the $m_{1}$ term, either
couple and eliminate the pair of crossings on high-symmetry lines or separate
them away from these lines, depending on the $\mathbb{Z}_{2}$ charge of the
pair defined in Eq.(16).

Figure 5: The $\pi$-phase jump in Pfaffian of the $D$ matrix induced by Fermi
level crossings. (a), (d) Fermi level crossings occurring on high-symmetry
lines when $m_{1}=0$ leads to an abrupt $\pi$-phase jump in
$\text{Pf}[D_{a}(0,\lambda)]$. For cases with finite $m_{1}$, there is no such
discontinuity. (b)-(c), (e)-(f) Energy spectra along $K=0$. No Fermi level
crossing exists on the high-symmetry line when $m_{1}$ is finite. In these
plots, $N_{a}$, the number of unit cells along the $a$ direction, is chosen to
be 40. The Pfaffian is calculated using the code from Ref.[51].

Given that we focus on boundary terminations commensurate with unit cells, the
boundary Hamiltonian is directly derived from the bulk Hamiltonian. In the
boundary basis, which in this model includes only the first and last unit
cells — denoted as $\\{\psi_{a,N_{a}}(k),\psi_{a,1}(k)\\}$ with $N_{a}$
representing the number of unit cells in the $a$ direction and $\psi_{a,j}(k)$
encompassing all internal degrees of freedom $[c_{a}(k)]_{j\alpha}$ — the
resulting matrix $B_{a}$ takes the following simple form,

$\displaystyle B_{a}=\begin{pmatrix}0&h_{a}\\\ h_{a}^{\dagger}&0\end{pmatrix}$
(28)

where $h_{a}$ includes all inter-cell hopping and pairing terms, with
$h_{x}=-t\tau_{z}-I\Delta\tau_{y}$ and
$h_{y}=-t\tau_{z}+I\Delta\tau_{x}s_{z}\sigma_{z}$. The unitary matrix that
sends the bulk Hamiltonian $H(\bm{k})$ into block off-diagonal form is

$\displaystyle U_{S}=\frac{1}{\sqrt{2}}\begin{pmatrix}s_{0}&-is_{y}\\\
s_{y}&is_{0}\end{pmatrix}\otimes\sigma_{0},$ (29)

where $s_{0}$ and $\sigma_{0}$ are identity matrices in the spin and orbital
space, respectively. The Fermi level crossings can then be identified
according to Eqs.(23) and (24).

In Figs. 4(a) and (e), we plot the positions of Fermi level crossings in the
parameter space for various values of $\delta m$ and $m_{1}$. When $m_{1}=0$,
the system supports a pair of Fermi level crossings on the high-symmetry line
$K=0$. The introduction of the $m_{1}$ term drives the pair of crossings away
from this line along opposite directions. The arrows in Figs. 4(a) and (e)
indicate where the Fermi level crossings flow with increasing $m_{1}$.

In this model, each Fermi level crossing at $(k_{F},\lambda_{F})$ is also
paired with another at $(k_{F},-\lambda_{F})$, the latter not explicitly shown
here as we focus only on crossings with $\lambda_{F}\in[0,1]$. As $\delta m$
increases (see Figs. 4(a) and (e)), the crossings at $\pm\lambda_{F}$ migrate
towards $\lambda=0$, eventually annihilating. A straightforward way to
illustrate this transition is through the $\lambda$ profile, which displays
all roots of Eq.(24) for $k\in[-\pi,\pi)$. Fermi level crossings correspond to
real roots for $\lambda$. As depicted in Figs. 4(b)-(d) and (f)-(h), the Fermi
level crossings vanish at sufficiently large $\delta m$, signaling edge phase
transitions where Majorana zero modes at different corners may either change
positions or couple. Our previous study in the D-class model [46] showed that,
in the absence of $m_{1}$, the energy gap on edges along the $x$ direction
closes at $|\delta m|=|m\sin\theta|$, and along the $y$ direction at $|\delta
m|=|m\cos\theta|$. The presence of $m_{1}$ term has a negligible impact on
these gap-closing points. Pinpointing the phase transition position for finite
$m_{1}$ is not our focus here. Instead our next step is to examine the
topological charges of these Fermi level crossings in order to derive the
higher-order topological invariants.

Figure 6: Topological charges of Fermi level crossings off high-symmetry
lines. Nontrivial windings in the phase of $\det[\tilde{A}_{a}(k,\lambda)]$
are observed along the loops $l_{1}$ and $l_{2}$, each of which encircles a
single Fermi level crossing. The two time-reversal partners at $\pm k_{F}$
carry opposite charges. In these plots, $m_{1}=0.3t$ and $\delta m=0$. Figure
7: Distributions of Majorana Kramers pairs in nontrivial higher-order phases.
(a) The probability distributions of Majorana zero modes. The inset shows the
energy spectrum for an $80\times 80$ open-boundary system, with red dots
representing Majorana zero modes. The two invariants, $\nu_{x}$ and $\nu_{y}$,
both take nontrivial values, as indicated by $\lambda$ profiles in (b) and
(c). This results in two Majorana Kramers pairs sitting at opposite corners.
(d) Only one of the invariants is nonzero, as revealed by $\lambda$ profiles
in (e) and (f), and the two Majorana pairs therefore reside on adjacent
corners. In these plots, $m_{1}=0.3t$ and $\delta m=0.1t$.

Let’s first consider Fermi level crossings on high-symmetry lines, denoted as
$(K,\lambda_{F})$. We already know that $\lambda_{F}$ is a root of Eq.(24) at
$k=K$. According to Eq.(12), $\eta_{a}(K,\lambda_{F})$, which equals to half
the root’s multiplicity, determines whether the Pfaffian of the $D$ matrix
experiences a $\pi$-phase jump as $\lambda$ crosses this point. The plots
corresponding to $m_{1}=0$ in Figs. 5 (a) and (c) indeed exhibit such a jump,
suggesting an odd $\eta_{a}$. By solving Eq.(24), we obtain a real root with
multiplicity 2, implying $\eta_{a}(K,\lambda_{F})=1$. This is consistent with
the energy spectrum along the high-symmetry line in Figs. 5(b) and (e), where
each level is doubly degenerate due to Kramers degeneracy, resulting in
exactly two linear band crossings at the zero energy for $m_{1}=0$. Turning on
$m_{1}$ drive the crossings away from the high-symmetry line, and hence the
Pfaffian of $D$ matrix varies continuously along the line, as is seen in Figs.
5(a) and (d) for finite $m_{1}$. This leads us to the first case investigated
in Sec. IIIB, where Fermi level crossings appear off the high-symmetry lines.

Since the pair of crossings on the high-symmetry line has a nontrivial
$\mathbb{Z}_{2}$ charge, we can expect finite $m_{1}$ doesn’t couple them. To
verify this argument, we may first locate the position of each crossing,
followed by the determination of their topological charges according to
Eq.(26). Note that each crossing at $(k_{F},\lambda_{F})$ has its time-
reversal partner at $(-k_{F},\lambda_{F})$, with opposite charges. In Figs.
6(b) and (c), we show the phase of $\det[\tilde{A}_{a}(k,\lambda)]$ along a
contour that encircles each crossing for $m_{1}=0.3$, which clearly shows
nontrivial windings, with topological charges being $\pm 1$. This suggests the
$\mathbb{Z}_{2}$ charge associated with the pair of crossings remains to be
$1$. Due to this nontrivial $\mathbb{Z}_{2}$ charge, finite $m_{1}$ can at
most separate the two crossings away from each other.

Through the analysis above, we are able to obtain the pair of topological
invariants $(\nu_{x},\nu_{y})$ as defined in Eq.(15). For the examples shown
in Figs. 5 and 6, we have $(\nu_{x},\nu_{y})=(1,1)$, which indicates that two
Majorana Kramers pairs sit at diagonals, as depicted in Fig. 7(a). Without
inversion symmetry, the two Majorana pairs may reside at any two corners. In
any of the cases, at least one of the invariants is nonzero. As showcased in
Fig. 7(d) for $\theta=0$, Fermi level crossings only occur in
$(k_{y},\lambda_{x})$ parameter space and hence $(\nu_{x},\nu_{y})=(1,0)$,
suggesting two Majorana Kramers pairs emerging at neighboring corners. It’s
important to note that when both the two invariants are zero, the system may
still reside in nontrivial higher-order phases. For instance, if all four
corners host Majorana Kramers pairs, both invariants are zero. Hence, when at
least one of the two invariants is nonzero, we immediately know the system
resides in the nontrivial higher-order phase, not the other way around. This
argument also applies to the BDI class which we shall deal with in the
following.

##  IV BDI class

###  IV.1 $\mathbb{Z}$ invariants

The Hamiltonian in the BDI class also satisfies the three intrinsic symmetries
listed in Eq.(3), except that the time-reversal operator has the property
$T^{2}=1$. This key difference determines that Fermi level crossings on high-
symmetry lines in the parameter space do not necessarily come in pairs, due to
the absence of Kramers degeneracy, as illustrated schematically in Fig. 8(b).
However, time-reversal symmetry still ensures that a crossing off the high-
symmetry lines, say at $(k_{F},\lambda_{F})$, has a partner at
$(-k_{F},\lambda_{F})$, only the two may not have opposite topological
charges, as the $D$ matrix doesn’t satisfy the relation
$D_{a}(k,\lambda)=-D^{T}_{a}(-k,\lambda)$. Hence each crossing within the pair
needs to be examined individually. Consequently, the higher-order topological
invariants in this case are of $\mathbb{Z}$ type.

In 2D, the BDI class does not have stable first-order topology but can exhibit
nontrivial higher-order topology [2]. Here, we continue to consider systems
with gapped bulk and edges, meaning that Hamiltonian $H_{a}(k,\lambda)$ is
gapped for both $\lambda=1$ and $\lambda=0$. We may associate a topological
invariant to the boundary-modulated Hamiltonian for generic $\lambda$, akin to
the $\tilde{\nu}_{a}(\lambda)$ in Eq.(6) for the DIII class. The topological
invariant in the BDI class is the winding number of $D$ matrix along the loop
with constant $k$ and takes the form

$\tilde{\omega}_{a}(\lambda)=\frac{1}{2\pi
i}\int_{-\pi}^{\pi}dk\partial_{k}\ln\det[D_{a}(k,\lambda)].$ (30)

Similar to the DIII class, the higher-order topological invariant $\omega_{a}$
here is defined to be the first-order topological invariant of the quasi-1D
cylindrical Hamiltonian $H_{a}(k,0)$, i.e.,

$\omega_{a}=\tilde{\omega}_{a}(0).$ (31)

Meanwhile, $\tilde{\omega}_{a}(1)$ indicates the first-order topology of a
toroidal system and is hence always zero.

Following from Eqs.(30) and (31), $\omega_{a}$ may take any integer value. In
a 1D system, this invariant reflects the difference in the count of Majorana
zero modes with distinct chiralities ($\pm 1$), namely the eigenvalues of
chiral operator, at one end. For the Hamiltonian $H_{a}(k,\lambda)$ in
consideration, we may denote $N_{a\pm}^{L(R)}$ as the number of Majorana zero
modes at left (right) end with chirality $\pm 1$. Then we have

$\omega_{a}=N_{a+}^{L}-N_{a-}^{L}=N_{a-}^{R}-N_{a+}^{R}.$ (32)

In a truly 1D system, Majorana zero modes of different chiralities can couple,
resulting in only one type of Majorana modes residing at each end if no other
symmetries are enforced. However, in the quasi-1D system here, each end,
comprising two spatially separated corners, can accommodate stable Majorana
corner modes of different types, as depicted in Fig. 8(a). The system enters
nontrivial higher-order phases when one or both of the invariants
$(\omega_{x},\omega_{y})$ take nonzero values.

Figure 8: (a) The relation between higher-order topological invariants and
Majorana corner modes in the BDI class. The invariant $\omega_{y}$ corresponds
to the difference in the number of Majorana zero modes with opposite
chiralities at two corners on the left, while $\omega_{x}$ counts the two
bottom corners. Blue and red dots represent Majorana modes with chiralities
$+1$ and $-1$, respectively. (b) Fermi level crossings in the BDI class. The
crossings off high-symmetry lines always come in pairs. On the line, the
crossing may appear alone. In contrast to the DIII class, signs of topological
charges for time-reversal partners, as indicated by colors of the crosses, can
be identical. Figure 9: Fermi level crossings and higher-order topological
invariants in three distinct phases of the BDI model. The pair of higher-order
invariants, $(\omega_{x},\omega_{y})$, is determined by topological charges of
Fermi level crossings. In (a-b), (d)-(e) and (g)-(h), the crosses mark the
relative positions of Fermi level crossings, and the numbers adjacent to each
cross indicate their respective topological charges, which are obtained from
winding numbers of the $D$ matrix along loops, $l_{1}$ and $l_{2}$, enclosing
each crossing. (c), (f), (i) The distributions of Majorana corner modes. The
solid blue and red circles represent Majorana zero modes with chiralities $+1$
and $-1$, respectively. Insets showcase energy spectra for an $80\times 80$
open-boundary system. In these plots, $m_{1}=0.3t$ and $\delta m=0.1t$.

To establish the relationship between higher-order invariants and Fermi level
crossings, we consider a large loop, $\mathcal{C}_{3}$, along the boundary of
the parameter space. This loop consists of four line segments, $l_{1}$ to
$l_{4}$, as depicted in Fig. 8(b). The winding number of the $D$ matrix along
this giant loop is given by

$\displaystyle\Omega$ $\displaystyle=\frac{1}{2\pi
i}\sum_{i=1}^{4}\int_{l_{i}}dl\nabla_{l}\ln\det[D_{a}(k,\lambda)]$
$\displaystyle=\sum_{(k_{F},\lambda_{F})}n_{a}(k_{F},\lambda_{F}),$ (33)

where the summation in the last equality runs over all Fermi level crossings
in the parameter space. The line integrals along $l_{1}$ and $l_{3}$ exactly
cancel each other. We are thus left with integrals along $l_{2}$ and $l_{4}$,
which essentially represent the winding number $\tilde{\omega}(1)$ and
$\tilde{\omega}(0)$ as defined in Eq.(30), up to a sign difference. Hence we
have

$\Omega=\tilde{\omega}_{a}(0)-\tilde{\omega}_{a}(1).$ (34)

Given that $\tilde{\omega}_{a}(1)=0$, we now arrive at the relation between
higher-order topological invariants and Fermi level crossings, which reads

$\omega_{a}=\sum_{(k_{F},\lambda_{F})}n_{a}(k_{F},\lambda_{F}).$ (35)

It follows from Eq.(35) that the sign of the topological charge for each
crossing matters, and that all the Fermi level crossings together determine
the overall higher-order topological invariants.

###  IV.2 An example

As a demonstration, we consider a 2D Hamiltonian in the BDI class, which is
based on the model described in Eq.(27). The Hamiltonian takes the form

$H_{\text{BDI}}(\bm{k})=H_{\text{DIII}}(\bm{k})-m_{2}\sigma_{y}.$ (36)

Due to the additional $\sigma_{y}$ term, this Hamiltonian does not preserve
time-reversal symmetry of the DIII-class Hamiltonian, which is
$T=-is_{y}\mathcal{K}$. Instead, the model exhibits a different time-reversal
symmetry with $T=\sigma_{x}\mathcal{K}$. Since the latter time-reversal
operator satisfies $T^{2}=1$, the Hamiltonian falls into BDI classification.
The chiral symmetry operator is now represented as $S=\tau_{x}\sigma_{x}$.

In the case where $m_{1}=m_{2}=0$, the system could possibly host a pair of
Fermi level crossings on high-symmetry lines, depending on the relative
strength of $m$ and $\delta m$. We already know that $m_{1}$ term doesn’t
couple the two crossings but instead separates them and drives them away from
the high-symmetry lines. The $m_{2}$ term doesn’t couple them either, albeit
for a different reason. Due to the redefined chiral symmetry, the two
crossings within the pair now carry the same topological charge, implying they
cannot annihilate with each other. At the onset of $m_{2}$, the pair moves
along the high-symmetry line in opposite directions. In this situation,
turning on $m_{1}$ cannot immediately move each crossing off the line because
of the time-reversal symmetry. However, for sufficiently large $m_{1}$, the
two separated crossings may meet again. If we continue to increase $m_{1}$,
the pair would move off the line. Therefore, the positions of the two
crossings vary according to the relative strength of $m_{1}$ and $m_{2}$. When
the $m_{1}$ term is dominant, one might anticipate that the pair of Fermi
level crossings will occur off the high-symmetry lines at $(\pm
k_{F},\lambda_{F})$. On the contrary, if the $m_{2}$ term is large compared to
$m_{1}$, they may both reside on the high-symmetry lines, but at distinct
points.

To illustrate these arguments, we present three representative cases with
varying $m_{2}$ in Fig. 9. It’s evident that the number and positions of Fermi
level crossings, along with their topological charges, differ in each case,
leading to distinct topological invariants. For smaller $m_{2}$, two crossings
appear off the high-symmetry lines and have identical topological charges, as
shown in Figs. 9(a)-(c). Following Eq.(35), we determine the topological
invariants to be $(\omega_{x},\omega_{y})=(2,-2)$. As a result, four Majorana
zero modes are found distributed across two opposite corners. At each corner,
the two zero modes share the same chirality, as indicated by colors of the
solid circles in Fig. 9(c). These topological invariants are further
corroborated using Eq.(32), by noting that $N_{x-}^{L}=N_{y+}^{L}=0$ and
$N_{x+}^{L}=N_{y-}^{L}=2$.

As $m_{2}$ increases, the two crossings are driven towards the high-symmetry
line $K=0$ where they will meet with each other. Subsequently, the pair is
pushed in opposite directions along the line. One of the crossings gradually
approaches $\lambda=0$, and annihilate with its counterpart at negative
$\lambda$. As depicted in Figs. 9(d) and (e) for $m=0.4$, two crossings are
located on the line $K=0$ in the $(k_{x},\lambda_{y})$ space, whereas only one
crossing are left in the $(k_{y},\lambda_{x})$ space as the other one has
annihilated with the crossing with $-\lambda_{x}$ when they meet at
$\lambda_{x}=0$. In this case, the topological invariants are
$(\omega_{x},\omega_{y})=(1,-2)$. There are still four Majorana zero modes at
the corners, but their distribution varies from the previous case. As shown in
Fig. 9(f), two Majorana zero modes with chirality $-1$ group in one corner,
while the others with chirality $+1$ occupy two separate corners, i.e.,
$N_{x+}^{L}=1$, $N_{x-}^{L}=N_{y+}^{L}=0$ and $N_{y-}^{L}=2$.

With a further increase in $m_{2}$, the two crossings in the
$(k_{x},\lambda_{y})$ space remain. In the $(k_{y},\lambda_{x})$ space,
however, another crossing emerges on the high-symmetry line carrying an
opposite charge. We then have $(\omega_{x},\omega_{y})=(0,-2)$. Four Majorana
modes are still present, but their locations differ from the previous cases.
Now, each corner hosts one Majorana mode, with $N_{x+}^{L}=1$, $N_{x-}^{L}=1$,
$N_{y+}^{L}=0$ and $N_{y-}^{L}=2$.

In all three cases, at least one of the two invariants is nonzero, indicating
the nontrivial higher-order topology. It is possible that two systems with
identical topological invariants can have distinct distributions, or even
different numbers of Majorana corner modes. Nonetheless, the presence of
nontrivial higher-order topology in BDI class is assured whenever at least one
of the invariants assumes a nonzero value.

##  V Discussions and conclusions

We have established that Fermi level crossings can act as indicators for
higher-order topology in DIII and BDI classes. With variations in a bulk
parameter, topological phase transitions may occur when either the number or
the topological charges of these crossings change, even if the bulk gap
remains open. The close relationship between Fermi level crossings and the
bulk Hamiltonian directly illustrates the bulk-boundary correspondence in
higher-order topological superconductors, irrespective of crystalline
symmetries. The fact that Fermi level crossings could change when either the
bulk or edge gap closes, highlights a striking contrast to conventional bulk-
boundary correspondence, where topological phase transitions typically occur
only when the bulk gap closes.

In the DIII class, the $\mathbb{Z}_{2}$ charge for crossings on high-symmetry
lines is obtained from the root’s multiplicity by solving Eq.(24) in the case
of linear band crossings. For other types of crossings, this argument still
holds. Let’s first consider a pair of band crossings at $(K,\lambda_{F})$,
which means only $m_{a,1}$ in Eq.(13) equals zero at the crossing, and assume
$m_{a,1}\propto(\lambda-\lambda_{F})^{r}$ near the crossing. The exponent $r$
indicates the type of crossings. Apparently, only if $r$ is odd will $m_{a,1}$
change its sign after crossing $\lambda_{F}$. Hence, $\eta_{a}=1$ for odd $r$
and $\eta_{a}=0$ for even $r$. For instance, the quadratic band crossing
corresponds to $r=2$, in which case $\eta_{a}$=0, suggesting that the Pfaffian
of the $D$ matrix doesn’t experience a $\pi$-phase jump across $\lambda_{F}$.
Given that $\det[D_{a}(K,\lambda)]\propto(\lambda-\lambda_{F})^{2r}$ near the
crossing, we immediately obtain the root’s multiplicity as $p_{a}=2r$. The
parity of $p_{a}/2$ is exactly equal to that of $\eta_{a}$. For a generic
case, $\eta_{a}$ is given by the number of $m_{a,j}$ that changes sign after
crossing $\lambda_{F}$, which equals the number of crossing pairs with odd
$r$. The root’s multiplicity is expressed as $p_{a}=2\rho_{r}r$, where
$\rho_{r}$ represents the number of crossing pairs with exponent $r$. One can
verify that $p_{a}/2\text{ mod }2=\eta_{a}\text{ mod }2$. Therefore, the
$\mathbb{Z}_{2}$ charge of Fermi level crossings on the high-symmetry lines
can be determined by the root’s multiplicity, independent of the crossing
type. For a crossing away from the high-symmetry lines, it is necessary to
evaluate its topological charge defined in Eq.(26), which carries opposite
charge to its time-reversal partner. In this case, the $\mathbb{Z}_{2}$ charge
defined for the time-reversal pair is nonzero when the topological charge of
either crossing within the pair is odd. An interesting question arises about
whether one may determine the $\mathbb{Z}_{2}$ charge based solely on the
root’s multiplicity, akin to the crossings on high-symmetry lines. For
example, consider two crossings off high-symmetry lines at $(\pm
k_{F},\lambda_{F})$, each with multiplicity $p_{a}$. If we move these
crossings towards the high-symmetry lines, which doesn’t affect the
$\mathbb{Z}_{2}$ charge, the multiplicity at the meeting point
$(K,\lambda_{F})$ on the high-symmetry lines becomes $2p_{a}$ unless the
crossing type changes in this process. This implies that under these
conditions, the $\mathbb{Z}_{2}$ charge could be inferred from the root’s
multiplicity. However, it’s challenging to conclusively prove that the
crossing type remains the same, which means we don’t know whether $p_{a}$
remains constant before the crossing pair arrive at the high-symmetry lines,
though it seems unlikely to change in most cases. In the model we studied in
Sec. IIID, $p_{a}$ for crossings off the high-symmetry lines doesn’t change
while the crossings move. In BDI classes, the distinction between Fermi level
crossings on and off high-symmetry lines becomes irrelevant, and the root’s
multiplicity seems also to lose its significance. Here, the higher-order
invariants rely on the cumulative topological charges of all crossings.

In this paper, we’ve characterized higher-order topology through pairs of
topological invariants: $(\nu_{x},\nu_{y})$ for the DIII class, as defined in
Eq.(15), and $(\omega_{x},\omega_{y})$ for the BDI class, introduced in
Eq.(35). A system enters nontrivial higher-order phases when at least one of
these invariants becomes nonzero. However, the reverse isn’t necessarily true.
Nontrivial higher-order topology can occur even if both topological invariants
are zero. For example, in the DIII class, if Majorana Kramers pairs are
present at all four corners, the invariants may still be zero,
$(\nu_{x},\nu_{y})=0$. This happens because these invariants represent the
first-order topology of a quasi-1D cylindrical system, where each end
corresponds to two corners in a square sample. When each corner supports a
Majorana pair in the DIII class, their combined effect results in trivial
invariants, due to the $\mathbb{Z}_{2}$ classification. A similar scenario
occurs in the BDI class, where trivial invariants can arise if the two corners
at each end host Majorana zero modes with opposing chiralities. This can also
be viewed from the perspective of edge topology, where these invariants
actually indicate the topological differences between opposite edges, as we
have established in our previous study of D-class systems [46]. However, in
higher-order topological phases, topological differences may only exist
between adjacent edges, rendering these invariants inadequate. For a
comprehensive bulk-boundary correspondence in such cases, other general
twisted boundary conditions, involving phase modulations of the boundary
Hamiltonian, might be necessary [52, 53, 54, 37, 55, 56]. Our research
underscores that boundary condition modulation can reveal a direct connection
between bulk and corners. Future explorations within this framework may
encompass a wider array of higher-order phases.

In summary, our study delves into 2D higher-order topological superconductors
in the DIII and BDI classes, both characterized by the presence of chiral
symmetry. We demonstrate that Fermi level crossings, occurring when boundary
conditions vary, can act as reliable indicators of nontrivial higher-order
topology. These crossings also provide a direct and intuitive understanding of
how bulk properties influence the emergence of Majorana zero modes at the
corners. A pair of higher-order topological invariants are introduced, which
are associated with topological charges of Fermi level crossings. The
invariants indicate nontrivial higher-order phases when at least one of them
assumes a nonzero value.

###### Acknowledgements.

This work was supported by National Science Foundation of China (NSFC) under
Grant No. 11704305, and the Innovation Program for Quantum Science and
Technology (2021ZD0302400).

##  Appendix A The matrix $\tilde{A}_{a}$

In this appendix, we will present the details in deriving the matrix
$\tilde{A}_{a}(k,\lambda)$ introduced in Eq.(23), whose determinant is shown
to be equal to that of $A_{a}(k,\lambda)$ in Eq.(21).

Firstly, we note that the $D$ matrix of bulk Hamiltonian in momentum space is
given by $D(\bm{k})=u^{\dagger}H(\bm{k})v$, where $U_{S}=(u\ v)$ is the
unitary matrix that diagonalizes the chiral symmetry operator $S$. For the
boundary-modulated Hamiltonian $H_{a}(k,\lambda)$ represented in the basis
$\Psi_{a}(k)=\\{\psi_{a,1}(k),...,\psi_{a,j}(k),...,\psi_{a,N_{a}}(k)\\}$,
with $\psi_{a,j}(k)$ including all internal degrees of freedom, denoted as
$[c_{a}(k)]_{j\alpha}$, on the site $j$, the corresponding matrices $u$ and
$v$ are replaced by $u_{a}=\mathbb{I}_{N_{a}\times N_{a}}\otimes u$ and
$v_{a}=\mathbb{I}_{N_{a}\times N_{a}}\otimes v$, respectively. Here, $N_{a}$
represents the number of unit cells along the $a$ direction. As such, the $D$
matrix of boundary-modulated Hamiltonian is given by
$D_{a}(k,\lambda)=u_{a}^{\dagger}H_{a}(k,\lambda)v_{a}$, and similarly for the
boundary Hamiltonian $D_{a}^{b}(k)=u_{a}^{\dagger}B_{a}(k)v_{a}$. The entries
of $A_{a}$ can be written as

$\displaystyle[A_{a}(k,\lambda)]_{i\bar{\alpha},j\bar{\beta}}$
$\displaystyle=\delta_{ij}\delta_{\bar{\alpha}\bar{\beta}}-(1-\lambda)$ (37)
$\displaystyle\times[D_{a}^{-1}(k,1)]_{i\bar{\alpha},l\bar{\gamma}}[D^{b}_{a}(k)]_{l\bar{\gamma},j\bar{\beta}},$

with Einstein summation assumed. In Eq.(37) we distinguish $\bar{\alpha}$ and
$\bar{\beta}$ from $\alpha$ and $\beta$ as the former indices run over only
half the inner degrees of freedom. The two representations of bulk
Hamiltonian, $D_{a}(k_{\bar{a}},1)$ and $D(\bm{k})$, are related with each
other through Fourier transformation, represented in matrix form as

$[D_{a}(k_{\bar{a}},1)]_{i\bar{\alpha},j\bar{\beta}}=\sum_{k_{a}}[U_{F}^{\dagger}]_{i,k_{a}}[D(\bm{k})]_{\bar{\alpha},\bar{\beta}}[U_{F}]_{k_{a},j}.$
(38)

Here, the subscript of $k$ is explicitly shown for clarity, and the Fourier
matrix $[U_{F}]_{k_{a},j}=\frac{1}{\sqrt{N_{a}}}e^{-ijk_{a}}$. Since the
boundary Hamiltonian operates only on boundary cells, it follows that the
entries
$[A_{a}(k,\lambda)]_{i\bar{\alpha},j\bar{\beta}}=\delta_{ij}\delta_{\bar{\alpha}\bar{\beta}}$
if $j$ does not represent a boundary cell. Consequently, the determinant of
$A_{a}$ can be calculated using only the block corresponding to boundary
cells, denoted by $\bar{A}_{a}(k,\lambda)$, with

$\det[A_{a}(k,\lambda)]=\det[\bar{A}_{a}(k,\lambda)].$ (39)

Combing Eqs.(37) and (38), we arrive at

$\displaystyle[\bar{A}_{a}$
$\displaystyle(k_{\bar{a}},\lambda_{a})]_{i_{b}\bar{\alpha},j_{b}\bar{\beta}}=\delta_{i_{b},j_{b}}\delta_{\bar{\alpha},\bar{\beta}}-\frac{1}{2\pi}(1-\lambda_{a})$
(40) $\displaystyle\times\int
dk_{a}[F(k_{a})]_{i_{b},l_{b}}[D^{-1}(\bm{k})]_{\bar{\alpha},\bar{\gamma}}[\bar{D}^{b}_{a}(k_{\bar{a}})]_{l_{b}\bar{\gamma},j_{b}\bar{\beta}},$

where $\bar{D}^{b}_{a}$ represents the block of $D^{b}_{a}$ corresponding to
boundary cells labelled by $i_{b}$, $j_{b}$,
$[F(k_{a})]_{i_{b},l_{b}}=e^{ik_{a}(i_{b}-l_{b})}$, and the summation over
$k_{a}$ is replaced by integration.

In certain cases, $\bar{D}^{b}_{a}(k)$ may not be a full-rank matrix. Its
singular value decomposition is represented as

$[\bar{D}^{b}_{a}(k)]_{i_{b}\bar{\alpha},j_{b}\bar{\beta}}=[\bar{U}_{a}^{b}(k)]_{i_{b}\bar{\alpha},n}[\bar{\Lambda}_{a}^{b}(k)]_{n,n}[\bar{V}_{a}^{b\dagger}(k)]_{n,j_{b}\bar{\beta}},$
(41)

where $\bar{\Lambda}_{a}^{b}$ is a diagonal matrix with $r_{a}$ nonzero
entries, with $r_{a}$ representing the rank of $\bar{D}^{b}_{a}(k)$. In
addition, $\bar{\Lambda}_{a}^{b}$ can be arranged such that
$[\bar{\Lambda}_{a}^{b}]_{n,n}\neq 0$ for $n\leq r_{a}$. Given that
$\bar{U}_{a}^{b}$ and $\bar{V}_{a}^{b}$ are unitary matrices, it follows that

$\det[\bar{A}_{a}(k,\lambda)]=\det[\bar{V}_{a}^{b\dagger}(k)\bar{A}_{a}(k,\lambda)\bar{V}_{a}^{b}(k)].$
(42)

From Eqs.(40) and (41) we know
$[\bar{V}_{a}^{b\dagger}\bar{A}_{a}\bar{V}_{a}^{b}]_{m,n}=\delta_{m,n}$ if
$n>r_{a}$. Hence we only need to consider the entries where $m,n\leq r_{a}$
while calculating its determinant. For convenience, we group the first $r_{a}$
columns of $\bar{V}_{a}^{b}$ into a new matrix, denoted as
$\tilde{V}_{a}^{b}$. Define

$\displaystyle\tilde{A}_{a}(k,\lambda)=\tilde{V}_{a}^{b\dagger}(k)\bar{A}_{a}(k,\lambda)\tilde{V}_{a}^{b}(k),$
(43)

and thus

$\det[\bar{A}_{a}(k,\lambda)]=\det[\tilde{A}_{a}(k,\lambda)].$ (44)

Combining Eqs.(40) and (43) we arrive at Eq.(23) in the main text, which shows
the specific form of this $r_{a}\times r_{a}$ matrix
$\tilde{A}_{a}(k,\lambda)$.

## References

  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Topological insulators and superconductors, Rev. Mod. Phys. 83, 1057 (2011).
  * Chiu  _et al._ [2016] C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Classification of topological quantum matter with symmetries, Rev. Mod. Phys. 88, 035005 (2016).
  * Read and Green [2000] N. Read and D. Green, Paired states of fermions in two dimensions with breaking of parity and time-reversal symmetries and the fractional quantum Hall effect, Phys. Rev. B 61, 10267 (2000).
  * Kitaev [2001] A. Y. Kitaev, Unpaired Majorana fermions in quantum wires, Phys.-Usp. 44, 131 (2001).
  * Wilczek [2009] F. Wilczek, Majorana returns, Nature Phys 5, 614 (2009).
  * Alicea [2012] J. Alicea, New directions in the pursuit of Majorana fermions in solid state systems, Rep. Prog. Phys. 75, 076501 (2012).
  * Stanescu and Tewari [2013] T. D. Stanescu and S. Tewari, Majorana fermions in semiconductor nanowires: Fundamentals, modeling, and experiment, J. Phys.: Condens. Matter 25, 233201 (2013).
  * Beenakker [2013] C. Beenakker, Search for Majorana Fermions in Superconductors, Annu. Rev. Condens. Matter Phys. 4, 113 (2013).
  * Elliott and Franz [2015] S. R. Elliott and M. Franz,  _Colloquium_ : Majorana fermions in nuclear, particle, and solid-state physics, Rev. Mod. Phys. 87, 137 (2015).
  * Aguado [2017] R. Aguado, Majorana quasiparticles in condensed matter, Riv. Nuovo Cimento 40, 523 (2017).
  * Benalcazar  _et al._ [2017a] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Quantized electric multipole insulators, Science 357, 61 (2017a).
  * Benalcazar  _et al._ [2017b] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Electric multipole moments, topological multipole moment pumping, and chiral hinge states in crystalline insulators, Phys. Rev. B 96, 245115 (2017b).
  * Langbehn  _et al._ [2017] J. Langbehn, Y. Peng, L. Trifunovic, F. von Oppen, and P. W. Brouwer, Reflection-Symmetric Second-Order Topological Insulators and Superconductors, Phys. Rev. Lett. 119, 246401 (2017).
  * Song  _et al._ [2017] Z. Song, Z. Fang, and C. Fang, ( d - 2 ) -Dimensional Edge States of Rotation Symmetry Protected Topological States, Phys. Rev. Lett. 119, 246402 (2017).
  * Schindler  _et al._ [2018a] F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. P. Parkin, B. A. Bernevig, and T. Neupert, Higher-order topological insulators, Sci. Adv. 4, eaat0346 (2018a).
  * Schindler  _et al._ [2018b] F. Schindler, Z. Wang, M. G. Vergniory, A. M. Cook, A. Murani, S. Sengupta, A. Y. Kasumov, R. Deblock, S. Jeon, I. Drozdov, H. Bouchiat, S. Guéron, A. Yazdani, B. A. Bernevig, and T. Neupert, Higher-order topology in bismuth, Nature Phys 14, 918 (2018b).
  * Benalcazar and Cerjan [2022] W. A. Benalcazar and A. Cerjan, Chiral-Symmetric Higher-Order Topological Phases of Matter, Phys. Rev. Lett. 128, 127601 (2022).
  * Khalaf [2018] E. Khalaf, Higher-order topological insulators and superconductors protected by inversion symmetry, Phys. Rev. B 97, 205136 (2018).
  * Zhu [2018] X. Zhu, Tunable Majorana corner states in a two-dimensional second-order topological superconductor induced by magnetic fields, Phys. Rev. B 97, 205134 (2018).
  * Yan  _et al._ [2018] Z. Yan, F. Song, and Z. Wang, Majorana Corner Modes in a High-Temperature Platform, Phys. Rev. Lett. 121, 096803 (2018).
  * Wang  _et al._ [2018a] Q. Wang, C.-C. Liu, Y.-M. Lu, and F. Zhang, High-Temperature Majorana Corner States, Phys. Rev. Lett. 121, 186801 (2018a).
  * Wang  _et al._ [2018b] Y. Wang, M. Lin, and T. L. Hughes, Weak-pairing higher order topological superconductors, Phys. Rev. B 98, 165144 (2018b).
  * Yan [2019] Z. Yan, Higher-Order Topological Odd-Parity Superconductors, Phys. Rev. Lett. 123, 177001 (2019).
  * Zhang  _et al._ [2019a] R.-X. Zhang, W. S. Cole, and S. Das Sarma, Helical Hinge Majorana Modes in Iron-Based Superconductors, Phys. Rev. Lett. 122, 187001 (2019a).
  * Ghosh  _et al._ [2021] A. K. Ghosh, T. Nag, and A. Saha, Hierarchy of higher-order topological superconductors in three dimensions, Phys. Rev. B 104, 134508 (2021).
  * Scammell  _et al._ [2022] H. D. Scammell, J. Ingham, M. Geier, and T. Li, Intrinsic first- and higher-order topological superconductivity in a doped topological insulator, Phys. Rev. B 105, 195149 (2022).
  * Geier  _et al._ [2018] M. Geier, L. Trifunovic, M. Hoskam, and P. W. Brouwer, Second-order topological insulators and superconductors with an order-two crystalline symmetry, Phys. Rev. B 97, 205135 (2018).
  * Trifunovic and Brouwer [2019] L. Trifunovic and P. W. Brouwer, Higher-Order Bulk-Boundary Correspondence for Topological Crystalline Phases, Phys. Rev. X 9, 011012 (2019).
  * Skurativska  _et al._ [2020] A. Skurativska, T. Neupert, and M. H. Fischer, Atomic limit and inversion-symmetry indicators for topological superconductors, Phys. Rev. Research 2, 013064 (2020).
  * Ono  _et al._ [2020] S. Ono, H. C. Po, and H. Watanabe, Refined symmetry indicators for topological superconductors in all space groups, Sci. Adv. 6, eaaz8367 (2020).
  * Takahashi  _et al._ [2020] R. Takahashi, Y. Tanaka, and S. Murakami, Bulk-edge and bulk-hinge correspondence in inversion-symmetric insulators, Phys. Rev. Research 2, 013300 (2020).
  * Hsu  _et al._ [2020] Y.-T. Hsu, W. S. Cole, R.-X. Zhang, and J. D. Sau, Inversion-protected Higher-order Topological Superconductivity in Monolayer WTe 2, Phys. Rev. Lett. 125, 097001 (2020).
  * Kruthoff  _et al._ [2017] J. Kruthoff, J. de Boer, J. van Wezel, C. L. Kane, and R.-J. Slager, Topological Classification of Crystalline Insulators through Band Structure Combinatorics, Phys. Rev. X 7, 041069 (2017).
  * Tang  _et al._ [2019] F. Tang, H. C. Po, A. Vishwanath, and X. Wan, Comprehensive search for topological materials using symmetry indicators, Nature 566, 486 (2019).
  * Zhang  _et al._ [2019b] T. Zhang, Y. Jiang, Z. Song, H. Huang, Y. He, Z. Fang, H. Weng, and C. Fang, Catalogue of topological electronic materials, Nature 566, 475 (2019b).
  * Vergniory  _et al._ [2019] M. G. Vergniory, L. Elcoro, C. Felser, N. Regnault, B. A. Bernevig, and Z. Wang, A complete catalogue of high-quality topological materials, Nature 566, 480 (2019).
  * Zhang [2022] R.-X. Zhang, Bulk-Vortex Correspondence of Higher-Order Topological Superconductors (2022), arxiv:2208.01652 .
  * Roberts  _et al._ [2020] E. Roberts, J. Behrends, and B. Béri, Second-order bulk-boundary correspondence in rotationally symmetric topological superconductors from stacked Dirac Hamiltonians, Phys. Rev. B 101, 155133 (2020).
  * Huang and Hsu [2021] S.-J. Huang and Y.-T. Hsu, Faithful derivation of symmetry indicators: A case study for topological superconductors with time-reversal and inversion symmetries, Phys. Rev. Research 3, 013243 (2021).
  * Tang  _et al._ [2022] F. Tang, S. Ono, X. Wan, and H. Watanabe, High-Throughput Investigations of Topological and Nodal Superconductors, Phys. Rev. Lett. 129, 027001 (2022).
  * Zhu [2019] X. Zhu, Second-Order Topological Superconductors with Mixed Pairing, Phys. Rev. Lett. 122, 236401 (2019).
  * Khalaf  _et al._ [2021] E. Khalaf, W. A. Benalcazar, T. L. Hughes, and R. Queiroz, Boundary-obstructed topological phases, Phys. Rev. Research 3, 013239 (2021).
  * Tiwari  _et al._ [2020] A. Tiwari, A. Jahin, and Y. Wang, Chiral Dirac superconductors: Second-order and boundary-obstructed topology, Phys. Rev. Research 2, 043300 (2020).
  * Volpez  _et al._ [2019] Y. Volpez, D. Loss, and J. Klinovaja, Second-Order Topological Superconductivity in $\pi$ -Junction Rashba Layers, Phys. Rev. Lett. 122, 126402 (2019).
  * Wu  _et al._ [2020] X. Wu, W. A. Benalcazar, Y. Li, R. Thomale, C.-X. Liu, and J. Hu, Boundary-Obstructed Topological High- T c Superconductivity in Iron Pnictides, Phys. Rev. X 10, 041014 (2020).
  * Wang and Zhu [2023] H. Wang and X. Zhu, Higher-order topological superconductors characterized by Fermi level crossings, Phys. Rev. B 108, 125426 (2023).
  * Tanaka  _et al._ [2020] Y. Tanaka, R. Takahashi, and S. Murakami, Appearance of hinge states in second-order topological insulators via the cutting procedure, Phys. Rev. B 101, 115120 (2020).
  * Schnyder and Ryu [2011] A. P. Schnyder and S. Ryu, Topological phases and surface flat bands in superconductors without inversion symmetry, Phys. Rev. B 84, 060504(R) (2011).
  * Qi  _et al._ [2010] X.-L. Qi, T. L. Hughes, and S.-C. Zhang, Topological invariants for the Fermi surface of a time-reversal-invariant superconductor, Phys. Rev. B 81, 134508 (2010).
  * Rhim  _et al._ [2018] J.-W. Rhim, J. H. Bardarson, and R.-J. Slager, Unified bulk-boundary correspondence for band insulators, Phys. Rev. B 97, 115143 (2018).
  * Wimmer [2012] M. Wimmer, Algorithm 923: Efficient Numerical Computation of the Pfaffian for Dense and Banded Skew-Symmetric Matrices, ACM Trans. Math. Softw. 38, 1 (2012).
  * Niu  _et al._ [1985] Q. Niu, D. J. Thouless, and Y.-S. Wu, Quantized Hall conductance as a topological invariant, Phys. Rev. B 31, 3372 (1985).
  * Qi  _et al._ [2006] X.-L. Qi, Y.-S. Wu, and S.-C. Zhang, General theorem relating the bulk topological number to edge states in two-dimensional insulators, Phys. Rev. B 74, 045125 (2006).
  * Song  _et al._ [2020] Z.-D. Song, L. Elcoro, and B. A. Bernevig, Twisted bulk-boundary correspondence of fragile topology, Science 367, 794 (2020).
  * Jahin  _et al._ [2024] A. Jahin, Y.-M. Lu, and Y. Wang, Many-body higher-order topological invariant for $C_n$-symmetric insulators (2024), arxiv:2401.00050 .
  * Lin and Lee [2024] L. Lin and C. Lee, Probing Chiral-Symmetric Higher-Order Topological Insulators with Multipole Winding Number (2024), arxiv:2401.03699 .


m O O

O O

O O O

O O

O O

O O

O O

A Mixed Integer Quadratic Program for Valuing the Impact of Price and
Forecast Uncertainty for Wind Generators

The information, data, or work
presented herein was supported by the National Science Foundation Graduate
Research Fellowship under Grant No. 2141064 as well as in part by the Advanced
Research Projects Agency-Energy (ARPA-E), U.S. Department of Energy under Award
Number DE-AR0001277. The views and opinions of authors expressed herein do not
necessarily state or reflect those of the United States Government or any
agency thereof.

Daniel Shen and Marija Ilic
Department of Electrical Engineering and Computer Science

Massachusetts Institute of Technology

Cambridge, USA

Owners of wind power plants are exposed to financial risk in wholesale
electricity markets due to the uncertain nature of wind forecasts and price
volatility. In the event of a wind shortfall, the plant may have to
repurchase power at a higher price in the real-time market. However, reducing
the power offered in the day-ahead market may also be interpreted by
regulators as physical withholding. We formulate and solve a mixed-integer
quadratic program (MIQP) that prices the uncertain portion of a wind
generator’s forecast to hedge against uncertainties and which addresses
concerns around withholding. We exploit the structure of the MIQP inputs to
introduce additional constraints to improve computation time. Additionally,
we provide a qualitative approach for generators and regulators to interpret
the results of the MIQP. Finally, we simulate a real-world application for a
wind farm in New York using past wind forecasts and NYISO prices.

renewable energy, electricity markets, uncertainty management

[A, 01]Set of price and wind forecasts.
[A, 02]Index over forecast scenarios.
[A, 03]Time indexer.
[A, 04]Offer curve segment indexer.
[A, 05]Total number of offer curve segments.

[B, 01]$\varsegmentpower$
Power offered in segment .

[B, 02]$\varsegmentprice$
Price of the power offered in segment .

[C, 01]
Indicator if segment $\indexsegment$ clears the market.

[C, 01]
Power cleared in day-ahead under scenario .

[C, 01]
Shortfall between the day-ahead clearing amount and the real-time forecast
realization; negative if there is a repurchase.

[D, 01]DA,RT
Day-ahead and real-time prices in scenario .

[D, 02]
Max. dispatchable wind power in scenario .

[E, 01]Risk factor for conditional value at risk.
[E, 02]$M$
Arbitrarily large constant for indicator constraints.

[E, 03]$\epsilon$
Arbitrarily small constant for indicator constraints.

§ INTRODUCTION

Significant amounts of variable renewable energy (VRE) must be installed on the
grid to achieve decarbonization targets. Such resources incur no fuel costs to
produce electricity, but their generation output is uncertain and the asset
owners face financial risk from underdelivery of power in the real-time market.
It is typically assumed in restructured electricity markets that generators
should offer power at a price that reflects their variable costs. While the
fuel portion of variable cost for VRE resources is zero, VRE generators do
incur nonzero costs in potential shortfall penalties as the volume of their
power sold on the day-ahead market increases. In the scenario where the
real-time electricity price is higher than the day-ahead electricity price, a
wind generator which always offers and clears at its full rated capacity will
incur costs from consistently repurchasing its real-time shortfall.

For work done previously on price-taker behavior offer/bid curve creation for
producers and consumers, we refer to <cit.> and
<cit.> <cit.>, respectively. Our work
builds upon the approach described in <cit.>, but for
VRE deviations instead of dispatch deviations by conventional generators.
This work also extends previous work on wind generator offer creation in
<cit.> to a piecewise linear context. We refer to
<cit.> for a detailed treatment of strategic
wind power offering with conditional value at risk.

This work extends some of the aforementioned work by combining the concept of
breaking the problem down into discrete block price-quantity segments, as is
required in nearly all restructured markets, and the concept of optimizing
those segments under the conditional value at risk (CVaR) risk measure. We
consider the problem of a price-taking wind generator pricing its variable
costs associated with uncertainty in a day-ahead offer curve. We propose a
mixed-integer quadratic program (MIQP) that considers CVaR to reflect the
variable costs of shortfall penalties. The MIQP takes into account the
distribution of day-ahead prices, real-time prices, and wind forecasts to
create an offer curve of energy blocks at specific marginal prices.

We wish to note the following contributions of this paper:

* Formulating a MIQP program which optimizes both the price and quantity
dimensions of the generator offer curve to reflect the variable costs of

* Reducing the number of unconstrained integer decision variables in the
MIQP by exploiting structure in the optimization problem to improve
computational speed.

* Presenting an approach for qualitatively interpreting the offer curve
results by vizualizing the distributions of the “active” price and wind
scenarios. Such a method could prove useful for generation operators to
justify their nonzero marginal cost offers to other stakeholders or market

§ METHODOLOGY

§.§ Problem Setup & Background

In restructured electricity markets, generators submit offer curves of blocks
of power offered at discrete marginal prices. We make the following
considerations for a wind generator:

* The agent takes the day-ahead and real-time price scenarios as
predicted (price-taker behavior) and also has a set of discrete forecast
scenarios of its real-time wind production.

* All power quantities offered at a price below the day-ahead price
clear the market and are sold at that day-ahead price.

* In the day-ahead optimization, the agent should not
consider any profits from selling or repurchasing electricity in the real-time
market to counteract claims that the generator is engaging in physical
withholding <cit.> and could potentially be exercising
market power.

* If the realized wind power is less in real-time than was contracted
for dispatch in the day-ahead market, the generator must buyback the
shortfall at the real-time price.

* The generator is downwards dispatchable to 0 MW from any realized wind
scenario and such downwards dispatch incurs no operating cost.

Electricity prices are generally more volatile than that of other commodities,
and generators may have a preference to reduce the risk of incurring financial
losses from price fluctuations. We integrate a popular risk measure,
conditional value at risk (CVaR) <cit.>, into our objective to allow
for modeling of the agent's aversion to risk. In using CVaR, we only consider
the most undesired scenarios as determined by the choice of offer curve. Fig.
<ref> gives a graphical representation of CVaR. The CVaR
parameter $\parambeta$ represents the agent's aversion to risk; as $\parambeta$
increases up to 1, the agent's aversion to risk also increases. Using a
sample-average approximation approach with 100 scenarios, the CVaR is equal to
the $(1-\parambeta) * 100$ lowest- profit scenarios. For $\parambeta=0$, the
problem is equivalent to optimizing the expectation.

Graphical representation of value at risk (VaR) and conditional
value at risk (CVaR).

§.§ Mixed-Integer Quadratic Program

Market operators require the offer curve of the generator to be a piecewise
linear function of total production costs ($/h) vs. output (MW) for the market
clearing software. This is equivalent to a step function that specifies the
marginal cost ($/MWh) as a function of output for $\nsegments$ steps. For the
remainder of this paper, the term “offer” will refer to the latter step

For a single period $\indextime$, the optimal offer based on considerations
in Section <ref> is found by solving the following MIQP:

\begin{equation}
\max_{\bm{x}} \quad CVaR_{\parambeta}
\left(
\rvprice{DA}[] \: \auxcleared[] + \rvprice{RT}[] \auxshortfall[]
\right) \label{eq:objective}
\end{equation}

\begin{align}
\text{subject to}\quad
& \auxcleared = \sum_{\indexsegment}^{\nsegments}
\varsegmentactive \varsegmentpower \label{eq:cleared} \\
& max(\auxcleared) \leq max(\rvwind) \label{eq:pmax} \\
& \auxshortfall = min(0, \rvwind - \auxcleared) \label{eq:shortfall} \\
& \varsegmentactive \leq 1 - \frac{\varsegmentprice - \rvprice{DA}}{M}
\label{eq:active_binary_zero} \\
& \varsegmentactive \geq \epsilon +
\frac{\rvprice{DA} - \varsegmentprice}{M}
\label{eq:active_binary_one}\\
& \varsegmentactive \in \{ 0, 1 \} \label{eq:active_binary_binary} \\
& \sum_i^{\nsegments} \varsegmentactive \leq
\sum_i^{\nsegments} \varsegmentactive[\indextime][(\indexscenario + 1)]
\label{eq:indicator-rowsum} \\
\varsegmentprice[\indexsegment] \leq \varsegmentprice[(\indexsegment + 1)] \label{eq:increasingprice} \\
& \varsegmentactive[\indextime][\indexscenario][(\indexsegment + 1)] \leq
\varsegmentactive \label{eq:indicator-decreasing}
\end{align}

\begin{equation*}
\text{where} \quad \bm{x} =
\varsegmentactive, \varsegmentprice, \varsegmentpower
\end{equation*}

Variables $\varsegmentprice$ and $\varsegmentpower$ specify the price and
quantity characteristics of each block in the offer curve
(Fig. <ref>). $\varsegmentactive$ is 1 if a block
$\indexsegment$ clears in a specific scenario and 0 otherwise.

The objective is to maximize the revenue cleared in the day-ahead market
subject to the $CVaR_\parambeta$ risk measure. Each scenario $\indexscenario$
consists of a hourly day-ahead price sample $\rvprice{DA}$, real-time price
sample $\rvprice{RT}$, and a forecast for the maximum wind power $\rvwind$.

The constraints can be interpreted as follows:

* (<ref>) gives the quantity of power cleared in
the day-ahead market under scenario $\indexscenario$.

* (<ref>) restricts the generator from submitting a total
quantity offer that exceeds the maximum forecast scenario.

* (<ref>) restricts the optimization from considering
net profit obtained from withholding power for the real-time market; the
variable $\auxshortfall$ is negative if the generator is buying back power to
cover a shortfall.

* (<ref>), (<ref>), and
(<ref>) reflect which of the offered blocks clear
for a given scenario in the day-ahead market; is is assume that if the market
price is equal to the offer price, the entire quantity associated with the
block clears. For a graphical representation of this constraint, see
Fig. <ref>.

* (<ref>) follows from assuming the scenarios are
presorted by increasing day-ahead prices.

* (<ref>) follows from
(<ref>), since lower-priced blocks will always clear
if a higher-priced block clears as well.

Note that constraints (<ref>)-(<ref>)
are not necessary to solve the optimization problem, but are added to reduce
the feasible solution search space.

Practically, when solving (<ref>) on a commercial solver, we also
introduce an $\ell_2$ regularization term for $\varsegmentprice$ and
$\varsegmentpower$ to the objective to improve numerical stability.

Representation of offer curve parameters. Observe that based on the
indicated forecast price in scenario , only segments 1 and 2
would clear and
$[\varsegmentactive[1], \varsegmentactive[2], \varsegmentactive[3]] = [1, 1, 0]$

§ RESULTS

§.§ Curve Interpretation for Multivariate Gaussian Scenarios

In this section, we show the results of running the MIQP described in Section
<ref> on two different synthetic Gaussian distributions of scenarios
to examine the behavior of the program and demonstrate a qualitative method
of interpreting the curve shapes. Table
<ref> describes the simulation parameters
for the multivariate Gaussian used to generate the two cases. We denote the
case where wind and real-time prices are negatively correlated as Case 1,
and the case where they are positively correlated as Case 2. For both curves,
we set the program to generate six offer curve segments and then remove
segments that are smaller than 1 MW in width to improve readability.

Fig. <ref> shows the offer curves created by running the
MIQP for Cases 1 and 2. We observe that in both cases, the right side of the
offer curve reflects the nonzero variable costs associated with buyback
uncertainty, as would be expected.

For Case 1, increasing $\parambeta$ (more risk averse) reduces
the amount of total power offered in the day-ahead; the total amount offered
is larger than would come from just truncating the marginal wind distribution
at the corresponding $\parambeta$ alone, since there are scenarios where a
power shortfall may not result in significant financial loss if the real-time
price is also not high.

For Case 2, increasing $\parambeta$ decreases the price offered at the
right side of the curve when compared to a lower risk aversity. This result
is initially counterintuitive; one would generally expect that a higher
risk aversity would result in the uncertain wind portion in the 100-125MW
range to be priced higher, since that portion will be more likely to be
subject to buyback. However, in Case 2, real-time prices and wind are
positively correlated, meaning that any wind shortfall will actually be
repurchased a lower price in real-time than was generated through clearing
at a higher price in the day-ahead market; put differently, in expectation
there is no financial penalty for underdelivery, and the bulk of the loss
which is targeted by optimizing CVaR is due to insufficient revenue clearing in
the day-ahead market. We can qualitatively justify this claim by examining the
joint distribution plot in Fig. <ref>, which shows the
“active” samples of day-ahead price, real-time price, and wind that
constitute the tail of the profit distribution that CVaR is calculated from.
At the higher $\parambeta$ values, the active samples are clustered where
day-ahead prices, real-time prices, and wind outputs are all lower than the
mean. Since the positive covariance means that there is less penalty for
underdelivery, the offer curve is constructed to clear power that may
not be deliverable in the real-time market to gain the additional day-ahead

Simulation parameters for synthetic Gaussian data cases

Parameter     Symbol     Value

Mean day-ahead price     $\mu_{DA}$     30 $/MWh

Mean real-time price     $\mu_{RT}$     30 $/MWh

Mean wind forecast     $\mu_{wind}$     100 MW

Variance, day-ahead price     $\sigma^2_{DA}$     100

Variance, real-time price     $\sigma^2_{RT}$     100

Variance, wind forecast     $\sigma^2_{wind}$     100

Covariance*, RT price & wind     $cov_{RT, wind}$     $\pm$ 80

Number of scenarios generatred     $N_{\setscenario}$     250

3l*All other covariances are zero.

Offer curve for Case 1.

Offer curve for Case 2.

Active samples for Case 1 optimal solution.

Active samples for Case 2 optimal solution.

Offer curves and pairwise visualizations of distributions.

§ NYISO CASE STUDY

We demonstrate the practicality of our program on real-world price data
from the New York Independent System Operator and synthetic wind forecast data
corresponding to a wind farm in the NYC zone. We compare the profits generated
when the MIQP is run for varying $\parambeta$ against a naive offer strategy
that offers at zero marginal price in the day-ahead market, up to a percentile
quantity of the forecast wind powers $\rvwind$. Hourly-level simulations are run
for the period of October 1st – 31st, 2019.

§.§ Scenario Generation

For each operating day, we generate price samples by taking the historical
day-ahead and real-time prices in the corresponding hour over the
preceeding 50 calendar days. These price scenarios are paired with 100 wind
forecast scenarios generated from historical forecast data using PGscen
<cit.> to create 100 total scenarios for the hour of interest.

§.§ Results

For both the “MIQP CVaR” offer strategy and the “forecast percentile”
strategy, the offer strategy is used to create the day-ahead offer curve. After
the day-ahead clearing quantity is determined from the actual day-ahead price,
the shortfall (or excess) wind production relative to the day-ahead clearing
quantity is bought-back (sold) at the real-time price.

To account for daily variations in revenue, we subtract the total profit from
the day-ahead and real-time markets under these two strategies from an
idealized profit number calculated by selling all the realized wind production
in whichever market is more profitable in hindsight. We refer to this
difference as the regret. For example, if the actual day-ahead hourly price was
$10/MWh, the actual real-time price was $20/MWh, and the actual maximum
wind production was 100 MW for that hour, the ideal profit would be $200 for
that hour. If in that hour the MIQP CVaR strategy yielded a total profit of
$150, the regret would be $\$50 = \$200 - \$150$.

Fig. <ref> shows the total regret over different choices
of risk preference $\parambeta$. There is no directional influence of
$\parambeta$ on the total regret. The MIQP CVaR approach also underperforms
the percentile-based offer strategy in the day-ahead market.

Fig. <ref> shows the average of the standard
deviations of the daily regret. Increasing the $\parambeta$ value (more risk
averse) decreases the variance of the daily regret, as expected. However, the
MIQP CVaR method still underperforms the naive 25% and 50% forecast percentile
methods in both total regret and variance of regret, indicating that there
is still room for improvement in the MIQP CVaR approach, possibly through
improve scenario creation.

Strategy performance considering the sum of regret.

Strategy performance considering the standard deviation of regret.

§ COMPUTATION & SOFTWARE

All computations were performed on a laptop with an Intel Core i5-8300H CPU
(2.30GHz) and 31.2 GB of memory. The model was formulated using JuMP v1.15.1
in Julia v1.9.2. Gurobi v10.0.3 was used to solve the MIQP. Table
<ref> describes the time to solve the model for a range
of scenario sample sizes; this time excludes model creation time. All models
were solved to optimality under the default Gurobi MIP gap of $10^{-4}$.

§ CONCLUDING REMARKS

In this work, we present an MIQP for valuing the CVaR-calculated risk of
different portions of a wind generator's offer curve, along with a pairwise
distribution plot to interpret the shapes of individual offer curves by
examining the set of scenarios that influence the calculation of those curves.
We introduce additional constraints to reduce the integer search space to
improve computation time. Empirical results using synthetic wind forecasts and
real-world NYISO price data demonstrate that adjusting the $\parambeta$ factor
of the MIQP CVaR strategy can adjust the variances in daily profits with little
tradeoff in total profits, although for the example dataset, such a strategy
underperformed a naive forecast percentile offer method in both total
profits and reduced variance of those profits. Such shortcomings could be
because there is insufficient price volatility in today's electricity markets
to reward hedging uncertainties. Further work should examine the price forecast
improvements which can be made to improve the quality of the CVaR-calculated
offer curve in reducing both the total and variance of regret.

MIQP simulation times, Six curve segments

Number of scenarios     Median solve time (s)

50       6.8

100       22.7

250       143.2

500       720.3